Eigenvalues of $A^T A$

I have a real, positive-definite, symmetric matrix $M$, calculated from a matrix $A$ which is known:

$M = A^\top A$

Given the nature of $M$, is there any specialised (read: fast code and/or simple code) way to find the largest eigenvalue of $M$? It's years since I did linear algebra, so I'm a bit rusty. I've had a look in Numerical Recipes, Wikipedia and the Matrix cookbook.

I found the following in the Matrix cookbook and got excited too early:

For a symmetric, positive-definite matrix $A$:

$eig\left(A A^\top\right) = eig\left(A^\top A\right) = eig\left(A\right) \cdot eig\left(A\right)$

Of course, it is $M$ that is positive-definite and symmetric, not $A$. $A$ isn't even square. M will probably be 1024x1024 or larger. The size of $M$ (i.e. width of $A$) can be constrained by the algorithm if needed, i.e. I don't mind if your suggestions require that it be a multiple of 4 or a power of 2.

[Edit:] I'm mainly looking for something that I can code natively (taking advantage of SSE/AVX where possible), although your MATLAB/Octave suggestions will definitely be useful to aid my understanding of the algorithms!

• SVD of $A$ might help. I don't know however if there are specialized algorithms which calculate only the largest singular values. Sep 5, 2012 at 16:06
• – user2468
Sep 5, 2012 at 16:42
• I expect any mathematical software package will have specialized routines for this, as it is a quite common requirement. For example, in Mathematica, Eigenvalues[Transpose[A].A, 1] or SingularValueList[A, 1]. Sep 5, 2012 at 16:50

A simple iterative algorithm to find the largest eigenvalue of a symmetric matrix $B$ is to compute repeatedly $x_{n+1} = B \, x_n$ starting from a random vector $x_0$, normalizing the length in each iteration. When it converges, you've found an eigenvector, and then you get the eigenvalue. There are some complications (repeated eigenvalues, event that $x_0$ is very near another eigenvector, etc), but that's the idea. Eg in Octave:

>>> A = rand(A,4); B = A'*A;  x = rand(1,4);
>>> x = x * B  /norm(x)
x =   1.7215   2.4537   2.3203   1.5036
>>> x = x * B  /norm(x)
x =   1.9544   3.0235   2.8072   1.6825
>>> x = x * B  /norm(x)
x =   1.9484   3.0362   2.8112   1.6717
>>> x = x * B  /norm(x)
x =   1.9477   3.0372   2.8113   1.6705
>>> x = x * B  /norm(x)
x =   1.9476   3.0373   2.8113   1.6704
>>> x = x * B  /norm(x)
x =   1.9476   3.0373   2.8113   1.6704
>>> norm(x * B) / norm(x)
ans =  4.8695
>>> eig(B)
ans =  0.074894   0.130166   0.485889   4.869534

• Recalling the classroom of eigenvalues+vectors as measurements of stretching and the corresponding directions, your method makes loads of sense (in addition to a nice simplicity about it). It will also play nicely with SSE or GPU optimisation, I'll profile it against Ed Gorcenski's suggestion. Sep 5, 2012 at 22:59

There are algorithms to compute the largest eigenvalues; in fact, many eigenvalue algorithms compute the eigenvalues in descending order (according to magnitude).

However, to compute just the single largest eigenvalue, the algorithm in this paper should work, but it is too complex to describe completely here.

• Thanks, I was looking at Lanczos on wikipedia before I decided to ask on here - I'll give that paper a proper read tomorrow! Sep 5, 2012 at 22:55