时间戳

首页 > 解决方案 > 使用 scipy.sparse.linalg.eigs 计算复杂特征值

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

问题描述

鉴于以下输入 numpy 2d-arrayA可以通过文件通过以下链接hill_mat.npy检索，如果我可以使用像scipy.sparse.linalg.eigs这样的迭代求解器仅计算其特征值的一个子集，那就太好了。

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

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

A_imshow.png

以及以下特点：

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 教程中，当寻找像 的小特征值时，应该通过指定kwargwhich = 'SI'来使用所谓的移位反转模式，即为了让算法知道它可以在哪里找到这些特征值。尽管如此，我所有的尝试都没有产生任何结果......sigma

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

下面是一个完整的代码片段：

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)

标签： pythonnumpyscipyeigenvaluearpack

解决方案

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

... 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]

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

推荐阅读