Hot-path Gaussian arithmetic round-trips through sqrt; convergence deltas computed in moment space #15

Closed
opened 2026-08-04 19:16:54 +00:00 by logaritmisk · 1 comment
Owner

The whole point of storing Gaussian in natural parameters (src/gaussian.rs:5-12) is that message passing becomes adds and subtracts with "no sqrt or reciprocal in the hot path". Mul and Div deliver on that. Add, Sub, and every convergence check do not.

1. Add/Sub compute a square root and then immediately square it away

src/gaussian.rs:130-144:

fn sub(self, rhs: Gaussian) -> Self::Output {
    let mu  = self.mu() + rhs.mu();                          // 2 divisions
    let var = self.sigma().powi(2) + rhs.sigma().powi(2);    // 2 sqrt, then squared back
    Self::from_ms(mu, var.sqrt())                            // 1 sqrt, then from_ms squares it
}

sigma() is 1.0 / pi.sqrt(), and .powi(2) undoes it. from_ms then does 1.0 / (sigma * sigma) — undoing var.sqrt(). Three square roots and several divisions to compute what is, in natural parameters:

pi'  = pi_a * pi_b / (pi_a + pi_b)
tau' = pi' * (tau_a/pi_a ± tau_b/pi_b)

This is the hot path. RankDiffFactor::propagate is a - b (src/factor/rank_diff.rs:26) and runs for every adjacent team pair on every forward and backward sweep, inside every EP iteration, inside every slice iteration, inside every history iteration. TeamSumFactor::propagate folds Add once per player (src/factor/team_sum.rs:20).

exclude and forget (src/gaussian.rs:88-105) have the same shape, and forget is called per player in the likelihood extraction (src/game.rs:383) and in Rating::performance.

Adding a from_mv(mu, var) constructor (skipping the sqrt→square round trip) is the minimal fix; a full natural-parameter formulation of Add/Sub removes the moment-space detour entirely. Care is needed around the existing pi <= 0 improper-Gaussian guards and the N00/N_INF identities — benches/gaussian.rs already benchmarks these ops, so the win is directly measurable.

2. Convergence is measured in moment space, against the spec

Factor::propagate returns (|Δmu|, |Δsigma|) (src/factor/mod.rs:62-64), computed by Gaussian::delta (src/gaussian.rs:81), which calls mu() and sigma() on both operands — 2 sqrts and 2 divisions per delta, once per factor per propagation, plus again per skill per slice iteration in History::iteration (src/history.rs:237, 256).

Spec §5 ("Stopping in natural-param space") explicitly calls for this to be (|Δpi|, |Δtau|):

  • mu and sigma are on different scales; one tolerance is wrong for both
  • We store in nat-params anyway — checking convergence in mu/sigma costs free sqrts
  • Nat-param delta is the natural geometry of the EP fixed point

None of that was implemented. Switching removes the sqrts and fixes the scale mismatch, at the cost of re-tuning test epsilons (the spec anticipates this: "tests' epsilons re-tuned mechanically") and a slightly different exact stopping threshold.

While changing this, note that the NaN-blindness of the comparison needs fixing regardless — see #8.

3. Team performance is computed twice per game

run_chain folds each team's weighted performance into arena.team_prior (src/game.rs:290-295), then the likelihood extraction re-folds the identical sum from scratch for every team (src/game.rs:374-377):

let performance = players.iter().zip(weights.iter())
    .fold(N00, |p, (player, &w)| p + (player.performance() * w));

The value is already sitting in arena.team_prior[si], and si is in hand on the line above. Each performance() call is itself a forget (another moment-space round trip), so this doubles that cost across every player in the game.

Suggested order

  1. from_mv + natural-parameter Add/Sub/exclude/forget — pure win, no semantic change beyond floating-point association.
  2. Reuse arena.team_prior in the likelihood loop — pure win.
  3. Natural-parameter convergence deltas — behaviour-affecting, needs epsilon re-tuning.

