Measuring code speed

Comparing code speed

The base R language is very, very slow. However, a lot of the packages and functions in R have been written in C/C++/Fortran and are very, very fast. So depending on what functions we use to do things, we can get significantly faster performance. Of course, this depends on us being able to measure how fast our code is in the first place.

To measure how fast a bit of code is, we can benchmark it. The easiest way of benchmarking our code is with the microbenchmark package. microbenchmark provides just one function - microbenchmark().

Let’s try an example. There are two ways of doing a square root - which is faster?

  • sqrt(number)
  • number ^ 0.5

microbenchmark() will run each set of code 100 times (can be changed with the times argument), then return a summary table of how long it took our code to run. We are usually interested in the median value (the mean is more affected by crazy outliers in execution speed).

library(microbenchmark)
microbenchmark(
               sqrt(10),
               10 ^ 0.5
               )
## Unit: nanoseconds
##      expr min    lq   mean median    uq  max neval
##  sqrt(10) 109 117.0 159.03  128.0 149.0 2377   100
##    10^0.5 232 245.5 325.82  251.5 270.5 6359   100

In this case, the sqrt() function was quite a bit faster, likely because of some optimization specific to doing a square root (as opposed to exponentiation operations in general). Now we can measure how fast two pieces of code are.

The R compiler

R has a JIT (just in time) compiler that is enabled by default. You don’t need to worry about compiling your code with the compiler::cmpfun() function, because it is already done for you automatically!

Profiling code

microbenchmark() is an extremely useful tool for determining which several blocks of several approaches is faster. But it does not highlight areas that need to be improved - which elements of a chunk of code take the longest amount of time to run. This second technique is called profiling. It is used to identify the slow parts of our code so we can focus on fixing them.

We’ll use the profvis package for this task. profvis executes a chunk of code (everything between the {}), stopping it at regular interals to determine which functions are currently executing, then displays a nice interactive graph that shows where time was spent in each line of our code. It also steps into the functions that get called by our code - there is no need to profvis() our function, then profvis() the code that it calls, etc., etc.

library(profvis)
library(nycflights13)
library(tidyverse)
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Conflicts with tidy packages ----------------------------------------------
## filter(): dplyr, stats
## lag():    dplyr, stats
library(stringr)

profvis({
  # quick and dirty aggregation by plane and manufacturer
  # not using pipes to break up the code a bit
  planes_simple <- planes %>% 
    mutate(model=str_replace(model, "-.*$", "")) %>% 
    mutate(manufacturer=map_chr(str_split(manufacturer, " "), ~.x[[1]])) %>% 
    group_by(model, manufacturer) %>% 
    summarize(number = n()) %>% 
    filter(number > 3)
  
  ggplot(planes_simple, aes(x=model, y=number, fill=manufacturer)) +
    geom_bar(stat = "identity") +
    theme_classic() + theme(axis.text.x = element_text(angle=45, hjust=1))
}, interval = 0.001)
## Intervals smaller than ~5ms will probably not result in accurate timings.

In this case, everything run through the pipe appears as one giant line of code. This can make it very tricky to identify what’s going on and when. In cases like this, it’s a good idea to break apart our pipe and have each as a separate expression. This lets us easily see how much time was spent in each line of code individual dplyr function.

profvis({
  planes_simple <- planes
  planes_simple <- mutate(planes_simple, model=str_replace(model, "-.*$", ""))
  planes_simple <- mutate(planes_simple, manufacturer=map_chr(str_split(manufacturer, " "), ~.x[[1]]))
  planes_simple <- group_by(planes_simple, model, manufacturer)
  planes_simple <- summarize(planes_simple, number = n())
  planes_simple <- filter(planes_simple, number > 3)
}, interval = 0.001)
## Intervals smaller than ~5ms will probably not result in accurate timings.