An accept-reject sampler using RcppArmadillo::sample()

Jonathan Olmsted — written May 8, 2013 — source

The recently added RcppArmadillo::sample() functionality provides the same algorithm used in R’s sample() to Rcpp-level code. Because R’s own sample() is written in C with minimal work done in R, writing a wrapper around RcppArmadillo::sample() to then call in R won’t get you much of a performance boost. However, if you need to repeatedly call sample(), then calling a single function which performs everything in Rcpp-land (including multiple calls to sample()) before returning to R can produce a noticeable speedup over a purely R-based solution.

Accept-Reject Sampler Example

One place where this situation arises is in an accept-reject sampler where the candidate “draw” is the output of a call to sample(). Concretely, let’s suppose we want to sample 20 integers (without replacement) from 1 to 50 such that the sum of the 20 integers is less than 400. Far fewer than 10% of randomly drawn samples will meet this constraint.

require(RcppArmadillo)

Loading required package: RcppArmadillo

Loading required package: Rcpp

require(rbenchmark)

Loading required package: rbenchmark

The R code is straightforward enough. It has been written to mirror the logic of the C++ code, although that doesn’t come at the cost of much performance.

r_getInts <- function(samples) {
    thresh <- 400
    results <- matrix(0, 20, samples) ;
    cnt <-  0

    while(cnt < samples) {
        candidate = sample(1:50, 20)

        if (sum(candidate) < thresh) {
            results[, cnt + 1] <- candidate
            cnt <- cnt + 1
        }
    }

    return(results)
}

Although it is a bit longer, the logic of the C++ code is similar.

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

using namespace Rcpp ;

// [[Rcpp::export]]
IntegerMatrix cpp_getInts(int samples
                          ) {
    int cnt = 0 ;
    IntegerMatrix results(20, samples) ;
    IntegerVector frame = seq_len(50) ;
    IntegerVector candidate(20) ;
    int thresh = 400 ;
  
    while (cnt < samples) {
        candidate = RcppArmadillo::sample(frame, 
                                          20, 
                                          FALSE, NumericVector::create()
                                          ) ;
        double sum = std::accumulate(candidate.begin(), candidate.end(), 0.0) ;
    
        if (sum < thresh) {
            results(_, cnt) = candidate ;
            cnt++ ;
        }
    }
  
    return results ;
}

Performance

The Rcpp code tends to be about 7-9 times faster and this boost increases as the constraint becomes more complicated (and necessarily more costly in R).

benchmark(r = {set.seed(1); r_getInts(50)},
          cpp = {set.seed(1); cpp_getInts(50)},
          replications = 10,
          order = 'relative',
          columns = c("test", "replications", "relative", "elapsed")
          )

  test replications relative elapsed
2  cpp           10     1.00   0.036
1    r           10    11.97   0.431

In the Real World …

Where might the structure in this problem arise in practice? One set of instances are those where “space” matters:

sampling US cities such that no more than two are in any one state
sampling cellphone towers such that no two are closer than X miles apart
sampling nodes in a graph/network such that no one has more than K edges

In these situations, R code to check the acceptance condition will likely be less efficient relative to the corresponding C++ code and so even larger speed-ups are realized.

tags: rng

Extending R with C++ and Fortran — Dirk Eddelbuettel and JBrandon Duck-Mayr
Benchmarking Rcpp code with RcppClock — Zach DeBruine
Simulation Smoother using RcppArmadillo — Tomasz Woźniak