Better benchmarking (pt. 1)

We can do better (better precision) using the microbenchmark package

install.packages("microbenchmark")

library(microbenchmark)

d = abs(rnorm(1000))
r = microbenchmark(
      exp(log(d)/2),
      d^0.5,
      sqrt(d),
      times = 1000
    )
print(r)

## Unit: microseconds
##           expr    min      lq      mean  median      uq     max neval
##  exp(log(d)/2) 18.885 22.2520 27.099327 22.9090 25.0270 204.077  1000
##          d^0.5 29.842 32.9425 40.330387 34.0590 36.0955 624.320  1000
##        sqrt(d)  2.935  5.8760  7.983377  6.2765  7.1700 315.885  1000

Better benchmarking (pt. 2)

We can also do better using the rbenchmark package

install.packages("rbenchmark")

library(rbenchmark)

d = abs(rnorm(1000))
benchmark(
  exp(log(d)/2),
  d^0.5,
  sqrt(d),
  replications = 1000,
  order = "relative"
)

##            test replications elapsed relative user.self sys.self user.child sys.child
## 3       sqrt(d)         1000   0.008    1.000     0.007    0.001          0         0
## 1 exp(log(d)/2)         1000   0.025    3.125     0.023    0.001          0         0
## 2         d^0.5         1000   0.040    5.000     0.037    0.002          0         0

Example 1

Earlier we mentioned that growing a vector as you collect results is bad, just how bad is it? Benchmark the following three functions and compare their performance.

good = function() {
    res = rep(NA, 1e4)
    for(i in seq_along(res))
        res[i] = sqrt(i)
}

bad = function() {
    res = numeric()
    for(i in 1:1e4)
        res = c(res,sqrt(i))

}

best = function() {
    sqrt(1:1e4)
}

Example 2

Lets compare looping vs. the map (apply) function vs dplyr.

• First we will construct a large data frame

set.seed(523)
d = data_frame(
  A = rnorm(1e5),
  B = runif(1e5),
  C = rexp(1e5),
  D = sample(letters, 1e5, replace=TRUE)               
)

• Implement functions that will return a data frame containing the maximum value of each column

  • A map function
  • A single for loop
  • dplyr

• Benchmark all of your preceding functions using data frame d, which is the fastest, why do you think this is the case? 10 replicates per function is sufficient.

• What would happen if d was smaller? (e.g. take only the first 1000 rows)

mclapply

Parallelized version of lapply

system.time(rnorm(1e6))

##   user  system elapsed 
##  0.101   0.007   0.107 

system.time(unlist(mclapply(1:10, function(x) rnorm(1e5), mc.cores = 2)))

##   user  system elapsed 
##  0.148   0.136   0.106 

system.time(unlist(mclapply(1:10, function(x) rnorm(1e5), mc.cores = 4)))

##   user  system elapsed 
##  0.242   0.061   0.052 ```

system.time(unlist(mclapply(1:10, function(x) rnorm(1e5), mc.cores = 4)))

##   user  system elapsed 
##  0.097   0.047   0.079 

system.time(unlist(mclapply(1:10, function(x) rnorm(1e5), mc.cores = 8)))

##   user  system elapsed 
##  0.193   0.076   0.040 

system.time(unlist(mclapply(1:10, function(x) rnorm(1e5), mc.cores = 10)))

##   user  system elapsed 
##  0.162   0.083   0.041 

system.time(unlist(mclapply(1:10, function(x) rnorm(1e5), mc.cores = 12)))

##   user  system elapsed 
##  0.098   0.065   0.037 

mcparallel

Asynchronously evaluation of an R expression in a separate process

m = mcparallel(rnorm(1e6))
n = mcparallel(rbeta(1e6,1,1))
o = mcparallel(rgamma(1e6,1,1))

str(m)

## List of 2
##  $ pid: int 17004
##  $ fd : int [1:2] 4 7
##  - attr(*, "class")= chr [1:3] "parallelJob" "childProcess" "process"

str(n)

## List of 2
##  $ pid: int 17005
##  $ fd : int [1:2] 5 9
##  - attr(*, "class")= chr [1:3] "parallelJob" "childProcess" "process"

mccollect - waiting

p = mcparallel(mean(rnorm(1e5)))
mccollect(p, wait = FALSE, 10) # will retrieve the result (since it's fast)

## $`17007`
## [1] -0.00102628

mccollect(p, wait = FALSE)     # will signal the job as terminating

## $`17007`
## NULL

mccollect(p, wait = FALSE)     # there is no longer such a job

## NULL

registerDoMC

Primarily used to set the number of cores used by foreach, by default uses options("cores") or half the number of cores found by detectCores from the parallel package.

options("cores")

## $cores
## NULL

detectCores()

## [1] 24

getDoParWorkers()

## [1] 1

registerDoMC(4)
getDoParWorkers()

## [1] 4

foreach

