I'll start this post by a very brief introduction to parallelization using the doMC package in R, in which I introduce a toy example with matrix inversions. I will however not go into detail about how to use the doMC package. If you already are familiar with doMC, you'll probably still want to glance at the matrix inversion functions before you skip ahead to the four points about what I've learned from my experiments, as these functions will occur again later in the text.
Parallel programming in R
Modern processors have several cores, which essentially are independent CPUs that can work simultaneously. When your code is parallelized several computations are carried out in parallel on the different cores. This allows you to, for instance, run several steps of a for-loop in parallel. In R it is simple to do this using e.g. the doMC package.
As an example, imagine that you wish to generate and invert several random 200 by 200 matrices and save the inverted matrices in a list. Using a standard for loop, you can do this using the following function:
Using the doMC package, we can write a parallelized version of this function as follows:# A standard function: myFunc<-function(B,size) { a<-vector("list", B) for(i in 1:B) { myMatrix<-matrix(runif(size^2),size,size) a[[i]]<-solve(myMatrix) } return(a) }
library(doMC) # Use doMC package
registerDoMC() # Determine number of cores to use, default = #cores/2
# A parallelized function:
myParFunc<-function(B,size)
{
foreach(i = 1:B) %dopar%
{
myMatrix<-matrix(runif(size^2),size,size)
solve(myMatrix)
}
}
Comparing the two functions, we find that the parallelized function is substantially faster:
The improvement is even greater for larger matrices. Parallelization can save you a lot of time! Unfortunately, there are situations where it actually makes your code slower, as we will see next.> system.time(myFunc(1000,200)) user system elapsed 19.026 0.148 19.186 > system.time(myParFunc(1000,200)) user system elapsed 24.982 8.448 5.576
1. Parallelization is only useful when you have computation-heavy loop steps
Each parallelization has an overhead cost in that it takes some time to assign the different loop steps to the different cores. If each step of your loops is quickly executed, the overhead cost might actually be larger than the gain from the parallelization! As an example, if you wish to generate and invert several small matrices, say of size 25 by 25, the overhead cost is too large for parallelization to be of any use:
In this case, the parallelized function is much slower than the standard unparallelized version. Parallelization won't magically make all your loops faster!> system.time(myFunc(1000,25)) user system elapsed 0.264 0.004 0.266 > system.time(myParFunc(1000,25)) user system elapsed 0.972 1.480 1.137
2. Only parallelize the outer loops
If you are using nested loops, you should probably refuse the urge to replace all your for-loops by parallelized loops. If you have parallelized loops within parallelized loops, you tend to lose much of the gain from the outer parallelization, as overhead costs increase. This is a special case of point 1 above: it is probably the entire inner loop that makes each step of the outer loop computation-heavy, rather than a single step in the inner loop.
Returning to the matrix inversion example, let's now assume that we wish to use an inner loop to generate a new type of random matrices by multiplying with other matrices with uniformly distributed elements. The unparallelized and parallelized versions of this are as follows:
# No parallelization:
myFunc<-function(B,size)
{
a<-vector("list", B)
for(i in 1:B)
{
myMatrix<-matrix(runif(size^2),size,size)
for(j in 1:B)
{
myMatrix<-myMatrix %*% matrix(runif(size^2),size,size)
}
a[[i]]<-solve(myMatrix)
}
}
# Outer loop parallelized:
myParFunc<-function(B,size)
{
foreach(i = 1:B) %dopar%
{
myMatrix<-matrix(runif(size^2),size,size)
for(j in 1:B)
{
myMatrix<-myMatrix %*% matrix(runif(size^2),size,size)
}
solve(myMatrix)
}
}
# Both loops parallelized:
myParFunc2<-function(B,size)
{
foreach(i = 1:B) %dopar%
{
myMatrix<-matrix(runif(size^2),size,size)
myMatrix<-myMatrix %*% Reduce("%*%",foreach(j = 1:B) %dopar%
{
matrix(runif(size^2),size,size)
})
solve(myMatrix)
}
}> system.time(myFunc(1000,25))
user system elapsed
89.305 0.012 89.381
> system.time(myParFunc(1000,25))
user system elapsed
46.767 0.780 5.239 > system.time(myParFunc2(1000,25))
user system elapsed
970.809 428.371 51.302 3. More is not always merrier when it comes to cores
Using a larger number of cores means spending more time on assigning tasks to different cores. When I tested the code for a large-scale simulation study by using different number of cores on our 64-core number cruncher, it turned out that the code executed much faster when I used only 16 cores instead of 32 or 64. Using 16 cores, I might add, was however definitely faster than using 8 cores. Then again, in another simulation, using 32 cores turned out to be faster than using 16. Less is not always merrier either.
An example where a larger number of cores slows down the execution is the following, which is a modified version of myParFunc2 above, in which the multiplication has been replaced by addition to avoid issues with singular matrices:
myParFunc3<-function(B,size)
{
foreach(i = 1:B) %dopar%
{
myMatrix<-matrix(runif(size^2),size,size)
myMatrix<-myMatrix %*% Reduce("+",foreach(j = 1:B) %dopar%
{
matrix(runif(size^2),size,size)
})
solve(myMatrix)
}
}
> registerDoMC(cores=16)
> system.time(myParFunc3(1000,75))
user system elapsed
2095.571 1340.748 137.397
> registerDoMC(cores=32)
> system.time(myParFunc3(1000,75))
user system elapsed
2767.289 2336.450 104.773
> registerDoMC(cores=64)
> system.time(myParFunc3(1000,75))
user system elapsed
2873.555 3522.984 112.809In this example, using 32 cores if faster than using either 16 or 64 cores. It pays off to make a small trial in order to investigate how many cores you should use for a particular parallelized function!
4. Make fair comparisons to unparallelized code
The doMC package includes the option of running foreach-loops without parallelization by simply typing %do% instead of %dopar%. This is nice as it lets you compare the performance of parallelized and unparallelized loops without having to rewrite the code too much. It is however important to know that the %do% loop behaves differently than the standard for loop that you probably would use if you weren't going to parallelize your code. An example is given by the matrix-inversion functions from point 1 above. The %do% function for this is
But this function is noticeably slower than the standard function described above:# A standard (?) function: myFunc2<-function(B,size) { foreach(i = 1:B) %do% { myMatrix<-matrix(runif(size^2),size,size) solve(myMatrix) } }
> system.time(myFunc(1000,25)) user system elapsed 0.264 0.004 0.266> system.time(myParFunc2(1000,25)) user system elapsed 0.952 0.044 0.994
When comparing parallelized and unparallelized loops, avoid %do% loops, as these make the parallelized loops seem better than they actually are.
I'm honestly not sure why there is such a large discrepancy between the two unparallelized for-loops. I posted a question on this matter on Stack Overflow last year, but to date it has not received any satisfactory answers. If you happen to know wherein the difference lies, I'd much appreciate if you would let me know.
That's all for now. These are my four insights from playing around with parallelization in R. It is entirely possible that they contradict some well-known principles of parallel programming - I confess that I am almost completely ignorant about the literature on parallelization. Please let me know in the comments if you disagree with any of the points above or if you have counterexamples!
For more on parallelization in R, check out R bloggers!
doMc is only one of many many packages for parallelization in R.
