加速矩陣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操作

回答

相關問題