I have a dense real symmetric square matrix. The dimension is about 1000x1000. I need to compute the first principal component and wonder what the best algorithm to do this might be.

It seems that MATLAB uses the Arnoldi/Lanczos algorithms (for eigs). But from reading about them I'm not sure whether they have any advantages over simple power iteration, since my matrix is not sparse and I'm only interested in the first eigenvector.

Any recommendations what's the fastest algorithm in this case?

  • 1
    $\begingroup$ On my computer, on a randomly-generated 1000 X 1000 symmetric matrix, the "eigen" function in R took about one second to compute all the eigenvalues and vectors, rounding up. Your mileage may vary, but I doubt your algorithm choice makes any difference at timings like that. $\endgroup$
    – jbowman
    Jan 28 '12 at 2:34
  • $\begingroup$ Yes, that's true of course. I'm not really concerned in making my program run faster. I'm just curious whether the mentioned more complicated techniques are also considered superior in this use-case (dense, only first eigenvector), or whether there are different techniques for dense matrices. $\endgroup$
    – Mika Fischer
    Jan 28 '12 at 2:57
  • $\begingroup$ Do you mean the eigenvector corresponding to the largest or smallest eigenvalue? It sounds like you want the former. $\endgroup$ Jan 28 '12 at 15:45
  • $\begingroup$ Yes, the eigenvector corresonding to the eigenvalue with largest magnitude. $\endgroup$ Jan 28 '12 at 16:00

The fastest method will likely depend upon the spectrum and normality of your matrix, but in all cases Krylov algorithms should be strictly better than power iteration. G.W. Stewart has a nice discussion of this issue in Chapter 4, Section 3 of Matrix Algorithms, Volume II: Eigensystems:

The power method is based on the observation that if $A$ has a dominant eigenpair then under mild restrictions on $u$ the vectors $A^k u$ produce increasingly accurate approximations to the dominant eigenvector. However, at each step the power method considers only the single vector $A^k u$, which amounts to throwing away information contained in the previously generated vectors. It turns out that this information is valuable..."

and he goes on to show that, for $100 \times 100$ diagonal matrix with the $i$'th diagonal value set to $.95^i$ (counting from $i=0$), after 25 iterations the Krylov subspace captures the dominant eigenvector eight orders of magnitude better than power iteration.

  • $\begingroup$ Hmm, I'd have thought MRRR was now the standard method when one just wants a few eigenvectors... $\endgroup$
    – J. M.
    May 2 '13 at 14:18
  • $\begingroup$ MRRR runs on symmetric tridiagonal matrices, and so to use it for dense matrices one must first reduce the dense matrix to tridiagonal form through similarity transformation(s), which takes cubic work in general. If only $k$ matrix-vector multiplications are required for a Krylov method, then $O(k n^2 + k^2 n + k^3)$ work is required. If $k$ is very small relative to $n$, Krylov methods should win. $\endgroup$ May 2 '13 at 21:47
  • $\begingroup$ I see; somehow I had the impression that you needed to tridiagonalize first before doing Krylov. Thanks! $\endgroup$
    – J. M.
    May 3 '13 at 0:50
  • $\begingroup$ Lanczos is actually gradually building up said tridiagonal matrix. $\endgroup$ May 3 '13 at 1:20

Power iteration is the simplest, but as mentioned above it would likely converge very slowly if the matrix is very non-normal. You get a "hump" phenomenon where the sequence appears to diverge for many iterations before asymptotic behavior kicks in.

Since your matrix is symmetric you could consider RQI iterations, which in the symmetric case yields cubic convergence: http://en.wikipedia.org/wiki/Rayleigh_quotient_iteration .

What makes Arnoldi or Lanczos iterations very nice (at least in my opinion, but I do not research numerical linear algebra) is that they are very versatile. It is usually possible to control which eigenvalues they give you, and how many you get. This is especially true in the symmetric case (and even better if your matrix is definite). For symmetric problems they are very robust. As a black box they work well, but they are also very receptive to new problem information, such as the ability to solve systems involving the matrix.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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