import numpy as np
import numpy.linalg as la
np.set_printoptions(precision=4, linewidth=100)
A = np.random.randn(5, 3)
Now compute the SVD of $A$. Note that numpy.linalg.svd
returns $V^T$:
U, singval, VT = la.svd(A)
V = VT.T
Let's first understand the shapes of these arrays:
print(U.shape)
print(singval.shape)
print(V.shape)
(5, 5) (3,) (3, 3)
Check the orthogonality of $U$ and $V$:
U.T.dot(U)
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]])
V.T.dot(V)
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$:
Sigma = np.zeros(A.shape)
Sigma[:3, :3] = np.diag(singval)
Sigma
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:
U.dot(Sigma).dot(V.T) - A
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:
SigmaInv = np.zeros((3,5))
SigmaInv[:3, :3] = np.diag(1/singval)
SigmaInv
array([[ 0.4042, 0. , 0. , 0. , 0. ], [ 0. , 0.7402, 0. , 0. , 0. ], [ 0. , 0. , 2.0971, 0. , 0. ]])
A_pinv = V.dot(SigmaInv).dot(U.T)
Now use the pseudoinverse to "solve" $Ax=b$ for our tall-and-skinny $A$:
b = np.random.randn(5)
A_pinv.dot(b)
array([ 0.2425, -3.5347, 0.9537])
Compare with the solution from QR-based Least Squares:
Q, R = la.qr(A, "complete")
la.solve(R[:3], Q.T.dot(b)[:3])
array([ 0.2425, -3.5347, 0.9537])