SpaNorm() gains a BPPARAM argument to parallelise normalisation across workers via BiocParallel, accelerating large datasets. It defaults to BiocParallel::SerialParam() (no parallelisation), and results are identical regardless of the backend used.
SpaNorm() now normalises DelayedArray-backed count assays (e.g. disk-backed via HDF5Array) block-wise, so out-of-core datasets are processed without ever loading the full matrix into memory. Results match the in-memory path.
Exported fitNB(), which fits a per-gene negative binomial GLM over an arbitrary design matrix using SpaNorm's IRLS engine (with optional ridge regularisation and adjustable outlier winsorisation). This exposes the model-fitting machinery for reuse independently of SpaNorm's spatial model.
fitNB() and calculateMu() gain an offset argument: a genes x cells matrix added to the linear predictor with its coefficient fixed at 1, so log(mu) = gmean + tcrossprod(alpha, W) + offset. Use it for an effect that is already known rather than adding a column to the design, whose coefficient would float and could absorb signal correlated with that effect. The offset is subset alongside the counts, applied to the dispersion estimation as well as the mean, and sliced in step with the counts when the fit is carved into gene-blocks. It defaults to NULL, which is a strict no-op, so all existing behaviour is unchanged. The GPU path is implemented backend-agnostically but has not yet been executed on an accelerator.
fitNB() also gains a psi argument for supplying per-gene NB dispersions instead of estimating them. Supplied dispersions are used as-is (no re-estimation, no winsorisation) and the outer dispersion loop is bypassed: the coefficients come from a single IRLS fit at the given psi, so edgeR::estimateDisp is never called and the returned sampling factor carries no "dispersion" level. The dispersions need not come from an identical design -- values estimated on a design nested in (or equal to) the fitting design are appropriate, e.g. pooled across a coarser model, which errs conservative. psi composes with offset, and psi = NULL (the default) is a strict no-op.
Exported the quasi-likelihood dispersion machinery of edgeR v4 as generic NB-GLM functions with a CPU and a torch backend: nbUnitDeviance() (the NB unit deviance, elementwise over genes x cells), nbDevianceMoments() (the unit deviance's first two moments under the fitted NB by direct summation over the pmf, the definition edgeR's Chebyshev tables approximate) and qlDispersion() (per-gene adjusted deviance, effective residual df and their ratio). The moments depend only on (mu, phi), so qlDispersion() evaluates them once on a shared (log mu, log phi) table and interpolates onto every gene and cell; the per-gene work is then elementwise and runs on the accelerator when the counts and means are torch tensors. Oracle-tested against edgeR::glmQLFit(); used by spiDE for the standard-error scale of its niche tests.
Exported polishNB(), nbProfilePsi() and the Newton solvers behind them (nbNewtonSolver(), nbNewtonSolverBatch(), nbGramBatch(), nbAbsorbGramBatch(), nbMuFloor()), moved here from spiDE's per-gene convergence stage. polishNB() takes an already-fitted fitNB()-style model (design, coefficients, dispersion) and converges every gene separately to its own penalised negative binomial optimum by damped Newton, gene-blocked and dispatched via BiocParallel; nbProfilePsi() re-estimates just the dispersion by profile maximum likelihood at fixed coefficients. Both take an offset (a known effect added to the linear predictor with its coefficient fixed at 1, rather than a floating design column that could absorb correlated signal); polishNB() also takes absorb/absorb.batch (a wide indicator block -- such as per-sample random-effect columns -- folded into the per-gene solve by a Schur complement, at the cost of one dense-column gram per iteration regardless of block width). spiDE calls these for its own per-gene convergence stage.
Added polishSpaNorm() and isPolished(). polishSpaNorm() maps a fitted SpaNorm() (or SpaNormSVG() null) model onto polishNB()'s generic problem -- the shared library-size coefficient a1 pulled out of the design as a per-cell offset, the per-gene mean gmean an unpenalised intercept column, and the biology/library-size/batch columns ridge-penalised as SpaNorm() penalises them (lambda.a times the number of cells/spots, by column wtype) -- converges every gene, and rewrites the normalised assay from the polished fit. psi.method sets how each gene's dispersion is handled at the converged mean ("fixed", the default, keeps SpaNorm()'s value; "profile" re-estimates it by profile ML and re-polishes); ls sets whether the shared a1 is held fixed ("fixed", the default) or re-estimated jointly with the per-gene coefficients at the optimum of the total penalised likelihood over the polished genes ("joint"); cells chooses all cells/spots or only the ones SpaNorm() sampled to fit; a gene with no counts in the polished cells, or with none in the polished cells of one batch level (a level of the unpenalised 0/1 batch columns), is held out and keeps its input fit (polished = FALSE, the reason in the diagnostics' held_out column, and a warning counting the batch-level cases) rather than driven to an unbounded coefficient; a warning also counts any gene the Newton engine could not polish, and an unpenalised block that is collinear on the polished cells (a batch matrix given with every level's indicator) is caught up front by one QR, so that no gene is polished along its null direction and the outcome does not depend on the platform's rounding; and overwrite allows re-polishing an already-polished fit from the unpolished one kept alongside it. isPolished() reports whether a stored fit has been polished. Measured on four real cores (YTMA CosMx WTA, 948-16,350 genes x 395-2,645 cells, psi.method = "fixed" only): ls = "joint" moved the shared a1 by 7.9-79.0 of its own profiled standard error and took the SVG calls (FDR < 0.05) from 137 to 123 on the 395-cell core, and from 5 to 5, 180 to 178 and 545 to 551 on the other three. On two of the cores (10,738 x 395 and 948 x 646), over a random 50 genes each, it shifted the genes' fitted log-means by a median 0.09 and 0.03 of their own standard error (max 0.68 and 0.17), concentrated in the lowest-library-size cells. The joint polish took 4-7x the wall time of a fixed one in that comparison, which ran with an earlier stop on a1 (tol.ls = 1e-6) under which every joint fit took all 10 steps; the shipped stop (1e-3) ends sooner, and the cost has not been re-measured. ls = "fixed" is the default because SpaNorm()'s own a1 is an unweighted mean over genes where the joint estimate is information-weighted and so dominated by the brightest ones, changing the estimand rather than only its precision. That comparison covers four cores only, with no ground truth for which a1 is correct and no independent validation of the resulting SVG-call changes.
Added nbBlockDesign(), a compact form for a design with a per-group block -- [X | Z_1 | ... | Z_G], each Z_g non-zero only on group g's cells, such as a per-patient intercept plus a per-patient spatial library-size spline -- given as the dense columns X, one cells x q matrix Z of each cell's own block values, and the group of every cell. polishNB(), nbProfilePsi() and nbNewtonSolver() accept it as W and never build the dense n x (p_x + G q) matrix: the linear predictor and the score cost O(n (p_x + q)), the per-gene gram O(n (p_x + q)^2), and every group's q x q block is absorbed by a Schur complement, whatever the number of groups. The result is the dense grouped result (as.matrix() of the design with absorb grouping the block columns by group), tested to 1e-8 on coefficients, log-likelihoods, dispersions, the Schur complement, the Newton step and xcov(). A group whose block is rank-deficient (fewer cells than columns, a level with no cells, an unpenalised column constant over the group; detected on the unit-diagonal scale at rank.tol = 1e-10) is solved by a generalised inverse, which fits exactly as the design with the aliased columns removed. CPU only. Measured on one dedicated core with one BLAS thread (benchmarks/block_design_timing.R, simulated data, 45 patients, 40 dense columns, a 10-column block per patient -- intercept plus a 3 x 3 spatial library-size spline -- 3 genes, profile dispersion, SLURM job 29107530, results in benchmarks/block_design_timing_29107530.csv): against the dense grouped path the per-gene factorisation was 8.0-13.2x faster and the whole per-gene polish 3.8-5.6x faster over 20,000-300,000 cells (13.9 s against 3.7 s per gene at 300,000 cells), with the same coefficients, log-likelihoods, dispersions and Newton iteration counts; and at every size the compact spline model took 0.85-0.97x the per-gene time of the intercept-only model (the patient intercepts alone, absorbed as 1x1 blocks from a dense design), although it took more Newton iterations (10.7-19 against 6-7 per gene). The intercept-only model itself, passed as nbBlockDesign(X, matrix(1, n, 1), patient), gave the same fit in 0.62-0.93x the time of the dense design with a logical absorb. Peak resident memory of the process running all four arms at 300,000 cells was 12.5 GiB, set by the dense grouped arm (its design matrix alone is 1.1 GiB, against 0.26 GiB for the compact design object). At the shape of spiDE's depth = "spatial_spline" design for one index type (benchmarks/block_design_spide_shape.R, built exactly as spiDE's own design functions build it: 30,000 cells, 50 patients with 1-3 sections each, a per-patient block of an intercept plus [l | l * B] per section -- 759 block columns, 6-31 per patient -- and 10 dense columns; 1,000 genes, profile dispersion, one dedicated core, SLURM job 29111836, results in benchmarks/block_design_spide_shape_29111836.csv), the compact design (each patient's block padded with zero columns to the widest, which spiDE's 1e-3 ridge holds at 0) took 0.48 s per gene against 1.57 s on the dense grouped path (3.3x; 8.0 against 26.2 minutes for the 1,000 genes), and the intercept-only model 0.35 s; the two spline fits agreed to 2.3e-13 on every coefficient of every gene.
Exported tpsBasis(), the tensor-product natural-spline basis over spatial coordinates that SpaNorm's fit is built from (identical to the internal basis on its own coordinates). It can be built on reference coordinates (such as a whole section) and evaluated at any subset, with the columns centred by the reference means, so every subset of a section shares one basis; df may be given per axis.
The optional GPU backend now uses the torch package instead of TensorFlow, adding native support for NVIDIA CUDA and Apple Silicon (Metal/MPS) devices and removing the Python/reticulate dependency. Users of backend = "gpu" should install torch in place of tensorflow.
The GPU backend now automatically detects available accelerator memory and fits large datasets in gene-blocks so peak GPU memory stays bounded, avoiding out-of-memory failures on GPUs with limited VRAM. This requires no additional arguments; the detected budget can be overridden via the new gpu.mem.budget parameter, and results match the CPU backend within a small numerical tolerance.
The dispersion winsorisation used during normalisation now clamps at 4 MAD (previously 3), matching the coefficient and mean winsorisation, and is configurable via the winsorisation controls on the fitting/normalisation helpers.
Exported the generic GPU device/memory-budget/tensor-conversion layer (checkGPU(), getBackendDevice(), getBackendDtype(), gpuDtypeBytes(), getGPUMemoryBudget(), resetGPUCache(), is_torch_tensor(), toGPUMatrix(), toGPUVector(), toRMatrix(), diag_mat(), tcrossprod_gpu(), matmul_gpu(), add_vec_mat_gpu(), mult_vec_mat_gpu(), rowSums_gpu(), dnbinom_gpu(), hasBadValues()) and added a new invert_mat_batched(), so downstream packages (e.g. spiDE) can build their own GPU-blocked per-gene computation on the same device/dtype/memory-budget machinery, without depending on unexported internals. calculateMu() now accepts a backend argument ("cpu" by default, unchanged behaviour) and dispatches to the accelerator when requested.
Fixed getGPUMemoryBudget()'s CUDA auto-detection reporting the whole physical GPU's free memory rather than the process's assigned MIG instance's, causing a many-fold budget overestimate (and downstream out-of-memory errors) on MIG-partitioned GPUs; it now resolves the correct device via CUDA_VISIBLE_DEVICES. Added setGPUMemoryBudget() to explicitly set (and cache) the budget for the session, for cases where auto-detection remains unreliable.
SpaNormSVG() now fits, retrieves or re-polishes the null (technical-only) model so that it is paired with the full fit's polish state AND settings (psi.method, ls, cells, and under psi.method = "profile" the dispersion's search interval psi.range), not only whether each is isPolished(): two fits can both be "polished" at different settings, which changes the objective just as much as one being unpolished. svgTest() carries the same check and refuses (rather than silently scoring) a full/null pair that is not paired this way. Either polish the full model before calling SpaNormSVG() (which then polishes the null to match) or call SpaNormSVG() directly and let it pair the null itself.
plotSpatial() and plotCovariate() now find genes with non-syntactic names such as HLA-B. The feature names taken from the aesthetics carried backticks, so no feature matched: plotSpatial() failed with "object 'HLA-B' not found" and plotCovariate() with "no genes selected for the plot".exp(median(log(psi)) + winsor * mad(log(psi)))) is now computed over the finite log-dispersions only. Previously a single NaN dispersion made the threshold NaN, so pmin() turned every gene's logPAC, Pearson and median-bio output into NaN; and a dispersion of exactly 0 (log -Inf) dragged the median to -Inf. Now only the gene with the NaN dispersion is NaN, genes with psi = 0 are excluded from the median/MAD, and if no dispersion is finite the threshold is Inf (no winsorisation). winsorisePsi() and subsetFitGenes() share one internal helper, .psiMax().polishNB() no longer stops when the grouped (absorb with multi-column blocks) solver cannot factor one gene's block -- an unpenalised block with no cells, or with fewer cells than columns. The per-gene engine stopped in sqrt(NULL) and the batch engine on a subscript, taking every gene's polish with it; the gene now falls back to its input fit, flagged singular.SpaNormSVG() being penalised less than the full model. fitSpaNorm() scales the ridge penalty lambda.a by the number of cells/spots before fitting, but fitSpaNormTechnical() rebuilt the penalty from the stored, unscaled lambda.a and did not rescale it, so the nested null was under-penalised by a factor of ncol(Y) on the library-size terms. The two nested fits now use the same penalty, making the likelihood-ratio test consistent. Because the less-penalised null generally fitted slightly better, svg.F was slightly deflated for most genes; SVG statistics will typically increase slightly after refitting the null model, although individual genes can move in either direction (a cached SpaNormNull fit in metadata(spe) is reused as-is, so remove it to refit). fitSpaNormTechnical() now also stops if the data and the full fit have different numbers of cells.plotSpatial to visualise colData, gene expression, and reduced dimensions.plotCovariate to visualise the biolgy, batch, and library size functions estimated by SpaNorm.SpaNormPCA function. The null model is considered to consist of the library size effects, batch effects, and the gene mean.