logaritmisk/trueskill-tt
1
0
Fork 0
Code Issues 1 Pull Requests Packages Projects Releases Activity

The joint precision matrix is 0.19% dense and we factorise it densely #52

New Issue
Closed
opened 2026-09-08 05:54:35 +00:00 by logaritmisk · 2 comments
logaritmisk commented 2026-09-08 05:54:35 +00:00
Owner
Copy Link
Copy Source

Split out of #51, which asked for caching. Caching is done. This is the other half, and it is the one that makes the first query expensive.

Measured

Instrumenting time_expanded_joint on a synthetic history shaped like ustat's hole-grain model — 76 slices, 988 scored duels, 200 competitors:

n=1976  nnz=7504  density=0.1922%  bandwidth=207
dense Cholesky: 745 ms

7,504 nonzeros in a 3,904,576-entry matrix. We allocate 31 MB, fill 99.81% of it with zeros, and then run an O(n^3) factorisation over the whole thing. At ustat's ~4,000 appearances that is 128 MB and ~5.9 s.

The structure is not accidental — it is what a time-expanded factor graph is. A row couples to its own previous and next appearance through the drift link, and to whoever co-appeared in its slice. Nothing else. The density falls as the history grows, so this gets worse with scale, not better.

Why bandwidth is not the answer

207 on this fixture, which would make a banded Cholesky ~15x faster. But bandwidth here is data-dependent and degenerates: one competitor who appears in slice 0 and then not again until slice 75 creates a drift link spanning nearly the whole matrix, and the band swallows everything between. A real history has those. Fill-reducing ordering is the robust version of the same idea.

Scouting, and what it does not establish

Registry metadata and docs only — nothing here was compiled or benchmarked, so treat all of it as a shortlist rather than a finding.

The minimal-dependency path looks unusually good. feral-amd is only an AMD ordering — 2 crates total, its core has zero dependencies, it is #![forbid(unsafe_code)] like we are, MIT, no build script, and its test suite asserts byte-for-byte equality with the SuiteSparse AMD reference including the internal counters. Its API is amd_order(&CscPattern) -> Vec<usize>. That would let us keep our own ~40-line Cholesky and just permute, which is by far the smallest change that could work. amd 0.2.2 is the older BSD-3 direct SuiteSparse port it was validated against — 2M downloads but stale since 2022.

The full-solver option is faer 0.24, which has exactly the right API shape (symbolic analysis with AMD once, numeric factorisation once, unlimited solves) and Par as an explicit per-call argument rather than an ambient thread pool. Against it: gemm and pulp are required rather than optional, so it is ~30-40 crates replacing our current four, and raw-cpuid in the required tree means runtime CPU-feature dispatch. That is a real problem for us specifically. Our just determinism job checks bit-identity across thread counts on one machine; it would not catch a factorisation whose last bits differ between an AVX-512 host and an AVX2 one. We already took the libm-over-std dependency precisely to avoid that class of drift, so accepting it here would undo that decision.

Two to avoid: nalgebra-sparse (603K downloads) disqualifies itself in its own docs — "performs no fill-in reduction... not currently recommended to use this implementation for serious projects" — which is strictly worse than what we have. And sprs-ldl, the top crates.io hit for "sparse cholesky", is LGPL-2.1 and stale since 2022.

Proposed order

  1. Reorder with AMD, keep our own Cholesky, and measure. If that closes most of the gap it costs 2 dependencies and no unsafe.
  2. Only if it does not: exploit sparsity in the factorisation itself, which is a real sparse-Cholesky implementation and a much larger change.
  3. faer last, and only with a bitwise cross-architecture check first.

Not urgent. With #51 done, a batch of queries pays this once rather than per query, so the practical cost is already bounded.

