8
$\begingroup$

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?

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ 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? $\endgroup$
    – Federico Poloni
    Aug 25, 2012 at 13:34
  • 2
    $\begingroup$ You might consider also asking this on scicomp.stackexchange.com $\endgroup$
    – Andrew T. Barker
    Aug 25, 2012 at 15:26
  • $\begingroup$ @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. $\endgroup$
    – Jeff
    Aug 26, 2012 at 23:29

2 Answers 2

Reset to default
4
$\begingroup$

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

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – Jeff
    Aug 24, 2012 at 23:52
  • 2
    $\begingroup$ 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. $\endgroup$
    – Igor Rivin
    Aug 25, 2012 at 1:52
  • $\begingroup$ @Jeff have you tried this? I am having the same issue. Is it better Lanczos? It also seems that the Gershgorin bound would be a quickly-computable always-valid value for $c$. $\endgroup$
    – VF1
    Dec 20, 2016 at 18:49
  • 1
    $\begingroup$ Yes, adding $cI$ and doing inverse iterations works well. $\endgroup$
    – Jeff
    Dec 21, 2016 at 21:39
3
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Federico Poloni
    Aug 25, 2012 at 13:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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