r中嵌套for循环的替代方法或可能需要的Rcpp？答案

【问题标题】：an alternative to nested for loops in r OR possible Rcpp needed?r中嵌套for循环的替代方法或可能需要的Rcpp？
【发布时间】：2021-04-19 10:45:34
【问题描述】：

我有一个嵌套循环

X <- matrix(c(0.5,0,0,0.75), nrow = 2)
k = nrow(X)

ans1 <- 0
ans2 <- 0
for (aa in 1:k) {
  for (bb in 1:k) {
    for (cc in 1:k) {
      for (dd in 1:k) {
        ans1 = ans1 + (0.45 * X[aa,bb] * X[cc,dd])
        for (xx in 1:k) {
          for (yy in 1:k){ 
            ans2 = ans2 + (1.7*X[aa,bb]*X[xx,yy]*X[cc,dd] + 0.2*X[aa,xx]*X[bb,yy]*X[cc,dd])  
          }
        }
      }
    }
  }
}

但是必须是方阵的矩阵X 可以是非常高的维度。因此，这将大大减慢循环。例如X = matrix(rnorm(10000),nrow = 100,byrow = T)

我想知道是否有更短的压缩方法。这将是可读的，最重要的是稍微快一点。我试过expand.grid，但这并没有多大帮助。

例如ans1

library(tidyverse)
an1 <- expand.grid(rep(list(seq(length(X))),2)) %>% arrange_all()
an11 <- t(apply(an1, 1, function(x) as.vector(t(X))[x]))

但正如我所提到的，这并不能提高速度。有什么建议么？我也认为 Rcpp 可能会有所帮助，但我不确定，我还没有尝试过（c++ 语法不太好）。

【问题讨论】：

尽量分解。例如，X[aa,bb] 可以从 cc 和 dd 循环中删除。
ans1 不等于0.45*sum(X)^2 和ans2 不等于1.9*sum(X)^3？这些将很快评估。

标签： r loops for-loop rcpp

【解决方案1】：

您根本不需要使用循环。由于 ans1 和 ans2 的代码只是术语的总和，并且这些术语根本不交互，因此表达式简化为

ans1simple <- 0.45*sum(X)^2
ans2simple <- 1.9*sum(X)^3

您可以对随机数据进行测试。如果您不确定，请更改 X 的种子或大小：

set.seed(123)

X <- matrix(rnorm(9), nrow = 3)
k = nrow(X)

ans1 <- 0
ans2 <- 0
for (aa in 1:k) {
  for (bb in 1:k) {
    for (cc in 1:k) {
      for (dd in 1:k) {
        ans1 = ans1 + (0.45 * X[aa,bb] * X[cc,dd])
        for (xx in 1:k) {
          for (yy in 1:k){ 
            ans2 = ans2 + (1.7*X[aa,bb]*X[xx,yy]*X[cc,dd] + 0.2*X[aa,xx]*X[bb,yy]*X[cc,dd])  
          }
        }
      }
    }
  }
}

ans1simple <- 0.45*sum(X)^2
ans2simple <- 1.9*sum(X)^3
ans1 - ans1simple
#> [1] 2.220446e-16
ans2 - ans2simple
#> [1] -7.993606e-15

^{由reprex package (v1.0.0) 于 2021-04-19 创建}

差异只是舍入误差。

【讨论】：

【解决方案2】：

与 C++ 中的for 循环相比，R 中的for 循环非常慢。

C++ for 循环语法与 R 的某些风格没有太大区别。

我高度怀疑您可以显着压缩您的代码。但只是按照你非常嵌套的语法：

Rcpp 函数：

//[[Rcpp::export]]
Rcpp::NumericVector foo(Rcpp::NumericMatrix& X) {
    Rcpp::NumericVector ans(2);
    int k = X.rows();
    for (int aa = 0; aa < k; ++aa) {
        for (int bb = 0; bb < k; ++bb) {
            for (int cc = 0; cc < k; ++cc) {
                for (int dd = 0; dd < k; ++dd) {
                    ans[0] += 0.45 * X[aa, bb] * X[cc, dd]l;
                    for (int xx = 0; xx < k; ++xx) {
                        for (int yy = 0; yy < k; ++yy) {
                            ans[1] += (1.7 * X[aa, bb] * X[xx, yy] * X[cc, dd] + 0.2 * X[aa, xx] * X[bb, yy] * X[cc, dd]);
                        }
                    }
                }
            }
        }
    }
    return ans;
}

在 R 端：

X <- matrix(c(0.5,0,0,0.75), nrow = 2)
ans <- foo(X)
ans1 <- ans[1]
ans2 <- ans[2]

反复使用上述代码并不是不优化代码的借口。再次，减少循环次数。你不应该都需要它们。

在看到@user2554330 的回答后（我怀疑是这种情况，但懒得解决），Rcpp 实现不会比 R 实现快得多（我怀疑你担心这样的边际收益）

【讨论】：