python - 寻找大型稀疏矩阵的最小特征向量，在 SciPy 中比在 Octave 中慢 100 倍以上

Question

我正在尝试计算几个（5-500）个特征向量，这些特征向量对应于大型对称方形稀疏矩阵（高达 30000x30000）的最小特征值，其中小于 0.1% 的值非零。

我目前在 shift-invert 模式（sigma=0.0）中使用 scipy.sparse.linalg.eigsh，我通过有关该主题的各种帖子发现这是首选的解决方案。但是，在大多数情况下，最多需要 1 小时才能解决问题。另一方面，如果我要求文档中预期的最大特征值（我系统上的亚秒），该函数非常快。

由于我在工作中更熟悉 Matlab，因此我尝试在 Octave 中解决问题，使用 eigs (sigma=0) 在短短几秒钟内（低于 10 秒）给了我相同的结果。由于我想对包括特征向量计算在内的算法进行参数扫描，因此在 python 中也能获得这种时间增益。

我首先更改了参数（尤其是容差），但在时间尺度上并没有太大变化。

我在 Windows 上使用 Anaconda，但尝试将 scipy 使用的 LAPACK/BLAS（这是一个巨大的痛苦）从 mkl（默认 Anaconda）切换到 OpenBlas（根据文档由 Octave 使用），但看不到变化表现。

我无法弄清楚，使用的 ARPACK 是否有什么需要改变的（以及如何改变）？

我将以下代码的测试用例上传到以下保管箱文件夹： https ://www.dropbox.com/sh/l6aa6izufzyzqr3/AABqij95hZOvRpnnjRaETQmka?dl=0

在 Python 中

import numpy as np
from scipy.sparse import csr_matrix, csc_matrix, linalg, load_npz   
M = load_npz('M.npz')
evals, evecs = linalg.eigsh(M,k=6,sigma=0.0)

在八度：

M=dlmread('M.txt');
M=spconvert(M);
[evecs,evals] = eigs(M,6,0);

任何帮助都会得到帮助！

我根据评论和建议尝试了一些其他选项：

八度： eigs(M,6,0)给eigs(M,6,'sm')我同样的结果：

[1.8725e-05 1.0189e-05 7.5622e-06 7.5420e-07 -1.2239e-18 -2.5674e-16]

同时eigs(M,6,'sa',struct('tol',2))收敛到

[1.0423 2.7604 6.1548 11.1310 18.0207 25.3933] 

快得多，但前提是容差值大于 2，否则它根本不会收敛并且值差异很大。

Python： eigsh(M,k=6,which='SA')并且eigsh(M,k=6,which='SM')两者都不收敛（未达到收敛时出现 ARPACK 错误）。仅eigsh(M,k=6,sigma=0.0)给出一些特征值（几乎一个小时后），这些特征值与最小的八度音阶不同（甚至发现了 1 个额外的小值）：

[3.82923317e-17 3.32269886e-16 2.78039665e-10 7.54202273e-07 7.56251500e-06 1.01893934e-05]

如果容差足够高，我也会从中得到结果eigsh(M,k=6,which='SA',tol='1')，这接近于其他获得的值

[4.28732218e-14 7.54194948e-07 7.56220703e-06 1.01889544e-05, 1.87247350e-05 2.02652719e-05]

再次使用不同数量的小特征值。计算时间仍然接近 30 分钟。虽然不同的非常小的值可能是可以理解的，因为它们可能代表 0 的倍数，但不同的多重性让我感到困惑。

此外，SciPy 和 Octave 之间似乎存在一些根本差异，我还无法弄清楚。

标签： pythonscipyoctavesparse-matrixeigenvector

Answer 1

5 月 19 日添加：Cholesky 内部求解器：

scipy eigsh的文档 说

shift-invert 模式 ... 需要一个算子来计算线性系统的解 (A - sigma * I) x = b... 这通过一个显式矩阵的稀疏 LU 分解 (splu) 在内部计算，或者通过一个一般线性算子的迭代求解器计算。

ARPACK 多次调用此“内部求解器”，具体取决于tol等；显然，慢内部求解器 => 慢eigs。对于Aposdef， sksparse.cholmod 比 slu 快得多。

Matlab eig 也使用 cholesky：

如果 A 是 Hermitian 而 B 是 Hermitian 正定的，那么算法的默认值是 'chol'

---

