Parallelized loops with R

What is a parallel backend?

When running tasks in parallel, there should be a director node that tells a group of workers what to do with a given set of data and functions. The workers execute the iterations, and the director manages execution and gathers the results provided by the workers. A parallel backend provides the means for the director and workers to communicate, while allocating and managing the required computing resources (processors, RAM memory, and network bandwidth among others).

There are two types of parallel backends that can be used with foreach, FORK and PSOCK.

Setup for a single computer

Setting up a cluster in a single computer requires first to find out how many cores we want to use from the ones we have available. It is recommended to leave one free core for other tasks.

parallel::detectCores()

## [1] 8

n.cores <- parallel::detectCores() - 1

Now we need to define the cluster with parallel::makeCluster() and register it so it can be used by %dopar% with doParallel::registerDoParallel(my.cluster). The type argument of parallel::makeCluster() accepts the strings “PSOCK” and “FORK” to define the type of parallel backend to be used.

#create the cluster
my.cluster <- parallel::makeCluster(
  n.cores, 
  type = "PSOCK"
  )

#check cluster definition (optional)
print(my.cluster)

## socket cluster with 7 nodes on host 'localhost'

#register it to be used by %dopar%
doParallel::registerDoParallel(cl = my.cluster)

#check if it is registered (optional)
foreach::getDoParRegistered()

## [1] TRUE

#how many workers are available? (optional)
foreach::getDoParWorkers()

## [1] 7

Now we can run a set of tasks in parallel!

x <- foreach(
  i = 1:10, 
  .combine = 'c'
) %dopar% {
    sqrt(i)
  }
x

##  [1] 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751 2.828427
##  [9] 3.000000 3.162278

If everything went well, now %dopar% should not be throwing the warning executing %dopar% sequentially: no parallel backend registered, meaning that the parallel execution is working as it should. In this little example there is no gain in execution speed, because the operation being executed is extremely fast, but this will change when the operations running inside of the loop take longer times to run.

Finally, it is always recommendable to stop the cluster when we are done working with it.

parallel::stopCluster(cl = my.cluster)

Setup for a Beowulf cluster

This setup is a bit more complex, because it requires to open a port in every computer of the cluster. Ports are virtual communication channels, and are identified by a number.

First, lets tell R what port we want to use:

#define port
Sys.setenv(R_PARALLEL_PORT = 11000)

#check that it
Sys.getenv("R_PARALLEL_PORT")

Now, we need to open the selected port in every computer of the network. In Linux we need to setup the firewall to allow connections from the network 10.42.1.0/24 (replace this with your network range if different!) to the port 11000 by splitting the window of the Terminator console in as many computers available in your network (the figure below shows three, one for my PC and two for my Intel NUCs), opening an ssh session on each remote machine, and setting Terminator with Grouping equal to Broadcast all so we only need to type the commands once.

Opening port 11000 in three computers at once with Terminator

Now we have to create an object defining the IPs of the computers in the network, the number of cores to use from each computer, the user name, and the identity of the director. This will be the spec argument required by parallel::makeCluster() to create the cluster throughtout the machines in the network. It is a list of lists, with as many lists as nodes are defined. Each sub-list has a slot named host with the IP of the computer where the given node is, and user, with the name of the user in each computer.

The code below shows how this would be done, step by step. Yes, this is CUMBERSOME.

#main parameters
director <- '10.42.0.1'
nuc2 <- '10.42.0.34'
nuc1 <- '10.42.0.104'
user <- "blas"

#list of machines, user names, and cores
spec <- list(
  list(
    host = director, 
    user = user,
    ncore = 7
  ), 
  list(
    host = nuc1, 
    user = user,
    ncore = 4
  ),
  list(
    host = nuc2, 
    user = user,
    ncore = 4
  )
)

#generating nodes from the list of machines
spec <- lapply(
  spec, 
  function(spec.i) rep(
    list(
      list(
        host = spec.i$host, 
        user = spec.i$user)
      ), 
    spec.i$ncore
    )
)

#formating into a list of lists
spec <- unlist(
  spec, 
  recursive = FALSE
)

Generating the spec definition is a bit easier with the function below.

