single-svdlib

Sparse SVD and PCA over sprs matrices. Built for matrices too big to densify — you can mask columns and mean-center without copying or touching the original.

[dependencies]
single-svdlib = "2.0"

Quick start

use single_svdlib::{sprs::TriMatI, SvdMat};

let mut tri = TriMatI::<f64, u32>::new((4, 3));
tri.add_triplet(0, 0, 1.0);
tri.add_triplet(1, 1, 2.0);
tri.add_triplet(2, 2, 3.0);
tri.add_triplet(3, 0, 4.0);
let a: SvdMat<f64> = tri.to_csr::<u64>();

let svd = single_svdlib::svd(&a, 2)?;

// A ≈ u · diag(s) · vt
assert_eq!(svd.u.dim(),  (4, 2));  // left vectors are columns
assert_eq!(svd.vt.dim(), (2, 3));  // right vectors are rows
# Ok::<(), single_svdlib::SvdLibError>(())

Picking a solver

single_svdlib::svd calls irlba.

use single_svdlib::{irlba, randomized};
# use single_svdlib::{sprs::TriMatI, SvdMat};
# let mut t = TriMatI::<f64, u32>::new((40, 20));
# for i in 0..40 { for j in 0..20 { t.add_triplet(i, j, ((i * 7 + j * 3) % 11) as f64); } }
# let a: SvdMat<f64> = t.to_csr::<u64>();

let exact = irlba::svd_seed(&a, 10, 42)?;            // reproducible
let pca = irlba::svd_centered(&a, 10, Some(42))?;    // mean-centered, no densifying
let approx = randomized::svd_seed(&a, 10, 42)?;      // faster, less accurate
let better = randomized::svd_block_krylov(&a, 10, 3, Some(42))?;  // slow spectra
# Ok::<(), single_svdlib::SvdLibError>(())

IrlbaConfig and RandomizedConfig give you the knobs:

use single_svdlib::randomized::{svd_with, Normalizer, RandomizedConfig};
# use single_svdlib::{sprs::TriMatI, SvdMat};
# let mut t = TriMatI::<f64, u32>::new((40, 20));
# for i in 0..40 { for j in 0..20 { t.add_triplet(i, j, ((i * 5 + j) % 7) as f64); } }
# let a: SvdMat<f64> = t.to_csr::<u64>();

let cfg = RandomizedConfig::new(10)
    .oversamples(15)
    .power_iterations(4)
    .normalizer(Normalizer::Tsqr)
    .mean_center(true)
    .seed(42);

// Last argument is an optional progress callback, one call per stage.
let svd = svd_with(&a, &cfg, Some(&|stage: &str| eprintln!("{stage}")))?;
# Ok::<(), single_svdlib::SvdLibError>(())

Memory

SvdMat<T> is CsMatI<T, u32, u64> — 32-bit column indices, 64-bit row pointers. That's 12 bytes per non-zero for f64 instead of the 16 you get with usize everywhere, and 8 instead of 16 for f32. The wider pointer still handles more than u32::MAX non-zeros. Need bigger indices? Name them: SvdMat<f64, u64, u64>.

irlba restarts, so it holds work + 1 right vectors and work left vectors (work = rank + 7 by default) no matter how many restarts it takes. Peak is known up front:

(work + 1) · cols + work · rows   scalars

200 000 × 30 000, 4.9M non-zeros, rank 50, 10 cores:

The matrix is 57.8 MiB, against 76.6 MiB with usize indices. Reproduce with cargo run --release --example scale irlba.

Block Krylov's basis is blocks × (rank + oversamples) wide, so its memory grows with blocks. That's what its accuracy costs.

Big matrices, column subsets, PCA

The thing this crate exists for:

