Parallelization in r

There are a few ways to parallelize within R, these generally broken up into 3 different approaches the capture different “for-loop” frameworks in R. They include:

1. parallelization in straight for-loop-style with foeach
2. parallelization in apply-style loops with mclapply
3. parallelization with tidyverse nesting-style loops with multidplyr (not discussed see: R and python parallelization blog for discussion)

The final option the above list (as of 8/2021) is still too experimental and limited to truly encourage.

General starting point

To start parallelizing within R we need to have some understanding of how many cores we have. To do figures with out we use the following lines of code:

 

library(parallel)
numCores <- detectCores()
numCores

This will tell you the total number of cores on your computer. We often recommend using at least we less than the total number of cores when working on your personal computer or on the department servers - to leave room for other jobs / use.

Extra resources:

• Matt Jones
• R and python more options

for-loop style: foreach and doParallel library (foreach () %dopar%)

The combination of foreach and doParallel can convert for-loop style coding into parallel processing and provide a lot of additional parameters that can help reduce unnecessary memory usage, provide progress reports and more.

In the following presentation, we’ll demonstrate the speed obtained with different approaches using a simulation function defined below

 

my_k <- rep(1:14, each = 20)

my_simulation <- function(k) {
  x <- runif(30,-3, 3)
  y <- 2*x + 2* rnorm(n = 30)
  x0 <- runif(10000,-3, 3)
  y0 <- 2*x0 + 2*rnorm(n = 10000)

  my_lm <- lm(y ~ poly(x,k))
  y0_hat <- predict(my_lm, data.frame(x = x0))
  return(mean((y0-y0_hat)^2))
}

single core analysis

We naturally want to compare a parallelization approach with the base single core analysis, which can be done with the code as follows:

 

#base for-loop
system.time({
  storage_f <- c()
  for (k in my_k) {
    storage_f <- c(storage_f, my_simulation(k))
  }
})

##    user  system elapsed
##   1.811   0.633   2.449

Additionally, in the foreach package we can also use just a single core in the following way:

 

#single-core `foreach`
library(foreach)
system.time({
  storage_fe <- foreach (k=my_k, .combine=c) %do% {
    my_simulation(k)
  }
})

##    user  system elapsed
##   1.876   0.640   2.522

This structure in the above code using foreach and %do% are very similar to the actual code we use, so it’s useful to understand. It should be clear that the foreach function works very similar to the base for-loop but we need to combined it with %do%. Additionally, notice that the foreach loops of k for values in my_k and we express that differently than in a for loop (k in my_k vs k = my_k).

multi-core analysis

To multiple cores with foreach we just change the above code to

 

library(doParallel)
registerDoParallel(numCores-2)  # use multicore, set to the number of our cores

system.time({
  storage_p <- foreach (k=my_k, .combine=c) %dopar% {
     my_simulation(k)
  }
})
stopImplicitCluster()

##    user  system elapsed
##   3.513   2.520   0.601

The last number under “elapsed” suggest that although we used more core + user time overall we completed the overall process faster than without parallelization.

additional parameters

1. .combine: this parameter tells the system how you’d like to return of each step of the for-loop combined at the end. Above we used .combine = c to get a vector back, but we could also use things like rbind, cbind. Note that the default is a list.
2. progress bars: To have a running progress bar one can use the package doSNOW. For the previous process we can add a progress bar with the following set of code

   library(doSNOW)
   cl <- makeCluster(numCores-4)
   registerDoSNOW(cl)
   pb <- txtProgressBar(max = length(my_k), style = 3)
   progress <- function(n) setTxtProgressBar(pb, n)
   opts <- list(progress = progress)
   
   storage_p <- foreach (k=my_k, .combine=c,
                         .options.snow = opts) %dopar% {
      my_simulation(k)
   }
   parallel::stopCluster(cl)

|==============================================================        | 83%

1. memory and library passing: Sometimes you have a lot of functions or memory in your RAM before your run for parallelized for-loop. The foreach function tries to reduce the amount memory used for each parallel process by constraining which functions are passed in. To specify which packages you’d like to use use the .packages command in foreach, for example: .packages = c("dplyr", "tidyr"). Additionally, you can constraint what data is passed into each process with the parameter .export (which also uses strings relative to your global envirnoment).

 

additional resources

doParallel + foreach vignette

lapply style: parallel library (mclapply)

The second package we will look at is parallel, which includes the function mclapply which can be used very similar to an lapply.

single core analysis

For a single-core analysis of the same analysis seen in the previous section, we could run

 

system.time({
  results <- lapply(my_k, my_simulation)
  my_data <- data.frame(k = my_k,
                        values = unlist(results))
})

##    user  system elapsed
##   1.819   0.641   2.467

multi core analysis

Now, with the mclapply function we can do basically the same thing:

 

system.time({
  results <- mclapply(my_k, my_simulation, mc.cores = numCores)
  my_data <- data.frame(k = my_k,
                        values = unlist(results))
})

##    user  system elapsed
##   2.208   1.636   0.522

additional resources

• parallel manual

The fastest one as of February 2024: BiocParallel

All of the above are good solutions, but from practical experience we have seen that sometimes they can incur in high overheads. A consistently faster solution comes from a package developed and maintainted by the (huge) bioinformatics community that works with R:  BioConductor. The package is named BiocParallel and its interface is almost exactly the same as the one for the Parallel package seen in the previous section. The only difference is that we need another package.

First of all, install the bioconductor package manager via

install.packages("BiocManager")

This will make available an install command, which we use as follows:

install("BiocParallel")

We can then use a simple command. Assuming the function to run multiple times in parallel is named foo, we type

library('BiocParallel')
system.time( 
     y <- bplapply(seq_len(100), foo) 
)

See this article for more information.