Chapter 5 Performance

Some resources used here or for further reading:

The people who say that “R is just always slow” are usually not great R programmers. It is true that writing inefficient R code is easy, yet writing efficient R code is also possible when you know what you’re doing. In this chapter, you will learn how to write R(cpp) code that is fast.

5.1 R’s memory management

Read more with this chapter of Advanced R.

5.1.1 Understanding binding basics

x <- c(1, 2, 3)
  • It’s creating an object, a vector of values, c(1, 2, 3).
  • And it’s binding that object to a name, x.


y <- x

There are now two names for the same object in memory.

5.1.2 Copy-on-modify

x <- c(1, 2, 3)
y <- x
y[3] <- 4
x
#> [1] 1 2 3

The object in memory is copied before being modified, so that x is not modified.

5.1.3 Copy-on-modify: what about inside functions?

f <- function(a) {
  a
}
x <- c(1, 2, 3)
z <- f(x)

f2 <- function(a) {
  a[1] <- 10
  a
}
z2 <- f2(x)
cbind(x, z2)
#>      x z2
#> [1,] 1 10
#> [2,] 2  2
#> [3,] 3  3

The input parameter is not modified; you operate on a local copy of a = x in f2().

5.1.4 Lists

It’s not just names (i.e. variables) that point to values; elements of lists do too.

l1 <- list(1, 2, 3)

l2 <- l1

5.1.5 Copy-on-modify for lists?

l2[[3]] <- 4

Only the third element needs to be copied.

5.1.6 Data frames

Data frames are lists of vectors.

d1 <- data.frame(x = c(1, 5, 6), y = c(2, 4, 3))

d2 <- d1
d2[, 2] <- d2[, 2] * 2  # modify one column

d3 <- d1
d3[1, ] <- d3[1, ] * 3  # modify one row

By modifying the first row, you’re modifying the first element of all vectors, therefore the full data frame is copied..

5.2 Early advice

5.2.1 NEVER GROW A VECTOR

Example computing the cumulative sums of a vector:

x <- rnorm(2e4)  # Try also with n = 1e5
system.time({
  current_sum <- 0
  res <- c()
  for (x_i in x) {
    current_sum <- current_sum + x_i
    res <- c(res, current_sum)
  }
})
#>    user  system elapsed 
#>    0.46    0.43    0.89

Here, at each iteration, you are reallocating a vector (of increasing size). Not only computations take time, memory allocations do too. This makes your code quadratic with the size of x (if you multiply the size of x by 2, you can expect the execution time to be multiplied by 4, for large sample sizes), whereas it should be only linear.

What happens is similar to if you would like to climb these stairs, you climb one stair, go to the bottom, then climb two stairs, go to bottom, climb three, and so on. That takes way more time than just climbing all stairs at once.

A good solution is to always pre-allocate your results (if you know the size):

system.time({
  current_sum <- 0
  res2 <- double(length(x))
  for (i in seq_along(x)) {
    current_sum <- current_sum + x[i]
    res2[i] <- current_sum
  }
})
#>    user  system elapsed 
#>       0       0       0
all.equal(res2, res)
#> [1] TRUE

If you don’t know the size of the results, you can store them in a list and merge them afterwards:

system.time({
  current_sum <- 0
  res3 <- list()
  for (i in seq_along(x)) {
    current_sum <- current_sum + x[i]
    res3[[i]] <- current_sum
  }
})
#>    user  system elapsed 
#>    0.01    0.00    0.01
all.equal(unlist(res3), res)
#> [1] TRUE

With recent versions of R (>= 3.4), you can efficiently grow a vector using

system.time({
  current_sum <- 0
  res4 <- c()
  for (i in seq_along(x)) {
    current_sum <- current_sum + x[i]
    res4[i] <- current_sum
  }
})
#>    user  system elapsed 
#>    0.02    0.00    0.01
all.equal(res4, res)
#> [1] TRUE

Assigning to an element of a vector beyond the current length now over-allocates by a small fraction. The new vector is marked internally as growable, and the true length of the new vector is stored in the truelength field. This makes building up a vector result by assigning to the next element beyond the current length more efficient, though pre-allocating is still preferred. The implementation is subject to change and not intended to be used in packages at this time. (NEWS)

