Performance drop when including code exported from other Rcpp package

Viewed 238

I have recently created a package and would like to recycle many of the "under-the-hood" functions I wrote for it in a new package. However, in a first try I have found that there's a significant performance decrease when importing cpp code into the new package. I will clarify below.

I have package1, created via RcppArmadillo::RcppArmadillo.package.skeleton(). The package's only source file is package1/src/shared.cpp, which includes a function that calculates column sums of a matrix using RcppArmadillo. The source code for shared.cpp is thus:

//[[Rcpp::depends(RcppArmadillo)]]
//[[Rcpp::interfaces(r, cpp)]]

#include "RcppArmadillo.h"

// [[Rcpp::export]]
arma::vec col_sums(const arma::mat& matty){
  return arma::sum(matty, 0).t();
}

Now suppose I want to recycle this function in another package, called package2. I do so by editing Imports and LinkingTo in DESCRIPTION, adding package1. Then, the only source file for this new package is package2/src/testimport.cpp

//[[Rcpp::depends(RcppArmadillo, package1)]]

#include "RcppArmadillo.h"
#include "package1.h"

//[[Rcpp::export]]
arma::vec col_sums(const arma::mat& test){
  return arma::sum(test,0).t();
}

//[[Rcpp::export]]
arma::vec col_sums_imported(const arma::mat& test){
  return package1::col_sums(test);
}

Now if I compile both packages, and benchmark the 3+1 functions, I get

library(magrittr)
library(rbenchmark)

nr <- 100
p <- 800

testmat <- rnorm(nr * p) %>% matrix(ncol=p)

benchmark(package2::col_sums(testmat),
          package2::col_sums_imported(testmat), 
          colSums(testmat),
          package1::col_sums(testmat),
          replications=1000)

I would expect there to be no difference at all between package1::col_sums and package2::col_sums, but also a minimal or small difference between those two and package2::col_sums_imported, which calls package1::col_sums from package2 using the cpp interface.

