1 Introduction

The aim of this document is to measure the performance of the HDF5Array package for normalization and PCA (Principal Component Analysis) of on-disk single cell data, two computationally intensive operations at the core of single cell analysis.

The benchmarks presented in the document were specifically designed to observe the impact of two critical parameters on performance:

data representation (i.e. sparse vs dense);
size of the blocks used for block-processed operations.

Hopefully these benchmarks will also facilitate comparing performance of single cell analysis workflows based on HDF5Array with workflows based on other tools like Seurat or Scanpy.

2 Install and load the required packages

Let’s install and load HDF5Array as well as the other packages used in this vignette:

if (!require("BiocManager", quietly=TRUE))
    install.packages("BiocManager")

pkgs <- c("HDF5Array", "ExperimentHub", "DelayedMatrixStats", "RSpectra")
BiocManager::install(pkgs)

Load the packages:

library(HDF5Array)
library(ExperimentHub)
library(DelayedMatrixStats)
library(RSpectra)

3 The test datasets

3.1 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:

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.
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.

3.2 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.

3.2.1 Bring the sparse dataset in R

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

brain_s is a 27,998 x 1,306,127 TENxMatrix object:

class(brain_s)

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

dim(brain_s)

## [1]   27998 1306127

is_sparse(brain_s)

## [1] TRUE

See ?TENxMatrix in the HDF5Array package for more information about TENxMatrix objects.

3.2.2 Bring the dense dataset in R

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

brain_D is a 27,998 x 1,306,127 HDF5Matrix object that contains the same data as brain_s:

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

See ?HDF5Matrix in the HDF5Array package for more information about HDF5Matrix objects.

Even though the data in brain_D_path is stored in a dense format, we can flag it as quantitatively sparse. This is done by calling the HDF5Array() constructor function with as.sparse=TRUE:

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

The only difference between brain_D and brain_Ds is that the latter is now seen as a sparse object, and will be treated as such:

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

Concretely this means that, when blocks of data are loaded from the dense HDF5 file to memory during block-processed operations, they end up directly in an in-memory sparse representation without going thru an in-memory dense representation first. This is expected to reduce memory footprint and (hopefully) will help with overall performance.

Finally note that the dense HDF5 file does not contain the dimnames of the matrix, so we manually add them to brain_s and brain_Ds:

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

3.3 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]

4 Block-processed normalization and PCA

4.1 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))
}

4.2 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:

	NORMALIZATION	ON-DISK REALIZATION	PCA
TENxMatrix (sparse)	250 Mb	100 Mb	40 Mb
HDF5Matrix (dense)	16 Mb	100 Mb	100 Mb
HDF5Matrix as sparse	250 Mb	100 Mb	40 Mb

4.3 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().

5 Benchmarks

5.1 Normalization

In this section we run simple_normalize() on the two smaller test datasets only (27,998 x 12,500 and 27,998 x 25,000), on the three different representations (TENxMatrix, HDF5Matrix, and “HDF5Matrix as sparse”) for each of them. We report the time and memory usage for each run.

See the Comprehensive timings obtained on various machines section at the end of 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.

5.1.1 Normalization of the 27,998 x 12,500 dataset

5.1.1.1 TENxMatrix (sparse)

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 
##  30.934   1.945  32.856

dim(norm_brain1_s)

## [1]  1000 12500

5.1.1.2 HDF5Matrix (dense)

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 
##  38.882   0.594  39.477

dim(norm_brain1_D)

## [1]  1000 12500

5.1.1.3 HDF5Matrix as sparse

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 
##  32.840   1.671  34.495

dim(norm_brain1_Ds)

## [1]  1000 12500

5.1.2 Normalization of the 27,998 x 25,000 dataset

5.1.2.1 TENxMatrix (sparse)

dim(brain2_s)

## [1] 27998 25000

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

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

system.time(norm_brain2_s <- simple_normalize(brain2_s))

##    user  system elapsed 
##  62.547   2.882  65.392

dim(norm_brain2_s)

## [1]  1000 25000

5.1.2.2 HDF5Matrix (dense)

dim(brain2_D)

## [1] 27998 25000

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_brain2_D <- simple_normalize(brain2_D))

##    user  system elapsed 
##  85.198   1.262  86.461

dim(norm_brain2_D)

## [1]  1000 25000

5.1.2.3 HDF5Matrix as sparse

dim(brain2_Ds)

## [1] 27998 25000

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_brain2_Ds <- simple_normalize(brain2_Ds))

##    user  system elapsed 
##  71.670   3.206  74.839

dim(norm_brain2_Ds)

## [1]  1000 25000

5.2 On-disk realization of the normalized datasets

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

class(norm_brain2_s)

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

showtree(norm_brain2_s)

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

All the other norm_brain* 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.

5.2.1 On-disk realization of the 1000 x 12,500 normalized dataset

5.2.1.1 TENxMatrix (sparse)

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

## 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 
##   2.405   0.156   2.563

The new norm_brain1_s object is a pristine TENxMatrix object:

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

5.2.1.2 HDF5Matrix (dense)

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

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

dim(norm_brain1_D)

## [1]  1000 12500

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

##    user  system elapsed 
##    2.41    0.18    2.59

The new norm_brain1_D object is a pristine HDF5Matrix object:

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

5.2.1.3 HDF5Matrix as sparse

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

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

dim(norm_brain1_Ds)

## [1]  1000 12500

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

##    user  system elapsed 
##   3.592   0.175   3.767

The new norm_brain1_Ds object is a pristine sparse HDF5Matrix object:

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

5.2.2 On-disk realization of the 1000 x 25,000 normalized dataset

5.2.2.1 TENxMatrix (sparse)

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

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

dim(norm_brain2_s)

## [1]  1000 25000

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

##    user  system elapsed 
##   6.571   0.240   6.811

The new norm_brain2_s object is a pristine TENxMatrix object:

class(norm_brain2_s)

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

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

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

5.2.2.2 HDF5Matrix (dense)

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

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

dim(norm_brain2_D)

## [1]  1000 25000

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

##    user  system elapsed 
##   4.035   0.323   4.358

The new norm_brain2_D object is a pristine HDF5Matrix object:

class(norm_brain2_D)

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

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

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

5.2.2.3 HDF5Matrix as sparse

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

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

