R Optimisation Patterns

Subcategory: None · Core

Don't Grow data.frame(s)

When you grow a data.frame by one, for example using rbind() or cbind() to add a row or column, R will allocate a new data.frame one larger and copy across all the old data. If you’re doing this regularly, that’s a lot of redundant expensive memory movement. Either pre-allocate your data.frame to the correct size at the start, or use a data.table instead which provides faster growth (and many other performance improvements).

Read More

When you grow a data.frame by one, for example using rbind() or cbind() to add a row or column, R will allocate a new data.frame one larger and copy across all the old data. If you’re doing this regularly, that’s a lot of redundant expensive memory movement. Either pre-allocate your data.frame to the correct size at the start, or use a data.table instead which provides faster growth (and many other performance improvements).

Prior to R 3.4, vectors suffered from the same overhead during growth. However, the newer implementation ensures this overhead of redundant copies is not incurred most of the time. It may still occur when there are multiple references to a growing vector, so beware.

Example Code

It might feel natural to grow a data.frame row-by-row as you collect data.

# Allocate a data.frame with 3 columns
df <- data.frame(x = numeric(), y = numeric(), z = numeric())
# Append 10,000 rows
for (i in 1:10000) {
  new_row <- data.frame(x = i, y = i * 2, z = i * 3)
  df <- rbind(df, new_row)
}

This took 4.42 seconds (141.91 seconds for 100,000 rows!).

But if we initialise the data-frame with the correct number of rows, and instead update the elements in-place.

# Allocate a data.frame with 3 columns and preallocate 10,000 rows.
df <- data.frame(x = numeric(10000), y = numeric(10000), z = numeric(10000))
# Append 10,000 rows
for (i in 1:10000) {
  df$x[i] <- i
  df$y[i] <- i * 2
  df$z[i] <- i * 3
}

This now executes in 0.81 seconds (26.56 seconds for 100,000 rows).

Using lapply followed by do.call to rbind the rows is slightly faster:

list_of_rows <- lapply(1:100000, function(i){
    data.frame(x = i, y = i * 2, z = i * 3)
  })
df <- do.call(rbind, list_of_rows)

This executes in ~0.75 seconds (20.78 seconds for 100,000 rows).

We can however improve this even further, by taking advantage of vectorisation.

This can be done using the data.table library:

# Import the data.table library
library(data.table)
df <- data.table(x = numeric(10000), y = numeric(10000), z = numeric(10000))
rows <- 1:10000
df[rows, `:=`(
  x = rows,
  y = rows * 2,
  z = rows * 3
)]

This now executes in 0.05 seconds (and doesn’t get slower for 100,000 elements).

If you need to convert it back to a data.frame, you can call as.data.frame() passing your data.table as a parameter.

The dplyr library is even faster in this case:

# Import the dplyr package
library(dplyr)
df <- data.frame(x = 1:10000)
dplyr::mutate(df, y = x * 2, z = x * 3)

This executes in 0.001 seconds (and doesn’t get slower for 100,000 elements).

The Technical Detail

R data-frames are column major, and many operations force the columns to be re-allocated. In particular, if you add a row each column of the data-table will be re-allocated.

data.table is almost always equal if not faster in performance to the core data.frame, so there’s little reason not to use it.

Loop ordering

In Fortran, MATLAB, and R, arrays are stored in column-major order, meaning that the entries from the left-most index are stored contiguously in memory. This is often referred to as the left-most index varying quickest. This is different from many other popular scientific programming languages such as C, C++, and Python, which use a row-major ordering. The implication for writing loops is that it is always best to order loops by which indices vary quickest. In this way, the loop nest will traverse contiguous memory, which can be done efficiently.

Read More

In Fortran, MATLAB, and R, arrays are stored in column-major order, meaning that the entries from the left-most index are stored contiguously in memory. This is often referred to as the left-most index varying quickest. This is different from many other popular scientific programming languages such as C, C++, and Python, which use a row-major ordering. The implication for writing loops is that it is always best to order loops by which indices vary quickest. In this way, the loop nest will traverse contiguous memory, which can be done efficiently.

Example Code

Consider the following Fortran example involving a loop over a 3D array, which uses the recommended column-major ordering:

