三重対称で疎な行列をPythonで対角化する

Question

私はPythonコードで計算されたNxN対称および三重対角行列を持っています。対角化したいのです。

具体的にはN = 6000を扱っていますが、マトリックスが大きくなることがあります。それは疎であるので、対角化する最良の方法は、アルゴリズムscipy.sparse.linalg.eigsh()を使用することであると仮定しました。これは、他の疎な対称行列（三重対角行列ではありませんが）でもうまくいきました。特に、スペクトルの低い部分だけが必要なので、私はk=2とwhich='SM'を関数に指定しています。

ただし、この場合には、このアルゴリズムが動作していないようです、計算の約20分後、私は次のエラーを取得するため、：

ArpackNoConvergence: ARPACK error -1: No convergence (60001 iterations, 0/2 eigenvectors converged)

なぜこの出来事はありますか？三重対角行列のいくつかの性質に関連する問題ですか？効率的な方法で行列を対角化するために、どのPython（そしてPythonのみ！）ルーチンを使用できますか？対角化が報告されたエラーを与えながら、行列の建設が364msを取る私のPC上で

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') 

：

は、ここに私のエラーを再現するために要求された最小限のコードです。

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