使用 Rcpp 的 R 代码的效率和速度答案

【问题标题】：Efficiency and Speed of R code using Rcpp使用 Rcpp 的 R 代码的效率和速度
【发布时间】：2015-04-16 22:32:26
【问题描述】：

这篇文章是关于使用 Rcpp 包加速 R 代码以避免递归循环。

我的输入由以下示例（长度 7）定义，它是我使用的 data.frame（长度 51673）的一部分：

S=c(906.65,906.65,906.65,906.65,906.65,906.65,906.65)
T=c(0.1371253,0.1457896,0.1248953,0.1261278,0.1156931,0.0985253,0.1332596)
r=c(0.013975,0.013975,0.013975,0.013975,0.013975,0.013975,0.013975)             
h=c(0.001332596,0.001248470,0.001251458,0.001242143,0.001257921,0.001235755,0.001238440)       
P=c(3,1,5,2,1,4,2)
A= data.frame(S=S,T=T,r=r,h=h,P=P)

         S         T        r           h   Per 
1   906.65 0.1971253 0.013975 0.001332596   3   
2   906.65 0.1971253 0.013975 0.001248470   1   
3   906.65 0.1971253 0.013975 0.001251458   5   
4   906.65 0.1971253 0.013975 0.001242143   2   
5   906.65 0.1971253 0.013975 0.001257921   1   
6   906.65 0.1971253 0.013975 0.001235755   4  
7   906.65 0.1971253 0.013975 0.001238440   2

参数为：

w=0.001; b=0.2; a=0.0154; c=0.0000052; neta=-0.70

我有以下我想使用的函数代码：

F<-function(x,w,b,a,c,neta,S,T,r,P){
    u=1i*x 
    nu=(1/(neta^2))*(((1-2*neta)^(1/2))-1)
    # Recursion back to time t
    # Terminal condition for the A and B 
    A_Q=0
    B_Q=0
    steps<-round(T*250,0)  

    for (j in 1:steps){
      A_Q= A_Q+ r*u + w*B_Q-(1/2)*log(1-2*a*(neta^4)*B_Q)
      B_Q= b*B_Q+u*nu+ (1/neta^2)*(1-sqrt((1-2*a*(neta^4)*B_Q)*( 1- 2*c*B_Q - 2*u*neta)))
    }
    F= exp(log(S)*u + A_Q + B_Q*h[P])
  return(F)
}

S = A$S ; r= A$r ; T= A$T ; P=A$P;  h= A$h

然后我想使用我的 Data.set 应用前面的函数，将长度为 N= 100000 的向量：

  Z=length(S); N=100000  ; alpha=2 ; delta= 0.25     
  lambda=(2*pi)/(N*delta) 

res = matrix(nrow=N, ncol=Z)
  for (i in 1:N){
    for (j in 1:Z){
      res[i,j]= Re(F(((delta*(i-1))-(alpha+1)*1i),w,b,a,c,neta,S[j],T[j],r[j],P[j]))
    }
  }

但这需要很多时间：执行这行代码需要 20 秒 N=100 但我想执行 N=100000 次，整体运行时间可能需要几个小时。如何使用 Rcpp 对上述代码进行微调，以减少执行时间并获得高效的程序？

是否可以减少执行时间，如果可以，即使没有 Rcpp，也请建议我一个解决方案。

谢谢。

【问题讨论】：