A slightly more powerful version of base for loops (think for with an lapply flavor). Combined with %do% or %dopar% for single or multicore execution.

for(i in 1:10) sqrt(i)

foreach(i = 1:5) %do% sqrt(i)   

## [[1]]
## [1] 1
## 
## [[2]]
## [1] 1.414214
## 
## [[3]]
## [1] 1.732051
## 
## [[4]]
## [1] 2
## 
## [[5]]
## [1] 2.236068

foreach - iterators

foreach can iterate across more than one value, but it doesn’t do length coercion

foreach(i = 1:5, j = 1:5) %do% sqrt(i^2+j^2)   

## [[1]]
## [1] 1.414214
## 
## [[2]]
## [1] 2.828427
## 
## [[3]]
## [1] 4.242641
## 
## [[4]]
## [1] 5.656854
## 
## [[5]]
## [1] 7.071068

foreach(i = 1:5, j = 1:2) %do% sqrt(i^2+j^2)   

## [[1]]
## [1] 1.414214
## 
## [[2]]
## [1] 2.828427

foreach - combining results

foreach(i = 1:5, .combine='c') %do% sqrt(i)   

## [1] 1.000000 1.414214 1.732051 2.000000 2.236068

foreach(i = 1:5, .combine='cbind') %do% sqrt(i)   

##      result.1 result.2 result.3 result.4 result.5
## [1,]        1 1.414214 1.732051        2 2.236068

foreach(i = 1:5, .combine='+') %do% sqrt(i)   

## [1] 8.382332

foreach - parallelization

Swapping out %do% for %dopar% will use the parallel backend.

registerDoMC(4)
system.time(foreach(i = 1:10) %dopar% mean(rnorm(1e6)))

##    user  system elapsed 
##   0.428   0.058   0.604

registerDoMC(8)
system.time(foreach(i = 1:10) %dopar% mean(rnorm(1e6)))

##    user  system elapsed 
##   1.516   0.182   0.646

registerDoMC(12)
system.time(foreach(i = 1:10) %dopar% mean(rnorm(1e6)))

##    user  system elapsed 
##   0.937   0.123   0.360

Example 3 - Bootstraping

Bootstrapping is a resampling scheme where the original data is repeatedly reconstructed by taking a sample (with replacement) of the same size as the original data, and using that to conduct whatever analysis procedure is of interest. Below is an example of fitting a local regression (loess) to some synthetic data, we will construct a bootstrap prediction interval for this model.

set.seed(3212016)
d = data.frame(x = 1:120) %>%
    mutate(y = sin(2*pi*x/120) + runif(length(x),-1,1))

l = loess(y ~ x, data=d)
d = d %>% mutate(
  pred_y = predict(l),
  pred_y_se = predict(l,se=TRUE)$se.fit
) %>% mutate(
  pred_low  = pred_y - 1.96 * pred_y_se,
  pred_high = pred_y + 1.96 * pred_y_se
)

ggplot(d, aes(x,y)) +
  geom_point(color="darkgrey") +
  geom_ribbon(aes(ymin=pred_low, ymax=pred_high), fill="red", alpha=0.25) +
  geom_line(aes(y=pred_y)) +
  theme_bw()

Example 3 - Cont.

We will now re-implement the code below using one of the parallelization techniques we have just discussed and will then check the performance with 1, 2, and 4 cores.

n_rep = 5000
d_xy = select(d, x, y)

res = map(1:n_rep, function(i) {
  d_xy %>% 
    select(x,y) %>% 
    sample_n(nrow(d), replace=TRUE) %>%
    loess(y ~ x, data=.) %>%
    predict(newdata=d) %>%
    setNames(NULL)
}) %>% do.call(cbind, .)

d = d %>% mutate(
  bs_low = apply(res,1,quantile,probs=c(0.025), na.rm=TRUE),
  bs_high  = apply(res,1,quantile,probs=c(0.975), na.rm=TRUE)
)

ggplot(d, aes(x,y)) +
  geom_point(color="gray50") +
  geom_ribbon(aes(ymin=pred_low, ymax=pred_high), fill="red", alpha=0.25) +
  geom_ribbon(aes(ymin=bs_low, ymax=bs_high), fill="blue", alpha=0.25) +
  geom_line(aes(y=pred_y)) +
  theme_bw()

BLAS and LAPACK

Low level algorithms for common linear algebra operations

BLAS

• Basic Linear Algebra Subprograms

• Copying, scaling, multiplying vectors and matrices

• Origins go back to 1979, written in Fortran

LAPACK

• Linear Algebra Package

• Higher level functionality building on BLAS.

• Linear solvers, eigenvalues, and matrix decompositions

• Origins go back to 1992, mostly Fortran (expanded on LINPACK, EISPACK)

OpenBLAS DGEMM (Matrix Multiply) Performance

GPU vs OpenBLAS (NVidia K20X - Kepler)

GPU vs OpenBLAS (NVidia P100 - Pascal)