How to optimize my R script in order to use “multicore”

I am using GNU R at a Ubuntu-Lucid PC which has 4 CPUs. In order to use all 4 CPUs, I installed the "r-cran-multicore" package. As the package's manual lacks of practical examples that I understand, I need advice in how to optimize my script in order to make use of all 4 CPUs.

My dataset is a data.frame (called P1) that has 50,000 rows and 1600 cols. For each row, I'd like to calc the maximun, sum and mean. My script looks as follows:

p1max <- 0
p1mean <- 0
p1sum <-0
plength <- length(P1[,1])
for(i in 1:plength){
p1max <- c(p1max, max(P1[i,]))
p1mean <- c(p1mean, mean(P1[i,]))
p1sum <- c(p1sum, sum(P1[i,]))
}


Could anyone please tell me how to modify and run the script in order to use all 4 CPUs?

-
there's an error in the above program: the line should be "for(i in 1:plength)" –  Simon Byrne Jan 19 '11 at 13:38
you are rigth, thx! –  Produnis Jan 19 '11 at 13:40
doesn't this belong on StackOverflow? –  R_Coholic Jan 20 '11 at 2:14
This does belong on StackOverflow. There's no statistics question here at all. Only a general programming question. –  JD Long Jan 24 '11 at 8:31

Use foreach and doMC. The detailed explanation can be found here. Your script will change very little, the line

for(i in 1:plength){


should be changed to

foreach(i=1:plength) %dopar% {


The prerequisites for any multitasking script using these packages are

library(foreach)
library(doMC)
registerDoMC()


Note of caution. According to the documentation you cannot use this in GUI.

As for your problem, do you really need multitasking? Your data.frame takes about 1.2GB of RAM, so it should fit into your memory. So you can simply use apply:

p1smry <- apply(P1,1,summary)


The result will be a matrix with summaries of each row.

You can also use function mclapply which is in the package multicore. Then your script might look like this:

loopfun <- function(i) {
summary(P1[i,])
}

res <- mclapply(1:nrow(P1),loopfun)


This will return the list, where i-th element will be the summary of i-th row. You can convert it to matrix using sapply

mres <- sapply(res,function(x)x)

-
thank you very much. You are right, that with "apply" the script could be optimized. I just used my script as a minimal-example in order to get the message through... Thx a lot, your answer is exactly what I was looking for!! –  Produnis Jan 19 '11 at 14:27

You've already got an answer as to how to use more than one core, but the real problem is with the way you have written your loops. Never extend your result vector/object at each iteration of a loop. If you do this, you force R to copy your result vector/object and extend it which all takes time. Instead, preallocate enough storage space before you start the loop and fill in as you go along. Here is an example:

set.seed(1)
p1 <- matrix(rnorm(10000), ncol=100)
system.time({
p1max <- p1mean <- p1sum <- numeric(length = 100)
for(i in seq_along(p1max)){
p1max[i] <- max(p1[i,])
p1mean[i] <- mean(p1[i,])
p1sum[i ]<- sum(p1[i,])
}
})

user  system elapsed
0.005   0.000   0.005


Or you can do these things via apply():

system.time({
p1max2 <- apply(p1, 1, max)
p1mean2 <- apply(p1, 1, mean)
p1sum2 <- apply(p1, 1, sum)
})
user  system elapsed
0.007   0.000   0.006


But note that this is no faster than doing the loop properly and sometimes slower.

However, always be on the lookout for vectorised code. You can do row sums and means using rowSums() and rowMeans() which are quicker than either the loop or apply versions:

system.time({
p1max3 <- apply(p1, 1, max)
p1mean3 <- rowMeans(p1)
p1sum3 <- rowSums(p1)
})

user  system elapsed
0.001   0.000   0.002


If I were a betting man, I would have money on the third approach I mention beating foreach() or the other multi-core options in a speed test on your matrix because they would have to speed things up considerably to justify the overhead incurred in setting up the separate processes that are farmed out the the different CPU cores.

Update: Following the comment from @shabbychef is it faster to do the sums once and reuse in the computation of the mean?

system.time({
p1max4 <- apply(p1, 1, max)
p1sum4 <- rowSums(p1)
p1mean4 <- p1sum4 / ncol(p1)
})

user  system elapsed
0.002   0.000   0.002


Not in this test run, but this is far from exhaustive...

-
thanks Gavin, this is really helpfull!!! –  Produnis Jan 19 '11 at 16:44
FWIW, Matlab has the same issues regarding preallocation and expanding vectors, and is a classic code 'blooper'. In addition to your wager, it is probably faster to use the results of rowSums to compute the row means (unless I am missing something regarding e.g. Na or NaN). The code in your third approach sums each column twice. –  shabbychef Jan 19 '11 at 20:16
@shabbychef you'll be surprised (see my edited answer). Yes the sums are notionally computed twice, but rowSums and rowMeans are highly optimised compiled code and what we gain in only computing the sums once, we loose again in doing the mean computation in interpreted code. –  Gavin Simpson Jan 19 '11 at 21:34
@Gavin Simpson: not so fast: try instead system.time({ for (iii in c(1:1000)) { p1max3 <- apply(p1, 1, max) p1mean3 <- rowMeans(p1) p1sum3 <- rowSums(p1) } }) and similarly system.time({ for (iii in c(1:1000)) { p1max4 <- apply(p1, 1, max) p1sum4 <- rowSums(p1) p1mean4 <- p1sum4 / ncol(p1) } }); the version which does not recompute the sum takes 1.368 seconds on my computer; the one which does takes 1.396. again, far from exhaustive, but more compelling... –  shabbychef Jan 20 '11 at 0:45
@shabbychef we must have different ideas on what is or is not compelling ;-) In fact, your more rigorous simulations reinforce my main point, that as rowMeans and rowSums are implemented in efficient, optimised compiled code they are going to be difficult to beat. –  Gavin Simpson Jan 20 '11 at 8:14

Have a look at the snow and snowfall packages. Plenty of examples with those...

If you want to speed up that specific code rather than learning about R and parallelism you should do this

P1 = matrix(rnorm(1000), ncol=10, nrow=10
apply(P1, 1, max)
apply(P1, 1, mean)
apply(P1, 1, sum)

-