机器算法验证 - 为什么要通过数据的 SVD 对数据进行 PCA？ - 吾爱随笔录

为什么要通过数据的 SVD 对数据进行 PCA？

机器算法验证主成分分析算法 svd 矩阵分解

2022-02-08 07:40:28

这个问题是关于计算主成分的有效方法。

许多关于线性 PCA 的文章都提倡使用个案数据的奇异值分解。也就是说，如果我们有数据并想用主成分替换变量（它的列），我们做 SVD：，奇异值（特征值的平方根）占据主对角线，右特征向量是轴变量到轴分量的正交旋转矩阵，左特征向量类似于，仅适用于情况。然后我们可以将组件值计算为。 $\bf X$ $\bf X=USV'$ $\bf S$ $\bf V$ $\bf U$ $\bf V$ $\bf C=XV=US$
对变量进行 PCA 的另一种方法是通过分解方阵（即可以是变量之间的相关性或协方差等）。分解可以是特征分解或奇异值分解：对于方形对称正半定矩阵，它们将给出相同的结果，其特征值作为的对角线，并且如前所述。组件值将是。 $\bf R=X'X$ $\bf R$ $\bf R=VLV'$ $\bf L$ $\bf V$ $\bf C=XV$

现在，我的问题：如果数据是一个大矩阵，并且案例数量（通常是案例）远大于变量数量，那么方式（1）预计比方式（2 ）慢得多)，因为方式 (1) 将非常昂贵的算法（例如 SVD）应用于大矩阵；它计算并存储在我们的案例中我们真正不需要的巨大矩阵如果是这样，那么为什么这么多教科书似乎提倡或仅提及方式（1）？也许它是有效的，我错过了什么？ $\bf X$ $\bf U$

2个回答

SVD 速度较慢，但通常被认为是首选方法，因为它具有更高的数值精度。

正如您在问题中所说，主成分分析 (PCA) 可以通过中心数据矩阵的 SVD （有关更多详细信息，请参阅此 Q&A 线程）或协方差矩阵（或者，如果，请参阅此处了解更多详细信息）。 $\mathbf X$ $\frac{1}{n-1}\mathbf X^\top \mathbf X$ $\mathbf{XX}^\top$ $n\ll p$

这是 MATLAB 的pca()函数帮助中写的内容：

pca用于执行主成分分析的主成分算法[...]：

'svd' -- 默认值。X 的奇异值分解 (SVD)。

'eig' - 协方差矩阵的特征值分解 (EIG)。EIG 算法比 SVD 更快，当观察次数， $n$ , 超过变量的数量, $p$ , 但不太准确，因为协方差的条件数是 X 的条件数的平方。

最后一句话强调了关键的速度-准确性权衡。

您正确地观察到协方差矩阵的特征分解通常比数据矩阵的 SVD 更快。这是 Matlab 中的一个简短基准，随机 $1000\times 100$ 数据矩阵：

X = randn([1000 100]);

