r - 在没有嵌套 for 循环的情况下使用排序计算矩阵以加快计算速度

Question

我正在从 Excel 转换一些代码，我们在其中根据它之前的元素计算矩阵中的值。这在 Excel 中简单明了。但是在 R 中，我定义了矩阵的第一行，随后的每一行都是根据前面的行在嵌套的 for 循环中使用以下等式计算的。

step1 <- c(0.0013807009, 0.0005997510, 0.0011314072, 0.0016246001, 0.0014240778)
A <- c( 34.648458,  1.705335,  0.000010, 11.312707,  9.167534)
n <- 10

tau <- matrix(0,nrow=n+1,ncol=5)
tau[1,] <- A
for(j in 1:5){
  for(i in 2:nrow(tau)){
    tau[i,j] <- tau[i-1,j] + step1[j]*1.0025^(i-2)
  }
}

我的矩阵非常大，有数千行和几列，所以我猜这不是进行这些计算的非常有效的方法。我查看了 sapply 和 vapply，但不明白如何执行基于前一行计算每一行的顺序步骤。

Answer 1

只需在 Rcpp 中实现您的代码:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix to_col_cumsum(const NumericVector& step1,
                            const NumericVector& A,
                            int n) {

  int m = step1.length();
  NumericMatrix tau(n + 1, m);
  int i, j;

  // precomputing this is important
  NumericVector pows(n + 1);
  for (i = 1; i < (n + 1); i++) pows[i] = pow(1.0025, i - 1);

  for (j = 0; j < m; j++) {
    tau(0, j) = A[j];
    for (i = 1; i < (n + 1); i++) {
      tau(i, j) = tau(i - 1, j) + step1[j] * pows[i];
    }
  }

  return tau;
}

---

验证:

step1 <- c(0.0013807009, 0.0005997510, 0.0011314072, 0.0016246001, 0.0014240778)
A <- c( 34.648458,  1.705335,  0.000010, 11.312707,  9.167534)
n <- 10

# OP
f1 <- function(step1, A, n) {
  m <- length(step1)
  tau <- matrix(0,nrow=n+1,ncol=m)
  tau[1,] <- A
  for(j in 1:m){
    for(i in 2:nrow(tau)){
      tau[i,j] <- tau[i-1,j] + step1[j]*1.0025^(i-2)
    }
  }
  tau
}

# Hayden
f2 <- function(step1, A, n) {
  calc_next_row <- function(tau, row_idx) {
    tau + step1 * 1.0025 ^ row_idx
  }
  do.call(rbind, Reduce(calc_next_row, 
                        init = A, 
                        x = 0:(n - 1), 
                        accumulate = TRUE))
}
all.equal(f2(step1, A, n), f1(step1, A, n))
all.equal(to_col_cumsum(step1, A, n), f1(step1, A, n))

---

基准:

step1 <- runif(1000)
A <- rnorm(1000)
n <- 2000
microbenchmark::microbenchmark(
  HR = f2(step1, A, n), 
  FP = to_col_cumsum(step1, A, n), 
  times = 100
)

结果:

Unit: milliseconds
 expr       min        lq      mean    median        uq       max neval cld
   HR 10.907345 13.127121 18.337656 14.680584 16.419786 131.97709   100   b
   FP  6.516132  7.308756  9.140994  9.139504  9.841078  17.28872   100  a 

Hayden Rabel 的 R 代码相当快!