Pseudo-inverse of A

Here's what happens if some of the singular values are zero.

Suppose A = (v v v ) (s ) (u ) 1 2 3 1 0 0 1 ... 0 2 0 2 u 0 0 0 3 where {u_1, u_2, u_3} and {v_1, v_2, v_3} are orthonormal bases for ^3 and s_1 and s_2 are non-zero.

Carry out the three steps of solving A x = y using a matrix M_i to carry out each step.

Step 1:  Find v_i coords of y .  
Since {v_1, v_2, v_3} is an orthonormal basis the v_i coords are c_i = y · v_i .
But since s_3 = 0 you aren't interested in c_3.
Choose M_1 = (v ) 1 v 2.  Then M_1y   = (v ) y = (v · y ) = (c ) 1 1 1 v v · y c 2 2 2.

Step 2:  Divide the numbers (c ) 1 c 2 by the singular values.  
Choose M_2 = (1/s ) 1 0 1/s 0 2.  
Then M_2(c ) = (1/s ) (c ) = (c /s ) 1 1 0 1 1 1 c 1/s c c /s 2 0 2 2 2 2.

Step 3:  Use the entries of (c /s ) 1 1 c /s 2 2 as the u_i coords of x.
Choose M_3 = (u u ) 1 2.
Then M_3(c /s ) = (u u ) (c /s ) = c_1/s_1u_1 + c_2/s_2u_2 1 1 1 2 1 1 c /s c /s 2 2 2 2.

Combining all three steps, a least squares solution to A x = y  is
    x = Underscript[Underscript[(u u ) (1/s ) (v ), ︸], A^+] y ... 1/s v 0 2 2.

The matrix product A^+= (u u ) (1/s ) (v ) 1 2 1 0 1 1/s v 0 2 2 is the pseudo-inverse of A.
You get the pseudo-inverse by keeping the parts of A^(-1) that correspond to non-zero singular values.  

Computing A^+ by hand is a hard job so you usually use a machine.  In Mathematica you use PseudoInverse[A] and the only hard thing is remembering how to spell pseudo.

One big advantage of using A^+ to solve A x = y is that you don't need to know in advance whether A x = y has solutions.  
Without any worries you compute x = A^+y .  
If A x = y has solutions, x = A^+y will be one of them.
If A x = y has NO solutions, then x = A^+y will give you a least squares solution.

(Q44) Suppose A = (v v v ) (5 0 0 0) (u ) 1 2 3 1 ... 3 u 4.
Express A^+ as a product of three matrices.

Up

Copyright © 2007 Todd Will
Created by Mathematica  (April 15, 2007)