An even better solution would be to avoid the loop by using a vectorized function:

system.time(res5 <- cumsum(x))
#>    user  system elapsed 
#>       0       0       0
all.equal(res5, res)
#> [1] TRUE
x <- rnorm(1e7)
system.time(cumsum(x))
#>    user  system elapsed 
#>    0.08    0.00    0.08

As a second example, let us generate a matrix of uniform values (max changing for every column):

n <- 1e3
max <- 1:1000
system.time({
  mat <- NULL
  for (m in max) {
    mat <- cbind(mat, runif(n, max = m))
  }
})
#>    user  system elapsed 
#>    0.68    0.53    1.25
apply(mat, 2, max)[1:10]
#>  [1] 0.999412 1.998023 2.996562 3.997387 4.992261 5.993579 6.990147 7.997312 8.998179
#> [10] 9.994124

Instead, we should pre-allocate a matrix of the right size:

system.time({
  mat3 <- matrix(0, n, length(max))
  for (i in seq_along(max)) {
    mat3[, i] <- runif(n, max = max[i])
  }
})
#>    user  system elapsed 
#>    0.01    0.01    0.04
apply(mat3, 2, max)[1:10]
#>  [1] 0.9999812 1.9995295 2.9990553 3.9976881 4.9972887 5.9925429 6.9920594 7.9991727
#>  [9] 8.9987023 9.9951227

Or we could use a list instead. What is nice with using a list is that you don’t need to pre-allocate. Indeed, as opposed to atomic vectors, each element of a list is in different places in memory so that you don’t have to reallocate all the data when you add an element to a list.

system.time({
  l <- list()
  for (i in seq_along(max)) {
    l[[i]] <- runif(n, max = max[i])
  }
  mat4 <- do.call("cbind", l)
})
#>    user  system elapsed 
#>    0.03    0.00    0.03
apply(mat4, 2, max)[1:10]
#>  [1] 0.9989994 1.9984834 2.9996546 3.9997688 4.9988993 5.9976521 6.9924608 7.9990048
#>  [9] 8.9998035 9.9953478

Instead of pre-allocating yourself, you can use sapply (or lapply and calling do.call() after, as previously done):

system.time(
  mat4 <- sapply(max, function(m) runif(n, max = m))
)
#>    user  system elapsed 
#>    0.03    0.00    0.03
apply(mat4, 2, max)[1:10]
#>  [1] 0.9998305 1.9984741 2.9992208 3.9964728 4.9992998 5.9916205 6.9962307 7.9888878
#>  [9] 8.9998210 9.9917989

Don’t listen to people telling you that sapply() is a vectorized operation that is so much faster than loops. That’s false, and for-loops can actually be much faster than sapply() when using just-in-time (JIT) compilation. You can learn more with this blog post.

5.2.2 Use the right function

Often, in order to optimize your code, you can simply find the right function to do what you need to do.

For example, rowMeans(x) is much faster than apply(x, 1, mean). Similarly, if you want more efficient functions that apply to rows and columns of matrices, you can check package {matrixStats}.

Another example is when reading large text files; in such cases, prefer using data.table::fread() rather than read.table().

Generally, packages that uses C/Rcpp are efficient.

5.2.3 Do not try to optimize everything

“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.”

— Donald Knuth.

If you try to optimize each and every part of your code, you will end up losing a lot of time writing it, and it will probably make your code less readable.

R is great at prototyping quickly because you can write code in a concise and easy way. Start by doing just that. If performance matters, then profile your code to see which part of your code is taking too much time and optimize only this part!

Learn more on how to profile your code in RStudio in this article.

5.3 Vectorization

I call vectorized a function that takes vectors as arguments and operate on each element of these vectors in another (compiled) language (such as C, C++ and Fortran). Again, sapply() is not a vectorized function (cf. above).

Take this code:

N <- 10e3; x <- runif(N); y <- rnorm(N)

