# 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 posibive-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. –  Fabian Sep 5 '12 at 16:06
–  user2468 Sep 5 '12 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]. –  David Speyer Sep 5 '12 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. –  Mark K Cowan Sep 5 '12 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! –  Mark K Cowan Sep 5 '12 at 22:55