避免在嵌套循环中使用太多 "declare same elements of Rcpp List"

Avoid too many "declare same elements of Rcpp List" in nested loops

我想问一下如何在嵌套循环中访问List来节省计算时间。这是一个包含两个函数的 Rcpp 示例:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
#include <Rcpp.h>

// [[Rcpp::export]]
void first(const Rcpp::List ETA
          ,const int N
          ,const int I
                      ) {   
    for(int i = 0; i < I; ++i) {
        arma::rowvec eta = ETA[i];
        for(int n = 0; n < N; ++n) {
            // use eta to do some calculation
        }
    }
}

// [[Rcpp::export]]
void second(const Rcpp::List ETA
           ,const int N
           ,const int I
                      ) {   
    for(int n = 0; n < N; ++n) {
        for(int i = 0; i < I; ++i) {
            arma::rowvec eta = ETA[i];
            // use eta to do some calculation
        }
    }
}

比较R中的时间:

Rcpp::sourceCpp("test.cpp") # save the above code as "test.cpp"
ETA = list()
N = 10^8
I = 10
for (i in 1:I) ETA[[i]] = rep(0,i) # just an example.
ptm <- proc.time(); first(ETA,N,I); print((proc.time() - ptm))
#    user  system elapsed 
#       0       0       0 
ptm <- proc.time(); second(ETA,N,I); print((proc.time() - ptm))
#    user  system elapsed 
#   16.69    0.00   16.69

这里的 ETA 是一个列表,其每个元素可以具有动态长度(向量)或动态维度(矩阵)。在这段代码中,第一种方式比第二种方式快得多。但出于实际需要,当有其他变量迭代n.

时,第二种方式可以减少计算时间

问题:

对于第一种方式或第二种方式,我们能否在循环外(之前)声明 eta,这样我们就不需要多次声明相同的 eta

您在每个循环中都在进行隐式转换和深度复制,因此编译器无法优化它也就不足为奇了。

您可以做的是预先计算所有的 etas 并将它们存储在一个向量中。 (我在循环中填充了一些工作,因为可以完全优化空循环)。

// [[Rcpp::export]]
double third(const Rcpp::List ETA
              ,const int N
              ,const int I) {   
  double output=0;
  
  // pre-create eta rowvecs
  std::vector<arma::rowvec> eta_vec(Rf_xlength(ETA));
  for(int i = 0; i < I; ++i) {
    eta_vec[i] = arma::rowvec(Rcpp::NumericVector(ETA[i]));
  }
  
  for(int n = 0; n < N; ++n) {
    for(int i = 0; i < I; ++i) {
      output += sum(eta_vec[i]);
    }
  }
  return output;
}

结果:

> ptm <- proc.time(); first(ETA,N,I); print((proc.time() - ptm))
[1] 0
   user  system elapsed 
  2.761   0.000   2.761 
> #    user  system elapsed 
> #       0       0       0 
> ptm <- proc.time(); second(ETA,N,I); print((proc.time() - ptm))
[1] 0
   user  system elapsed 
 29.935   0.000  29.930 
> 
> ptm <- proc.time(); third(ETA,N,I); print((proc.time() - ptm))
[1] 0
   user  system elapsed 
  2.599   0.000   2.598