How do we use the SVD to invert a matrix \( \boldsymbol{X}^\boldsymbol{X} \) which is singular or near singular? The simple answer is to use the linear algebra function for the computation of the pseudoinverse of a given matrix \( \boldsymbol{X} \), that is
import numpy as np
X = np.array( [ [1,2,3],[2,4,5],[3,5,6]])
Xinv = np.linlag.pinv(X)
Let us first look at a matrix which does not causes problems and write our own function where we just use the SVD.
import numpy as np
# SVD inversion
def SVDinv(A):
''' Takes as input a numpy matrix A and returns inv(A) based on singular value decomposition (SVD).
SVD is numerically more stable than the inversion algorithms provided by
numpy and scipy.linalg at the cost of being slower.
'''
U, s, VT = np.linalg.svd(A)
print('test U')
print( (np.transpose(U) @ U - U @np.transpose(U)))
print('test VT')
print( (np.transpose(VT) @ VT - VT @np.transpose(VT)))
D = np.zeros((len(U),len(VT)))
D = np.diag(s)
UT = np.transpose(U); V = np.transpose(VT); invD = np.linalg.inv(D)
return np.matmul(V,np.matmul(invD,UT))
#X = np.array([ [1.0, -1.0, 2.0], [1.0, 0.0, 1.0], [1.0, 2.0, -1.0], [1.0, 1.0, 0.0] ])
# Non-singular square matrix
X = np.array( [ [1,2,3],[2,4,5],[3,5,6]])
print(X)
A = np.transpose(X) @ X
# Brute force inversion
B = np.linalg.inv(A) # here we could use np.linalg.pinv(A)
C = SVDinv(A)
print(np.abs(B-C))