From 884eae89814b0cb28c2ec8e09318ed45cb19ddf3 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Tue, 22 Sep 2020 16:56:33 -0700 Subject: [PATCH 1/6] remove BCE reference --- efficient-R.Rmd | 4 ++-- efficient-R.html | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/efficient-R.Rmd b/efficient-R.Rmd index 8956f7c..4fda01f 100644 --- a/efficient-R.Rmd +++ b/efficient-R.Rmd @@ -15,7 +15,7 @@ This tutorial covers strategies for writing efficient R code by taking advantage 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)](http://bce.berkeley.edu). 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. @@ -24,7 +24,7 @@ Materials for this tutorial, including the R markdown file and associated code f 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. ```{r, build-html, eval=FALSE} Rscript -e "library(knitr); knit2html('efficient-R.Rmd')" ``` diff --git a/efficient-R.html b/efficient-R.html index 078b341..82c571c 100644 --- a/efficient-R.html +++ b/efficient-R.html @@ -214,7 +214,7 @@

0) This Tutorial

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 @@

0) This Tutorial

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')"
From 75e99ebc1122b3d252d2b92b43bca174eaf31ee6 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Wed, 23 Sep 2020 12:40:22 -0700 Subject: [PATCH 2/6] update benchmark/profiling/timing of name lookup --- efficient-R.R | 129 +++++++++++------- efficient-R.Rmd | 73 +++++----- efficient-R.html | 345 +++++++++++++++++++++++++++-------------------- 3 files changed, 314 insertions(+), 233 deletions(-) diff --git a/efficient-R.R b/efficient-R.R index 11ed1f2..722759d 100644 --- a/efficient-R.R +++ b/efficient-R.R @@ -17,44 +17,84 @@ microbenchmark( df[2, 'vals'] ) +## @knitr microbenchmark2 +library(microbenchmark) +n <- 1000 +x <- matrix(rnorm(n^2), n) +microbenchmark( + t(x) %*% x, + crossprod(x), + times = 10) + ## @knitr benchmark library(rbenchmark) # speed of one calculation -n <- 1000 -x <- matrix(rnorm(n^2), n) -benchmark(crossprod(x), replications = 10, +benchmark(t(x) %*% x, + crossprod(x), + replications = 10, columns=c('test', 'elapsed', 'replications')) -# 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')) - - -## @knitr Rprof -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)) + + + +## @knitr Rprof-fun +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) +} + +## @knitr Rprof-run1 + +## generate random observations and random matrix of predictors +y <- rnorm(5000) +x <- matrix(rnorm(5000*1000), nrow = 5000) + +library(proftools) + +pd1 <- profileExpr(lr_slow(y, x)) +hotPaths(pd1) +hotPaths(pd1, value = 'time') + + +## @knitr Rprof-run2 +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) +hotPaths(pd2, value = 'time') + + +## @knitr Rprof-run3 +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) +hotPaths(pd3, value = 'time') + + +## @knitr Rprof-old + ## 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") -} -## using proftools instead: -library(proftools) -pd <- profileExpr(makeTS(0.1, 3000)) -hotPaths(pd) +Rprof("makeRegr.prof", interval = 0.005, line.profiling = TRUE) +out <- lr_slow(y, x) +Rprof(NULL) +summaryRprof("makeRegr.prof") + + ## @knitr preallocate @@ -223,36 +263,29 @@ select <- c("h", "h", "a", "c") vals[select] ## @knitr index-lookup -n <- 1000000 +n <- 1000 x <- 1:n xL <- as.list(x) nms <- 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')) +microbenchmark( + x[500], # index lookup in vector + x["500"], # name lookup in vector + xL[[500]], # index lookup in list + xL[["500"]]) # name lookup in list ## @knitr env-lookup xEnv <- as.environment(xL) # convert from a named list -xEnv$"500000" +xEnv$"500" # I need quotes above because numeric; otherwise xEnv$nameOfObject is fine -xEnv[["500000"]] -benchmark( - x[500000], - xL[[500000]], - xEnv[["500000"]], - xEnv$"500000", - replications = 10000, columns=c('test', 'elapsed', 'replications')) +xEnv[["500"]] microbenchmark( - x[500000], - xL[[500000]], - xEnv[["500000"]], - xEnv$"500000") + x[500], + xL[[500]], + xEnv[["500"]], + xEnv$"500") ## @knitr cache-aware diff --git a/efficient-R.Rmd b/efficient-R.Rmd index 4fda01f..9e8a79c 100644 --- a/efficient-R.Rmd +++ b/efficient-R.Rmd @@ -87,15 +87,18 @@ implementations. Here's a basic comparison of the time to calculate the row mean 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. ```{r, microbenchmark} ``` +*microbenchmark* is also a good way to compare slower calculations, e.g., doing a crossproduct via *crossproduct()* compared to "manually": -The *rbenchmark* package provides a nice wrapper function, *benchmark*, -that automates timings and comparisons. +```{r, microbenchmark2} +``` + +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. ```{r, benchmark} ``` @@ -106,58 +109,46 @@ You might also checkout the *tictoc* package. ## 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. +will probably 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. +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. -```{r, Rprof, eval=FALSE} -``` -Here's the result for the *makeTS* function: +```{r, Rprof-fun} ``` - 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 + +Let's run the function with profiling turned on. + +```{r, Rprof-run1} ``` +The first call to *hotPaths* shows the percentage of time spent in each call, while the second shows the actual time. + +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 a generic method that immediately calls *solve.default*. -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. +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. -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). -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). +```{r, Rprof-run2} +``` -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. +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.) 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. -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. +```{r, Rprof-run3} +``` + +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. @@ -410,7 +401,7 @@ names of matrices/arrays, row and column names of dataframes, and 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. @@ -421,15 +412,13 @@ 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. ```{r, env-lookup, cache=TRUE} ``` -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 diff --git a/efficient-R.html b/efficient-R.html index 82c571c..4e8a87e 100644 --- a/efficient-R.html +++ b/efficient-R.html @@ -214,7 +214,7 @@

0) This Tutorial

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. However, in many cases the R you are using may not be linked to a fast linear algebra package.

+

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.

@@ -291,20 +291,20 @@

3.1) Benchmarking

##    user  system elapsed 
-##   0.232   0.016   0.246
+##   0.231   0.020   0.251
 
system.time(rowMeans(x))
 
##    user  system elapsed 
-##   0.024   0.000   0.024
+##   0.025   0.000   0.024

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'))
@@ -316,41 +316,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
+##           expr    min      lq     mean  median      uq    max neval cld
+##       df[2, 1] 13.570 14.0245 14.55557 14.1945 14.4675 29.628   100   b
+##     df$vals[2]  1.033  1.4335  1.85085  1.9300  2.0540  7.692   100  a 
+##  df[2, "vals"] 13.529 13.8460 14.96311 14.1025 14.3445 79.596   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
+##    t(x) %*% x 38.5048 44.10347 52.56555 46.44601 57.27781 74.89820    10
+##  crossprod(x) 19.0293 19.07104 23.59746 20.53333 28.25155 32.17322    10
+##  cld
+##    b
+##   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'))
+
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.238           10
+## 1   t(x) %*% x   0.564           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 +366,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.

 
-## 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")
+
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.

+
+
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)
+

+
+
##  path               total.pct self.pct
+##  lr_slow            100.00     0.00   
+##  . %*% (<text>:2)    44.32    44.32   
+##  . t (<text>:3)      29.55     0.00   
+##  . . t.default       29.55    29.55   
+##  . solve (<text>:4)  22.73     0.00   
+##  . . solve.default   22.73    21.59   
+##  . . . diag           1.14     1.14   
+##  . t (<text>:2)       2.27     0.00   
+##  . . t.default        2.27     2.27   
+##  . %*% (<text>:3)     1.14     1.14
 

 
-
Here's the result for the makeTS function:

+
hotPaths(pd1, value = 'time')
+

 
-
 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   
+
##  path               total.time self.time
+##  lr_slow            1.76       0.00     
+##  . %*% (<text>:2)   0.78       0.78     
+##  . t (<text>:3)     0.52       0.00     
+##  . . t.default      0.52       0.52     
+##  . solve (<text>:4) 0.40       0.00     
+##  . . solve.default  0.40       0.38     
+##  . . . diag         0.02       0.02     
+##  . t (<text>:2)     0.04       0.00     
+##  . . t.default      0.04       0.04     
+##  . %*% (<text>:3)   0.02       0.02
 

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

+
The first call to hotPaths shows the percentage of time spent in each call, while the second shows the actual time.

 
-
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). 