dim(norm_brain2_Ds)

## [1]  1000 25000

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

##    user  system elapsed 
##   8.552   0.246   8.797

The new norm_brain2_Ds object is a pristine sparse HDF5Matrix object:

class(norm_brain2_Ds)

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

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

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

5.3 PCA

In this section we run simple_pca() on the normalized datasets obtained in the previous section (1000 x 12,500 and 1000 x 25,000), on the three different representations (TENxMatrix, HDF5Matrix, and “HDF5Matrix as sparse”) for each dataset. We report the time and memory usage for each run.

See the Comprehensive timings obtained on various machines section at the end of 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.

5.3.1 PCA on the 1000 x 12,500 normalized dataset

5.3.1.1 TENxMatrix (sparse)

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 
##  32.720   2.170  33.418

5.3.1.2 HDF5Matrix (dense)

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 
##  34.332   8.815  43.150

Sanity check:

stopifnot(all.equal(pca1s, pca1D))

5.3.1.3 HDF5Matrix as sparse

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 
## 112.941   5.684 117.104

Sanity check:

stopifnot(all.equal(pca1s, pca1Ds))

5.3.2 PCA on the 1000 x 25,000 normalized dataset

5.3.2.1 TENxMatrix (sparse)

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

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

dim(norm_brain2_s)

## [1]  1000 25000

system.time(pca2s <- simple_PCA(norm_brain2_s))

##    user  system elapsed 
##  70.938   4.141  72.055

5.3.2.2 HDF5Matrix (dense)

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

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

dim(norm_brain2_D)

## [1]  1000 25000

system.time(pca2D <- simple_PCA(norm_brain2_D))

##    user  system elapsed 
##  56.868  16.676  73.548

Sanity check:

stopifnot(all.equal(pca2s, pca2D))

5.3.2.3 HDF5Matrix as sparse

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

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

dim(norm_brain2_Ds)

## [1]  1000 25000

system.time(pca2Ds <- simple_PCA(norm_brain2_Ds))

##    user  system elapsed 
## 220.943  11.191 229.004

Sanity check:

stopifnot(all.equal(pca2s, pca2Ds))

5.4 Comprehensive timings obtained on various machines

Here we report timings (and memory usage) observed on various machines. For each machine, the results are presented in a table that shows the normalization & realization & PCA timings obtained for all our test datasets and using four different block sizes: 16 Mb, 40 Mb, 100 Mb, and 250 Mb. For each operation, the best time across the four different block sizes is displayed in bold.

All the timings (and memory usage) were produced by running the run_benchmarks.sh script located in the HDF5Array/inst/scripts/ folder of the package, using R 4.5 and HDF5Array 1.35.12 (Bioconductor 3.21).

5.4.1 Timings for DELL XPS 15 laptop

Machine specs:

DELL XPS 15 laptop (model 9520)
OS	Linux Ubuntu 24.04	Disk performance	Output of `sudo hdparm -Tt <device>`: `Timing cached reads: 35188 MB in 2.00 seconds = 17620.75 MB/sec Timing buffered disk reads: 7842 MB in 3.00 seconds = 2613.57 MB/sec`
RAM	32GB
Disk	1TB SSD

Timings:

