计算R中前两个主成分的最快方法是什么？
问题描述 投票：17回答：8
我在R中使用princomp来执行PCA。我的数据矩阵很大（10K x 10K，每个值最多4个小数点）。在Xeon 2.27 GHz处理器上需要大约3.5小时和~6.5 GB的物理内存。
由于我只想要前两个组件，有没有更快的方法呢？
更新：
除了速度之外，还有一种内存有效的方法吗？
使用svd(,2,)计算前两个组件需要大约2小时和~6.3 GB的物理内存。
8个回答
19
投票
投票
您有时可以访问所谓的“经济”分解，这些分解允许您限制特征值/特征向量的数量。看起来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()重建princomp()时，相对于svd()的因子可能是值得的，因为svd()允许你在两个值之后停止。
6
投票
投票
'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
4
投票
投票
我尝试了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
1
投票
投票
power method可能是你想要的。如果你用R编写它，这根本不是很难，我想你可能会发现它并不比其他答案中建议的SVD方法快，后者利用了LAPACK编译的例程。
0
投票
投票
您可以使用神经网络方法来查找主要组件。这里给出了基本的描述.. http://www.heikohoffmann.de/htmlthesis/node26.html
第一主分量y = w1 * x1 + w2 * x2和第二正交分量可以被计算为q = w2 * x1-w1 * x2。
0
投票
投票
“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
0
投票
投票
我很惊讶没人提到irlba包：
它甚至比svd的propack.svd快一点，为方便起见提供了irlba::prcomp_irlba(X, n=2)类似stats::prcomp的接口，并且在以下基准测试中不需要参数调整以适应不同大小的矩形矩阵（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
当数据实际上具有结构时，这可能会更好。
-1
投票
投票
您可以自己编写函数并停止在2个组件上。这不是太困难。我把它放在某个地方，如果我发现它我会发布它。
最新问题
- 如何将方程符号作为变量传递
- 根据正确值列表更正输入字符串
- ActiveMQ Artemis 代理连接到 ActiveMQ Classic
- 更改 Metronic Angular 以与 Firebase 配合使用
- 成功获取请求后，React useState 变量未更新
- Angular 12 数据路径“”不得具有其他属性（styleext）
- 重复的工作表及其脚本
- FastAPI TypeError：“ModelMetaclass”类型的对象不可 JSON 序列化
- 如何使用 Google Drive 作为一种 CDN
- Shopify Webhooks 未重定向到 Laravel Route
- Async/Await 异常处理混乱
- 通过 StoreKit 进行应用内购买可以在 macOS 上的 Developer-ID 签名应用中使用吗？
- NQueens - 递归回溯问题
- 透视数据帧值及其单位
- C# 的每个字符串颜色是什么，例如“<color=lime></color>”
- 如何在firebase genkit中将媒体文件发送到流http端点？
- csharp 的每个字符串颜色是什么，例如“<color=lime></color>”
- 为什么在函数原型中使用“[*]”而不是“[]”？
- 从数组中选择数组？
- PDB 文件中配体的 Bond Type 全部显示为 SINGLE
© www.soinside.com 2019 - 2024. All rights reserved.