[SOLVED] SVD on a non-square matrix

SVD on a non-square matrix

I'm using numpy.linalg.svd() to get the singular value decomposition of matrices. But I can't get back from the decomposition to the original matrix for a non-square matrix.

For example, for a square matrix :

import numpy as np
n=5
# make a random (n,n) matrix
A= np.reshape( np.random.random_integers(0, 9, size= n**2), (n, n))
#SVD 
U,S,Vh = np.linalg.svd(A) 
# to get A from the SVD back
A_svd = U@np.diag(S)@Vh
#check if its the same 
print(np.allclose(A,A_svd))

I get : >>> True

Now for a non-square matrix, for example A of shape (m,n), then the shape U is (m,m), the shape V is (n,n) and S is a diagonal matrix of (of length k), with k = min(m,n). For example :

import numpy as np
n=5
m= 8
# make a random (m,n) matrix
A= np.reshape( np.random.random_integers(0, 9, size= m*n), (m, n))
#SVD 
U,S,Vh = np.linalg.svd(A)

With the following shapes :

>>> U.shape
(8, 8)
>>> S.shape
(5,)
>>> Vh.shape
(5, 5)

I dont know how to get the matrix A back if I have the svd decomposition though. I cant do a simple multiplication because of the shapes difference. U@np.diag(S)@Vh or np.matmul(U,S,Vh) or with np.dot. So I tried to reshape S and fill it with zeroes.

S_m = np.diag(S) 
S_m.resize((U.shape[1], Vh.shape[0]))
#check if its the same 
print(np.allclose(A,U @S_m@ Vh))

>>> False

Solution

I found an answer here, using diagsvd from scipy.linalg.

import scipy.linalg as la    
A_svd = U@la.diagsvd(S,*A.shape)@Vh