加速矩阵rowMeans操作

考虑以下矩阵，加速矩阵rowMeans操作

nc <- 5000 
nr <- 1024 
m <- matrix(rnorm(nc*nr), ncol=nc)

我希望采取两个rowMeans组相同大小的随机在该矩阵中采取之间的差。

n <- 1000 # group size 

system.time(replicate(100, { 
    ind1 <- sample(seq.int(nc), n) 
    ind2 <- sample(seq.int(nc), n) 
    rowMeans(m[, ind1]) - rowMeans(m[, ind2]) 
}))

这是很慢的，可惜我听不懂Rprof的输出（它似乎大部分的时间用在is.data.frame？）

建议的东西更有效率？

我已经考虑了以下几点：

Rcpp：从我在线阅读，我相信的r rowMeans是相当有效的，所以目前还不清楚这将有助于在这一步。我想确信瓶颈的真正起点在哪里，也许我的整个设计并不理想。如果大部分时间都花在为每个较小的矩阵制作副本上，Rcpp会表现得更好吗？
更新到R-devel，似乎有一个新的.rowMeans功能更有效。有人试过吗？

谢谢。

来源

2012-02-28 baptiste

如果你这样做了采样，子集和差异都在犰狳，我会怀疑你获得一点点。应该足够快以通过RcppArmadillo尝试，不是吗？ – 2012-02-28 01:16:21

这很容易，是的，但希望我可以摆脱纯粹的R.本质上，我会尝试何时/如果所有R方法失败。另外，我没有在Rcpp中管理随机数的经验。 – baptiste 2012-02-28 03:58:40

Rcpp sugar为您提供了相同的数据流R使用:-) – 2012-02-28 03:59:24

每个rowSums()呼叫上的列从m的子集可以被看作是与m之间的矩阵乘法的0或1指示所选择的列的向量。如果并列所有这些载体，你结束了两个矩阵之间的乘法（这是更有效的）：

ind1 <- replicate(100, seq.int(nc) %in% sample(seq.int(nc), n)) 
ind2 <- replicate(100, seq.int(nc) %in% sample(seq.int(nc), n)) 
output <- m %*% (ind1 - ind2)

来源

2012-02-28 02:15:53 flodel

这听起来很有希望，谢谢！我需要说服自己，它做的是正确的事情，但它确实快速而优雅。 – baptiste 2012-02-28 03:16:37

您不需要拨打电话rowMeans。您可以先进行减法，并在结果上调用rowMeans。

x1 <- rowMeans(m[,ind1])-rowMeans(m[,ind2]) 
x2 <- rowMeans(m[,ind1]-m[,ind2]) 
all.equal(x1,x2) 
# [1] TRUE

is.data.frame是在rowMeans完成检查的一部分。

更新：关于R-devel中的.rowMeans，它看起来像是直接调用内部代码（假设do_colsum没有改变）。它的定义为：

.rowMeans <- function(X, m, n, na.rm = FALSE) 
    .Internal(rowMeans(X, m, n, na.rm))

在你的情况，m=1024和n=1000。

来源

2012-02-28 01:19:33

事实上，这比你说的更好，因为OP有200个调用（2 * 100个重复）到'rowMeans'，可以减少到1个。 ..'rm < - rowMeans（m）; system.time（replicate（100，{rm [sample（seq.int（nc），n）] - rm [sample（seq.int（nc），n）]}））'经过0.1秒... – 2012-02-28 01:36:53

@Joshua，你确定采用两个矩阵的差异不会像计算其中一个矩阵的行数那么昂贵吗？毕竟这是相同数量的操作。 – flodel 2012-02-28 01:43:49

@BenBolker。这也是我最初的猜测，rowMeans（m）'可能被存储在'replicate'调用之外，但它不能解决同样的问题。 OP的输出是1024×10;你和我都认为会是1000×10 ... – flodel 2012-02-28 01:53:02

加速矩阵rowMeans操作

回答

相关问题