+
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 a generic method that immediately calls solve.default.

 
-
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).

+
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.

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

+
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)
+}                   
 
-
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.

+pd2 <- profileExpr(lr_medium(y, x))
+hotPaths(pd2)
+

+
+
##  path                   total.pct self.pct
+##  lr_medium              100.00     0.00   
+##  . crossprod (<text>:2)  50.94    50.94   
+##  . solve (<text>:4)      35.85     0.00   
+##  . . solve.default       35.85    33.96   
+##  . . . diag               1.89     1.89   
+##  . %*% (<text>:5)         9.43     9.43   
+##  . crossprod (<text>:3)   3.77     3.77
+

+
+
hotPaths(pd2, value = 'time')
+

+
+
##  path                   total.time self.time
+##  lr_medium              1.06       0.00     
+##  . crossprod (<text>:2) 0.54       0.54     
+##  . solve (<text>:4)     0.38       0.00     
+##  . . solve.default      0.38       0.36     
+##  . . . diag             0.02       0.02     
+##  . %*% (<text>:5)       0.10       0.10     
+##  . crossprod (<text>:3) 0.04       0.04
+

+
+
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.) 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.00     0.00   
+##  . crossprod (<text>:2)  78.79    78.79   
+##  . chol (<text>:4)       12.12     0.00   
+##  . . chol.default        12.12    12.12   
+##  . crossprod (<text>:3)   6.06     6.06   
+##  . backsolve (<text>:6)   3.03     3.03
+

+
+
hotPaths(pd3, value = 'time')
+

+
+
##  path                   total.time self.time
+##  lr_fast                0.66       0.00     
+##  . crossprod (<text>:2) 0.52       0.52     
+##  . chol (<text>:4)      0.08       0.00     
+##  . . chol.default       0.08       0.08     
+##  . crossprod (<text>:3) 0.04       0.04     
+##  . backsolve (<text>:6) 0.02       0.02
+

+
+
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. 

@@ -471,8 +545,8 @@ 
4.1) Pre-allocate memory

 

 
 
##      test elapsed replications
-## 1 fun1(z)   2.347           20
-## 2 fun2(z)   0.021           20
+## 1 fun1(z)   2.104           20
+## 2 fun2(z)   0.026           20
 ## 3 fun3(z)   0.006           20
 

 
@@ -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.158           10
 

 
 
For lists, we can do this

@@ -586,7 +660,7 @@ 
4.2) Vectorized calculations

 ##     }
 ##     .Internal(mean(x))
 ## }
-## <bytecode: 0x2cb7598>
+## <bytecode: 0x55bf04f01630>
 ## <environment: namespace:base>
 

 
@@ -599,7 +673,7 @@ 
4.2) Vectorized calculations

 ##         stop("complex matrices not permitted at present")
 ##     .Internal(La_chol(as.matrix(x), pivot, tol))
 ## }
-## <bytecode: 0x6a3f188>
+## <bytecode: 0x55bf051f0d18>
 ## <environment: namespace:base>
 

 
@@ -697,8 +771,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.880           10
+## 2       out <- rowMeans(x)   0.212           10
 

 
 
We can 'sweep' out a summary statistic, such as subtracting
@@ -708,7 +782,7 @@ 
4.3) Using apply and specialized functions

 
##    user  system elapsed 
-##   0.124   0.040   0.162
+##   0.121   0.024   0.144

Here's a trick for doing the sweep based on vectorized calculations, remembering @@ -720,7 +794,7 @@

4.3) Using apply and specialized functions

##    user  system elapsed 
-##   0.276   0.048   0.324
+##   0.304   0.021   0.325
 
identical(out, out2)
@@ -753,7 +827,7 @@ 
Are apply, lapply, sapply, etc. faster than loops
 
##    user  system elapsed 
-##   0.288   0.016   0.302
+##   0.085   0.000   0.084
 
system.time({
@@ -765,7 +839,7 @@ 
Are apply, lapply, sapply, etc. faster than loops
 
##    user  system elapsed 
-##   0.300   0.016   0.312
+##   0.085   0.000   0.084

And here's an example where sapply is much faster, because the core function evaluation at each iteration is very fast:

@@ -787,8 +861,8 @@

Are apply, lapply, sapply, etc. faster than loops 
##      test elapsed replications
-## 1 fun2(z)   2.302           10
-## 2 fun4(z)   0.032           10
+## 1 fun2(z)   2.258           10
+## 2 fun4(z)   0.034           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.

@@ -803,7 +877,7 @@

Are apply, lapply, sapply, etc. faster than loops ## X <- as.list(X) ## .Internal(lapply(X, FUN)) ## } -## <bytecode: 0xa2fd48> +## <bytecode: 0x55bf0309a0d0> ## <environment: namespace:base> @@ -820,8 +894,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.032           10
+## 2 mat %*% rep(1, ncol(mat))   0.003           10
 ## 3              rowSums(mat)   0.006           10
@@ -943,83 +1017,68 @@

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)
 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'))
+microbenchmark(
+    x[500],  # index lookup in vector
+    x["500"], # name lookup in vector
+    xL[[500]], # index lookup in list
+    xL[["500"]]) # name lookup in list
 

 
-
##             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
+
## Unit: nanoseconds
+##         expr   min      lq     mean  median      uq    max neval cld
+##       x[500]   475   523.5   587.95   552.5   605.5   1278   100 a  
+##     x["500"]  6202  6269.5  7407.16  6305.5  6395.0 105028   100  b 
+##    xL[[500]]   185   219.5   394.80   251.0   273.5   8283   100 a  
+##  xL[["500"]] 12043 12108.5 12432.33 12149.5 12213.0  20153   100   c
 

 
 
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. 

 
 
xEnv <- as.environment(xL)  # convert from a named list
-xEnv$"500000"  
+xEnv$"500"  
 

 
-
## [1] 500000
+
## [1] 500
 

 
 
# I need quotes above because numeric; otherwise xEnv$nameOfObject is fine
-xEnv[["500000"]]
-

-
-
## [1] 500000
+xEnv[["500"]]
 

 
-
benchmark(
-  x[500000],
-  xL[[500000]],
-  xEnv[["500000"]],
-  xEnv$"500000",
-  replications = 10000, columns=c('test', 'elapsed', 'replications'))
-

-
-
##               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[["500"]],
+  xEnv$"500")
 

 
 
## 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] 466 520.5 687.27  540.5 563.0 13909   100   b
+##      xL[[500]] 168 190.5 252.75  218.5 245.0  3511   100  a 
+##  xEnv[["500"]] 159 170.5 228.14  195.5 219.5  3054   100  a 
+##     xEnv$"500" 199 222.5 297.29  277.0 298.0  2763   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

 
 
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. 

@@ -1040,7 +1099,7 @@ 
Cache-aware programming

 

 
 
##    user  system elapsed 
-##   1.116   0.316   1.430
+##   0.856   0.224   1.081
 

 
 
A = t(A)
@@ -1048,7 +1107,7 @@ 
Cache-aware programming

 

 
 
##    user  system elapsed 
-##   1.804   0.272   2.079
+##   1.751   0.296   2.047
 

 
 
nr = 800
@@ -1062,8 +1121,8 @@ 
Cache-aware programming

 

 
 
##                 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.012           10
+## 2 apply(tA, 1, mean)   0.012           10
 

 
 

From 17f75d2a382db0ce2ff6ca865803b264104c8dd6 Mon Sep 17 00:00:00 2001
From: Christopher Paciorek 
Date: Thu, 24 Sep 2020 07:50:03 -0700
Subject: [PATCH 3/6] add missing section number

---
 efficient-R.Rmd  |   2 +-
 efficient-R.html | 148 ++++++++++++++++++++++-------------------------
 2 files changed, 70 insertions(+), 80 deletions(-)

