[Math] Eigenvalues of $A^T A$

eigenvalues-eigenvectorsmatrices

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 1024×1024 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!

Best Answer

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
Related Question