R：forループとスピードコードを取り除く

Question

ループを関数mを使用せずに計算をスピードアップして結果を得たいと思います。再現性の例：我々はないベクトル行列を扱っているので、どのようにこれをベクトル化することができforループの各反復で

N <- 2500 
n <- 500 
r <- replicate(1000, sample(N, n)) 

m <- function(r, N) { 
    ic <- matrix(0, nrow = N, ncol = N) 
    for (i in 1:ncol(r)) { 
    p <- r[, i] 
    ic[p, p] <- ic[p, p] + 1 
    } 
    ic 
} 

system.time(ic <- m(r, N)) 
# user system elapsed 
# 6.25 0.51 6.76 
isSymmetric(ic) 
# [1] TRUE 

？

@ joel.wilsonこの関数の目的は、要素のペアごとの頻度を計算することです。それ以来、ペアワイズ包摂確率を推定することができました。

@Khashaaと@alexis_lazに感謝します。ベンチマーク：それは二重のインデックスの操作を避けるよう

> require(rbenchmark) 
> benchmark(m(r, N), 
+   m1(r, N), 
+   mvec(r, N), 
+   alexis(r, N), 
+   replications = 10, order = "elapsed") 
      test replications elapsed relative user.self sys.self user.child sys.child 
4 alexis(r, N)   10 4.73 1.000  4.63  0.11   NA  NA 
3 mvec(r, N)   10 5.36 1.133  5.18  0.18   NA  NA 
2  m1(r, N)   10 5.48 1.159  5.29  0.19   NA  NA 
1  m(r, N)   10 61.41 12.983  60.43  0.90   NA  NA 

Answer 1

これは大幅に高速化する必要があります

m1 <- function(r, N) { 
    ic <- matrix(0, nrow = N, ncol=ncol(r)) 
    for (i in 1:ncol(r)) { 
    p <- r[, i] 
    ic[, i][p] <- 1 
    } 
    tcrossprod(ic) 
} 

system.time(ic1 <- m1(r, N)) 
# user system elapsed 
# 0.53 0.01 0.55 

all.equal(ic, ic1) 
# [1] TRUE 

シンプルな「カウント/追加」の操作はほとんど常にコアが何であるかを

mvec <- function(r, N) { 
    ic <- matrix(0, nrow = N, ncol=ncol(r)) 
    i <- rep(1:ncol(r), each=nrow(r)) 
    ic[cbind(as.vector(r), i)] <- 1 
    tcrossprod(ic) 
}