实对称矩阵本征值问题的理论与实践

18 Jul 2019

引言
Lapack 几种实对称矩阵本征值算法的异同
调用 Lapack 的几种不同姿势
大矩阵对角化编程实践中的问题
一些更深度的分析讨论
References

引言

本文重点关注求大矩阵的全部本征值（或包括本征向量）这一数值问题（如果你对稀疏矩阵本征值问题感兴趣，或是不熟悉本文提到的各种线性代数库的关系，可以参考我之前的文章）。这一问题乍一看当然很简单，无论是 numpy 的 np.linalg.eigh()，抑或是 mathematica 的 Eigensystems[] 都可以一行轻松解决。但是当我们想把结果 push to the limit，想算的矩阵越大越好时，是否内存耗尽就是尽头了，还有没有办法占用更小的内存，有更快的速度呢。读者的第一反应可能是无法提升了，原因很简单，numpy matlab 或者 mathematica 都内置或链接了优化过的 Lapack 库，比如 MKL，因此理论上本征值问题，调用这些高级语言和直接调用 MKL，应该一样快，占的内存一样多。然而现实并没有那么简单，这些底层的 Lapack routine 有大量的变种和可调参数，这些在高级语言甚至 C++ 的各种线性代数库中，都被封装得不见了，然而这些正是在现有硬件条件下可以对角化更大矩阵的精髓。例如，通过本文我们可以看到，我们通过直接调用 lapack routine，不需要太多 hack，就能把本征值问题所需的内存压缩到 numpy 同样规模问题所需内存的四分之一。

Lapack 几种实对称矩阵本征值算法的异同

Lapack 对于实矩阵的本征值问题提供的 routine，可以参考¹。更多关于本征值问题和三对角化的 routine 也可以参考²最后的表格。其中第一个字母的 d s 对应 double 和 float32，表格中已经省略了。第二个字母的 s 是对称矩阵的意思（h 代表 Hermitian matrix）。第三个字母的 y 代表对称矩阵的正常存储方式，p 代表 packed storage，b 代表 band matrix，而 t 代表 tridiag matrix。下面我们会看到，对于所有普通的实对称矩阵，Lapack 内部的对角化算法都是两步走，第一步是通过 Household transformation 将对称矩阵 tridiagolization，第二步是选择适当的算法，对三对角化矩阵求本征值。因此我们只关注一般实对称矩阵的本征值 routine，其他的都是同理的，只是矩阵的排布和存储方式的区别。当然，如果 state-of-the-art 对角化的瓶颈是内存的话，可以考虑专门使用 packed 排列的矩阵来对角化，这样可以节约一半的内存，但是速度可能有略微的损失³。后三个字母则代表了 routine 利用的算法。也就是说 Lapack 提供的本征值问题的算法共四种，分别是 ev, evd,evr,evx, 这也是这部分要讨论的问题：究竟选择哪个 routine，他们各自的特点是什么。这里的特点又分为两个方面，一是功能性上的区别，比如四个算法中有些可以计算部分本征值，有些只能算全谱。另一个是性能和准确性上的区别，不同方法需要的时空复杂度，和最后结果的精度不尽相同。关于这些算法最详细的对比，可以参考⁴。

简单来讲，这四个 routine 对应的算法分别为 ev, QR iteration，simple driver，计算全部本征值和本征态（可选）；evx，Bisection and Inverse Iteration, expert driver, 计算全部或部分本征值和本征态（可选）；evd，divide and conquer，计算全部本征值和本征态（可选），对于大矩阵明显比 simple driver 更快，内存空间需求更大；evr，Multiple Relatively Robust Representations，计算全部或部分本征值和本征态（可选），在大多数情况最快，且需要最少的内存空间。其中，最早的 Lapack 只有前两种算法，后两种是随着研究的进展而增加的，并且 Lapack 有用 evr 算法统一的想法⁵。当然随着数值计算研究的新进展，很可能有更多的算法或者算法改进进入 Lapack。换句话说，Lapack 听起来很像祖传老代码，但是实际上一直保持在人类数值算法知识的前沿，在线性代数实践可用的稳定算法上，几乎是无法打败的。MKL 等实现做的事，只不过是考虑如何将 Lapack 算法的实现更多考虑硬件的特性，比如 AVX 指令集的利用或者 L3 Cache 的更好利用等等。