integer, parameter :: n = 500
real, dimension(n, n, n) :: a, b, c
integer :: i, j, k

! Initialisation of a and b omitted

! Hand-coded c = a + b using column-major
do k = 1, n
  do j = 1, n
    do i = 1, n
      c(i,j,k) = a(i,j,k) + b(i,j,k)
    end do
  end do
end do

Without compiler optimisations, this will likely run several times faster than the (not recommended) equivalent row-major version:

integer, parameter :: n = 500
real, dimension(n, n, n) :: a, b, c
integer :: i, j, k

! Initialisation of a and b omitted

! Hand-coded c = a + b using row-major
do i = 1, n
  do j = 1, n
    do k = 1, n
      c(i,j,k) = a(i,j,k) + b(i,j,k)
    end do
  end do
end do

With compiler optimisations, simple loops such as this where the result of the sum are unused would however likely be optimised out by the compiler. In practice you wouldn’t normally have a redundant loop in your code, so this shouldn’t be a concern.

A further comment on the simple loops above is that the same could be achieved with an array operation, which has a more compact notation:

integer, parameter :: n = 500
real, dimension(n, n, n) :: a, b, c

! Initialisation of a and b omitted

! Compute c = a + b with an array operation
c(:,:,:) = a + b

Benchmark

Below these examples have been extended into a full benchmark, to demonstrate the impact. An additional checksum is used to ensure the code we’re benchmarking isn’t optimised away by the compiler.

Benchmark Code

```fortran program benchmark_loops implicit none integer, parameter :: n = 1000 real, allocatable :: a(:,:,:), b(:,:,:), c(:,:,:) real :: t1, t2 integer :: i real :: checksum allocate(a(n,n,n), b(n,n,n), c(n,n,n)) ! Initialize data a = 1.0 b = 2.0 c = 0.0 ! Warm-up runs (optional but recommended) call col_major_add(a, b, c, n) call row_major_add(a, b, c, n) ! Benchmark column-major loop order call cpu_time(t1) call col_major_add(a, b, c, n) call cpu_time(t2) print *, "Column-major time (s): ", t2 - t1 checksum = sum(c) print *," checksum: ", checksum ! Reset result array to avoid reuse artifacts c = 0.0 ! Benchmark row-major loop order call cpu_time(t1) call row_major_add(a, b, c, n) call cpu_time(t2) print *, "Row-major time (s): ", t2 - t1 checksum = sum(c) print *," checksum: ", checksum contains subroutine col_major_add(a, b, c, n) integer, intent(in) :: n real, intent(in) :: a(n,n,n), b(n,n,n) real, intent(out) :: c(n,n,n) integer :: i, j, k do k = 1, n do j = 1, n do i = 1, n c(i,j,k) = a(i,j,k) + b(i,j,k) end do end do end do end subroutine col_major_add subroutine row_major_add(a, b, c, n) integer, intent(in) :: n real, intent(in) :: a(n,n,n), b(n,n,n) real, intent(out) :: c(n,n,n) integer :: i, j, k do i = 1, n do j = 1, n do k = 1, n c(i,j,k) = a(i,j,k) + b(i,j,k) end do end do end do end subroutine row_major_add end program benchmark_loops ```

Compiled with GFortran and compiler optimisations enabled, column-major accesses were 2.4x faster under Ubuntu!

$ gfortran -O3 benchmark.f90
$ ./a.out
 Column-major time (s):    3.98375320    
  checksum:    67108864.0    
 Row-major time (s):       9.72134018    
  checksum:    67108864.0

Under WSL on different hardware a greater 7.1x speedup was seen, so results are likely to vary, but should remain positive.

$ gfortran -O3 benchmark.f90
$ ./a.out
 Column-major time (s):   0.440809250
  checksum:    67108864.0
 Row-major time (s):       3.11280060
  checksum:    67108864.0

Technical Details

The underlying reason that order has such a big impact here is due to how computer memory operates. When a processor loads a variable from RAM into it’s caches, it doesn’t load an individual variable (likely 4 bytes), it loads a full cache line (likely 64 bytes).

