Merge perf/sparse-joint (#52)

745 ms -> 1.11 ms on the fixture #52 was opened about.

The joint precision matrix is 0.19% dense at scale and gets sparser as
the history grows. We allocated all n^2 entries — 31 MB at n = 1976,
128 MB at ustat's ~4000 appearances — filled 99.8% of it with zeros, and
ran an O(n^3) factorisation over the whole thing.

Two measurements shaped the fix, and the first killed the plan #52
proposed.

**Ordering alone does nothing to a dense factorisation.** Its inner
loops run over every k whether the entry is zero or not, so a
permutation changes which entries are zero and not how many
multiplications happen. A 700x700 banded matrix at 0.43% density:
30.196 ms in band order, 29.544 ms under a scramble that destroyed the
band. Identical, as the flop count says it must be. #52's step 1 —
"reorder with AMD, keep our own Cholesky, and measure" — could not have
worked, and measuring said so before any of it was written.

**Sparsity and AMD together are worth four orders of magnitude.**
Symbolic factorisation on the n = 1976 fixture, against 2.572e9 dense
flops: sparse in natural order needs 5.597e7 (46x), sparse under AMD
needs 8.656e4 — 29,710x. AMD is worth 646x on top of sparsity and
nothing without it. Natural order fills in badly for exactly the reason
#52 predicted about bandwidth: nnz(L) is 292,437 against A's 7,504,
because a competitor idle from slice 0 to slice 75 links across the
whole matrix.

Measured end to end, factorising through `History::joint`:

    n =  480     215 us   (bench: 9.11 ms -> 167 us, 54x)
    n = 1976    1.112 ms  (was ~745 ms, 670x)
    n = 7800    4.616 ms  (dense would be 1.58e11 flops)

Scaling is near-linear now rather than cubic: 16x the variables costs
21x the time, where dense would cost 4096x.

The factorisation is the up-looking sparse Cholesky of Davis's *Direct
Methods for Sparse Linear Systems*, written here rather than taken from
a crate. The scouting in #52 still holds and got one addition: `feral`
itself pulls `pulp`, so it has the same runtime CPU-dispatch problem
that ruled out `faer` — results could differ between an AVX-512 host and
an AVX2 one, the drift the libm-over-std decision was made to avoid.
`sprs-ldl` is still LGPL and `nalgebra-sparse` still disclaims
fill-reduction in its own docs. Only the ordering is a dependency:
`feral-amd`, two crates, both `#![forbid(unsafe_code)]`.

The matrix is accumulated into a `BTreeMap`, not a hash map: the
iteration order becomes the summation order, and a hash map's varies per
process. `tests/cross_process_determinism.rs` exists because that has
bitten before.

`whiten` returns its result in the permuted order and leaves it there —
a dot product does not care, as long as both operands were permuted the
same way — so `bilinear` is unchanged.

Correctness: the existing analytic goldens are 2x2 and 3x3, too small to
permute or fill in, so they could not have caught a symbolic-pass bug.
`agrees_with_a_dense_reference_on_random_sparse_systems` checks every
bilinear form against a deliberately naive dense factorisation that
shares no code with the thing it is checking, on chain-plus-long-range
matrices up to n = 60.

Closes #52.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ

Cargo.toml

@@ -47,12 +47,15 @@ harness = false
 
 [dependencies]
 approx = { version = "0.5.1", optional = true }
+feral-amd = "0.2"
 libm = "0.2.16"
 rayon = { version = "1", optional = true }
 smallvec = "1"
 
 [features]
 approx = ["dep:approx"]
+# Exposes the joint sparsity pattern for the #52 measurement. Test-only.
+measure-sparsity = []
 rayon = ["dep:rayon"]
 
 [dev-dependencies]

src/history.rs

@@ -386,8 +386,8 @@ pub(crate) struct CompetitorConfig {
 /// The joint precision over a history's appearances, with the maps needed to
 /// address a competitor either at their latest appearance or at a given slice.
 struct TimeExpanded {
-    /// Row-major precision matrix over appearances.
-    lambda: Vec<f64>,
+    /// Precision matrix over appearances, accumulated sparsely.
+    lambda: crate::joint::SymmetricBuilder,
     /// `(row, slice)` of each competitor's latest appearance.
     latest: HashMap<Index, (usize, usize)>,
     /// Row of each `(competitor, slice)` appearance.
@@ -1415,6 +1415,23 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<K, T, D
         prior_variance * SQRT_EPSILON
     }
 
+    /// Test-only: the joint's sparsity pattern, for the #52 measurement.
+    ///
+    /// Returns `(n, adjacency)` where `adjacency[i]` holds the off-diagonal
+    /// nonzero columns of row `i`.
+    #[cfg(feature = "measure-sparsity")]
+    pub fn joint_pattern_for_measurement(&self) -> (usize, Vec<std::collections::HashSet<usize>>) {
+        let te = self.time_expanded_joint();
+        let n = te.width;
+        let mut adj = vec![std::collections::HashSet::new(); n];
+        for (i, j) in te.lambda.pattern() {
+            if i != j {
+                adj[i].insert(j);
+            }
+        }
+        (n, adj)
+    }
+
     fn time_expanded_joint(&self) -> TimeExpanded {
         let mut latest: HashMap<Index, (usize, usize)> = HashMap::new();
         let mut at_slice: HashMap<(Index, usize), usize> = HashMap::new();
@@ -1455,17 +1472,23 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<K, T, D
             }
         }
 
-        let mut lambda = vec![0.0; n * n];
+        // Accumulated sparsely: this matrix is ~0.19% dense at scale, and the
+        // dense form was 31 MB at n = 1976 and 128 MB at ustat's ~4000
+        // appearances, of which 99.8% was zeros. See `joint.rs` and #52.
+        let mut lambda = crate::joint::SymmetricBuilder::new();
 
         for (row, competitor) in first_rows {
-            lambda[row * n + row] +=
-                1.0 / self.competitors[competitor].rating.prior.sigma().powi(2);
+            lambda.add(
+                row,
+                row,
+                1.0 / self.competitors[competitor].rating.prior.sigma().powi(2),
+            );
         }
         for (a, b, drift) in drift_links {
-            lambda[a * n + a] += 1.0 / drift;
-            lambda[b * n + b] += 1.0 / drift;
-            lambda[a * n + b] -= 1.0 / drift;
-            lambda[b * n + a] -= 1.0 / drift;
+            lambda.add(a, a, 1.0 / drift);
+            lambda.add(b, b, 1.0 / drift);
+            lambda.add(a, b, -1.0 / drift);
+            lambda.add(b, a, -1.0 / drift);
         }
         for (slice_idx, slice) in self.time_slices.iter().enumerate() {
             for (contrast, noise) in slice.scored_contrasts(&self.competitors) {
@@ -1473,7 +1496,7 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<K, T, D
                     let ra = at_slice[&(*ia, slice_idx)];
                     for (ib, cb) in &contrast {
                         let rb = at_slice[&(*ib, slice_idx)];
-                        lambda[ra * n + rb] += ca * cb / noise;
+                        lambda.add(ra, rb, ca * cb / noise);
                     }
                 }
             }

src/joint.rs

@@ -1,4 +1,4 @@
-//! Cholesky factorisation of a joint precision matrix.
+//! Sparse Cholesky factorisation of a joint precision matrix.
 //!
 //! Every question the joint answers is a *bilinear form* in the precision
 //! matrix's inverse — the variance of a contrast is `c^T L^-1 c`, and the
@@ -18,72 +18,324 @@
 //! negative number, where the same quantity as `|L^-1 c|^2` is a sum of
 //! squares and cannot.
 //!
-//! Factorising is `O(n^3)` and whitening is `O(n^2)`, so the split also
-//! matters structurally: the expensive half depends only on the fit, and is
-//! shared across every query a [`Joint`](crate::Joint) answers.
+//! # Why this is sparse (#52)
+//!
+//! A time-expanded joint is *extremely* sparse and gets sparser as the history
+//! grows: a row couples only to its own previous and next appearance through
+//! the drift link, and to whoever co-appeared in its slice. Measured on a
+//! 76-slice, 988-duel, 200-competitor history: `n = 1976`, `nnz = 7504`,
+//! **0.19% dense**.
+//!
+//! This used to store all `n^2` entries and run a dense `O(n^3)` factorisation
+//! over them. Two measurements decided the replacement:
+//!
+//! - **Ordering alone does nothing to a dense factorisation.** Its inner loops
+//!   run over every `k` whether or not the entry is zero. A 700x700 banded
+//!   matrix at 0.43% density factorised in 30.196 ms in band order and
+//!   29.544 ms under a scramble that destroyed the band — identical, as the
+//!   flop count says it must be. Fill-reducing order is worth nothing until
+//!   the factorisation skips zeros.
+//! - **Together they are worth four orders of magnitude.** On that `n = 1976`
+//!   fixture, against `n^3/3 = 2.572e9` flops dense: sparse in the natural
+//!   order needs `5.597e7` (46x better), and sparse under an AMD fill-reducing
+//!   order needs `8.656e4` — **29,710x**. AMD is worth 646x *on top of*
+//!   sparsity and nothing without it.
+//!
+//! Natural ordering fills in badly here for the reason #52 predicted: a
+//! competitor who appears in slice 0 and not again until slice 75 creates a
+//! drift link spanning nearly the whole matrix. `nnz(L)` is 292,437 under the
+//! natural order against 11,583 under AMD, from an `A` with 7,504.
+//!
+//! The ordering comes from `feral-amd`. The factorisation is the up-looking
+//! sparse Cholesky of Davis's *Direct Methods for Sparse Linear Systems*,
+//! written here rather than taken from a crate: the sparse solvers on
+//! crates.io either pull SIMD dispatch (`faer`, and `feral` itself, both
+//! through `pulp`), which would make results differ between an AVX-512 host
+//! and an AVX2 one — the same class of drift the `libm`-over-`std` decision
+//! was made to avoid — or are LGPL, or disclaim fill-reduction in their own
+//! docs.
+
+use std::collections::BTreeMap;
+
+/// A symmetric matrix accumulated entry by entry, before factorisation.
+///
+/// A `BTreeMap` rather than a hash map because the iteration order becomes the
+/// factorisation's summation order, and a hash map's order varies per process.
+/// `tests/cross_process_determinism.rs` exists because that has bitten before.
+#[derive(Default)]
+pub(crate) struct SymmetricBuilder {
+    entries: BTreeMap<(usize, usize), f64>,
+}
+
+impl SymmetricBuilder {
+    pub(crate) fn new() -> Self {
+        Self::default()
+    }
+
+    /// Add `value` to entry `(row, col)`. Both triangles must be supplied.
+    pub(crate) fn add(&mut self, row: usize, col: usize, value: f64) {
+        *self.entries.entry((row, col)).or_insert(0.0) += value;
+    }
+
+    /// The `(row, col)` positions that hold a nonzero. For the #52 measurement.
+    #[cfg(feature = "measure-sparsity")]
+    pub(crate) fn pattern(&self) -> impl Iterator<Item = (usize, usize)> + '_ {
+        self.entries
+            .iter()
+            .filter(|(_, v)| **v != 0.0)
+            .map(|(&rc, _)| rc)
+    }
+}
 
 /// A factorised symmetric positive-definite matrix, reusable across queries.
 pub(crate) struct Cholesky {
-    /// Lower triangle of `L`, row-major `n * n`. The upper triangle is
-    /// leftover scratch from the factorisation and is never read.
-    l: Vec<f64>,
     n: usize,
+    /// `inv[old] = new`: where each original row sits after the AMD reorder.
+    inv: Vec<usize>,
+    /// `L` in compressed-column form, permuted. Within a column the diagonal
+    /// is first and the rest ascend by row.
+    col_ptr: Vec<usize>,
+    row_idx: Vec<usize>,
+    val: Vec<f64>,
 }
 
 impl Cholesky {
-    /// Factorise `a` (row-major, `n * n`, symmetric) into `L L^T`.
-    ///
-    /// `a` is consumed as scratch.
+    /// Factorise the accumulated matrix into `L L^T`, under a fill-reducing
+    /// permutation.
     ///
     /// Returns `None` if the matrix is not positive-definite, which for a
     /// precision matrix means the model is improper — a competitor with
-    /// neither a proper prior nor any evidence.
-    pub(crate) fn factor(mut a: Vec<f64>, n: usize) -> Option<Self> {
-        debug_assert_eq!(a.len(), n * n);
+    /// neither a proper prior nor any evidence — or if the ordering fails.
+    pub(crate) fn factor(built: SymmetricBuilder, n: usize) -> Option<Self> {
+        if n == 0 {
+            return Some(Self {
+                n: 0,
+                inv: Vec::new(),
+                col_ptr: vec![0],
+                row_idx: Vec::new(),
+                val: Vec::new(),
+            });
+        }
 
-        for j in 0..n {
-            let mut d = a[j * n + j];
-            for k in 0..j {
-                d -= a[j * n + k] * a[j * n + k];
+        let inv = Self::amd_permutation(n, &built)?;
+
+        // Upper triangle of the permuted matrix, column-major: column `c`
+        // holds the rows `r <= c`. Exactly one of a symmetric pair survives
+        // the `r <= c` filter, so nothing is double-counted.
+        let mut cols: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
+        for (&(old_r, old_c), &v) in &built.entries {
+            if v == 0.0 {
+                continue;
             }
+            let (r, c) = (inv[old_r], inv[old_c]);
+            if r <= c {
+                cols[c].push((r, v));
+            }
+        }
+        let mut a_ptr = Vec::with_capacity(n + 1);
+        let mut a_row = Vec::new();
+        let mut a_val = Vec::new();
+        a_ptr.push(0usize);
+        for col in &mut cols {
+            col.sort_unstable_by_key(|&(r, _)| r);
+            for &(r, v) in col.iter() {
+                a_row.push(r);
+                a_val.push(v);
+            }
+            a_ptr.push(a_row.len());
+        }
+
+        let parent = Self::etree(n, &a_ptr, &a_row);
+
+        // Symbolic pass: how many entries each column of L will hold. Running
+        // `ereach` per column costs O(nnz(L)) in total, which is the same order
+        // as the numeric pass it sizes.
+        let mut counts = vec![0usize; n];
+        let mut stack = vec![0usize; n];
+        let mut mark = vec![false; n];
+        for k in 0..n {
+            let top = Self::ereach(k, &a_ptr, &a_row, &parent, &mut stack, &mut mark);
+            for &i in &stack[top..] {
+                counts[i] += 1;
+            }
+            counts[k] += 1; // the diagonal
+        }
+
+        let mut col_ptr = Vec::with_capacity(n + 1);
+        col_ptr.push(0usize);
+        for &c in &counts {
+            col_ptr.push(col_ptr[col_ptr.len() - 1] + c);
+        }
+        let nnz = col_ptr[n];
+        let mut row_idx = vec![0usize; nnz];
+        let mut val = vec![0.0f64; nnz];
+
+        // `next[i]` is the slot column `i` will fill next. Column `i`'s
+        // diagonal lands first, at `col_ptr[i]`, because nothing is written to
+        // a column before its own iteration.
+        let mut next: Vec<usize> = col_ptr[..n].to_vec();
+        let mut x = vec![0.0f64; n];
+
+        for k in 0..n {
+            let top = Self::ereach(k, &a_ptr, &a_row, &parent, &mut stack, &mut mark);
+
+            for p in a_ptr[k]..a_ptr[k + 1] {
+                if a_row[p] <= k {
+                    x[a_row[p]] = a_val[p];
+                }
+            }
+
+            let mut d = x[k];
+            x[k] = 0.0;
+
+            for &i in &stack[top..] {
+                let lki = x[i] / val[col_ptr[i]];
+                x[i] = 0.0;
+                for p in col_ptr[i] + 1..next[i] {
+                    x[row_idx[p]] -= val[p] * lki;
+                }
+                d -= lki * lki;
+                let p = next[i];
+                next[i] += 1;
+                row_idx[p] = k;
+                val[p] = lki;
+            }
+
             // Explicit rather than `!(d > 0.0)`: a NaN pivot must fail here
             // too, and a negated comparison would let it through as "not
             // positive".
             if d.is_nan() || d <= 0.0 {
                 return None;
             }
-            let d = d.sqrt();
-            a[j * n + j] = d;
-
-            for i in j + 1..n {
-                let mut s = a[i * n + j];
-                for k in 0..j {
-                    s -= a[i * n + k] * a[j * n + k];
-                }
-                a[i * n + j] = s / d;
-            }
+            let p = next[k];
+            next[k] += 1;
+            row_idx[p] = k;
+            val[p] = d.sqrt();
         }
 
-        Some(Self { l: a, n })
+        Some(Self {
+            n,
+            inv,
+            col_ptr,
+            row_idx,
+            val,
+        })
     }
 
-    /// Whiten a contrast: `y = L^-1 b`.
+    /// AMD fill-reducing order, as `inv[old] = new`.
+    fn amd_permutation(n: usize, built: &SymmetricBuilder) -> Option<Vec<usize>> {
+        let mut cols: Vec<Vec<i32>> = vec![Vec::new(); n];
+        for (&(r, c), &v) in &built.entries {
+            if v != 0.0 {
+                cols[c].push(i32::try_from(r).ok()?);
+            }
+        }
+        let mut col_ptr = Vec::with_capacity(n + 1);
+        let mut row_idx = Vec::new();
+        col_ptr.push(0i32);
+        for (j, col) in cols.iter_mut().enumerate() {
+            col.push(i32::try_from(j).ok()?);
+            col.sort_unstable();
+            col.dedup();
+            row_idx.extend_from_slice(col);
+            col_ptr.push(i32::try_from(row_idx.len()).ok()?);
+        }
+
+        let pattern = feral_amd::CscPattern::new(n, &col_ptr, &row_idx)?;
+        // `perm[new] = old`; we want the inverse.
+        let perm = feral_amd::amd_order(&pattern).ok()?;
+        let mut inv = vec![0usize; n];
+        for (new, &old) in perm.iter().enumerate() {
+            inv[usize::try_from(old).ok()?] = new;
+        }
+        Some(inv)
+    }
+
+    /// Elimination tree of the upper-triangular pattern. `usize::MAX` is "no
+    /// parent", i.e. a root.
+    fn etree(n: usize, col_ptr: &[usize], row_idx: &[usize]) -> Vec<usize> {
+        let mut parent = vec![usize::MAX; n];
+        let mut ancestor = vec![usize::MAX; n];
+        for k in 0..n {
+            for &row in &row_idx[col_ptr[k]..col_ptr[k + 1]] {
+                let mut i = row;
+                while i != usize::MAX && i < k {
+                    let next = ancestor[i];
+                    ancestor[i] = k;
+                    if next == usize::MAX {
+                        parent[i] = k;
+                    }
+                    i = next;
+                }
+            }
+        }
+        parent
+    }
+
+    /// Nonzero pattern of row `k` of `L`, written into `stack[top..n]` in
+    /// topological order. Returns `top`.
+    ///
+    /// `stack` is used from both ends — a scratch region from `0` while walking
+    /// each path up the tree, and the result from `n` downwards. They cannot
+    /// collide because every node is pushed at most once across the whole call.
+    fn ereach(
+        k: usize,
+        col_ptr: &[usize],
+        row_idx: &[usize],
+        parent: &[usize],
+        stack: &mut [usize],
+        mark: &mut [bool],
+    ) -> usize {
+        let n = mark.len();
+        let mut top = n;
+        mark[k] = true;
+        for &row in &row_idx[col_ptr[k]..col_ptr[k + 1]] {
+            let mut i = row;
+            if i > k {
+                continue;
+            }
+            let mut len = 0usize;
+            while i != usize::MAX && !mark[i] {
+                stack[len] = i;
+                len += 1;
+                mark[i] = true;
+                i = parent[i];
+            }
+            // Reverse the path onto the output end, so the result stays in
+            // topological order overall.
+            while len > 0 {
+                len -= 1;
+                top -= 1;
+                stack[top] = stack[len];
+            }
+        }
+        for &i in &stack[top..] {
+            mark[i] = false;
+        }
+        mark[k] = false;
+        top
+    }
+
+    /// Whiten a contrast: `y = L^-1 P b`.
     ///
     /// The point of the result is the dot product, not the vector: for two
     /// contrasts `b` and `b'`, `y . y'` is `b^T A^-1 b'`. See the module docs.
+    ///
+    /// The result is in the permuted order, and stays there — a dot product
+    /// does not care, as long as both operands were permuted the same way.
     pub(crate) fn whiten(&self, b: &[f64]) -> Vec<f64> {
         debug_assert_eq!(b.len(), self.n);
         let n = self.n;
-        let mut y = b.to_vec();
-        for i in 0..n {
-            // Folded from `y[i]` rather than summed and subtracted once, so the
-            // accumulation order matches a plain substitution loop exactly.
-            let row = &self.l[i * n..i * n + i];
-            let s = row
-                .iter()
-                .zip(&y[..i])
-                .fold(y[i], |acc, (l, v)| acc - l * v);
-            y[i] = s / self.l[i * n + i];
+        let mut y = vec![0.0f64; n];
+        for (old, &v) in b.iter().enumerate() {
+            y[self.inv[old]] = v;
+        }
+        for j in 0..n {
+            y[j] /= self.val[self.col_ptr[j]];
+            let yj = y[j];
+            for p in self.col_ptr[j] + 1..self.col_ptr[j + 1] {
+                y[self.row_idx[p]] -= self.val[p] * yj;
+            }
         }
         y
     }
@@ -98,11 +350,24 @@ pub(crate) fn bilinear(y: &[f64], y_prime: &[f64]) -> f64 {
 mod tests {
     use super::*;
 
+    /// Factorise a dense row-major matrix, for the goldens below.
+    fn dense(a: &[f64], n: usize) -> Option<Cholesky> {
+        let mut b = SymmetricBuilder::new();
+        for i in 0..n {
+            for j in 0..n {
+                if a[i * n + j] != 0.0 {
+                    b.add(i, j, a[i * n + j]);
+                }
+            }
+        }
+        Cholesky::factor(b, n)
+    }
+
     /// `[[4, 1], [1, 3]] z = [1, 2]` has `z = [1/11, 7/11]`, so the quadratic
     /// form `b^T A^-1 b` is `1 * 1/11 + 2 * 7/11 = 15/11`.
     #[test]
     fn reproduces_a_known_quadratic_form() {
-        let c = Cholesky::factor(vec![4.0, 1.0, 1.0, 3.0], 2).unwrap();
+        let c = dense(&[4.0, 1.0, 1.0, 3.0], 2).unwrap();
         let y = c.whiten(&[1.0, 2.0]);
         assert!((bilinear(&y, &y) - 15.0 / 11.0).abs() < 1e-12);
     }
@@ -113,8 +378,8 @@ mod tests {
     fn recovers_the_inverse_diagonal() {
         // A = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]; inverse diagonal is
         // [0.75, 1.0, 0.75].
-        let a = vec![2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0];
-        let c = Cholesky::factor(a, 3).unwrap();
+        let a = [2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0];
+        let c = dense(&a, 3).unwrap();
         for (i, expected) in [0.75, 1.0, 0.75].into_iter().enumerate() {
             let mut e = vec![0.0; 3];
             e[i] = 1.0;
@@ -127,8 +392,8 @@ mod tests {
     #[test]
     fn recovers_an_off_diagonal_covariance() {
         // Same A; (A^-1)_{0,1} = 0.5.
-        let a = vec![2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0];
-        let c = Cholesky::factor(a, 3).unwrap();
+        let a = [2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0];
+        let c = dense(&a, 3).unwrap();
         let y0 = c.whiten(&[1.0, 0.0, 0.0]);
         let y1 = c.whiten(&[0.0, 1.0, 0.0]);
         assert!((bilinear(&y0, &y1) - 0.5).abs() < 1e-12);
@@ -138,15 +403,108 @@ mod tests {
     /// A variance can never come out negative, because it is a sum of squares.
     #[test]
     fn a_quadratic_form_is_never_negative() {
-        let a = vec![1e12, 1e12 - 1.0, 1e12 - 1.0, 1e12];
-        let c = Cholesky::factor(a, 2).unwrap();
+        let a = [1e12, 1e12 - 1.0, 1e12 - 1.0, 1e12];
+        let c = dense(&a, 2).unwrap();
         let y = c.whiten(&[1.0, -1.0]);
         assert!(bilinear(&y, &y) >= 0.0);
     }
 
+    /// Against an independent dense reference, on random sparse SPD matrices.
+    ///
+    /// The goldens above are 2x2 and 3x3 — small enough that AMD does nothing
+    /// and no fill-in occurs, so they cannot catch a symbolic-pass bug. This
+    /// builds matrices big enough to permute and fill in, and checks every
+    /// bilinear form against a textbook dense factorisation of the *same*
+    /// matrix in its original order.
+    #[test]
+    fn agrees_with_a_dense_reference_on_random_sparse_systems() {
+        /// Dense Cholesky and quadratic form, deliberately naive: this is the
+        /// reference, so it must not share code with what it is checking.
+        fn dense_quadratic_form(a: &[f64], n: usize, b: &[f64], c: &[f64]) -> f64 {
+            let mut l = a.to_vec();
+            for j in 0..n {
+                let mut d = l[j * n + j];
+                for k in 0..j {
+                    d -= l[j * n + k] * l[j * n + k];
+                }
+                let d = d.sqrt();
+                l[j * n + j] = d;
+                for i in j + 1..n {
+                    let mut sum = l[i * n + j];
+                    for k in 0..j {
+                        sum -= l[i * n + k] * l[j * n + k];
+                    }
+                    l[i * n + j] = sum / d;
+                }
+            }
+            let solve = |rhs: &[f64]| -> Vec<f64> {
+                let mut y = rhs.to_vec();
+                for i in 0..n {
+                    for k in 0..i {
+                        y[i] -= l[i * n + k] * y[k];
+                    }
+                    y[i] /= l[i * n + i];
+                }
+                y
+            };
+            let (yb, yc) = (solve(b), solve(c));
+            yb.iter().zip(&yc).map(|(x, y)| x * y).sum()
+        }
+
+        // A cheap deterministic generator; no dependency, and reproducible.
+        let mut seed = 0x2545_F491_4F6C_DD1Du64;
+        let mut rand = move || {
+            seed ^= seed << 13;
+            seed ^= seed >> 7;
+            seed ^= seed << 17;
+            (seed >> 11) as f64 / (1u64 << 53) as f64
+        };
+
+        for n in [7usize, 23, 60] {
+            let mut a = vec![0.0f64; n * n];
+            // A chain plus scattered long-range couplings: the shape of a
+            // time-expanded joint, where a competitor's drift link can span
+            // the whole matrix.
+            for i in 0..n {
+                a[i * n + i] = 4.0 + rand();
+                if i + 1 < n {
+                    let v = -(0.5 + rand() * 0.5);
+                    a[i * n + i + 1] = v;
+                    a[(i + 1) * n + i] = v;
+                }
+            }
+            for step in 0..n / 3 {
+                let i = (step * 7) % n;
+                let j = (step * 29 + 3) % n;
+                if i != j {
+                    let v = -(0.1 + rand() * 0.2);
+                    a[i * n + j] = v;
+                    a[j * n + i] = v;
+                    // Keep it diagonally dominant, hence positive-definite.
+                    a[i * n + i] += 0.6;
+                    a[j * n + j] += 0.6;
+                }
+            }
+
+            let sparse = dense(&a, n).expect("spd");
+
+            for trial in 0..8 {
+                let b: Vec<f64> = (0..n).map(|_| rand() * 2.0 - 1.0).collect();
+                let c: Vec<f64> = (0..n).map(|_| rand() * 2.0 - 1.0).collect();
+                let got = bilinear(&sparse.whiten(&b), &sparse.whiten(&c));
+                let want = dense_quadratic_form(&a, n, &b, &c);
+                assert!(
+                    (got - want).abs() <= 1e-10 * want.abs().max(1.0),
+                    "n={n} trial={trial}: sparse {got} vs dense {want}"
+                );
+            }
+        }
+    }
+
+    /// A permutation must not change which matrices are rejected.
     #[test]
     fn rejects_a_non_positive_definite_matrix() {
         // Singular: the second row is a multiple of the first.
-        assert!(Cholesky::factor(vec![1.0, 2.0, 2.0, 4.0], 2).is_none());
+        assert!(dense(&[1.0, 2.0, 2.0, 4.0], 2).is_none());
     }
 }

tests/sparsity_measurement.rs

@@ -0,0 +1,178 @@
+//! What a sparse factorisation of the joint would actually buy (#52).
+//!
+//! Run explicitly:
+//!
+//! ```text
+//! cargo test --release --features approx,measure-sparsity \
+//!     --test sparsity_measurement -- --ignored --nocapture
+//! ```
+//!
+//! The whole file is gated: it reaches for the joint's sparsity pattern, which
+//! is exposed only under `measure-sparsity`.
+#![cfg(feature = "measure-sparsity")]
+
+use std::collections::HashSet;
+
+use trueskill_tt::{ConvergenceOptions, History};
+
+/// A history shaped like the issue's fixture: many slices, scored duels,
+/// competitors reappearing across slices so the drift links are long.
+fn fitted(slices: i64, duels: usize, competitors: usize) -> History<String> {
+    let mut h: History<String> = History::builder()
+        .key_type::<String>()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .gamma(0.05)
+        .convergence(ConvergenceOptions {
+            max_iter: trueskill_tt::ITERATIONS,
+            epsilon: 1e-8,
+            alpha: 1.0,
+        })
+        .build();
+
+    let mut k = 0usize;
+    for t in 0..slices {
+        for _ in 0..duels {
+            k += 1;
+            h.event(t)
+                .team([format!("p{}", k % competitors)])
+                .team([format!("p{}", (k + 37) % competitors)])
+                .scores([
+                    (k as f64 * 0.3).sin().abs() * 20.0,
+                    (k as f64 * 0.3).cos().abs() * 20.0,
+                ])
+                .commit()
+                .expect("ingests");
+        }
+    }
+    h.converge().expect("converges");
+    h
+}
+
+/// Symbolic Cholesky by row-merge: returns (nnz(L), flops).
+///
+/// Fill-in is simulated directly — for each column, the set of rows below the
+/// diagonal that are nonzero — which is exact and easily checked, at the cost
+/// of being O(n * nnz(L)) rather than the linear elimination-tree method.
+fn symbolic(n: usize, adj: &[HashSet<usize>], perm_of: &[usize]) -> (usize, f64) {
+    // `perm_of[old] = new`. Build the permuted lower-triangle pattern.
+    let mut cols: Vec<HashSet<usize>> = vec![HashSet::new(); n];
+    for (old, nbrs) in adj.iter().enumerate() {
+        let i = perm_of[old];
+        for &old_j in nbrs {
+            let j = perm_of[old_j];
+            if j < i {
+                cols[j].insert(i);
+            }
+        }
+    }
+
+    let mut nnz = 0usize;
+    let mut flops = 0.0f64;
+    for j in 0..n {
+        // Column j's pattern is final once every earlier column has merged in.
+        let rows: Vec<usize> = cols[j].iter().copied().collect();
+        let c = rows.len();
+        nnz += c + 1; // below-diagonal entries plus the diagonal
+        // Cholesky work for this column: one outer product over its pattern.
+        flops += (c as f64 + 1.0) * (c as f64 + 1.0);
+        // Fill-in: every pair in column j becomes an edge in the remaining graph.
+        for (a_idx, &a) in rows.iter().enumerate() {
+            for &b in &rows[a_idx + 1..] {
+                let (lo, hi) = if a < b { (a, b) } else { (b, a) };
+                cols[lo].insert(hi);
+            }
+        }
+    }
+    (nnz, flops)
+}
+
+#[test]
+#[ignore = "measurement, run explicitly"]
+fn what_sparsity_would_buy() {
+    for (slices, duels, competitors) in [(30, 8, 100), (76, 13, 200)] {
+        let h = fitted(slices, duels, competitors);
+        let (n, pattern) = h.joint_pattern_for_measurement();
+
+        let nnz_a: usize = pattern.iter().map(HashSet::len).sum::<usize>() + n;
+        let dense_flops = (n as f64).powi(3) / 3.0;
+
+        let natural: Vec<usize> = (0..n).collect();
+        let (nnz_nat, flops_nat) = symbolic(n, &pattern, &natural);
+
+        // AMD returns `perm[new] = old`; invert it.
+        let (col_ptr, row_idx) = csc(n, &pattern);
+        let p = feral_amd::amd_order(
+            &feral_amd::CscPattern::new(n, &col_ptr, &row_idx).expect("valid pattern"),
+        )
+        .expect("amd");
+        let mut perm_of = vec![0usize; n];
+        for (new, &old) in p.iter().enumerate() {
+            perm_of[old as usize] = new;
+        }
+        let (nnz_amd, flops_amd) = symbolic(n, &pattern, &perm_of);
+
+        println!(
+            "\n=== {slices} slices x {duels} duels, {competitors} competitors ===\n\
+             n              = {n}\n\
+             nnz(A)         = {nnz_a}  ({:.4}% dense)\n\
+             dense flops    = {:.3e}\n\
+             nnz(L) natural = {nnz_nat}   flops = {:.3e}  ({:.1}x vs dense)\n\
+             nnz(L) AMD     = {nnz_amd}   flops = {:.3e}  ({:.1}x vs dense)",
+            100.0 * nnz_a as f64 / (n * n) as f64,
+            dense_flops,
+            flops_nat,
+            dense_flops / flops_nat,
+            flops_amd,
+            dense_flops / flops_amd,
+        );
+    }
+}
+
+/// Full symmetric pattern to CSC, as `feral-amd` wants it.
+fn csc(n: usize, adj: &[HashSet<usize>]) -> (Vec<i32>, Vec<i32>) {
+    let mut col_ptr = Vec::with_capacity(n + 1);
+    let mut row_idx = Vec::new();
+    col_ptr.push(0i32);
+    for (j, nbrs) in adj.iter().enumerate() {
+        let mut rows: Vec<i32> = nbrs.iter().map(|&i| i as i32).collect();
+        rows.push(j as i32);
+        rows.sort_unstable();
+        rows.dedup();
+        row_idx.extend_from_slice(&rows);
+        col_ptr.push(row_idx.len() as i32);
+    }
+    (col_ptr, row_idx)
+}
+
+/// End-to-end factorisation time at the scale #52 was opened about.
+#[test]
+#[ignore = "measurement, run explicitly"]
+fn factorisation_time_at_scale() {
+    use std::time::Instant;
+
+    for (slices, duels, competitors) in [(30, 8, 100), (76, 13, 200), (150, 26, 400)] {
+        let h = fitted(slices, duels, competitors);
+        let (n, _) = h.joint_pattern_for_measurement();
+
+        // Warm, then time.
+        let _ = h.joint().expect("scored history");
+        let t = Instant::now();
+        let joint = h.joint().expect("scored history");
+        let factor = t.elapsed();
+
+        let a = "p0".to_string();
+        let b = "p1".to_string();
+        let t = Instant::now();
+        let _ = joint.posterior_of(&[(&a, 1.0), (&b, -1.0)]).expect("known");
+        let query = t.elapsed();
+
+        println!(
+            "n = {n:5}  factorise = {factor:>12?}   query = {query:>10?}   \
+             (dense was O(n^3): {:.3e} flops)",
+            (n as f64).powi(3) / 3.0
+        );
+    }
+}