Как ускорить эту функцию 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 только используя цикл:

#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, дает следующее:

 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++,