Rcpp packages benchmarks

Sum of vector

std::accumulate

#include <Rcpp.h>

// [[Rcpp::export]]
double sumRcpp(const Rcpp::NumericVector& x) {
    return std::accumulate(x.begin(), x.end(), 0.0);
}

Loop

#include <Rcpp.h>

// [[Rcpp::export]]
double sumRcpp2(const Rcpp::NumericVector& x) {
    std::size_t n = x.size();
    double total = 0;
    if (n > 0) {
        for (std::size_t i = 0, k = n - 1; i < k; ++i, --k)
            total += x[i] + x[k];
        if (n % 2 == 1)
            total += x[n / 2];
    }
    return total;
}

RcppArmadillo

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
double sumArmadillo(const arma::vec& x) {
    return arma::accu(x);
}

RcppEigen

// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>

typedef Eigen::Map<Eigen::VectorXd> MapNumVec;

// [[Rcpp::export]]
double sumEigen(const MapNumVec& x) {
    return x.sum();
}

RcppParalell

// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
#include <Rcpp.h>

using namespace RcppParallel;

struct Sum : Worker {
    const RVector<double> input;
    double value;
    
    Sum(const Rcpp::NumericVector input) : input(input), value(0) {}
    Sum(const Sum& sum, Split) : input(sum.input), value(0) {}
    
    void operator()(std::size_t begin, std::size_t end) {
        value += std::accumulate(input.begin() + begin, input.begin() + end, 0.0);
    }
    
    void join(const Sum& rhs) {
        value += rhs.value;
    }
};

// [[Rcpp::export]]
double sumRcppParallel(const Rcpp::NumericVector& x) {
    Sum sum(x);
    parallelReduce(0, x.size(), sum);
    return sum.value;
}

Benchmark function

library(microbenchmark)
library(ggplot2)
theme_set(theme_bw())
bench <- function(..., n = 10e2 * 1:10e2, times = 100) {
    exprs <- as.list(match.call(expand.dots = FALSE)$...)
    nexprs <- length(exprs)
    nout <- nexprs * times
    out <- list(n = integer(nout),
                expr = factor(integer(nout), levels = as.character(exprs)),
                time = numeric(nout))
    for (i in seq_along(n)) {
        start <- (i - 1) * nexprs + 1
        end <- start + nexprs - 1
        x <- runif(n[i])
        res <- microbenchmark(list = exprs, times = times)
        res <- stack(lapply(split(res$time, res$expr), median.default))
        out[["n"]][start:end] <- n[i]
        out[["expr"]][start:end] <- res$ind
        out[["time"]][start:end] <- res$values
    }
    class(out) <- "data.frame"
    attr(out, "row.names") <- .set_row_names(length(out$n))
    gp <- ggplot(out, aes(x = n, y = time, linetype = expr, color = expr)) +
        geom_smooth(se = FALSE)
    print(gp)
    invisible(out)
}

Benchmark results

bench(sum(x), sumRcpp(x), sumRcpp2(x),
      sumArmadillo(x), sumEigen(x), sumRcppParallel(x))