Test Dataset	F o r m a t	block size = 16 Mb		block size = 40 Mb		block size = 100 Mb		block size = 250 Mb		Normalized Test Dataset	F o r m a t	block size = 16 Mb		block size = 40 Mb		block size = 100 Mb		block size = 250 Mb		block size = 16 Mb		block size = 40 Mb		block size = 100 Mb		block size = 250 Mb
nrow x ncol (# genes x # cells)	F o r m a t	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	nrow x ncol (# sel. genes x # cells)	F o r m a t	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used

		1. NORMALIZATION & selection of 1000 most variable genes										2. ON-DISK REALIZATION of the normalized dataset								3. PCA of the normalized dataset

27998 x 12500	[s]	74	1.2Gb	57	1.2Gb	40	1.6Gb	34	1.9Gb	1000 x 12500	[s]	3	1.1Gb	3	1.3Gb	4	1.6Gb	4	1.7Gb	41	1.2Gb	39	1.3Gb	51	1.5Gb	39	1.6Gb
	[D]	51	1.3Gb	50	1.4Gb	52	1.6Gb	46	2.2Gb		[D]	3	1.3Gb	2	1.2Gb	2	1.7Gb	2	2.2Gb	49	1.4Gb	41	1.4Gb	55	1.3Gb	35	1.7Gb
	[Ds]	46	1.2Gb	47	1.3Gb	37	1.6Gb	36	2.2Gb		[Ds]	6	1.2Gb	5	1.3Gb	4	1.6Gb	4	2.1Gb	135	1.4Gb	119	1.7Gb	159	1.8Gb	132	2.0Gb

27998 x 25000	[s]	152	1.2Gb	122	1.3Gb	90	1.7Gb	63	2.3Gb	1000 x 25000	[s]	6	1.2Gb	6	1.3Gb	7	1.9Gb	7	2.1Gb	84	1.2Gb	76	1.5Gb	94	1.7Gb	76	2.0Gb
	[D]	94	1.3Gb	98	1.5Gb	114	1.7Gb	99	2.2Gb		[D]	7	1.3Gb	4	1.3Gb	5	1.5Gb	4	1.8Gb	94	1.4Gb	78	1.5Gb	81	1.7Gb	130	1.5Gb
	[Ds]	91	1.3Gb	89	1.4Gb	88	1.9Gb	75	2.3Gb		[Ds]	11	1.3Gb	8	1.4Gb	8	1.8Gb	9	2.3Gb	241	1.5Gb	221	2.0Gb	285	2.0Gb	244	2.6Gb

27998 x 50000	[s]	328	1.2Gb	238	1.5Gb	178	2.0Gb	140	2.5Gb	1000 x 50000	[s]	13	1.2Gb	15	1.4Gb	13	1.8Gb	12	2.5Gb	165	1.2Gb	156	1.4Gb	193	1.8Gb	176	2.2Gb
	[D]	183	1.4Gb	198	1.5Gb	219	1.7Gb	209	2.4Gb		[D]	17	1.3Gb	13	1.3Gb	9	1.5Gb	9	1.6Gb	241	1.4Gb	215	1.5Gb	189	1.6Gb	167	2.2Gb
	[Ds]	183	1.2Gb	174	1.5Gb	167	1.9Gb	165	2.5Gb		[Ds]	27	1.2Gb	20	1.5Gb	16	1.7Gb	15	2.4Gb	535	1.6Gb	474	1.9Gb	620	1.9Gb	472	3.6Gb

27998 x 100000	[s]	642	1.2Gb	458	1.6Gb	365	2.0Gb	280	2.9Gb	1000 x 100000	[s]	23	1.3Gb	25	1.5Gb	25	1.8Gb	23	2.7Gb	304	1.3Gb	268	1.5Gb	305	1.7Gb	332	2.6Gb
	[D]	344	1.4Gb	392	1.6Gb	439	1.7Gb	417	2.3Gb		[D]	34	1.3Gb	25	1.3Gb	18	1.6Gb	17	1.9Gb	520	1.5Gb	424	1.6Gb	260	1.7Gb	310	2.2Gb
	[Ds]	354	1.2Gb	340	1.5Gb	323	2.0Gb	334	3.0Gb		[Ds]	52	1.3Gb	40	1.5Gb	34	1.9Gb	30	3.1Gb	1057	1.7Gb	1016	2.0Gb	1080	2.2Gb	944	3.9Gb

27998 x 200000	[s]	1284	1.3Gb	985	1.6Gb	738	2.0Gb	610	3.7Gb	1000 x 200000	[s]	46	1.3Gb	47	1.6Gb	52	2.0Gb	46	3.4Gb	585	1.3Gb	535	1.6Gb	703	1.8Gb	665	2.9Gb
	[D]	694	1.4Gb	796	1.9Gb	886	1.7Gb	855	3.0Gb		[D]	78	1.3Gb	62	1.8Gb	49	1.6Gb	36	2.8Gb	1588	1.6Gb	1088	2.1Gb	892	1.8Gb	676	3.0Gb
	[Ds]	775	1.3Gb	680	1.6Gb	642	2.4Gb	677	3.6Gb		[Ds]	113	1.4Gb	93	1.6Gb	74	2.2Gb	62	3.4Gb	2756	2.0Gb	2327	2.1Gb	2525	2.6Gb	2103	3.8Gb

		1. NORMALIZATION & selection of 2000 most variable genes										2. ON-DISK REALIZATION of the normalized dataset								3. PCA of the normalized dataset

27998 x 12500	[s]	81	1.1Gb	58	1.3Gb	45	1.6Gb	34	2.0Gb	2000 x 12500	[s]	4	1.1Gb	4	1.3Gb	4	1.5Gb	5	1.8Gb	76	1.1Gb	55	1.3Gb	70	1.4Gb	121	1.7Gb
	[D]	55	1.3Gb	50	1.4Gb	56	1.7Gb	48	2.3Gb		[D]	5	1.3Gb	3	1.3Gb	4	1.6Gb	3	1.8Gb	91	1.4Gb	84	1.5Gb	86	1.6Gb	75	2.0Gb
	[Ds]	54	1.1Gb	49	1.3Gb	45	1.6Gb	38	2.2Gb		[Ds]	9	1.1Gb	6	1.4Gb	6	1.7Gb	6	2.2Gb	198	1.5Gb	183	1.8Gb	233	1.9Gb	224	2.5Gb

27998 x 25000	[s]	162	1.2Gb	127	1.3Gb	95	1.6Gb	75	2.2Gb	2000 x 25000	[s]	8	1.2Gb	9	1.3Gb	9	1.7Gb	8	2.1Gb	128	1.2Gb	127	1.3Gb	127	1.5Gb	155	2.0Gb
	[D]	103	1.3Gb	102	1.5Gb	117	1.6Gb	105	2.1Gb		[D]	13	1.3Gb	9	1.4Gb	7	1.5Gb	6	1.6Gb	217	1.4Gb	186	1.5Gb	184	1.6Gb	150	2.1Gb
	[Ds]	102	1.3Gb	94	1.4Gb	95	1.7Gb	91	2.6Gb		[Ds]	18	1.3Gb	15	1.5Gb	11	1.5Gb	11	2.4Gb	401	1.6Gb	398	1.8Gb	426	1.8Gb	435	3.1Gb

27998 x 50000	[s]	341	1.1Gb	250	1.3Gb	181	1.8Gb	146	2.7Gb	2000 x 50000	[s]	18	1.1Gb	18	1.3Gb	18	1.7Gb	17	2.7Gb	238	1.1Gb	205	1.4Gb	256	1.6Gb	229	2.3Gb
	[D]	194	1.4Gb	204	1.5Gb	226	1.7Gb	214	2.2Gb		[D]	27	1.4Gb	18	1.3Gb	13	1.4Gb	13	1.8Gb	358	1.4Gb	341	1.5Gb	403	1.3Gb	360	2.0Gb
	[Ds]	210	1.2Gb	183	1.5Gb	176	2.0Gb	172	2.5Gb		[Ds]	40	1.2Gb	28	1.5Gb	22	1.9Gb	21	2.6Gb	701	1.7Gb	653	2.0Gb	732	2.3Gb	684	3.4Gb

27998 x 100000	[s]	684	1.2Gb	471	1.4Gb	374	1.9Gb	292	2.9Gb	2000 x 100000	[s]	33	1.2Gb	30	1.4Gb	32	1.9Gb	29	2.7Gb	434	1.2Gb	351	1.4Gb	496	1.7Gb	410	2.5Gb
	[D]	376	1.4Gb	407	1.5Gb	452	1.7Gb	433	2.3Gb		[D]	61	1.4Gb	44	1.4Gb	33	1.6Gb	25	2.0Gb	990	1.5Gb	788	1.6Gb	656	1.9Gb	859	1.4Gb
	[Ds]	405	1.2Gb	359	1.5Gb	341	2.1Gb	349	2.8Gb		[Ds]	88	1.3Gb	68	1.5Gb	51	2.1Gb	44	2.7Gb	1621	1.8Gb	1425	2.0Gb	1591	2.4Gb	1451	3.4Gb

27998 x 200000	[s]	1373	1.3Gb	1033	1.6Gb	774	1.9Gb	643	2.9Gb	2000 x 200000	[s]	62	1.4Gb	62	1.6Gb	66	1.9Gb	69	2.8Gb	787	1.4Gb	692	1.6Gb	929	1.7Gb	944	2.4Gb
	[D]	762	1.4Gb	820	1.9Gb	925	1.7Gb	930	3.3Gb		[D]	138	1.3Gb	102	1.9Gb	75	1.5Gb	57	2.5Gb	2487	1.6Gb	1725	2.1Gb	1361	1.8Gb	1285	2.2Gb
	[Ds]	860	1.3Gb	718	1.8Gb	686	2.3Gb	748	3.4Gb		[Ds]	185	1.4Gb	141	1.8Gb	119	2.1Gb	92	2.9Gb	3745	1.9Gb	3048	2.1Gb	3196	2.9Gb	3187	3.3Gb

Table 1: Timings for DELL XPS 15 laptop Formats: [s] TENxMatrix (sparse) — [D] HDF5Matrix (dense) — [Ds] HDF5Matrix as sparse. For each operation, the best time across the four different block sizes is displayed in bold. In addition, if it’s also the best time across the three formats ([s],[D], and [Ds]), then we box it (only for Normalization and PCA). The “max. mem. used” is the max RSS (Resident Set Size) value obtained by running `ps u -p <PID>` every second while performing a given operation.

5.4.2 Timings for Supermicro SuperServer 1029GQ-TRT

Machine specs:

Supermicro SuperServer 1029GQ-TRT
OS	Linux Ubuntu 22.04	Disk performance	Output of `sudo hdparm -Tt <device>`: `Timing cached reads: 20592 MB in 1.99 seconds = 10361.41 MB/sec Timing buffered disk reads: 1440 MB in 3.00 seconds = 479.66 MB/sec`
RAM	384GB
Disk	1.3TB ATA Disk

Timings:

Test Dataset	F o r m a t	block size = 16 Mb		block size = 40 Mb		block size = 100 Mb		block size = 250 Mb		Normalized Test Dataset	F o r m a t	block size = 16 Mb		block size = 40 Mb		block size = 100 Mb		block size = 250 Mb		block size = 16 Mb		block size = 40 Mb		block size = 100 Mb		block size = 250 Mb
nrow x ncol (# genes x # cells)	F o r m a t	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	nrow x ncol (# sel. genes x # cells)	F o r m a t	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used

		1. NORMALIZATION & selection of 1000 most variable genes										2. ON-DISK REALIZATION of the normalized dataset								3. PCA of the normalized dataset

27998 x 12500	[s]	131	1.1Gb	92	1.3Gb	70	1.5Gb	59	2.0Gb	1000 x 12500	[s]	5	1.1Gb	5	1.2Gb	6	1.4Gb	6	1.6Gb	91	1.1Gb	85	1.2Gb	92	1.4Gb	172	1.5Gb
	[D]	79	1.2Gb	84	1.3Gb	84	1.6Gb	77	2.2Gb		[D]	5	1.2Gb	4	1.2Gb	3	1.0Gb	3	1.6Gb	75	1.3Gb	67	1.3Gb	80	1.5Gb	62	1.6Gb
	[Ds]	73	1.1Gb	72	1.2Gb	59	1.6Gb	59	2.1Gb		[Ds]	9	1.1Gb	7	1.2Gb	7	1.6Gb	8	1.6Gb	250	1.4Gb	217	1.6Gb	274	1.7Gb	252	1.8Gb

27998 x 25000	[s]	250	1.2Gb	193	1.3Gb	156	1.6Gb	115	1.9Gb	1000 x 25000	[s]	11	1.1Gb	12	1.3Gb	11	1.5Gb	12	1.9Gb	180	1.1Gb	147	1.2Gb	178	1.5Gb	168	1.8Gb
	[D]	155	1.2Gb	165	1.4Gb	186	1.6Gb	165	2.1Gb		[D]	11	1.2Gb	6	1.3Gb	7	1.4Gb	7	1.4Gb	170	1.3Gb	125	1.3Gb	148	1.6Gb	221	1.4Gb
	[Ds]	146	1.1Gb	148	1.2Gb	131	1.9Gb	118	2.1Gb		[Ds]	21	1.1Gb	15	1.3Gb	14	1.6Gb	13	2.2Gb	432	1.5Gb	449	1.8Gb	521	1.8Gb	452	2.4Gb

27998 x 50000	[s]	508	1.1Gb	367	1.3Gb	310	2.1Gb	241	2.4Gb	1000 x 50000	[s]	22	1.1Gb	21	1.4Gb	20	1.9Gb	21	2.6Gb	335	1.1Gb	306	1.4Gb	275	1.8Gb	312	2.3Gb
	[D]	292	1.3Gb	341	1.4Gb	358	1.6Gb	329	2.2Gb		[D]	25	1.3Gb	19	1.2Gb	15	1.3Gb	15	1.6Gb	396	1.3Gb	369	1.4Gb	320	1.2Gb	283	2.1Gb
	[Ds]	286	1.1Gb	287	1.3Gb	261	1.8Gb	261	2.6Gb		[Ds]	40	1.1Gb	32	1.4Gb	27	1.6Gb	26	2.2Gb	948	1.5Gb	885	1.9Gb	1138	1.8Gb	939	3.5Gb

27998 x 100000	[s]	1047	1.1Gb	729	1.4Gb	613	1.9Gb	468	3.1Gb	1000 x 100000	[s]	43	1.2Gb	42	1.4Gb	43	1.7Gb	40	2.8Gb	629	1.2Gb	568	1.3Gb	664	1.8Gb	594	2.7Gb
	[D]	595	1.3Gb	662	1.4Gb	734	1.6Gb	739	2.3Gb		[D]	56	1.2Gb	37	1.3Gb	29	1.3Gb	28	1.9Gb	870	1.4Gb	594	1.5Gb	617	1.6Gb	568	2.1Gb
	[Ds]	564	1.2Gb	534	1.5Gb	491	2.1Gb	506	3.3Gb		[Ds]	81	1.2Gb	68	1.4Gb	59	1.7Gb	50	2.9Gb	1972	1.6Gb	1957	2.0Gb	1943	2.3Gb	1739	3.6Gb

27998 x 200000	[s]	2192	1.2Gb	1610	1.4Gb	1252	2.2Gb	1062	3.4Gb	1000 x 200000	[s]	81	1.2Gb	92	1.4Gb	88	2.0Gb	83	3.0Gb	1192	1.2Gb	1098	1.4Gb	1462	1.9Gb	1287	2.9Gb
	[D]	1145	1.3Gb	1354	1.6Gb	1494	1.6Gb	1382	2.9Gb		[D]	119	1.3Gb	96	1.4Gb	78	1.3Gb	57	2.9Gb	2562	1.4Gb	1805	1.7Gb	1581	1.8Gb	1140	2.9Gb
	[Ds]	1175	1.2Gb	1108	1.5Gb	996	2.2Gb	NA	NA		[Ds]	176	1.3Gb	153	1.4Gb	124	2.1Gb	NA	NA	4664	1.8Gb	4202	1.9Gb	4107	2.5Gb	NA	NA

		1. NORMALIZATION & selection of 2000 most variable genes										2. ON-DISK REALIZATION of the normalized dataset								3. PCA of the normalized dataset

27998 x 12500	[s]	134	1.1Gb	93	1.2Gb	78	1.5Gb	59	2.0Gb	2000 x 12500	[s]	7	1.0Gb	7	1.2Gb	9	1.5Gb	8	1.7Gb	135	1.1Gb	127	1.2Gb	161	1.4Gb	123	1.7Gb
	[D]	85	1.2Gb	85	1.4Gb	90	1.7Gb	77	2.2Gb		[D]	7	1.1Gb	5	1.3Gb	5	1.4Gb	5	1.8Gb	146	1.3Gb	133	1.1Gb	144	1.5Gb	126	1.9Gb
	[Ds]	81	1.1Gb	77	1.3Gb	70	1.7Gb	63	2.1Gb		[Ds]	15	1.1Gb	9	1.3Gb	9	1.5Gb	11	2.0Gb	357	1.5Gb	336	1.5Gb	426	1.8Gb	360	2.3Gb

27998 x 25000	[s]	289	1.2Gb	198	1.3Gb	161	1.6Gb	130	2.0Gb	2000 x 25000	[s]	15	1.1Gb	15	1.3Gb	19	1.5Gb	13	2.1Gb	284	1.1Gb	248	1.3Gb	283	1.4Gb	269	1.9Gb
	[D]	164	1.2Gb	174	1.4Gb	195	1.6Gb	180	2.1Gb		[D]	20	1.2Gb	14	1.3Gb	12	1.4Gb	11	1.6Gb	344	1.3Gb	314	1.4Gb	312	1.2Gb	464	1.3Gb
	[Ds]	162	1.1Gb	168	1.3Gb	142	1.8Gb	144	2.3Gb		[Ds]	33	1.1Gb	30	1.3Gb	19	1.6Gb	18	2.2Gb	787	1.4Gb	731	1.7Gb	721	1.9Gb	723	3.3Gb

27998 x 50000	[s]	526	1.1Gb	376	1.3Gb	320	1.7Gb	258	2.3Gb	2000 x 50000	[s]	28	1.1Gb	30	1.3Gb	28	1.6Gb	29	2.2Gb	423	1.1Gb	388	1.3Gb	434	1.5Gb	477	1.9Gb
	[D]	316	1.3Gb	356	1.4Gb	364	1.6Gb	348	2.2Gb		[D]	40	1.2Gb	30	1.2Gb	22	1.4Gb	21	2.0Gb	593	1.3Gb	595	1.4Gb	690	1.1Gb	467	2.1Gb
	[Ds]	324	1.1Gb	302	1.4Gb	277	1.9Gb	285	2.6Gb		[Ds]	64	1.2Gb	50	1.4Gb	38	1.8Gb	39	2.5Gb	1272	1.5Gb	1181	1.9Gb	1350	2.1Gb	1250	3.2Gb

27998 x 100000	[s]	1066	1.2Gb	NA	NA	752	2.1Gb	519	2.8Gb	2000 x 100000	[s]	54	1.1Gb	NA	NA	92	1.8Gb	58	2.7Gb	823	1.2Gb	NA	NA	883	1.7Gb	794	2.4Gb
	[D]	617	1.3Gb	675	1.4Gb	767	1.6Gb	746	2.2Gb		[D]	93	1.2Gb	66	1.3Gb	56	1.5Gb	44	2.1Gb	1549	1.4Gb	1329	1.5Gb	1194	1.7Gb	954	2.2Gb
	[Ds]	646	1.2Gb	581	1.4Gb	532	2.0Gb	566	2.9Gb		[Ds]	135	1.3Gb	104	1.5Gb	85	1.9Gb	77	2.7Gb	2834	1.6Gb	2620	1.9Gb	2880	2.3Gb	2572	3.3Gb

27998 x 200000	[s]	NA	NA	NA	NA	NA	NA	NA	NA	2000 x 200000	[s]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
	[D]	NA	NA	NA	NA	NA	NA	NA	NA		[D]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
	[Ds]	NA	NA	NA	NA	NA	NA	NA	NA		[Ds]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA

Table 2: Timings for Supermicro SuperServer 1029GQ-TRT Formats: [s] TENxMatrix (sparse) — [D] HDF5Matrix (dense) — [Ds] HDF5Matrix as sparse. For each operation, the best time across the four different block sizes is displayed in bold. In addition, if it’s also the best time across the three formats ([s],[D], and [Ds]), then we box it (only for Normalization and PCA). The “max. mem. used” is the max RSS (Resident Set Size) value obtained by running `ps u -p <PID>` every second while performing a given operation.

5.4.3 Timings for Apple Silicon Mac Pro

Machine specs:

Apple Silicon Mac Pro (Apple M2 Ultra)
OS	macOS 13.7.1	Disk performance	N/A
RAM	192GB
Disk	2TB SSD

Timings:

Test Dataset	F o r m a t	block size = 16 Mb		block size = 40 Mb		block size = 100 Mb		block size = 250 Mb		Normalized Test Dataset	F o r m a t	block size = 16 Mb		block size = 40 Mb		block size = 100 Mb		block size = 250 Mb		block size = 16 Mb		block size = 40 Mb		block size = 100 Mb		block size = 250 Mb
nrow x ncol (# genes x # cells)	F o r m a t	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	nrow x ncol (# sel. genes x # cells)	F o r m a t	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used	time in sec.	max. mem. used

		1. NORMALIZATION & selection of 1000 most variable genes										2. ON-DISK REALIZATION of the normalized dataset								3. PCA of the normalized dataset

27998 x 12500	[s]	66	1.4Gb	48	2.2Gb	35	2.5Gb	28	3.2Gb	1000 x 12500	[s]	3	1.6Gb	2	2.2Gb	4	2.5Gb	3	3.4Gb	40	1.7Gb	29	2.0Gb	29	2.7Gb	29	3.4Gb
	[D]	41	1.6Gb	36	2.1Gb	35	2.5Gb	31	3.6Gb		[D]	3	1.7Gb	2	2.1Gb	2	2.5Gb	2	3.6Gb	38	1.9Gb	34	2.2Gb	30	2.8Gb	23	3.7Gb
	[Ds]	36	1.6Gb	35	2.1Gb	30	2.5Gb	29	3.0Gb		[Ds]	5	1.6Gb	3	2.2Gb	4	2.6Gb	4	3.1Gb	104	1.9Gb	101	2.4Gb	125	3.2Gb	104	3.6Gb

27998 x 25000	[s]	138	1.8Gb	102	2.3Gb	76	3.1Gb	54	3.6Gb	1000 x 25000	[s]	5	1.8Gb	5	2.5Gb	6	3.1Gb	6	3.7Gb	61	1.8Gb	60	2.5Gb	76	3.1Gb	62	4.0Gb
	[D]	73	1.9Gb	70	2.4Gb	73	2.4Gb	65	3.8Gb		[D]	6	1.7Gb	3	2.4Gb	4	2.4Gb	3	3.8Gb	74	1.8Gb	53	2.4Gb	54	2.4Gb	50	4.0Gb
	[Ds]	71	1.6Gb	71	2.1Gb	68	3.2Gb	58	4.1Gb		[Ds]	10	1.6Gb	8	2.2Gb	7	3.3Gb	7	4.3Gb	204	2.1Gb	207	2.6Gb	218	3.6Gb	178	4.8Gb

27998 x 50000	[s]	257	1.8Gb	199	2.7Gb	150	3.2Gb	116	5.0Gb	1000 x 50000	[s]	11	1.9Gb	12	2.7Gb	12	3.2Gb	11	4.7Gb	128	2.1Gb	124	2.7Gb	162	3.1Gb	147	4.9Gb
	[D]	144	1.8Gb	140	2.9Gb	144	2.3Gb	139	3.5Gb		[D]	13	1.9Gb	11	2.9Gb	7	2.4Gb	6	3.5Gb	189	2.0Gb	157	2.7Gb	140	2.5Gb	123	3.7Gb
	[Ds]	144	1.7Gb	138	2.2Gb	132	3.2Gb	129	5.5Gb		[Ds]	21	1.9Gb	16	2.2Gb	14	3.4Gb	13	5.9Gb	418	2.6Gb	400	3.2Gb	479	4.1Gb	365	7.1Gb

27998 x 100000	[s]	575	1.9Gb	398	2.9Gb	308	3.5Gb	237	4.9Gb	1000 x 100000	[s]	22	1.9Gb	22	3.0Gb	22	3.5Gb	21	4.6Gb	256	1.9Gb	237	3.0Gb	249	3.5Gb	260	5.1Gb
	[D]	281	2.2Gb	276	2.5Gb	295	3.0Gb	289	3.8Gb		[D]	28	2.0Gb	20	2.5Gb	14	3.0Gb	14	3.9Gb	407	2.1Gb	297	2.6Gb	345	3.0Gb	223	3.9Gb
	[Ds]	285	2.0Gb	269	2.2Gb	259	4.1Gb	255	6.0Gb		[Ds]	43	2.1Gb	33	2.7Gb	27	3.9Gb	25	6.3Gb	887	3.0Gb	935	3.9Gb	910	5.0Gb	730	7.0Gb

27998 x 200000	[s]	1122	2.6Gb	834	3.3Gb	631	4.1Gb	514	7.6Gb	1000 x 200000	[s]	41	2.6Gb	41	3.4Gb	44	4.1Gb	41	7.5Gb	466	2.7Gb	495	3.4Gb	655	4.2Gb	547	7.4Gb
	[D]	561	2.2Gb	554	3.3Gb	588	2.7Gb	595	4.4Gb		[D]	62	2.2Gb	48	3.3Gb	37	2.8Gb	28	4.6Gb	1244	2.5Gb	809	3.4Gb	652	2.9Gb	453	4.7Gb
	[Ds]	583	2.0Gb	529	2.8Gb	501	4.2Gb	525	7.5Gb		[Ds]	91	2.2Gb	78	2.9Gb	62	4.2Gb	52	7.2Gb	2221	3.2Gb	1977	4.2Gb	1992	5.3Gb	1718	7.9Gb

		1. NORMALIZATION & selection of 2000 most variable genes										2. ON-DISK REALIZATION of the normalized dataset								3. PCA of the normalized dataset

27998 x 12500	[s]	68	1.5Gb	50	1.9Gb	39	2.3Gb	29	3.3Gb	2000 x 12500	[s]	4	1.7Gb	4	2.0Gb	3	2.3Gb	4	3.3Gb	58	1.7Gb	54	2.1Gb	48	2.2Gb	52	3.7Gb
	[D]	44	1.8Gb	37	2.1Gb	38	2.5Gb	32	3.8Gb		[D]	4	1.8Gb	3	2.2Gb	2	2.5Gb	2	3.8Gb	72	2.1Gb	55	2.3Gb	56	2.6Gb	51	3.8Gb
	[Ds]	42	1.6Gb	39	2.1Gb	35	2.5Gb	31	3.3Gb		[Ds]	7	1.6Gb	5	2.1Gb	5	2.5Gb	5	3.4Gb	149	2.0Gb	151	2.8Gb	166	3.3Gb	159	4.4Gb

27998 x 25000	[s]	148	1.6Gb	104	2.0Gb	80	2.7Gb	62	4.0Gb	2000 x 25000	[s]	8	1.8Gb	7	2.1Gb	7	2.7Gb	7	4.2Gb	119	1.8Gb	94	2.3Gb	107	2.7Gb	127	4.5Gb
	[D]	80	2.1Gb	72	2.4Gb	76	2.5Gb	71	3.8Gb		[D]	10	2.1Gb	7	2.6Gb	5	2.5Gb	5	3.8Gb	174	2.3Gb	138	2.7Gb	108	2.5Gb	107	4.4Gb
	[Ds]	82	1.6Gb	77	2.1Gb	74	3.1Gb	69	5.1Gb		[Ds]	16	1.7Gb	13	2.3Gb	9	3.1Gb	9	5.2Gb	327	2.5Gb	368	3.0Gb	336	3.6Gb	334	6.2Gb

27998 x 50000	[s]	272	1.8Gb	208	2.2Gb	155	3.4Gb	122	4.8Gb	2000 x 50000	[s]	14	2.0Gb	14	2.4Gb	13	3.1Gb	14	4.8Gb	175	2.2Gb	158	2.3Gb	242	3.1Gb	172	5.0Gb
	[D]	157	2.2Gb	144	3.0Gb	149	2.3Gb	144	3.6Gb		[D]	22	2.2Gb	14	2.9Gb	10	2.3Gb	9	3.7Gb	296	2.3Gb	240	3.0Gb	157	2.3Gb	178	3.7Gb
	[Ds]	163	1.8Gb	147	2.2Gb	141	3.5Gb	137	5.4Gb		[Ds]	33	1.9Gb	25	2.4Gb	20	3.4Gb	17	5.4Gb	563	2.8Gb	561	3.2Gb	591	4.2Gb	526	7.0Gb

27998 x 100000	[s]	596	1.9Gb	414	2.5Gb	318	3.6Gb	247	5.6Gb	2000 x 100000	[s]	30	2.1Gb	27	2.5Gb	28	3.8Gb	28	5.5Gb	355	2.6Gb	299	2.5Gb	385	3.8Gb	341	5.1Gb
	[D]	309	2.3Gb	284	2.6Gb	305	2.9Gb	299	3.7Gb		[D]	52	2.4Gb	35	2.7Gb	27	2.9Gb	18	3.7Gb	771	2.4Gb	621	3.0Gb	525	3.1Gb	676	3.9Gb
	[Ds]	323	1.9Gb	282	2.6Gb	276	4.0Gb	269	6.4Gb		[Ds]	73	1.9Gb	54	2.7Gb	43	4.1Gb	34	6.3Gb	1337	3.0Gb	1241	3.8Gb	1297	5.2Gb	1120	6.8Gb

27998 x 200000	[s]	NA	NA	NA	NA	NA	NA	NA	NA	2000 x 200000	[s]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
	[D]	NA	NA	NA	NA	NA	NA	NA	NA		[D]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
	[Ds]	NA	NA	NA	NA	NA	NA	NA	NA		[Ds]	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA

Table 3: Timings for Apple Silicon Mac Pro Formats: [s] TENxMatrix (sparse) — [D] HDF5Matrix (dense) — [Ds] HDF5Matrix as sparse. For each operation, the best time across the four different block sizes is displayed in bold. In addition, if it’s also the best time across the three formats ([s],[D], and [Ds]), then we box it (only for Normalization and PCA). The “max. mem. used” is the max RSS (Resident Set Size) value obtained by running `ps u -p <PID>` every second while performing a given operation. “max. mem. used” values > 4Gb are displayed in light red.

5.4.4 Timings for Intel Mac Pro

Machine specs:

Intel Mac Pro (24-Core Intel Xeon W)
OS	macOS 12.7.6	Disk performance	N/A
RAM	96GB
Disk	1TB SSD

Timings: Timings for this machine will be added soon.

5.5 Discussion

5.5.1 TENxMatrix (sparse) vs HDF5Matrix (dense)

The “[Ds] HDF5Matrix as sparse” representation didn’t live up to its promise so we leave it alone for now. Note that the code for loading blocks of data from the dense HDF5 file to memory does not currently take full advantage of the SVT_SparseArray representation, a new efficient data structure for multidimensional sparse data implemented in the SparseArray package that overcomes some of the limitations of the dgCMatrix representation from the Matrix package. This will need to be addressed.

5.5.1.1 Normalization

There’s no obvious best choice between the “[s] TENxMatrix (sparse)” and “[D] HDF5Matrix (dense)” representations. More precisely, for normalization, the former tends to give the best times when using bigger blocks (e.g. 250 Mb), whereas the latter tends to give the best times when using smaller blocks (e.g. 16 Mb).

Therefore, on a machine with enough memory to support big block sizes, one will get the best results with the [s] representation and blocks of 250 Mb or more. However, on a machine with not enough memory to support such big blocks, one should instead use the [D] representation with blocks of 16 Mb.

[TODO: Add some plots to help vizualize the above observations.]

5.5.1.2 PCA

For PCA, choosing the “[s] TENxMatrix (sparse)” representation and using small block sizes (40 Mb) tends to give the best performance.

[TODO: Add some plots to help vizualize this observation.]

5.5.1.3 Hybrid approach

Note that, on a machine where using blocks of 250 Mb or more for normalization is not an option, one should use the following hybrid approach:

Start with the “[D] HDF5Matrix (dense)” representation.
Perform normalization, using very small blocks (16 Mb).
Switch to the “[s] TENxMatrix (sparse)” format when writing the normalized dataset to disk, that is, use writeTENxMatrix() instead of writeHDF5Array() for on-disk realization of the intermediate dataset.
Perform PCA on the object returned by the on-disk realization step (writeTENxMatrix()), using small blocks (40 Mb).

5.5.1.4 A note about memory usage

Overall, for a given choice of block size, memory usage doesn’t seem to be affected too much by the number of cells in the dataset, that is, all operations seem to perform at almost constant memory.

However, the “[D] HDF5Matrix (dense)” representation appears to be better than the “[s] TENxMatrix (sparse)” representation at keeping memory usage (mostly) flat as the number of cells in the dataset increases.

5.5.2 Main takeaways

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 this table:

Machine	NORMALIZATION time block size = 250 Mb	REALIZATION time block size = 250 Mb	PCA time block size = 40 Mb	TOTAL time	Max. mem. used
DELL XPS 15 laptop	643	69	692	1404	2.9Gb
Supermicro SuperServer 1029GQ-TRT	NA	NA	NA	NA	NA
Apple Silicon Mac Pro	NA	NA	NA	NA	NA
Comparing times across machines For each machine, we show the normalization, realization, and PCA times (plus total time) obtained on the 27998 x 200000 dataset, using the “[s] TENxMatrix (sparse)” format, and selecting the 2000 most variable genes during the normalization step. All times are in seconds.

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.

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 PCA prefers small blocks.

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

6 Session information

sessionInfo()

## R Under development (unstable) (2025-01-08 r87545)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /home/hpages/R/R-4.5.r87545/lib/libRblas.so 
## LAPACK: /home/hpages/R/R-4.5.r87545/lib/libRlapack.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/Los_Angeles
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] TENxBrainData_1.27.0        SingleCellExperiment_1.29.1
##  [3] SummarizedExperiment_1.37.0 Biobase_2.67.0             
##  [5] GenomicRanges_1.59.1        GenomeInfoDb_1.43.4        
##  [7] RSpectra_0.16-2             DelayedMatrixStats_1.29.1  
##  [9] ExperimentHub_2.15.0        AnnotationHub_3.15.0       
## [11] BiocFileCache_2.15.1        dbplyr_2.5.0               
## [13] HDF5Array_1.35.12           h5mread_0.99.4             
## [15] rhdf5_2.51.2                DelayedArray_0.33.5        
## [17] SparseArray_1.7.5           S4Arrays_1.7.1             
## [19] IRanges_2.41.2              abind_1.4-8                
## [21] S4Vectors_0.45.2            MatrixGenerics_1.19.1      
## [23] matrixStats_1.5.0           BiocGenerics_0.53.6        
## [25] generics_0.1.3              Matrix_1.7-2               
## [27] BiocStyle_2.35.0           
## 
## loaded via a namespace (and not attached):
##  [1] KEGGREST_1.47.0          xfun_0.50                bslib_0.9.0             
##  [4] lattice_0.22-6           rhdf5filters_1.19.0      vctrs_0.6.5             
##  [7] tools_4.5.0              curl_6.2.0               tibble_3.2.1            
## [10] AnnotationDbi_1.69.0     RSQLite_2.3.9            blob_1.2.4              
## [13] pkgconfig_2.0.3          sparseMatrixStats_1.19.0 lifecycle_1.0.4         
## [16] GenomeInfoDbData_1.2.13  compiler_4.5.0           Biostrings_2.75.3       
## [19] htmltools_0.5.8.1        sass_0.4.9               yaml_2.3.10             
## [22] pillar_1.10.1            crayon_1.5.3             jquerylib_0.1.4         
## [25] cachem_1.1.0             mime_0.12                tidyselect_1.2.1        
## [28] digest_0.6.37            purrr_1.0.2              dplyr_1.1.4             
## [31] bookdown_0.42            BiocVersion_3.21.1       fastmap_1.2.0           
## [34] grid_4.5.0               cli_3.6.3                magrittr_2.0.3          
## [37] withr_3.0.2              filelock_1.0.3           UCSC.utils_1.3.1        
## [40] rappdirs_0.3.3           bit64_4.6.0-1            rmarkdown_2.29          
## [43] XVector_0.47.2           httr_1.4.7               bit_4.5.0.1             
## [46] png_0.1-8                memoise_2.0.1            evaluate_1.0.3          
## [49] knitr_1.49               rlang_1.1.5              Rcpp_1.0.14             
## [52] glue_1.8.0               DBI_1.2.3                BiocManager_1.30.25     
## [55] jsonlite_1.8.9           R6_2.5.1                 Rhdf5lib_1.29.0

HDF5Array performance

Compiled 5 February 2025; Modified 5 February 2025

Package

Contents

1 Introduction

2 Install and load the required packages

3 The test datasets

3.1 Sparse vs dense representation

3.2 TENxMatrix vs HDF5Matrix objects

3.2.1 Bring the sparse dataset in R

3.2.2 Bring the dense dataset in R

3.3 Create the test datasets

4 Block-processed normalization and PCA

4.1 Code used for normalization and PCA

4.2 Block processing and block size

4.3 Monitoring memory usage

5 Benchmarks

5.1 Normalization

5.1.1 Normalization of the 27,998 x 12,500 dataset

5.1.1.1 TENxMatrix (sparse)

5.1.1.2 HDF5Matrix (dense)

5.1.1.3 HDF5Matrix as sparse

5.1.2 Normalization of the 27,998 x 25,000 dataset

5.1.2.1 TENxMatrix (sparse)

5.1.2.2 HDF5Matrix (dense)

5.1.2.3 HDF5Matrix as sparse

5.2 On-disk realization of the normalized datasets

5.2.1 On-disk realization of the 1000 x 12,500 normalized dataset

5.2.1.1 TENxMatrix (sparse)

5.2.1.2 HDF5Matrix (dense)

5.2.1.3 HDF5Matrix as sparse

5.2.2 On-disk realization of the 1000 x 25,000 normalized dataset

5.2.2.1 TENxMatrix (sparse)

5.2.2.2 HDF5Matrix (dense)

5.2.2.3 HDF5Matrix as sparse

5.3 PCA

5.3.1 PCA on the 1000 x 12,500 normalized dataset

5.3.1.1 TENxMatrix (sparse)

5.3.1.2 HDF5Matrix (dense)

5.3.1.3 HDF5Matrix as sparse

5.3.2 PCA on the 1000 x 25,000 normalized dataset

5.3.2.1 TENxMatrix (sparse)

5.3.2.2 HDF5Matrix (dense)

5.3.2.3 HDF5Matrix as sparse

5.4 Comprehensive timings obtained on various machines

5.4.1 Timings for DELL XPS 15 laptop

5.4.2 Timings for Supermicro SuperServer 1029GQ-TRT

5.4.3 Timings for Apple Silicon Mac Pro

5.4.4 Timings for Intel Mac Pro

5.5 Discussion

5.5.1 TENxMatrix (sparse) vs HDF5Matrix (dense)

5.5.1.1 Normalization

5.5.1.2 PCA

5.5.1.3 Hybrid approach

5.5.1.4 A note about memory usage

5.5.2 Main takeaways

6 Session information