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

Question

我正在阅读Abdi和Williams（2010年）的《主成分分析》，我想重新做一次奇异值分解（SVD），以便获得进一步进行主成分分析所需的数值。

文章中提到，经过奇异值分解后，可以得到：

X = P D Q^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

虽然这篇帖子有点旧，但我觉得有必要更新一下。在之前的回答中，提到的右奇异向量（通常放在矩阵V的列里）是直接从np.linalg.svd()得到的。然而，这个说法是不对的。np.linalg.svd()返回的矩阵是Vh，也就是V的共轭转置，因此右奇异向量实际上是在Vh的行里。要注意这一点，因为这个矩阵是方形的，所以你不能仅仅通过形状来判断，但你可以通过重构来测试你是否正确理解了这个矩阵。

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）时，还有一些重要的注意事项。首先，这个链接是一个很好的参考，可以帮助你理解SVD的计算函数。

在进行SVD计算时，可以把它看作是 A = U D (V^T)。当你用代码 U, D, V = np.linalg.svd(A) 来计算时，这个函数返回的 V 已经是 V^T 的形式了。而且 D 只包含特征值，所以需要把它转换成矩阵的形式。这样你就可以用下面的方式来重构矩阵：

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

需要注意的是，如果 A 矩阵不是方阵，而是一个长方形的矩阵，这种方法就不适用了，你可以用下面的方式来处理：

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

Answer 4

简而言之：numpy的SVD计算的是X = PDQ，所以Q已经是转置的了。

SVD可以把矩阵X有效地分解成旋转矩阵P和Q，还有一个对角矩阵D。我用的linalg.svd()这个版本返回的是P和Q的正向旋转。在计算X_a的时候，你不需要对Q进行变换。

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, 1.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())