Permute matrix columns in parallel with R

Question

I am trying to permute a matrix by columns using R. However, it is taking ages (the matrix is 68k x 32k integers).

I would like to do it in parallel (since each column is permuted independently). How can I achieve it? It should be something related to embarrassingly parallel for in R but I didn't find a solution.

Currently, my function is the following:

permMTX <- function(x) {
    nr <- nrow(x)
    nc <- ncol(x)
    # I'd like to parallelize this for, since each
    # column can be permuted independently
    for (i in 1:nc) {
        x[,i] <- x[sample(nr),i]
    }
    x
} 

Answer 1

Solution

First, I'd make use of vectorization, which should make it more efficient.

permMTX = function(x) apply(x, 2L, sample)

Then we can use library parallel to parallelize that function:

library(parallel)

parPermMTX = function(x, cluster) parApply(cl = cluster, X = x, MARGIN = 2L, FUN = sample)

---

Usage

With parallel you have to register a cluster before usage. Here's an example:

cl = makeCluster(detectCores(logical = FALSE))
parPermMTX(diag(10), cl)
#     [,1] [,2] [,3] [,4] [,5]
#[1,]    0    1    0    0    0
#[2,]    0    0    0    0    0
#[3,]    0    0    0    0    0
#[4,]    1    0    0    1    1
#[5,]    0    0    1    0    0

The way parallel works (spawning multiple R processes) you have to assure that you have enough memory to fit multiple copies of your data as well.

I think it's recommended to export the data to the processes as well, you can do that simply calling

clusterExport(cl, varlist = "exampleData")

While it does run in parallel on my end, it's not faster at all than simply employing apply, but I couldn't test with data at the same scale as yours, so I can't be sure it'll work.

This is due to the fact sample is heavily optimized already, so the overhead of spawning processes is bigger than simply calling sample. See Why is the parallel package slower than just using apply?

On my system, sampling 68E3 integers 32E3 times takes roughly 40 seconds:

microbenchmark(sample(68E3), times = 32E3)
#Unit: milliseconds
#          expr      min       lq     mean   median       uq      max neval
# sample(68000) 1.132273 1.192923 1.290838 1.227912 1.286229 7.880191 32000

Perhaps you're running out of memory, and using the hard disk cache, which is really slow.

---

2nd Solution

So, what if we tried to allocate as many calls to sample sequentially to a single process? This is what I tried here:

parPermMTX2 = function(x, cluster) do.call(cbind, parLapply(cl = cluster, X = list(x[,seq(floor(ncol(x)/2))], x[,floor(ncol(x)/2)+seq(ceiling(ncol(x)/2))]), fun = permMTX))

We divide x into two halves, then call permMTX in each, then recombine with cbind.

Sadly, neither this way I could achieve better performance. So, while I answered your question, I'm not sure it's any help at all.

Answer 2

Disclaimer: I'm the author of package bigstatsr.

You can use shared memory (matrices stored on disk) and do this:

# devtools::install_github("privefl/bigstatsr")
library(bigstatsr)

# matrix on disk
mat <- FBM(68e3, 32e2, backingfile = "test")
# inialize with 1:nrow(mat) for each column
system.time(
  big_apply(mat, a.FUN = function(X, ind) {
    print(min(ind))
    X[, ind] <- rep(rows_along(X), length(ind))
    NULL
  }, a.combine = 'c')
) # 15 sec

# permute each column, in parallel
system.time(
  big_apply(mat, a.FUN = function(X, ind) {
    print(min(ind))
    X[, ind] <- apply(X[, ind], 2, sample)
    NULL
  }, a.combine = 'c', ncores = nb_cores())
) # 27 sec

This takes 27 secs on a tenth of the data and 378 sec on the whole dataset (on a laptop computer with only 2 physical cores and 8GB of RAM).