tic; svd(X); toc         %// Elapsed time is 0.004075 seconds.
tic; svd(X'); toc        %// Elapsed time is 0.011194 seconds.
tic; eig(X'*X); toc      %// Elapsed time is 0.001620 seconds.
tic; eig(X*X'); toc;     %// Elapsed time is 0.126723 seconds.

在这种情况下，最快的方法是通过协方差矩阵（第三行）。当然，如果 $n \ll p$ （而不是相反）那么这将是最慢的方式，但在这种情况下使用 Gram 矩阵 $\mathbf{XX}^\top$ （第四行）将是最快的方式。无论哪种方式，数据矩阵本身的 SVD 都会变慢。

但是，它会更准确，因为相乘 $\mathbf X$ 本身会导致数值精度损失。这是一个示例，改编自@JM 对Why SVD on的回答 $X$ 优于特征分解 $XX^⊤$ 在Math.SE 上的 PCA中。

考虑一个数据矩阵

X = (\begin{matrix} 1 & 1 & 1 \\ ϵ & 0 & 0 \\ 0 & ϵ & 0 \\ 0 & 0 & ϵ \end{matrix}),

$\mathbf X = \begin{pmatrix}1&1&1\\\epsilon & 0 & 0\\ 0 & \epsilon & 0 \\ 0 & 0 & \epsilon\end{pmatrix},$ 有时称为 Läuchli 矩阵（让我们在此示例中省略居中）。它的平方奇异值为

3 + ϵ^{2}

$3+\epsilon^2$ ,

ϵ^{2}

$\epsilon^2$ ，和

ϵ^{2}

$\epsilon^2$ . 服用

ϵ = 10^{- 5}

$\epsilon = 10^{-5}$ ，我们可以使用 SVD 和 EIG 来计算这些值：

eps = 1e-5;
X = [1 1 1; eye(3)*eps];
display(['Squared sing. values of X: ' num2str(sort(svd(X),'descend').^2')])
display(['Eigenvalues of X''*X:       ' num2str(sort(eig(X'*X),'descend')')])

获得相同的结果：

Squared sing. values of X: 3       1e-10       1e-10
Eigenvalues of X'*X:       3       1e-10       1e-10

但现在服用 $\epsilon = 10^{-10}$ 我们可以观察到 SVD 如何仍然表现良好，但 EIG 出现故障：

Squared sing. values of X: 3       1e-20       1e-20
Eigenvalues of X'*X:       3           0 -3.3307e-16

这里发生的是，协方差矩阵的计算使条件数平方 $\mathbf X$ , 所以特别是在 $\mathbf X$ 有一些几乎共线的列（即一些非常小的奇异值），首先计算协方差矩阵，然后计算其特征分解与直接 SVD 相比会导致精度损失。

我应该补充一点，人们通常很乐意忽略这种潜在的[微小]精度损失，而是使用更快的方法。

这是我关于这个话题的 2ct

我第一次学习 PCA 的化学计量学讲座使用了解决方案 (2)，但它不是面向数值的，而且我的数值讲座只是一个介绍，据我记得没有讨论 SVD。
如果我正确理解Holmes: Fast SVD for Large-Scale Matrices，那么您的想法已被用于获得长矩阵的计算快速 SVD。
这意味着如果遇到合适的矩阵，一个好的 SVD 实现可能会在内部遵循 (2)（我不知道是否还有更好的可能性）。这意味着对于高级实现，最好使用 SVD (1) 并将其留给 BLAS 来处理内部使用的算法。

快速实际检查：OpenBLAS 的 svd 似乎没有做出这种区分，在 5e4 x 100 的矩阵svd (X, nu = 0)上，中位数为 3.5 秒，而svd (crossprod (X), nu = 0)需要 54 毫秒（从 R 调用microbenchmark）。
特征值的平方当然很快，并且两个调用的结果是等价的。

timing  <- microbenchmark (svd (X, nu = 0), svd (crossprod (X), nu = 0), times = 10)
timing
# Unit: milliseconds
#                      expr        min         lq    median         uq        max neval
#            svd(X, nu = 0) 3383.77710 3422.68455 3507.2597 3542.91083 3724.24130    10
# svd(crossprod(X), nu = 0)   48.49297   50.16464   53.6881   56.28776   59.21218    10

更新：看看Wu，W。Massart, D. & de Jong, S.：宽数据的内核 PCA 算法。第一部分：理论和算法，化学计量学和智能实验室系统，36, 165 - 172 (1997)。DOI：http://dx.doi.org/10.1016/S0169-7439(97)00010-5

本文讨论了 PCA 的 4 种不同算法的数值和计算特性：SVD、特征分解 (EVD)、NIPALS 和 POWER。

它们的关系如下：

computes on      extract all PCs at once       sequential extraction    
X                SVD                           NIPALS    
X'X              EVD                           POWER

论文的上下文很宽，它们在（内核 PCA）上工作 - 这与您询问的情况正好相反。因此，要回答有关长矩阵行为的问题，您需要交换“内核”和“经典”的含义。 $\mathbf X^{(30 \times 500)}$ $\mathbf{XX'}$

性能比较

毫不奇怪，EVD 和 SVD 会根据使用的是经典算法还是内核算法而改变位置。在这个问题的上下文中，这意味着根据矩阵的形状，一个或另一个可能更好。

但从他们对“经典”SVD 和 EVD 的讨论中可以清楚地看出，的分解是计算 PCA 的一种非常常用的方法。但是，除了使用 Matlab 的函数外，他们没有指定使用哪种 SVD 算法。 $\mathbf{X'X}$ svd ()

    > sessionInfo ()
    R version 3.0.2 (2013-09-25)
    Platform: x86_64-pc-linux-gnu (64-bit)

    locale:
     [1] LC_CTYPE=de_DE.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8        LC_COLLATE=de_DE.UTF-8     LC_MONETARY=de_DE.UTF-8   
     [6] LC_MESSAGES=de_DE.UTF-8    LC_PAPER=de_DE.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
    [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       

    attached base packages:
    [1] stats     graphics  grDevices utils     datasets  methods   base     

    other attached packages:
    [1] microbenchmark_1.3-0

loaded via a namespace (and not attached):
[1] tools_3.0.2

$ dpkg --list libopenblas*
[...]
ii  libopenblas-base              0.1alpha2.2-3                 Optimized BLAS (linear algebra) library based on GotoBLAS2
ii  libopenblas-dev               0.1alpha2.2-3                 Optimized BLAS (linear algebra) library based on GotoBLAS2

其它你可能感兴趣的问题

上一篇文本分类的词袋：为什么不只使用词频而不是 TFIDF？下一篇回归中的 Wald 检验（OLS 和 GLM）：t 分布与 z 分布