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])