Using Numpy (np.linalg.svd) for Singular Value Decomposition
TL;DR: numpy's SVD computes X = PDQ, so the Q is already transposed.
SVD decomposes the matrix X
effectively into rotations P
and Q
and the diagonal matrix D
. The version of linalg.svd()
I have returns forward rotations for P
and Q
. You don't want to transform Q
when you calculate X_a
.
import numpy as npX = np.random.normal(size=[20,18])P, D, Q = np.linalg.svd(X, full_matrices=False)X_a = np.matmul(np.matmul(P, np.diag(D)), Q)print(np.std(X), np.std(X_a), np.std(X - X_a))
I get: 1.02, 1.02, 1.8e-15, showing that X_a
very accurately reconstructs X
.
If you are using Python 3, the @
operator implements matrix multiplication and makes the code easier to follow:
import numpy as npX = np.random.normal(size=[20,18])P, D, Q = np.linalg.svd(X, full_matrices=False)X_a = P @ diag(D) @ Qprint(np.std(X), np.std(X_a), np.std(X - X_a))print('Is X close to X_a?', np.isclose(X, X_a).all())
From the scipy.linalg.svd docstring, where (M,N) is the shape of the input matrix, and K is the lesser of the two:
Returns-------U : ndarray Unitary matrix having left singular vectors as columns. Of shape ``(M,M)`` or ``(M,K)``, depending on `full_matrices`.s : ndarray The singular values, sorted in non-increasing order. Of shape (K,), with ``K = min(M, N)``.Vh : ndarray Unitary matrix having right singular vectors as rows. Of shape ``(N,N)`` or ``(K,N)`` depending on `full_matrices`.
Vh, as described, is the transpose of the Q used in the Abdi and Williams paper. So just
X_a = P.dot(D).dot(Q)
should give you your answer.
I think there are still some important points for those who use SVD in Python/linalg library. Firstly, https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.svd.html is a good reference for SVD computation function.
Taking SVD computation as A= U D (V^T), For U, D, V = np.linalg.svd(A), this function returns V in V^T form already. Also D contains eigenvalues only, hence it has to be shaped into matrix form. Hence the reconstruction can be formed with
import numpy as npU, D, V = np.linalg.svd(A)A_reconstructed = U @ np.diag(D) @ V
The point is that, If A matrix is not a square but rectangular matrix, this won't work, you can use this instead
import numpy as npU, D, V = np.linalg.svd(A)m, n = A.shapeA_reconstructed = U[:,:n] @ np.diag(D) @ V[:m,:]
or you may use 'full_matrices=False' option in the SVD function;
import numpy as npU, D, V = np.linalg.svd(A,full_matrices=False)A_reconstructed = U @ np.diag(D) @ V