# Finding the smallest eigenvalues of a large, but structured, matrix

I'm trying to find the eigenvector corresponding to the second smallest eigenvalue of a large $(4,000,000 \times 4,000,000)$ matrix $M$. $M$ is a Laplacian matrix, and it has the following structure: $M = D - AA'$, where $D$ is a diagonal matrix and $A$ is a large sparse matrix. $A$ has dimension $4,000,000 \times 10,000,000$, but only about $40,000,000$ non-zero entries. So I can rapidly perform matrix-vector multiplication: $Mv = Dv - A(A'v)$.

Currently I'm using Scipy (which calls ARPACK) to find the smallest eigenvalues and corresponding eigenvectors of $M$. The implementation is a variant of Lanczos method, as far as I can tell. Unfortunately the implementation fails to converge with some frequency. Even when it does converge, since I need fairly high accuracy, it takes many iterations to converge. Any suggestions?

Lanczos method is much more reliable when finding the largest eigenvalues and the corresponding eigenvectors, as I understand it. So I was thinking that I could transform the problem, perhaps instead finding the largest eigenvalues of M's left inverse. But, even if I somehow manage to find a good approximation of M's left inverse (most methods don't converge), the largest eigenvalue is infinite (the smallest eigenvalue of $M$ -- a Laplacian matrix -- is $0$), so it seems there would still be convergence issues, since I can't solve for the second largest eigenvalue without also solving for the largest, as I understand it. Is there an easier way?

-
Arpack should already use shift-and-invert if you set it to find the smallest eigenvalues of a matrix. Are you sure you are not trying to reinvent what is already implemented there? – Federico Poloni Aug 25 '12 at 13:34
You might consider also asking this on scicomp.stackexchange.com – Andrew T. Barker Aug 25 '12 at 15:26
@Federico -- It doesn't use shift-and-invert by default, at least when accessed through Scipy, but you're basically right -- once sigma is specified, shift-and-invert is used. However, with the default settings, I was encountering an error messages ("Error in inverting [A-sigma*M]: function gmres did not converge info=40"), even when run with a tiny (and invertible) 4 x 4 matrix. So I asked this question to figure out how to do it manually, though with the answers below, I'm starting to see how to work around the error, by overriding the defaults and substituting conjugate gradient for gmres. – Jeff Aug 26 '12 at 23:29

Well, if you add $c I$ to your matrix, for some reasonable value of $c,$ it will become nonsingular. As for the left inverse, since your matrix is sparse, to compute the backward iteration you can use the conjugate gradient method, which will be very fast (there are fancier "Krylov space" methods, but you need not go there, unless conjugate gradient fails).
Thanks for the response. So you're proposing I 1) set $M=cI + D - AA'$, 2) solve for $M^{-1}$ with conjugate gradient, and 3) find the eigenvector corresponding to the second largest eigenvalue of $M^{-1}$? I like the use of ridging -- I had forgotten about that technique. But, one concern: Since $M$ isn't sparse, won't conjugate gradient produce a representation of $M^{-1}$ that is too large to fit into memory? – Jeff Aug 24 '12 at 23:52
You don't need a representation of $M^{-1},$ you just need to iterate $M^{-1}$ on some random starting vector, which is what you do with conjugate gradient, so this requires no storage. – Igor Rivin Aug 25 '12 at 1:52
In other words, use ARPACK on $M-\alpha ww^T$, where $w$ is the eigenvector corresponding to zero and $\alpha$ is a constant. This "moves away" the zero eigenvalue, but leaves every other one in place. – Federico Poloni Aug 25 '12 at 13:37