FastR

Introduction

This document provides ways to speed up your code before parallelization on your own machine or cluster. There are many ways to use R in an optimal fashion, but often these are not taught in applied statistics courses or perhaps are not otherwise presented in a cohesive way. In other cases, one may simply be new to R, or coming at it from a programming language that works differently.

The ‘R is slow’ mantra is an old one, and has never been entirely accurate in my opinion for two reasons. First, even if it is slower, and it is relative to many languages for many situations, the time saved in programming will far exceed any computational time lost. R often can do things in 1 or 2 lines of code that will require custom functions to be written in other languages, and in that sense R is one of the fastest languages out there for statistical analysis. This is usually the case for data processing/manipulation as well. R was designed expressly for statistical analysis, not for computational speed. For more on some of this point see Hadley Wickham’s chapter.

Secondly, there are simply faster ways to use R, and they often involve a programming approach that is different from what one would use in say, Python or other languages. That is the focus of this document.

For the timings we’ll be using the microbenchmark package, which can show speed differences for even extremely fast operations. This way we don’t have to engage in necessarily slow operations to show the differences well. You’ll need to install it if you want to make your own comparisons. Note also that timings will be in microseconds, milliseconds, or seconds.

• To be added
  • data structures/representations (alternatives to data.frames)
    • data.table, dplyr, dtplyr, multidplyr
    • pooling
    • sparklyr

Color coding:

emphasis

package

function

object/class

link

Faster Data Input/Output

The base R functions read.csv and similar work fine for smaller data sets. Given their flexibility they may be preferred when useful. However with larger data sets they can be too slow and/or use too much memory. One has plenty of alternatives though.

readr

readr is a newer package and provides many complements of the base R functions, substituting an _ for the . of commonly used file reading functions. For example:

• read_delim: read delimited files
  • read_csv: comma delimited
  • read_tsv: tab delimited
• read_lines: read a file line by line
• read_file: read a whole file in as a string (e.g. for text processing)

To show the difference, we have a simple data set of two columns, one numeric (random normal), and one string (letters). It is large enough for demonstration, but still not very big (less than 6MB on disk).

'data.frame':   260000 obs. of  2 variables:
 $ nums: num  -0.157 -1.935 2.439 0.127 -0.233 ...
 $ lets: Factor w/ 26 levels "a","b","c","d",..: 1 1 1 1 1 1 1 1 1 1 ...

Here is a comparison reading the the file in with both readr and read.csv.

Unit: milliseconds
  expr      min       lq     mean    median        uq       max neval
 readr  71.2180  71.8998  72.6430  72.10712  72.65593  81.51621    20
 baser 335.3346 336.4661 341.5536 341.49618 344.54989 353.95927    20

The readr package makes assumptions about the column type based on the first 1000 observations, but one can override this. As mentioned, for small data the standard function is fine, and even faster by default.

Unit: milliseconds
  expr      min       lq     mean   median       uq      max neval
 readr 3.651792 3.740206 3.858049 3.798122 3.844015 8.953091   100
 baser 1.211106 1.250107 1.328205 1.272541 1.322686 5.587508   100

Chunked

The readr package has chunked versions of those same functions, so that not all the data has to be read in at the same time. This is extremely important for data processing. While one may not have as much choice once the analysis stage is reached, many data processing tasks can be chunked or be done line by line, thereby opening doors for parallelization and memory control.

Unit: milliseconds
  expr      min       lq     mean   median       uq      max neval
 readr 710.0352 722.2291 737.2585 724.8806 727.8937 850.4921    20
 chunk 695.6169 707.1263 726.6579 711.2742 763.2438 769.3471    20

data.table

The data.table package also has a read function that will amount to faster read speeds. The function fread works much like the previous functions to detect column types quickly and otherwise has optimizations in place to speed the process. Both read_* and fread might be viable depending on the situation. Wickham notes the following comparison:

• Compared to fread, readr:

