Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.$$

share|improve this question
    
Do you have $\mathbf A$, or just $\mathbf M$? –  J. M. Aug 15 '12 at 1:22
1  
Your $\mathbf B$ is in fact expressible in terms of the Moore-Penrose pseudoinverse, of which much has been written about... –  J. M. 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. –  J. M. 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
show 2 more comments

3 Answers

up vote 3 down vote accepted

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.

share|improve this answer
    
Just to let you know that Cholesky worked a charm, thank you! –  Mark K Cowan Jul 3 '13 at 18:42
add comment

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$.

share|improve this answer
1  
"...the pseudo-inverse has inversed singular values..." - if none of the singular values are zero, of course. –  J. M. 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
add comment

$ 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} $

share|improve this answer
1  
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
1  
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.