wpsze

Open Menu 实对称矩阵本征值问题的理论与实践 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，可以参考 1 。更多关于本征值问题和三对角化的 routine 也可以参考 2 最后的表格。其中第一个字母的 d s 对应 double 和 float32，表格中已经省略了。第二个字母的 s 是对称矩阵的意思（h 代表 Hermitian matrix）。第三个字母的 y 代表对称矩阵的正常存储方式，p 代表 packed storage，b 代表 band matrix，而 t 代表 tridiag matrix。下面我们会看到，对于所有普通的实对称矩阵，Lapack 内部的对角化算法都是两步走