计算R中前两个主成分的最快方法是啥？

Posted 2023-02-16

【中文标题】计算R中前两个主成分的最快方法是啥？

【英文标题】：What is the fastest way to calculate first two principal components in R?计算R中前两个主成分的最快方法是什么？

【发布时间】：2012-01-08 03:00:54

【问题描述】：

我在 R 中使用 princomp 来执行 PCA。我的数据矩阵很大（10K x 10K，每个值最多 4 个小数点）。在 Xeon 2.27 GHz 处理器上需要约 3.5 小时和约 6.5 GB 的物理内存。

由于我只想要前两个组件，有没有更快的方法来做到这一点？

更新：

除了速度之外，有没有一种内存有效的方法来做到这一点？

使用svd(,2,) 计算前两个组件需要大约 2 小时和大约 6.3 GB 的物理内存。

【问题讨论】：

可以使用 NIPALS 算法。搜索 R 包。

【参考方案1】：

您可以自己编写函数并在 2 个组件处停止。这不是太难。我把它放在某个地方，如果我找到它，我会发布它。

【讨论】：

可能你能给出函数的逻辑，我可以尝试自己编码！

作为 PCA 的入门者，我写了一篇博客文章，试图用 OLS 来解释这一点：cerebralmastication.com/2010/09/… 在底部有一个链接指向我找到的 Lindsay I Smith 的文章真的很有帮助。史密斯 PDF 链接：cs.otago.ac.nz/cosc453/student_tutorials/…

@JD Long：这是一篇有趣的文章。让我试试！

Bioc 项目中的 pcaMethods 包可能值得一看。我不知道它有多快，但它是另一个参考点。 bioconductor.org/packages/release/bioc/html/pcaMethods.html

【参考方案2】：

您有时可以访问所谓的“经济”分解，它允许您限制特征值/特征向量的数量。看起来 eigen() 和 prcomp() 不提供此功能，但 svd() 允许您指定要计算的最大数量。

在小矩阵上，收益似乎不大：

R> set.seed(42); N <- 10; M <- matrix(rnorm(N*N), N, N)
R> library(rbenchmark)
R> benchmark(eigen(M), svd(M,2,0), prcomp(M), princomp(M), order="relative")
          test replications elapsed relative user.self sys.self user.child
2 svd(M, 2, 0)          100   0.021  1.00000      0.02        0          0
3    prcomp(M)          100   0.043  2.04762      0.04        0          0
1     eigen(M)          100   0.050  2.38095      0.05        0          0
4  princomp(M)          100   0.065  3.09524      0.06        0          0
R> 

但是在将 princomp() 从 svd() 重构为 svd() 时，相对于 princomp() 的三倍可能值得您在两个值之后停止。

【讨论】：

当 N=200 时，我的机器 princomp 最快（不是很多，基本上等于 svd(,2,)，因此结果可能会因处理器和缩放而异。

在 rbenchmark 包中。还有一个 microbenchmark 包。

corpcor 包中的

fast.svd 非常快。

【参考方案3】：

您可以使用神经网络方法来找到主成分。 这里给出基本描述.. http://www.heikohoffmann.de/htmlthesis/node26.html

第一个主成分，y= w1*x1+w2*x2 第二正交分量可以计算为q = w2*x1-w1*x2。

【参考方案4】：

power method 可能是您想要的。如果您在 R 中编写代码，这一点也不难，我想您可能会发现它并不比其他答案中建议的 SVD 方法快，它使用 LAPACK 编译例程。

【讨论】：

我建议不要这样做，因为幂法收敛速度极慢。

这在很多情况下都是正确的。速度取决于最大特征值与下一个特征值的相对大小；所以这将取决于问题。尽管如此，我认为如果只寻找两个特征向量并且矩阵非常大，该方法可能具有竞争力。不尝试就无法知道。

【参考方案5】：

'svd' 包提供了通过 Lanczos 算法进行截断 SVD/特征分解的例程。您可以使用它来计算前两个主成分。

这里有：

> library(svd)
> set.seed(42); N <- 1000; M <- matrix(rnorm(N*N), N, N)
> system.time(svd(M, 2, 0))
   user  system elapsed 
  7.355   0.069   7.501 
> system.time(princomp(M))
   user  system elapsed 
  5.985   0.055   6.085 
> system.time(prcomp(M))
   user  system elapsed 
  9.267   0.060   9.368 