Fwiw，在我的旧 4 核 imac 上在一小时内np.linalg.eigh找到7 Gb 密集矩阵的 所有特征值和特征向量——哇。A.A它的光谱是这样的：

---

2020 年 2 月，TL;DR

一个猜想和一些评论，因为我没有 Matlab / Octave：

要找到特征值 >= 0 的对称矩阵的小特征值，就像您的一样，以下方法比移位反转快得多：

# flip eigenvalues e.g.
# A:     0 0 0 ... 200 463
# Aflip: 0 163 ... 463 463 463
maxeval = eigsh( A, k=1 )[0]  # biggest, fast
Aflip = maxeval * sparse.eye(n) - A
bigevals, evecs = eigsh( Aflip, which="LM", sigma=None ... )  # biggest, near 463
evals = maxeval - bigevals  # flip back, near 463 -> near 0
# evecs are the same

eigsh( Aflip )大特征对比小特征对的移位反转更快，因为A * x它比solve()移位反转必须做的要快。Aflip在使用 Cholesky 进行正定快速测试后，可以想象 Matlab / Octave 可以自动执行此操作。
你可以eigsh( Aflip )在 Matlab / Octave 中运行吗？

可能影响准确性/速度的其他因素：

Arpack 的默认起始向量v0是随机向量。我使用v0 = np.ones(n)，这对某些人来说可能很糟糕，A但可以重现 :)

这个A矩阵几乎完全是奇异的，A * ones~ 0。

多核：带有 openblas / Lapack 的 scipy-arpack 在我的 iMac 上使用了 4 个内核中的 3.9 个；Matlab / Octave 使用所有内核吗？

---

以下是几个kand的 scipy-Arpack 特征值tol，从gist.github下的日志文件中提取：

k 10  tol 1e-05:    8 sec  eigvals [0 8.5e-05 0.00043 0.0014 0.0026 0.0047 0.0071 0.0097 0.013 0.018] 
k 10  tol 1e-06:   44 sec  eigvals [0 3.4e-06 2.8e-05 8.1e-05 0.00015 0.00025 0.00044 0.00058 0.00079 0.0011] 
k 10  tol 1e-07:  348 sec  eigvals [0 3e-10 7.5e-07 7.6e-06 1.2e-05 1.9e-05 2.1e-05 4.2e-05 5.7e-05 6.4e-05] 

k 20  tol 1e-05:   18 sec  eigvals [0 5.1e-06 4.5e-05 0.00014 0.00023 0.00042 0.00056 0.00079 0.0011 0.0015 0.0017 0.0021 0.0026 0.003 0.0037 0.0042 0.0047 0.0054 0.006
k 20  tol 1e-06:   73 sec  eigvals [0 5.5e-07 7.4e-06 2e-05 3.5e-05 5.1e-05 6.8e-05 0.00011 0.00014 0.00016 0.0002 0.00025 0.00027 0.0004 0.00045 0.00051 0.00057 0.00066
k 20  tol 1e-07:  267 sec  eigvals [-4.8e-11 0 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 5.8e-05 6.4e-05 6.9e-05 8.3e-05 0.00011 0.00012 0.00013 0.00015

k 50  tol 1e-05:   82 sec  eigvals [-4e-13 9.7e-07 1e-05 2.8e-05 5.9e-05 0.00011 0.00015 0.00019 0.00026 0.00039 ... 0.0079 0.0083 0.0087 0.0092 0.0096 0.01 0.011 0.011 0.012
k 50  tol 1e-06:  432 sec  eigvals [-1.4e-11 -4e-13 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 ... 0.00081 0.00087 0.00089 0.00096 0.001 0.001 0.0011 0.0011
k 50  tol 1e-07: 3711 sec  eigvals [-5.2e-10 -4e-13 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 ... 0.00058 0.0006 0.00063 0.00066 0.00069 0.00071 0.00075

versions: numpy 1.18.1  scipy 1.4.1  umfpack 0.3.2  python 3.7.6  mac 10.10.5 

Matlab / Octave 差不多吗？如果不是，所有的赌注都被取消——首先检查正确性，然后检查速度。

为什么特征值摆动如此之大？对于一个所谓的非负定矩阵来说，微小的 < 0 是 舍入误差的标志，但通常的微小移位技巧A += n * eps * sparse.eye(n), 并没有帮助。

---

这A来自哪里，什么问题领域？你能生成相似A的、更小的或更稀疏的吗？

希望这可以帮助。