Menu
I'd like to find the N smallest eigenvalues of a sparse matrix in Python. I've tried using the scipy.sparse.linalg.eigen.arpack package, but it is very slow at computing the smallest eigenvalues. I read somewhere that there is a shift-invert mode, but when I try using it, I receive an error message telling me that the shift-invert mode is not yet supported. Any ideas as to how I should proceed?
421
scipy/sparse/linalg/eigsh can output the k smallest(largest) eigenvalues and eigenvectors; scipy/linalg/eigh also provides the option to select subset of eigenvalues; numpy/linalg/eigvalsh outputs all the eigenvalues.
If you know that A is symmetric positive-definite, then the spectral shift B=A−λmaxI will work. Use the power method on B, then add λmax to the result to get the smallest eigenvalue of A. The reason this shift works is that a positive-definite matrix has all positive eigenvalues.
ARPACK is a Fortran package which provides routines for quickly finding a few eigenvalues/eigenvectors of large sparse matrices. In order to find these solutions, it requires only left-multiplication by the matrix in question. This operation is performed through a reverse-communication interface.
Sparse Matrices in Python SciPy provides tools for creating sparse matrices using multiple data structures, as well as tools for converting a dense matrix to a sparse matrix. Many linear algebra NumPy and SciPy functions that operate on NumPy arrays can transparently operate on SciPy sparse arrays.
Comparing the documentation of scipy.sparse.linalg.eigs from SciPy v0.9 with the documentation of scipy.sparse.linalg.eigs from SciPy v0.10 it appears that the shift-invert mode is implemented and working since v0.10. Specifically, the explanation of the sigma parameter in the v0.9 documentation states it is not implemented, but the v0.10 documentation does not indicate that is the case.
If you do not have SciPy v0.10, or later, installing the latest should enable you to utilize shift-invert mode with the sparse eigensolver.
As mentioned in the question, it is possible to use the ARPACK interface to find small-magnitude eigenvalues. This is done by passing which='SM' when calling scipy.sparse.linalg.eigs. It is, however, as stated in the question, slow. This is confirmed in the SciPy Tutorial's section on Sparse Eigenvalue Problems with ARPACK, where it states:
Note that ARPACK is generally better at finding extremal eigenvalues: that is, eigenvalues with large magnitudes. In particular, using
which = 'SM'may lead to slow execution time and/or anomalous results. A better approach is to use shift-invert mode.
Let's look at some code trying to use shift-invert with both v0.9 and v0.10 of SciPy. In both cases, we will use the following code.
from scipy.sparse import identity
from scipy.sparse.linalg import eigs
A = identity(10, format='csc')
A.setdiag(range(1, 11))
eigs(A, 3, sigma=0) # find three eigenvalues near zero using shift-invert mode
Running the code in SciPy v0.9 results in an exception being raised.
NotImplementedError: shifted eigenproblem not supported yet
Running the code in SciPy 0.10 produces expected results.
(array([ 1.+0.j, 2.+0.j, 3.+0.j]),
array([[ -1.00000000e+00+0.j, 5.96300068e-17+0.j, 9.95488924e-17+0.j],
[ 3.55591776e-17+0.j, 1.00000000e+00+0.j, -4.88997616e-16+0.j],
[ -3.79110898e-17+0.j, 1.16635626e-16+0.j, 1.00000000e+00+0.j],
[ -1.08397454e-17+0.j, 1.23544164e-17+0.j, 1.78854096e-15+0.j],
[ 1.68486368e-17+0.j, -9.37965967e-18+0.j, 2.05571432e-16+0.j],
[ -2.97859557e-19+0.j, -3.43100887e-18+0.j, 3.35947574e-17+0.j],
[ 1.89565432e-17+0.j, -3.61479402e-17+0.j, -1.33021453e-17+0.j],
[ -1.40925577e-18+0.j, 3.16953070e-18+0.j, 7.91193025e-17+0.j],
[ 6.76947854e-19+0.j, -3.75674631e-19+0.j, 3.61821551e-17+0.j],
[ -3.07505146e-17+0.j, -6.52050102e-17+0.j, -8.57423599e-16+0.j]]))
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
