14

I've got a program that computes the largest eigenvalue of many real symmetric 50x50 matrices by performing singular-value decompositions on all of them. The SVD is a bottleneck in the program.

Are there algorithms that are much faster in finding the largest eigenvalue, or would optimizing this part not give much return on investment?

David Ketcheson
  • 16,522
  • 4
  • 54
  • 105
Anna
  • 243
  • 1
  • 4
  • Could you give some more information on your matrices, e.g. if anything is known about their structure, the range of their eigenvalues or their similarity to each other? – Pedro Jun 18 '12 at 15:26
  • It's a covariance matrix ($XX^T$). Testing shows that all but the 5 or so largest eigenvalues are close to zero, and that the largest eigenvalue is at least ~20% larger than the second largest. Since there's a lot of eigenvalues close to zero, I suppose the range isn't important? It could be rescaled to any range. The scale I'm using currently gives me a range of 150~200. – Anna Jun 18 '12 at 15:31
  • Also, the matrix is not very closely singular, so the SVD problem is well-conditioned. – Anna Jun 18 '12 at 15:33
  • Since $XX^T$ is symmetric and positive (semi) definite you could use the Cholesky factorization instead of the SVD. The Cholesky factorization takes a lot fewer flops to compute than the SVD but being an exact method still takes $O(n^3)$ flops. – Ken Jun 19 '12 at 05:44
  • @Anna: Have you tried out any of the many approaches proposed here? I'd be quite curious to know what worked best in practice for you... – Pedro Jun 25 '12 at 16:09
  • @Pedro: I'm not done with it yet. When I am, I'll be sure to post the results for my specific case and accept an answer. Sorry about the delay. – Anna Jun 26 '12 at 09:04
  • @Anna I am working on a similar problem. Have you found anything that you can post as your own answer? – adam W Mar 12 '13 at 18:32

3 Answers3

13

Depending on the precision you require for the largest eigenvalue, you could try using the Power Iteration.

For your specific example, I would go as far as to not form $A=XX^\mathsf{T}$ explicitly, but compute $x \leftarrow X(X^\mathsf{T}x)$ in each iteration. Computing $A$ would require $\mathcal O(n^3)$ operations whereas the matrix-vector product requires only $\mathcal O(n^2)$.

The convergence rate depends on the separation between the largest two eigenvalues, so this may not be a good solution in all cases,

Pedro
  • 9,573
  • 1
  • 36
  • 45
  • 1
    If the largest eigenvalue is 20% larger than the next, the power iteration should converge pretty rapidly (all other eigenvalues get damped by a factor of 5/6 in each iteration, so you get one digit for every 13 iterations. – Wolfgang Bangerth Jun 21 '12 at 06:45
  • 2
    Krylov subspace methods are strictly better than power methods, as they contain the vector from power iteration with the same number of iterations. – Jack Poulson Jun 25 '12 at 13:08
  • 1
    @JackPoulson: Yes, but each iteration is more expensive to compute... Would it really be worth it for such a small problem? – Pedro Jun 25 '12 at 13:14
  • @Pedro: of course, the matvecs require quadratic work and the Rayleigh quotient eigensolve and subsequent expansion are trivial in comparison. – Jack Poulson Jun 25 '12 at 15:43
  • 1
    Code expense? Since @JackPoulson broached the issue, B. Parlett et al (1982) ("On Estimating the Largest Eigenvalue with the Lanczos Algorithm") compare power method, power method + Aitken acceleration, and an application of Lanczos targetting the largest eigenvalue of a real symmetric (or Hermitian) pos. def. matrix. They conclude the Lanczos method is more efficient if even modest accuracy (of first eigenvalue relative to second) is needed, and better at avoiding misconvergence. – hardmath Jun 27 '12 at 12:02
5

If only 5 eigenvalues are very significant, the Lanczsos algorithm with $X(X^Tx)$ as matrix-vector multiply should give fast linear convergence after 5 initial steps, hence a fairly accurate largest eigenvalue with few iterations.

Arnold Neumaier
  • 11,318
  • 20
  • 47
  • Are you ( @ArnoldNeumaier ) thinking of something like this, suitably simplified ($B = T = I$)? It's interesting that it gives a different approximation from Lanczos if a third vector is retained, over the same Krylov subspace. – hardmath Jun 27 '12 at 12:24
  • No; I meant the standard Lanczsos algorithm but had in a hurry written CG. Now corrected. – Arnold Neumaier Jun 27 '12 at 14:07
5

For a positive semi-definite matrix such as $A = XX^T$ it may be worth the effort to accelerate convergence with a spectrum shift. That is, a suitable scalar $\mu$ is chosen and the power method is applied to $A - \mu I$ instead of $A$.

A few iterations of the basic power method should give you a rough estimate $||Ax||/||x||$ of the largest eigenvalue $\lambda_1$. Assuming the dominant eigenvalue has multiplicity 1 and that all the others are in $[0,\frac{5}{6} \lambda_1]$, then $A - \frac{5}{12} \lambda_1 I$ would have a largest eigenvalue $\frac{7}{12} \lambda_1$ and the rest in $[\frac{-5}{12} \lambda_1, \frac{5}{12} \lambda_1]$.

In other words you would increase the dominance of the largest eigenvalue from 20% over the next largest to 40% over the next largest (absolute value of an) eigenvalue. Geometric convergence of the power method would accelerate accordingly. Once the largest eigenvalue of $A - \mu I$ is found to sufficient accuracy, $\lambda_1$ is estimated by adding back the shift $\mu$ that had been taken away.

Note that you need not explicitly form $A - \mu I$ because $(A - \mu I)x = X(X^Tx) - \mu x$ can still be computed with $O(n^2)$ effort.

Community
  • 1
hardmath
  • 3,359
  • 2
  • 21
  • 40
  • This would seem to require having a good idea of what the magnitude of the second-largest eigenvalue is. How would you approximate it in such a case? – Pedro Jun 21 '12 at 22:05
  • @Pedro: Applying the shifted power iteration only requires an estimate of the largest eigenvalue $\lambda_1$, but as with the unshifted power method, analysis shows the rate of convergence depends on $|\lambda_2|/|\lambda_1|$ (throwing in abs. values that are unnecessary for the pos. semi-definite case at hand). In turn observed rates of convergence can be used to estimate $|\lambda_2|/|\lambda_1|$, and hence the size of $\lambda_2$ relative to $\lambda_1$ if desired. I was suggesting what benefit you'd see in a case such as Anna describes in her comments below the Question. – hardmath Jun 21 '12 at 22:26