R性能幂函数

Question

有人有如何加快下面代码的提示吗？特别是避免 for 循环？

J <- 10000
I <- 10000
Y <- matrix(0,J,I)
X <- runif(I,0,1)
P <- runif(I,0,1)
Z <- matrix(runif(n = J*I,0,1),J,I)
K <- matrix(runif(n = J*I,0,1),J,I)

for(j in 1:J){
  for (i in 1:I){
    Y[j,i] <- X[i]^(Z[j,i])*P[i]^(K[j,i])       
  }
}

谢谢!

Answer 1

我认为t(X^t(Z)*P^t(K))会导致相同的结果并且更快。这是一个具有 5 X 5 矩阵和性能评估的可重现示例。

set.seed(543)

### Original Code
J <- 5
I <- 5
Y <- matrix(0,J,I)
X <- runif(I,0,1)
P <- runif(I,0,1)
Z <- matrix(runif(n = J*I,0,1),J,I)
K <- matrix(runif(n = J*I,0,1),J,I)

for(j in 1:J){
  for (i in 1:I){
    Y[j,i] <- X[i]^(Z[j,i])*P[i]^(K[j,i])       
  }
}

# View the result
Y
#           [,1]      [,2]      [,3]       [,4]      [,5]
# [1,] 0.8244760 0.7717289 0.3884273 0.30937614 0.6807137
# [2,] 0.8878758 0.3761184 0.2819624 0.08388850 0.5299624
# [3,] 0.9559749 0.7813653 0.2048310 0.05117558 0.4069641
# [4,] 0.9317235 0.6614524 0.1619824 0.08777542 0.3037913
# [5,] 0.9507279 0.5434549 0.3950076 0.08050582 0.3244810

### A solution without for loop
Y2 <- t(X^t(Z)*P^t(K))

# View the result
Y2
#           [,1]      [,2]      [,3]       [,4]      [,5]
# [1,] 0.8244760 0.7717289 0.3884273 0.30937614 0.6807137
# [2,] 0.8878758 0.3761184 0.2819624 0.08388850 0.5299624
# [3,] 0.9559749 0.7813653 0.2048310 0.05117558 0.4069641
# [4,] 0.9317235 0.6614524 0.1619824 0.08777542 0.3037913
# [5,] 0.9507279 0.5434549 0.3950076 0.08050582 0.3244810

identical(Y, Y2)
# [1] TRUE

### Performance evaluation
library(microbenchmark)

perf <- microbenchmark(
  m1 = { Y <- matrix(0,J,I)
  for(j in 1:J){
    for (i in 1:I){
      Y[j,i] <- X[i]^(Z[j,i])*P[i]^(K[j,i])       
    }
  }},
  m2 = {Y2 <- t(X^t(Z)*P^t(K))},
  times = 100L
)

# View the result
perf
# Unit: microseconds
# expr      min       lq       mean   median       uq      max neval cld
#   m1 3649.287 3858.250 4107.31032 3932.017 4112.965 6240.644   100   b
#   m2   13.365   14.907   21.66753   15.422   26.731   60.658   100  a