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!
 A: 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

A: 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.
