4

Suppose I have a symmetric linear operator $A:\mathbb{R}^k \rightarrow \mathbb{R}^k$ where $k$ is "small" (eg., $k=100$), and I want to find it's first few eigenvectors, (eg., $10$ eigenvectors).

If we had a matrix representation of $A$ this would be a standard problem in dense numerical linear algebra. However, the catch is that we don't have the entries of $A$, only a code that can evaluate the action of $A$ on a vector. Furthermore, evaluating $Av$ is moderately costly - it involves solving a PDE behind the scenes.

What is the best way to find the first few eigenvectors of $A$ in this situation? My first thought was the Lanczos algorithm, but that's usually good for large sparse systems, whereas my system is small and dense.

Nick Alger
  • 3,143
  • 15
  • 25

2 Answers2

6

Is the matrix also positive, such that eigenvalues are equal to singular values? You can estimate $k$ singular values of a general matrix (possibly nonsquare) with several digits of accuracy provided they are at least somewhat isolated in $8k$ matrix multiplies. The algorithm, as described in the excellent survey by Halko, Martinsson, and Tropp (2012) is:

  1. choose a random $n\times 2k$ matrix $X$
  2. compute $Y = (A^* A)^q A X$ where $q=1$ is typically sufficient
  3. orthogonalize $Y = QR$ with $Q$ orthogonal
  4. compute $B = Q^* A$
  5. compute an SVD of the small matrix: $\tilde U \Sigma V^* = B$
  6. If you want the left singular vectors, let $U = Q \tilde U$

The first $k$ singular values and singular vectors typically provide several digits of accuracy. This algorithm is often more accurate than Lanczos or Arnoldi iteration and has fewer synchronization points. The advantage is more significant for large problems and if applying the operator to several vectors at once is significantly cheaper than applying it to one at a time. This could be a good choice if you only want a couple values.

Jed Brown
  • 25,650
  • 3
  • 72
  • 130
  • Thanks Jed! Yes the matrix is positive, but I don't have a fixed bound on the condition number. I actually was thinking along these lines earlier today - take $k$ initial guesses for the eigenvectors (from a previous timestep, for example), append $k$ more random vectors, and then apply $A$ one or two times and orthogonalize like a simultaneous power method. – Nick Alger Jul 14 '12 at 03:09
  • 1
    A large condition number is not "bad" at all. This method is more accurate if the largest eigenvalues are well-isolated.

    Note that to find 10 out of 100 eigenvalues, Arnold Neumaier's method will cost about the same and accurately gives you all the eigenvalues. This method is better if you only want, say, 3 eigenvalues.

     – Jed Brown Jul 14 '12 at 03:40
  • Yeah, it seems these methods are ideally suited for finding finite rank approximations of compact operators. – Nick Alger Sep 08 '12 at 05:07
4

To find with reasonable accuracy the 10 dominant eigenvectors of a general 100 x 100 matrix you can hardly do better than to build the matrix using 100 matrix vector multiplies and then find its eigenvalues.

An exception is when you are confident that the rank is much smaller than $n$. In that case, Lanczsos with full reorthogonalization is the way to go.

Arnold Neumaier
  • 11,318
  • 20
  • 47