Benchmarking and Profiling

R Implementation

gibbs_r <- function(N, thin) {

  mat <- matrix(NA, nrow = N, ncol = 2)
  x <- 0
  y <- 0

  for(i in 1:N) {
    j <- 0
    for(j in 1:thin) {
      x <- rgamma(1, 3.0, 1.0 / (y * y + 4));
      y <- rnorm(1, 1.0 / (x + 1), 1.0 / sqrt(2 * x + 2));
    }

    mat[i, 1] <- x; # Note the difference in Indexing
    mat[i, 2] <- y;

  }

  return(mat);
}

This is a straightforward implementation of a Gibbs sampler which involves sampling from conditional distributions repeatedly to generate samples from the target distribution. Why might this perform poorly as a pure R implementation? Profile the code to see what could be going on.

C++ Implementation

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix gibbs_cpp(int N, int thin) {

  NumericMatrix mat(N, 2);
  double x = 0;
  double y = 0;

  for(int i = 0; i < N; i++) {
    for(int j = 0; j < thin; j++) {
      x = R::rgamma(3.0, 1.0 / (y * y + 4)); // Rcpp::rgamma is the vectorized version, returns a vector
      y = R::rnorm(1.0 / (x + 1), 1.0 / sqrt(2 * x + 2)); // Rcpp::rnorm is the vectorized version, returns a vector
    }

    mat(i, 0) = x;  // Note the difference in Indexing
    mat(i, 1) = y;
  }

  return(mat);
}

The C++ implementation looks almost exactly the same, except for a few differences in indexing and initialization. Now, how do they each perform? Let’s go ahead and run a benchmark comparing the two.

results <- microbenchmark(
  gibbs_r   = gibbs_r(100, 1000),
  gibbs_cpp = gibbs_cpp(100, 1000),
  times = 10
)

results

## Unit: milliseconds
##       expr       min        lq      mean    median        uq       max
##    gibbs_r 472.37406 493.12392 541.83436 499.00573 519.95170 854.55663
##  gibbs_cpp  15.23462  15.77395  16.10817  15.91789  16.04142  18.60962
##  neval cld
##     10   b
##     10  a

The difference is staggering, the C++ code ends up nearly 40x faster than the pure R implementation. This is an example of why you should avoid for loops in R.