Instead I get (I also added R's colSums for a comparison)

                                  test replications elapsed relative user.self sys.self user.child sys.child
3                     colSums(testmat)         1000   0.050    1.429     0.052    0.000          0         0
4          package1::col_sums(testmat)         1000   0.035    1.000     0.036    0.000          0         0
1        package2::col_sums(testmat)         1000   0.038    1.086     0.036    0.000          0         0
2 package2::col_sums_imported(testmat)         1000   0.214    6.114     0.100    0.108          0         0

This 6x slow-down puzzles me, because I was not expecting there to be such a difference. Is it preferrable to just copy the source of the "shared" functions to the new package, and why? I feel like having only one source for col_sums would allow me to more easily propagate changes across both packages. Or is there another reason why my code slows down so much?

EDIT: in addition to @duckmayr's answer below, I have updated my minimal github package example to show how a user-created function can be used in package1, exported to other packages, imported to package2. Code can be found at https://github.com/mkln/rcppeztest

2 Answers

As mentioned by others, allowing other packages to call your C++ code from C++ requires the use of header files in inst/include/. Rcpp::interfaces allows you to automate the creation of such files. However, as I demonstrate below, creating your own headers manually can result in faster execution time. I believe it is because relying on Rcpp::interfaces to create your headers for you may result in more complicated header code.

Before I go further and demonstrate a "simpler" approach that results in faster execution time, I need to note that while this works for me (and I have used the approach I will demonstrate below several times without issue), the more "complex" approach taken by Rcpp::interfaces is used in part to comport with statements in Section 5.4.3. of the Writing R Extensions manual. (Specifically, the bits having to do with R_GetCCallable you'll see below). So, improve your execution time with the code I offer below at your own peril.1,2

A simple header to share the code for col_sums might look like this:

#ifndef RCPP_package3
#define RCPP_package3

#include <RcppArmadillo.h>

namespace package3 {
    inline arma::vec col_sums(const arma::mat& test){
      return arma::sum(test,0).t();
    }
}

#endif

However, the header created by Rcpp::interfaces looks like this:

// Generated by using Rcpp::compileAttributes() -> do not edit by hand
// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#ifndef RCPP_package1_RCPPEXPORTS_H_GEN_
#define RCPP_package1_RCPPEXPORTS_H_GEN_

#include <RcppArmadillo.h>
#include <Rcpp.h>

namespace package1 {

    using namespace Rcpp;

    namespace {
        void validateSignature(const char* sig) {
            Rcpp::Function require = Rcpp::Environment::base_env()["require"];
            require("package1", Rcpp::Named("quietly") = true);
            typedef int(*Ptr_validate)(const char*);
            static Ptr_validate p_validate = (Ptr_validate)
                R_GetCCallable("package1", "_package1_RcppExport_validate");
            if (!p_validate(sig)) {
                throw Rcpp::function_not_exported(
                    "C++ function with signature '" + std::string(sig) + "' not found in package1");
            }
        }
    }

    inline arma::vec col_sums(const arma::mat& matty) {
        typedef SEXP(*Ptr_col_sums)(SEXP);
        static Ptr_col_sums p_col_sums = NULL;
        if (p_col_sums == NULL) {
            validateSignature("arma::vec(*col_sums)(const arma::mat&)");
            p_col_sums = (Ptr_col_sums)R_GetCCallable("package1", "_package1_col_sums");
        }
        RObject rcpp_result_gen;
        {
            RNGScope RCPP_rngScope_gen;
            rcpp_result_gen = p_col_sums(Shield<SEXP>(Rcpp::wrap(matty)));
        }
        if (rcpp_result_gen.inherits("interrupted-error"))
            throw Rcpp::internal::InterruptedException();
        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
            throw Rcpp::LongjumpException(rcpp_result_gen);
        if (rcpp_result_gen.inherits("try-error"))
            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
        return Rcpp::as<arma::vec >(rcpp_result_gen);
    }

}

#endif // RCPP_package1_RCPPEXPORTS_H_GEN_

So, I created two additional packages via

library(RcppArmadillo)
RcppArmadillo.package.skeleton(name = "package3", example_code = FALSE)
RcppArmadillo.package.skeleton(name = "package4", example_code = FALSE)

Then in package3/inst/include, I added package3.h containing the "simple header" code above (I also added a throwaway "Hello World" cpp file in src/). In package4/src/ I added the following:

#include <package3.h>

// [[Rcpp::export]]
arma::vec col_sums(const arma::mat& test){
  return arma::sum(test,0).t();
}

// [[Rcpp::export]]
arma::vec simple_header_import(const arma::mat& test){
  return package3::col_sums(test);
}

as well as adding package3 to LinkingTo in the DESCRIPTION file.

Then, after installing the new packages, I benchmarked all the functions against each other:

library(rbenchmark)

set.seed(1)
nr <- 100
p <- 800
testmat <- matrix(rnorm(nr * p), ncol = p)

benchmark(original = package1::col_sums(testmat),
          first_copy = package2::col_sums(testmat),
          complicated_import = package2::col_sums_imported(testmat),
          second_copy = package4::col_sums(testmat),
          simple_import = package4::simple_header_import(testmat),
          replications = 1e3,
          columns = c("test", "relative", "elapsed", "user.self", "sys.self"),
          order = "relative")


                test relative elapsed user.self sys.self
2         first_copy    1.000   0.174     0.174    0.000
4        second_copy    1.000   0.174     0.173    0.000
5      simple_import    1.000   0.174     0.174    0.000
1           original    1.126   0.196     0.197    0.000
3 complicated_import    6.690   1.164     0.544    0.613

While the more "complicated" header function was 6x slower, the "simpler" one was not.


1. However, the automated code generated by Rcpp::interfaces does in fact include some features that may be superfluous for you beside the R_GetCCallable issue, though they may be desirable and in some other contexts necessary.

2. Registering functions is always portable, and package authors are instructed to do so by the Writing R Extensions manual, but for internal/organizational/etc. use I believe the approach featured here should work if all packages involved are built from source. See this section of Hadley Wickham's R Packages for some discussion as well as the section of the Writing R Extensions manual linked above.

Three things come to my mind:

  1. Does rbenchmark do a "warm up" cycle? If not, then the first invocation of package1::col_sums is paying the price of calling an R function. That might account for the 0.1s of system time.

  2. The function returns an Armadillo object. But when called via R, this has to be converted to an R object and back. I am not sure how light-weight these conversions are, or if data copies are made in (some) cases.

  3. The function might be to simple. Execution time is about 36 µs per function call. It seems plausible that doing this via R adds some significant overhead.

Overall, if you want to share such short running functions, you should convert them into "header only" and place them in inst/include/, as suggested by F. Privé in the comments. However, you will only share the source code but not the object code this way, i.e. package2 will have to be recompiled when the function in package1 changes.

I was curious to see how efficient it is to call a function exported via R. I therefore added a simple test function to a sample package:

//[[Rcpp::interfaces(r, cpp)]]
#include <thread>
#include <chrono>

#include <Rcpp.h>

// [[Rcpp::export]]
int mysleep(int msec) {
  std::this_thread::sleep_for (std::chrono::microseconds(msec));
  return msec;
}

I then compared calling this function either directly or indirectly as exported function for sleep times 50, 500 and 5000 µs. The median execution times as reported by bench::mark:

            50µs  500µs     5ms   mem_alloc
direct     153µs  688µs  5.37ms      2.47KB 
indirect   163µs  705µs  5.39ms      4.95KB

To me this looks as if calling such a simple function indirectly adds only a few 10µs of overhead on this rather slow machine. However, we already see that the amount of memory allocated is doubled. If we look at your function that returns a more complicated structure, we get:

  expression   min  mean median      max `itr/sec` mem_alloc  n_gc n_itr
  <chr>      <bch> <bch> <bch:> <bch:tm>     <dbl> <bch:byt> <dbl> <int>
1 direct     141µs 148µs  145µs 830.14µs     6737.    10.4KB     0  3342
2 imported   344µs 703µs  832µs   1.17ms     1423.   644.2KB     7   628 

The amount of memory allocated in the indirect call is more than 60 times greater! To me this explains the performance drop.

Related