Split out of #51, which asked for caching. Caching is done. This is the other half, and it is the one that makes the *first* query expensive. ## Measured Instrumenting `time_expanded_joint` on a synthetic history shaped like ustat's hole-grain model — 76 slices, 988 scored duels, 200 competitors: ``` n=1976 nnz=7504 density=0.1922% bandwidth=207 dense Cholesky: 745 ms ``` 7,504 nonzeros in a 3,904,576-entry matrix. We allocate 31 MB, fill 99.81% of it with zeros, and then run an O(n^3) factorisation over the whole thing. At ustat's ~4,000 appearances that is 128 MB and ~5.9 s. The structure is not accidental — it is what a time-expanded factor graph *is*. A row couples to its own previous and next appearance through the drift link, and to whoever co-appeared in its slice. Nothing else. The density falls as the history grows, so this gets worse with scale, not better. ## Why bandwidth is not the answer 207 on this fixture, which would make a banded Cholesky ~15x faster. But bandwidth here is data-dependent and degenerates: one competitor who appears in slice 0 and then not again until slice 75 creates a drift link spanning nearly the whole matrix, and the band swallows everything between. A real history has those. Fill-reducing ordering is the robust version of the same idea. ## Scouting, and what it does not establish Registry metadata and docs only — **nothing here was compiled or benchmarked**, so treat all of it as a shortlist rather than a finding. The minimal-dependency path looks unusually good. `feral-amd` is *only* an AMD ordering — 2 crates total, its core has zero dependencies, it is `#![forbid(unsafe_code)]` like we are, MIT, no build script, and its test suite asserts byte-for-byte equality with the SuiteSparse AMD reference including the internal counters. Its API is `amd_order(&CscPattern) -> Vec<usize>`. That would let us keep our own ~40-line Cholesky and just permute, which is by far the smallest change that could work. `amd` 0.2.2 is the older BSD-3 direct SuiteSparse port it was validated against — 2M downloads but stale since 2022. The full-solver option is `faer` 0.24, which has exactly the right API shape (symbolic analysis with AMD once, numeric factorisation once, unlimited solves) and `Par` as an explicit per-call argument rather than an ambient thread pool. Against it: `gemm` and `pulp` are required rather than optional, so it is ~30-40 crates replacing our current four, and `raw-cpuid` in the required tree means runtime CPU-feature dispatch. That is a real problem for us specifically. Our `just determinism` job checks bit-identity across thread counts *on one machine*; it would not catch a factorisation whose last bits differ between an AVX-512 host and an AVX2 one. We already took the libm-over-std dependency precisely to avoid that class of drift, so accepting it here would undo that decision. Two to avoid: `nalgebra-sparse` (603K downloads) disqualifies itself in its own docs — *"performs no fill-in reduction... not currently recommended to use this implementation for serious projects"* — which is strictly worse than what we have. And `sprs-ldl`, the top crates.io hit for "sparse cholesky", is **LGPL-2.1** and stale since 2022. ## Proposed order 1. Reorder with AMD, keep our own Cholesky, and measure. If that closes most of the gap it costs 2 dependencies and no unsafe. 2. Only if it does not: exploit sparsity in the factorisation itself, which is a real sparse-Cholesky implementation and a much larger change. 3. `faer` last, and only with a bitwise cross-architecture check first. Not urgent. With #51 done, a batch of queries pays this once rather than per query, so the practical cost is already bounded.
logaritmisk added the enhancementnumericsperformance labels 2026-09-08 08:26:38 +00:00
logaritmisk commented 2026-09-10 04:59:49 +00:00
Author
Owner
Copy Link
Copy Source

Step 1 as written cannot work, and measuring says so before any code is written. The plan needs reordering — not the matrix, the plan.

AMD on its own buys exactly nothing

Cholesky::factor is dense. Its inner loops are for k in 0..j { d -= a[j*n+k] * a[j*n+k] } — unconditional, every k, zero or not. A permutation changes which entries are zero; it does not change how many multiplications happen. Measured on a 700×700 banded matrix at 0.43% density, factorised in band order and then under a scramble that destroys the band:

banded    = 30.196 ms
scrambled = 29.544 ms

