Pseudoinverse and Least Squares

In [1]: 
import numpy as np
import numpy.linalg as la

np.set_printoptions(precision=4, linewidth=100)
In [2]: 
A = np.random.randn(5, 3)

Now compute the SVD of $A$. Note that numpy.linalg.svd returns $V^T$:

In [3]: 
U, singval, VT = la.svd(A)
V = VT.T

Let's first understand the shapes of these arrays:

In [4]: 
print(U.shape)
print(singval.shape)
print(V.shape) 

(5, 5)
(3,)
(3, 3)

Check the orthogonality of $U$ and $V$:

In [5]: 
U.T.dot(U)
Out[5]:
array([[  1.0000e+00,   0.0000e+00,  -6.9389e-17,   5.5511e-17,   5.5511e-17],
       [  0.0000e+00,   1.0000e+00,  -1.9429e-16,  -6.9389e-17,  -5.5511e-17],
       [ -6.9389e-17,  -1.9429e-16,   1.0000e+00,  -1.6653e-16,   5.5511e-17],
       [  5.5511e-17,  -6.9389e-17,  -1.6653e-16,   1.0000e+00,   0.0000e+00],
       [  5.5511e-17,  -5.5511e-17,   5.5511e-17,   0.0000e+00,   1.0000e+00]])
In [6]: 
V.T.dot(V)
Out[6]:
array([[  1.0000e+00,   3.8858e-16,   5.5511e-17],
       [  3.8858e-16,   1.0000e+00,   1.9429e-16],
       [  5.5511e-17,   1.9429e-16,   1.0000e+00]])

Now build the matrix $\Sigma$:

In [7]: 
Sigma = np.zeros(A.shape)
Sigma[:3, :3] = np.diag(singval)
Sigma
Out[7]:
array([[ 2.4741,  0.    ,  0.    ],
       [ 0.    ,  1.3509,  0.    ],
       [ 0.    ,  0.    ,  0.4768],
       [ 0.    ,  0.    ,  0.    ],
       [ 0.    ,  0.    ,  0.    ]])

Now piece $A$ back together from the factorization:

In [8]: 
U.dot(Sigma).dot(V.T) - A
Out[8]:
array([[ -1.6653e-16,  -2.0817e-16,   0.0000e+00],
       [ -2.2204e-16,   1.2837e-16,   3.1919e-16],
       [ -8.3267e-17,  -2.2204e-16,  -2.4980e-16],
       [  6.6613e-16,  -3.0531e-16,  -4.4409e-16],
       [ -2.2204e-16,   1.3878e-16,   1.6653e-16]])

Next, compute the pseudoinverse:

In [9]: 
SigmaInv = np.zeros((3,5))
SigmaInv[:3, :3] = np.diag(1/singval)
SigmaInv
Out[9]:
array([[ 0.4042,  0.    ,  0.    ,  0.    ,  0.    ],
       [ 0.    ,  0.7402,  0.    ,  0.    ,  0.    ],
       [ 0.    ,  0.    ,  2.0971,  0.    ,  0.    ]])
In [10]: 
A_pinv = V.dot(SigmaInv).dot(U.T)

Now use the pseudoinverse to "solve" $Ax=b$ for our tall-and-skinny $A$:

In [11]: 
b = np.random.randn(5)
In [12]: 
A_pinv.dot(b)
Out[12]:
array([ 0.2425, -3.5347,  0.9537])

Compare with the solution from QR-based Least Squares:

In [13]: 
Q, R = la.qr(A, "complete")
la.solve(R[:3], Q.T.dot(b)[:3])
Out[13]:
array([ 0.2425, -3.5347,  0.9537])