稍微总结一下这四个 routine 在时空复杂度和准确度上的区别。

evd 和 evr 大矩阵时远远比其他两个进程快，因此另两个算法，无特殊情况就不需要考虑了。有趣的是，注意到 Lapack 以前只有这两个快被淘汰的 routine，这充分说明了这一数值领域还是比较 active，一直有新的算法和改进出来，并且 Lapack 也能跟进，基本保持在实践上可行的学术前沿。evd 和 evr 具体谁更快，取决于矩阵的具体形式。
evr 对角化所需的浮点数操作最少（不代表其最快，因为其 FLOP rate 也很低）。其浮点数操作吞吐低，可能是因为除法较多。
ev 和 evd 可以给出最准确的结果，准确性可以通过不同本征向量之间的正交性，和 $\sum_{x}\vert Ax-\lambda x\vert/\vert A\vert$ 这一残差来刻画。evr 残差的经验上限大概在 $80n\epsilon$ 左右，其中 n 是矩阵的维度，$\epsilon$ 是数值精度。而 evd 和 ev，对应的残差则不到一个 $n\epsilon$。
对于需要计算全部本征向量的问题，evd 需要的工作空间是 2 个额外同大小矩阵的数量级，而其他进程则只需要 1 个额外矩阵的数量级。
值得一提的是，虽然这些求本征值算法 worst case 的时间复杂度是 $O(n^3)$，但实践中，对于比较“正常”的矩阵，其 scaling 指数很可能小于 3。对于非常稀疏的矩阵，evd 给出的指数甚至只有 2.5 左右。对于 evr 算法，其对角化 tridiagonal 矩阵的步骤，理论复杂度只有 $O(n^2)$ （前边常数取决于本征值谱的“稠密”程度），不过其实践上很多时候并不比 evd 快。这里讨论的时间，都包括了第一步三对角的部分，其复杂度为 $O(n^3)$。上述第二步对角化步骤的复杂度分析指的是求本征向量的情形，求本征值的话，这一步基本上是 $O(n^2)$。

基本建议就是，尽量选择 evd 或者 evr 做对角化，具体用哪个需要根据问题的矩阵，进行真实的 benchmark 来测量时间和准确程度是否能够接受（evr 对于精度的影响是否会影响到解决问题，一定要提前试验清楚）。

调用 Lapack 的几种不同姿势

首先需要指出的是 Lapack 的代码基本上是 Fortran 写成的，包括 MKL 提供的 Lapack 实现也是 Fortran。因此 Fortran 使用 Lapack 就没什么好说的了，直接无脑调用，记得编译时连接 mkl 库就行。

在 C 里调用 Lapack 就有几种不同的途径。一种是直接调用，只要确保 C 中定义的函数原型头文件与 Fortran routine 相同，那么连接 Fortran 编译成的库，则可直接在 C 代码中调用 Fortran 函数（需要注意的时，根据编译器的不同，name mangling 的机制也不同，由此导致 C 中调用的 Fortran 函数的名字也和 Fortran 中定义的名字有所不同，可能会涉及大小写的变化和下划线的添加等等）。当然使用 mkl 直接调用 Fortran Lapack routine 更简单，因为脏活都被 Intel 封装好了，我们只需要也按照名字，比如 dsyevd 在 C 中直接调用即可。需要注意的是 Fortran 的函数原型，所有变量总是传址的，因此在 C 中要注意所有函数变量都是地址而非值。