Identical within noise, as the flop count says it must be. AMD reduces fill-in, and fill-in only costs anything once the factorisation is skipping zeros. Permuting a dense factorisation is pure overhead.

What sparsity is worth, and what AMD is worth on top of it

Symbolic factorisation on two fixtures, the second reproducing this issue's numbers exactly (n=1976, nnz=7504, 0.1922%):

=== 30 slices x 8 duels, 100 competitors ===
n              = 480
nnz(A)         = 1720  (0.7465% dense)
dense flops    = 3.686e7
nnz(L) natural = 24428   flops = 2.160e6  (17.1x vs dense)
nnz(L) AMD     = 1913   flops = 8.111e3  (4544.9x vs dense)

=== 76 slices x 13 duels, 200 competitors ===
n              = 1976
nnz(A)         = 7504  (0.1922% dense)
dense flops    = 2.572e9
nnz(L) natural = 292437   flops = 5.597e7  (46.0x vs dense)
nnz(L) AMD     = 11583   flops = 8.656e4  (29709.6x vs dense)

So, on the fixture this issue was opened about:

speedup vs the dense factorisation we have
AMD alone (step 1 as proposed) 1x — measured
sparse, natural ordering 46x
sparse + AMD 29,710x

AMD is worth 646x on top of sparsity and nothing without it. The ordering in the issue is inverted: sparsity is the enabler, AMD is the multiplier. Note also how badly natural ordering fills in — 292,437 nonzeros in L against 7,504 in A, which is the drift links spanning the matrix exactly as the "why bandwidth is not the answer" section predicted.

Revised plan

Steps 1 and 2 are one step. Doing sparse Cholesky and AMD together, keeping our own factorisation (the scouting conclusion still holds — feral itself pulls pulp, so it has the same runtime CPU-dispatch problem that ruled out faer, and sprs-ldl is still LGPL).

feral-amd checks out on inspection: #![forbid(unsafe_code)] in both it and feral-ordering-core, two crates total, amd_order(&CscPattern) -> Result<Vec<i32>, OrderingError>.

One consequence worth flagging now: reordering changes the summation order, so joint goldens will move in their last bits. That is expected, not a regression — but it means the goldens need checking against the analytic reference rather than against their current values.

Measurement lives in tests/sparsity_measurement.rs behind a measure-sparsity feature.

**Step 1 as written cannot work, and measuring says so before any code is written.** The plan needs reordering — not the matrix, the plan. ## AMD on its own buys exactly nothing `Cholesky::factor` is dense. Its inner loops are `for k in 0..j { d -= a[j*n+k] * a[j*n+k] }` — unconditional, every `k`, zero or not. A permutation changes which entries are zero; it does not change how many multiplications happen. Measured on a 700×700 banded matrix at 0.43% density, factorised in band order and then under a scramble that destroys the band: ``` banded = 30.196 ms scrambled = 29.544 ms ``` Identical within noise, as the flop count says it must be. AMD reduces *fill-in*, and fill-in only costs anything once the factorisation is skipping zeros. Permuting a dense factorisation is pure overhead. ## What sparsity is worth, and what AMD is worth on top of it Symbolic factorisation on two fixtures, the second reproducing this issue's numbers exactly (n=1976, nnz=7504, 0.1922%): ``` === 30 slices x 8 duels, 100 competitors === n = 480 nnz(A) = 1720 (0.7465% dense) dense flops = 3.686e7 nnz(L) natural = 24428 flops = 2.160e6 (17.1x vs dense) nnz(L) AMD = 1913 flops = 8.111e3 (4544.9x vs dense) === 76 slices x 13 duels, 200 competitors === n = 1976 nnz(A) = 7504 (0.1922% dense) dense flops = 2.572e9 nnz(L) natural = 292437 flops = 5.597e7 (46.0x vs dense) nnz(L) AMD = 11583 flops = 8.656e4 (29709.6x vs dense) ``` So, on the fixture this issue was opened about: | | speedup vs the dense factorisation we have | |---|---| | AMD alone (step 1 as proposed) | **1x** — measured | | sparse, natural ordering | 46x | | sparse + AMD | **29,710x** | AMD is worth **646x on top of** sparsity and **nothing without it**. The ordering in the issue is inverted: sparsity is the enabler, AMD is the multiplier. Note also how badly natural ordering fills in — 292,437 nonzeros in L against 7,504 in A, which is the drift links spanning the matrix exactly as the "why bandwidth is not the answer" section predicted. ## Revised plan Steps 1 and 2 are one step. Doing sparse Cholesky **and** AMD together, keeping our own factorisation (the scouting conclusion still holds — `feral` itself pulls `pulp`, so it has the same runtime CPU-dispatch problem that ruled out `faer`, and `sprs-ldl` is still LGPL). `feral-amd` checks out on inspection: `#![forbid(unsafe_code)]` in both it and `feral-ordering-core`, two crates total, `amd_order(&CscPattern) -> Result<Vec<i32>, OrderingError>`. One consequence worth flagging now: reordering changes the summation order, so joint goldens will move in their last bits. That is expected, not a regression — but it means the goldens need checking against the analytic reference rather than against their current values. Measurement lives in `tests/sparsity_measurement.rs` behind a `measure-sparsity` feature.
logaritmisk commented 2026-09-10 05:08:38 +00:00
Author
Owner
Copy Link
Copy Source