res <- double(length(x))
for (i in seq_along(x)) {
  res[i] <- x[i] + y[i]
}

As an interpreted language, for each iteration res[i] <- x[i] + y[i], R has to ask:

  • what is the type of x[i] and y[i]?

  • can I add these two types?

  • what is the type of x[i] + y[i] then?

  • can I store this result in res or do I need to convert it?

These questions must be answered for each iteration, which takes time.
Some of this is alleviated by JIT compilation.

On the contrary, for vectorized functions, these questions must be answered only once, which saves a lot of time. Read more with Noam Ross’s blog post on vectorization.

5.3.1 Exercise

Monte-Carlo integration (example from book Efficient R programming)

Suppose we wish to estimate the integral \(\int_0^1 x^2 dx\) using a Monte-Carlo method. Essentially, we throw darts (at random) and count the proportion of darts that fall below the curve (as in the following figure).

Naively implementing this Monte-Carlo algorithm in R would typically lead to something like:

monte_carlo <- function(N) {
  
  hits <- 0
  
  for (i in seq_len(N)) {
    x <- runif(1)
    y <- runif(1)
    if (y < x^2) {
      hits <- hits + 1
    }
  }
  
  hits / N
}

This takes a few seconds for N = 1e6:

N <- 2e5
system.time(res <- monte_carlo(N))
#>    user  system elapsed 
#>    0.52    0.30    0.83
res
#> [1] 0.33386

Your task: find a vectorized solution for this problem:

system.time(res2 <- monte_carlo_vec(N))
#>    user  system elapsed 
#>    0.01    0.00    0.02
res2
#> [1] 0.33295

5.4 Rcpp

See this presentation.

You have this data and this working code (a loop) that is slow

mydf <- readRDS(system.file("extdata/one-million.rds", package = "advr38pkg"))

QRA_3Dmatrix <- array(0, dim = c(max(mydf$ID), max(mydf$Volume), 2))
transform_energy <- function(x) {
  1 - 1.358 / (1 + exp( (1000 * x - 129000) / 120300 ))
}

for (i in seq_len(nrow(mydf))) {
  # Row corresponds to the ID class
  row    <- mydf$ID[i]
  # Column corresponds to the volume class
  column <- mydf$Volume[i]
  # Number of events, initially zero, then increment
  QRA_3Dmatrix[row, column, 1] <- QRA_3Dmatrix[row, column, 1] + 1  
  # Sum energy 
  QRA_3Dmatrix[row, column, 2] <- QRA_3Dmatrix[row, column, 2] + 
    transform_energy(mydf$Energy[i])
}

Rewrite this for-loop with Rcpp.

You can also try to use {dplyr} for this problem.

5.5 Linear algebra

It is faster to use crossprod(X) and tcrossprod(X) instead of t(X) %*% X and X %*% t(X). Moreover, using A %*% (B %*% y) is faster than A %*% B %*% y, and solve(A, y) is faster than solve(A) %*% y.

Don’t re-implement linear algebra operations (such as matrix products) yourself. There exist some highly optimized libraries for this. If you want to use linear algebra in Rcpp, try RcppArmadillo or RcppEigen.

If you want to use some optimized multi-threaded linear library, you can try Microsoft R Open.

5.5.1 Exercises

Compute the Euclidean distances between each of row of X and each row of Y:

set.seed(1)
X <- matrix(rnorm(5000), ncol = 5)
Y <- matrix(rnorm(2000), ncol = 5)

A naive implementation would be:

system.time({
  dist <- matrix(NA_real_, nrow(X), nrow(Y))
  for (i in seq_len(nrow(X))) {
    for (j in seq_len(nrow(Y))) {
      dist[i, j] <- sqrt(sum((Y[j, ] - X[i, ])^2))
    }
  }
})
#>    user  system elapsed 
#>    0.45    0.05    0.50

Try first to remove one of the two loops using sweep() instead. Which loop should you choose to remove ideally? A solution with sweep() can take

#>    user  system elapsed 
#>    0.08    0.01    0.09

