mean in RCPP slower than R mean

Question

Interested in Rcpp, I copied a simple example from Hadley Wickham's "Advanced R":

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double meanC(NumericVector x) {
  int n = x.size();
  double total = 0;

  for(int i = 0; i < n; ++i) {
    total += x[i];
  }
  return total / n;
}

/*** R
library(microbenchmark)
x <- runif(1e5)
microbenchmark(
  mean(x),
  meanC(x)
)
*/

which gives me:

Unit: microseconds
     expr     min       lq     mean   median       uq      max neval cld
  mean(x) 149.412 161.4115 181.1470 180.3395 204.2910  216.656   100  a 
 meanC(x) 394.605 400.4335 489.2311 481.6755 539.6835 1425.628   100   b

It seems that meanC() is substantially slower than mean()! Why? Can I do anything to speed up meanC?

Tested on macOS Catalina 64bit.

Answer 1

Because the code in the main loop for a (manual) mean() is so simple, optmization settings matter greatly.

If I enforce -O0 (and note that -g is also used):

R> microbenchmark(mean(x), meanC(x), meanS(x)
+ )
Unit: microseconds
     expr      min       lq     mean   median       uq     max neval cld
  mean(x)  653.089  654.093  693.971  670.952  708.419 1090.22   100 a  
 meanC(x) 1922.536 1951.835 2067.521 1980.786 2058.981 3078.64   100  b 
 meanS(x) 3409.202 3467.219 3660.131 3520.522 3618.264 5999.65   100   c
R> 

If I use -O1, or the -O3 default value I commonly use I get essentially identical results. Here is -O3:

R> microbenchmark(mean(x), meanC(x), meanS(x)
+ )
Unit: microseconds
     expr     min      lq    mean  median      uq      max neval cld
  mean(x) 653.006 653.400 683.852 668.616 699.988  869.978   100   b
 meanC(x) 435.107 435.435 460.909 438.860 465.111 1078.962   100  a 
 meanS(x) 652.505 652.873 689.620 660.695 693.213 1270.513   100   b
R> 

If I try -O6 -march=native I get about he same. There is not too much one can do, and the compiler apparently is good enough to add something worthwhile even at easiest settings.

Code below

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double meanC(NumericVector x) {
  int n = x.size();
  double total = 0;

  for(int i = 0; i < n; ++i) {
    total += x[i];
  }
  return total / n;
}

// [[Rcpp::export]]
double meanS(const Rcpp::NumericVector& x) {
  return Rcpp::mean(x);
}

/*** R
library(microbenchmark)
x <- runif(5e5)
microbenchmark(mean(x), meanC(x), meanS(x)
)
*/