diff --git a/efficient-R.Rmd b/efficient-R.Rmd
index 9e8a79c..b660e68 100644
--- a/efficient-R.Rmd
+++ b/efficient-R.Rmd
@@ -420,7 +420,7 @@ In contrast, we can look up by name in an environment very quickly, because envi
 ```
 
 
-## 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. 
 
diff --git a/efficient-R.html b/efficient-R.html
index 4e8a87e..63b6e0d 100644
--- a/efficient-R.html
+++ b/efficient-R.html
@@ -291,14 +291,14 @@ 
3.1) Benchmarking

 

 
 
##    user  system elapsed 
-##   0.231   0.020   0.251
+##   0.236   0.012   0.248
 

 
 
system.time(rowMeans(x))
 

 
 
##    user  system elapsed 
-##   0.025   0.000   0.024
+##   0.023   0.000   0.023
 

 
 
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.

@@ -315,11 +315,11 @@ 
3.1) Benchmarking

 )
 

 
-
## Unit: microseconds
-##           expr    min      lq     mean  median      uq    max neval cld
-##       df[2, 1] 13.570 14.0245 14.55557 14.1945 14.4675 29.628   100   b
-##     df$vals[2]  1.033  1.4335  1.85085  1.9300  2.0540  7.692   100  a 
-##  df[2, "vals"] 13.529 13.8460 14.96311 14.1025 14.3445 79.596   100   b
+
## Unit: nanoseconds
+##           expr   min      lq     mean median      uq   max neval cld
+##       df[2, 1] 13869 14563.0 15817.02  14748 14964.0 93003   100   b
+##     df$vals[2]   997  1625.5  2306.03   2049  2162.5 27105   100  a 
+##  df[2, "vals"] 13792 14334.5 14845.36  14501 14864.0 24041   100   b
 

 
 
microbenchmark is also a good way to compare slower calculations, e.g., doing a crossproduct via crossproduct() compared to “manually”:

@@ -334,9 +334,9 @@ 
3.1) Benchmarking

 

 
 
## Unit: milliseconds
-##          expr     min       lq     mean   median       uq      max neval
-##    t(x) %*% x 38.5048 44.10347 52.56555 46.44601 57.27781 74.89820    10
-##  crossprod(x) 19.0293 19.07104 23.59746 20.53333 28.25155 32.17322    10
+##          expr      min       lq     mean   median       uq      max neval
+##    t(x) %*% x 87.71254 88.20979 89.55519 88.74442 89.52175 97.70361    10
+##  crossprod(x) 43.63842 44.12943 44.76383 44.61579 45.39003 46.04143    10
 ##  cld
 ##    b
 ##   a
@@ -353,8 +353,8 @@ 
3.1) Benchmarking

 

 
 
##           test elapsed replications
-## 2 crossprod(x)   0.238           10
-## 1   t(x) %*% x   0.564           10
+## 2 crossprod(x)   0.428           10
+## 1   t(x) %*% x   0.863           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.

@@ -394,31 +394,27 @@ 
3.2) Profiling

 
 
##  path               total.pct self.pct
 ##  lr_slow            100.00     0.00   
-##  . %*% (<text>:2)    44.32    44.32   
-##  . t (<text>:3)      29.55     0.00   
-##  . . t.default       29.55    29.55   
-##  . solve (<text>:4)  22.73     0.00   
-##  . . solve.default   22.73    21.59   
-##  . . . diag           1.14     1.14   
-##  . t (<text>:2)       2.27     0.00   
-##  . . t.default        2.27     2.27   
-##  . %*% (<text>:3)     1.14     1.14
+##  . %*% (<text>:2)    54.05    54.05   
+##  . t (<text>:3)      21.62     0.00   
+##  . . t.default       21.62    21.62   
+##  . solve (<text>:4)  18.92     0.00   
+##  . . solve.default   18.92    18.92   
+##  . t (<text>:2)       5.41     0.00   
+##  . . t.default        5.41     5.41
 

 
 
hotPaths(pd1, value = 'time')
 

 
 
##  path               total.time self.time
-##  lr_slow            1.76       0.00     
-##  . %*% (<text>:2)   0.78       0.78     
-##  . t (<text>:3)     0.52       0.00     
-##  . . t.default      0.52       0.52     
-##  . solve (<text>:4) 0.40       0.00     
-##  . . solve.default  0.40       0.38     
-##  . . . diag         0.02       0.02     
+##  lr_slow            0.74       0.00     
+##  . %*% (<text>:2)   0.40       0.40     
+##  . t (<text>:3)     0.16       0.00     
+##  . . t.default      0.16       0.16     
+##  . solve (<text>:4) 0.14       0.00     
+##  . . solve.default  0.14       0.14     
 ##  . t (<text>:2)     0.04       0.00     
-##  . . t.default      0.04       0.04     
-##  . %*% (<text>:3)   0.02       0.02
+##  . . t.default      0.04       0.04
 

 
 
The first call to hotPaths shows the percentage of time spent in each call, while the second shows the actual time.

@@ -440,25 +436,21 @@ 
3.2) Profiling

 
 
##  path                   total.pct self.pct
 ##  lr_medium              100.00     0.00   
-##  . crossprod (<text>:2)  50.94    50.94   
-##  . solve (<text>:4)      35.85     0.00   
-##  . . solve.default       35.85    33.96   
-##  . . . diag               1.89     1.89   
-##  . %*% (<text>:5)         9.43     9.43   
-##  . crossprod (<text>:3)   3.77     3.77
+##  . crossprod (<text>:2)  56.25    56.25   
+##  . solve (<text>:4)      37.50     0.00   
+##  . . solve.default       37.50    37.50   
+##  . crossprod (<text>:3)   6.25     6.25
 

 
 
hotPaths(pd2, value = 'time')
 

 
 
##  path                   total.time self.time
-##  lr_medium              1.06       0.00     
-##  . crossprod (<text>:2) 0.54       0.54     
-##  . solve (<text>:4)     0.38       0.00     
-##  . . solve.default      0.38       0.36     
-##  . . . diag             0.02       0.02     
-##  . %*% (<text>:5)       0.10       0.10     
-##  . crossprod (<text>:3) 0.04       0.04
+##  lr_medium              0.32       0.00     
+##  . crossprod (<text>:2) 0.18       0.18     
+##  . solve (<text>:4)     0.12       0.00     
+##  . . solve.default      0.12       0.12     
+##  . crossprod (<text>:3) 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.) 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.

@@ -477,23 +469,21 @@ 
3.2) Profiling

 
 
##  path                   total.pct self.pct
 ##  lr_fast                100.00     0.00   
-##  . crossprod (<text>:2)  78.79    78.79   
-##  . chol (<text>:4)       12.12     0.00   
-##  . . chol.default        12.12    12.12   
-##  . crossprod (<text>:3)   6.06     6.06   
-##  . backsolve (<text>:6)   3.03     3.03
+##  . crossprod (<text>:2)  81.82    81.82   
+##  . chol (<text>:4)        9.09     0.00   
+##  . . chol.default         9.09     9.09   
+##  . crossprod (<text>:3)   9.09     9.09
 

 
 
hotPaths(pd3, value = 'time')
 

 
 
##  path                   total.time self.time
-##  lr_fast                0.66       0.00     
-##  . crossprod (<text>:2) 0.52       0.52     
-##  . chol (<text>:4)      0.08       0.00     
-##  . . chol.default       0.08       0.08     
-##  . crossprod (<text>:3) 0.04       0.04     
-##  . backsolve (<text>:6) 0.02       0.02
+##  lr_fast                0.22       0.00     
+##  . crossprod (<text>:2) 0.18       0.18     
+##  . chol (<text>:4)      0.02       0.00     
+##  . . chol.default       0.02       0.02     
+##  . crossprod (<text>:3) 0.02       0.02
 

 
 
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.)

@@ -545,9 +535,9 @@ 
4.1) Pre-allocate memory

 
##      test elapsed replications
-## 1 fun1(z)   2.104           20
+## 1 fun1(z)   2.027           20
 ## 2 fun2(z)   0.026           20
-## 3 fun3(z)   0.006           20
+## 3 fun3(z)   0.005           20

It's not necessary to use as.numeric above though it saves @@ -571,8 +561,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.055 10 -## 1 0.158 10 +## 2 0.053 10 +## 1 0.147 10

For lists, we can do this

@@ -660,7 +650,7 @@

4.2) Vectorized calculations

## } ## .Internal(mean(x)) ## } -## <bytecode: 0x55bf04f01630> +## <bytecode: 0x559af71c2bd0> ## <environment: namespace:base> @@ -673,7 +663,7 @@

4.2) Vectorized calculations

## stop("complex matrices not permitted at present") ## .Internal(La_chol(as.matrix(x), pivot, tol)) ## } -## <bytecode: 0x55bf051f0d18> +## <bytecode: 0x559afb2c2ad0> ## <environment: namespace:base> @@ -771,8 +761,8 @@

4.3) Using apply and specialized functions

##                       test elapsed replications
-## 1 out <- apply(x, 1, mean)   1.880           10
-## 2       out <- rowMeans(x)   0.212           10
+## 1 out <- apply(x, 1, mean)   1.829           10
+## 2       out <- rowMeans(x)   0.209           10

We can 'sweep' out a summary statistic, such as subtracting @@ -782,7 +772,7 @@

4.3) Using apply and specialized functions

##    user  system elapsed 
-##   0.121   0.024   0.144
+##   0.111   0.032   0.143

Here's a trick for doing the sweep based on vectorized calculations, remembering @@ -794,7 +784,7 @@

4.3) Using apply and specialized functions

##    user  system elapsed 
-##   0.304   0.021   0.325
+##   0.281   0.024   0.305
 
identical(out, out2)
@@ -827,7 +817,7 @@ 
Are apply, lapply, sapply, etc. faster than loops
 
##    user  system elapsed 
-##   0.085   0.000   0.084
+##   0.086   0.000   0.086
 
system.time({
@@ -839,7 +829,7 @@ 
Are apply, lapply, sapply, etc. faster than loops
 
##    user  system elapsed 
-##   0.085   0.000   0.084
+##   0.082   0.000   0.082

And here's an example where sapply is much faster, because the core function evaluation at each iteration is very fast:

@@ -861,8 +851,8 @@

Are apply, lapply, sapply, etc. faster than loops 
##      test elapsed replications
-## 1 fun2(z)   2.258           10
-## 2 fun4(z)   0.034           10
+## 1 fun2(z)   2.255           10
+## 2 fun4(z)   0.035           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.

@@ -877,7 +867,7 @@

Are apply, lapply, sapply, etc. faster than loops ## X <- as.list(X) ## .Internal(lapply(X, FUN)) ## } -## <bytecode: 0x55bf0309a0d0> +## <bytecode: 0x559af5368f40> ## <environment: namespace:base> @@ -894,9 +884,9 @@

4.4) Matrix algebra efficiency

##                        test elapsed replications
-## 1        apply(mat, 1, sum)   0.032           10
-## 2 mat %*% rep(1, ncol(mat))   0.003           10
-## 3              rowSums(mat)   0.006           10
+## 1        apply(mat, 1, sum)   0.031           10
+## 2 mat %*% rep(1, ncol(mat))   0.002           10
+## 3              rowSums(mat)   0.007           10

On the other hand, big matrix operations can be slow. Challenge: Suppose you @@ -1036,10 +1026,10 @@

4.5) Fast mapping/lookup tables

## Unit: nanoseconds
 ##         expr   min      lq     mean  median      uq    max neval cld
-##       x[500]   475   523.5   587.95   552.5   605.5   1278   100 a  
-##     x["500"]  6202  6269.5  7407.16  6305.5  6395.0 105028   100  b 
-##    xL[[500]]   185   219.5   394.80   251.0   273.5   8283   100 a  
-##  xL[["500"]] 12043 12108.5 12432.33 12149.5 12213.0  20153   100   c
+##       x[500]  1161  1452.0  1715.65  1631.5  1942.5   2767   100 a  
+##     x["500"] 11117 11800.5 13812.84 12073.5 12433.5 168811   100  b 
+##    xL[[500]]   451   698.5   916.10   833.5  1012.5   7390   100 a  
+##  xL[["500"]] 22140 23256.5 23789.66 23495.0 23731.5  39917   100   c

Lookup by name is slow because R needs to scan through the objects @@ -1079,7 +1069,7 @@

4.5) Fast mapping/lookup tables

## xEnv$"500" 199 222.5 297.29 277.0 298.0 2763 100 a -

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.

@@ -1099,7 +1089,7 @@

Cache-aware programming

##    user  system elapsed 
-##   0.856   0.224   1.081
+##   0.883   0.256   1.139
 
A = t(A)
@@ -1107,7 +1097,7 @@ 
Cache-aware programming

 
##    user  system elapsed 
-##   1.751   0.296   2.047
+##   1.830   0.339   2.169
 
nr = 800

From 7dc8f8ed2fa27a9d48d199e76f95f59d81bd3b37 Mon Sep 17 00:00:00 2001
From: Christopher Paciorek 
Date: Fri, 25 Sep 2020 16:22:15 -0700
Subject: [PATCH 4/6] fix a few things in apply vs loo

---
 efficient-R.R    | 14 ++++++--
 efficient-R.html | 91 ++++++++++++++++++++++++++----------------------
 2 files changed, 60 insertions(+), 45 deletions(-)

diff --git a/efficient-R.R b/efficient-R.R
index 722759d..d110721 100644
--- a/efficient-R.R
+++ b/efficient-R.R
@@ -104,6 +104,7 @@ z <- rnorm(n)
 
 fun1 <- function(vals) {
    x <- exp(vals[1])
+   n <- length(vals) 
    for(i in 2:n) x <- c(x, exp(vals[i]))
    return(x)
 }
@@ -201,16 +202,23 @@ system.time({
 z <- rnorm(10000)
 fun2 <- 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) {
     x <- sapply(vals, exp)
     return(x)
 }
-    
-benchmark(fun2(z), fun4(z),
+
+fun3 <- function(vals) {
+    x <- exp(vals)
+    return(x)
+}
+
+benchmark(fun2(z), fun4(z), fun3(z),
 replications = 10, columns=c('test', 'elapsed', 'replications'))
 
 
diff --git a/efficient-R.html b/efficient-R.html
index 63b6e0d..b0a7183 100644
--- a/efficient-R.html
+++ b/efficient-R.html
@@ -291,7 +291,7 @@ 
3.1) Benchmarking

 
##    user  system elapsed 
-##   0.236   0.012   0.248
+##   0.231   0.020   0.251
 
system.time(rowMeans(x))
@@ -316,10 +316,10 @@ 
3.1) Benchmarking

 
## Unit: nanoseconds
-##           expr   min      lq     mean median      uq   max neval cld
-##       df[2, 1] 13869 14563.0 15817.02  14748 14964.0 93003   100   b
-##     df$vals[2]   997  1625.5  2306.03   2049  2162.5 27105   100  a 
-##  df[2, "vals"] 13792 14334.5 14845.36  14501 14864.0 24041   100   b
+##           expr   min      lq     mean  median      uq   max neval cld
+##       df[2, 1] 13692 14349.5 15099.59 14595.0 14924.5 31647   100   b
+##     df$vals[2]   959  1595.0  2139.88  2072.5  2174.5  9987   100  a 
+##  df[2, "vals"] 13761 14216.0 15452.07 14507.5 14808.0 75373   100   b

microbenchmark is also a good way to compare slower calculations, e.g., doing a crossproduct via crossproduct() compared to “manually”:

@@ -335,8 +335,8 @@

3.1) Benchmarking

## Unit: milliseconds
 ##          expr      min       lq     mean   median       uq      max neval
-##    t(x) %*% x 87.71254 88.20979 89.55519 88.74442 89.52175 97.70361    10
-##  crossprod(x) 43.63842 44.12943 44.76383 44.61579 45.39003 46.04143    10
+##    t(x) %*% x 88.57042 88.62179 89.81577 88.72816 89.26367 96.27504    10
+##  crossprod(x) 44.20568 44.36209 45.53898 44.60348 44.82433 53.05143    10
 ##  cld
 ##    b
 ##   a
@@ -353,8 +353,8 @@ 
3.1) Benchmarking

 
##           test elapsed replications
-## 2 crossprod(x)   0.428           10
-## 1   t(x) %*% x   0.863           10
+## 2 crossprod(x)   0.424           10
+## 1   t(x) %*% x   0.866           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.

@@ -469,10 +469,9 @@

3.2) Profiling

##  path                   total.pct self.pct
 ##  lr_fast                100.00     0.00   
-##  . crossprod (<text>:2)  81.82    81.82   
+##  . crossprod (<text>:2)  90.91    90.91   
 ##  . chol (<text>:4)        9.09     0.00   
-##  . . chol.default         9.09     9.09   
-##  . crossprod (<text>:3)   9.09     9.09
+##  . . chol.default         9.09     9.09
 
hotPaths(pd3, value = 'time')
@@ -480,10 +479,9 @@ 
3.2) Profiling

 
 
##  path                   total.time self.time
 ##  lr_fast                0.22       0.00     
-##  . crossprod (<text>:2) 0.18       0.18     
+##  . crossprod (<text>:2) 0.20       0.20     
 ##  . chol (<text>:4)      0.02       0.00     
-##  . . chol.default       0.02       0.02     
-##  . crossprod (<text>:3) 0.02       0.02
+##  . . chol.default       0.02       0.02
 

 
 
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.)

@@ -517,6 +515,7 @@ 
4.1) Pre-allocate memory

 
 fun1 <- function(vals) {
    x <- exp(vals[1])
+   n <- length(vals) 
    for(i in 2:n) x <- c(x, exp(vals[i]))
    return(x)
 }
@@ -535,9 +534,9 @@ 
4.1) Pre-allocate memory

 
##      test elapsed replications
-## 1 fun1(z)   2.027           20
-## 2 fun2(z)   0.026           20
-## 3 fun3(z)   0.005           20
+## 1 fun1(z)   2.061           20
+## 2 fun2(z)   0.023           20
+## 3 fun3(z)   0.006           20

It's not necessary to use as.numeric above though it saves @@ -561,8 +560,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.053 10 -## 1 0.147 10 +## 2 0.056 10 +## 1 0.150 10

For lists, we can do this

@@ -650,7 +649,7 @@

4.2) Vectorized calculations

## } ## .Internal(mean(x)) ## } -## <bytecode: 0x559af71c2bd0> +## <bytecode: 0x55d05e0e2d70> ## <environment: namespace:base> @@ -663,7 +662,7 @@

4.2) Vectorized calculations

## stop("complex matrices not permitted at present") ## .Internal(La_chol(as.matrix(x), pivot, tol)) ## } -## <bytecode: 0x559afb2c2ad0> +## <bytecode: 0x55d062194bf0> ## <environment: namespace:base> @@ -761,8 +760,8 @@

4.3) Using apply and specialized functions

##                       test elapsed replications
-## 1 out <- apply(x, 1, mean)   1.829           10
-## 2       out <- rowMeans(x)   0.209           10
+## 1 out <- apply(x, 1, mean)   1.849           10
+## 2       out <- rowMeans(x)   0.207           10

We can 'sweep' out a summary statistic, such as subtracting @@ -772,7 +771,7 @@

4.3) Using apply and specialized functions

##    user  system elapsed 
-##   0.111   0.032   0.143
+##   0.141   0.003   0.145

Here's a trick for doing the sweep based on vectorized calculations, remembering @@ -784,7 +783,7 @@

4.3) Using apply and specialized functions

##    user  system elapsed 
-##   0.281   0.024   0.305
+##   0.282   0.029   0.310
 
identical(out, out2)
@@ -817,7 +816,7 @@ 
Are apply, lapply, sapply, etc. faster than loops
 
##    user  system elapsed 
-##   0.086   0.000   0.086
+##   0.084   0.000   0.084
 
system.time({
@@ -829,7 +828,7 @@ 
Are apply, lapply, sapply, etc. faster than loops
 
##    user  system elapsed 
-##   0.082   0.000   0.082
+##   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:

@@ -837,21 +836,29 @@

Are apply, lapply, sapply, etc. faster than loops 
z <- rnorm(10000)
 fun2 <- 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) {
     x <- sapply(vals, exp)
     return(x)
 }
 
-benchmark(fun2(z), fun4(z),
+fun3 <- function(vals) {
+    x <- exp(vals)
+    return(x)
+}
+
+benchmark(fun2(z), fun4(z), fun3(z),
 replications = 10, columns=c('test', 'elapsed', 'replications'))
 
##      test elapsed replications
-## 1 fun2(z)   2.255           10
+## 1 fun2(z)   0.011           10
+## 3 fun3(z)   0.004           10
 ## 2 fun4(z)   0.035           10
@@ -867,7 +874,7 @@

Are apply, lapply, sapply, etc. faster than loops ## X <- as.list(X) ## .Internal(lapply(X, FUN)) ## } -## <bytecode: 0x559af5368f40> +## <bytecode: 0x55d05c27b820> ## <environment: namespace:base> @@ -884,9 +891,9 @@

4.4) Matrix algebra efficiency

##                        test elapsed replications
-## 1        apply(mat, 1, sum)   0.031           10
+## 1        apply(mat, 1, sum)   0.037           10
 ## 2 mat %*% rep(1, ncol(mat))   0.002           10
-## 3              rowSums(mat)   0.007           10
+## 3              rowSums(mat)   0.006           10

On the other hand, big matrix operations can be slow. Challenge: Suppose you @@ -1025,11 +1032,11 @@

4.5) Fast mapping/lookup tables

## Unit: nanoseconds
-##         expr   min      lq     mean  median      uq    max neval cld
-##       x[500]  1161  1452.0  1715.65  1631.5  1942.5   2767   100 a  
-##     x["500"] 11117 11800.5 13812.84 12073.5 12433.5 168811   100  b 
-##    xL[[500]]   451   698.5   916.10   833.5  1012.5   7390   100 a  
-##  xL[["500"]] 22140 23256.5 23789.66 23495.0 23731.5  39917   100   c
+##         expr   min      lq     mean  median      uq   max neval cld
+##       x[500]   439   475.5   522.18   499.0   533.5  1265   100 a  
+##     x["500"]  5702  5772.0  6784.64  5800.0  5843.5 92263   100  b 
+##    xL[[500]]   156   193.0   269.90   213.0   238.5  5488   100 a  
+##  xL[["500"]] 11298 11354.5 11847.55 11390.5 11433.0 19982   100   c

Lookup by name is slow because R needs to scan through the objects @@ -1089,7 +1096,7 @@

4.6 Cache-aware programming

##    user  system elapsed 
-##   0.883   0.256   1.139
+##   0.872   0.248   1.120
 
A = t(A)
@@ -1097,7 +1104,7 @@ 
4.6 Cache-aware programming

 
##    user  system elapsed 
-##   1.830   0.339   2.169
+##   1.688   0.340   2.029
 
nr = 800
@@ -1111,7 +1118,7 @@ 
4.6 Cache-aware programming

 
##                 test elapsed replications
-## 1  apply(A, 2, mean)   0.012           10
+## 1  apply(A, 2, mean)   0.011           10
 ## 2 apply(tA, 1, mean)   0.012           10
From 19744eb4ce711617a1b09f4d00e13e396a2a5a78 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Fri, 25 Sep 2020 16:47:45 -0700 Subject: [PATCH 5/6] updaet for vs apply --- efficient-R.R | 2 +- efficient-R.Rmd | 3 +- efficient-R.html | 107 +++++++++++++++++++++++------------------------ 3 files changed, 56 insertions(+), 56 deletions(-) diff --git a/efficient-R.R b/efficient-R.R index d110721..5caf973 100644 --- a/efficient-R.R +++ b/efficient-R.R @@ -199,7 +199,7 @@ system.time({ ## @knitr apply-vs-for-part2 -z <- rnorm(10000) +z <- rnorm(1e6) fun2 <- function(vals) { x <- as.numeric(NA) n <- length(vals) diff --git a/efficient-R.Rmd b/efficient-R.Rmd index b660e68..91f0f1b 100644 --- a/efficient-R.Rmd +++ b/efficient-R.Rmd @@ -314,7 +314,8 @@ than writing a loop. ```{r, apply-vs-for} ``` -And here's an example where *sapply* is much faster, because the core function evaluation at each iteration is very fast: +I believe that in old version of the R the *sapply* in the next example was faster than 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. ```{r, apply-vs-for-part2} ``` diff --git a/efficient-R.html b/efficient-R.html index b0a7183..aa40900 100644 --- a/efficient-R.html +++ b/efficient-R.html @@ -291,14 +291,14 @@

3.1) Benchmarking

##    user  system elapsed 
-##   0.231   0.020   0.251
+##   0.238   0.016   0.255
 
system.time(rowMeans(x))
 
##    user  system elapsed 
-##   0.023   0.000   0.023
+##   0.024   0.000   0.023

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.

@@ -315,11 +315,11 @@

3.1) Benchmarking

) -
## Unit: nanoseconds
-##           expr   min      lq     mean  median      uq   max neval cld
-##       df[2, 1] 13692 14349.5 15099.59 14595.0 14924.5 31647   100   b
-##     df$vals[2]   959  1595.0  2139.88  2072.5  2174.5  9987   100  a 
-##  df[2, "vals"] 13761 14216.0 15452.07 14507.5 14808.0 75373   100   b
+
## Unit: microseconds
+##           expr    min      lq     mean  median      uq    max neval cld
+##       df[2, 1] 13.982 14.5520 15.08594 14.8910 15.1485 22.980   100   b
+##     df$vals[2]  1.014  1.9675  2.40240  2.1650  2.3220 23.304   100  a 
+##  df[2, "vals"] 13.553 14.3555 15.52685 14.6075 14.8860 62.408   100   b
 

 
 
microbenchmark is also a good way to compare slower calculations, e.g., doing a crossproduct via crossproduct() compared to “manually”:

@@ -335,8 +335,8 @@ 
3.1) Benchmarking

 
 
## Unit: milliseconds
 ##          expr      min       lq     mean   median       uq      max neval
-##    t(x) %*% x 88.57042 88.62179 89.81577 88.72816 89.26367 96.27504    10
-##  crossprod(x) 44.20568 44.36209 45.53898 44.60348 44.82433 53.05143    10
+##    t(x) %*% x 83.37811 84.06533 87.49997 87.84960 89.24589 95.10511    10
+##  crossprod(x) 41.63207 42.61395 43.96466 43.93525 45.33973 45.89729    10
 ##  cld
 ##    b
 ##   a
@@ -354,7 +354,7 @@ 
3.1) Benchmarking

 
 
##           test elapsed replications
 ## 2 crossprod(x)   0.424           10
-## 1   t(x) %*% x   0.866           10
+## 1   t(x) %*% x   0.868           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.

@@ -394,25 +394,25 @@ 
3.2) Profiling

 
 
##  path               total.pct self.pct
 ##  lr_slow            100.00     0.00   
-##  . %*% (<text>:2)    54.05    54.05   
-##  . t (<text>:3)      21.62     0.00   
-##  . . t.default       21.62    21.62   
-##  . solve (<text>:4)  18.92     0.00   
-##  . . solve.default   18.92    18.92   
-##  . t (<text>:2)       5.41     0.00   
-##  . . t.default        5.41     5.41
+##  . %*% (<text>:2)    55.56    55.56   
+##  . t (<text>:3)      22.22     0.00   
+##  . . t.default       22.22    22.22   
+##  . solve (<text>:4)  16.67     0.00   
+##  . . solve.default   16.67    16.67   
+##  . t (<text>:2)       5.56     0.00   
+##  . . t.default        5.56     5.56
 

 
 
hotPaths(pd1, value = 'time')
 

 
 
##  path               total.time self.time
-##  lr_slow            0.74       0.00     
+##  lr_slow            0.72       0.00     
 ##  . %*% (<text>:2)   0.40       0.40     
 ##  . t (<text>:3)     0.16       0.00     
 ##  . . t.default      0.16       0.16     
-##  . solve (<text>:4) 0.14       0.00     
-##  . . solve.default  0.14       0.14     
+##  . solve (<text>:4) 0.12       0.00     
+##  . . solve.default  0.12       0.12     
 ##  . t (<text>:2)     0.04       0.00     
 ##  . . t.default      0.04       0.04
 

@@ -436,21 +436,19 @@ 
3.2) Profiling

 
 
##  path                   total.pct self.pct
 ##  lr_medium              100.00     0.00   
-##  . crossprod (<text>:2)  56.25    56.25   
-##  . solve (<text>:4)      37.50     0.00   
-##  . . solve.default       37.50    37.50   
-##  . crossprod (<text>:3)   6.25     6.25
+##  . crossprod (<text>:2)  58.82    58.82   
+##  . solve (<text>:4)      41.18     0.00   
+##  . . solve.default       41.18    41.18
 

 
 
hotPaths(pd2, value = 'time')
 

 
 
##  path                   total.time self.time
-##  lr_medium              0.32       0.00     
-##  . crossprod (<text>:2) 0.18       0.18     
-##  . solve (<text>:4)     0.12       0.00     
-##  . . solve.default      0.12       0.12     
-##  . crossprod (<text>:3) 0.02       0.02
+##  lr_medium              0.34       0.00     
+##  . crossprod (<text>:2) 0.20       0.20     
+##  . solve (<text>:4)     0.14       0.00     
+##  . . solve.default      0.14       0.14
 

 
 
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.) 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.

@@ -534,8 +532,8 @@ 
4.1) Pre-allocate memory

 

 
 
##      test elapsed replications
-## 1 fun1(z)   2.061           20
-## 2 fun2(z)   0.023           20
+## 1 fun1(z)   2.058           20
+## 2 fun2(z)   0.022           20
 ## 3 fun3(z)   0.006           20
 

 
@@ -649,7 +647,7 @@ 
4.2) Vectorized calculations

 ##     }
 ##     .Internal(mean(x))
 ## }
-## <bytecode: 0x55d05e0e2d70>
+## <bytecode: 0x55ebdd7b91d0>
 ## <environment: namespace:base>
@@ -662,7 +660,7 @@

4.2) Vectorized calculations

## stop("complex matrices not permitted at present") ## .Internal(La_chol(as.matrix(x), pivot, tol)) ## } -## <bytecode: 0x55d062194bf0> +## <bytecode: 0x55ebdf073f20> ## <environment: namespace:base> @@ -760,8 +758,8 @@

4.3) Using apply and specialized functions

##                       test elapsed replications
-## 1 out <- apply(x, 1, mean)   1.849           10
-## 2       out <- rowMeans(x)   0.207           10
+## 1 out <- apply(x, 1, mean)   1.823           10
+## 2       out <- rowMeans(x)   0.214           10

We can 'sweep' out a summary statistic, such as subtracting @@ -771,7 +769,7 @@

4.3) Using apply and specialized functions

##    user  system elapsed 
-##   0.141   0.003   0.145
+##   0.138   0.008   0.146

Here's a trick for doing the sweep based on vectorized calculations, remembering @@ -783,7 +781,7 @@

4.3) Using apply and specialized functions

##    user  system elapsed 
-##   0.282   0.029   0.310
+##   0.278   0.024   0.303
 
identical(out, out2)
@@ -828,12 +826,13 @@ 
Are apply, lapply, sapply, etc. faster than loops
 
##    user  system elapsed 
-##   0.087   0.000   0.087
+##   0.083   0.000   0.082
-

And here's an example where sapply is much faster, because the core function evaluation at each iteration is very fast:

+

I believe that in old version of the R the sapply in the next example was faster than 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)
+
z <- rnorm(1e6)
 fun2 <- function(vals) {
     x <- as.numeric(NA)
     n <- length(vals)
@@ -857,9 +856,9 @@ 
Are apply, lapply, sapply, etc. faster than loops
 

 
 
##      test elapsed replications
-## 1 fun2(z)   0.011           10
-## 3 fun3(z)   0.004           10
-## 2 fun4(z)   0.035           10
+## 1 fun2(z)   1.123           10
+## 3 fun3(z)   0.201           10
+## 2 fun4(z)   6.739           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.

@@ -874,7 +873,7 @@ 
Are apply, lapply, sapply, etc. faster than loops
 ##         X <- as.list(X)
 ##     .Internal(lapply(X, FUN))
 ## }
-## <bytecode: 0x55d05c27b820>
+## <bytecode: 0x55ebdb955820>
 ## <environment: namespace:base>
@@ -891,9 +890,9 @@

4.4) Matrix algebra efficiency

##                        test elapsed replications
-## 1        apply(mat, 1, sum)   0.037           10
+## 1        apply(mat, 1, sum)   0.032           10
 ## 2 mat %*% rep(1, ncol(mat))   0.002           10
-## 3              rowSums(mat)   0.006           10
+## 3              rowSums(mat)   0.005           10

On the other hand, big matrix operations can be slow. Challenge: Suppose you @@ -1032,11 +1031,11 @@

4.5) Fast mapping/lookup tables

## Unit: nanoseconds
-##         expr   min      lq     mean  median      uq   max neval cld
-##       x[500]   439   475.5   522.18   499.0   533.5  1265   100 a  
-##     x["500"]  5702  5772.0  6784.64  5800.0  5843.5 92263   100  b 
-##    xL[[500]]   156   193.0   269.90   213.0   238.5  5488   100 a  
-##  xL[["500"]] 11298 11354.5 11847.55 11390.5 11433.0 19982   100   c
+##         expr   min    lq     mean  median      uq   max neval cld
+##       x[500]   439   476   618.02   510.5   547.0  8385   100 a  
+##     x["500"]  5720  5788  7163.52  5829.5  5881.0 94053   100  b 
+##    xL[[500]]   149   189   298.38   230.0   248.0  6067   100 a  
+##  xL[["500"]] 11649 11727 12529.01 11769.5 11848.5 27445   100   c

Lookup by name is slow because R needs to scan through the objects @@ -1096,7 +1095,7 @@

4.6 Cache-aware programming

##    user  system elapsed 
-##   0.872   0.248   1.120
+##   0.833   0.140   0.974
 
A = t(A)
@@ -1104,7 +1103,7 @@ 
4.6 Cache-aware programming

 
##    user  system elapsed 
-##   1.688   0.340   2.029
+##   1.766   0.228   1.994
 
nr = 800
@@ -1118,7 +1117,7 @@ 
4.6 Cache-aware programming

 
##                 test elapsed replications
-## 1  apply(A, 2, mean)   0.011           10
+## 1  apply(A, 2, mean)   0.012           10
 ## 2 apply(tA, 1, mean)   0.012           10
From a7ce4b803ac60f8621d07d439989c267fde32dc9 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Mon, 4 Oct 2021 16:51:22 -0700 Subject: [PATCH 6/6] updates to improve clarity --- efficient-R.R | 73 +++++----- efficient-R.Rmd | 59 +++++--- efficient-R.html | 350 ++++++++++++++++++++++++++++------------------- 3 files changed, 288 insertions(+), 194 deletions(-) diff --git a/efficient-R.R b/efficient-R.R index 5caf973..3632313 100644 --- a/efficient-R.R +++ b/efficient-R.R @@ -102,23 +102,23 @@ summaryRprof("makeRegr.prof") 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')) ## @knitr init-matrix @@ -200,7 +200,7 @@ system.time({ ## @knitr apply-vs-for-part2 z <- rnorm(1e6) -fun2 <- function(vals) { +fun_loop <- function(vals) { x <- as.numeric(NA) n <- length(vals) length(x) <- n @@ -208,17 +208,17 @@ fun2 <- function(vals) { return(x) } -fun4 <- function(vals) { +fun_sapply <- function(vals) { x <- sapply(vals, exp) return(x) } -fun3 <- function(vals) { +fun_vec <- function(vals) { x <- exp(vals) return(x) } -benchmark(fun2(z), fun4(z), fun3(z), +benchmark(fun_loop(z), fun_sapply(z), fun_vec(z), replications = 10, columns=c('test', 'elapsed', 'replications')) @@ -261,56 +261,65 @@ info <- data.frame( 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 ## @knitr name-lookup -vals <- rnorm(10) -names(vals) <- letters[1:10] -select <- c("h", "h", "a", "c") -vals[select] +info2 <- info$numGrade +names(info2) <- info$grade +info2 +info2[df$clusterLabel] + ## @knitr index-lookup 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 +x[1:3] +xL[1:3] microbenchmark( x[500], # index lookup in vector - x["500"], # name lookup in vector + x["var500"], # name lookup in vector xL[[500]], # index lookup in list - xL[["500"]]) # name lookup in list + xL[["var500"]]) # name lookup in list ## @knitr env-lookup xEnv <- as.environment(xL) # convert from a named list -xEnv$"500" -# I need quotes above because numeric; otherwise xEnv$nameOfObject is fine -xEnv[["500"]] +xEnv$var500 microbenchmark( x[500], xL[[500]], - xEnv[["500"]], - xEnv$"500") - + xEnv[["var500"]], + xEnv$var500 +) ## @knitr cache-aware -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 = t(A) -system.time(apply(A, 1, mean)) ## same calculation, but by row +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')) + +## @knitr cache-aware2 -nr = 800 -nc = 100 +nr <- 800 +nc <- 100 ## small matrix that should fit in cache -A = matrix(rnorm(nr * nc), nrow = nr) -tA = t(A) +A <- matrix(rnorm(nr * nc), nrow = nr) +## Yep, the size is less than the L3 cache: +object.size(A) +memuse::Sys.cachesize() + +tA <- t(A) benchmark(apply(A, 2, mean), ## by column apply(tA, 1, mean), ## by row replications = 10, columns=c('test', 'elapsed', 'replications')) diff --git a/efficient-R.Rmd b/efficient-R.Rmd index 91f0f1b..c90db87 100644 --- a/efficient-R.Rmd +++ b/efficient-R.Rmd @@ -71,7 +71,7 @@ To use an optimized BLAS, talk to your systems adminstrator, see [Section A.3 of 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](https://github.com/berkeley-scf/tutorial-parallel-basics). # 3) Tools for assessing efficiency @@ -98,6 +98,8 @@ package. Of course one would generally only care about accurately timing quick c ```{r, microbenchmark2} ``` +**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. ```{r, benchmark} @@ -129,7 +131,7 @@ Let's run the function with profiling turned on. The first call to *hotPaths* shows the percentage of time spent in each call, while the second shows the actual time. -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 a generic method that immediately calls *solve.default*. +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*. 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. @@ -137,7 +139,7 @@ We can see that a lot of time is spent in doing the two crossproducts. So let's ```{r, Rprof-run2} ``` -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.) 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. +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. ```{r, Rprof-run3} ``` @@ -175,7 +177,7 @@ like this using loops because we would in practice use vectorized calculations. 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. @@ -233,6 +235,7 @@ 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. +Here we see that `nchar` is vectorized and that various arguments to `substring` can be vectors. ```{r, vectorized} ``` @@ -277,8 +280,11 @@ then compare the more efficient code to make sure the results are the same. It's ## 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}. ```{r, apply} @@ -300,9 +306,10 @@ need to transpose twice. ### 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. @@ -314,13 +321,13 @@ than writing a loop. ```{r, apply-vs-for} ``` -I believe that in old version of the R the *sapply* in the next example was faster than in R, but that doesn't seem to be the case currently. +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. ```{r, apply-vs-for-part2} ``` -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. ```{r, lapply-callout} print(lapply) @@ -329,7 +336,9 @@ print(lapply) ## 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. ```{r, matrix-calc} ``` @@ -343,6 +352,8 @@ by another matrix that contains 0s, 1s, and -1s in appropriate places. *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? @@ -351,6 +362,8 @@ be critical for efficiency. How should I order the following calculation? 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 @@ -375,13 +388,13 @@ such as C and Fortran packages. ## 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*. @@ -391,8 +404,11 @@ Here's how you would create an index vector, *grp*, if it doesn't already exist. ```{r, match-lookup} ``` -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: ```{r, name-lookup} ``` @@ -415,7 +431,7 @@ In contrast, to look up by index, R can just go directly to the position of inte 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. +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. ```{r, env-lookup, cache=TRUE} ``` @@ -436,3 +452,8 @@ If you know the size of the cache, you can try to design your code so that in a ```{r, cache-aware} ``` +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. + +```{r, cache-aware2} +``` + diff --git a/efficient-R.html b/efficient-R.html index aa40900..b5ab98b 100644 --- a/efficient-R.html +++ b/efficient-R.html @@ -272,7 +272,7 @@

2) Fast linear algebra

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.

3) Tools for assessing efficiency

@@ -291,14 +291,14 @@

3.1) Benchmarking

##    user  system elapsed 
-##   0.238   0.016   0.255
+##   0.251   0.011   0.262
 
system.time(rowMeans(x))
 
##    user  system elapsed 
-##   0.024   0.000   0.023
+##   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.

@@ -315,11 +315,11 @@

3.1) Benchmarking

) -
## Unit: microseconds
-##           expr    min      lq     mean  median      uq    max neval cld
-##       df[2, 1] 13.982 14.5520 15.08594 14.8910 15.1485 22.980   100   b
-##     df$vals[2]  1.014  1.9675  2.40240  2.1650  2.3220 23.304   100  a 
-##  df[2, "vals"] 13.553 14.3555 15.52685 14.6075 14.8860 62.408   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
 

 
 
microbenchmark is also a good way to compare slower calculations, e.g., doing a crossproduct via crossproduct() compared to “manually”:

@@ -334,14 +334,13 @@ 
3.1) Benchmarking

 
## Unit: milliseconds
-##          expr      min       lq     mean   median       uq      max neval
-##    t(x) %*% x 83.37811 84.06533 87.49997 87.84960 89.24589 95.10511    10
-##  crossprod(x) 41.63207 42.61395 43.96466 43.93525 45.33973 45.89729    10
-##  cld
-##    b
-##   a
+##          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
+

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)
@@ -353,8 +352,8 @@ 
3.1) Benchmarking

 
##           test elapsed replications
-## 2 crossprod(x)   0.424           10
-## 1   t(x) %*% x   0.868           10
+## 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.

@@ -394,32 +393,38 @@

3.2) Profiling

##  path               total.pct self.pct
 ##  lr_slow            100.00     0.00   
-##  . %*% (<text>:2)    55.56    55.56   
-##  . t (<text>:3)      22.22     0.00   
-##  . . t.default       22.22    22.22   
-##  . solve (<text>:4)  16.67     0.00   
-##  . . solve.default   16.67    16.67   
-##  . t (<text>:2)       5.56     0.00   
-##  . . t.default        5.56     5.56
+##  . %*% (<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
 
hotPaths(pd1, value = 'time')
 
##  path               total.time self.time
-##  lr_slow            0.72       0.00     
-##  . %*% (<text>:2)   0.40       0.40     
-##  . t (<text>:3)     0.16       0.00     
-##  . . t.default      0.16       0.16     
-##  . solve (<text>:4) 0.12       0.00     
-##  . . solve.default  0.12       0.12     
+##  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.default      0.04       0.04     
+##  . t (<text>:3)     0.48       0.00     
+##  . . t.default      0.48       0.48

The first call to hotPaths shows the percentage of time spent in each call, while the second shows the actual time.

-

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 a generic method that immediately calls solve.default.

+

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.

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.

@@ -436,22 +441,26 @@

3.2) Profiling

##  path                   total.pct self.pct
 ##  lr_medium              100.00     0.00   
-##  . crossprod (<text>:2)  58.82    58.82   
-##  . solve (<text>:4)      41.18     0.00   
-##  . . solve.default       41.18    41.18
+##  . 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.34       0.00     
-##  . crossprod (<text>:2) 0.20       0.20     
-##  . solve (<text>:4)     0.14       0.00     
-##  . . solve.default      0.14       0.14
+##  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.) 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.

+

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)
@@ -466,20 +475,24 @@ 
3.2) Profiling

 
##  path                   total.pct self.pct
-##  lr_fast                100.00     0.00   
-##  . crossprod (<text>:2)  90.91    90.91   
-##  . chol (<text>:4)        9.09     0.00   
-##  . . chol.default         9.09     9.09
+##  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                0.22       0.00     
-##  . crossprod (<text>:2) 0.20       0.20     
-##  . chol (<text>:4)      0.02       0.00     
-##  . . chol.default       0.02       0.02
+##  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.)

@@ -511,35 +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.058           20
-## 2 fun2(z)   0.022           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.

@@ -558,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.056           10
-## 1   0.150           10
+## 2   0.055           10
+## 1   0.147           10

For lists, we can do this

@@ -647,7 +660,7 @@

4.2) Vectorized calculations

## } ## .Internal(mean(x)) ## } -## <bytecode: 0x55ebdd7b91d0> +## <bytecode: 0x5646e07c69f0> ## <environment: namespace:base> @@ -656,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: 0x55ebdf073f20>
+## <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, ",
@@ -689,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
@@ -746,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)
@@ -758,8 +777,8 @@ 
4.3) Using apply and specialized functions

 

 
 
##                       test elapsed replications
-## 1 out <- apply(x, 1, mean)   1.823           10
-## 2       out <- rowMeans(x)   0.214           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
@@ -769,7 +788,7 @@ 
4.3) Using apply and specialized functions

 

 
 
##    user  system elapsed 
-##   0.138   0.008   0.146
+##   0.138   0.009   0.147
 

 
 
Here's a trick for doing the sweep based on vectorized calculations, remembering
@@ -781,7 +800,7 @@ 
4.3) Using apply and specialized functions

 
##    user  system elapsed 
-##   0.278   0.024   0.303
+##   0.293   0.011   0.303
 
identical(out, out2)
@@ -792,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.

@@ -814,7 +834,7 @@ 
Are apply, lapply, sapply, etc. faster than loops
 
##    user  system elapsed 
-##   0.084   0.000   0.084
+##   0.087   0.000   0.086
 
system.time({
@@ -826,14 +846,14 @@ 
Are apply, lapply, sapply, etc. faster than loops
 
##    user  system elapsed 
-##   0.083   0.000   0.082
+##   0.087   0.000   0.087
-

I believe that in old version of the R the sapply in the next example was faster than in R, but that doesn't seem to be the case currently. +

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(1e6)
-fun2 <- function(vals) {
+fun_loop <- function(vals) {
     x <- as.numeric(NA)
     n <- length(vals)
     length(x) <- n
@@ -841,27 +861,27 @@ 
Are apply, lapply, sapply, etc. faster than loops
     return(x)
 }
 
-fun4 <- function(vals) {
+fun_sapply <- function(vals) {
     x <- sapply(vals, exp)
     return(x)
 }
 
-fun3 <- function(vals) {
+fun_vec <- function(vals) {
     x <- exp(vals)
     return(x)
 }
 
-benchmark(fun2(z), fun4(z), fun3(z),
+benchmark(fun_loop(z), fun_sapply(z), fun_vec(z),
 replications = 10, columns=c('test', 'elapsed', 'replications'))
-
##      test elapsed replications
-## 1 fun2(z)   1.123           10
-## 3 fun3(z)   0.201           10
-## 2 fun4(z)   6.739           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)
 

@@ -873,14 +893,16 @@ 
Are apply, lapply, sapply, etc. faster than loops
 ##         X <- as.list(X)
 ##     .Internal(lapply(X, FUN))
 ## }
-## <bytecode: 0x55ebdb955820>
+## <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),
@@ -890,9 +912,9 @@ 
4.4) Matrix algebra efficiency

 
##                        test elapsed replications
-## 1        apply(mat, 1, sum)   0.032           10
-## 2 mat %*% rep(1, ncol(mat))   0.002           10
-## 3              rowSums(mat)   0.005           10
+## 1        apply(mat, 1, sum)   0.033           10
+## 2 mat %*% rep(1, ncol(mat))   0.003           10
+## 3              rowSums(mat)   0.006           10

On the other hand, big matrix operations can be slow. Challenge: Suppose you @@ -904,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?

@@ -931,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 @@ -962,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.

@@ -983,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 @@ -995,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.

-
vals <- rnorm(10)
-names(vals) <- letters[1:10]
-select <- c("h", "h", "a", "c")
-vals[select]
+
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
 

 
-
##          h          h          a          c 
-##  0.2511293  0.2511293 -0.4466439  0.2898673
+
##  A  B  C 
+## 95 85 75
+

+
+
info2[df$clusterLabel]
+

+
+
##  C  B  B  A  C 
+## 75 85 85 95 75
 

 
 
You can do similar things in terms of looking up by name with dimension
@@ -1020,22 +1055,42 @@ 
4.5) Fast mapping/lookup tables

 
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
-microbenchmark(
+x[1:3]
+

+
+
## var1 var2 var3 
+##    1    2    3
+

+
+
xL[1:3]
+

+
+
## $var1
+## [1] 1
+## 
+## $var2
+## [1] 2
+## 
+## $var3
+## [1] 3
+

+
+
microbenchmark(
     x[500],  # index lookup in vector
-    x["500"], # name lookup in vector
+    x["var500"], # name lookup in vector
     xL[[500]], # index lookup in list
-    xL[["500"]]) # name lookup in list
+    xL[["var500"]]) # name lookup in list
 

 
 
## Unit: nanoseconds
-##         expr   min    lq     mean  median      uq   max neval cld
-##       x[500]   439   476   618.02   510.5   547.0  8385   100 a  
-##     x["500"]  5720  5788  7163.52  5829.5  5881.0 94053   100  b 
-##    xL[[500]]   149   189   298.38   230.0   248.0  6067   100 a  
-##  xL[["500"]] 11649 11727 12529.01 11769.5 11848.5 27445   100   c
+##            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
 

 
 
Lookup by name is slow because R needs to scan through the objects
@@ -1044,17 +1099,10 @@ 
4.5) Fast mapping/lookup tables

 
 
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. 

+
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$"500"  
-

-
-
## [1] 500
-

-
-
# I need quotes above because numeric; otherwise xEnv$nameOfObject is fine
-xEnv[["500"]]
+xEnv$var500  
 

 
 
## [1] 500
@@ -1063,16 +1111,17 @@ 
4.5) Fast mapping/lookup tables

 
microbenchmark(
   x[500],
   xL[[500]],
-  xEnv[["500"]],
-  xEnv$"500")
+  xEnv[["var500"]],
+  xEnv$var500
+)
 

 
 
## Unit: nanoseconds
-##           expr min    lq   mean median    uq   max neval cld
-##         x[500] 466 520.5 687.27  540.5 563.0 13909   100   b
-##      xL[[500]] 168 190.5 252.75  218.5 245.0  3511   100  a 
-##  xEnv[["500"]] 159 170.5 228.14  195.5 219.5  3054   100  a 
-##     xEnv$"500" 199 222.5 297.29  277.0 298.0  2763   100  a
+##              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
 

 
 
4.6 Cache-aware programming

@@ -1087,37 +1136,52 @@ 
4.6 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 
-##   0.833   0.140   0.974
+
##                 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.766   0.228   1.994
+
## 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.012           10
+## 1  apply(A, 2, mean)   0.013           10
 ## 2 apply(tA, 1, mean)   0.012           10