#function to generate cluster specifications from a vector of IPs, a vector with the number of cores to use on each IP, and a user name
cluster_spec <- function(
  ips,
  cores,
  user
){
  
  #creating initial list
  spec <- list()
  
  for(i in 1:length(ips)){
    spec[[i]] <- list()
    spec[[i]]$host <- ips[i]
    spec[[i]]$user <- user
    spec[[i]]$ncore <- cores[i]
  }

  #generating nodes from the list of machines
  spec <- lapply(
    spec, 
    function(spec.i) rep(
      list(
        list(
          host = spec.i$host, 
          user = spec.i$user)
        ), 
      spec.i$ncore
      )
  )

  #formating into a list of lists
  spec <- unlist(
    spec, 
    recursive = FALSE
  )
  
  return(spec)
  
}

This function is also available in this GitHub Gist, so you can load it into your R environment by executing:

source("https://gist.githubusercontent.com/BlasBenito/93ee54d3a98d101754aaff0d658dccca/raw/de57b23740ca90bc02fbd0d5cd3551106ff2fb6d/cluster_spec.R")

Below I use it to generate the input to the spec argument to start the cluster with parallel::makeCluster(). Notice that I have added several arguments.

• The argument outfile determines where the workers write a log. In this case it is set to nowhere with the double quotes, but the path to a text file in the director could be provided here.
• The argument homogeneous = TRUE indicates that all machines have the Rscript in the same location. In this case all three machines have it at “/usr/lib/R/bin/Rscript”. Otherwise, set it up to FALSE.

#generate cluster specification
spec <- cluster_spec(
  ips = c('10.42.0.1', '10.42.0.34', '10.42.0.104'),
  cores = c(7, 4, 4),
  user = "blas"
)

#setting up cluster
my.cluster <- parallel::makeCluster(
  master = '10.42.0.1', 
  spec = spec,
  port = Sys.getenv("R_PARALLEL_PORT"),
  outfile = "",
  homogeneous = TRUE
)

#check cluster definition (optional)
print(my.cluster)

#register cluster
doParallel::registerDoParallel(cl = my.cluster)

#how many workers are available? (optional)
foreach::getDoParWorkers()

Now we can use the cluster to execute a dummy operation in parallel using all machines in the network.

x <- foreach(
  i = 1:20, 
  .combine = 'c'
) %dopar% {
    sqrt(i)
  }
x

Once everything is done, remember to close the cluster.

parallel::stopCluster(cl = my.cluster)

Tuning random forest hyperparameters

Random forest has several hyperparameters that influence model fit:

• num.trees is the total number of trees to fit. The default value is 500.
• mtry is the number of variables selected by chance (from the total pool of variables) as candidates for a tree split. The minimum is 2, and the maximum is the total number of predictors.
• min.node.size is the minimum number of cases that shall go together in the terminal nodes of each tree. For classification models as the ones we are going to fit, 1 is the minimum.

Here we are going to explore how combinations of these values increase or decrease the prediction error of the model (percentage of misclassified cases) on the out-of-bag data (not used to train each decision tree). This operation is usually named grid search for hyperparameter optimization.

To create these combinations of hyperparameters we use expand.grid().

sensitivity.df <- expand.grid(
  num.trees = c(500, 1000, 1500),
  mtry = 2:4,
  min.node.size = c(1, 10, 20)
)

Each row in sensitivity.df corresponds to a combination of parameters to test, so there are 27 models to fit. The code below prepares the cluster, and uses the ability of foreach to work with several iterators at once to easily introduce the right set of hyperparameters to each fitted model.

Notice how in the foreach definition I use the .packages argument to export the ranger package to the environments of the workers.

#create and register cluster
my.cluster <- parallel::makeCluster(n.cores)
doParallel::registerDoParallel(cl = my.cluster)
  
#fitting each rf model with different hyperparameters
prediction.error <- foreach(
  num.trees = sensitivity.df$num.trees,
  mtry = sensitivity.df$mtry,
  min.node.size = sensitivity.df$min.node.size,
  .combine = 'c', 
  .packages = "ranger"
) %dopar% {
  
  #fit model
  m.i <- ranger::ranger(
    data = penguins,
    dependent.variable.name = "species",
    num.trees = num.trees,
    mtry = mtry,
    min.node.size = min.node.size
  )
  
  #returning prediction error as percentage
  return(m.i$prediction.error * 100)
  
}

#adding the prediction error column
sensitivity.df$prediction.error <- prediction.error

To plot the results:

ggplot2::ggplot(data = sensitivity.df) + 
  ggplot2::aes(
    x = mtry,
    y = as.factor(min.node.size),
    fill = prediction.error
  ) + 
  ggplot2::facet_wrap(as.factor(num.trees)) +
  ggplot2::geom_tile() + 
  ggplot2::scale_y_discrete(breaks = c(1, 10, 20)) +
  ggplot2::scale_fill_viridis_c() + 
  ggplot2::ylab("min.node.size")

