6
$\begingroup$

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

|cite|improve this question
$\endgroup$
  • $\begingroup$ Do you have $\mathbf A$, or just $\mathbf M$? $\endgroup$ – J. M. is a poor mathematician Aug 15 '12 at 1:22
  • 1
    $\begingroup$ Your $\mathbf B$ is in fact expressible in terms of the Moore-Penrose pseudoinverse, of which much has been written about... $\endgroup$ – J. M. is a poor mathematician Aug 15 '12 at 1:27
  • $\begingroup$ @J.M. - I have the "A" matrix $\endgroup$ – Mark K Cowan Aug 15 '12 at 10:35
  • $\begingroup$ Then chao's solution would be most appropriate to your circumstances. $\endgroup$ – J. M. is a poor mathematician Aug 15 '12 at 11:44
  • $\begingroup$ 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. $\endgroup$ – Mark K Cowan Aug 15 '12 at 12:44
4
$\begingroup$

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.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Just to let you know that Cholesky worked a charm, thank you! $\endgroup$ – Mark K Cowan Jul 3 '13 at 18:42
3
$\begingroup$

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

|cite|improve this answer
$\endgroup$
  • 1
    $\begingroup$ "...the pseudo-inverse has inversed singular values..." - if none of the singular values are zero, of course. $\endgroup$ – J. M. is a poor mathematician Aug 15 '12 at 2:37
  • $\begingroup$ @J.M. you are right. Answer edited. $\endgroup$ – chaohuang Aug 15 '12 at 3:08
  • $\begingroup$ 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! $\endgroup$ – Mark K Cowan Aug 15 '12 at 10:44
1
$\begingroup$

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

|cite|improve this answer
$\endgroup$
  • 1
    $\begingroup$ Please add some english prose to your answer and consider formatting the mathematics in LaTeX. $\endgroup$ – Dan Rust Jul 2 '13 at 18:11
  • $\begingroup$ @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. $\endgroup$ – rschwieb Jul 2 '13 at 18:46
  • 1
    $\begingroup$ 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}$ $\endgroup$ – Mark K Cowan Jul 3 '13 at 18:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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