Usage¶

VCSCArray/VCSRArray are standalone array types built on top of SciPy sparse arrays – see Converting from/to scipy sparse to use them without AnnData at all. The sections below also cover vsparse’s optional AnnData integration.

Converting an AnnData object¶

import anndata as ad
import vsparse

adata = ad.read_h5ad("data.h5ad")

# From adata.X, as a VCSCArray (column-compressed)
v = vsparse.from_anndata(adata)

# From adata.raw.X, as a VCSRArray (row-compressed)
v = vsparse.from_anndata(adata, use_raw=True, format="csr")

# From a specific layer
v = vsparse.from_anndata(adata, layer="counts")

Converting from/to scipy sparse¶

import scipy.sparse as sp
from vsparse import VCSCArray

csc = sp.csc_array(dense_or_sparse_matrix)
v = VCSCArray.from_scipy(csc)

back = v.to_scipy()  # -> scipy.sparse.csc_array
dense = v.toarray()  # -> numpy.ndarray

Supported operations¶

v.T  # transpose (free: returns a VCSRArray sharing buffers)
v * 2.0  # scalar multiplication
v / 2.0  # scalar division
v.log1p()  # elementwise log1p
v @ x  # matrix-vector product, x: 1-D array of length n_cols
x @ v  # vector-matrix product, x: 1-D array of length n_rows
v @ B  # matrix-matrix product, B: 2-D array with B.shape[0] == n_cols
B @ v  # matrix-matrix product, B: 2-D array with B.shape[1] == n_rows
v[:, [1, 3, 5]]  # major-axis (column, for VCSCArray) indexing
v[0:2, 0:2]  # general 2-D indexing (falls back to scipy internally)

Attach a decompressed result back onto an AnnData object as a layer:

vsparse.to_layer(adata, v, key="vcsc_roundtrip")

`VCSCAnnData`: X backed directly by a VCSCArray¶

vsparse.VCSCAnnData is an AnnData subclass whose X (and, separately, raw_X) is a VCSCArray/VCSRArray, not a scipy array:

import vsparse

va = vsparse.VCSCAnnData.from_anndata(adata)  # compresses X and raw.X
va.X  # a VCSCArray
va.raw_X  # a VCSCArray (kept separately from anndata's own `.raw`)

# Persist -- read back with the matching classmethod, not anndata.read_h5ad/read_zarr
va.write_h5ad("compressed.h5ad")
va2 = vsparse.VCSCAnnData.read_h5ad("compressed.h5ad")

va.write_zarr("compressed.zarr")
va3 = vsparse.VCSCAnnData.read_zarr("compressed.zarr")

# Escape hatch: decompress to a normal, fully-featured AnnData
plain = va.to_anndata()

Standard AnnData validates every array assigned to X/layers/etc. against a fixed allowlist of types, so a plain AnnData cannot hold a VCSCArray in X.

On-disk compression and the IVCSC/IVCSR storage format¶

write_h5ad/write_zarr compress every numeric array with zstd at level 1 by default (h5py/hdf5plugin for .h5ad, numcodecs/zarr’s native zstd codec for zarr). On read, .h5ad arrays are decompressed chunk-by-chunk in parallel (numba.get_num_threads() threads) rather than through HDF5’s single-threaded filter pipeline. Files written with the earlier Blosc2+LZ4 default still read, through the ordinary h5py path. Pass dataset_kwargs={} to write uncompressed, or your own dataset_kwargs to use a different codec.

X/raw_X are always stored as IVCSC/IVCSR – the VCSC/VCSR layout, but with the minor-axis indices array delta+varint byte-packed (inspired by IVSparse’s IVCSC). This is a file-storage-only format: it trades extra CPU on write/read for a smaller file, and there is no in-memory IVCSC array type to compute with – read_h5ad/read_zarr always hand back an ordinary VCSCArray/VCSRArray, decompressed from IVCSC/IVCSR immediately on load.

va.write_h5ad("archived.h5ad")
va4 = vsparse.VCSCAnnData.read_h5ad("archived.h5ad")  # va4.X is a plain VCSCArray

VCSCAnnData works around this by overriding the X property; as a consequence, operations that need anndata’s normal per-element type dispatch on X – slicing into views, concatenation, most of the scanpy/anndata ecosystem – are not supported while X is VCSC/VCSR-backed. Call .to_anndata() first if you need those.

VCSCArray/VCSRArray are also registered with anndata’s on-disk IO registry directly, so they round-trip through write_elem/read_elem even when nested inside a plain AnnData’s .uns, independent of VCSCAnnData.

Loading and normalizing IVCSR files¶

vsparse.load_packed decodes an IVCSR/IVCSC-stored .h5ad into a VCSCAnnData. Loading, subsetting, and normalization are separate steps; normalization is a lazy view (.normalized("parafac2")) over the loaded array, so the steps can be mixed and matched rather than fused into one loader:

import vsparse

va = vsparse.load_packed("archived.ivcsr.h5ad")
norm = va.X.normalized("parafac2")

Subsetting while loading¶

load_packed(path, obs_idx=..., var_idx=...) returns exactly load_packed(path)[obs_idx, var_idx], but subsets X as it decodes rather than after. Each selection is a boolean mask over the axis, an array of integer positions (negative positions count from the end; unsorted or repeated positions come back in that order), or a slice; select by label by converting first, e.g. with obs.index.get_indexer(labels).

For an IVCSR file only the selected cells’ indices are decoded, and a gene selection drops and renumbers indices in that same pass (IVCSC: the same, with cells and genes swapped). The full matrix is never allocated. Every packed byte is still read from disk, though: with no stored per-row byte offsets, finding where each selected cell starts takes a parallel scan of the stream. The subset doesn’t carry the file’s statistics sidecar (see below), which describes the whole matrix; its own stats are computed on first use.

import anndata, h5py

with h5py.File("archived.ivcsr.h5ad") as f:  # read only obs to build the mask
    mask = anndata.io.read_elem(f["obs"])["condition"] == "control"
controls = vsparse.load_packed("archived.ivcsr.h5ad", obs_idx=mask, var_idx=hvg_positions)

Statistics sidecar¶

Each VCSCArray/VCSRArray carries a vsparse.VCSStats as arr.stats: per-row row_sums/row_nnz (total counts, n-genes-by-counts), per-column col_sums/col_nnz (total counts, n-cells-by-counts), and col_depth_sums, the per-column sum of depth-normalized values sum_i x[i, j] / row_sums[i]. They are computed on first access (one pass over the data) and kept.

write_h5ad/write_zarr always store them, as a versioned stats subgroup of X/raw_X, and readers (including load_packed) attach them, so QC filters on these columns need no read of X. Normalization takes every recipe’s per-cell and per-gene scales (a/b in y = (g(x * a * b) - c) * s) from the sidecar; the centering/scaling terms c/s depend on the nonlinear g, so they still take one pass per recipe. A file without a usable sidecar – written by an earlier version, or with an unrecognized sidecar version – reads as before, and the statistics are recomputed when first needed.

CUDA¶

A normalized view can be moved onto an NVIDIA GPU with .to_gpu(), which needs CuPy (pip install vsparse[cuda]):

import numpy as np
import vsparse

arr = vsparse.VCSRArray.from_scipy(counts)
gpu = arr.normalized("parafac2").to_gpu()

B = np.random.default_rng(0).normal(size=(arr.shape[1], 16))
out = gpu @ B  # a float32 CuPy array
out = B.T @ gpu  # the other direction, likewise

The transfer keeps the value-compressed layout: to_gpu() uploads the same major_ptr/values/value_ptr/indices arrays the host holds, plus the \(O(n_{\text{rows}} + n_{\text{cols}})\) statistics, never an expanded one-float-per-nonzero copy. That is the point on a device, where memory is the binding constraint: the benchmarks measure the device-resident view at 0.51x to 0.62x the bytes of the equivalent CuPy CSR.

Keeping the layout is also why vsparse ships its own kernels. cuSPARSE reads CSR and CSC only, so nothing in CuPy can walk the two-level major_ptr -> (values, value_ptr) -> indices structure. Both directions of the product, in both formats, are handled by the four kernels in vsparse._cuda.

Precision¶

Everything floating-point on the device is float32, without exception – the values, the statistics, the dense operand and the result. The statistics themselves are still computed on the host in float64 and narrowed on transfer, since the \(O(nnz)\) passes that derive them run once where the matmul runs many times.

So a device product will not match the host’s float64 one to better than roughly 1e-6 relative. In the two misaligned directions (VCSC @ B and B @ VCSR) the kernels also accumulate with atomicAdd, so repeated runs can differ in the last bits; the other two are deterministic.

Materializing instead¶

CudaNormalizedView.to_cupy_sparse() is the device twin of to_scipy_sparse(): it expands the runs into a CuPy CSR/CSC of the uncentered Delta term, leaving means to subtract externally. Useful for handing the matrix to a library that wants real cuSPARSE input, and it is the baseline the kernels are benchmarked against.

Two traps it handles, both worth knowing about if you build such a matrix yourself:

A VCS slice stores its indices grouped by the value they share, never ascending, so the expanded matrix is not canonical. cuSPARSE measured ~6x slower on the unsorted form, which is why to_cupy_sparse sorts by default (pass sort_indices=False to skip it).
cuSPARSE is much faster with a CSR operand in both directions – B @ csc measured ~5x B @ csr – so pass format="csr" if you will multiply repeatedly, even from a VCSC-backed view.