Then, try to implement a fully vectorized solution based on this hint: \(\text{dist}(X_i, Y_j)^2 = (X_i - Y_j)^T (X_i - Y_j) = X_i^T X_i + Y_j^T Y_j - 2 X_i^T Y_j\). A faster solution using outer() and tcrossprod() takes

#>    user  system elapsed 
#>    0.02    0.00    0.01

5.6 Algorithms & data structures

Sometimes, getting the right data structure (e.g. using a matrix instead of a data frame or integers instead of characters) can save you some computation time.

Is your algorithm doing some redundant computations making it e.g. quadratic instead of linear with respect to the dimension of your data?

See exercises (section 5.7) for some insights.

You can also find a detailed example in this blog post.

5.7 Exercises

Generate \(10^7\) (start with \(10^5\)) steps of the process described by the formula:\[X(0)=0\]\[X(t+1)=X(t)+Y(t)\] where \(Y(t)\) are independent random variables with the distribution \(N(0,1)\). Then, calculate the percentage of \(X(t)\) that are negative. You do not need to store all values of \(X\).

A naive implementation with a for-loop could be:

set.seed(1)

system.time({
  N <- 1e5
  x <- 0
  count <- 0
  for (i in seq_len(N)) {
    y <- rnorm(1)
    x <- x + y
    if (x < 0) count <- count + 1
  }
  p <- count / N
})
#>    user  system elapsed 
#>    0.19    0.05    0.23
p
#> [1] 0.88454

Try to vectorize this after having written the value of X(0), X(1), X(2), and X(3). What would be the benefit of writing an Rcpp function over a simple vectorized R function?

set.seed(1)
system.time(p2 <- advr38pkg::random_walk_neg_prop(1e5))
#>    user  system elapsed 
#>    0.01    0.00    0.02
p2
#> [1] 0.88454
set.seed(1)
system.time(p3 <- advr38pkg::random_walk_neg_prop(1e7))
#>    user  system elapsed 
#>    0.41    0.00    0.40
p3
#> [1] 0.3400444
mat <- as.matrix(mtcars)
ind <- seq_len(nrow(mat))
mat_big <- mat[rep(ind, 1000), ]  ## 1000 times bigger dataset
last_row <- mat_big[nrow(mat_big), ]

Speed up these loops (vectorize):

system.time({
  for (j in 1:ncol(mat_big)) {
    for (i in 1:nrow(mat_big)) {
      mat_big[i, j] <- 10 * mat_big[i, j] * last_row[j]
    }
  }
})
#>    user  system elapsed 
#>    0.43    0.00    0.42

Why colSums() on a whole matrix is faster than on only half of it?

m0 <- matrix(rnorm(1e6), 1e3, 1e3)
microbenchmark::microbenchmark(
  colSums(m0[, 1:500]), 
  colSums(m0)
)
#> Unit: microseconds
#>                  expr    min     lq     mean  median     uq    max neval
#>  colSums(m0[, 1:500]) 1710.6 1841.2 2320.028 2002.25 2310.6 9423.6   100
#>           colSums(m0)  756.4  853.1  959.932  912.60 1042.6 1522.0   100

Try to speed up this code by vectorizing it first, and/or by precomputing. Then, recode it in Rcpp and benchmark all the solutions you came up with.

M <- 50
step1 <- runif(M)
A <- rnorm(M)
N <- 1e4

tau <- matrix(0, N + 1, M)
tau[1, ] <- A
for (j in 1:M) {
  for (i in 2:nrow(tau)) {
    tau[i, j] <- tau[i - 1, j] + step1[j] * 1.0025^(i - 2)
  }
}

Make a fast function that counts the number of elements between a sequence of breaks. Can you do it in base R? Try also implementing it in Rcpp. How can you implement a solution whose computation time doesn’t depend on the number of breaks? [Which are the special cases that you should consider?]

