Bug in rowSums,DelayedMatrix-method

[PeteHaitch]

Calling rowSums() on a large-ish DelayedMatrix leads to a serialization/forking/memory issue on macOS and Linux

library(DelayedArray)
x <- DelayedArray(matrix(1L, nrow = 10000000, ncol = 100))

# Errors
rowSums(x)

On macOS (16GB RAM) the error is:

Error: vector memory exhausted (limit reached?)

On Linux (20GB RAM) the error is:

Error in serialize(data, node$con, xdr = FALSE) : ignoring SIGPIPE signal
Error in serialize(data, node$con, xdr = FALSE) : ignoring SIGPIPE signal
Error: failed to stop ‘SOCKcluster’ cluster: error writing to connection
Error in serialize(data, node$con, xdr = FALSE) : ignoring SIGPIPE signal
# Then "Error in serialize(data, node$con, xdr = FALSE) : ignoring SIGPIPE signal" repeats indefinitely

I thought it might be a more general issue with blockApply() and its use of BiocParallel, but I haven't been able to trigger the problem in some brief testing. For example, using colSums() or blockApply()-ing max() over individual columns or rows of x worked fine.

[ttriche]

nb. I see this problem even on Linux boxes with 384GB of RAM.

It does not appear to be a small-memory issue per se.

[PeteHaitch]

@ttriche can you please try DelayedMatrixStats::rowSums2() and let me know if you still have the issue? It's implemented slightly differently

[ttriche]

I use rowSums2() in my code, but the trigger for this behavior seems to be inside of the `dmrseq` package, and I haven't been able to track it down yet. Running some separate tests to try and figure it out now.

…

--t

On Tue, Apr 24, 2018 at 1:31 PM, Peter Hickey ***@***.***> wrote: @ttriche <https://github.com/ttriche> can you please try DelayedMatrixStats::rowSums2() and let me know if you still have the issue? It's implemented slightly differently — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#16 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAARIu1AEpx1kOPsQnTpbFp-T4kA7gtfks5tr2GDgaJpZM4TR-hQ> .

[kdkorthauer]

@PeteHaitch Thanks for looking into this. Should the problem be fixed in hansenlab/bsseq with the most recent commit?

@ttriche If you try installing the latest bsseq from hansenlab/bsseq, does dmrseq run properly again?

[PeteHaitch]

@kdkorthauer Yes, it should be fixed now in bsseq. However, it's not clear to me if @ttriche's dmrseq example is triggering this issue from within bsseq or elsewhere in dmrseq. I'm happy to help debug this further

[kdkorthauer]

@PeteHaitch awesome, thanks! I'll investigate whether the issue persists within dmrseq, and if so whether a similar patch might work.

[mtmorgan]

DelayedArray does this

    bplapply(seq_len(nblock),
        function(b) {
            if (get_verbose_block_processing())
                message("Processing block ", b, "/", nblock, " ... ",
                        appendLF=FALSE)
            viewport <- grid[[b]]
            block <- extract_block(x, viewport)
            if (!is.array(block))
                block <- .as_array_or_matrix(block)
            attr(block, "from_grid") <- grid
            attr(block, "block_id") <- b
            block_ans <- FUN(block, ...)
            if (get_verbose_block_processing())
                message("OK")
            block_ans
        },
        BPREDO=BPREDO,
        BPPARAM=BPPARAM
    )

where block <- extract_block(x, viewport) is done on the worker. This means that x needs to be made available (serialized to, even in the case of MulticoreParam()) on the worker. A different implementation is to use bpiterate() and a generator function ITER to produce blocks on the manager.

    ITER <- local({
        b <- 0L
        function() {
            b <<- b + 1L
            if (b > nblock)
                return(NULL)
            if (get_verbose_block_processing())
                message("Processing block ", b, "/", nblock, " ... ",
                        appendLF=FALSE)
            viewport <- grid[[b]]
            block <- extract_block(x, viewport)
            if (!is.array(block))
                block <- .as_array_or_matrix(block)
            attr(block, "from_grid") <- grid
            attr(block, "block_id") <- b
            block
        }
    })
            
    bpiterate(ITER, FUN, BPPARAM = BPPARAM)

This 'works' but is incredibly slow (use set_verbose_block_processing(TRUE) to convince yourself that it's chugging away) because the chunks are still serialized to each worker, and because the garbage collector is being called often; I'll explore a better solution for the common multicore data transfer problem in BiocParallel. Also I'm not sure where, if x were something like an HDF5Array, the object is actually realized in memory as a matrix; one would like to do that step on the worker.

I say 'works', but actually on my laptop (after setting BiocParallel::register(BiocParallel::SerialParam()) for better speed) in the rowSums,DelayedArray method there is

    block_results <- blockApply(x, rowSums, na.rm=na.rm)

    ans <- rowSums(matrix(unlist(block_results, use.names=FALSE), nrow=nrow(x)))

and on the last line I end up with Error: cannot allocate vector of size 7.5 Gb -- the original object x is consuming 3.7G, block_results consumes 7.5G (these are doubles, rather than ints, which is a little surprising, I would have thought rowSums() would have returned ints if possible @lawremi), unlist(block_results, use.names=FALSE) is another 7.5G, and the matrix(...) is another 7.5G so 3.7G + 3 * 7.5G. At this point the unlist() result is available for garbage collection (is it?), but then ans takes its place using another 7.5G! If this has been written as

    result <- blockApply(x, rowSums, na.rm=na.rm)
    result <- unlist(result, use.names = FALSE)
    result <- matrix(result, nrow = nrow(x))
    result <- rowSums(result)

there would only ever need to be 3.7G + 2 x 7.5G in memory.

Hmm, but now I'm confused! unlist(blockApply()) is as.vector(x) (!), which we reshape into a matrix equal to x (except with numeric rather than integer type) and then we calculate rowSums locally! I guess DelayedArray has chosen to block into chunks where each chunk has a single column...

[hpages]

unlist(blockApply()) is as.vector(x) because blocApply() is still using the old default block grid where the blocks go "along the columns". This is not optimal in most cases and needs to change. Will do ASAP.