computeSumFactors converts some sparse matrices to dense

[scottgigante]

Some formats of sparse matrix are causing massive allocation of memory. dgCMatrix works perfectly, using only a total of 2MB:

``` > library(scran) > library(lineprof) > > nrow = 5000 > ncol = 100000 > nnonzero = nrow*100 > > entry_index <- sample(nrow*ncol-1, nnonzero, replace=FALSE) > X <- Matrix::sparseMatrix(i=(entry_index %/% ncol)+1, j=(entry_index %% ncol)+1, x=1) > sce <- SingleCellExperiment(list(counts=as(X, "dgCMatrix"))) > > prof <- lineprof(computeSumFactors(sce)) dimnames(.) <- NULL: translated to dimnames(.) <- list(NULL,NULL) <==> unname(.) > prof Reducing depth to 2 (from 34) time alloc release dups ref src 1 0.046 2.245 0 639 c(".calculate_sum_factors", "ncol") .calculate_sum_factors/ncol 2 0.009 0.224 0 4 c(".calculate_sum_factors", ".limit_cluster_size") .calculate_sum_factors/.limit_cluster_size 3 0.006 0.273 0 1 ".calculate_sum_factors" .calculate_sum_factors 4 0.021 0.347 0 6 c(".calculate_sum_factors", ".limit_cluster_size") .calculate_sum_factors/.limit_cluster_size 5 0.013 0.103 0 3 c(".calculate_sum_factors", "split") .calculate_sum_factors/split 6 0.040 0.986 0 152 c(".calculate_sum_factors", ".guess_min_mean") .calculate_sum_factors/.guess_min_mean 7 0.041 0.616 0 530 c(".calculate_sum_factors", "bplapply") .calculate_sum_factors/bplapply ```

while dgTMatrix allocates a whopping 65MB for the same operation

``` > library(scran) > library(lineprof) > > nrow = 5000 > ncol = 100000 > nnonzero = nrow*100 > > entry_index <- sample(nrow*ncol-1, nnonzero, replace=FALSE) > X <- Matrix::sparseMatrix(i=(entry_index %/% ncol)+1, j=(entry_index %% ncol)+1, x=1) > sce <- SingleCellExperiment(list(counts=as(X, "dgTMatrix"))) > > prof <- lineprof(computeSumFactors(sce)) dimnames(.) <- NULL: translated to dimnames(.) <- list(NULL,NULL) <==> unname(.) > shine(prof) Starting interactive profile explorer. Press Escape / Ctrl + C to exit > prof Reducing depth to 2 (from 49) time alloc release dups ref src 1 0.031 2.516 0.000 659 c(".calculate_sum_factors", "ncol") .calculate_sum_factors/ncol 2 0.001 0.049 0.000 14 ".calculate_sum_factors" .calculate_sum_factors 3 0.006 0.176 0.000 1 c(".calculate_sum_factors", ".limit_cluster_size") .calculate_sum_factors/.limit_cluster_size 4 0.001 0.048 0.000 0 ".calculate_sum_factors" .calculate_sum_factors 5 0.001 0.049 0.000 1 c(".calculate_sum_factors", ".limit_cluster_size") .calculate_sum_factors/.limit_cluster_size 6 0.003 0.224 0.000 0 ".calculate_sum_factors" .calculate_sum_factors 7 0.019 0.347 0.000 6 c(".calculate_sum_factors", ".limit_cluster_size") .calculate_sum_factors/.limit_cluster_size 8 0.011 0.124 0.000 3 c(".calculate_sum_factors", "split") .calculate_sum_factors/split 9 0.094 1.874 0.000 183 c(".calculate_sum_factors", ".guess_min_mean") .calculate_sum_factors/.guess_min_mean 10 4.964 65.128 19.824 1874 c(".calculate_sum_factors", "bplapply") .calculate_sum_factors/bplapply ```

and the same issue for dgRMatrix, allocating ~40MB.

``` library(scran) library(lineprof) nrow = 100000 ncol = 5000 nnonzero = ncol*100 X <- Matrix::sparseMatrix(i=(entry_index %/% ncol)+1, j=(entry_index %% ncol)+1, x=1) X <- new("dgRMatrix", j=X@i, p=X@p, x=rep(1, nnonzero), Dim=as.integer(c(ncol, nrow))) sce <- SingleCellExperiment(list(counts=X)) prof <- lineprof(computeSumFactors(sce)) prof time alloc release dups ref src 1 0.011 2.304 0 659 c(".calculate_sum_factors", "ncol") .calculate_sum_factors/ncol 2 0.001 0.049 0 14 ".calculate_sum_factors" .calculate_sum_factors 3 0.011 0.224 0 1 c(".calculate_sum_factors", ".limit_cluster_size") .calculate_sum_factors/.limit_cluster_size 4 0.001 0.000 0 0 ".calculate_sum_factors" .calculate_sum_factors 5 0.001 0.049 0 1 c(".calculate_sum_factors", ".limit_cluster_size") .calculate_sum_factors/.limit_cluster_size 6 0.006 0.224 0 0 ".calculate_sum_factors" .calculate_sum_factors 7 0.028 0.347 0 6 c(".calculate_sum_factors", ".limit_cluster_size") .calculate_sum_factors/.limit_cluster_size 8 0.022 0.123 0 3 c(".calculate_sum_factors", "split") .calculate_sum_factors/split 9 0.156 1.462 0 288 c(".calculate_sum_factors", ".guess_min_mean") .calculate_sum_factors/.guess_min_mean 10 0.001 0.067 0 10 c(".calculate_sum_factors", ".subset2index") .calculate_sum_factors/.subset2index 11 0.698 38.735 0 1080 c(".calculate_sum_factors", "bplapply") .calculate_sum_factors/bplapply 12 0.001 0.000 0 16 character(0) ```

