Have more things happen at once: Parallel Programming

Parallel processing is about using multiple cores of your computer’s CPU to run multiple tasks simultaneously. This enables you to complete the same task multiple times quicker!

In R, usually computations run sequentially. When we initiate multiple tasks they are performed one after the other, new task starts only after the previous one is completed.

This might become the bottleneck when you come across a computationally heavy process. Such as running monte carlo simulations or fitting multiple machine learning models.Typically, if I have a process which runs longer than 3 minutes, I consider using parallel processing.

It might sound complicated but parallel computing is simple. When you have a repetitive task which takes too much of your valuable time why not implement it and save time? Even if you have a single task, you can benefit from parallel processing by dividing the task into smaller pieces.

Two widely used packages for parallel processing are parallel and foreach.

Let’s do some parallel processing..

What is our goal?

I recently analyzed data from horror movies released in the world between 2012 and 2017. What caught my attention was that movies are more likely to be released at the 13th of each month. You can visit my recent post if you want to know more about this data.

library(tidyverse) # ggplot2, dplyr, tidyr, readr, 
                   # purrr, tibble, stringr, forcats
horror_movies <- readr::read_csv("https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2019/2019-10-22/horror_movies.csv")

dim(horror_movies)
## [1] 3328   12
horror_date <-  horror_movies %>% 
                separate(
                  release_date, 
                  c("day", "month", "year"),
                  sep = "-")

horror_date$day <- as.numeric(horror_date$day)

# Remove rows without Date of the month

horror_date <- horror_date %>% filter(day < 32) 

# I am excluding Day 1 from the analysis (Most aggreements starts at 1st)

horror_date_table <- horror_date %>% filter(day > 1)

# let's check what is the most common day in the month for a horror movie release
data <-  horror_date_table %>%
  group_by(day) %>% 
  count() %>% 
  arrange(desc(n))

data %>% print(n=30)
## # A tibble: 30 x 2
## # Groups:   day [30]
##      day     n
##    <dbl> <int>
##  1    13   124
##  2    18   119
##  3    25   119
##  4    21   110
##  5    31   107
##  6    28   102
##  7    10   100
##  8    20   100
##  9     5    99
## 10     7    98
## 11     2    95
## 12    29    94
## 13    14    91
## 14    12    90
## 15    24    90
## 16     8    88
## 17    15    88
## 18    17    88
## 19    16    87
## 20     3    86
## 21     4    86
## 22     6    86
## 23     9    85
## 24    22    85
## 25    26    84
## 26    11    83
## 27    27    82
## 28    19    77
## 29    30    74
## 30    23    65

multfun <- function(x){ 
  a <- list()
  avc <- rmultinom(1000000, size=2782, prob=rep(1/30, 30))  # Create probability matrix
  for (i in 1:30){ a[[i]] <- mean(avc[i,1:1000000] >=x) }   # Find probabilities for each day
  ifelse(do.call(sum,a)/30 <= 0.5, p_val <-  2 * do.call(sum,a)/30, 
                                   p_val <- 2*(1 - do.call(sum,a)/30))
}

# Run an example
multfun(124)
## [1] 0.001844267

system.time({
multfun(124)
multfun(119)
multfun(110)
multfun(107)
multfun(102)
})
##    user  system elapsed 
##   15.72    0.21   15.94

system.time({
lapply(c(124,119,119,197,102), multfun)
})
##    user  system elapsed 
##   15.64    0.29   15.92

Using parallel package functions: parLapply()

If you have ever used any of lapply() or sapply function you are halfway in parallel processing.

library(parallel)
system.time({
n_cores <- detectCores(logical=FALSE)
cl <- makeCluster(n_cores-1)
parLapply(cl , c(124,119,119,197,102), multfun)
})
##    user  system elapsed 
##    0.02    0.01    4.54
stopCluster(cl)