此外，还有更多底层的 Lapack 的简单 wrapper 可以使用，最著名的就是 LAPACKE，mkl 中也支持通过 LAPACKE 的界面在 C 中调用 Lapack。其函数名称为 LAPACKE_dsyevd 的格式。关于 Lapack 谜一样的函数命名的规则，可以参考⁶。这种封装，自动会完成最优 lwork 的查询和在堆上开出所需的数组空间，而不需要自己两次调用 lapack routine。同时这种封装可以说完全没有任何性能上的牺牲，因此推荐直接使用 LAPACKE 的 interface，在 C 中使用。需要注意 Lapacke 和 Lapack 的 API 不相同，不要混淆。

除了这种底层的 Lapack wrapper，还有很多 C++ 的库提供封装比较好的线性代数库，比如 armadillo 等，其底层也可以设置调用相应的 Lapack 后端，算是从 C++ 中使用 Lapack。当然，C++ 也可以直接调用 Lapack routine，不过由于 C++ 函数重载的原因，其编译器又有自己的 name mangling，因此要对外部 Lapack 函数记得声明 extern "C"。

对于 Python，线性代数的处理几乎全部依赖于 numpy 和 scipy 两个库。对于 dense matrix 全部本征值和本征向量的求解，numpy.linalg.eigh 对应的底层 Lapack routine 是 ?syevd (调用源码)，而 scipy.linalg.eigh 的底层 routine 是 ?syevr(调用源码)。scipy 之所以选用了 syevr routine，是因为该算法，可以有选择的指定求部分本征谱。关于两个 routine 的区别，上一节已经简要论述过。

至于 matlab 或者 mathematica 这种更高级的语言，其实也是内置了 MKL 或者类似的 Lapack 优化实现，只不过这一链接不需要手动配置，用起来更傻瓜式而已。这也是为什么会有 “matlab 矩阵和线性代数计算更快” 的无稽之谈，这都源于对于数值计算和计算机不熟悉的人的口口相传。事实上，只要你有足够的经验，连接 MKL 的任何语言的 interface 都可以有不输 Matlab 的速度，只不过这可能需要你了解很多编译器优化，多线程，动态库连接等知识。换句话说，大部分人无法顺利和高效的连接到 MKL 库使用，这就使得 Matlab 这种提前配置连接好 MKL 的软件赢得了“算矩阵比较快”的虚名。这种无知简直是用自己对编译连接的不熟悉在抹杀大量算法科学家对于 Lapack 和大量 Intel 工程师对于 MKL 的贡献。在线性代数底层算法方面，Matlab 做的仅仅是正确的配置和链接了 MKL 而已。

大矩阵对角化编程实践中的问题

你以为学过些 numpy 或者 lapack，算过 100*100 矩阵的本征值，再算更大的矩阵的本征值就是改个 size 的简单问题？完全不是，当你尝试计算更大维度的矩阵对角化问题的时候，无数的软件问题会在真实内存耗尽之前阻止你能够计算的矩阵的大小。而这些问题，基本都需要扎实的计算机体系结构，数值方法，和 linux 下 C 开发的基本功，才能完美的定位分析和解决问题，试看下面几例。

C 直接调用 MKL lapack routine 计算本征值时，dimension 大于 1020 的矩阵出现 Segmentation fault (syslog 显示 segfault error 6)。

原因：由于 linux 上默认的 stack size limit。linux 上 ulimit -s 可以查看，给每个用户默认添加了 8192k 的 soft limit。这主要是为了防止用户代码没写好，函数栈无限迭代下去。然而对于 Lapack 计算，需要在 C 中开一个数组 double a[N*N] 来存放矩阵。$1020^2\times 8/1024=8128k$，再加上其他一些变量，就触及了 stack size 的限制，因此会产生 segfault。

解决：命令行直接 ulimit -s unlimited 之后在运行即可。这样做就取消了 stack size 的 soft limit。~~2019 年了，操作系统默认的 stack size limit 还这么小，我是无法理解的~~。
C 直接调用 MKL LAPACKE 封装或直接调用 lapack routine dsyevd 计算本征值时， dimension 大于 32766 的矩阵出现 Not enough memory to allocate work array in LAPACKE_dsyevd。

