PCA 如何计算 `sklearn` 中的转换版本？

Question

我对 sklearn 的 PCA(here is the documentation) 及其与奇异值分解 (SVD) 的关系感到困惑。

在 Wikipedia 我们有，

The full principal components decomposition of X can, therefore, be given as T=WX, where W is a p-by-p matrix of weights whose columns are the eigenvectors of $X^T X$. The transpose of W is sometimes called the whitening or sphering transformation.

后面一解释与SVD的关系，我们有：

X=U $\Sigma W^T$

所以我假设矩阵 W，将样本嵌入到潜在的 space 中（注意矩阵的维度是有意义的）并使用 transform 模块sklearn 中的 class PCA 应该给出与我将观察矩阵乘以 W 相同的结果。但是，我检查了它们，它们不匹配。

我是否遗漏了任何错误或代码中存在错误？

import numpy as np
from sklearn.decomposition import PCA

x = np.random.rand(200).reshape(20,10)
x = x-x.mean(axis=0)
u, s, vh = np.linalg.svd(x, full_matrices=False)
pca = PCA().fit(x)

# transformed version based on WIKI: t = X@vh.T = u@np.diag(s)
t_svd1= x@vh.T
t_svd2= u@np.diag(s)
# the pca transform
t_pca = pca.transform(x)

print(np.abs(t_svd1-t_pca).max()) # should be a small value, but it's not :(
print(np.abs(t_svd2-t_pca).max()) # should be a small value, but it's not :(

Answer 1

理论上的维基百科描述与实际 sklearn 实现之间存在差异，但这不是错误，只是稳定性和可重复性增强。

您几乎已经确定了 PCA 的确切实现，但是为了能够完全重现计算，sklearn 开发人员在他们的实现中又增加了一项实施。问题源于 SVD 的不确定性，即 SVD 没有唯一解。通过设置 U_s = -U 和 W_s = -W，也可以很容易地从您的等式中看出这一点，然后 U_s 和 W_s 也满足：

X=U_s $\Sigma W_s^T$

更重要的是，这在切换 U 和 W 列的符号时也适用。如果我们只是将U和W的第k列的符号取反，等式仍然成立。您可以阅读有关此问题的更多信息 f.e。这里 https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2007/076422.pdf.

PCA 的实现通过强制绝对值中的最高加载值始终为正来解决此问题，具体而言，正在使用方法 sklearn.utils.extmath.svd_flip。这样，无论不确定性方法 np.linalg.svd 生成的向量具有哪个符号，绝对值的加载值都将保持不变，即矩阵的符号将保持不变。

因此，为了使您的代码具有与 PCA 实现相同的结果：

import numpy as np
from sklearn.decomposition import PCA

np.random.seed(41)
x = np.random.rand(200).reshape(20,10)
x = x-x.mean(axis=0)
u, s, vh = np.linalg.svd(x, full_matrices=False)
max_abs_cols = np.argmax(np.abs(u), axis=0)
signs = np.sign(u[max_abs_cols, range(u.shape[1])])
u *= signs
vh *= signs.reshape(-1,1)
pca = PCA().fit(x)

# transformed version based on WIKI: t = X@vh.T = u@np.diag(s)
t_svd1= x@vh.T
t_svd2= u@np.diag(s)
# the pca transform
t_pca = pca.transform(x)

print(np.abs(t_svd1-t_pca).max()) # pretty small value :)
print(np.abs(t_svd2-t_pca).max()) # pretty small value :)

PCA 如何计算 `sklearn` 中的转换版本？

How PCA computes the transformed version in `sklearn`?

python

pca

scikit-learn