Single-cell performance and memory

What changed

The development version keeps selected features sparse while preparing native CUDA PCA. It transposes the sparse matrix and uses cudaverse’s existing sparse PCA interface. The dense PCA input is created on the GPU. Inputs whose dense representation is smaller than 1 MiB retain the existing dense host path to avoid sparse setup overhead. CPU and torch execution also retain dense host preparation.

PCA and sparse arithmetic remain float64. Function arguments, result objects, feature ordering, cell identifiers, duplicate-name handling, exact neighbours and explicit TF32 policy are unchanged. GPU PCA still needs dense workspace; this is not an out-of-core decomposition.

Measured workflows

The October 19, 2026 comparison used one NVIDIA A100-SXM4-40GB host, Linux, R 4.3.3, cudaverse 0.4.2, native CUDA, and OMP_NUM_THREADS=2 / OPENBLAS_NUM_THREADS=2. The host exposed 32 virtual CPUs on AMD EPYC 7713; R used the system reference BLAS/LAPACK. The baseline was the unchanged submitted cudacellr 0.4.0 archive; the candidate runtime was commit b1b2262. Each workflow had five warmups and ten synchronized timed runs in fresh processes for each variant. All cases used k = 15, batch_size = 256, and scaled PCA.

Input	Cells	Genes	Selected genes / PCs	Median seconds, baseline → candidate	Peak host GiB, baseline → candidate
Historical small synthetic	1,000	50	40 / 30	0.032 → 0.030	0.429 → 0.427
Historical medium synthetic	10,000	100	80 / 30	0.176 → 0.166	0.459 → 0.456
Historical large synthetic	50,000	128	102 / 30	1.394 → 1.164	0.698 → 0.611
Sparse stress fixture	50,000	10,000	2,000 / 50	9.435 → 3.367	3.473 → 0.663
10x PBMC v3	11,769	33,538	2,000 / 50	5.465 → 4.567	2.096 → 1.570

The sparse stress fixture has 2,549,762 nonzeros; PBMC has 24,825,783. PBMC comes from 10x Genomics, 10k PBMCs from a Healthy Donor (v3 chemistry), licensed CC BY 4.0. The input manifest pins the download and derived inputs by SHA-256 and records synthetic seeds and all workflow parameters.

The host was available for shared use, without an exclusive GPU reservation. The 200 before/after compute-process observations in the accepted comparison reported no other compute process. These observations do not rule out activity between samples or other host contention. The table reports this run’s observations; it does not promise the same speedup on other hardware.

Peak host memory is the whole R process high-water mark, including runtime, input and returned results. Timings include the public workflow and a final CUDA synchronization, with garbage collection, serialization and diagnostic queries outside the timed call. The initial all-sparse candidate failed the small-input runtime gate; its evidence is retained. The 1 MiB dispatch threshold resolved that regression. The unchanged baseline measurements were reused for the final candidate on the same host, R/core versions and pinned inputs.

All five cases passed exact normalization, HVG, identifier and neighbour-index comparisons, PCA reconstruction/subspace tolerances, distance tolerances and native allocation cleanup. The large sparse fixture exceeded the required 25% host-memory reduction (80.9%), and every case passed the maximum 10% end-to-end slowdown gate. PBMC host memory fell 25.1%.

Where the time goes

Separate single-pass stage profiles locate work; they are not warmed medians and should not be added to the warmed workflow timings above.

Candidate stage, seconds	Sparse stress fixture	PBMC
Normalization	0.473	1.811
Variable-feature selection	0.071	0.785
PCA preparation	0.069	0.207
PCA including upload	2.527	2.022
Exact neighbours	0.856	0.097

The next profiling target is dense PCA workspace and decomposition for the large sparse case. For PBMC, CPU normalization and feature selection together also merit attention. These observations do not justify changing numerical algorithms or precision without a separate correctness study.

Transfers and GPU memory

Separate Nsight Systems 2022.4.2 traces measured CUDA memory operations over the complete profile process, including capability probes and result copies. Profiler runs are excluded from the timing gate. MB below means 1,000,000 bytes.

Trace	Host → GPU MB, baseline → candidate	GPU → host MB, both	GPU → GPU MB, both
Sparse stress fixture	800.016 → 8.176	29.865	800.000
PBMC	188.320 → 9.930	7.691	188.304

The additional native sparse row-index upload is included in these traces. The benchmark JSON’s smaller CSR payload estimate excludes that extra copy and is not a measured transfer total. Native PCA retains GPU scores, and kNN provenance confirms reuse of those scores without another upload.

Peak allocations tracked by cudaverse were 4.646 GiB for the candidate sparse stress workflow and 1.192 GiB for PBMC. The corresponding baseline peaks were 4.638 and 1.183 GiB: reducing host memory does not reduce dense PCA’s GPU workspace. These are package-owned allocation peaks, excluding CUDA context, library-internal allocations and other applications. Allow additional headroom.

Sizing and batch guidance

The table defines tested configurations, not maximum supported sizes. In particular, 50,000 cells with 2,000 selected genes was validated at about 0.5% input density; a denser 50,000-cell matrix can need substantially more host memory. The full sparse input and normalized output still reside in host memory.

One float64 selected dense matrix needs approximately 8 * cells * n_hvg bytes, before decomposition workspace. At 50,000 cells and 2,000 selected genes that is 0.745 GiB per copy. Lower n_hvg when dense GPU memory is the constraint. Reducing n_components alone does not shrink that input matrix.

Exact kNN’s distance block needs about 8 * min(batch_size, cells) * cells bytes, plus scores, outputs and workspace. For 50,000 cells, batch sizes 128, 256 and 512 imply approximately 48.8, 97.7 and 195.3 MiB for that block. Start with the measured value of 256; reduce it to 128 if available memory is tight. The smaller batch’s throughput is not measured here. Batch size preserves exact neighbour definitions and original row-index tie breaking. It does not change quadratic search time or PCA memory.

Remove unused results and call gc() before checking available memory with cudaverse::cuda_memory_info("cuda"). Keep GPU scores with their PCA model when continuing directly into cuda_cell_neighbors().

Reproduce and inspect

The repository’s benchmark evidence contains the input manifest, all timing samples, stage profiles, numerical gate results and transfer summaries. Maintainer scripts under tools/ describe the protocol and preserve failed attempts. Public derivatives remove unrelated process identities; original reports, result matrices and raw CUDA traces are retained with their hashes in the publication/development evidence.