Diagonalization of a tridiagonal, symmetric sparse matrix with Python

user_2337272 0
Viewed 1471

I have an NxN symmetric and tridiagonal matrix computed by a Python code and I want to diagonalize it.

In the specific case I'm dealing with N = 6000, but the matrix can become larger. Since it is sparse, I assumed the best way to diagonalize it was to use the algorithm scipy.sparse.linalg.eigsh(), which performed extremely good with other sparse and symmetric matrices (not tridiagonal ones, though) I worked with. In particular, since I need only the low lying part of the spectrum, I'm specifying k=2 and which='SM' in the function.

However, in this case this algorithm seems not to work, since after approximately 20 minutes of computation I get the following error:

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

Why is this happening? Is it a problem related to some properties of tridiagonal matrices? Which Python (and please, only Python!) routine can I use in order to diagonalize my matrix in an efficient way?

Here's the requested minimal code to reproduce my error:

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

On my pc the construction of the matrix takes 364ms, while the diagonalization gives the reported error.

scipy linear-algebra matrix numpy python
2 Answers

ARPACK is good at finding the large-magnitude eigenvalues but can struggle to find the small ones. Fortunately, you can work around this quite easily by using the shift-invert options built into eigsh. See, for example, 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)

For this particular example problem, you can verify that the above is finding the correct eigenvalues since the eigenvalues happen to have an 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])
user_19711870

Using eigsh may not be a good idea with regard to performance if you want more than a few eigenvalues. Instead, consider using eigh_tridiagonal, which uses the appropriate LAPACK routine (instead of ARPACK). Note that with eigh_tridiagonal, I do not think you can compute the smallest magnitude eigenvalues. You can compute all eigenvalues, eigenvalues corresponding to particular indexes, or eigenvalues within a particular range (see the documentation page). I think the last option may be useful for your present purposes though.

Some example code is below:

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

You can then make a comparison with the sparse ARPACK approach:

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)

Some results:

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
user_80798700
Related