Parallel Design-Matrix, Parameter-Vector multiplication in ML-Estimation

Question

I do maximum-likelihood-optimization with R's optim-procedure (using BFGS).

Each time the target function is evaluated a couple of vector-matrix-multiplications will take place where the same design-matrix will be postmultilpied with the vector of changing parameters.

If I simply use the parallel package to distribute the job among cores, the distribution time within each iteration basically kills the computation time such it even takes longer with the parallel version if compared to the ordinary matrix vector product.

What I would like to do is to distribute the pieces of the matrix among cores once and then perform the multiplication on the pieces since the matrix will not change between iterations. Basically I do not want that the same object will be distributed every iteration.

What Ive done so far is

nc <- detectCores()
cl <- makeCluster(rep("localhost", nc))
matprod.par <- function(cl, A, B){
    idx <- splitIndices(nrow(A), length(cl))
    Alist <- lapply(idx, function(ii) A[ii,,drop=FALSE])
    ans <- clusterApply(cl, Alist, get("%*%"), B)
    do.call(rbind, ans)
}

Here, the clusterApply-Function distributes the pieces of A, i.e., Alistamong cores. Is there a possibility to distribute Alist among cores once and then perform the multiplication on the distributed pieces and put them back together via clusterApply?

/edit

I compared the clusterCall-approach of Steve Weston with a simple foreach-%dopar%-approch:

matprod.par1 <-function(Alist,B,nc){
    par <- foreach(i=1:nc, .combine=rbind) %dopar%{
        Alist[[i]]%*%B
    }
}

and the clusterExport-approach which, unlike the clusterCall-approach, copies everything to every spawn.

matprod.par2 <-function(obj1,obj2,nc){
    return( do.call(rbind, sfClusterApplyLB(1:nc, function(i) eval(as.name(obj1))[[i]]%*%eval(as.name(obj2)))) )
}

I will run the matrix multiplication of two 1000x1000 matrices, 100times on a 8core cluster with 30GB RAM

The setting is

nr <- 1000
A <- matrix(round(rnorm(nr^2),1),nr=nr)
B <- t(A) + 4
ordinary <- A %*% B

This is what I see enter image description here

On the y-axis we see the time in seconds.

Clearly the sfCluster-Approach performs the best but is not feasible if the matrices are very huge (which in my case they are). So is it wise to go for the clusterCall-approach even though the foreach one is better? I'm not sure bot the foreach one does not copy everything to right?

Steve Weston Steve Weston · Accepted Answer · 2014-07-10T22:42:12

I would use clusterApply to distribute submatrices of A to the workers, and then use clusterCall to perform operations repeatedly on those submatrices. For example:

A <- matrix(rnorm(16), 4)
idx <- splitIndices(nrow(A), length(cl))
Alist <- lapply(idx, function(ii) A[ii,,drop=FALSE])
clusterApply(cl, Alist, function(a) { subA <<- a; NULL })

matprod.par <- function(cl, B) {
  do.call(rbind, clusterCall(cl, function(b) subA %*% b, B))
}

AB <- matprod.par(cl, matrix(1:16, 4))
AC <- matprod.par(cl, matrix(rnorm(16), 4))

Unlike clusterExport, clusterApply can export a different value of the variable "subA" for each of cluster workers.

The clusterCall function is very useful in this context since it allows you to iterate over data that has already been distributed to the workers, and you can still pass "B" along with each of the tasks.

Update

First, I should note that my example makes three assumptions:

The real problem involves more computations than a simple matrix multiply, since that is much better done sequentially on the master;
The "matprod.par" function will be executed multiple times making it worthwhile to pre-distribute the matrix "A" since it will be reused;
The "B" matrix is different on each call to "matprod.par", so no reuse is possible.

The "foreach" example that you present doesn't pre-distribute "A", so it can't reuse "A" if it is called multiple times. Both the "foreach" and "sfClusterApplyLB" examples copy all of "A" to all workers, which uses more memory, as you point out.

If you aren't going to call "matprod.par" multiple times, you could use:

matprod.par <- function(Alist, B) {
    foreach(a=Alist, .combine=rbind) %dopar% {
        a %*% B
    }
}

which avoids copying all of "A" to all the workers. This is equivalent to:

matprod.par <- function(cl, Alist, B) {
  do.call(rbind, clusterApply(cl, Alist, get('%*%'), B))
}

which will run somewhat faster than the foreach version since it has less overhead.

If you are going to call matprod.par many times, the pre-distribution time becomes insignificant, and the "clusterCall" example will be faster since none of "A" is sent to the workers again. Where the cross-over occurs depends on the number of times matprod.par is called, the size of "A", and the number of workers.

Parallel Design-Matrix, Parameter-Vector multiplication in ML-Estimation

1 Answers