使用Numpy(np.linalg.svd)进行奇异值分解

Question

我正在阅读Abdi&Williams(2010)的"主成分分析",我正在尝试重做SVD以获得进一步PCA的值.

文章指出在SVD之后:

X = PDQ ^ t

我在np.array X中加载我的数据.

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

但是在检查时我没有得到上述相同

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

X_a和X是相同的尺寸,但值不相同.我是否遗漏了某些内容,或者np.linalg.svd函数的功能是否与文中的等式不兼容？

Answer 1

TL; DR:numpy的SVD计算X = PDQ,因此Q已经转置.

SVD分解矩阵X有效成旋转P和Q和对角矩阵D.linalg.svd()我的版本为P和返回前向旋转Q.你不想Q在计算时进行变换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))

我得到:1.02,1.02,2.8e-15,表明X_a非常准确地重建X.

如果您使用的是Python 3,则@运算符将实现矩阵乘法并使代码更易于遵循:

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

从 scipy.linalg.svd 文档字符串，其中 (M,N) 是输入矩阵的形状，K 是两者中的较小者：

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 是 Abdi 和 Williams 论文中使用的 Q 的转置。所以就

X_a = P.dot(D).dot(Q)

应该给你答案。

Answer 3

我认为对于那些在 Python/linalg 库中使用 SVD 的人来说，仍然有一些重要的点。首先，https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.svd.html是 SVD 计算函数的一个很好的参考。

将 SVD 计算为 A=UD (V^T)，对于 U，D，V = np.linalg.svd(A)，该函数已经以 V^T 形式返回 V。此外，D 仅包含特征值，因此必须将其整形为矩阵形式。因此重建可以用

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

关键是，如果矩阵不是正方形而是矩形矩阵，这将不起作用，您可以改用它

import numpy as np
U, D, V = np.linalg.svd(A)
m, n = A.shape
A_reconstructed = U[:,:n] @ np.diag(D) @ V[:m,:]

或者您可以在 SVD 函数中使用 'full_matrices=False' 选项；

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