An example of accelerating calculations in R by multithreading

Introduction

As follows from Wikipedia:

R is a programming language for statistical data processing and graphics, as well as a free open-source computing environment within the GNU project.

This language, at present, has found wide application in many practical and purely scientific fields. However, historically, the speed of many demanding computing leaves much to be desired. The topic of parallel computing in R for habrahabr has already been raised . In this article I will try to show the application of this approach to a specific example using the parallel package for multi-threaded computing.

Formulation of the problem

The need to speed up computations arose when studying the connection of genetic markers with human phenotypic traits using the neuralnet package for building neural networks. However, the proposed solution can be used for other tasks.

Materials and methods

R 2.15.2 was used in the work, neuralnet - 1.32, parallel - 2.15.2, rbenchmark - 1.0.0. Computer running Windows 7, processor Intell Core i5-2550K CPU @ 3.40GHz (4 cores) and 8GB of RAM. The built-in data set (infert, package = "datasets").
')

Results and discussion

First, let's look at the processor load in the usual neural network calculation. The example is taken from the description of the neuralnet package with a slight modification: threshold = 0.0001 is one of the criteria for stopping, and rep = 20–20 repetitions of the computation, in order for the calculations to be more “complex”.

nn<-function() { print(" ") neuralnet(case~parity+induced+spontaneous, infert, err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=20, threshold=0.0001) } nn() 

From the Windows Task Manager, you can see that when executing, only one of the processor cores is used and the CPU load is 25%.
Now measure the execution time of this code. For this, I used the rbenchmark package. In general, the code for measurement will look like this:

 within(benchmark(test.name=test.function(), #      replications=c(3), #     (+1   "") columns=c('test', 'replications', 'elapsed'), #       order=c('elapsed', 'test')), #    { average = elapsed/replications }) #         

The result is:

  test replications elapsed average 1 _ 3 47.83 15.94333333 

Now let's compare this time with performing 20 calculations not through the rep = 20 option, but through, for example, sapply.

 nn.s<-function() { print(" ") nets<-sapply(1:20, function(X) neuralnet(case~parity+induced+spontaneous, infert, err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=1, threshold=0.0001)) } 

It turns out:

  test replications elapsed average 1 _ 3 46.05 15.35 2 _ 3 47.52 15.84 

The execution time is about the same. Apparently there is no optimization of the execution of repetitions inside the neuralnet function. To use the rest of the kernels use the parallel package. It is included in R c version 2.14.0 and allows you to run code in several independent threads. Each thread will perform the neuralnet function.

 nn.p<-function() { print(" ") cl <- makeCluster(getOption("cl.cores", 4)) #       clusterExport(cl,"infert") #     clusterEvalQ(cl,library(neuralnet)) #   neuralnet   parSapply(cl, 1:20, function(X) #   sapply neuralnet(case~parity+induced+spontaneous, infert, err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=1, threshold=0.0001) ) stopCluster(cl) } 

Compare lead time:

  test replications elapsed average 3 _ 3 17.38 5.793333333 2 _ 3 45.88 15.293333333 1 _ 3 46.61 15.536666667 

Thus, the parallel version runs more than 2.5 times faster. In this case, the processor load is 100%.

For convenience, you can create your own wrapper for the neuralnet function.

 pneuralnet <- function(formula, data, rep=1, ..., cl) { clusterExport(cl,"data") clusterEvalQ(cl,library(neuralnet)) nets <- parLapply(cl, 1:rep, function(X) neuralnet(formula, data, rep=1, ...) ) #       reached.threshold nets <- nets[order(sapply(1:rep,function(i){nets[[i]]$result.matrix["reached.threshold", ]}))] return(nets) } cl <- makeCluster(getOption("cl.cores", 4)) nets <- pneuralnet(case~parity+induced+spontaneous, infert, err.fct="ce", linear.output=FALSE, likelihood=TRUE,rep=4, threshold=0.0001, cl=cl) stopCluster(cl)