使用 Python 对三对角对称稀疏矩阵进行对角化

Question

我有一个由 Python 代码计算的 NxN 对称和三对角矩阵，我想将它对角化。

在特定情况下，我正在处理N = 6000，但矩阵可以变得更大。由于它是稀疏的，我认为对角化它的最佳方法是使用算法scipy.sparse.linalg.eigsh()，它在我使用的其他稀疏和对称矩阵（但不是三对角矩阵）中表现得非常好。特别是，由于我只需要频谱的低位部分，因此我在函数中指定了 k=2 和 which='SM'。

但是，在这种情况下，该算法似乎不起作用，因为经过大约 20 分钟的计算后，我收到以下错误：

ArpackNoConvergence：ARPACK 错误 -1：没有收敛（60001 次迭代，0/2 特征向量收敛）

为什么会这样？这是与三对角矩阵的某些属性有关的问题吗？为了以有效的方式对角化我的矩阵，我可以使用哪个 Python（请只使用 Python！）例程？

这是重现我的错误所需的最少代码：

import scipy.sparse.linalg as sl
import numpy as np

dim = 6000
a = np.empty( dim - 1 )
a.fill( 1. )
diag_up = np.diag( a, 1 )
diag_bot = np.diag( a, -1 )

b = np.empty( dim )
b.fill( 1. )

mat = np.diag( b ) + diag_up + diag_bot
v, w = sl.eigsh(mat, 2, which = 'SM')

在我的电脑上，矩阵的构建需要 364 毫秒，而对角化给出了报告的错误。

Answer 1

ARPACK 擅长找到较大的特征值，但很难找到较小的特征值。幸运的是，您可以使用 eigsh 中内置的 shift-invert 选项轻松解决此问题。例如，请参阅here。

import scipy.sparse.linalg as sl
import scipy.sparse as spr
import numpy as np

dim = 6000
diag = np.empty( dim )
diag.fill( 1. )

# construct the matrix in sparse format and cast to CSC which is preferred by the shift-invert algorithm
M = spr.dia_matrix((np.array([diag, diag, diag]), [0,-1, 1]), shape=(dim,dim)).tocsc()

# Set sigma=0 to find eigenvalues closest to zero, i.e. those with smallest magnitude. 
# Note: under shift-invert the small magnitued eigenvalues in the original problem become the large magnitue eigenvalue
# so 'which' parameter needs to be 'LM'
v, w = sl.eigsh(M, 2, sigma=0, which='LM')
print(v)

对于这个特定的示例问题，您可以验证上述是否找到了正确的特征值，因为特征值恰好有一个 explicit formula：

from math import sqrt, cos, pi
eigs = [abs(1-2*cos(i*pi/(1+dim))) for i in range(1, dim+1)]
print(sorted(eigs)[0:2])

Answer 2

如果您需要多个特征值，则使用eigsh 可能不是一个好主意。相反，请考虑使用eigh_tridiagonal，它使用适当的 LAPACK 例程（而不是 ARPACK）。请注意，使用eigh_tridiagonal，我认为您无法计算最小幅度特征值。您可以计算所有特征值、与特定索引对应的特征值或特定范围内的特征值（请参阅documentation page）。我认为最后一个选项可能对您目前的目的有用。

一些示例代码如下：

import numpy as np
import scipy.linalg as la
from timeit import default_timer as timer

dim = 6000
diag = np.ones(dim)
offdiag = np.ones(dim - 1)

eig_range = (-0.005, 0.005) # center it on zero for the comparison with eigsh to work properly!

lapack_time = timer()
evals_lapack, evecs_lapack = la.eigh_tridiagonal(diag, offdiag, select = 'v', select_range = eig_range)
lapack_time = timer() - lapack_time

然后您可以与稀疏 ARPACK 方法进行比较：

import scipy.sparse as spa
import scipy.sparse.linalg as sla

n_eigvals = np.size(evals_lapack)
evals_true = 1.0 - 2.0 * np.cos(np.arange(1, dim + 1) * np.pi / (dim + 1))
sm_idxs = np.argpartition(np.abs(evals_true), (0, n_eigvals))[:n_eigvals]
evals_true_sm = np.sort(evals_true[sm_idxs])

arpack_time = timer()
sparse_csc_matrix = spa.dia_matrix((np.array([diag, diag, diag]), [0, -1, 1]), shape = (dim, dim)).tocsc()
evals_arpack, evecs_arpack = sla.eigsh(sparse_csc_matrix, n_eigvals, sigma = 0, which = 'LM')
arpack_time = timer() - arpack_time

lapack_accuracy = np.sqrt(np.mean(np.square((evals_lapack - evals_true_sm) / evals_true_sm)))
arpack_accuracy = np.sqrt(np.mean(np.square((evals_arpack - evals_true_sm) / evals_true_sm)))

print("for", n_eigvals, "eigenvalues...")
print("are results the same?", np.allclose(evals_lapack, evals_arpack))
print("lapack time: ", lapack_time)
print("arpack time: ", arpack_time)
print("lapack error:", lapack_accuracy)
print("arpack error:", arpack_accuracy)

一些结果：

for 11 eigenvalues...
are results the same? True
lapack time:  0.027329199999996945
arpack time:  0.1122953999999936
lapack error: 1.0656446384601465e-13
arpack error: 3.0175238399729594e-13

for 1137 eigenvalues...
are results the same? True
lapack time:  10.724007
arpack time:  62.3984447
lapack error: 4.6350742410143475e-14
arpack error: 3.0063051257502015e-14