title: How to speed up this Rcpp function?
body: <p>I wish to implement a simple <code>split-apply-combine</code> routine in <code>Rcpp</code> where a dataset (matrix) is split up into groups, and then the groupwise column sums are returned. This is a procedure easily implemented in <code>R</code>, but often takes quite some time. I have managed to implement an <code>Rcpp</code> solution that beats the performance of <code>R</code>, but I wonder if I can further improve upon it. To illustrate, here some code, first for the use of <code>R</code>:</p>
  
  <pre><code>n &lt;- 50000
  k &lt;- 50
  set.seed(42)
  X &lt;- matrix(rnorm(n*k), nrow=n)
  g=rep(1:8,length.out=n )
  
  use.for &lt;- function(mat, ind){
    sums &lt;- matrix(NA, nrow=length(unique(ind)), ncol=ncol(mat))
    for(i in seq_along(unique(ind))){
      sums[i,] &lt;- colSums(mat[ind==i,])
    }
    return(sums)
  }
  
  use.apply &lt;- function(mat, ind){
    apply(mat,2, function(x) tapply(x, ind, sum))
  }
  
  use.dt &lt;- function(mat, ind){ # based on Roland's answer
     dt &lt;- as.data.table(mat)
     dt[, cvar := ind]
     dt2 &lt;- dt[,lapply(.SD, sum), by=cvar]
     as.matrix(dt2[,cvar:=NULL])
  }
  </code></pre>
  
  <p>It turns out that the <code>for</code>-loops is actually quite fast and is the easiest (for me) to implement with <code>Rcpp</code>. It works by creating a submatrix for each group and then calling <code>colSums</code> on the matrix. This is implemented using <code>RcppArmadillo</code>:</p>
  
  <pre><code>#include &lt;RcppArmadillo.h&gt;
  // [[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&lt;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;
  }
  </code></pre>
  
  <p>However, based on <a href="https://stackoverflow.com/questions/24933290/elementwise-matrix-multiplication-r-versus-rcpp-how-to-speed-this-code-up">answers to this question</a>, I learned that creating copies is expensive in <code>C++</code> (just as in <code>R</code>), but that loops are not as bad as they are in <code>R</code>. Since the <code>arma</code>-solution relies on creating matrixes (<code>submat</code> in the code) for each group, my guess is that avoiding this will speed up the process even further. Hence, here a second implementation based on <code>Rcpp</code> only using a loop:</p>
  
  <pre><code>#include &lt;Rcpp.h&gt;
  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&lt;gr_n; g++){
       int g_id = gr(g);
  
        for (int j = 0; j &lt; ncol; j++) {
        double total = 0;
          for (int i = 0; i &lt; 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;
  }
  </code></pre>
  
  <p>Benchmarking these solutions, including the <code>use_dt</code> version provided by @Roland (my previous version discriminted unfairly against <code>data.table</code>), as well as the <code>dplyr</code>-solution suggested by @beginneR, yields the following:</p>
  
  <pre><code> 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
  </code></pre>
  
  <p>My intution (<code>use_Rcpp</code> better than <code>use_arma</code>) did not turn out right. Having said that, I guess that the line <code>if (G(i) != g_id) continue;</code> in my <code>use_Rcpp</code> function slows down everything. I am happy to learn about alternatives to set this up. </p>
  
  <p>I am happy that I have achieved the same task in half the time it takes <code>R</code> to do it, but maybe the several <code>Rcpp is much faster than R</code>-examples have messed with my expectations, and I am wondering if I can speed this up even more. Does anyone have an idea? I also welcome any programming / coding comments in general since I am relatively new to <code>Rcpp</code> and <code>C++</code>. </p>