So I am trying to find the generalised inverse A+ of a matrix A using svd
this is my code for svd
A4 = np.array([[3, 1,2,3,4],
[5,2,3,4,5],
[7,4,5,6,7]])
print(A4)
U4, s4, V4 = np.linalg.svd(A4)
D4 = sp.linalg.diagsvd(s4, 3, 5)
print('UDVt =\n', U4@D4@V4 )
Which works fine and gives me the correct matrix
And this is the code that I am trying to use to get the generalised inverse
D4p = np.array([[D4[0,0]**-1,0,0],
[0,D4[1,1]**-1,0],
[0,0,D4[2,2]**-1],
[0,0,0],
[0,0,0]])
A4plus = np.transpose(V4) @ D4p @ np.transpose(U4)
print(np.round(A4plus @ A4,2))
print(np.round(A4 @ A4plus,2))
I have tried changing the transposed matrixed.
Assuming that the rows are not linear dependent (otherwise the S-matrix has elements of value zero and I cannot invert them), You can say
U * S * Vt = A
=> A * (Vt.T * 1/S * U.T) = I
In code:
from numpy import array, round
from numpy.linalg import svd
from scipy.linalg import diagsvd
A = array([[3, 1, 2, 3, 4],
[5, 2, 3, 4, 5],
[7, 4, 5, 6, 7]])
U, S, Vt = svd(A)
S_inv = diagsvd(1/S, A.shape[1], A.shape[0])
A_inv = Vt.T @ S_inv @ U.T
print(round(A @ A_inv, 5))
This gives as a result the pseudo-inverse of A
[[ 1. 0. -0.] [ 0. 1. 0.] [-0. -0. 1.]]
Here is a relevant paragraph about this: https://inst.eecs.berkeley.edu/~ee127/sp21/livebook/def_pseudo_inv.html