The figure shows that combinations of lower values of min.node.size and mtry generally lead to models with a lower prediction error across different numbers of trees. Retrieving the first line of sensitivity.df ordered by ascending prediction.error will give us the values of the hyperparameters we need to use to reduce the prediction error as much as possible.

best.hyperparameters <- sensitivity.df %>% 
  dplyr::arrange(prediction.error) %>% 
  dplyr::slice(1)

Confidence intervals of variable importance scores

Random forest has an important stochastic component during model fitting, and as consequence, the same model will return slightly different results in different runs (unless set.seed() or the seed argument of ranger are used). This variability also affects the importance scores of the predictors, and can be use to our advantage to assess whether the importance scores of different variables do really overlap or not.

I have written a little function to transform the vector of importance scores returned by ranger into a data frame (of one row). It helps arranging the importance scores of different runs into a long format, which helps a lot to plot a boxplot with ggplot2 right away. This function could have been just some code thrown inside the foreach loop, but I want to illustrate how foreach automatically transfers functions available in the R environment into the environments of the workers when required, without the intervention of the user. The same will happen with the best.hyperparameters tiny data frame we created in the previous section.

importance_to_df <- function(model){
  x <- as.data.frame(model$variable.importance)
  x$variable <- rownames(x)
  colnames(x)[1] <- "importance"
  rownames(x) <- NULL
  return(x)
}

The code chunk below setups the cluster and runs 1000 random forest models in parallel (using the best hyperparameters computed in the previous section) while using system.time() to assess running time.

#we don't need to create the cluster, it is still up
print(my.cluster)

## socket cluster with 7 nodes on host 'localhost'

#assessing execution time
system.time(
  
  #performing 1000 iterations in parallel
  importance.scores <- foreach(
    i = 1:1000, 
    .combine = 'rbind', 
    .packages = "ranger"
  ) %dopar% {
    
    #fit model
    m.i <- ranger::ranger(
      data = penguins,
      dependent.variable.name = "species",
      importance = "permutation",
      mtry = best.hyperparameters$mtry,
      num.trees = best.hyperparameters$num.trees,
      min.node.size = best.hyperparameters$min.node.size
    )
    
    #format importance
    m.importance.i <- importance_to_df(model = m.i)
    
    #returning output
    return(m.importance.i)
    
  }
  
)

##    user  system elapsed 
##   0.267   0.027   6.556

The output of system.time() goes as follows:

• user: seconds the R session has been using the CPU.
• system: seconds the operating system has been using the CPU.
• elapsed: the total execution time experienced by the user.

This will make sense in a minute. In the meantime, let’s plot our results!

ggplot2::ggplot(data = importance.scores) + 
  ggplot2::aes(
    y = reorder(variable, importance), 
    x = importance
  ) +
  ggplot2::geom_boxplot() + 
  ggplot2::ylab("")

The figure shows that the variable bill_length_mm is the most important in helping the model classifying penguin species, with no overlap with any other variable. In this particular case, since the distributions of the importance scores do not overlap, this analysis isn’t truly helpful, but now you know how to do it!

I assessed the running time with system.time() because ranger() can run in parallel by itself just by setting the num.threads argument to the number of cores available in the machine. This capability cannot be used when executing ranger() inside a parallelized foreach loop though, and it is only useful inside classic for loops.

What option is more efficient then? The code below executes a regular for loop running the function sequentially to evaluate whether it is more efficient to run ranger() in parallel using one core per model, as we did above, or sequentially while using several cores per model on each iteration.

#list to save results
importance.scores.list <- list()

#performing 1000 iterations sequentially
system.time(
  
  for(i in 1:1000){
    
    #fit model
    m.i <- ranger::ranger(
      data = penguins,
      dependent.variable.name = "species",
      importance = "permutation",
      seed = i,
      num.threads = parallel::detectCores() - 1
    )
    
    #format importance
    importance.scores.list[[i]] <- importance_to_df(model = m.i)
    
  }
)

##    user  system elapsed 
##  43.663   2.815  12.948

As you can see, ranger() takes longer to execute in a regular for loop using several cores at once than in a parallel foreach loop using one core at once. That’s a win for the parallelized loop!

We can stop our cluster now, we are done with it.

parallel::stopCluster(cl = my.cluster)