Session info:

``` > sessionInfo() R version 4.0.2 (2020-06-22) Platform: x86_64-pc-linux-gnu (64-bit) Running under: Arch Linux Matrix products: default BLAS/LAPACK: /usr/lib/libopenblas_haswellp-r0.3.10.so Random number generation: RNG: Mersenne-Twister Normal: Inversion Sample: Rounding locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] parallel stats4 stats graphics grDevices utils datasets methods base other attached packages: [1] lineprof_0.1.9001 scran_1.16.0 SingleCellExperiment_1.10.1 SummarizedExperiment_1.18.2 DelayedArray_0.14.1 matrixStats_0.56.0 [7] Biobase_2.48.0 GenomicRanges_1.40.0 GenomeInfoDb_1.24.2 IRanges_2.22.2 S4Vectors_0.26.1 BiocGenerics_0.34.0 loaded via a namespace (and not attached): [1] Rcpp_1.0.5 rsvd_1.0.3 locfit_1.5-9.4 lattice_0.20-41 R6_2.4.1 ggplot2_3.3.2 [7] pillar_1.4.6 zlibbioc_1.34.0 rlang_0.4.7 rstudioapi_0.11 irlba_2.3.3 Matrix_1.2-18 [13] BiocNeighbors_1.6.0 BiocParallel_1.22.0 statmod_1.4.34 stringr_1.4.0 igraph_1.2.5 RCurl_1.98-1.2 [19] munsell_0.5.0 beachmat_2.4.0 compiler_4.0.2 vipor_0.4.5 BiocSingular_1.4.0 pkgconfig_2.0.3 [25] ggbeeswarm_0.6.0 tidyselect_1.1.0 tibble_3.0.3 gridExtra_2.3 GenomeInfoDbData_1.2.3 edgeR_3.30.3 [31] viridisLite_0.3.0 crayon_1.3.4 dplyr_1.0.2 bitops_1.0-6 grid_4.0.2 gtable_0.3.0 [37] lifecycle_0.2.0 magrittr_1.5 scales_1.1.1 dqrng_0.2.1 stringi_1.5.3 XVector_0.28.0 [43] viridis_0.5.1 limma_3.44.3 scater_1.16.2 DelayedMatrixStats_1.10.1 ellipsis_0.3.1 generics_0.0.2 [49] vctrs_0.3.4 tools_4.0.2 glue_1.4.2 beeswarm_0.2.3 purrr_0.3.4 colorspace_1.4-1 ```

[LTLA]

Probably it's hitting beachmat code in C++, which doesn't know anything about dgTMatrix and dgRMatrix. All such "unknown" matrices are instead handled via block processing, which realizes blocks of the matrix as ordinary arrays with an upper memory usage defined by DelayedArray::getAutoBlockSize(). This defaults to 100 MB.

There are no plans to provide native support for dgTMatrix and dgRMatrix. The former can't be accessed efficiently and R-level operations convert them to dgCMatrix anyway. As for the latter, I have never seen it used in real analyses.

[scottgigante]

You may have never seen it but it's the reason why my jobs are failing, hence this issue report :) Seems like a fairly trivial fix to include a

if (is(X, "sparseMatrix")) X <- as(X, "CsparseMatrix")

rather than silently coerce the matrix to dense, hope the user notices and fixes it themself.

[LTLA]

It's not like the entire matrix is being coerced to dense form. It's bounded by a 100 MB limit, as defined by getAutoBlockSize; your jobs should be able to handle that. The same code is used to handle all unknown matrices, e.g., RleMatrix, ResidualMatrix, various other forms of DelayedMatrix representations: I don't see why dgTMatrixes and dgRMatrixes should get special treatment here.

[scottgigante]

In this toy example, the dgTMatrix also takes ~100x longer than the dgCMatrix. Seems like a good reason to me.

[LTLA]

One could say that about many matrix formats.

For example, I could store a sparse matrix in a HDF5Matrix, and under this reasoning, the function would be expected to convert it to a dgCMatrix for further processing. This is not a decision that the function should be allowed to make - for example, it would not be aware of memory constraints that motivated the use of the HDF5Matrix in the first place.

In the specific case of the dgTMatrix, an automated conversion is unwise if there are specific reasons for storing it as a dgTMatrix instead of a dgCMatrix. The most obvious is that the latter fails due to integer overflow in its p vector when there are more than ~2e9 non-zero elements in the matrix. The former stays operational at the cost of speed and memory.

Given that you're talking about toy examples, the easiest solution would just be to convert your matrix to the desired format. scran won't do any conversion, that is not its decision to make. For real analyses, all ingestion pipelines that I use will produce a dgCMatrix directly so I don't think this lack of automated conversion has much practical consequence.