c++ - 如何加速这个 Rcpp 函数？

Question

我希望在 Rcpp 中实现一个简单的 split-apply-combine 例程，其中数据集(矩阵)被分成几组，然后分组列总和回来。这是一个在 R 中很容易实现的过程，但通常需要相当长的时间。我已经设法实现了一个性能优于 R 的 Rcpp 解决方案，但我想知道我是否可以进一步改进它。为了说明，这里有一些代码，首先用于 R 的使用:

n <- 50000
k <- 50
set.seed(42)
X <- matrix(rnorm(n*k), nrow=n)
g=rep(1:8,length.out=n )

use.for <- function(mat, ind){
  sums <- matrix(NA, nrow=length(unique(ind)), ncol=ncol(mat))
  for(i in seq_along(unique(ind))){
    sums[i,] <- colSums(mat[ind==i,])
  }
  return(sums)
}

use.apply <- function(mat, ind){
  apply(mat,2, function(x) tapply(x, ind, sum))
}

use.dt <- function(mat, ind){ # based on Roland's answer
   dt <- as.data.table(mat)
   dt[, cvar := ind]
   dt2 <- dt[,lapply(.SD, sum), by=cvar]
   as.matrix(dt2[,cvar:=NULL])
}

事实证明，for 循环实际上非常快，而且(对我而言)是使用 Rcpp 最容易实现的。它的工作原理是为每个组创建一个子矩阵，然后在矩阵上调用 colSums。这是使用 RcppArmadillo 实现的:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
arma::mat use_arma(arma::mat X, arma::colvec G){

  arma::colvec gr = arma::unique(G);
  int gr_n = gr.n_rows;
  int ncol = X.n_cols;

  arma::mat out = zeros(gr_n, ncol); 

  for(int g=0; g<gr_n; g++){
   int g_id = gr(g);
   arma::uvec subvec = find(G==g_id);
   arma::mat submat = X.rows(subvec);
   arma::rowvec res = sum(submat,0);
   out.row(g) = res;     
  }
 return out;
}

但是，基于answers to this question ，我了解到在 C++ 中创建拷贝的成本很高(就像在 R 中一样)，但是循环并不像在 R 中那么糟糕.由于 arma 解决方案依赖于为每个组创建矩阵(代码中的 submat)，我的猜测是避免这种情况会进一步加快该过程。因此，这里是基于 Rcpp 的第二个实现，仅使用循环:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix use_Rcpp(NumericMatrix X, IntegerVector G){

  IntegerVector gr = unique(G);
  std::sort(gr.begin(), gr.end());
  int gr_n = gr.size();
  int nrow = X.nrow(), ncol = X.ncol();

  NumericMatrix out(gr_n, ncol);

  for(int g=0; g<gr_n; g++){
     int g_id = gr(g);

      for (int j = 0; j < ncol; j++) {
      double total = 0;
        for (int i = 0; i < nrow; i++) {

          if (G(i) != g_id) continue;    // not sure how else to do this
          total += X(i, j);
        }
        out(g,j) = total;
      }
  }
      return out;
}

对这些解决方案进行基准测试，包括@Roland 提供的use_dt 版本(我之前的版本不公平地歧视了data.table)，以及dplyr-@beginneR 建议的解决方案产生以下结果:

 library(rbenchmark)
 benchmark(use.for(X,g), use.apply(X,g), use.dt(X,g), use.dplyr(X,g), use_arma(X,g), use_Rcpp(X,g), 
+           columns = c("test", "replications", "elapsed", "relative"), order = "relative", replications = 1000)
             test replications elapsed relative
# 5  use_arma(X, g)         1000   29.65    1.000
# 4 use.dplyr(X, g)         1000   42.05    1.418
# 3    use.dt(X, g)         1000   56.94    1.920
# 1   use.for(X, g)         1000   60.97    2.056
# 6  use_Rcpp(X, g)         1000  113.96    3.844
# 2 use.apply(X, g)         1000  301.14   10.156

我的直觉(use_Rcpp 比 use_arma 更好)结果并不正确。话虽如此，我猜 if (G(i) != g_id) continue; 在我的 use_Rcpp 函数中会减慢一切。我很高兴了解设置它的替代方法。

我很高兴我用 R 一半的时间完成了相同的任务，但也许几个 Rcpp 比 R 快得多 -examples打乱了我的期望，我想知道我是否可以进一步加快速度。有人有想法吗？我也欢迎任何一般的编程/编码评论，因为我是 Rcpp 和 C++ 的新手。

Answer 1

不，您需要击败的不是 for 循环:

library(data.table)
#it doesn't seem fair to include calls to library in benchmarks
#you need to do that only once in your session after all

use.dt2 <- function(mat, ind){
  dt <- as.data.table(mat)
  dt[, cvar := ind]
  dt2 <- dt[,lapply(.SD, sum), by=cvar]
  as.matrix(dt2[,cvar:=NULL])
}

all.equal(use.dt(X,g), use.dt2(X,g))
#TRUE

benchmark(use.for(X,g), use.apply(X,g), use.dt(X,g), use.dt2(X,g),
          columns = c("test", "replications", "elapsed", "relative"), 
          order = "relative", replications = 50)

#             test replications elapsed relative
#4   use.dt2(X, g)           50    3.12    1.000
#1   use.for(X, g)           50    4.67    1.497
#3    use.dt(X, g)           50    7.53    2.413
#2 use.apply(X, g)           50   17.46    5.596