Sparse vs dense representation

The datasets that we will use for our benchmarks are subsets of the 1.3 Million Brain Cell Dataset from 10x Genomics. This is a sparse 27,998 x 1,306,127 matrix of counts, with one gene per row and one cell per column. Around 7% of the matrix values are nonzero counts. The dataset is available via the ExperimentHub package in two forms:

1. As a sparse HDF5 file: This is the original HDF5 file provided by 10x Genomics. It uses the CSR/CSC/Yale representation to store the sparse data.

2. As a dense HDF5 file: The same data as the above but stored as a regular HDF5 dataset with (compressed) chunks of dimensions 100 x 100.

The two files are hosted on ExperimentHub under resource ids EH1039 and EH1040:

hub <- ExperimentHub()
hub["EH1039"]$description  # sparse representation

## [1] "Single-cell RNA-seq data for 1.3 million brain cells from E18 mice. 'HDF5-based 10X Genomics' format originally provided by TENx Genomics"

hub["EH1040"]$description  # dense representation

## [1] "Single-cell RNA-seq data for 1.3 million brain cells from E18 mice. Full rectangular, block-compressed format, 1GB block size."

Let’s download them to the local ExperimentHub cache if they are not there yet:

## Note that this will be quick if the HDF5 files are already in the
## local ExperimentHub cache. Otherwise, it will take a while!
brain_s_path <- hub[["EH1039"]]
brain_D_path <- hub[["EH1040"]]

brain_s_path and brain_D_path are the paths to the downloaded files.

TENxMatrix vs HDF5Matrix objects

We use the TENxMatrix() and HDF5Array() constructors to bring the sparse and dense datasets in R, as DelayedArray derivatives. Note that this does not load the matrix data in memory.

## Use 'h5ls(brain_s_path)' to find out the group.
brain_s <- TENxMatrix(brain_s_path, group="mm10")

class(brain_s)

## [1] "TENxMatrix"
## attr(,"package")
## [1] "HDF5Array"

dim(brain_s)

## [1]   27998 1306127

is_sparse(brain_s)

## [1] TRUE

## Use 'h5ls(brain_D_path)' to find out the name of the dataset.
brain_D <- HDF5Array(brain_D_path, name="counts")

class(brain_D)

## [1] "HDF5Matrix"
## attr(,"package")
## [1] "HDF5Array"

dim(brain_D)

## [1]   27998 1306127

chunkdim(brain_D)

## [1] 100 100

is_sparse(brain_D)

## [1] FALSE

brain_Ds <- HDF5Array(brain_D_path, name="counts", as.sparse=TRUE)

class(brain_Ds)

## [1] "HDF5Matrix"
## attr(,"package")
## [1] "HDF5Array"

dim(brain_Ds)

## [1]   27998 1306127

chunkdim(brain_Ds)

## [1] 100 100

is_sparse(brain_Ds)

## [1] TRUE

dimnames(brain_Ds) <- dimnames(brain_D) <- dimnames(brain_s)

Create the test datasets

