Inverse of symmetric matrix $M = A A^\top$

I have a matrix, generated by the product of a non-square matrix with its own transpose:

$$M = A A^\top.$$

I need the inverse of $M$, assuming $\det(M) \neq 0$.

Given the nature of the matrix $M$, are there any specialised computational methods to find its inverse, prioritising precision over speed? Gauss-Jordan is prone to error, and I hope to find something nicer than and with comparable precision to the adj/det method.

$M$ will either have be around $70 \times 70$, or $1000 \times 1000$ in size.

I've had a quick read of the Matrix Cookbook and of this page, but (at the present time of 1 am) I'm struggling to see how it could help me.

In case it helps, I'm actually trying to calculate:

$$B = (A A^\top)^{-1} A.$$

-
Do you have $\mathbf A$, or just $\mathbf M$? –  Guess who it is. Aug 15 '12 at 1:22
Your $\mathbf B$ is in fact expressible in terms of the Moore-Penrose pseudoinverse, of which much has been written about... –  Guess who it is. Aug 15 '12 at 1:27
@J.M. - I have the "A" matrix –  Mark K Cowan Aug 15 '12 at 10:35
Then chao's solution would be most appropriate to your circumstances. –  Guess who it is. Aug 15 '12 at 11:44
Thanks. I've already started implementing Cholesky (which vectorises easily with SSE), I'll give SVD a go later though since I'll need to implement it anyway soon in another project. –  Mark K Cowan Aug 15 '12 at 12:44

Cholesky decomposition! It's a very stable method and it works on sparse matrices too.

But if you are just trying to solve the normal equation, I may suggest conjugate gradient or SOR.

-
Just to let you know that Cholesky worked a charm, thank you! –  Mark K Cowan Jul 3 '13 at 18:42

Actually, you don't need to calculate $(A A^T)^{-1}$ to compute B. What you are trying to calculate is the left inverse of $A^T$, and it is give by

$$B = (AA^T)^{-1} A = {A^{-T}}_L = ((A^T)^T A^T)^{-1} (A^T)^T$$

Since in your case the left inverse exists, it's equivalent to the pseudo-inverse of $A^T$, which can be computed by SVD. And since the left (right) singular vectors of a matrix are right (left) singular vectors of its transpose or left/pseudo inverse, and a matrix and its transpose have the same singular values, but the pseudo-inverse has inversed singular values (non-zeros ones), the pseudo-inverse of $A^T$ has the same singular vectors with A, and has the inversed singular values of $A$.

So, all you have to do is calculating the SVD of $A$, and inversing it's non-zero singular values, then you will get $B$.

-
"...the pseudo-inverse has inversed singular values..." - if none of the singular values are zero, of course. –  Guess who it is. Aug 15 '12 at 2:37
@J.M. you are right. Answer edited. –  chaohuang Aug 15 '12 at 3:08
Thanks - I've not done much linear algebra since I started university (6 years ago...) so it'll take me a while to research and digest your suggestions - hopefully by the time I've tried and compared them I'll have enough rep to give your answers the appropriate votes! –  Mark K Cowan Aug 15 '12 at 10:44

$B = \left(A A^\top\right)^{-1}A = \left(A A^\top\right)^{-1}A A^\top \left(A^\top\right)^{-1} = \left(A A^\top\right)^{-1}A A^\top \left(A^\top\right)^{-1} = I \left(A^\top\right)^{-1} = A^{-\top}$

-
Please add some english prose to your answer and consider formatting the mathematics in LaTeX. –  Daniel Rust Jul 2 '13 at 18:11
@Xiang you can wrap the math in dollar signs. The carats do work for powers, but if you have anything more complex than a single character in a power, you need to gather it all together with {} so that it will all become a superscript. –  rschwieb Jul 2 '13 at 18:46
This only applies if $A$ is a square matrix, which I explicitly stated that it is not. Also, a tidier proof (for a square matrix $A$) would be: $B = \left(A A^\top\right)^{-1}A = A^{-\top} A^{-1} A = A^{-\top}$ –  Mark K Cowan Jul 3 '13 at 18:39