Therefore, variables stored consecutively in memory will be loaded at the same time if they fall within the same cache line. This means the processor does not need to go to RAM to access the variable, which is orders of magnitude higher latency. Due to how Fortran lays out memory, column-major ordering achieves this.

In contrast, if you operate as though memory is row-major, consecutively accessed variables will be stored a long way apart from each other in memory. Hence each individual access will need to load a full cache line from RAM. Due to this high turnover, it’s likely by the time a second variable would be accessed from any cache line, that the cache line has been evicted to be replaced with a fresher load.

If you’re working with 4-byte types you could be doing 16x the RAM accesses, with 8-byte types it’s still 8x!

Furthermore, in order for the compiler to perform vectorisation, it needs to be able to apply vector instructions to a full cache line of variables. If you are not operating consecutively on variables within a cache line, it can be much harder for the compiler to detect that vectorisation is appropriate.

Vectorisation

R has many functions which operate on entire vectors and matrices, which have been implemented in faster C. Vectorisation can perform mathematical operations many times faster than a for loop and even take advantage of conditional logic.

Read More

R has many functions which operate on entire vectors and matrices, which have been implemented in faster C. Vectorisation can perform mathematical operations many times faster than a for loop and even take advantage of conditional logic.

Note, when R users talk about vectorisation they’re typically referring to functions which operate on whole vectors, not advanced SIMD instructions.

Example Code

To take advantage of vectorisation within R, operations should be performed at an array or matrix level, rather than by looping over the individual elements or rows.

The below example sums the rows within a 1000x1000 matrix to produce an output array with length 1000.

# Allocate an 1000x1000 matrix
mat <- matrix(1:1e6, nrow = 1000)
# Allocate an output array length 1000
x <- integer(length = nrow(mat))
# Calculate the sum of each row into the output array
for (i in 1:nrow(mat)) { 
  x[i] <- sum(mat[i,]) 
}

To use vectorisation, it can be rewritten as a single line of code using rowSums().

# Allocate an 1000x1000 matrix
mat <- matrix(1:1e6, nrow = 1000)
tic()
# Call the rowSums() method
x <- rowSums(mat)
toc()

In this example vectorisation also reduces the runtime from 0.05 to 0.02 seconds.

Most R functions support vectorisation, so it’s important to try and use purpose built functions to process your arrays and matrices, rather than manually looping.

Not all of your R code will be this simple though, maybe you have conditional logic which restricts which elements the operation should be performed on.

# Allocate an array of 1 million elements
x <- 1:1e6
# Divide even elements by 2, multiply odd elements by 3

for (i in x) {
  if (i %% 2 == 0) {
    i <- i / 2         # even: divide by 2
  } else {
    i <- i * 3         # odd: multiply by 3
  }
}

The above example can be written using ifelse() as

# Allocate an array of 1 million elements
x <- 1:1e6
# Divide even elements by 2, multiply odd elements by 3
result <- ifelse(x %% 2 == 0,         # condition: even numbers
                 x / 2,               # if TRUE: divide even numbers by 2
                 x * 3)               # if FALSE: multiply odd numbers by 3

The ifelse() approach reduced the runtime from 0.35 seconds to 0.06 seconds.

Using vectorisation is both more concise and more performant.

You might be tempted to use sapply() or vapply(), which allow you to apply a transformation function to each element within an array, but these were slower than the original loop in testing.

This is because the apply family of functions, in contrast to popular believe, are not more performant than for loops in R. In fact, the apply implementation contains a regular R for loop. The reason they may sometimes be more performant than for loops is because they take care of the common pitfalls with loops, like the growing dataframes problem. However, using an apply function is not vectorisation, it is loop-hiding (see also the Fourth Circle of The R Inferno).

The Technical Detail

Many of R’s mathematical vector functions rely on the scientific libraries BLAS and LAPACK, which provide efficient implementations of the operations. These are faster when executed in a single call to C, rather than an R loop over a vector which makes many individual calls.

BLAS and LAPACK can take advantage of even faster SIMD (Single Instruction, Multiple Data) instructions to operate on several values in parallel which are packed into a register. However, most R installations are provided with a reference implementation of these libraries which is unlikely to take advantage of these. It’s currently non-trivial to install R with an accelerated version of BLAS, but it’s likely your will find one setup on HPC systems.