While some of the strategies covered here are specific to R, many are built on principles that can guide your coding in other languages.
-You should be able to work through this tutorial in any working R installation, including through RStudio. To work through it using R linked to a fast linear algebra package, you may want to use a virtual machine developed here at Berkeley, the Berkeley Common Environment (BCE). BCE is a virtual Linux machine - basically it is a Linux computer that you can run within your own computer, regardless of whether you are using Windows, Mac, or Linux. This provides a common environment so that things behave the same for all of us. However, note that BCE has not been updated in a while.
+You should be able to work through this tutorial in any working R installation, including through RStudio. However, in many cases the R you are using may not be linked to a fast linear algebra package.
This tutorial assumes you have a working knowledge of R.
@@ -223,7 +223,7 @@git clone https://github.com/berkeley-scf/tutorial-efficient-R
To create this HTML document, simply compile the corresponding R Markdown file in R as follows (the following will work from within BCE after cloning the repository as above).
+To create this HTML document, simply compile the corresponding R Markdown file in R as follows.
Rscript -e "library(knitr); knit2html('efficient-R.Rmd')"
Any calls to BLAS or to the LAPACK libraries that use BLAS to do higher-level linear algebra calculations will be nearly as fast as if you used C/C++ or Matlab, because R is using the compiled code from the BLAS and LAPACK libraries.
-In addition, the BLAS libraries above are threaded – they can use more than one core, and often will do so by default. More details in the tutorial on parallel programming.
+In addition, the BLAS libraries above are threaded – they can use more than one core, and often will do so by default. More details in the SCF tutorial on parallel programming.
## user system elapsed
-## 0.232 0.016 0.246
+## 0.251 0.011 0.262
system.time(rowMeans(x))
## user system elapsed
-## 0.024 0.000 0.024
+## 0.026 0.000 0.026
In general, user time gives the CPU time spent by R and system time gives the CPU time spent by the kernel (the operating system) on behalf of R. Operations that fall under system time include opening files, doing input or output, starting other processes, etc.
-To time code that runs very quickly, you may want to use the microbenchmark -package. Of course one would generally only care about such timing if a larger operation does the quick calculation very many times. Here's a comparison of different ways of accessing an element of a dataframe.
+To time code that runs very quickly, you should use the microbenchmark +package. Of course one would generally only care about accurately timing quick calculations if a larger operation does the quick calculation very many times. Here's a comparison of different ways of accessing an element of a dataframe.
library(microbenchmark)
df <- data.frame(vals = 1:3, labs = c('a','b','c'))
@@ -315,42 +315,45 @@ 3.1) Benchmarking
)
## Unit: microseconds
-## expr min lq mean median uq max neval cld
-## df[2, 1] 12.532 13.2260 13.79943 13.6360 13.946 18.980 100 b
-## df$vals[2] 6.731 7.9160 8.71348 8.3565 8.641 52.525 100 a
-## df[2, "vals"] 12.670 13.5355 14.36814 13.7850 14.136 44.447 100 b
+## Unit: nanoseconds
+## expr min lq mean median uq max neval cld
+## df[2, 1] 13976 14559.0 15169.09 14823 15000.0 27324 100 b
+## df$vals[2] 748 1063.0 1723.97 1591 1697.5 13343 100 a
+## df[2, "vals"] 13834 14413.5 15216.43 14621 14791.0 45620 100 b
-The rbenchmark package provides a nice wrapper function, benchmark,
-that automates timings and comparisons.
+microbenchmark is also a good way to compare slower calculations, e.g., doing a crossproduct via crossproduct() compared to “manually”:
-library(rbenchmark)
-# speed of one calculation
+library(microbenchmark)
n <- 1000
x <- matrix(rnorm(n^2), n)
-benchmark(crossprod(x), replications = 10,
- columns=c('test', 'elapsed', 'replications'))
+microbenchmark(
+ t(x) %*% x,
+ crossprod(x),
+ times = 10)
-## test elapsed replications
-## 1 crossprod(x) 0.481 10
+## Unit: milliseconds
+## expr min lq mean median uq max neval cld
+## t(x) %*% x 43.75747 43.79796 52.49518 48.89571 57.60102 78.42165 10 b
+## crossprod(x) 19.58690 19.67473 23.19242 19.82202 27.82647 33.67125 10 a
-# comparing different approaches to a task
-benchmark(
- {mns <- rep(NA, n); for(i in 1:n) mns[i] <- mean(x[i , ])},
- rowMeans(x),
- replications = 10,
- columns=c('test', 'elapsed', 'replications'))
+Challenge: What is inefficient about the manual approach above?
+
+An alternative that also automates timings and comparisons is the benchmark function from the rbenchmark package, but there's not really a reason to use it when microbenchmark is available.
+
+library(rbenchmark)
+# speed of one calculation
+benchmark(t(x) %*% x,
+ crossprod(x),
+ replications = 10,
+ columns=c('test', 'elapsed', 'replications'))
-## test
-## 1 {\n mns <- rep(NA, n)\n for (i in 1:n) mns[i] <- mean(x[i, ])\n}
-## 2 rowMeans(x)
-## elapsed replications
-## 1 0.196 10
-## 2 0.024 10
+## test elapsed replications
+## 2 crossprod(x) 0.275 10
+## 1 t(x) %*% x 0.382 10
In general, it's a good idea to repeat (replicate) your timing, as there is some stochasticity in how fast your computer will run a piece of code at any given moment.
@@ -362,73 +365,143 @@ 3.2) Profiling
The Rprof function will show you how much time is spent in
different functions, which can help you pinpoint bottlenecks in your
code. The output from Rprof can be hard to decipher, so you
-may want to use the proftools package functions, which make use of
-Rprof under the hood.
-
-Here's a function that works with a correlation matrix such as one
-might have for time series data. Basically, it creates a matrix of time lags (dd) and
-computes the correlation between the outcome for all pairs of times, based on the time
-lag between the pair of times. Then it computes the Cholesky factor of the correlation
-matrix so that it can generate a random time series in the last line. The question
-one might ask is which part(s) of the code take the most time.
-
-makeTS <- function(param, len){
- times <- seq(0, 1, length = len)
- dd <- rdist(times)
- C <- exp(-dd/param)
- U <- chol(C)
- white <- rnorm(len)
- return(crossprod(U, white))
-}
+will probably want to use the proftools package functions, which make use of
+Rprof under the hood.
+
+Here's a function that does the linear algebra to find the least squares solution in a linear regression, assuming x
+is the matrix of predictors, including a column for the intercept.
-## old approach:
-library(fields)
-if(FALSE) { # not running this, just for illustration
- Rprof("makeTS.prof", interval = 0.005, line.profiling = TRUE)
- out <- makeTS(0.1, 3000)
- Rprof(NULL)
- summaryRprof("makeTS.prof")
+lr_slow <- function(y, x) {
+ xtx <- t(x) %*% x
+ xty <- t(x) %*% y
+ inv <- solve(xtx) ## explicit matrix inverse is slow and generally a bad idea numerically
+ return(inv %*% xty)
}
+
+
+Let's run the function with profiling turned on.
+
+## generate random observations and random matrix of predictors
+y <- rnorm(5000)
+x <- matrix(rnorm(5000*1000), nrow = 5000)
-## using proftools instead:
library(proftools)
-pd <- profileExpr(makeTS(0.1, 3000))
-hotPaths(pd)
+
+pd1 <- profileExpr(lr_slow(y, x))
+hotPaths(pd1)
-Here's the result for the makeTS function:
+## path total.pct self.pct
+## lr_slow 100.00 0.00
+## . %*% (<text>:2) 33.65 33.65
+## . %*% (<text>:3) 1.92 1.92
+## . %*% (<text>:5) 0.96 0.96
+## . solve (<text>:4) 38.46 0.00
+## . . solve.default 38.46 37.50
+## . . . diag 0.96 0.96
+## . t (<text>:2) 1.92 0.00
+## . . t.default 1.92 1.92
+## . t (<text>:3) 23.08 0.00
+## . . t.default 23.08 23.08
+
- path total.pct self.pct
- makeTS 100.00 0.00
- . ?? (#File 1: :4) 56.06 56.06
- . chol (#File 1: :5) 30.30 0.00
- . . standardGeneric 30.30 0.00
- . . . chol 30.30 0.00
- . . . . chol.default 30.30 30.30
- . rdist (#File 1: :3) 12.12 0.00
- . . .Call 12.12 12.12
- . crossprod (#File 1: :7) 1.52 0.00
- . . crossprod 1.52 0.00
- . . . base::crossprod 1.52 1.52
+hotPaths(pd1, value = 'time')
-Note the nestedness of the results. For example, 12 percent of the time was spent in the call to rdist, of which essentially all of that was spent in .Call, which is a call out to C code.
+## path total.time self.time
+## lr_slow 2.08 0.00
+## . %*% (<text>:2) 0.70 0.70
+## . %*% (<text>:3) 0.04 0.04
+## . %*% (<text>:5) 0.02 0.02
+## . solve (<text>:4) 0.80 0.00
+## . . solve.default 0.80 0.78
+## . . . diag 0.02 0.02
+## . t (<text>:2) 0.04 0.00
+## . . t.default 0.04 0.04
+## . t (<text>:3) 0.48 0.00
+## . . t.default 0.48 0.48
+
-In this case, the results are not fully helpful, as 56 percent of the time is spent in other computations within makeTS that are not shown individually (see the “??” line).
+The first call to hotPaths shows the percentage of time spent in each call, while the second shows the actual time.
-As we increase the number of time points,
-the time taken up by the Cholesky would increase since that calculation
-is order of \(n^{3}\) while the others are order \(n^{2}\) (more in the Linear Algebra unit).
+Note the nestedness of the results. For example, essentially all the time spent in solve is actually spent in solve.default. In this case solve is just an S3 generic method that immediately calls the specific method solve.default.
-In this case, since the Cholesky and the main calculations in rdist, as well as exp,
-are all done in compiled C or Fortran code, there is probably not much we can do
-to speed this up (apart from using an optimized BLAS, which is essential). But in other cases profiling may reveal the slow steps in a piece of code.
+We can see that a lot of time is spent in doing the two crossproducts. So let's try using crossprod to make those steps faster.
-Note that Rprof works by sampling - every little while (the interval argument) during a calculation it finds out what function R is in and saves that information to the file given as the argument to Rprof. So if you try to profile code that finishes really quickly, there's not enough opportunity for the sampling to represent the calculation accurately and you may get spurious results.
+lr_medium <- function(y, x) {
+ xtx <- crossprod(x)
+ xty <- crossprod(x, y)
+ inv <- solve(xtx) ## explicit matrix inverse is slow and generally a bad idea numerically
+ return(inv %*% xty)
+}
+
+pd2 <- profileExpr(lr_medium(y, x))
+hotPaths(pd2)
+
+
+## path total.pct self.pct
+## lr_medium 100.00 0.00
+## . crossprod (<text>:2) 46.51 46.51
+## . crossprod (<text>:3) 4.65 4.65
+## . solve (<text>:4) 48.84 0.00
+## . . solve.default 48.84 46.51
+## . . . diag 2.33 2.33
+
+
+hotPaths(pd2, value = 'time')
+
+
+## path total.time self.time
+## lr_medium 0.86 0.00
+## . crossprod (<text>:2) 0.40 0.40
+## . crossprod (<text>:3) 0.04 0.04
+## . solve (<text>:4) 0.42 0.00
+## . . solve.default 0.42 0.40
+## . . . diag 0.02 0.02
+
+
+First note that this version takes about half the time of the previous one. Second note that a fair amount of time is spent computing the explicit matrix inverse using solve. (There's not much we can do to speed up the crossprod operations, apart from making sure we are using a fast BLAS and potentially a parallelized BLAS.) It's well known that one should avoid computing the explicit inverse if one can avoid it. Here's a faster version that avoids it.
+
+lr_fast <- function(y, x) {
+ xtx <- crossprod(x)
+ xty <- crossprod(x, y)
+ U <- chol(xtx)
+ tmp <- backsolve(U, xty, transpose = TRUE)
+ return(backsolve(U, tmp))
+}
+
+pd3 <- profileExpr(lr_fast(y, x))
+hotPaths(pd3)
+
+
+## path total.pct self.pct
+## lr_fast 100 0
+## . backsolve (<text>:6) 2 2
+## . chol (<text>:4) 12 0
+## . . chol.default 12 12
+## . crossprod (<text>:2) 80 80
+## . crossprod (<text>:3) 6 6
+
+
+hotPaths(pd3, value = 'time')
+
+
+## path total.time self.time
+## lr_fast 1.00 0.00
+## . backsolve (<text>:6) 0.02 0.02
+## . chol (<text>:4) 0.12 0.00
+## . . chol.default 0.12 0.12
+## . crossprod (<text>:2) 0.80 0.80
+## . crossprod (<text>:3) 0.06 0.06
+
+
+We can see we get another speedup from that final version of the code. (But beware my earlier caution that if comparing times between implementations, we should have replication.)
You might also check out profvis for an alternative to displaying profiling information
generated by Rprof.
+Note that Rprof works by sampling - every little while (the interval argument) during a calculation it finds out what function R is in and saves that information to the file given as the argument to Rprof. So if you try to profile code that finishes really quickly, there's not enough opportunity for the sampling to represent the calculation accurately and you may get spurious results.
+
Warning: Rprof conflicts with threaded linear algebra,
so you may need to set OMP_NUM_THREADS to 1 to disable threaded
linear algebra if you profile code that involves linear algebra.
@@ -451,34 +524,35 @@ 4.1) Pre-allocate memory
n <- 10000
z <- rnorm(n)
-fun1 <- function(vals) {
+fun_append <- function(vals) {
x <- exp(vals[1])
+ n <- length(vals)
for(i in 2:n) x <- c(x, exp(vals[i]))
return(x)
}
-fun2 <- function(vals) {
+fun_prealloc <- function(vals) {
n <- length(vals)
x <- rep(as.numeric(NA), n)
for(i in 1:n) x[i] <- exp(vals[i])
return(x)
}
-fun3 <- function(vals) {
+fun_vec <- function(vals) {
x <- exp(vals)
return(x)
}
-benchmark(fun1(z), fun2(z), fun3(z),
+benchmark(fun_append(z), fun_prealloc(z), fun_vec(z),
replications = 20, columns=c('test', 'elapsed', 'replications'))
-## test elapsed replications
-## 1 fun1(z) 2.347 20
-## 2 fun2(z) 0.021 20
-## 3 fun3(z) 0.006 20
+## test elapsed replications
+## 1 fun_append(z) 2.047 20
+## 2 fun_prealloc(z) 0.019 20
+## 3 fun_vec(z) 0.004 20
It's not necessary to use as.numeric above though it saves
a bit of time. Challenge: figure out why I have as.numeric(NA)
-and not just NA.
+and not just NA. Hint: what is the type of NA?
In some cases, we can speed up the initialization by initializing a vector of length one and then changing its length and/or dimension, although in many practical
circumstances this would be overkill.
@@ -497,8 +571,8 @@ 4.1) Pre-allocate memory
## 2 {\n x <- as.numeric(NA)\n length(x) <- nr * nc\n dim(x) <- c(nr, nc)\n}
## 1 x <- matrix(as.numeric(NA), nr, nc)
## elapsed replications
-## 2 0.130 10
-## 1 0.275 10
+## 2 0.055 10
+## 1 0.147 10
For lists, we can do this
@@ -586,7 +660,7 @@ 4.2) Vectorized calculations
## }
## .Internal(mean(x))
## }
-## <bytecode: 0x2cb7598>
+## <bytecode: 0x5646e07c69f0>
## <environment: namespace:base>
@@ -595,18 +669,21 @@ 4.2) Vectorized calculations
## function (x, pivot = FALSE, LINPACK = FALSE, tol = -1, ...)
## {
+## if (!missing(LINPACK))
+## stop("the LINPACK argument has been defunct since R 3.1.0")
## if (is.complex(x))
## stop("complex matrices not permitted at present")
## .Internal(La_chol(as.matrix(x), pivot, tol))
## }
-## <bytecode: 0x6a3f188>
+## <bytecode: 0x5646e34aa168>
## <environment: namespace:base>
Many R functions allow you to pass in vectors, and operate on those
vectors in vectorized fashion. So before writing a for loop, look
at the help information on the relevant function(s) to see if they
-operate in a vectorized fashion. Functions might take vectors for one or more of their arguments.
+operate in a vectorized fashion. Functions might take vectors for one or more of their arguments.
+Here we see that nchar is vectorized and that various arguments to substring can be vectors.
address <- c("Four score and seven years ago our fathers brought forth",
" on this continent, a new nation, conceived in Liberty, ",
@@ -628,8 +705,8 @@ 4.2) Vectorized calculations
substring(address[1], startIndices, startIndices + 1)
-## [1] "Fo" "r " "co" "e " "nd" "se" "en" "ye" "rs" "ag" " o" "r " "at" "er"
-## [15] " b" "ou" "ht" "fo" "th"
+## [1] "Fo" "r " "co" "e " "nd" "se" "en" "ye" "rs" "ag" " o" "r " "at" "er" " b"
+## [16] "ou" "ht" "fo" "th"
Challenge: Consider the chi-squared statistic involved in
@@ -685,8 +762,11 @@
4.2) Vectorized calculations
4.3) Using apply and specialized functions
-Another core efficiency strategy is to use the apply functionality.
-Even better than apply for calculating sums or means of columns
+
Historically, another core efficiency strategy in R has been to use the apply functionality (e.g., apply, sapply, lapply, mapply, etc.).
+
+Some faster alternatives to apply
+
+Note that even better than apply for calculating sums or means of columns
or rows (it also can be used for arrays) is {row,col}{Sums,Means}.
n <- 3000; x <- matrix(rnorm(n * n), nr = n)
@@ -697,8 +777,8 @@ 4.3) Using apply and specialized functions
## test elapsed replications
-## 1 out <- apply(x, 1, mean) 2.615 10
-## 2 out <- rowMeans(x) 0.220 10
+## 1 out <- apply(x, 1, mean) 1.562 10
+## 2 out <- rowMeans(x) 0.232 10
We can 'sweep' out a summary statistic, such as subtracting
@@ -708,7 +788,7 @@
4.3) Using apply and specialized functions
## user system elapsed
-## 0.124 0.040 0.162
+## 0.138 0.009 0.147
Here's a trick for doing the sweep based on vectorized calculations, remembering
@@ -720,7 +800,7 @@
4.3) Using apply and specialized functions
## user system elapsed
-## 0.276 0.048 0.324
+## 0.293 0.011 0.303
identical(out, out2)
@@ -731,9 +811,10 @@ 4.3) Using apply and specialized functions
Are apply, lapply, sapply, etc. faster than loops?
-Using apply with matrices and versions of apply with lists may or may not be faster
-than looping but generally produces cleaner code. Whether looping
-is slower will depend on whether a substantial part of the work is
+
Using apply with matrices and versions of apply with lists or vectors (e.g., lapply, sapply) may or may not be faster
+than looping but generally produces cleaner code.
+
+Whether looping and use of apply variants is slow will depend in part on whether a substantial part of the work is
in the overhead involved in the looping or in the time required by the function
evaluation on each of the elements. If you're worried about speed,
it's a good idea to benchmark the apply variant against looping.
@@ -753,7 +834,7 @@ Are apply, lapply, sapply, etc. faster than loops
## user system elapsed
-## 0.288 0.016 0.302
+## 0.087 0.000 0.086
system.time({
@@ -765,33 +846,42 @@ Are apply, lapply, sapply, etc. faster than loops
## user system elapsed
-## 0.300 0.016 0.312
+## 0.087 0.000 0.087
-And here's an example where sapply is much faster, because the core function evaluation at each iteration is very fast:
+And here's an example, where (unlike the previous example) the core computation is very fast, so we might expect the overhead of looping to be important. I believe that in old versions of R the sapply in this example was faster than looping in R, but that doesn't seem to be the case currently.
+I think this may be related to various somewhat recent improvements in R's handling of loops, possibly including the use of the byte compiler.
-z <- rnorm(10000)
-fun2 <- function(vals) {
+z <- rnorm(1e6)
+fun_loop <- function(vals) {
x <- as.numeric(NA)
- length(x) <- length(vals)
+ n <- length(vals)
+ length(x) <- n
for(i in 1:n) x[i] <- exp(vals[i])
return(x)
}
-fun4 <- function(vals) {
+
+fun_sapply <- function(vals) {
x <- sapply(vals, exp)
return(x)
}
-benchmark(fun2(z), fun4(z),
+fun_vec <- function(vals) {
+ x <- exp(vals)
+ return(x)
+}
+
+benchmark(fun_loop(z), fun_sapply(z), fun_vec(z),
replications = 10, columns=c('test', 'elapsed', 'replications'))
-## test elapsed replications
-## 1 fun2(z) 2.302 10
-## 2 fun4(z) 0.032 10
+## test elapsed replications
+## 1 fun_loop(z) 0.946 10
+## 2 fun_sapply(z) 4.792 10
+## 3 fun_vec(z) 0.113 10
-You'll notice if you look at the R code for lapply (sapply just calls lapply) that it calls directly out to C code, so the for loop is executed in compiled code.
+You'll notice if you look at the R code for lapply (sapply just calls lapply) that it calls directly out to C code, so the for loop is executed in compiled code. However, the code being executed at each iteration is still R code, so there is still all the overhead of the R interpreter.
print(lapply)
@@ -803,14 +893,16 @@ Are apply, lapply, sapply, etc. faster than loops
## X <- as.list(X)
## .Internal(lapply(X, FUN))
## }
-## <bytecode: 0xa2fd48>
+## <bytecode: 0x5646de5809b0>
## <environment: namespace:base>
4.4) Matrix algebra efficiency
Often calculations that are not explicitly linear algebra calculations
-can be done as matrix algebra. For example, we can sum the rows of a matrix by multiplying by a vector of ones. Given the extra computation involved in actually multiplying each number by one, it's surprising that this is faster than using R's heavily optimized rowSums function. It might be at least partially related to cache effects as colSums is more than twice as fast as rowSums in this case.
+can be done as matrix algebra. If our R installation has a fast (and possibly parallelized) BLAS, this allows our calculation to take advantage of it.
+
+For example, we can sum the rows of a matrix by multiplying by a vector of ones. Given the extra computation involved in actually multiplying each number by one, it's surprising that this is faster than using R's heavily optimized rowSums function.
mat <- matrix(rnorm(500*500), 500)
benchmark(apply(mat, 1, sum),
@@ -820,8 +912,8 @@ 4.4) Matrix algebra efficiency
## test elapsed replications
-## 1 apply(mat, 1, sum) 0.038 10
-## 2 mat %*% rep(1, ncol(mat)) 0.002 10
+## 1 apply(mat, 1, sum) 0.033 10
+## 2 mat %*% rep(1, ncol(mat)) 0.003 10
## 3 rowSums(mat) 0.006 10
@@ -834,6 +926,8 @@ 4.4) Matrix algebra efficiency
for loop is much faster than matrix multiplication. However,
there is a way to do it faster as matrix direct subtraction.
+Order of operations and efficiency
+
When doing matrix algebra, the order in which you do operations can
be critical for efficiency. How should I order the following calculation?
@@ -861,6 +955,8 @@ 4.4) Matrix algebra efficiency
Why is the second order much faster?
+Avoiding unnecessary operations
+
We can use the matrix direct product (i.e., A*B) to do
some manipulations much more quickly than using matrix multiplication.
Challenge: How can I use the direct product to find the trace
@@ -892,13 +988,13 @@
4.4) Matrix algebra efficiency
4.5) Fast mapping/lookup tables
Sometimes you need to map between two vectors. E.g.,
-\(y_{ij}\sim\mathcal{N}(\mu_{j},\sigma^{2})\)
-is a basic ANOVA type structure, where multiple observations in group \(j\)
-are associated with a common mean, \(\mu_j\).
+\(y_{i}\sim\mathcal{N}(\mu_{j[i]},\sigma^{2})\)
+is a basic ANOVA type structure, where multiple observations
+are associated with a common mean, \(\mu_j\), via the j[i] mapping.
How can we quickly look up the mean associated with each observation?
A good strategy is to create a vector, grp, that gives a numeric
-mapping of the observations to their cluster. Then you can access
+mapping of the observations to their cluster, playing the role of j[i] above. Then you can access
the \(\mu\) value relevant for each observation as: mus[grp]. This requires
that grp correctly map to the right elements of mus.
@@ -913,7 +1009,7 @@ 4.5) Fast mapping/lookup tables
numGrade = c(95, 85, 75),
fail = c(FALSE, FALSE, TRUE) )
grp <- match(df$clusterLabel, info$grade)
-df$numGrade = info$numGrade[grp]
+df$numGrade <- info$numGrade[grp]
df
@@ -925,17 +1021,26 @@ 4.5) Fast mapping/lookup tables
## 5 5 C 75
-R allows you to look up elements of vector by name.
-For example:
+Lookup by name versus index
+
+In the example above we looked up the mu values based on grp, which supplies the needed indexes as numeric indexes.
+
+R also allows you to look up elements of vector by name, as illustrated here by rearranging the code above a bit:
+
+info2 <- info$numGrade
+names(info2) <- info$grade
+info2
+
+
+## A B C
+## 95 85 75
+
-vals <- rnorm(10)
-names(vals) <- letters[1:10]
-select <- c("h", "h", "a", "c")
-vals[select]
+info2[df$clusterLabel]
-## h h a c
-## 0.2511293 0.2511293 -0.4466439 0.2898673
+## C B B A C
+## 75 85 85 95 75
You can do similar things in terms of looking up by name with dimension
@@ -943,84 +1048,83 @@
4.5) Fast mapping/lookup tables
named lists.
However, looking things up by name can be slow relative to looking up by index.
-Here's a toy example where we have a vector or list with a million elements and
+Here's a toy example where we have a vector or list with 1000 elements and
the character names of the elements are just the character versions of the
indices of the elements.
-n <- 1000000
+n <- 1000
x <- 1:n
xL <- as.list(x)
-nms <- as.character(x)
+nms <- paste0("var", as.character(x))
names(x) <- nms
names(xL) <- nms
-benchmark(
- x[500000], # index lookup in vector
- x["500000"], # name lookup in vector
- xL[[500000]], # index lookup in list
- xL[["500000"]], # name lookup in list
- replications = 10, columns=c('test', 'elapsed', 'replications'))
+x[1:3]
-## test elapsed replications
-## 2 x["500000"] 0.062 10
-## 1 x[5e+05] 0.000 10
-## 4 xL[["500000"]] 0.058 10
-## 3 xL[[5e+05]] 0.000 10
+## var1 var2 var3
+## 1 2 3
-Lookup by name is slow because R needs to scan through the objects
-one by one until it finds the one with the name it is looking for.
-In contrast, to look up by index, R can just go directly to the position of interest.
-
-In contrast, we can look up by name in an environment very quickly, because environments in R use hashing, which allows for fast lookup that does not require scanning through all of the names in the environment. In fact, this is how R itself looks for values when you specify variables in R code.
-
-xEnv <- as.environment(xL) # convert from a named list
-xEnv$"500000"
+xL[1:3]
-## [1] 500000
+## $var1
+## [1] 1
+##
+## $var2
+## [1] 2
+##
+## $var3
+## [1] 3
-# I need quotes above because numeric; otherwise xEnv$nameOfObject is fine
-xEnv[["500000"]]
+microbenchmark(
+ x[500], # index lookup in vector
+ x["var500"], # name lookup in vector
+ xL[[500]], # index lookup in list
+ xL[["var500"]]) # name lookup in list
-## [1] 500000
+## Unit: nanoseconds
+## expr min lq mean median uq max neval cld
+## x[500] 343 364 419.99 375.5 409.5 3035 100 a
+## x["var500"] 2655 2683 2954.80 2698.5 2750.0 17271 100 b
+## xL[[500]] 144 167 196.92 188.0 218.0 529 100 a
+## xL[["var500"]] 3696 3722 4074.27 3758.5 3788.5 11907 100 c
-benchmark(
- x[500000],
- xL[[500000]],
- xEnv[["500000"]],
- xEnv$"500000",
- replications = 10000, columns=c('test', 'elapsed', 'replications'))
+Lookup by name is slow because R needs to scan through the objects
+one by one until it finds the one with the name it is looking for.
+In contrast, to look up by index, R can just go directly to the position of interest.
+
+Side note: there is a lot of variability across the 100 replications shown above. This might have to do with cache effects (see next section).
+
+In contrast, we can look up by name in an environment very quickly, because environments in R use hashing, which allows for fast lookup that does not require scanning through all of the names in the environment. In fact, this is how R itself looks for values when you specify variables in R code, because the global environment, function frames, and package namespaces are all environments.
+
+xEnv <- as.environment(xL) # convert from a named list
+xEnv$var500
-## test elapsed replications
-## 1 x[5e+05] 0.026 10000
-## 3 xEnv[["500000"]] 0.030 10000
-## 4 xEnv$"500000" 0.031 10000
-## 2 xL[[5e+05]] 0.016 10000
+## [1] 500
microbenchmark(
- x[500000],
- xL[[500000]],
- xEnv[["500000"]],
- xEnv$"500000")
+ x[500],
+ xL[[500]],
+ xEnv[["var500"]],
+ xEnv$var500
+)
## Unit: nanoseconds
-## expr min lq mean median uq max neval cld
-## x[5e+05] 286 334.5 391.49 360.5 408 1881 100 a
-## xL[[5e+05]] 126 157.5 295.53 183.5 196 8172 100 a
-## xEnv[["500000"]] 1195 1221.5 1478.58 1239.0 1264 23034 100 b
-## xEnv$"500000" 1176 1218.5 1360.84 1258.5 1305 9762 100 b
+## expr min lq mean median uq max neval cld
+## x[500] 356 367.0 461.30 385.5 435.5 5604 100 b
+## xL[[500]] 155 170.5 226.74 197.0 220.0 2599 100 a
+## xEnv[["var500"]] 138 153.0 204.92 178.0 191.0 2490 100 a
+## xEnv$var500 201 218.0 263.41 232.5 259.5 1661 100 a
-The first benchmark indicates that lookup in an environment is nearly as fast as lookup by index in a vector or list, though the microbenchmark suggests that lookup in the environment is somewhat slower.
-
-Cache-aware programming
+4.6 Cache-aware programming
In addition to main memory (what we usually mean when we talk about RAM), computers also have memory caches, which are small amounts of fast memory that can be accessed very quickly by the processor. For example your computer might have L1, L2, and L3 caches, with L1 the smallest and fastest and L3 the largest and slowest. The idea is to try to have the data that is most used by the processor in the cache.
@@ -1032,38 +1136,53 @@ Cache-aware programming
If you know the size of the cache, you can try to design your code so that in a given part of your code you access data structures that will fit in the cache. This sort of thing is generally more relevant if you're coding in a language like C. But it can matter sometimes in R too. Here's an example:
-nr = 800000
-nc = 100
+nr <- 800000
+nc <- 100
## large matrix that won't fit in cache
-A = matrix(rnorm(nr * nc), nrow = nr)
-system.time(apply(A, 2, mean)) ## operate by column
+A <- matrix(rnorm(nr * nc), nrow = nr)
+tA <- t(A)
+benchmark(
+ apply(A, 2, mean), ## operate by column
+ apply(tA, 1, mean), ## exact same calculation, but by row
+ replications = 10, columns=c('test', 'elapsed', 'replications'))
-## user system elapsed
-## 1.116 0.316 1.430
+## test elapsed replications
+## 1 apply(A, 2, mean) 10.675 10
+## 2 apply(tA, 1, mean) 21.352 10
-A = t(A)
-system.time(apply(A, 1, mean)) ## same calculation, but by row
+Now let's compare things when we make the matrix small enough that it fits in the cache. In this case it fits into the largest (L3) cache but not the smaller (L2 and L1) caches, and we see that the difference in speed disappears.
+
+nr <- 800
+nc <- 100
+## small matrix that should fit in cache
+A <- matrix(rnorm(nr * nc), nrow = nr)
+## Yep, the size is less than the L3 cache:
+object.size(A)
-## user system elapsed
-## 1.804 0.272 2.079
+## 640216 bytes
-nr = 800
-nc = 100
-## small matrix that should fit in cache
-A = matrix(rnorm(nr * nc), nrow = nr)
-tA = t(A)
+memuse::Sys.cachesize()
+
+
+## L1I: 32.000 KiB
+## L1D: 32.000 KiB
+## L2: 256.000 KiB
+## L3: 8.000 MiB
+
+
+tA <- t(A)
benchmark(apply(A, 2, mean), ## by column
apply(tA, 1, mean), ## by row
replications = 10, columns=c('test', 'elapsed', 'replications'))
## test elapsed replications
-## 1 apply(A, 2, mean) 0.014 10
-## 2 apply(tA, 1, mean) 0.015 10
+## 1 apply(A, 2, mean) 0.013 10
+## 2 apply(tA, 1, mean) 0.012 10