Done — 695bb82 (merged as 1629176). 745 ms → 1.11 ms on this issue's fixture.

Sparse up-looking Cholesky with an AMD ordering, both together, since the measurement above showed neither is worth anything without the other.

Measured end to end, factorising through History::joint

n =  480     215 µs      query =  3.2 µs
n = 1976    1.112 ms     query = 11.9 µs      (was ~745 ms — 670x)
n = 7800    4.616 ms     query = 55.6 µs      (dense: 1.58e11 flops)

The benchmark agrees: joint_factorise_480_appearances went 9.11 ms → 167 µs (54x), and joint_query_480_appearances is 2.4 µs.

Scaling is near-linear now instead of cubic: 16x the variables costs 21x the time, where dense would have cost 4096x. The 128 MB allocation at ustat's ~4000 appearances is gone with it — storage is nnz(L), which AMD keeps at ~1.5x nnz(A).

On the scouting

Your shortlist held, with one addition: feral itself pulls pulp, so it has the same runtime CPU-dispatch problem that ruled out faer — a factorisation whose last bits differ between an AVX-512 host and an AVX2 one, which is the drift the libm-over-std decision was made to avoid. sprs-ldl is still LGPL, nalgebra-sparse still disclaims fill-reduction in its own docs, and the SuiteSparse bindings are all FFI, which #![forbid(unsafe_code)] rules out. So the factorisation is written here (Davis, Direct Methods for Sparse Linear Systems) and only the ordering is a dependency. feral-amd was exactly as advertised: two crates, both #![forbid(unsafe_code)], amd_order(&CscPattern) -> Result<Vec<i32>, _>.

Two things worth recording

The matrix accumulates into a BTreeMap, not a hash map. The iteration order becomes the factorisation's summation order, and a hash map's order varies per process — which is precisely what tests/cross_process_determinism.rs was written for.

whiten leaves its result in the permuted order. A dot product does not care which order its operands are in, as long as both were permuted the same way, so bilinear is untouched and the un-permute is skipped entirely.

Correctness

The existing analytic goldens are 2×2 and 3×3 — too small for AMD to do anything or for any fill-in to occur, so they could not have caught a symbolic-pass bug even though they all still pass. agrees_with_a_dense_reference_on_random_sparse_systems fills that gap: chain-plus-long-range matrices up to n=60, every bilinear form checked against a deliberately naive dense factorisation that shares no code with the thing it is checking.

The joint goldens elsewhere in the suite did not move, which surprised me — I expected last-bit drift from the changed summation order and flagged it above. The full suite, all four feature combinations plus release, and just determinism at 1/2/4/8 threads are green unchanged.

