使用 scipy.sparse.linalg.eigs

问题描述

鉴于以下输入 numpy 二维数组 A 可以通过文件 hill_mat.npy 使用 following link 检索，如果我只能计算其特征值的一个子集就好了使用像 scipy.sparse.linalg.eigs 这样的迭代求解器。

首先，介绍一下上下文。该矩阵 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 的结果特征值是 complex（以复共轭对的形式出现），我对 最小虚部 更感兴趣模数。

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

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 来使用所谓的 shift-invert 模式让算法知道在哪里可以找到这些特征值。尽管如此，我所有的尝试都没有产生任何结果......

有没有对这个功能更有经验的人帮我完成这项工作？

以下是完整的代码片段：

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,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),return_eigenvectors=False)

解决方法

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

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

这样，在寻找小的特征值时，为了帮助 LAPACK 完成工作，应该通过指定适当的 sigma 值来激活 shift-invert 模式，同时还反转 {{ 中指定的所需指定子集1}} 关键字参数。

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

which

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