---
title: "Benchmarking DuckDBArray"
author:
- name: Patrick Aboyoun
  email: aboyounp@gene.com
  affiliation: Genentech, Inc.
date: "Compiled: `r BiocStyle::doc_date()`; Modified: 6 July 2026"
package: DuckDBArray
vignette: >
  %\VignetteIndexEntry{2. Benchmarking DuckDBArray}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output:
  BiocStyle::html_document:
    number_sections: true
    toc: true
    toc_depth: 2
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      error = FALSE, warning = FALSE, message = FALSE)

.matrixstats_ok <- requireNamespace("sparseMatrixStats", quietly = TRUE) ||
    requireNamespace("DelayedMatrixStats", quietly = TRUE)
```

```{r matrixstats-unavailable-note, echo=FALSE, results='asis', eval=!.matrixstats_ok}
cat("> **Note:** Neither `sparseMatrixStats` nor `DelayedMatrixStats` is ",
    "installed in this build environment, so the `rowVars()` comparison below ",
    "is shown but not evaluated. Install one of them (with ",
    "`BiocManager::install()`) to run it.\n", sep = "")
```

# Introduction

The *Introduction to DuckDBArray* vignette shows how a `DuckDBMatrix` behaves
like an ordinary matrix while keeping its data on disk. This vignette asks the
follow-up question: **what does that cost, and what does it buy?** We compare the
DuckDB backend against three alternatives on real single-cell data:

1. an in-memory sparse matrix (`r CRANpkg("Matrix")`'s `dgCMatrix`): the
   fastest option when data fits in RAM, and our baseline;
2. `r Biocpkg("HDF5Array")`: the long-standing Bioconductor backend for
   out-of-memory single-cell analysis;
3. `r Biocpkg("TileDBArray")`: a tile-based on-disk backend.

The headline results below were produced **offline** on the full 10x Genomics
1.3 million brain-cell dataset (see [Benchmark setup]) and are rendered here from
a bundled results file, so this vignette builds quickly. The
[A small, live comparison] section runs a miniature version at build time so you
can see the mechanics end to end.

The short version: for the columnar aggregations that dominate single-cell QC and
feature selection, **DuckDB matches an in-memory sparse matrix while staying on
disk, and clearly beats the other disk-backed backends**, without block-size
tuning.

# What DuckDBArray optimizes

The reason is that `r Biocpkg("DuckDBArray")` turns array summaries into SQL
aggregations that DuckDB executes directly over the Parquet file. Existing
Bioconductor code picks these up automatically: calling `rowSums()` on a
`DuckDBMatrix` *is* a `SUM ... GROUP BY`.

Standard `r Biocpkg("MatrixGenerics")` row/column summaries:

| Function | SQL |
|----------|-----|
| `rowSums` / `colSums` | `SUM()` |
| `rowMeans` / `colMeans` | `AVG()` |
| `rowVars` / `colVars` | `VAR_SAMP()` |
| `rowSds` / `colSds` | `STDDEV_SAMP()` |
| `rowMins` / `rowMaxs` | `MIN()` / `MAX()` |

Plus two helpers for count-based feature selection:

| Function | Description | Use |
|----------|-------------|-----|
| `rowNnzs` / `colNnzs` | count non-zero entries | detection rate, filtering |
| `rowDeviances` | binomial/Poisson deviance | highly-variable-gene selection |

# A small, live comparison

To show the mechanics without a large download, we build a small sparse matrix
and run one summary on both an in-memory `dgCMatrix` and a `DuckDBMatrix`,
confirming the results agree.

```{r live-demo}
library(DuckDBArray)
library(Matrix)
library(MatrixGenerics)

set.seed(1L)
m <- Matrix(rpois(2000 * 400, lambda = 0.3), nrow = 2000, ncol = 400,
            sparse = TRUE)
rownames(m) <- paste0("Gene", seq_len(nrow(m)))
colnames(m) <- paste0("Cell", seq_len(ncol(m)))

path <- tempfile()
writeCoordArray(m, path)
mat <- DuckDBMatrix(path, datacol = "value",
                    keycols = list(index1 = setNames(seq_len(nrow(m)), rownames(m)),
                                   index2 = setNames(seq_len(ncol(m)), colnames(m))),
                    dimtbls = createDimTables(m))

## same answer, one in memory and one queried from disk;
## rowDeviances() is DuckDBArray's own generic and needs no optional package
all.equal(unname(rowDeviances(m, family = "binomial")),
          unname(rowDeviances(mat, family = "binomial")))