> system.time(trlan.svd(M, neig = 2))
   user  system elapsed 
  0.606   0.004   0.614 
> system.time(trlan.svd(M, neig = 20))
   user  system elapsed 
  1.894   0.009   1.910
> system.time(propack.svd(M, neig = 20))
   user  system elapsed 
  1.072   0.011   1.087 

【讨论】：

由于我的数据是方阵，有没有办法只输入上/下三角矩阵到任何函数（svd、pr​​incomp、prcomp）？这样可以节省将下三角形复制为上三角形的内存消耗步骤！

我认为这对于“常规”功能是不可能的。对于“svd”包中的内容，您可以使用所谓的“外部矩阵接口”，您只需定义如何将矩阵乘以向量，仅此而已。目前这个 API 只是 C 级的，但有传言说一切都将很快传播到普通的 R 级，所以人们可以在 R 中编写自己的例程（当然可以利用矩阵的对称性或稀疏性）。

【参考方案6】：

我尝试了 pcaMethods 包的 nipals 算法实现。默认情况下，它计算前 2 个主成分。结果比其他建议的方法慢。

set.seed(42); N <- 10; M <- matrix(rnorm(N*N), N, N)
library(pcaMethods)
library(rbenchmark)
m1 <- pca(M, method="nipals", nPcs=2)
benchmark(pca(M, method="nipals"),
          eigen(M), svd(M,2,0), prcomp(M), princomp(M), order="relative")

                       test replications elapsed relative user.self sys.self
3              svd(M, 2, 0)          100    0.02      1.0      0.02        0
2                  eigen(M)          100    0.03      1.5      0.03        0
4                 prcomp(M)          100    0.03      1.5      0.03        0
5               princomp(M)          100    0.05      2.5      0.05        0
1 pca(M, method = "nipals")          100    0.23     11.5      0.24        0

【参考方案7】：

“gmodels”和“corpcor”R 软件包提供了更快的 SVD 和 PCA 实现。这些执行类似于小型矩阵的核心版本：

> set.seed(42); N <- 10; M <- matrix(rnorm(N*N), N*N, N)
> library("rbenchmark")
> library("gmodels")    
> benchmark(svd(M,2,0), svd(M), gmodels::fast.svd(M), corpcor::fast.svd(M), prcomp(M), gmodels::fast.prcomp(M), princomp(M), order="relative")
                     test replications elapsed relative user.self sys.self user.child sys.child
1            svd(M, 2, 0)          100   0.005      1.0     0.005    0.000          0         0
2                  svd(M)          100   0.006      1.2     0.005    0.000          0         0
3    gmodels::fast.svd(M)          100   0.007      1.4     0.006    0.000          0         0
4    corpcor::fast.svd(M)          100   0.007      1.4     0.007    0.000          0         0
6 gmodels::fast.prcomp(M)          100   0.014      2.8     0.014    0.000          0         0
5               prcomp(M)          100   0.015      3.0     0.014    0.001          0         0
7             princomp(M)          100   0.030      6.0     0.029    0.001          0         0
> 

但是，它们为较大的矩阵（尤其是具有许多行的矩阵）提供了更快的结果。

> set.seed(42); N <- 10; M <- matrix(rnorm(N*N), N*N*N, N)
> library("rbenchmark")
> library("gmodels")
> benchmark(svd(M,2,0), svd(M), gmodels::fast.svd(M), corpcor::fast.svd(M), prcomp(M), gmodels::fast.prcomp(M), order="relative")

                     test replications elapsed relative user.self sys.self user.child sys.child
4    corpcor::fast.svd(M)          100   0.029    1.000     0.028    0.001          0         0
3    gmodels::fast.svd(M)          100   0.035    1.207     0.033    0.001          0         0
2                  svd(M)          100   0.037    1.276     0.035    0.002          0         0
1            svd(M, 2, 0)          100   0.039    1.345     0.037    0.001          0         0
5               prcomp(M)          100   0.068    2.345     0.061    0.006          0         0
6 gmodels::fast.prcomp(M)          100   0.068    2.345     0.060    0.007          0         0

【讨论】：

您的基准测试很好地表明 gmodels 函数实际上并没有更快。

这取决于您使用的是 PCA 还是 SVD。该问题还特别与大型矩阵的性能有关。

