Calculate a matrix with sequencing without a nested for loop for faster calculations

Question

I am converting some code over from Excel in which we calculate the values in a matrix based on the element that came before it. This is easy and straightforward in Excel. But in R, I define the first row of the matrix and each subsequent row is calculated based on the one before with the following equation in a nested for loop.

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

My matrices are very large, thousands of rows and columns, so my guess is this is not a very efficient way to make these calculations. I looked into sapply and vapply, but didn't understand how to perform the sequential step of calculating each row based on the previous row.

Answer 1

Just implementing your code in 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;
}

---

Verification:

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

---

Benchmark:

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

Results:

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 

The R code of Hayden Rabel is fairly fast!