原因：export MKL_VERBOSE=1 可以在 stdout 显示 MKL 库调用的详细信息。再次运行程序，我们可以看到，LAPACKE 第一次调用 lapack routine 做 size query 是成功了，但第二次调用做真实计算没有显示，因此调用失败了。这一问题源自 Lapack routine 函数原型的变量声明。很多 Lapack routine 有诸如 lwork 这种整数变量，用来指定 work space 的大小。但问题就出在这个整数变量上，用 C 的语言说就是，lwork 的变量类型是 int 还是 long 呢。也就是说 lwork 最大只能到 $2^{31}-1$ 还是可以到 $2^{63}-1$ 呢（lwork 默认的变量类型似乎是有符号的整形）。当矩阵维度为 32767 时，dsyevd 需要的 lwork 的值是 2N^2+6N+1=2^31-1+65534 已经超过了 int 的上限。那么该值传入 fortran 函数时会强制类型转换，因此对于 Lapack routine，其发现 lwork 传入了一个负数，因此很自然的报错。

解决：需要连接支持 64 位的整数类型的 intel 库，也即 ilp64 interface。具体的说就是：
```
icc  -DMKL_ILP64 -I${MKLROOT}/include src.c   -L${MKLROOT}/lib/intel64 -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm -ldl
```
对于其他的编译器或情形，可以参照⁷进行选择，需要注意 Select interface layer 中选择 64-bit integer。上述这行命令中，-D MKL_ILP64 定义了该宏变量，从而使得程序中的所有 MKL_INT 都被替换为 long 而非 int。同时还需注意要链接到 libmkl_intel_ilp64.so 这一动态库，从而使得 Lapack routine 接受的参数也是 long 型。此二者缺一不可。

本来对于 icc 的编译器，有默认的 -mkl flag 来实现连接到 MKL 的编译，但其并不适用于 64 bit 整数的界面。这一点可以通过 ldd a.out 查看输出二进制的动态库依赖来看出。直接用 mkl flag 编译的文件，依赖的是 mkl_intel_lp64 的动态库，而 lp64 对应的是 32bit 整数的 lapack 函数原型。因此必须手加 ilp64 的连接。icc mkl flag 这一行为，我认为是应该进行修改的。因为 mkl 早就提出了 libmkl_rt.so 这一动态库⁸，作为 one for all 的连接库使用，而通过 runtime 设置默认的环境变量，来改变整数定义的字节数，这一行为应为 icc -mkl 的默认行为才是更加合适的。

在突破了 32766 之后，我们还有个要注意的问题。通常 C style 是将矩阵的维度定义为宏变量，比如 #define N 50000，但这样做编译也是会报错的。理由很简单，定义矩阵数组的 $50000*50000$ 已经爆掉了 int （没错，int 就是这么小）。所以，我们需要强行注明 50000 是 long 型，也即 50000UL。这就是为什么说 C++ 的 const 比 C 的 macro 替换安全的多，macro 一个文本替换，somehow 还和变量类型纠缠在一起，比较危险。
使用 armadillo eig_sym() 间接调用 lapack dsyevd 计算本征值时，dimension 大于 23169 的矩阵出现 Intel MKL ERROR: Parameter 8 was incorrect on entry to DSYEVD 或者 error: arma::memory::acquire(): requested size is too large terminate called after throwing an instance of 'std::logic_error' what(): arma::memory::acquire(): requested size is too large Aborted (core dumped)。

原因：这一错误可以通过 gdb 很容易查到，并发现是上一个错误的变种。设置 ulimit -c unlimited 来允许转存 core 文件，并且打开 -g flag 编译源码。之后 segfault 后，我们就可以通过 gdb -c core a.out 来检查出问题的地方。通过 bt 来查看函数调用栈，并且在可疑的地方通过 b arma::aux_lib::eig_sym_dc 设置断点，从而 n 进行逐步检查，我们发现 armadillo 源码中写的 lwork space 比 lapack 需要的多了一倍，也即 blas_int lwork = 2 * (1 + 6*N + 2*(N*N)) 这句，使得 armadillo 在 23169 而不是 32767 维出问题。armadillo 的这种设计，可能是源自“更大的 workspace 可以使 Lapack 对角化计算更快”这一信念（但我在实践中无法发现不同 workspace 对角化时间的区别）。同时我们注意到 p sizeof(blas_int) 的结果为 4。因此当传入 lapack::syevd 时， lwork 变量面临 int overflow 的问题，从而报错。只不过 armadillo 的异常封装反而使得错误变了样子，改报内存不足了。

