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

Complex eigenvalues computation using scipy.sparse.linalg.eigs

给定以下输入 numpy 二维数组 A 可以使用 following link through the file hill_mat.npy, it would be great if I can compute only a subset of its eigenvalues using an iterative solver like 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 的最终特征值是 复数 （以复数共轭对形式出现），我对 最小的特征值更感兴趣模数中的虚部。

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

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 使用所谓的 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,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)

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

... ARPACK contains a mode that allows a quick determination of non-external eigenvalues: shift-invert mode. As mentioned above, this mode involves transforming the eigenvalue problem to an equivalent problem with different eigenvalues. In this case, we hope to find eigenvalues near zero, so we’ll choose sigma = 0. The transformed eigenvalues will then satisfy , so our small eigenvalues become large eigenvalues .

这样，在寻找小特征值时，为了帮助 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]

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