Straight, curly, or compiled?

September 7, 2010

(This article was first published on Thinking inside the box , and kindly contributed to R-bloggers)

Christian Robert, whose blog I
commented-on here once before, had
followed up on a recent set of posts by Radford Neal which had appeared both on Radford’s blog and on the
r-devel mailing list.

Now, let me prefix this by saying that I really enjoyed Radford’s posts. He obviously put a lot of time into finding a number of (all somewhat
small in isolation) inefficiencies in R which, when taken together, can make a difference in
performance. I already spotted one commit by Duncan in the SVN logs for R so this is being looked at.

Yet Christian, on the other hand, goes a little overboard in bemoaning performance differences somewhere between ten and fifteen percent — the
difference between curly and straight braces (as noticed in Radford’s first post). Maybe he spent too much time waiting for his MCMC runs to
finish to realize the obvious: compiled code is evidently much faster.

And before everybody goes and moans and groans that that is hard, allow me to just interject and note that it is not. It really
doesn’t have to be. Here is a quick
cleaned up version of Christian’s example code, with proper assigment operators and a second variable x. We then get to the
meat and potatoes and load our
Rcpp package as well as
inline to define the same little test function in C++. Throw in
rbenchmark which I am becoming increasingly fond of for these little timing tests,
et voila, we have ourselves a horserace:

# Xian's code, using <- for assignments and passing x down
f <- function(n, x=1) for (i in 1:n) x=1/(1+x)
g <- function(n, x=1) for (i in 1:n) x=(1/(1+x))
h <- function(n, x=1) for (i in 1:n) x=(1+x)^(-1)
j <- function(n, x=1) for (i in 1:n) x={1/{1+x}}
k <- function(n, x=1) for (i in 1:n) x=1/{1+x}

# now load some tools

# and define our version in C++
l <- cxxfunction(signature(ns="integer", xs="numeric"),
                 'int n = as(ns); double x=as(xs);
                  for (int i=0; i
                  return wrap(x); ',

# more tools

# now run the benchmark
N <- 1e6
benchmark(f(N, 1), g(N, 1), h(N, 1), j(N, 1), k(N, 1), l(N, 1),
          columns=c("test", "replications", "elapsed", "relative"),
          order="relative", replications=10)

And how does it do? Well, glad you asked. On my i7, which the other three cores standing around and watching, we get an
eighty-fold increase relative to the best interpreted version:

/tmp$ Rscript xian.R
Loading required package: methods
     test replications elapsed relative
6 l(N, 1)           10   0.122    1.000
5 k(N, 1)           10   9.880   80.984
1 f(N, 1)           10   9.978   81.787
4 j(N, 1)           10  11.293   92.566
2 g(N, 1)           10  12.027   98.582
3 h(N, 1)           10  15.372  126.000

So do we really want to spend time arguing about the ten and fifteen percent differences? Moore’s law gets you
those gains in a couple of weeks anyway. I’d much rather have a conversation about how we can get people speed increases that are orders of
magnitude, not fractions. Rcpp is one such tool. Let’s get more of them.

To leave a comment for the author, please follow the link and comment on their blog: Thinking inside the box . offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.


Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)