Efficient R code
The goal of this section is to highlight practices for writing correct, robust and efficient R code.
Priorities
- Correct: consistent with hand-worked examples (
identical(),all.equal()) - Robust: supports realistic inputs, e.g., 0-length vectors,
NAvalues, … - Simple: easy to understand next month; easy to describe what it does to a colleague; easy to spot logical errors; easy to enhance.
- Fast, or at least reasonable given the speed of modern computers.
Strategies
- Profile
- Look at the script to understand in general terms what it is doing.
- Step through the code to see how it is executed, and to gain an understanding of the speed of each line.
- Time evaluation of select lines or simple chunks of code with
system.time()or the microbenchmark package. - Profile the code with a tool that indicates how much time is spent in each function call or line – the built-in
Rprof()function, or packages such as lineprof or aprof
Vectorize – operate on vectors, rather than explicit loops
x <- 1:10 log(x) ## NOT for (i in seq_along) x[i] <- log(x[i])## [1] 0.0000000 0.6931472 1.0986123 1.3862944 1.6094379 1.7917595 1.9459101 ## [8] 2.0794415 2.1972246 2.3025851Pre-allocate memory, then fill in the result
result <- numeric(10) result[1] <- runif(1) for (i in 2:length(result)) result[i] <- runif(1) * result[i - 1] result## [1] 4.376985e-01 1.058143e-01 2.319046e-02 5.399646e-03 5.275794e-03 ## [6] 1.391109e-03 1.085354e-04 8.256891e-05 6.167133e-05 4.673066e-05- Hoist common sub-expressions outside of repeated calculations, so that the sub-expression is only calculated once
- Simple, e.g., ‘hoist’ constant multiplications from a
forloop - Higher-level, e.g., use
lm.fit()rather than repeatedly fitting the same design matrix.
- Re-use existing, tested code
- Efficient implementations of common operations –
tabulate(),rowSums()and friends,%in%, … - Efficient domain-specific implementations, e.g., snpStats for GWAS linear models; limma for microarray linear models; edgeR, DESeq2 for negative binomial GLMs relevant to RNASeq.
- Reuse others’ work – DESeq2, GenomicRanges, Biostrings, …, dplyr, data.table, Rcpp
Case Study: Pre-allocate and vectorize
Here’s an obviously inefficient function:
f0 <- function(n, a=2) {
## stopifnot(is.integer(n) && (length(n) == 1) &&
## !is.na(n) && (n > 0))
result <- numeric()
for (i in seq_len(n))
result[[i]] <- a * log(i)
result
}Use system.time() to investigate how this algorithm scales with n, focusing on elapsed time.
system.time(f0(10000))## user system elapsed
## 0.139 0.000 0.139n <- 1000 * seq(1, 20, 2)
t <- sapply(n, function(i) system.time(f0(i))[[3]])
plot(t ~ n, type="b")
Remember the current ‘correct’ value, and an approximate time
n <- 10000
system.time(expected <- f0(n))## user system elapsed
## 0.155 0.000 0.156head(expected)## [1] 0.000000 1.386294 2.197225 2.772589 3.218876 3.583519Revise the function to hoist the common multiplier, a, out of the loop. Make sure the result of the ‘optimization’ and the original calculation are the same. Use the microbenchmark package to compare the two versions
f1 <- function(n, a=2) {
result <- numeric()
for (i in seq_len(n))
result[[i]] <- log(i)
a * result
}
identical(expected, f1(n))## [1] TRUElibrary(microbenchmark)
microbenchmark(f0(n), f1(n), times=5)## Unit: milliseconds
## expr min lq mean median uq max neval
## f0(n) 118.6033 123.8989 137.1136 132.7924 132.8397 177.4337 5
## f1(n) 129.6826 175.5757 167.7663 175.9782 176.8241 180.7711 5Adopt a ‘pre-allocate and fill’ strategy
f2 <- function(n, a=2) {
result <- numeric(n)
for (i in seq_len(n))
result[[i]] <- log(i)
a * result
}
identical(expected, f2(n))## [1] TRUEmicrobenchmark(f0(n), f2(n), times=5)## Unit: milliseconds
## expr min lq mean median uq max neval
## f0(n) 136.82808 180.08703 172.27453 181.08084 181.26852 182.1082 5
## f2(n) 18.75233 18.85918 20.29981 20.28972 20.37353 23.2243 5Use an *apply() function to avoid having to explicitly pre-allocate, and make opportunities for vectorization more apparent.
f3 <- function(n, a=2)
a * sapply(seq_len(n), log)
identical(expected, f3(n))## [1] TRUEmicrobenchmark(f0(n), f2(n), f3(n), times=10)## Unit: milliseconds
## expr min lq mean median uq max
## f0(n) 134.451852 137.003841 171.20652 183.653188 188.31813 191.75448
## f2(n) 18.745654 18.833771 20.50144 20.254524 20.88116 23.44838
## f3(n) 9.467853 9.651863 10.05424 9.729817 10.56676 11.24802
## neval
## 10
## 10
## 10Now that the code is presented in a single line, it is apparent that it could be easily vectorized. Seize the opportunity to vectorize it:
f4 <- function(n, a=2)
a * log(seq_len(n))
identical(expected, f4(n))## [1] TRUEmicrobenchmark(f0(n), f3(n), f4(n), times=10)## Unit: microseconds
## expr min lq mean median uq max
## f0(n) 135787.889 137468.593 170961.713 184137.111 187689.288 188119.542
## f3(n) 8705.995 9402.749 9507.574 9554.713 9774.741 9841.998
## f4(n) 517.837 522.898 524.170 525.450 526.068 526.819
## neval
## 10
## 10
## 10f4() definitely seems to be the winner. How does it scale with n? (Repeat several times)
n <- 10 ^ (5:8) # 100x larger than f0
t <- sapply(n, function(i) system.time(f4(i))[[3]])
plot(t ~ n, log="xy", type="b")
Any explanations for the different pattern of response?
Lessons learned:
- Vectorizing offers a huge improvement over iteration
- Pre-allocate-and-fill is very helpful when explicit iteration is required.
*apply()functions help avoid need for explicit pre-allocation and make opportunities for vectorization more apparent. This may come at a small performance cost, but is generally worth it- Hoisting common sub-expressions can be helpful for improving performance when explicit iteration is required.