Take the 2-minute tour ×
MathOverflow is a question and answer site for professional mathematicians. It's 100% free, no registration required.

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?

share|improve this question
1  
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
2  
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
add comment

2 Answers 2

up vote 3 down vote accepted

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

share|improve this answer
    
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
2  
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
add comment

If it is a Laplacian then you not only know the smallest eigenvalue is zero, but you also know its corresponding eigenvector. You can use this information by essentially adding a single step to your procedure, it's something like adding back the mean of v to each element of Mv. Then use arnoldi to find the smallest (instead of second smallest) eigenpair.

share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.