• Is slower (currently ~1.2-2x slower. If you want absolutely the best performance, use data.table::fread().

• Readr has a slightly more sophisticated parser, recognising both doubled (“”“”) and backslash escapes (“"”). Readr allows you to read factors and date times directly from disk.

• fread() saves you work by automatically guessing the delimiter, whether or not the file has a header, how many lines to skip by default and more. Readr forces you to supply these parameters.

• The underlying designs are quite different. Readr is designed to be general, and dealing with new types of rectangular data just requires implementing a new tokenizer. fread() is designed to be as fast as possible. fread() is pure C, readr is C++ (and Rcpp).

I will note that readr was faster on the 54 Mb chunk test file, so you may need very large files before fread starts to shine.

Lines vs. Whole Data

Many do not realize that they don’t have to read in the whole data file. Many preprocessing steps can be done line by line, and this may save on memory and provide a means of (massively) parallel operations. Also, and here is a point further not considered by many, this provides a means to debug your code before going prime time. Your file may be huge, but you don’t have to work with the whole thing to get your code sorted out. See readLines in base R for details, but you’d want to use read_lines in readr.

Other

feather

Yet another alternative is to use a different data format altogether. Wickham and Wes McKinney (the latter is behind Python’s pandas module) are working on feather, a binary data format that would be read in R or Python, and allow for very fast read/write operations, as well as smaller file sizes, and more interoperability between the two languages. See this link for more details. The following compares times using feather vs. the readr counterparts.

Unit: milliseconds
          expr       min        lq      mean    median        uq      max neval
     write_csv 182.09487 193.72603 223.05187 214.96862 239.24572 350.7905    20
 write_feather  30.91076  31.46558  39.45274  32.33418  51.06321  64.0367    20
Unit: microseconds
         expr       min        lq      mean    median        uq       max neval
     read_csv 71054.955 71386.177 72744.536 71711.386 72984.514 77582.317    20
 read_feather   621.975  1192.631  2196.104  1469.016  1590.713  8013.532    20

Loops and Alternatives

Loops

Before getting to alternative methods, we’ll first do things with a standard loop, and then show others options later. Loops in and of themselves are not a problem in R, nor are they necessarily slow for some tasks. It’s just that they just often have more efficient alternatives, whether it be more succinct code, or possibly vastly quicker computation.

As one of our working examples, we’ll start with a function that does a random walk. The code is based on this document, co-authored by one of the R’s in R, with some slight additions and modifications.

rw2d_loop_nopre <- function(n) {
  xpos = ypos = 0
  xdir = c(T, F)
  pm1 = c(1,-1)
  
  for (i in 2:n)
  {
    if (sample(xdir, 1)) {
      xpos[i] = xpos[i-1] + sample(pm1, 1)
      ypos[i] = ypos[i-1]
    }
    else {
      xpos[i] = xpos[i-1]
      ypos[i] = ypos[i-1] + sample(pm1, 1)
    }
  }
  data.frame(x=xpos, y=ypos) 
}
# rw2d_loop_nopre(1e5) # for profiling
# ggplot2::qplot(x, y, data=rw2d_loop_nopre(1000), geom='path') # plot check

For better understanding of what’s going on, you can use the debug functionality on the functions used to see what they are doing line by line. Rstudio makes this especially easy, and to get started one can use code that looks like the following.

debugonce(rw2d_loop_nopre)
rw2d_loop_nopre(1000)

You will be able to see what’s produced and even manipulate those objects at each step.

Pre-allocation

One thing you’ll want to note, if you must use an explicit loop, first create an empty vessel to contain the elements of the loops output, rather than taking an append-as-you-go approach. The reasoning is similar to the reasoning behind vectorization noted later. By not having R do some (internal) operations every single iteration, you save time, and it adds up. Memory is also already made available to the object.

In the following we’ll compare the random walk both with and without pre-allocation. I have a line for each function so you can use them for profiling, to see where bottlenecks arise specifically. In (newer versions of) RStudio, you can simply highlight that line and hit Ctrl+Alt+Shft+P to do so. Also, for the first demonstration above I showed an example graph, but all others will look like some (random) version of that, so are subsequently commented out in the code.

# preallocation
rw2d_loop_pre <- function(n) {
  xpos = ypos = numeric(n)
  xdir = c(T, F)
  pm1 = c(1,-1)
  
  for (i in 2:n)
  {
    if (sample(xdir, 1)) {
      xpos[i] = xpos[i-1] + sample(pm1, 1)
      ypos[i] = ypos[i-1]
    }
    else {
      xpos[i] = xpos[i-1]
      ypos[i] = ypos[i-1] + sample(pm1, 1)
    }
  }
  data.frame(x=xpos, y=ypos)
}
# rw2d_loop_pre(1e5) # for profiling
# ggplot2::qplot(x, y, data=rw2d_loop_pre(1000), geom='path')

Unit: milliseconds
                  expr      min       lq     mean    median       uq      max neval
 rw2d_loop_nopre(1000) 9.336365 9.562017 12.40965 14.002640 14.19164 16.04319   100
   rw2d_loop_pre(1000) 7.402994 7.582314  8.90173  7.682311 12.02000 14.21070   100

As we can see, even with this quick example pre-allocation is already saving some time.

Smaller code != faster

Note that the goal is fast running code, but this doesn’t mean there is a one-to-one correlation between speed and code size. The following is a recursive version of the above and it’s three lines of code. It is a little faster than the non-pre-allocated loop, but slower than one with pre-allocation. I might also mention that R is not really that great for recursive functions in my experience, often ‘nesting too deeply’ where other languages have no issue.

rw2d2_recursive <- function(n, x=0, y=0) {
  if(n == 1) return(cumsum(data.frame(x=x, y=y)))
  step = sample(1:4, 1)
  rw2d2_recursive(n=n-1, 
                  x=c(x, c(-1, 1, 0, 0)[step]), 
                  y=c(y, c(0, 0, -1, 1)[step]))
}
# cannot be profiled
# ggplot2::qplot(x, y, data=rw2d2_recursive(1000), geom='path')

Unit: milliseconds
                  expr      min        lq     mean    median       uq      max neval
 rw2d_loop_nopre(1000) 8.816146 11.444509 13.67076 13.834903 14.90716 47.43077   100
   rw2d_loop_pre(1000) 7.135847  7.789052 10.40728  8.203848 10.92444 66.63750   100
 rw2d2_recursive(1000) 7.963387  8.481699 11.51162 10.673126 12.82716 69.51776   100

Looks can be deceiving!

Double loops

Double loops for independent operations are almost always an inefficient way to use R, especially for common data processing tasks, so if you find yourself using one, you should probably think very hard about what alternatives might be possible. In this example we create a distance matrix based on the Euclidean distance between rows in a given matrix m, which for timing demonstration will be a 100 x 100 matrix.

doubleLoopDist <- function(m) {
  n = nrow(m)
  dists = matrix(0, ncol=n, nrow=n)
  for (i in 1:n){
    for (j in 1:n){
      dists[i,j] = sqrt(sum((m[i,] - m[j,])^2))
    }
  }
  dists
}
d = matrix(rnorm(10000), ncol=100)
# doubleLoopDist(d)

As we’ll see in a moment, this is very inefficient, and I would again say there is rarely a need to a double loop. Even if you have to double loop, there might be an efficiency gain possible. Since the distance matrix is symmetric, we only need to compute half of it, then fill in the rest. This would at least cut the time in half.

doubleLoopDist2 <- function(m) {
  n = nrow(m)
  dists = matrix(0, ncol=n, nrow=n)
  for (i in 1:(n-1)){
    for (j in (i+1):n){
      dists[i,j] = sqrt(sum((m[i,] - m[j,])^2))
    }
  }
  
  dists = dists + t(dists) + diag(diag(dists))
  diag(dists) = 0
  dists
}
identical(doubleLoopDist(d), doubleLoopDist2(d))

[1] TRUE

Unit: milliseconds
   expr      min       lq     mean   median       uq      max neval
 dloop1 34.55346 37.79764 43.06150 45.81337 46.72830 98.36466   100
 dloop2 17.34432 18.49033 21.44364 19.33913 26.85913 29.12900   100

Later on we’ll see just how slow this is when we compare it to the base R dist function, which uses compiled C code.

Apply

A set of functions exists in R that typically make explicit loops unnecessary, and can reduce programming time, and potentially make code cleaner. Furthermore, they are parallelizable, and so potentially far faster. There is a misconception that these functions are faster than loops out-of-the-box, but this is incorrect, and they technically can even be slower. Again though, you automatically gain in coding efficiency and possibly gain a lot in computation time with the use of parallelization.

apply functions

• apply: arrays, matrices, data.frames
• lapply, sapply, vapply: lists, data.frames, vectors
• tapply: grouped operations (table apply)
• mapply: multivariate version of sapply
• replicate: similar to sapply

As a simple example, let’s estimate some column means. An explicit loop might look like the following.

means_loop = rep(0, ncol(d))  # initialize the means vector
for (n in 1:ncol(d)){
  means_loop[n] = mean(d[,n])
}
head(means_loop)

[1]  0.09693607  0.07491450  0.14641378  0.04092882  0.20245490 -0.09174778

Let’s use the apply function instead. The primary arguments are a matrix/data.frame (or any array), the margin to iterate over (e.g. rows or columns), and the function we want to apply, which in this case is the mean.

means_apply = apply(d, 2, mean)
identical(means_loop, means_apply)

[1] TRUE

And that’s all there is to it. A point of emphasis however: there is no guaranteed speed gain here, only more succinct code. The following demonstrates using a loop vs. apply on a 1000x1000 matrix. For a simple operation like this, the explicit loop is actually faster. However, most would not want to write a loop every time they want column or row means, and as we’ll see later, we’d use a vectorized approach rather than either of these.

Unit: milliseconds
        expr      min       lq     mean   median       uq      max neval
  means_loop 12.22220 12.55708 14.87569 12.95751 16.53907 69.19196   100
 means_apply 14.05909 15.30876 18.84277 19.11773 19.59968 71.32181   100

plyr alternates

The plyr package provides alternatives to the apply family of functions. These plyr functions work based on the first two letters. The first denotes the input, the second the output. For example, ldply takes a list and converts the output to a data.frame.

• Lists: ldply, llply, laply

• Arrays: a*ply

• Dataframes: d*ply

• Multiple arguments: m*ply

As an example will use adply on the column mean example. The d object is a matrix, and we’ll get a data.frame as output.

library(plyr)
adply(d, 2, mean)

Unfortunately, the flexibility does come at a price, and plyr functions are very slow (mostly because data.frames are slow). Using the 1000x1000 matrix again, and using aaply to be more comparable, we see the plyr function is much slower.

Unit: milliseconds
        expr      min       lq     mean   median       uq       max neval
  means_loop 12.25328 13.22891 20.46052 15.98190 16.15375  72.28072   100
 means_apply 16.02179 19.51361 30.64961 22.77069 26.35973  78.05620   100
       aaply 73.72027 78.07908 84.69652 79.26012 80.39513 134.93946   100

However, there is comparable functionality for data.frames in dplyr that would be much faster.

Parallel versions

Another thing that puts the apply approach ahead of loops regards parallelization. Most of the base R apply functions have par* counterparts, e.g. parApply, parLapply etc. This is also the case for plyr functions, which have a .parrallel=TRUE argument. Example code would like the following.

x = matrix(rnorm(1e7), ncol=1e6)
library(parallel)
cl = makeCluster(8)
clusterExport(cl, 'x')
means = parApply(cl, x, 2, mean)
stopCluster(cl)

For the plyr functions, one needs to use the foreach package and register the cluster, and supply the .parallel argument. It’s still exceedingly slow, so you might just stick with the parallel/foreach package.

library(foreach); library(doParallel)
library(parallel)
cl = makeCluster(12)
registerDoParallel(cl)
means = aaply(x, 2, mean, .parallel = T, .paropts = list(.export='x'))
stopCluster(cl)

For something simple like this, there isn’t going to be much gained in going parallel. And this is an important point to note: using more cores does not necessarily mean more speed, and this is especially the case when spreading operations across machines. Had I used 12 instead of 8 on my machine, the difference would have been negligible. I should also note that none of these approaches is nearly as fast as we could get with a built-in R function, colMeans.

Unit: seconds
          expr      min       lq     mean   median       uq      max neval
    means_loop 5.385455 5.446444 5.911896 5.624648 6.656046 6.671387    10
   means_apply 5.430653 5.491898 5.676836 5.575201 5.933212 6.021431    10
 parallelMeans 2.990517 3.152247 3.215829 3.222990 3.260839 3.424342    10

Using faster languages

Base R internals

Much of base R is not written in R, but much faster languages such as C, C++ and Fortran, and using those functions can lead to enormous speed gains. Relative to an explicit loop that does the same thing, there is no comparison to a base R function that’s calling internal C or other code. Not only will the built-in function run faster, but your code will be more succinct because you’re making a single function call as opposed to an explicit loop.

As an example let’s compare the previous double loop approach to the base R function dist.

Unit: milliseconds
   expr       min        lq      mean    median       uq       max neval
  dloop 36.542256 38.517853 43.065587 38.752303 47.46347 58.667367   100
 dloop2 17.985360 19.325787 20.986236 19.489124 19.65276 38.788812   100
   dist  1.523706  1.590126  1.614325  1.614613  1.63382  1.708158   100

If there are actually only two columns though, dist would be slower due to some preliminary checks and calculations that sum and sqrt functions do not have to make (note that these timings are in microseconds).

Unit: microseconds
    expr    min     lq     mean median     uq     max neval
    dist 65.687 67.740 72.81621 69.793 75.218 112.314   100
 sqrtsum  2.639  3.226  3.90039  3.666  4.106  17.301   100

Again, when there isn’t much to optimize, don’t waste time figuring out how to do so.

Premature optimization is the root of all evil.

~ Donald Knuth

As another example, let’s use the colMeans function on that large 10 x 1000000 matrix we used before.

Unit: milliseconds
          expr         min          lq        mean      median          uq        max neval
    means_loop 5385.455354 5446.444056 5911.895932 5624.647995 6656.045691 6671.38714    10
   means_apply 5430.653350 5491.897762 5676.835830 5575.201307 5933.212066 6021.43071    10
 parallelMeans 2990.517013 3152.246836 3215.829386 3222.989777 3260.838736 3424.34182    10
      colMeans    7.775123    7.855472    8.202892    7.900339    8.116022   12.03423   100

This shows that using the compiled C in the colMeans function is way faster than even parallelizing (in this case).

As a further illustration, let’s center a matrix, i.e. subtract the mean. We now have several options at this point, but we’ll keep it to the following:

• using the scale function with the argument scale = FALSE
• using the sweep function
• subtract the column means calculated via colMeans

Unit: milliseconds
              expr      min       lq      mean    median       uq      max neval
       means_sweep 89.46845 91.09817 119.72156 104.24922 167.7151 181.0814   100
             scale 95.17531 97.14974 119.93443 111.44093 121.7853 209.0363   100
 subtract_colMeans 33.78692 34.21637  49.67427  40.48495  54.2156 117.4253   100

The interesting thing to note here is that sweep is already as fast as the built in scale function, but we could have provided any statistic of the columns to ‘sweep out’. For example, maybe we prefer the median, or perhaps we want to add rather than subtract. The slowdown of sweep is that the statistic to be swept out may require a slower and/or separate operation (e.g. I could first create an object means that are the column means). In the end, one might have several options depending on their goals.

Rcpp

We can potentially take advantage of the power of C++ explicitly within R, by using the Rcpp package. It’s also worth noting that many R packages do so to pick up some speed.

library(Rcpp)
distC = cppFunction('
NumericMatrix distC(NumericMatrix x) {
  int n = x.nrow();
  NumericMatrix out(n, n);
  for(int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j) {
      out(i,j) = sqrt(sum(pow(x(i,_) - x(j,_), 2.0)));
    }
  }
  return out;
}
')
distC(d[1:3,])

         [,1]     [,2]     [,3]
[1,]  0.00000 15.68387 13.98549
[2,] 15.68387  0.00000 13.49691
[3,] 13.98549 13.49691  0.00000

dist(d[1:3,], diag = T, upper = T)

         1        2        3
1  0.00000 15.68387 13.98549
2 15.68387  0.00000 13.49691
3 13.98549 13.49691  0.00000

Now that we can see thse both produce the same thing, let’s time the dist function, the C++ version and the loop.

Unit: milliseconds
  expr       min        lq      mean    median        uq       max neval
  dist  1.523120  1.581476  1.604833  1.599950  1.633967  1.692030   100
 distC  1.196444  1.234420  1.245064  1.243363  1.261691  1.302305   100
 dloop 34.157286 34.931016 40.931563 38.261996 39.007134 68.798428   100

Note that in bigger matrices the R function would eventually win out, but the point here is that we’re at least in the ballpark. Furthermore, like most languages, it is easy enough to run another language script from within that program. So if you’re working within another language for the heavy lifting, you could call R from it, or vice versa, i.e. write what you need in another language and call it from R.

Vectorization

Vectorization is a key concept in R, though some coming to R from other languages (e.g. Matlab) may be familiar with it. The basic idea is that we pass whole vectors to operations rather than the elements of it, and as such R can treat the vector once and go from there, rather than treat each element of the vector as a unique snowflake that must be considered individually, and figuring things out every time the operation is run.

As an example let’s do a simple loop vs. a vectorized approach.

x = 1:3
y = 7:9
for (i in seq_along(x)) {
  print(x[i] + y[i])
}

[1] 8
[1] 10
[1] 12

Every iteration, R evaluates what type of object x[i] is, what y[i] is, whether adding them is possible etc. If we use a vectorized approach, R figures this out once and proceeds accordingly.

x + y

[1]  8 10 12

Vectorizing a sequential operation

Let’s return to our random walk example for some benchmarks. This is a nice example because we can see that sometimes even seemingly sequential operations might be put in a form that would not require a loop. The first vectorized example uses no explicit loop but gives us the same result. It uses ifelse, which is a vectorized version of the if... else construct. However, ifelse comes with some baggage, and dplyr provides a faster version, if_else.

# vectorize
rw2d2_vec_ifelse <- function(n) {
  steps = sample(c(-1, 1), n-1, replace=T)
  xdir = sample(c(T, F), n-1, replace=T)
  xpos = c(0, cumsum(ifelse(xdir, steps, 0)))
  ypos = c(0, cumsum(ifelse(xdir, 0, steps)))
  data.frame(x=xpos, y=ypos)
}
# rw2d2_vec_ifelse(1e7)
# ggplot2::qplot(x, y, data=rw2d2_vec_ifelse(1000), geom='path')
# dplyr if_else can be faster if there are many operations
rw2d2_vec_ifelse2 <- function(n) {
  steps = sample(c(-1, 1), n-1, replace=T)
  xdir = sample(c(T, F), n-1, replace=T)
  xpos = c(0, cumsum(if_else(xdir, steps, 0)))
  ypos = c(0, cumsum(if_else(xdir, 0, steps)))
  data.frame(x=xpos, y=ypos)
}
# rw2d2_vec_ifelse2(1e7)
# ggplot2::qplot(x, y, data=rw2d2_vec_ifelse2(1000), geom='path')

However, we can still get more efficient. Here we do all the necessary sampling for the x and y positions first, then subsequently fill out those positions. We also just use the fact that cumsum can work on the data.frame itself. As data.frames are lists, cumsum will work on its elements, which in this case are the columns.

rw2d2_vec_noifelse <- function(n) {
  dir = sample(1:4, n-1, replace=T)
  xpos = c(0, c(-1, 1, 0, 0)[dir])
  ypos = c(0, c(0, 0, -1, 1)[dir])
  cumsum(data.frame(x=xpos, y=ypos))
}
# rw2d2_vec_noifelse(1e7)
# ggplot2::qplot(x, y, data=rw2d2_vec_noifelse(1000), geom='path')

Now we compare all the looping choices considered thus far:

• loop no pre-allocation
• loop with pre-allocation
• vectorization (2 versions)
• more efficient vectorization

The last version is succinct in code and much faster as well.

Unit: microseconds
                      expr        min          lq        mean      median         uq        max neval
    rw2d_loop_nopre(10000) 213722.327 229259.2290 269065.1205 264013.2695 301059.176 424019.833   100
      rw2d_loop_pre(10000)  69603.975  79123.0310  82274.5622  81806.5240  84786.049 120327.309   100
   rw2d2_vec_ifelse(10000)   1456.846   1506.2575   1597.8680   1525.6125   1545.700   8314.695   100
  rw2d2_vec_ifelse2(10000)   1217.557   1255.3860   1274.2183   1270.6350   1285.884   1432.506   100
 rw2d2_vec_noifelse(10000)    551.010    584.4395    719.0892    597.9285    615.670   6764.598   100

For more examples and issues with vectorization see the following:

A very accessible intro to vectorization Noam Ross

When to leave loops alone Hadley Wickham

R Inferno (3rd Circle) Patrick Burns

Matrix operations

Consider a linear combination to produce model predictions as in a standard regression model. We multiply each column by a coefficient and sum them.

b = runif(11) # first value is the intercept
X = matrix(rnorm(100000), ncol=10)
lincom1 = b[1] + b[2]*X[,1] + b[3]*X[,2] + b[4]*X[,3] + b[5]*X[,4] + b[6]*X[,5] + 
  b[7]*X[,6] + b[8]*X[,7] + b[9]*X[,8] + b[10]*X[,9] + b[11]*X[,10]

While this is vectorized, we can do better. We can use matrix multiplication to produce the same result. In addition, the biggest gain here is programming time. Matrix multiplication takes far less to code up.

X = cbind(1, X) # add intercept column
lincom2  = X %*% b

However, even within matrix operations there might be faster ways to do it. In this case, crossprod is slightly faster.

lincom3 = tcrossprod(X, t(b))
identical(lincom1, c(lincom2))

[1] TRUE

identical(lincom2, lincom3)

[1] TRUE

Let’s time them.

Unit: microseconds
 expr     min      lq      mean   median       uq      max neval
  lc1 590.598 613.324 732.61942 627.9860 635.4640 4235.645   100
  lc2 138.999 140.465 147.76094 146.7695 154.5405  197.061   100
  lc3  56.303  57.770  61.37959  59.8220  62.9015   80.937   100

They all produce identical output, but using a matrix operation is over 4 times as fast as the simple sum, and crossprod over 11 times as fast. Programming time is quicker also.

Even how we do the matrix operations can have an effect. The following examples use various means to obtain coefficients from a standard linear regression model. We’ll start with just the base R lm function. The next is the lm.fit function, which is what lm uses, but after a lot of pre- and post-processing of the data. After that we use basic matrix multiplication, a grouped version of it, followed by use of the crossprod function.

set.seed(1234)
n = 10000000                   # sample size
x = rnorm(n)                   # covariate
y =  5 + 2*x + rnorm(n)        # target
X = cbind(1, x)                # model.matrix
# lm
lm_base  = lm(y ~ x)
# use the underlying fit function
lm_fit  = lm.fit(X, y)          # note: qr.solve isn't really a gain here but could also be used
# an even faster version
lm_fit2  = .lm.fit(X, y)
# normal equations to get coefficients
matrix_raw = solve(t(X) %*% X) %*% t(X) %*% y
matrix_grouped = solve(t(X)%*%X) %*% (t(X) %*% y)
# crossprod approach
matrix_crossprod = crossprod(solve(crossprod(X)), crossprod(X,y))
cbind(cbind(lm_base=coef(lm_base), 
      lm_fit=coef(lm_fit), 
      lm_fit2=coef(lm_fit2), 
      matrix_raw=matrix_raw[,1],
      matrix_grouped=matrix_grouped[,1],
      matrix_crossprod=matrix_crossprod[,1]))

             lm_base   lm_fit  lm_fit2 matrix_raw matrix_grouped matrix_crossprod
(Intercept) 4.999339 4.999339 4.999339   4.999339       4.999339         4.999339
x           1.999691 1.999691 1.999691   1.999691       1.999691         1.999691

Now for the timings.

Unit: milliseconds
             expr        min         lq       mean     median        uq        max neval
          lm_base 8518.96348 8553.59319 8785.89097 8652.52351 8846.2702 9598.13285    20
           lm_fit  608.17270  610.11326  725.07545  611.13302  846.8290 1431.15979    20
          lm_fit2  492.01975  495.77271  576.94008  496.67547  747.8333  810.55111    20
       matrix_raw  387.49752  390.13028  472.59917  392.63900  649.1642  676.84504    20
   matrix_grouped  246.16266  247.35969  380.93035  248.66786  371.0785 1225.53339    20
 matrix_crossprod   51.68899   52.91857   53.74915   53.54319   54.0288   57.83074    20

We can see above that these approaches all produce the same results, but take very different amounts of time to complete. We get a significant gain if we reduce the lm overhead, where lm.fit is almost 10x faster, and the .lm.fit, which is just a wrapper for the C code call, is faster still. We can still get faster by using matrix operations, and even get a slight gain if we group them a certain way. This reduces memory usage as well. The idea is that if we go from solve(t(X)%*%X) and post multiply by t(X), we’re dealing with a 2 x 10⁷ matrix. If we first t(X) %*% y, we’re now just multiplying a 2x2 matrix with a 2x1 matrix. In the end though, crossprod is the most optimized, and notably faster than any of the other approaches.

And finally, there are efficient methods, particularly for sparse matrices, in packages such as Matrix. For dense matrices it has the same speed as standard R, but would still have other within-package functionality that might be useful.

library(Matrix)
X2 = as(X, 'dgeMatrix')
y2 = as(y, 'dgeMatrix')

Unit: milliseconds
              expr      min       lq     mean   median       uq      max neval
  matrix_crossprod 49.23980 53.24276 53.59032 53.66591 54.01502 59.17703   100
 matrix_crossprod2 52.04411 53.28352 53.68203 53.81811 54.10270 55.42700   100

Here’s a comparison with sparse matrices, where memory as well as speed gains can be made. The base R object is 200 times larger.

m <- matrix(0, 500, 500)
m[runif(314, 0, length(m))] <- 1
mm <- as(m, "dgCMatrix")
object.size(m) / object.size(mm) # 

278.114571746385 bytes

The speed gains are notable as well.

microbenchmark(sparse_tcrossprod <- tcrossprod(mm),
               sparse_crossprod <- crossprod(mm),
               base_crossprod <- crossprod(m),
               base_tcrossprod <- tcrossprod(m))

Unit: microseconds
                                expr       min         lq       mean     median        uq       max neval
 sparse_tcrossprod <- tcrossprod(mm)   137.532   157.6205   183.0385   172.4285   198.381   332.248   100
   sparse_crossprod <- crossprod(mm)   141.932   157.7670   183.0298   168.7635   194.275   339.286   100
      base_crossprod <- crossprod(m) 49594.339 50380.0910 50557.3464 50560.1435 50721.575 51607.179   100
    base_tcrossprod <- tcrossprod(m)   929.296   974.6030  1016.5842  1009.3525  1052.460  1196.737   100

Other R builds

Other versions of R exist. For those using university high performance computing, the R build you are using will be different, for example. Here is a run down of some, partly taken from Wickham’s chapter linked in the introduction of this document. There are others, but I just note the ones that still appear to be in regular development.

Microsoft Open R: From Microsoft, specifically the former Revolution Analytics group. While the build is typically a point release behind the current version of R, that won’t affect most users, and it’s easy to install on one’s own machine just like any R.

pqR (pretty quick R): by Radford Neal. Fixes many obvious performance issues, and provides better memory management and some support for automatic multithreading.

Renjin: by BeDataDriven. Renjin uses the Java virtual machine, and has an extensive test suite.

Summary

To complete the partial quote referred to earlier…

“Programmers waste enormous amounts of time thinking about, or worrying about, the speed of noncritical parts of their programs, and these attempts at efficiency actually have a strong negative impact when debugging and maintenance are considered. We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil. Yet we should not pass up our opportunities in that critical 3%.”

~ Donald Knuth

There are no good benchmarks for the balance of programming vs. computational time for statistical analysis and related activities. However, the reason R will beat other languages hands down in many cases is because it is so much easier to program certain, often very common operations. Furthermore, the functionality available across a wide variety of modeling paradigms is unmatched at present. Knowing just a few ‘tricks’ such as those demonstrated here will help gain those small efficiencies in the 97% time. In the end, what approach works best will depend on the specifics of the task, but you will always have options with R, and often pretty fast ones.

References

Two excellent texts on R programming are:

The Art of R Programming by Norm Matloff

Advanced R by Hadley Wickham.