x <- sample(10, size = 1e4, replace = TRUE)
breaks <- c(1, 3, 8.5, 9.5, 10)
table(cut(x, breaks), exclude = NULL) # does not include first break (1)
#> 
#>     (1,3]   (3,8.5] (8.5,9.5]  (9.5,10]      <NA> 
#>      2006      4970      1027       944      1053
hist(x, breaks, plot = FALSE)$counts  # includes first break
#> [1] 3059 4970 1027  944
advr38pkg::count_by_breaks(x, breaks)
#> [1] 2006 4970 1027  944
advr38pkg::count_by_breaks_fast(x, breaks)
#> [1] 2006 4970 1027  944
microbenchmark::microbenchmark(
  table(cut(x, breaks)), 
  hist(x, breaks, plot = FALSE)$counts, 
  advr38pkg::count_by_breaks(x, breaks, use_outer = TRUE),
  advr38pkg::count_by_breaks(x, breaks, use_outer = FALSE),
  advr38pkg::count_by_breaks_fast(x, breaks)
)
#> Unit: microseconds
#>                                                      expr   min     lq     mean  median
#>                                     table(cut(x, breaks)) 830.4 970.10 1219.004 1039.30
#>                      hist(x, breaks, plot = FALSE)$counts 393.2 421.10  497.680  466.85
#>   advr38pkg::count_by_breaks(x, breaks, use_outer = TRUE) 415.3 488.45  559.892  529.00
#>  advr38pkg::count_by_breaks(x, breaks, use_outer = FALSE) 214.9 237.10  283.447  250.85
#>                advr38pkg::count_by_breaks_fast(x, breaks) 188.7 205.45  244.874  219.90
#>       uq     max neval
#>  1218.50 12819.2   100
#>   518.45   991.3   100
#>   580.20  1079.1   100
#>   305.95   755.7   100
#>   256.90   548.8   100
x2 <- sample(100, size = 1e5, replace = TRUE)
breaks2 <- breaks * 10
breaks3 <- seq(0, 100, length.out = 100)
microbenchmark::microbenchmark(
  advr38pkg::count_by_breaks(x2, breaks2),
  advr38pkg::count_by_breaks_fast(x2, breaks2),
  advr38pkg::count_by_breaks(x2, breaks3),
  advr38pkg::count_by_breaks_fast(x2, breaks3)
)
#> Unit: milliseconds
#>                                          expr     min       lq      mean   median
#>       advr38pkg::count_by_breaks(x2, breaks2)  1.8981  1.95885  2.365524  2.03270
#>  advr38pkg::count_by_breaks_fast(x2, breaks2)  1.4139  1.60140  1.964527  1.68520
#>       advr38pkg::count_by_breaks(x2, breaks3) 37.9734 40.82125 45.215242 43.49680
#>  advr38pkg::count_by_breaks_fast(x2, breaks3)  1.5561  1.63235  2.095755  1.69485
#>        uq     max neval
#>   2.23370 12.8919   100
#>   1.82425 13.0252   100
#>  49.93595 61.3280   100
#>   1.84390 13.1516   100

An R user wants to implement some sampling on a sparse matrix and provides this working code:

N <- 2000
system.time({
  m <- Matrix::Matrix(0, nrow = N, ncol = N)
  for (j in 1:N) {
    cols <- sample((1:N)[-j], 2)  # pick 2 columns that are not j 
    m[j, cols] <- 1
  }
})
#>    user  system elapsed 
#>    1.71    0.03    1.78

This code is slow; can you find two major reasons why?

How can you more efficiently assign 1s? A faster solution would take:

#>    user  system elapsed 
#>    0.08    0.00    0.07

Can you use sampling with replacement (to avoid unnecessarily allocating memory) in this example? A faster solution would take:

#>    user  system elapsed 
#>    0.03    0.00    0.03

It would be even faster using Rcpp (cf. this SO answer).

Make a fast function that returns all prime numbers up to a number N.

N <- 1e6
system.time(
  primes <- advr38pkg::AllPrimesUpTo(N)
)
#>    user  system elapsed 
#>    0.05    0.00    0.04
plot(primes, pch = 20, cex = 0.5)

5.8 Parallel computing

I basically always use foreach and recommend to do so. See my guide to parallelism in R with foreach.

Just remember to optimize your code before trying to parallelize it.

Try to parallelize some of your best solutions for the previous exercises.