解决：这一解决分为两方面。第一方面是，我们得让 armadillo 的 blas_int 用 long 替换，而非 int。这可以通过修改 armadillo include 中的 header config.hpp 实现，也即将其中 #define ARMA_BLAS_LONG 这句宏定义去掉注释。需要注意，这一改动，需要在 make install 之后，安装位置处的 armadillo 文件夹内直接进行。对源码进行改动的话，make install 之后位置的头文件并不会改变（这一行为的原因我还不理解）。这一宏替换可以在 runtime 轻松更换，从而使得 armadillo 的灵活性更好，而不需重新安装。

解决的第二方面是和上一个问题一样，我们的程序需要连接 ipl64 的 mkl 库。本来默认对于安装过的 armadillo 我们只需编译时加 -larmadillo 即可，不需要关心 lapack 库的动态链接。但现在，为了使用 ipl64，我们还是可以按照上面连接 mkl 的方式加 flag，多补上一句 -larmadillo 即可。其实我们也可以不这么麻烦，对于只加了 larmadillo 编译出的二进制，ldd 发现其依赖的是 libmkl_rt.so，也即是 one for all 的 mkl 动态库。那么我们还是只加 larmadillo 编译就好，为了使用 ipl64 界面，唯一需要做的就是⁹ export MKL_INTERFACE_LAYER=ILP64 即可。也即 mkl_rt 这一动态库的好处就是，可以把线程库和整数bit数等不同的选择，推迟到 runtime，只需要设定相应的环境变量即可。
numpy 利用 np.linalg.eigh() 计算本征值时（间接调用 dsyevd），dimension 大于 32766 的矩阵出现 ValueError: On entry to DORMTR parameter number 12 had an illegal value, 或 Segmentation fault (core dumped) 或快速返回全为零的错误结果。

原因：详细的前因后果可以参考我的这个 issue，hopefully 这个问题可以在不远的未来解决，（虽然 numpy 开发者已经拖了四五年了）。看到出问题的数字就知道，这一问题也来自于 lwork 参数的 int overflow。只不过 numpy 里的问题更难解决一些，并不能通过简单的编译连接到64位整数的 MKL 库就可以解决。这是因为，numpy 调用 lapack 的胶水层，有大量代码，并且写死了 int 类型。因此只能等未来 numpy 的代码变更，才能真的解决 numpy 里无法对角化大矩阵的问题。（其实我有点不可想象，2019年了，python 用途最广泛，甚至是几乎唯一的线性代数解决方案中，还不能对角化一个三万多维的矩阵，而且这一软肋还没有多少人知道！）

想在 python 中对角化更大的矩阵的一个并不是很理想的 workaround，是使用 scipy 的 scipy.linalg.eigh()。其之所以能够突破 32766 这一屏障，并不是因为 scipy 的胶水层支持 long 型整数传入 Lapack 进程，而是因为其调用的是 dsyevr 而非 dsyevd，后者所需的工作空间的参数 lwork 等，都在 N 的数量级，而非 $N^2$ 数量级，因此传入 Lapack 的参数并不会有 int overflow 的问题，从而实现对角化。之所以这个 workaround 不太理想，在于有些情况 dsyevr 计算本征值问题要比 dsyevd 慢一些。

一些更深度的分析讨论

benchmark 不同的 Lapack 算法，或是不同层级的调用，观察所需的内存空间大小时，需要注意以下细节。Linux 操作系统在计算程序占用的内存时（至少 htop 和 free 等监控工具如此），在堆上 malloc 开的内存区域，不会被立即计入到程序占用内存中。只有在对应的位置被赋值后，才会算入程序内存。这也就解释了观察调用 Lapack 求本征值的 routine 时，很多时候，只有在运行的最后阶段，内存才会突然增加到和分配的 workspace 一个数量级，因为只有那时候 malloc 的工作空间，才真的被赋值和使用上。