use single_svdlib::{irlba, MaskedCsMat};
# use single_svdlib::{sprs::TriMatI, SvdMat};
# let mut t = TriMatI::<f64, u32>::new((200, 60));
# for i in 0..200 { for j in 0..60 { t.add_triplet(i, j, ((i * 7 + j) % 13) as f64); } }
# let counts: SvdMat<f64> = t.to_csr::<u64>();
# let selected_genes: Vec<usize> = (0..60).step_by(3).collect();

// A view over the selected columns. No copy, `counts` untouched.
let view = MaskedCsMat::with_columns(&counts, &selected_genes);

let pca = irlba::svd_centered(&view, 10, Some(42))?;
// pca.u  — scores,   cells x components
// pca.vt — loadings, components x selected genes
# Ok::<(), single_svdlib::SvdLibError>(())

1 000 000 cells × 30 000 genes, 145M non-zeros, masked to 2238 genes, 50 components (cargo run --release --example pca_at_scale):

Three things to know:

A mask doesn't make products cheaper. Every product still walks all the source's non-zeros and checks each against the mask — only the output gets narrower. With a restrictive mask and hundreds of products, extract once instead. It stays sparse:

# use single_svdlib::{irlba, MaskedCsMat, sprs::TriMatI, SvdMat};
# let mut t = TriMatI::<f64, u32>::new((200, 60));
# for i in 0..200 { for j in 0..60 { t.add_triplet(i, j, ((i * 7 + j) % 13) as f64); } }
# let counts: SvdMat<f64> = t.to_csr::<u64>();
# let selected_genes: Vec<usize> = (0..60).step_by(3).collect();
let sub = MaskedCsMat::with_columns(&counts, &selected_genes).to_sparse();
let pca = irlba::svd_centered(&sub, 10, Some(42))?;
# Ok::<(), single_svdlib::SvdLibError>(())

At 1M × 30k that copy took 146 ms and made the PCA 18% faster (44 s → 36 s).

An unconverged result is an error, not a value. irlba won't hand back triplets that missed tol — too easy to skip the diagnostics and carry a bad decomposition downstream. The error tells you the residual and what to change. If you really want a best effort, call .allow_unconverged() and check converged() / max_residual() yourself.

Raise work on tall matrices. With many rows it's re-orthogonalization that dominates, not the sparse products. rank + 30 was 2.2× faster than the default at the same accuracy. cargo run --release --example tune prints the comparison for your shape.

Result

pub struct SvdRec<T> {
    pub d: usize,                 // triplets returned
    pub u: Array2<T>,             // m × d, left vectors are columns
    pub s: Array1<T>,             // d, descending
    pub vt: Array2<T>,            // d × n, right vectors are rows
    pub total_squared_norm: T,    // ‖A‖²_F of what was decomposed
    pub diagnostics: Diagnostics<T>,
}

A ≈ u · diag(s) · vt, same layout as numpy.linalg.svd. Diagnostics::matvecs counts sparse products (a block product against k columns counts as k), which is comparable across solvers. Diagnostics::detail has per-algorithm stuff like irlba's restart count.

Explained variance

# use single_svdlib::{irlba, sprs::TriMatI, SvdMat};
# let mut t = TriMatI::<f64, u32>::new((200, 60));
# for i in 0..200 { for j in 0..60 { t.add_triplet(i, j, ((i * 7 + j) % 13) as f64); } }
# let a: SvdMat<f64> = t.to_csr::<u64>();
let pca = irlba::svd_centered(&a, 10, Some(42))?;

let scores = pca.scores();                    // m × d, u · diag(s)
let var    = pca.explained_variance();        // sᵢ² / (m − 1), sklearn's convention
let ratio  = pca.explained_variance_ratio();  // fraction of total, sums to ≤ 1

assert!(ratio.iter().sum::<f64>() <= 1.0 + 1e-9);
# Ok::<(), single_svdlib::SvdLibError>(())

The ratio needs a denominator a truncated SVD can't give you — the tail you threw away — so every solver makes one extra pass over the non-zeros to record total_squared_norm. That's nothing next to the hundreds of products a solve does.

