带有Python的三角形，对称稀疏基质的对角线化

Question

i具有由python代码计算的对称和tridiagonal矩阵，我想将其对角度化。

在我处理N = 6000的具体情况下，矩阵可以变大。由于它很稀疏，因此我假设对角度化的最佳方法是使用算法scipy.sparse.linalg.eigsh()，尽管我与其他稀疏和对称的矩阵（不是Tridiagonal的矩阵）相对出色。特别是，由于我只需要频谱的低位置部分，因此我在功能中指定k=2和which='SM'。

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

arpacknoconervergence：arpack错误-1：无收敛（60001迭代，0/2 eigenVectors contresed）

为什么会发生这种情况？它是否与Tridiagonal矩阵的某些特性有关？我可以使用哪个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')

在我的PC上，矩阵的构造需要364ms，而对角线化给出了报告的错误。

Answer 1

arpack擅长找到大型特征值，但可能很难找到小型特征。幸运的是，您可以使用eigsh中内置的Shift Envert选项很容易地解决此问题。例如，请参阅此处。

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)

对于此特定示例问题，您可以验证以上是找到正确的特征值，因为特征值恰好具有明确的公式：

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可能不是一个好主意。相反，请考虑使用使用适当的Lapack例程（而不是Arpack）的eigh_tridiagonal。请注意，使用eigh_tridiagonal，我认为您不能计算最小的特征值。您可以计算与特定索引相对应的所有特征值，特征值或特定范围内的特征值（请参见文档页面）。我认为最后一个选项可能对您的当前目的有用。

某些示例代码如下：

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