このRcpp関数を高速化する方法は？

データセット（行列）がグループに分割され、グループごとの列の合計が返される単純なsplit-apply-combineルーチンを実装したいと思いRcppます。これはで簡単に実装できる手順ですが、R多くの場合、かなりの時間がかかります。Rcppのパフォーマンスを上回るソリューションを実装するRことができましたが、それをさらに改善できるかどうか疑問に思います。説明のために、ここではいくつかのコードを最初に使用します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-loopsは実際には非常に高速であるとして実装するために（私にとって）最も簡単です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;
}

ただし、この質問への回答に基づいて、コピーの作成にはC++（と同じようにR）コストがかかることを学びましたが、ループはの場合ほど悪くはありませんR。以来arma-溶液は、（マトリックスの作成に依存しているsubmatグループごとにコード内に）、私の推測では、これを回避することはさらにプロセスをスピードアップすることです。したがって、ここではRcpp、ループのみの使用に基づく2番目の実装を示します。

#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;
}

use_dt@Rolandによって提供されたバージョン（以前のバージョンは不当に差別さdata.tableれていました）、dplyrおよび@beginneRによって提案された-solutionを含むこれらのソリューションをベンチマークすると、次のようになります。

 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 is much faster than R例が私の期待を台無しにしている可能性があり、これをさらにスピードアップできるかどうか疑問に思っています。誰かアイデアがありますか？また、私はへの比較的新しいですので、一般的には任意のプログラミング/コーディングコメントを歓迎RcppしてC++。