The measurement lives in tests/sparsity_measurement.rs behind the measure-sparsity feature.

Done — 695bb82 (merged as 1629176). **745 ms → 1.11 ms** on this issue's fixture. Sparse up-looking Cholesky with an AMD ordering, both together, since the measurement above showed neither is worth anything without the other. ## Measured end to end, factorising through `History::joint` ``` n = 480 215 µs query = 3.2 µs n = 1976 1.112 ms query = 11.9 µs (was ~745 ms — 670x) n = 7800 4.616 ms query = 55.6 µs (dense: 1.58e11 flops) ``` The benchmark agrees: `joint_factorise_480_appearances` went 9.11 ms → 167 µs (54x), and `joint_query_480_appearances` is 2.4 µs. Scaling is near-linear now instead of cubic: 16x the variables costs 21x the time, where dense would have cost 4096x. The 128 MB allocation at ustat's ~4000 appearances is gone with it — storage is `nnz(L)`, which AMD keeps at ~1.5x `nnz(A)`. ## On the scouting Your shortlist held, with one addition: **`feral` itself pulls `pulp`**, so it has the same runtime CPU-dispatch problem that ruled out `faer` — a factorisation whose last bits differ between an AVX-512 host and an AVX2 one, which is the drift the libm-over-std decision was made to avoid. `sprs-ldl` is still LGPL, `nalgebra-sparse` still disclaims fill-reduction in its own docs, and the SuiteSparse bindings are all FFI, which `#![forbid(unsafe_code)]` rules out. So the factorisation is written here (Davis, *Direct Methods for Sparse Linear Systems*) and only the *ordering* is a dependency. `feral-amd` was exactly as advertised: two crates, both `#![forbid(unsafe_code)]`, `amd_order(&CscPattern) -> Result<Vec<i32>, _>`. ## Two things worth recording **The matrix accumulates into a `BTreeMap`, not a hash map.** The iteration order becomes the factorisation's summation order, and a hash map's order varies per process — which is precisely what `tests/cross_process_determinism.rs` was written for. **`whiten` leaves its result in the permuted order.** A dot product does not care which order its operands are in, as long as both were permuted the same way, so `bilinear` is untouched and the un-permute is skipped entirely. ## Correctness The existing analytic goldens are 2×2 and 3×3 — too small for AMD to do anything or for any fill-in to occur, so they could not have caught a symbolic-pass bug even though they all still pass. `agrees_with_a_dense_reference_on_random_sparse_systems` fills that gap: chain-plus-long-range matrices up to n=60, every bilinear form checked against a deliberately naive dense factorisation that shares no code with the thing it is checking. The joint goldens elsewhere in the suite did **not** move, which surprised me — I expected last-bit drift from the changed summation order and flagged it above. The full suite, all four feature combinations plus release, and `just determinism` at 1/2/4/8 threads are green unchanged. The measurement lives in `tests/sparsity_measurement.rs` behind the `measure-sparsity` feature.
Sign in to join this conversation.
No Branch/Tag Specified
Branches Tags
Labels
Clear labels
api
Shape or ergonomics of the public surface
 blocked
Waiting on another issue or an external decision
 breaking
Resolving it changes the public API; wants a major bump
 bug
Behaviour is wrong, including silently wrong
 decision
Needs a call to be made, not code to be written
 docs
Documentation, including doc comments and the README
 enhancement
New capability the crate does not have
 infra
CI, release process, packaging
 numerics
Floating-point accuracy, cancellation, convergence maths
 performance
Speed or memory
 refactor
Internal cleanup with no intended behaviour change
 testing
Test coverage and test infrastructure
No labels enhancement numerics performance
Milestone
No items
No Milestone
Projects
Clear projects
No projects
Assignees
Clear assignees
logaritmisk
No Assignees
1 Participants
Notifications
Due Date
No due date set.
Dependencies

No dependencies set.

Reference: logaritmisk/trueskill-tt#52