```

```{r live-demo-rowvars, eval=.matrixstats_ok}
## rowVars(m) dispatches to sparseMatrixStats/DelayedMatrixStats for the plain
## dgCMatrix; rowVars(mat) is DuckDBArray's own method and needs neither
all.equal(rowVars(m), rowVars(mat))
```

At this size the in-memory matrix is faster, there is nothing to gain from
going to disk. The advantage appears at scale, which is what the offline
benchmark measures.

# Benchmark setup

The full benchmark uses the 10x Genomics 1.3 million brain-cell dataset
(`27,998 x 1,306,127` UMI counts), the same dataset used by the
`r Biocpkg("HDF5Array")` performance vignette and available through
`r Biocpkg("ExperimentHub")` (accession `EH1039`). The headline table uses a
200,000-cell subset, the same size the `r Biocpkg("HDF5Array")` performance
vignette headlines, so the numbers are directly comparable; `inst/scripts/` can
run other cell counts.

## Giving every backend its best effort

The backends parallelize very differently, so a fair comparison configures each at
its best on the **same core budget**:

- **DuckDBMatrix** *autotunes*: it parallelizes inside the SQL engine across all
  available cores with no configuration.
- **HDF5Array** and **TileDBArray** need R-level finesse. Forked workers
  (`MulticoreParam`) clash with HDF5 file locking and TileDB's internal (TBB)
  threading, so the benchmark drives them with **`BiocParallel::SnowParam()`**
  (separate processes, no fork), initializes each worker (a fresh TileDB context;
  HDF5 file locking disabled), **budgets TileDB's internal threads across workers**
  to avoid oversubscription, and tunes the block size, the standard levers from
  the `r Biocpkg("HDF5Array")` and `r Biocpkg("TileDBArray")` playbooks.
- **dgCMatrix** is single-threaded (the in-memory baseline).

To keep this transparent we report **two regimes**: *single-threaded* (one core per
backend, per-core efficiency, matching the `r Biocpkg("HDF5Array")` vignette's
own methodology) and *best effort* (every backend given the full core budget,
configured as above). The exact per-backend configuration is recorded with the
results and shown beneath the tables.

# Results

Rendered from the bundled offline results (`inst/scripts/benchmark_results.rds`);
regenerate with `inst/scripts/run_vignette_benchmarks.R` (see that script's header).

```{r timings, echo=FALSE, results='asis'}
helper <- system.file("scripts", "make_timings_table.R", package = "DuckDBArray")
res <- if (nzchar(helper)) { source(helper); load_vignette_timings() } else NULL
if (is.null(res)) {
    cat("_Precomputed benchmark timings are not available in this build; ",
        "generate them with `inst/scripts/run_vignette_benchmarks.R`._\n", sep = "")
} else {
    cat("\n**Best effort (full core budget)**\n\n")
    cat(make_timings_table("parallel", results = res,
        caption = "Elapsed seconds, each backend using the full core budget."),
        sep = "\n")
    cat("\n\n**Single-threaded (one core)**\n\n")
    cat(make_timings_table("serial", results = res,
        caption = "Elapsed seconds, one core per backend."), sep = "\n")
    cat("\n\n")
    timings_config_note(res)
}
```

# Main takeaways

- **DuckDB's real advantage is zero configuration.** It reaches its performance by
  autotuning; matching it with the disk-backed backends takes `SnowParam` workers,
  per-worker contexts, thread budgeting, file-lock flags, and block-size tuning.
- **Basic statistics match in-memory, on disk.** For `colSums` and `rowVars`,
  DuckDB is on par with an in-memory `dgCMatrix` while never loading the matrix.
- **Feature selection can beat in-memory.** `rowDeviances` computes its aggregation
  entirely in SQL, avoiding the intermediate objects the `dgCMatrix` path builds.
- **Parallelism helps the disk backends unevenly.** With 8 `SnowParam` workers,
  `r Biocpkg("TileDBArray")` speeds up well (`rowVars` ~5x, `rowDeviances` ~4x);
  `r Biocpkg("HDF5Array")`'s gains are smaller and operation-dependent (~7x on
  `rowNnzs`, ~2x on `rowDeviances`, but essentially none on `rowVars`), reflecting
  contention on its single file. DuckDB shows no such unevenness, it scales
  ~11--13x across every operation, automatically. Compare the two regime tables.

# Backend comparison

| Feature | dgCMatrix | HDF5Array | TileDBArray | DuckDBMatrix |
|---------|-----------|-----------|-------------|--------------|
| Storage | Memory | HDF5 file | TileDB array | Parquet file |
| Memory footprint | Full data | Blocks | Blocks | Query results |
| Block tuning | N/A | Required | Required | Automatic |
| Parallelism | single-thread | `SnowParam` over blocks | `SnowParam` + TileDB threads | internal, automatic |
| Read by other languages | R only | R, Python | R, Python | R, Python, Julia, ... |

## When to use each

- **`dgCMatrix`**: data fits comfortably in RAM; the fastest choice for small
  and mid-size matrices.
- **`r Biocpkg("HDF5Array")`**: native 10x/HDF5 inputs; mature and widely used.
- **`r Biocpkg("TileDBArray")`**: when TileDB features (versioning, cloud
  arrays) are wanted.
- **`DuckDBMatrix`**: larger-than-memory data dominated by columnar
  aggregations, when you want near in-memory speed without block tuning, or when
  the Parquet files are shared with non-R tooling.

# Running your own benchmarks

`inst/scripts/` reproduces these numbers on your hardware and extends them to
larger cell counts:

- `run_vignette_benchmarks.R`: the operations timed above; writes
  `benchmark_results.rds`.
- `normalize_and_PCA.R`, `run_comparison.sh`: the fuller normalization/PCA
  sweep, comparable to the `r Biocpkg("HDF5Array")` performance vignette.

See `inst/scripts/README.md` for details. For higher-level scuttle/scran
benchmarks (QC, normalization, variance modelling, marker detection), see the
`r Biocpkg("BiocDuckDB")` package.

# Session information

```{r sessioninfo}
sessionInfo()
```