multfun <- function(x){ 
  a <- list()
  avc <- rmultinom(1000000, size=2782, prob=rep(1/30, 30))  # Create probability matrix
  for (i in 1:30){ a[[i]] <- mean(avc[i,1:1000000] >=x) }   # Find probabilities for each day
  ifelse(do.call(sum,a)/30 <= 0.5, p_val <-  2 * do.call(sum,a)/30, 
                                   p_val <- 2*(1 - do.call(sum,a)/30))
                       }

Improved multfun: Define the function without using a loop and list

Mean function in R can be vectorized instead of looping.

I can achieve the same result what the first multfun() was doing by including what for loop was doing inside the “[ ]”. Mean function will be applied find the total proportion of values equal or bigger than x.

multfun_opt <- function(x){ 
  avc <- rmultinom(1000000, size=2782, rep(1/30, 30))
  p <- mean(avc[1:30, 1:1000000] >=x)
  ifelse(p <= 0.5, p_val <-  2 * p, p_val <- 2*(1 - p))}

Let’s check whether improved multfun is faster.

system.time({
n_cores <- detectCores(logical=FALSE)
cl <- makeCluster(n_cores-1)
parLapply(cl , c(124,119,119,197,102), multfun_opt)
})
##    user  system elapsed 
##    0.00    0.06    4.32
stopCluster(cl)

Since rmultinom is the bottleneck here, we did not see very big improvement in overall function. We got 10-20% improvement. But it is good to know, what we can improve when we really need it.

Let’s compare only the part of the function which uses the vectorized operation and the one with for loop.

avc <- rmultinom(1000000, size=2782, prob=rep(1/30, 30)) 
a <- list()

system.time({
  for (i in 1:30){ a[[i]] <- mean(avc[i,1:1000000] >=124) }  
  ifelse(do.call(sum,a)/30 <= 0.5, p_val <-  2 * do.call(sum,a)/30, 
                                   p_val <- 2*(1 - do.call(sum,a)/30)) })
##    user  system elapsed 
##    0.41    0.01    0.42
system.time({
    p <- mean(avc[1:30, 1:1000000] >=124)
  ifelse(p <= 0.5, p_val <-  2 * p, p_val <- 2*(1 - p)) })
##    user  system elapsed 
##    0.22    0.02    0.24

Vectorised version is twice faster!!

Further Improvement

You can use Rprof() function to profile your process, it will return how much time each subprocess takes time. Let’s profile our function.

multfun_opt <- function(x){ 
  avc <- rmultinom(1000000, size=2782, rep(1/30, 30))
  p <- mean(avc[1:30, 1:1000000] >=x)
  ifelse(p <= 0.5, p_val <-  2 * p, p_val <- 2*(1 - p))
  }

Rprof()
multfun_opt(124)
## [1] 0.0018542
Rprof(NULL)

summaryRprof()$by.self
##                self.time self.pct total.time total.pct
## "rmultinom"         1.80    93.75       1.80     93.75
## ">="                0.08     4.17       0.08      4.17
## "mean"              0.02     1.04       0.12      6.25
## "mean.default"      0.02     1.04       0.02      1.04

In profile summary, total.pct shows the total percentage of the duration of each task inside the function. We can see that rmultinom() part of takes more than 90% of the processing time.

If we are looking further improvement we can focus on it. To improve rmultinom() part, I will isolate it by splitting the multfun into two functions.

Above, first improvement was using parLapply instead of lapply and second replacing a loop with a faster vectorized operation. Now, we will utilize another approach, dividing a task into smaller pieces.

I need to create a new rmultinom() function, which I call fastbinom. To keep things simple, it will take only one argument which I will use to divide 1.000.000 rmultinom simulations into 5 equal parts of 200.000 rmultinom simulations. Each will be handled by a different core in my CPU boosting my speed.

fastbinom <- function(x) {
  rmultinom(x, size=2782, rep(1/30,30)) 
}

multfun_fin <- function(x){ 
  p <- mean(avc[1:30, 1:1000000] >=x)
  ifelse(p <= 0.5, p_val <-  2 * p, p_val <- 2*(1 - p))
}

Here comes an additional overhead. Since, multfun_fin will be executed on different cores, we need to export the data to each cluster. Windows rely on clusterExport() for data to be copied to each cluster.

