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

$\endgroup$
7
  • $\begingroup$ Do you have $\mathbf A$, or just $\mathbf M$? $\endgroup$
    – J. M. ain't a mathematician
    Aug 15, 2012 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. ain't a mathematician
    Aug 15, 2012 at 1:27
  • $\begingroup$ @J.M. - I have the "A" matrix $\endgroup$
    – Mark K Cowan
    Aug 15, 2012 at 10:35
  • $\begingroup$ Then chao's solution would be most appropriate to your circumstances. $\endgroup$
    – J. M. ain't a mathematician
    Aug 15, 2012 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, 2012 at 12:44

3 Answers 3

Reset to default
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.

$\endgroup$
1
  • $\begingroup$ Just to let you know that Cholesky worked a charm, thank you! $\endgroup$
    – Mark K Cowan
    Jul 3, 2013 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$.

$\endgroup$
3
  • 1
    $\begingroup$ "...the pseudo-inverse has inversed singular values..." - if none of the singular values are zero, of course. $\endgroup$
    – J. M. ain't a mathematician
    Aug 15, 2012 at 2:37
  • $\begingroup$ @J.M. you are right. Answer edited. $\endgroup$
    – chaohuang
    Aug 15, 2012 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, 2012 at 10:44
0
$\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} $

$\endgroup$
3
  • 1
    $\begingroup$ Please add some english prose to your answer and consider formatting the mathematics in LaTeX. $\endgroup$
    – Dan Rust
    Jul 2, 2013 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, 2013 at 18:46
  • 2
    $\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, 2013 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.