Using Numpy (np.linalg.svd) for Singular Value Decomposition

Question

Im reading Abdi & Williams (2010) "Principal Component Analysis", and I'm trying to redo the SVD to attain values for further PCA.

The article states that following SVD:

X = P D Q^t

I load my data in a np.array X.

X = np.array(data)
P, D, Q = np.linalg.svd(X, full_matrices=False)
D = np.diag(D)

But i do not get the above equality when checking with

X_a = np.dot(np.dot(P, D), Q.T)

X_a and X are the same dimensions, but the values are not the same. Am I missing something, or is the functionality of the np.linalg.svd function not compatible somehow with the equation in the paper?

Answer 1

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 np
X = 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 np
X = np.random.normal(size=[20,18])
P, D, Q = np.linalg.svd(X, full_matrices=False)
X_a = P @ diag(D) @ Q
print(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())

Answer 2

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 np
U, 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 np
U, D, V = np.linalg.svd(A)
m, n = A.shape
A_reconstructed = U[:,:n] @ np.diag(D) @ V[:m,:]

or you may use 'full_matrices=False' option in the SVD function;

import numpy as np
U, D, V = np.linalg.svd(A,full_matrices=False)
A_reconstructed = U @ np.diag(D) @ V

Answer 3

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.

Answer 4

Though this post is quite old, I thought it deserves a crucial update. In the above answers, the right singular vectors (typically placed in columns of the matrix V) are said to be given directly as columns from np.linalg.svd(). However, this is incorrect. The matrix return from np.linalg.svd() is Vh, the hermitian or conjugate transpose of V, therefore the right singular vectors are in fact in the rows of Vh. Be careful with this as the matrix itself is square so you cannot determine this correctly using the shape, but you can use reconstruction to test if you are viewing the matrix correctly.