R中的循环加速

Question

我正在尝试在 R 中为 nrow=300,000 次模拟（在 ncol=30 个变量上）执行以下操作：

投反对票

接受

这是我的代码：

FS_DF <- read.csv("fs.csv", sep = ",")
Y_DF <- read.csv("Y.csv", sep = ",")
CALIBSCENS_DF <- read.csv("calib_scens.csv", sep = ",")
Y_DF$X <- NULL

X_mat <- matrix(1:1, nrow(CALIBSCENS_DF), nrow(FS_DF))

for (irow in 1:nrow(CALIBSCENS_DF)) { 
for (jrow in 1:nrow(FS_DF)) { 
for (krow in 1:ncol(FS_DF)) { 
    X_mat [irow, jrow] <- X_mat[irow, jrow] * (CALIBSCENS_DF[irow, krow] ^ FS_DF[jrow, krow])

}}}

fit <- .lm.fit(X_mat, as.matrix(sapply(Y_DF, as.numeric)))

填充我的 X 矩阵需要很长时间。有人可以建议一种更快的方法来执行此操作。 SCENS_DF、FS_DF 是数据帧。 X_mat 是一个矩阵。

Answer 1

如果此代码是您的瓶颈并且您使用循环，这始终是cpp 可能产生良好结果的好兆头。我们可以使用Rcpp 使其更容易，并在我们的代码中添加 cpp 函数。

您可以在下面找到我使用 Rcpp 的方法以及针对 minem 方法的一些基准测试，减少了大约 20% 的运行时间（很大程度上取决于矩阵的大小）。

library(Rcpp) # load the Rcpp library

# create some data...
CALIBSCENS_DF <- matrix(2:5, nrow = 2)
FS_DF <- matrix(2:5, nrow = 2)

# create the cpp-function, basically the same as yours, just adapted to cpp
cppFunction("
NumericMatrix cpp_fun(NumericMatrix A, NumericMatrix B) {
    NumericMatrix retMax(A.nrow(), B.nrow());

    long double mult;
    for (int irow = 0; irow < A.nrow(); irow++) {
        for (int jrow = 0; jrow < B.nrow(); jrow++) {
            mult = 1;
            for (int krow = 0; krow < B.ncol(); krow++) {
                mult *= pow(A(irow, krow), B(jrow, krow));
            }
            retMax(irow, jrow) = mult;
        }
    }
    return retMax;
}
")
# execute the function called 'cpp_fun' in R
cpp_mat <- cpp_fun(CALIBSCENS_DF, FS_DF)
cpp_mat
# [,1]  [,2]
# [1,] 1024  8192
# [2,] 5625 84375

将函数与 Minem 显示的结果进行比较

# for comparison, use Minems function
minem_fun <- function(A_mat, B_mat) {
  X <- matrix(1, ncol = nrow(B_mat), nrow = nrow(A_mat))
  for (irow in 1:nrow(A_mat)) {
    for (jrow in 1:nrow(B_mat)) {
      X [irow, jrow] <- prod(A_mat[irow, ] ^ B_mat[jrow, ])
    }
  }
  return(X)
}
minem_mat <- minem_fun(CALIBSCENS_DF, FS_DF)

identical(cpp_mat, minem_mat)
# [1] TRUE

速度基准

library(microbenchmark)
# small data
microbenchmark(
  minem = minem_fun(CALIBSCENS_DF, FS_DF),
  cpp = cpp_fun(CALIBSCENS_DF, FS_DF),
  times = 1000
)
# Unit: microseconds
# expr   min     lq      mean median     uq      max neval
# minem 9.386 10.239 11.198179  10.24 10.667   49.915  1000
# cpp 1.707  2.560  3.954538   2.56  2.987 1098.980  1000


# larger data
n <- 200
CALIB_large <- matrix(rnorm(n^2, mean = 100, sd = 10), nrow = n, ncol = n)
FS_large <- matrix(rnorm(n^2, mean = 2, sd = 0.5), nrow = n, ncol = n)

microbenchmark(
  minem = minem_fun(CALIB_large, FS_large),
  cpp = cpp_fun(CALIB_large, FS_large),
  times = 10
)
# Unit: seconds
# expr      min       lq     mean   median       uq      max neval
# minem 1.192011 1.197783 1.209692 1.201320 1.230812 1.238446    10
# cpp 1.009908 1.019727 1.023600 1.025791 1.028152 1.029427    10

这对你有帮助吗？

Answer 2

看起来我们可以这样删除一个循环：

CALIBSCENS_DF <- matrix(2:5, nrow = 2)
FS_DF <- matrix(2:5, nrow = 2)
X <- matrix(1, ncol = nrow(FS_DF), nrow = nrow(CALIBSCENS_DF))
for (irow in 1:nrow(CALIBSCENS_DF)) { 
  for (jrow in 1:nrow(FS_DF)) { 
      X [irow, jrow] <-
        X[irow, jrow] * prod(CALIBSCENS_DF[irow, ] ^ FS_DF[jrow, ])
    }}
X
#      [,1]  [,2]
# [1,] 1024  8192
# [2,] 5625 84375

Answer 3

这还不是您问题的真正答案，但不适合发表评论。我认为我们需要仔细查看您正在尝试执行的操作，并确定 for 循环是否正在执行您认为的操作。

让我们稍微简化一下代码。让我们有矩阵X、C 和F 并定义循环

for (i in 1:nrow(C)){
  for (j in 1:nrow(F)){
    for (k in 1:ncol(F)){
      X[i, j] <- X[i, j] * C[i, k] ^ F[j, k]
    }
  }
}

现在让我们逐步了解循环迭代时会发生什么

i = 1; j = 1; k = 1      X[1, 1] <- X[1, 1] * C[1, 1] ^ F[1, 1]
i = 1; j = 1; k = 2      X[1, 1] <- X[1, 1] * C[1, 2] ^ F[1, 2]
i = 1; j = 1; k = 3      X[1, 1] <- X[1, 1] * C[1, 3] ^ F[1, 3]
...
i = 1; j = 1; k = 30      X[1, 1] <- X[1, 1] * C[1, 30] ^ F[1, 30]

最终，X[1, 1] 依赖于 C[1, 30] 和 F[1, 30]。您已经完成了 29 次被覆盖的迭代。此时，循环将递增j，您将得到

i = 1; j = 2; k = 1      X[1, 2] <- X[1, 2] * C[1, 1] ^ F[2, 1]
i = 1; j = 2; k = 2      X[1, 2] <- X[1, 2] * C[1, 2] ^ F[2, 2]
i = 1; j = 2; k = 3      X[1, 2] <- X[1, 2] * C[1, 3] ^ F[2, 3]
...
i = 1; j = 2; k = 30      X[1, 2] <- X[1, 2] * C[1, 30] ^ F[2, 30]

所以X[1, 2] 依赖于C[1, 30] 和F[2, 30]。

这是您期望的行为吗？