Test the final optimized code.

# Final optimized code 
system.time({
cl <- makeCluster(5)
avc <- parLapply(cl , c(200000, 200000, 200000, 200000, 200000), fastbinom)
avc <- do.call(cbind, avc)
clusterExport(cl, varlist="avc", envir=environment())
parLapply(cl, c(124,119,119,197,102), multfun_fin)
})
##    user  system elapsed 
##    0.60    0.61    3.05
stopCluster(cl)

Quite a bit improvement!

Final test

Let’s now put things together, I will test different versions I created so far. But this time for more input values, I will use my function for each of the 30 categories and by using up to 1.000.000 simulations.

Let’s call the initial version: original, second parLapply, third parLapply, optimized and the last parLapply/fastbinom

multfun <- function(x){ 
  a <- list()
  avc <- rmultinom(1000000, size=2782, prob=rep(1/30, 30))  # Create probability matrix
  for (i in 1:30){ a[[i]] <- mean(avc[i,1:1000000] >=x) }   # Find probabilities for each day 
  ifelse(do.call(sum,a)/30 <= 0.5, p_val <-  2 * do.call(sum,a)/30, 
                                   p_val <- 2 * (1 - do.call(sum,a)/30))
}


# Original
system.time({
lapply(data$n, multfun)
})
##    user  system elapsed 
##   92.37    1.43   93.82
# parLapply
system.time({
n_cores <- detectCores(logical=FALSE)
cl <- makeCluster(n_cores-1)
parLapply(cl , data$n, multfun)
})
##    user  system elapsed 
##    0.03    0.03   24.97
stopCluster(cl)

multfun_opt <- function(x){ 
  avc <- rmultinom(1000000, size=2782, rep(1/30, 30))
  p <- mean(avc[1:30, 1:1000000] >=x)
  ifelse(p <= 0.5, p_val <-  2 * p, p_val <- 2*(1 - p))
  }

# parLapply, optimized
system.time({
n_cores <- detectCores(logical=FALSE)
cl <- makeCluster(n_cores-1)
parLapply(cl , data$n, multfun_opt)
})
##    user  system elapsed 
##    0.00    0.03   21.75
stopCluster(cl)

multfun_fin <- function(x){ 
  p <- mean(avc[1:30, 1:1000000] >=x)
  ifelse(p <= 0.5, p_val <-  2 * p, p_val <- 2*(1 - p))
}

# parLapply, fastbinom
system.time({
cl <- makeCluster(5)
avc <- parLapply(cl , c(200000, 200000, 200000, 200000, 200000), fastbinom)
avc <- do.call(cbind,avc)
clusterExport(cl, varlist="avc", envir=environment())
parLapply(cl, data$n, multfun_fin)
})
##    user  system elapsed 
##    0.56    0.60    4.10
stopCluster(cl)

Here is an additonal comparision, testing time differences between different versions and to be more precise I simulated up to 10.000.000 multinomial distributions.

library(readxl)
library(tidyverse)
library(viridis)
## Loading required package: viridisLite
speedtest <- read_excel("posts_data/parallelprocessing_speed_tests.xlsx", sheet =3)
## New names:
## * `` -> ...1
speedtest <- speedtest[-1]
spd <- gather(speedtest, version, seconds, -c(rmultinom,n))
spd$version <- factor(spd$version, labels = c("original", "parlapply", "parLapply/opt", "parLapply/fastbinom"))

spdplot <- spd %>% filter(n >5) %>% 
  ggplot(aes(x=rmultinom, y=seconds, color=version)) + 
  geom_line() + 
  theme(text = element_text(size=18)) +
   labs(  x = "Number of simulations", 
          y = "Time (seconds)")
spdplot

spd %>% filter(rmultinom ==10000000, n >5) %>%
  ggplot(aes(x=version, y=seconds)) + geom_col() + 
  theme(text = element_text(size=15)) + 
  geom_label(aes(label = seconds), alpha = 0.5)

By splitting our original function into two and then applying parallel processing for each subtask we achieved more than 30X speed gain.