使用 scipy.sparse.linalg.eigs 计算复杂特征值

Question

鉴于以下输入 numpy 2d-array A 可以通过文件 following link 检索到 hill_mat.npy，如果我可以使用像 @ 这样的迭代求解器仅计算其特征值的一个子集，那就太好了987654322@.

首先，有一点上下文。该矩阵A 是由大小为N 的二次特征值问题产生的，该问题已在双大小2*N 的等效特征值问题中线性化。 A 具有以下结构（蓝色为零）：

plt.imshow(np.where(A > 1e-15,1.,0), interpolation='None')

以及以下功能：

A shape = (748, 748)
A dtype = float64
A sparsity ratio = 77.64841716949297 %

A 的真实尺寸比这个可重现的小例子要大得多。对于这种情况，我预计真实的稀疏率和形状接近 95% 和 (5508, 5508)。

A 的结果特征值是复数（以复共轭对的形式出现），我对模数中虚部最小的那些更感兴趣。

问题：使用直接求解器时：

w_dense = np.linalg.eigvals(A) 
idx = np.argsort(abs(w_dense.imag))
w_dense = w_dense[idx]

计算时间很快变得令人望而却步。因此，我希望使用稀疏算法：

from scipy.sparse import csc_matrix, linalg as sla
w_sparse = sla.eigs(A, k=100, sigma=0+0j, which='SI', return_eigenvectors=False)

但似乎 ARPACK 没有以这种方式找到任何特征值。从scipy/arpack tutorial 中，当寻找像which = 'SI' 这样的小特征值时，应该使用所谓的移位反转模式，指定sigma kwarg，ie，以便算法知道在哪里它可以期望找到这些特征值。尽管如此，我所有的尝试都没有产生任何结果......

是否有对此功能更有经验的人帮我完成这项工作？

下面是完整的代码sn-p：

import numpy as np
from matplotlib import pyplot as plt
from scipy.sparse import csc_matrix, linalg as sla

A = np.load('hill_mat.npy')
print('A shape =', A.shape)
print('A dtype =', A.dtype) 
print('A sparsity ratio =',(np.product(A.shape) - np.count_nonzero(A)) / np.product(A.shape) *100, '%')

# quick look at the structure of A
plt.imshow(np.where(A > 1e-15,1.,0), interpolation='None')

# direct
w_dense = np.linalg.eigvals(A)
idx = np.argsort(abs(w_dense.imag))
w_dense = w_dense[idx]

# sparse
w_sparse = sla.eigs(csc_matrix(A), k=100, sigma=0+0j, which='SI', return_eigenvectors=False)

Answer 1

问题终于解决了，我想我应该更仔细地阅读文档，但是，以下内容非常违反直觉，我认为可以更好地强调：

... ARPACK 包含一种允许快速确定 非外部特征值：移位反转模式。如上所述，这 模式涉及将特征值问题转换为等价的 不同特征值的问题。在这种情况下，我们希望找到 特征值接近零，所以我们将选择sigma = 0。变身的 然后特征值将满足，因此我们的小特征值 变大 特征值。

这样，在寻找小特征值时，为了帮助 LAPACK 完成工作，应该通过指定适当的 sigma 值来激活移位反转模式，同时反转 which 关键字中指定的所需指定子集论据。

因此，这只是一个执行的问题：

w_sparse = sla.eigs(csc_matrix(A), k=100, sigma=0+0j, which='LM', return_eigenvectors=False, maxiter=2000)
idx = np.argsort(abs(w_sparse.imag))
w_sparse = w_sparse[idx]

因此，我只能希望这个错误对其他人有所帮助:)