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 

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 40446
##  $ fd : int [1:2] 5 8
##  - attr(*, "class")= chr [1:3] "parallelJob" "childProcess" "process"

str(n)

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

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

## [[1]]
## [1] 0.0002014624

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

## [[1]]
## NULL

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

## NULL

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

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 can iterate across more than one value

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

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.526   0.057   0.608

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

##    user  system elapsed 
##   1.773   0.168   0.415

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

##    user  system elapsed 
##   1.293   0.132   0.479

Bootstrapping is a resampling scheme where the original data is repeatedly reconstructed by taking a sample (with replacement) of the same size of 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$pred_y = predict(l)
d$pred_y_se = predict(l,se=TRUE)$se.fit

ggplot(d, aes(x,y)) +
  geom_point() +
  geom_line(aes(y=pred_y)) +
  geom_line(aes(y=pred_y + 1.96 * pred_y_se), color="red") +
  geom_line(aes(y=pred_y - 1.96 * pred_y_se), color="red")

Re-implement the code below using one of the parallelization techniques we have just discussed, check your performance in creating the bootstrap sample using for 1, 2, and 4 cores. (Work with a neighbor so you are not all running MC code at the same time)

n_rep = 10000
res = matrix(NA, ncol=n_rep, nrow=nrow(d))

for(i in 1:ncol(res))
{ 
  bootstrap_samp = d %>% select(x,y) %>% sample_n(nrow(d), replace=TRUE)
  res[,i] = predict(loess(y ~ x, data=bootstrap_samp), newdata=d)
}

# Calculate the 95% bootstrap prediction interval
d$bs_low = apply(res,1,quantile,probs=c(0.025), na.rm=TRUE)
d$bs_up  = apply(res,1,quantile,probs=c(0.975), na.rm=TRUE)

ggplot(d, aes(x,y)) +
  geom_point() +
  geom_line(aes(y=pred_y)) +
  geom_line(aes(y=pred_y + 1.96 * pred_y_se), color="red") +
  geom_line(aes(y=pred_y - 1.96 * pred_y_se), color="red") +
  geom_line(aes(y=bs_low), color="blue") +
  geom_line(aes(y=bs_up), color="blue")

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)

D - type double, GE - general matrix, MM - matrix / matrix multiplication.

dgemm(   character   TRANSA,
         character   TRANSB,
         integer     M,
         integer     N,
         integer     K,
         double precision    ALPHA,
         double precision, dimension(lda,*)  A,
         integer     LDA,
         double precision, dimension(ldb,*)  B,
         integer     LDB,
         double precision    BETA,
         double precision, dimension(ldc,*)  C,
         integer     LDC 
     )   

DGEMM performs one of the matrix-matrix operations

\[C = \alpha op( A ) \times op( B ) + \beta C\]

where \(op( X )\) is either \(op( X ) = X\) or \(op( X ) = X^T\), \(\alpha\) and \(\beta\) are scalars, and \(A\), \(B\) and \(C\) are matrices, with \(op( A )\) an \(m\) by \(k\) matrix, \(op( B )\) a \(k\) by \(n\) matrix and \(C\) an \(m\) by \(n\) matrix.

library(RhpcBLASctl)
x=matrix(runif(5000^2),ncol=5000)

sizes = c(100,500,1000,2000,3000,4000,5000)
cores = c(1,2,4,8)

sapply(
  cores, 
  function(n_cores) 
  {
    blas_set_num_threads(n_cores)
    sapply(
      sizes, 
      function(s) 
      {
           y = x[1:s,1:s]
           system.time(y %*% y)[3]
      }
    )
  }
)