我怎样才能加快这个Rcpp代码？

Question

我已经在R中实现了一个很长的函数。我在R中成功地进行了改进，但现在我想通过使用Rcpp包来加速它的速度。

我已经创建了以下Rcpp代码。不幸的是，与R代码运行大致相同。我想这样改进它。有没有人有关如何改进这段代码的想法？

非常感谢！

#include <math.h> 
#include <Rcpp.h> 
using namespace Rcpp; 


// [[Rcpp::export]] 
double kernelcpp(NumericVector a, NumericVector b, int N){ 
    int i; 
    double sum=0.0; 
    for (i=0;i<N;i++){ 
    if (a[i] > b[i]) 
     sum+= a[i] - b[i]; 
    else 
     sum+= b[i] - a[i]; 
    } 
    return(exp(- sum)); 
} 
// [[Rcpp::export]] 
NumericVector testFromontcpp(NumericMatrix z1, NumericMatrix z2, int Nbootstrap){ 
    // first element of TKeps = TK 
    int i,j,k,t; 
    int dim1 = z1.nrow(); 
    int dim2 = z2.nrow(); 
    double n1 = (double) dim1; 
    double n2 = (double) dim2; 
    int dimension = z1.ncol(); 
    int N = dim1 + dim2; 
    NumericVector TKeps(Nbootstrap+1); 
    Rcpp::NumericMatrix bb(N,N); 

    double cc = 1/(n1*n2*(n1+n2-2)); 
    double a = sqrt(1/(n1*n1-n1)-cc); 
    double b = - sqrt(1/(n2*n2-n2)-cc); 

    for (i=0 ; i<N ; i++){ 
    for (j=0 ; j<N ; j++){ 
    if (i != j){ 
     if (i < dim1) { 
     if (j < dim1){ 
      bb(i,j) = kernelcpp(z1(i,_),z1(j,_),dimension); 
     } else { 
      bb(i,j) = kernelcpp(z1(i,_),z2(j-dim1,_),dimension); 
     } 
     } 
     else{ 
     if (j < dim1){ 
      bb(i,j) = kernelcpp(z2(i-dim1,_),z1(j,_),dimension); 
     } else { 
      bb(i,j) = kernelcpp(z2(i-dim1,_),z2(j-dim1,_),dimension); 
     } 
     } 
    } 
    } 
    } 

    TKeps(0)=0.0; 
    for (i=0 ; i<N ; i++){ 
    for (j=0 ; j<N ; j++){ 
    if (i != j){ 
     if (i < dim1) { 
     if (j < dim1){ 
      TKeps(0) += bb(i,j)* (a*a + cc); 
     } else { 
      TKeps(0) += bb(i,j) * (a*b + cc); 
     } 
     } 
     else{ 
     if (j < dim1){ 
      TKeps(0) += bb(i,j) * (a*b + cc); 
     } else { 
      TKeps(0) += bb(i,j) * (b*b + cc); 
     } 
     } 
    } 
    } 
    } 


    for (k=1 ; k<=Nbootstrap ; k++){ 
    TKeps(k)=0; 
    int R[N]; 
    for (i = 0 ; i < N ; i++) 
     R[i] = i; 
    for (i = 0; i < N - 1 ; i++) { 
     int j = i + rand()/(RAND_MAX/(N - i) + 1); 
     t = R[j]; 
     R[j] = R[i]; 
     R[i] = t; 
    } 

    for (i=0 ; i<N ; i++){ 
     for (j=0 ; j<N ; j++){ 
     if (i != j){ 
      if (R[i] < n1) { 
      if (R[j] < n1){ 
       TKeps(k) += bb(i,j) * (a*a + cc); 
      } else { 
       TKeps(k) += bb(i,j) * (a*b + cc); 
      } 
      } else{ 
      if (R[j] < n1){ 
       TKeps(k) += bb(i,j) * (b*a + cc); 
      } else { 
       TKeps(k) += bb(i,j) * (b*b + cc); 
      } 
      } 
     } 
     } 
    } 
    } 
    return(TKeps); 
} 

Answer 1

因为我并不确切地知道你的代码做什么，我可以看到两件事情从暂存：

• 您从[R环境调用的函数是testFromontcpp（...）。我建议这个函数应该有SEXP值作为参数。那些S表达式是指向R的内存的指针。如果你不使用SEXP，那么两个矩阵都将被复制： 考虑一个1000x1000的矩阵，这意味着你有一百万个条目保存在R中，它们被复制到C++中。这样做写：

  testFromontcpp（SEXP的x，SEXP Y，SEXP Z）{

  NumericMatrix Z1（x）中，Z 2（Y）;

  int * Nbootstrap = INTEGER（z）;

  ... }

注意：在for循环中不能使用i<Nbootstrap。你必须写i<*Nbootstrap !!!

• 其次...更重要的是：由于R的矩阵被保存为指向列的指针，并且从列指针保存到行，因此C的矩阵被保存。我想说的是，跳入内存并跳回整个时间而不是跟随内存路径花费很多。我的建议是：切换for循环，所以首先遍历特定列的行，而不是相反。

最后一点：在大学的一项任务中，我也遇到了迭代矩阵的问题。在我的情况下，转置矩阵然后进行计算便宜得多。

我希望我能帮助你。

最佳， 迈克尔

PS：参照点1，......我只是基准你的代码与您的实现，并使用SEXP。使用SEXP，对于1到10之间的随机数字的100x100矩阵，其速度稍微快一些。