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)