Parallelization in r

There are a few ways to parallelize within R, 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

The final option in the above list 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.

library(parallel)
numCores <- detectCores()
numCores

This will tell you the total number of cores on your computer. We often recommend using less than the total number of cores when working on your personal computer or on the department servers to leave room for other jobs. You can always determine the current server system use by checking the hydra usage page.

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 using a standard for loop over various values of k :

#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

In the foreach package, we can also use just a single core in the following way which results in approximately the same run time as the previous code segment:

#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

It should be clear that the foreach function works very similar to the base for-loop; the difference in the syntax is mainly the fact that there is a %do% before the opening bracket. Additionally, notice that the loops of k for values in my_k is expressed differently than in a standard for loop (k in my_k vs k = my_k).

Multi-Core Analysis

To use 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. The %dopar% represents "do parallelization."

Additional Parameters

1. .combine: this parameter tells the system how you’d like to return each step of the for-loop when 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. Memory and library passing: Sometimes you have a lot of functions or memory in your RAM before your run the 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 constrain what data is passed into each process with the parameter .export which uses strings relative to your global environment.
3. 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%

 

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 simulation 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 as before. The difference here is the same as the difference between using for and lapply in base R. mclapply applies each of the values in my_k to the simulation, and parallelizes it across numCores cores.

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

Extra Resources:

• Matt Jones
• R and python more options
• parallel manual
• doParallel + foreach vignette

The Fastest One as of February 2024: BiocParallel

All of the above are good solutions, but from practical experience, we have seen that sometimes the previous methods can incur high overheads. A consistently faster solution comes from a package developed and maintained by the R bioinformatics community:  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 BiocParallel 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.