CODEWITHSUNDEEP
rperformanceloopsmatrix
minem
minem

Reputation: 3650

R: getting rid of for loop and speeding code

I would like to speed up my calculations and obtain results without using loop in function m. Reproducible example:

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

In every iteration of for loop we are dealing with matrix not vector, so how this could be Vectorized?

@joel.wilson The purpose of this function is to calculate pairwise frequencies of elements. So afterwards we could estimate pairwise inclusion probabilities.

Thanks to @Khashaa and @alexis_laz. Benchmarks:

> 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

Upvotes: 7

Views: 438

Answers (1)

Khashaa
Khashaa

Reputation: 7373

This should be significantly faster as it avoids operations on double indexing

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

Simple "counting/adding" operations can almost always be vectorized

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)
}

Upvotes: 7

Related Questions