For our benchmarks below, we create subsets of the 1.3 Million Brain Cell Dataset of increasing sizes: subsets with 12,500 cells, 25,000 cells, 50,000 cells, 100,000 cells, and 200,000 cells. Note that subsetting a TENxMatrix or HDF5Matrix object with [ is a delayed operation so has virtually no cost:

brain1_s  <- brain_s[  , 1:12500]
brain1_D  <- brain_D[  , 1:12500]
brain1_Ds <- brain_Ds[ , 1:12500]

brain2_s  <- brain_s[  , 1:25000]
brain2_D  <- brain_D[  , 1:25000]
brain2_Ds <- brain_Ds[ , 1:25000]

brain3_s  <- brain_s[  , 1:50000]
brain3_D  <- brain_D[  , 1:50000]
brain3_Ds <- brain_Ds[ , 1:50000]

brain4_s  <- brain_s[  , 1:100000]
brain4_D  <- brain_D[  , 1:100000]
brain4_Ds <- brain_Ds[ , 1:100000]

brain5_s  <- brain_s[  , 1:200000]
brain5_D  <- brain_D[  , 1:200000]
brain5_Ds <- brain_Ds[ , 1:200000]

Code used for normalization and PCA

We’ll use the following code for normalization:

## Also selects the most variable genes (1000 by default).
simple_normalize <- function(mat, num_var_genes=1000)
{
    stopifnot(length(dim(mat)) == 2, !is.null(rownames(mat)))
    mat <- mat[rowSums(mat) > 0, ]
    col_sums <- colSums(mat) / 10000
    mat <- t(t(mat) / col_sums)
    row_vars <- rowVars(mat)
    row_vars_order <- order(row_vars, decreasing=TRUE)
    variable_idx <- head(row_vars_order, n=num_var_genes)
    mat <- log1p(mat[variable_idx, ])
    mat / rowSds(mat)
}

and the following code for PCA:

simple_PCA <- function(mat, k=25)
{
    stopifnot(length(dim(mat)) == 2)
    row_means <- rowMeans(mat)
    Ax <- function(x, args)
        (as.numeric(mat %*% x) - row_means * sum(x))
    Atx <- function(x, args)
        (as.numeric(x %*% mat) - as.vector(row_means %*% x))
    RSpectra::svds(Ax, Atrans=Atx, k=k, dim=dim(mat))
}

Block processing and block size

Note that the implementations of simple_normalize() and simple_PCA() are expected to work on any matrix-like object regardless of its exact type/representation e.g. it can be an ordinary matrix, a SparseMatrix object from the SparseArray package, a dgCMatrix object from the Matrix package, a DelayedMatrix derivative (TENxMatrix, HDF5Matrix, TileDBMatrix), etc…

However, when the input is a DelayedMatrix object or derivative, it’s important to be aware that:

• Summarization methods like sum(), colSums(), rowVars(), or rowSds(), and matrix multiplication (%*%), are block-processed operations.

• The block size is 100 Mb by default. Increasing or decreasing the block size will typically increase or decrease the memory usage of block-processed operations. It will also impact performance, but sometimes in unexpected or counter-intuitive ways.

• The block size can be controlled with DelayedArray::getAutoBlockSize() and DelayedArray::setAutoBlockSize().

For our benchmarks below, we’ll use the following block sizes:

Monitoring memory usage

While manually running our benchmarks below on a Linux or macOS system, we will also monitor memory usage at the command line in a terminal with:

(while true; do ps u -p <PID>; sleep 1; done) >ps.log 2>&1 &

where <PID> is the process id of our R session. This will allow us to measure the maximum amount of memory used by the calls to simple_normalize() or simple_PCA().

Normalization

In this section we run simple_normalize() on the three different representations (TENxMatrix, HDF5Matrix, and “HDF5Matrix as sparse”) of the smaller test dataset only (27,998 x 12,500), and we report the time of each run.

See the Comprehensive timings obtained on various machines section below in this document for simple_normalize() and simple_pca() timings obtained on various machines on all our test datasets and using four different block sizes: 16 Mb, 40 Mb, 100 Mb, and 250 Mb.

dim(brain1_s)

## [1] 27998 12500

DelayedArray::setAutoBlockSize(250e6)  # set block size to 250 Mb

## automatic block size set to 2.5e+08 bytes (was 1e+08)

system.time(norm_brain1_s <- simple_normalize(brain1_s))

##    user  system elapsed 
##  46.243   0.993  48.820

dim(norm_brain1_s)

## [1]  1000 12500

dim(brain1_D)

## [1] 27998 12500

DelayedArray::setAutoBlockSize(16e6)   # set block size to 16 Mb

## automatic block size set to 1.6e+07 bytes (was 2.5e+08)

system.time(norm_brain1_D <- simple_normalize(brain1_D))

##    user  system elapsed 
##  57.787   0.609  58.421

dim(norm_brain1_D)

## [1]  1000 12500

dim(brain1_Ds)

## [1] 27998 12500

DelayedArray::setAutoBlockSize(250e6)  # set block size to 250 Mb

## automatic block size set to 2.5e+08 bytes (was 1.6e+07)

system.time(norm_brain1_Ds <- simple_normalize(brain1_Ds))

##    user  system elapsed 
##  45.028   0.881  45.922

dim(norm_brain1_Ds)

## [1]  1000 12500

On-disk realization of the normalized datasets

Note that the normalized datasets obtained in the previous section are DelayedMatrix objects that carry delayed operations. These operations can be displayed with showtree() e.g. for norm_brain1_s:

class(norm_brain1_s)

## [1] "DelayedMatrix"
## attr(,"package")
## [1] "DelayedArray"

showtree(norm_brain1_s)

## 1000x12500 double, sparse: DelayedMatrix object
## └─ 1000x12500 double, sparse: Unary iso op with args
##    └─ 1000x12500 double, sparse: Stack of 1 unary iso op(s)
##       └─ 1000x12500 double, sparse: Aperm (perm=c(2,1))
##          └─ 12500x1000 double, sparse: Unary iso op with args
##             └─ 12500x1000 integer, sparse: Aperm (perm=c(2,1))
##                └─ 1000x12500 integer, sparse: Subset
##                   └─ 27998x1306127 integer, sparse: [seed] TENxMatrixSeed object

The other norm_brain1_* objects carry similar operations.

Before we proceed with PCA, we’re going to write the normalized datasets to new HDF5 files. This introduces an additional cost, but, on the other hand, it has the benefit of triggering on-disk realization of the object. This means that all the delayed operations carried by the object will get realized on-the-fly before the matrix data actually lands on the disk, making the new object (TENxMatrix or HDF5Matrix) more efficient for PCA or whatever block-processed operations will come next.

We will use blocks of 100 Mb for all the writing operations.

DelayedArray::setAutoBlockSize(1e8)

## automatic block size set to 1e+08 bytes (was 2.5e+08)

dim(norm_brain1_s)

## [1]  1000 12500

system.time(norm_brain1_s <- writeTENxMatrix(norm_brain1_s, level=0))

##    user  system elapsed 
##   3.265   0.168   3.444

class(norm_brain1_s)

## [1] "TENxMatrix"
## attr(,"package")
## [1] "HDF5Array"

showtree(norm_brain1_s)  # "pristine" object (i.e. no more delayed operations)

## 1000x12500 double, sparse: TENxMatrix object
## └─ 1000x12500 double, sparse: [seed] TENxMatrixSeed object

dim(norm_brain1_D)

## [1]  1000 12500

system.time(norm_brain1_D <- writeHDF5Array(norm_brain1_D, level=0))

##    user  system elapsed 
##   2.883   0.227   3.110

class(norm_brain1_D)

## [1] "HDF5Matrix"
## attr(,"package")
## [1] "HDF5Array"

showtree(norm_brain1_D)  # "pristine" object (i.e. no more delayed operations)

## 1000x12500 double: HDF5Matrix object
## └─ 1000x12500 double: [seed] HDF5ArraySeed object

dim(norm_brain1_Ds)

## [1]  1000 12500

system.time(norm_brain1_Ds <- writeHDF5Array(norm_brain1_Ds, level=0))

##    user  system elapsed 
##   6.239   0.247   6.486

class(norm_brain1_Ds)

## [1] "HDF5Matrix"
## attr(,"package")
## [1] "HDF5Array"

showtree(norm_brain1_Ds)  # "pristine" object (i.e. no more delayed operations)

## 1000x12500 double, sparse: HDF5Matrix object
## └─ 1000x12500 double, sparse: [seed] HDF5ArraySeed object

PCA

In this section we run simple_pca() on the normalized datasets obtained in the previous section and report the time of each run.

See the Comprehensive timings obtained on various machines section below in this document for simple_normalize() and simple_pca() timings obtained on various machines on all our test datasets and using four different block sizes: 16 Mb, 40 Mb, 100 Mb, and 250 Mb.

DelayedArray::setAutoBlockSize(40e6)   # set block size to 40 Mb

## automatic block size set to 4e+07 bytes (was 1e+08)

dim(norm_brain1_s)

## [1]  1000 12500

system.time(pca1s <- simple_PCA(norm_brain1_s))

##    user  system elapsed 
##  53.483   1.817  55.304

DelayedArray::setAutoBlockSize(1e8)    # set block size to 100 Mb

## automatic block size set to 1e+08 bytes (was 4e+07)

dim(norm_brain1_D)

## [1]  1000 12500

system.time(pca1D <- simple_PCA(norm_brain1_D))

##    user  system elapsed 
##  69.504  84.533  59.402

stopifnot(all.equal(pca1D, pca1s))

DelayedArray::setAutoBlockSize(40e6)   # set block size to 40 Mb

## automatic block size set to 4e+07 bytes (was 1e+08)

dim(norm_brain1_Ds)

## [1]  1000 12500

system.time(pca1Ds <- simple_PCA(norm_brain1_Ds))

##    user  system elapsed 
## 158.321   3.779 162.195

stopifnot(all.equal(pca1Ds, pca1s))

Timings for DELL XPS 15 laptop

Timings for Supermicro SuperServer 1029GQ-TRT

Timings for Apple Silicon Mac Pro

Timings for Intel Mac Pro

TENxMatrix (sparse) vs HDF5Matrix (dense)

Main takeaways

1. Using the TENxMatrix representation (sparse format), one can perform normalization and PCA of the bigger dataset (200,000 cells) on an average consumer-grade laptop like the DELL XPS 15 laptop (model 9520) in less than 25 minutes and using less than 4Gb of memory, as shown in the table below. For comparison, normalization and PCA on an in-memory representation of the same dataset (e.g. on a SparseMatrix object) takes less then 4 minutes. However, that consumes about 18.5Gb of memory! This will be the subject of an upcoming vignette in the SparseArray package.

2. Normalization and PCA are roughly linear in time with respect to the number of cells in the dataset, regardless of representation (sparse or dense) or block size.

3. Block size matters. When using the TENxMatrix representation (sparse format), the bigger the blocks the faster normalization will be (at the cost of increased memory usage). On the other hand, it seems that PCA prefers small blocks, at least with our naive simple_PCA() implementation.

4. Disk performance is important (not surprisingly) as attested by the lower performance of the Supermicro SuperServer 1029GQ-TRT machine, likely due to its slower disk.