35 毫秒而不是 37 毫秒并不是真的更快。与 OP 的 10000 平方相比，1000x10 仍然非常小。您可能还打算在 rnorm 调用中添加 *Ns，目前您正在测试所有列都相同的 matices（R 的众多设计缺陷之一），这可能不是一个理想的测试用例。这两个软件包都声称仅对胖/宽矩阵有优势，但即使在那里，我也没有观察到快速测试的真正优势。如果您有时间解决这些问题，您的回答将与Kevin Wright's answer 一样有用。

是的，这不是一个理想的基准测试。在发布此内容时，我没有太多时间来运行大型矩阵。目的不是广泛测试或给出正确答案，而是提出更多选项（使用与该答案相同的基准测试）。我建议任何人认真应用它来尝试更大的测试作业，然后再将它应用到更大的矩阵，并考虑到由于开销，结果可能与更小的矩阵不同。

【参考方案8】：

我很惊讶还没有人提到irlba 包：

比svd的propack.svd还要快一点， 提供irlba::prcomp_irlba(X, n=2)stats::prcomp-like 接口，方便 对于不同大小的矩形矩阵 (2:1)，不需要在以下基准中调整参数。对于大小为 6000x3000 的矩阵，它比 stats::prcomp 快 50 倍。不过对于小于 100x50 的矩阵，stats::svd 仍然更快。

library(microbenchmark)
library(tidyverse)
#install.packages("svd","corpcor","irlba","rsvd")

exprs <- rlang::exprs(
  svd(M, 2, 2)$v,
  prcomp(M)$rotation[,1:2],
  irlba::prcomp_irlba(M, n=2)$rotation,
  irlba::svdr(M, k=2)$v,
  rsvd::rsvd(M, 2)$v,
  svd::propack.svd(M, neig=2, opts=list(maxiter=100))$v,
  corpcor::fast.svd(M)$v[,1:2]
)

set.seed(42)
tibble(N=c(10,30,100,300,1000,3000)) %>%
  group_by(N) %>%
  do(
    M <- scale(matrix(rnorm(.$N*.$N*2), .$N*2, .$N))
    microbenchmark(!!!exprs,
      times=min(100, ceiling(3000/.$N)))%>%
      as_tibble
  ) %>% 
ggplot(aes(x=N, y=time/1E9,color=expr)) +
  geom_jitter(width=0.05) +
  scale_x_log10("matrix size (2N x N)") +
  scale_y_log10("time [s]") +
  stat_summary(fun.y = median, geom="smooth") +
  scale_color_discrete(labels = partial(str_wrap, width=30))

rsvd 提供的随机 svd 甚至更快，但不幸的是，相当偏离：

set.seed(42)
N <- 1000
M <- scale(matrix(rnorm(N^2*2), N*2, N))
cor(set_colnames(sapply(exprs, function(x) eval(x)[,1]), sapply(exprs, deparse)))

                                                       svd(M, 2, 2)$v prcomp(M)$rotation[, 1:2] irlba::prcomp_irlba(M, n = 2)$rotation irlba::svdr(M, k = 2)$v rsvd::rsvd(M, 2)$v svd::propack.svd(M, neig = 2, opts = list(maxiter = 100))$v corpcor::fast.svd(M)$v[, 1:2]
svd(M, 2, 2)$v                                                   1.0000000                 1.0000000                             -1.0000000               0.9998748           0.286184                                                   1.0000000                     1.0000000
prcomp(M)$rotation[, 1:2]                                        1.0000000                 1.0000000                             -1.0000000               0.9998748           0.286184                                                   1.0000000                     1.0000000
irlba::prcomp_irlba(M, n = 2)$rotation                          -1.0000000                -1.0000000                              1.0000000              -0.9998748          -0.286184                                                  -1.0000000                    -1.0000000
irlba::svdr(M, k = 2)$v                                          0.9998748                 0.9998748                             -0.9998748               1.0000000           0.290397                                                   0.9998748                     0.9998748
rsvd::rsvd(M, 2)$v                                               0.2861840                 0.2861840                             -0.2861840               0.2903970           1.000000                                                   0.2861840                     0.2861840
svd::propack.svd(M, neig = 2, opts = list(maxiter = 100))$v      1.0000000                 1.0000000                             -1.0000000               0.9998748           0.286184                                                   1.0000000                     1.0000000
corpcor::fast.svd(M)$v[, 1:2]                                    1.0000000                 1.0000000                             -1.0000000               0.9998748           0.286184                                                   1.0000000                     1.0000000

如果数据实际上具有结构，这可能会更好。