If the solver was centering, that norm is of the centered matrix, summed per entry so it stays accurate on columns with a large offset. On a MaskedCsMat it covers only the selected entries, so a masked PCA's ratios are relative to the submatrix you decomposed.

Accuracy

The small reduced factors (IRLBA's B, randomized's R) decide the accuracy of the whole result, so they're factored with a one-sided Jacobi SVD rather than bidiagonal QR. Jacobi is accurate to the condition number after column scaling (Demmel & Veselić). On a 5 × 2 operand with κ ≈ 8·10⁶, nalgebra's Golub–Reinsch reconstructed to 1.3·10⁻⁹; Jacobi hits < 1·10⁻¹⁵. It's also why there's no linear-algebra backend in the dependency tree — sprs, ndarray, rayon, num-traits, rand and thiserror is the lot.

Three test suites:

• tests/cross_algorithm.rs — every solver against a dense LAPACK-grade reference.
• tests/robustness.rs — degenerate and hostile inputs (all-zero, rank-deficient, NaN, ∞, huge dynamic range, empty masks) must give a correct answer or a typed error. Never a panic, a hang, or a quietly wrong result.
• tests/properties.rs — invariants over shapes and distributions proptest picks: agreement with the reference, orthonormality, A·vᵢ = σᵢ·uᵢ, truncation error equal to the spectral tail, storage-order and index-width invariance, reproducibility, and the same set again in f32. Push harder with PROPTEST_CASES=100000 cargo test --release --test properties.

Storage order

A compressed matrix can only be walked along its outer dimension, so the two product directions differ. A · D on CSR writes output row i from sparse row i — threads own disjoint rows, no scratch, no reduction. Aᵀ · D scatters, so threads collide and have to accumulate. transpose_view() doesn't get you out of this; it relabels CSR as a CSC view without changing what's traversable. What it does mean is that a CSC-stored matrix gets the cheap kernel for Aᵀ · D for free. Store CSC if your workload is transpose-heavy.

For the scatter direction there's one accumulator per thread, and if threads × cols × k would blow past DEFAULT_SCRATCH_BUDGET (64 MiB) the dense columns are processed in blocks.

Migrating from 1.x

2.0 is a clean break. Full list in CHANGELOG.md.

Several 1.x results were wrong, not just differently shaped — randomized panicked for every built-in matrix type, mean centering computed the wrong quantity, and seed: None meant seed 0. If you have numbers from 1.x, re-check them. Details in the changelog.

LAS2 is deprecated

lanczos is still there so 2.0 doesn't silently drop the API, but it's #[deprecated], behind an off-by-default feature, and you shouldn't use it. Against a dense LAPACK reference it gets the largest singular value 18–100% wrong on every matrix class tested — including diag(n, n-1, …, 1), where asking for the full rank reports 32 when the answer is 40.

# Only to keep a 1.x caller compiling while it migrates.
single-svdlib = { version = "2.0", features = ["las2"] }

This came from 1.x, not the sprs port — running single-svdlib 1.0.9 from crates.io on the same fixtures reproduces the same wrong values. One cause is fixed: imtqlb hoisted its shift p = d[l] out of the iteration loop, where EISPACK's IMTQL1 assigns it inside (label 120), so every eigenvalue after the first used a stale shift. That's where 1.x's imtqlb had some convergence issues warnings came from. At least one more remains: ritvec reads s[k*js + i] — row k — while imtql2 stores eigenvectors as columns. Transposing roughly halves the error but doesn't close it.

Use irlba. It matches LAPACK to 1e-10 on the same fixtures, diagonal case included, and its memory is bounded. cargo test --release --all-features -- --ignored lanczos::tests::report_accuracy_vs_lapack prints the current error profile.

Licence

BSD, same as the original SVDLIBC. See SVDLIBC-LICENSE.txt.

Credits

• Original SVDLIBC by Doug Rohde
• Rust port of SVDLIBC by Dave Farnham
• Extensions and modifications by Ian F. Diks