Rcpp Function slower than Rf_eval

Question

I have been working on a package that uses Rcpp to apply arbitrary R code over a group of large medical imaging files. I noticed that my Rcpp implementation is considerably slower than the original pure C version. I traced the difference to calling a function via Function, vs the original Rf_eval. My question is why is there a close to 4x performance degradation, and is there a way to speed up the function call to be closer in performance to Rf_eval?

Example:

library(Rcpp)                                                                                                                                                          
library(inline)                                                                                                                                                        
library(microbenchmark)                                                                                                                                                

cpp_fun1 <-                                                                                                                                                            
  '                                                                                                                                                                    
Rcpp::List lots_of_calls(Function fun, NumericVector vec){                                                                                                             
  Rcpp::List output(1000);                                                                                                                                             
  for(int i = 0; i < 1000; ++i){                                                                                                                                       
    output[i] = fun(NumericVector(vec));                                                                                                                               
  }                                                                                                                                                                    
  return output;                                                                                                                                                       
}                                                                                                                                                                      
'                                                                                                                                                                      

cpp_fun2 <-                                                                                                                                                            
  '                                                                                                                                                                    
Rcpp::List lots_of_calls2(SEXP fun, SEXP env){                                                                                                                         
  Rcpp::List output(1000);                                                                                                                                             
  for(int i = 0; i < 1000; ++i){                                                                                                                                       
    output[i] = Rf_eval(fun, env);                                                                                                                                     
  }                                                                                                                                                                    
  return output;                                                                                                                                                       
}                                                                                                                                                                      
'                                                                                                                                                                      

lots_of_calls <- cppFunction(cpp_fun1)                                                                                                                                 
lots_of_calls2 <- cppFunction(cpp_fun2)                                                                                                                                

microbenchmark(lots_of_calls(mean, 1:1000),                                                                                                                            
               lots_of_calls2(quote(mean(1:1000)), .GlobalEnv))

Results

Unit: milliseconds
                                            expr      min       lq     mean   median       uq      max neval
                     lots_of_calls(mean, 1:1000) 38.23032 38.80177 40.84901 39.29197 41.62786 54.07380   100
 lots_of_calls2(quote(mean(1:1000)), .GlobalEnv) 10.53133 10.71938 11.08735 10.83436 11.03759 18.08466   100

Answer 1

Rcpp is great because it makes things look absurdly clean to the programmer. The cleanliness has a cost in the form of templated responses and a set of assumptions that weigh down the execution time. But, such is the case with a generalized vs. specific code setup.

Take for instance the call route for an Rcpp::Function. The initial construction and then outside call to a modified version of Rf_reval requires a special Rcpp specific eval function given in Rcpp_eval.h. In turn, this function is wrapped in protections to protect against a function error when calling into R via a Shield associated with it. And so on...

In comparison, Rf_eval has neither. If it fails, you will be up the creek without a paddle. (Unless, of course, you implement error catching via R_tryEval for it.)

With this being said, the best way to speed up the calculation is to simply write everything necessary for the computation in C++.

Answer 2

Besides the points made by @coatless, you aren't even comparing apples with apples. Your Rf_eval example does not pass the vector to the function, and, more importantly, plays tricks on the function via quote().

In short, it is all a little silly.

Below is a more complete example using the sugar function mean().

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List callFun(Function fun, NumericVector vec) {
  List output(1000);
  for(int i = 0; i < 1000; ++i){
    output[i] = fun(NumericVector(vec));
  }
  return output;
}

// [[Rcpp::export]]
List callRfEval(SEXP fun, SEXP env){
  List output(1000);
  for(int i = 0; i < 1000; ++i){
    output[i] = Rf_eval(fun, env);
  }
  return output;
}

// [[Rcpp::export]]
List callSugar(NumericVector vec) {
  List output(1000);
  for(int i = 0; i < 1000; ++i){
    double d = mean(vec);
    output[i] = d;
  }
  return output;
}

/*** R
library(microbenchmark)
microbenchmark(callFun(mean, 1:1000),   
               callRfEval(quote(mean(1:1000)), .GlobalEnv),
               callSugar(1:1000))
*/

You can just sourceCpp() this:

R> sourceCpp("/tmp/ch.cpp")

R> library(microbenchmark)

R> microbenchmark(callFun(mean, 1:1000), 
+                callRfEval(quote(mean(1:1000)), .GlobalEnv),
+                callSugar(1:1000))
Unit: milliseconds
                                        expr      min       lq     mean   median       uq       max neval
                       callFun(mean, 1:1000) 14.87451 15.54385 18.57635 17.78990 18.29127 114.77153   100
 callRfEval(quote(mean(1:1000)), .GlobalEnv)  3.35954  3.57554  3.97380  3.75122  4.16450   6.29339   100
                           callSugar(1:1000)  1.50061  1.50827  1.62204  1.51518  1.76683   1.84513   100
R>