通过 numpy 计算对角化问题，虽然底层调用了 Lapack 的 ?syevd routine，但是肯定无法满足 push to the limit 的情形。其阻碍是两重的，其本质问题来自 np.linalg.eigh(m) 之后，矩阵 m 还完整的存在，因此 numpy 在调用 Lapack routine 之前，需要完整的在内存中复制一遍 m 矩阵。因此时间上，numpy 要比 Fortran 直接调用 Lapack 多了一个矩阵复制的时间，当矩阵大小达到几十 GB 这个量级的时候，这一复制时间不可忽略。在空间上问题更明显，如果在原始矩阵上操作，那么 Lapack 的 dsyevd 只需要额外 $2N$ 的空间，就可以求得矩阵全部的本征谱，因此，求得全部本征谱的问题所需的内存为 $8(N^2+2N)$ （double 数据类型）。然而由于 numpy 为了保留原始矩阵，则需要的内存为 $8(2N^2+2N)$，相当于对角化同样的矩阵，python 比 Fortran 需要的内存整整大了一倍。因此，现在的 numpy 肯定是不适合 start-of-the-art 的对角化数值计算的。但这个锅本质上不是 python 的，我们完全可以设计一个不保留原始矩阵，直接将 python 的原始矩阵地址传入 Lapack routine 的 python interface，这样的新函数，将不会有时间和空间的 overhead，从而真正达到时间空间复杂度和底层 Lapack 一致。这一提议参见我在 numpy 提的 issue。注意到 scipy 虽然提供了所谓 low level 的 lapack wrapper，但并不能解决本段的问题。这是因为某种程度上来说，scipy 提供的 wrapper 还是有点 high level 了，其依旧保留了输入的矩阵，也即也有复制的过程，因此无助于解决该问题。该节的分析，可以完美解释矩阵对角化时 numpy 和 Lapack 对应的内存占用的区别。类似的，C++ 的 armadillo 之类的封装，也存在相应的问题。一个最简单的判别标准就是，任何抽象层级的一个线性代数库，当计算完本征值时，原矩阵还存在的话，那么该库一定是无法榨干硬件性能的，因为其浪费了一半内存。更进一步地，如果我们在 Lapack 中，选择 packed storage 的方式表示实对称矩阵，那么同样的问题所需的内存还可进一步减少一半。

最后是一些关于算法选择的讨论。为什么求全部能谱的时候，不能用稀疏矩阵表示？这样做貌似可以大量节约内存，毕竟我们关心的矩阵还是非常“稀疏”的。对于直接使用稀疏矩阵特有的本征值求解方法，比如 Lanczos 算法等等，大可以 argue 说，这些算法只适合求少量本征值，对于求大量本征值的问题，这些基于 subspace 的算法，速度上既没有优势，准确性上又堪忧，存在数值的不稳定性。那么，为什么不能采取 dense matrix 的本征值算法结合 sparse matrix 的矩阵表示呢，这样岂不是既省内存，又不会太慢？答案是，dense matrix 矩阵对角化算法的第一步，tridigonalization 利用了 Household transformation ¹⁰，而这种变换，会使得原本非零矩阵元很少的矩阵也变得非零矩阵元非常多，这就是所谓的 fill in （具体地说，由于每一步的变换，都是对各列向量做一个镜面变换，当该镜面的方向并不特殊的时候，原来只有少量非零分量的向量，将镜面映射到几乎所有分量都不为零的向量，由此使得变换后的矩阵 fill in，也即非零元素快速增加）。具有比较强 fill in 性质的算法，都无法应用于稀疏矩阵。因为经过这种算法，矩阵总是变得很“密”，稀疏表示将不会有内存和时间优势。也就是 dense matrix 的对角化算法无法适用于 sparse matrix。

References