Acceptance

  • benches/gaussian.rs shows a measurable drop on Gaussian::add/sub.
  • benches/batch.rs and benches/history_converge.rs show no regression, ideally an improvement.
  • Existing numerical goldens still pass (steps 1–2 within current epsilons; step 3 with re-tuned epsilons and a documented threshold change).
The whole point of storing `Gaussian` in natural parameters (`src/gaussian.rs:5-12`) is that message passing becomes adds and subtracts with "no `sqrt` or reciprocal in the hot path". `Mul` and `Div` deliver on that. `Add`, `Sub`, and every convergence check do not. ## 1. `Add`/`Sub` compute a square root and then immediately square it away `src/gaussian.rs:130-144`: ```rust fn sub(self, rhs: Gaussian) -> Self::Output { let mu = self.mu() + rhs.mu(); // 2 divisions let var = self.sigma().powi(2) + rhs.sigma().powi(2); // 2 sqrt, then squared back Self::from_ms(mu, var.sqrt()) // 1 sqrt, then from_ms squares it } ``` `sigma()` is `1.0 / pi.sqrt()`, and `.powi(2)` undoes it. `from_ms` then does `1.0 / (sigma * sigma)` — undoing `var.sqrt()`. Three square roots and several divisions to compute what is, in natural parameters: ``` pi' = pi_a * pi_b / (pi_a + pi_b) tau' = pi' * (tau_a/pi_a ± tau_b/pi_b) ``` **This is the hot path.** `RankDiffFactor::propagate` is `a - b` (`src/factor/rank_diff.rs:26`) and runs for every adjacent team pair on every forward and backward sweep, inside every EP iteration, inside every slice iteration, inside every history iteration. `TeamSumFactor::propagate` folds `Add` once per player (`src/factor/team_sum.rs:20`). `exclude` and `forget` (`src/gaussian.rs:88-105`) have the same shape, and `forget` is called per player in the likelihood extraction (`src/game.rs:383`) and in `Rating::performance`. Adding a `from_mv(mu, var)` constructor (skipping the sqrt→square round trip) is the minimal fix; a full natural-parameter formulation of `Add`/`Sub` removes the moment-space detour entirely. Care is needed around the existing `pi <= 0` improper-Gaussian guards and the `N00`/`N_INF` identities — `benches/gaussian.rs` already benchmarks these ops, so the win is directly measurable. ## 2. Convergence is measured in moment space, against the spec `Factor::propagate` returns `(|Δmu|, |Δsigma|)` (`src/factor/mod.rs:62-64`), computed by `Gaussian::delta` (`src/gaussian.rs:81`), which calls `mu()` and `sigma()` on both operands — 2 sqrts and 2 divisions per delta, once per factor per propagation, plus again per skill per slice iteration in `History::iteration` (`src/history.rs:237`, `256`). Spec §5 ("Stopping in natural-param space") explicitly calls for this to be `(|Δpi|, |Δtau|)`: > - `mu` and `sigma` are on different scales; one tolerance is wrong for both > - We store in nat-params anyway — checking convergence in mu/sigma costs free sqrts > - Nat-param delta is the natural geometry of the EP fixed point None of that was implemented. Switching removes the sqrts *and* fixes the scale mismatch, at the cost of re-tuning test epsilons (the spec anticipates this: "tests' epsilons re-tuned mechanically") and a slightly different exact stopping threshold. While changing this, note that the NaN-blindness of the comparison needs fixing regardless — see #8. ## 3. Team performance is computed twice per game `run_chain` folds each team's weighted performance into `arena.team_prior` (`src/game.rs:290-295`), then the likelihood extraction re-folds the identical sum from scratch for every team (`src/game.rs:374-377`): ```rust let performance = players.iter().zip(weights.iter()) .fold(N00, |p, (player, &w)| p + (player.performance() * w)); ``` The value is already sitting in `arena.team_prior[si]`, and `si` is in hand on the line above. Each `performance()` call is itself a `forget` (another moment-space round trip), so this doubles that cost across every player in the game. ## Suggested order 1. `from_mv` + natural-parameter `Add`/`Sub`/`exclude`/`forget` — pure win, no semantic change beyond floating-point association. 2. Reuse `arena.team_prior` in the likelihood loop — pure win. 3. Natural-parameter convergence deltas — behaviour-affecting, needs epsilon re-tuning. ## Acceptance - `benches/gaussian.rs` shows a measurable drop on `Gaussian::add`/`sub`. - `benches/batch.rs` and `benches/history_converge.rs` show no regression, ideally an improvement. - Existing numerical goldens still pass (steps 1–2 within current epsilons; step 3 with re-tuned epsilons and a documented threshold change).
Author
Owner

