Complex eigenvalues computation using scipy.sparse.linalg.eigs

Question

Given the following input numpy 2d-array A that may be retrieved with the 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.

First of all, a little bit of context. This matrix A results from a quadratic eigenvalue problem of size N which has been linearized in an equivalent eigenvalue problem of double size 2*N. A has the following structure (blue color being zeroes):

plt.imshow(np.where(A > 1e-15,1.,0), interpolation='None')

and the following features:

A shape = (748, 748)
A dtype = float64
A sparsity ratio = 77.64841716949297 %

The true dimensions of A are much bigger than this small reproducible example. I expect the real sparsity ratio and shape to be close to 95% and (5508, 5508) for this case.

The resulting eigenvalues of A are complex (which come in complex conjugate pairs) and I am more interested in the ones with the smallest imaginary part in modulus.

Problem: when using direct solver:

w_dense = np.linalg.eigvals(A) 
idx = np.argsort(abs(w_dense.imag))
w_dense = w_dense[idx]

calculation times become rapidly prohibitive. I am thus looking to use a sparse algorithm:

from scipy.sparse import csc_matrix, linalg as sla
w_sparse = sla.eigs(A, k=100, sigma=0+0j, which='SI', return_eigenvectors=False)

but it seems that ARPACK doesn't find any eigenvalues this way. From the scipy/arpack tutorial, when looking for small eigenvalues like which = 'SI', one should use the so-called shift-invert mode by specifying sigma kwarg, i.e. in order for the algorithm to know where it could expect to find these eigenvalues. Nonetheless, all of my attempts did not yield any results...

Could someone more experienced with this function give me a hand in order to make this work?

Here follows a whole code snippet:

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)

Answer 1

Problem finally solved, I guess I should have read the documentation more carefully, but yet, the following is quite counter-intuitive and could be better emphasized in my opinion:

... 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 .

This way, when looking for small eigenvalues, in order to help LAPACK do the work, one should activate shift-invert mode by specifying an appropriate sigma value while also reversing the desired specified subset specified in the which keyword argument.

Thus, it is simply a matter of executing:

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]

Therefore, I can only hope this mistake help someone else :)