对于计算 x= (((delta*(i-1))-(alpha+1)*1i) 其中 1i 是 r 中复数的定义，并且 i 取值之间1到N

标签： r performance for-loop vectorization rcpp

【解决方案1】：

您的函数F 可以通过利用Armadillo library（通过RcppArmadillo 包访问）中的vec 和cx_vec 类很容易地转换为C++ - 它对矢量化计算有很好的支持.

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

// [[Rcpp::export]]
arma::cx_vec Fcpp(const arma::cx_vec& x, double w, double b, double a, double c,
                  double neta, const arma::vec& S, const arma::vec& T, 
                  const arma::vec& r, Rcpp::IntegerVector P, Rcpp::NumericVector h) {

  arma::cx_vec u = x * arma::cx_double(0.0,1.0);

  double nu = (1.0/std::pow(neta,2.0)) * (std::sqrt(1.0-2.0*neta)-1.0);
  arma::cx_vec A_Q(r.size());
  arma::cx_vec B_Q(r.size()); 
  arma::vec steps = arma::round(T*250.0);

  for (size_t j = 0; j < steps.size(); j++) {
    for (size_t k = 0; k < steps[j]; k++) {
      A_Q = A_Q + r*u + w*B_Q - 
              0.5*arma::log(1.0 - 2.0*a*std::pow(neta,4.0)*B_Q);
      B_Q = b*B_Q + u*nu + (1.0/std::pow(neta,2.0)) * 
              (1.0 - arma::sqrt((1.0 - 2.0*a*std::pow(neta,4.0)*B_Q) *
                (1.0 - 2.0*c*B_Q - 2.0*u*neta)));
    }
  }

  arma::vec hP = Rcpp::as<arma::vec>(h[P-1]);
  arma::cx_vec F = arma::exp(arma::log(S)*u + A_Q + B_Q*hP);

  return F;
}

需要注意的几个小改动：

我正在使用arma::函数进行向量化计算，例如arma::log、arma::exp、arma::round、arma::sqrt，以及各种重载运算符（*、+、-） ;但使用std::pow 和std::sqrt 进行标量计算。在 R 中，这对我们来说是抽象的，但在这里我们必须区分这两种情况。
您的函数 F 有一个循环 - for (i in 1:steps) - 但 C++ 版本有两个，只是由于两种语言之间循环语义的差异。
大多数输入向量是arma:: 类（而不是使用Rcpp::NumericVector 和Rcpp::ComplexVector），例外是P 和h，因为Rcpp 向量提供类似R 的元素访问 - 例如。 h[P-1]。另请注意，P 需要偏移 1（C++ 中基于 0 的索引），然后使用 Rcpp::as<arma::vec> 转换为犰狳向量 (hP)，因为如果您尝试乘以 @，编译器会报错987654348@ 和 NumericVector (B_Q*hP)。
我添加了一个函数参数h - 依赖全局变量h 的存在不是一个好主意，您在F 中这样做。如果你需要在函数体中使用它，你应该将它传递给函数。

我将您的函数名称更改为Fr，并且为了使基准测试更容易一些，我只是将您将矩阵res 填充到函数Fr 和Fcpp 中的双循环包装起来：

loop_Fr <- function(mat = res) {
  for (i in 1:N) {
    for (j in 1:Z) {
      mat[i,j]= Re(Fr(((delta*(i-1))-(alpha+1)*1i),w,b,a,c,neta,S[j],T[j],r[j],P[j],h))
    }
  }
  return(mat)
}
loop_Fcpp <- function(mat = res) {
  for (i in 1:N) {
    for (j in 1:Z) {
      mat[i,j]= Re(Fcpp(((delta*(i-1))-(alpha+1)*1i),w,b,a,c,neta,S[j],T[j],r[j],P[j],h))
    }
  }
  return(mat)
}
##
R> all.equal(loop_Fr(),loop_Fcpp())
[1] TRUE

我比较了 N = 100、N = 1000 和 N = 100000 的两个函数（这需要很长时间） - 相应地调整 lambda 和 res，但保持其他所有内容相同。一般来说，Fcpp 在我的电脑上比Fr 快大约 10 倍：

N <- 100
lambda <- (2*pi)/(N*delta)
res <- matrix(nrow=N, ncol=Z)
##
R> microbenchmark::microbenchmark(loop_Fr(), loop_Fcpp(),times=50L)
Unit: milliseconds
        expr       min        lq    median        uq       max neval
   loop_Fr() 142.44694 146.62848 148.97571 151.86318 186.67296    50
 loop_Fcpp()  14.72357  15.26384  15.58604  15.85076  20.19576    50

N <- 1000
lambda <- (2*pi)/(N*delta) 
res <- matrix(nrow=N, ncol=Z)
##
R> microbenchmark::microbenchmark(loop_Fr(), loop_Fcpp(),times=50L)
Unit: milliseconds
        expr       min        lq    median        uq       max neval
   loop_Fr() 1440.8277 1472.4429 1491.5577 1512.5636 1565.6914    50
 loop_Fcpp()  150.6538  153.2687  155.4156  158.0857  181.8452    50

N <- 100000
lambda <- (2*pi)/(N*delta)
res <- matrix(nrow=N, ncol=Z)
##
R> microbenchmark::microbenchmark(loop_Fr(), loop_Fcpp(),times=2L)
Unit: seconds
        expr       min        lq    median        uq       max neval
   loop_Fr() 150.14978 150.14978 150.33752 150.52526 150.52526     2
 loop_Fcpp()  15.49946  15.49946  15.75321  16.00696  16.00696     2

其他变量，如您的问题所示：

S <- c(906.65,906.65,906.65,906.65,906.65,906.65,906.65)
T <- c(0.1371253,0.1457896,0.1248953,0.1261278,0.1156931,0.0985253,0.1332596)
r <- c(0.013975,0.013975,0.013975,0.013975,0.013975,0.013975,0.013975)             
h <- c(0.001332596,0.001248470,0.001251458,0.001242143,0.001257921,0.001235755,0.001238440)       
P <- c(3,1,5,2,1,4,2)
w <- 0.001; b <- 0.2; a <- 0.0154; c <- 0.0000052; neta <- (-0.70)
Z <- length(S)
alpha <- 2; delta <- 0.25

【讨论】：