Fixed in 355cdb7 — all three items.

  1. Add, Sub, exclude and forget go through new from_mv(mu, var) and variance(), both pure divisions. No square roots remain in the variance path.
  2. run_chain reads the team performance from arena.team_prior[si] instead of re-folding it.
  3. Natural-parameter convergence deltas were not done — see below.

Measured on this machine, same fixtures, before and after:

before after
Batch::iteration 23.57 µs 19.31 µs (−18%)
scored_history_60_events 1.071 ms 983 µs (−8%)

The Gaussian::add/sub microbenchmarks in benches/gaussian.rs cannot resolve this change — they sit at ~234 ps against a ~218 ps harness floor that mul/div also hit, so a single operation is lost in overhead. Worth knowing before reading those numbers as evidence of anything; only the end-to-end benches showed the effect.

One golden moved, and the new value is the correct one: two identical competitors drawing must land on their shared prior mean exactly by symmetry, and the root-free path returns 25.0 where the reference transcription recorded 24.999999 — that value rounded to six decimals. Asserting a six-decimal transcription at epsilon 1e-6 left no headroom.

Item 3 left open, tracked in #22 alongside the other convergence-reporting work: switching Factor::propagate to (|Δpi|, |Δtau|) changes the stopping threshold and needs every test epsilon re-tuned, which is a larger and more disruptive change than the two pure wins here. The NaN-blindness half of it is already fixed (#8).

Fixed in 355cdb7 — all three items. 1. `Add`, `Sub`, `exclude` and `forget` go through new `from_mv(mu, var)` and `variance()`, both pure divisions. No square roots remain in the variance path. 2. `run_chain` reads the team performance from `arena.team_prior[si]` instead of re-folding it. 3. Natural-parameter convergence deltas were **not** done — see below. Measured on this machine, same fixtures, before and after: | | before | after | |---|---:|---:| | `Batch::iteration` | 23.57 µs | 19.31 µs (−18%) | | `scored_history_60_events` | 1.071 ms | 983 µs (−8%) | The `Gaussian::add`/`sub` microbenchmarks in `benches/gaussian.rs` **cannot resolve this change** — they sit at ~234 ps against a ~218 ps harness floor that `mul`/`div` also hit, so a single operation is lost in overhead. Worth knowing before reading those numbers as evidence of anything; only the end-to-end benches showed the effect. One golden moved, and the new value is the correct one: two identical competitors drawing must land on their shared prior mean exactly by symmetry, and the root-free path returns 25.0 where the reference transcription recorded 24.999999 — that value rounded to six decimals. Asserting a six-decimal transcription at epsilon 1e-6 left no headroom. **Item 3 left open**, tracked in #22 alongside the other convergence-reporting work: switching `Factor::propagate` to `(|Δpi|, |Δtau|)` changes the stopping threshold and needs every test epsilon re-tuned, which is a larger and more disruptive change than the two pure wins here. The NaN-blindness half of it is already fixed (#8).
logaritmisk added the numericsperformance labels 2026-09-07 13:53:18 +00:00
Sign in to join this conversation.