# Pseudoinverse and Least Squares

In :
import numpy as np
import numpy.linalg as la

np.set_printoptions(precision=4, linewidth=100)

In :
A = np.random.randn(5, 3)


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

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


Let's first understand the shapes of these arrays:

In :
print(U.shape)
print(singval.shape)
print(V.shape)

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



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

In :
U.T.dot(U)

Out:
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 :
V.T.dot(V)

Out:
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 :
Sigma = np.zeros(A.shape)
Sigma[:3, :3] = np.diag(singval)
Sigma

Out:
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 :
U.dot(Sigma).dot(V.T) - A

Out:
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 :
SigmaInv = np.zeros((3,5))
SigmaInv[:3, :3] = np.diag(1/singval)
SigmaInv

Out:
array([[ 0.4042,  0.    ,  0.    ,  0.    ,  0.    ],
[ 0.    ,  0.7402,  0.    ,  0.    ,  0.    ],
[ 0.    ,  0.    ,  2.0971,  0.    ,  0.    ]])

In :
A_pinv = V.dot(SigmaInv).dot(U.T)


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

In :
b = np.random.randn(5)

In :
A_pinv.dot(b)

Out:
array([ 0.2425, -3.5347,  0.9537])


Compare with the solution from QR-based Least Squares:

In :
Q, R = la.qr(A, "complete")
la.solve(R[:3], Q.T.dot(b)[:3])

Out:
array([ 0.2425, -3.5347,  0.9537])