From 6139061740adc1edc35ed9dce267e34a1a3aeb8f Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 17:42:16 +0200 Subject: [PATCH] fix: keep the truncated variance representable in the far tail MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `v_w` returned `w` and let `trunc` form `1 - w`. `w` tends to 1 out in the tail, so that subtraction lost about log10(alpha^2) digits — and the quantity it was destroying is perfectly representable. Two separate cancellations, fixed separately. The non-tie half: `half_line_truncation` now returns `1 - w` computed symbolically rather than as `1 - v*gap`. With `alpha*gap = 1 - inv^2*b` the leading ones cancel on paper instead of in floating point. Measured against the exact truncated variance: alpha before after 1e6 8.9e-5 rel 0.0 rel (exact) 1e8 returns 0.0 0.0 rel (exact) At 1e8 the old form gave `sigma_trunc = 0`, and `from_ms(mu, 0.0)` is a point mass whose `mu()` is inf/inf = NaN. `beta(1e-8).sigma(1e-8)` with priors 1000 apart went from Err + NaN skills to a finite fit. The tie half is a different subtraction — `w = v^2 - u`, where both grow as alpha^2 while their difference stays O(1). The existing escape hatch could not cover it: it keys on `alpha * width >= HALF_LINE_WINDOW`, how many window-widths from the mean the window sits, and a NARROW window fails that however deep it is. Measured at alpha 1e6 with a 1e-6 window it kept four digits and returned `1 - w = -2.4e-4` where the truth is +2.8e-13. One step earlier it was quietly wrong instead: `1 - w = 1.0` exactly, a truncation reported as a no-op, where the truth was 5e-17. Over a narrow window the density is a truncated exponential in `s = (x - alpha)/width`, whose mean and variance are closed forms, so `v = alpha + width*m(t)` and `1 - w = width^2 * V(t)` with no large subtraction at all. Validated against high-precision quadrature: v exact to 4e-10, `1 - w` to 4e-10 across the region it is used in. The crossover is on `alpha / width` rather than on either alone, because that ratio is what says how many digits the subtraction has left — and the approximation is most accurate exactly where the subtraction is worst, since both improve as the window narrows. Defaults are bit-identical (pi 0.02398318151216503 before and after). Tests: the three reproductions from the issue, the narrow-window form against pinned quadrature values, and a continuity sweep across all three tie branches — a misplaced crossover is the real risk here, and a jump at a boundary is visible even without pinning absolute values. Closes #60 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- src/lib.rs | 181 ++++++++++++++++++++++++++++++++++-- tests/non_finite_results.rs | 54 ++++++++++- 2 files changed, 224 insertions(+), 11 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 88deb5a..1feca81 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -218,6 +218,12 @@ const HALF_LINE_WINDOW: f64 = 10.0; /// four-term series is good to ~1e-10 by here, so the two are at their closest /// agreement around this point. Below it the subtraction is exact; above it the /// series is. +/// `alpha / width` past which the tie branch's `v^2 - u` has lost too many +/// digits to trust, and the narrow-window form takes over. +/// +/// The subtraction retains about `(width / alpha)^2 / EPSILON` of its +/// precision, so this is the ratio at which that falls below roughly 1e-6. +const NARROW_WINDOW_RATIO: f64 = 2.0e4; const ASYMPTOTIC_MILLS_ALPHA: f64 = 100.0; pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0); @@ -476,10 +482,72 @@ fn pdf(x: f64, mu: f64, sigma: f64) -> f64 { fn half_line_truncation(alpha: f64) -> (f64, f64) { let inv = alpha.recip(); let inv_sq = inv * inv; - let gap = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq))); + let b = 2.0 - inv_sq * (10.0 - 74.0 * inv_sq); + let gap = inv * (1.0 - inv_sq * b); let v = alpha + gap; - (v, v * gap) + // Returns `1 - w`, not `w`, and that is the whole point of this shape. + // + // `w` tends to 1 out here, so a caller forming `1 - w` loses about + // `log10(alpha^2)` digits: measured against the exact truncated variance, + // `1 - w` came back with 8.9e-5 relative error at alpha = 1e6 and **0.0** + // from alpha = 1e8 — where the true value is 1e-16 and perfectly + // representable. `sigma * (1 - w).sqrt()` was then exactly zero, and + // `from_ms(mu, 0.0)` is a point mass whose `mu()` is `inf/inf = NaN`. + // + // Expanding `1 - v*gap` symbolically removes the subtraction: with + // `alpha*gap = 1 - inv^2*b`, the leading ones cancel on paper instead of in + // floating point, leaving `inv^2` times a bracket that tends to 1. Measured + // exact — 0.0 relative error — from alpha = 1e3 to 1e8. + let one_minus_w = inv_sq + * ((1.0 - inv_sq * (10.0 - 74.0 * inv_sq)) + 2.0 * inv_sq * b - inv_sq * inv_sq * b * b); + + (v, one_minus_w) +} + +/// Truncation to a *narrow* window `[alpha, alpha + d]`, as `(v, 1 - w)`. +/// +/// The tie branch forms `w` from `v^2 - u`, and both grow as `alpha^2` while +/// their difference stays `O(1)`. Far enough into the tail that subtraction has +/// nothing left: measured at `alpha = 1e6` with a window of `1e-6` it kept four +/// significant digits and returned `1 - w = -2.4e-4` where the truth is +/// `+2.8e-13`, so `sqrt` of it was NaN. One step earlier it was quietly wrong +/// instead — `1 - w = 1.0` exactly, a truncation reported as a no-op, where the +/// truth was `5e-17`. +/// +/// The existing half-line escape hatch does not cover it, because that keys on +/// `alpha * d >= HALF_LINE_WINDOW` — how many window-widths from the mean the +/// window sits — and a *narrow* window fails that however deep it is. +/// +/// Over a narrow window the density is `exp(-t*s - s^2 d^2 / 2)` in +/// `x = alpha + s*d`, with `t = alpha * d`. Dropping the `d^2` term leaves a +/// truncated exponential on `[0, 1]`, whose mean and variance are closed forms. +/// So `v = alpha + d*m(t)` and `1 - w = d^2 * V(t)`, with no subtraction of +/// large quantities anywhere. +/// +/// Measured against high-precision quadrature over `alpha` in `[1e2, 1e9]`: +/// `v` exact to 4e-10 or better, `1 - w` to 4e-10 across the region this is +/// used in. +fn narrow_window_truncation(alpha: f64, d: f64) -> (f64, f64) { + let t = alpha * d; + + // `m` and `V` are the mean and variance of a truncated exponential on + // [0, 1] with rate `t`, both of which cancel as `t -> 0`. The series is + // their limit (1/2 and 1/12, a uniform window) with the leading correction. + let (m, v_s) = if t < 1e-3 { + ( + 0.5 - t / 12.0 + t * t * t / 720.0, + 1.0 / 12.0 - t * t / 240.0, + ) + } else { + let em1 = libm::expm1(t); + ( + 1.0 / t - 1.0 / em1, + 1.0 / (t * t) - (em1 + 1.0) / (em1 * em1), + ) + }; + + (alpha + d * m, d * d * v_s) } fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { @@ -507,7 +575,7 @@ fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { (v, v - alpha) }; - (v, v * gap) + (v, 1.0 - v * gap) } else { // v is odd in mu and w is even, so fold to mu <= 0. Both truncation // points then sit in the upper tail, where the scaled form applies. @@ -523,9 +591,22 @@ fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { // Once the window sits many of its own widths into the tail it is // indistinguishable from a half-line, so the asymptotic covers it with // no subtraction at all. - if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * (beta - alpha) >= HALF_LINE_WINDOW { - let (v, w) = half_line_truncation(alpha); - return (if flipped { -v } else { v }, w); + let width = beta - alpha; + + if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * width >= HALF_LINE_WINDOW { + let (v, one_minus_w) = half_line_truncation(alpha); + return (if flipped { -v } else { v }, one_minus_w); + } + + // A narrow window deep in the tail: too narrow for the half-line above, + // too deep for the subtraction below. The direct form keeps roughly + // `1 / (alpha/width)^2` of its digits, so the crossover is on that + // ratio rather than on either quantity alone — and the approximation is + // most accurate exactly where the subtraction is worst, since both + // improve as the window narrows. + if alpha > 0.0 && alpha > NARROW_WINDOW_RATIO * width { + let (v, one_minus_w) = narrow_window_truncation(alpha, width); + return (if flipped { -v } else { v }, one_minus_w); } let (v, u) = if alpha > 0.0 { @@ -548,17 +629,23 @@ fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { ) }; - let w = -(u - v.powi(2)); + // `1 - w` where `w = v^2 - u`. Both `v^2` and `u` grow as alpha^2 while + // their difference stays O(1), so this subtraction is the one place the + // tie branch can still lose everything — see the escape hatch above, + // which is what keeps the far tail away from it. + let one_minus_w = 1.0 + u - v.powi(2); - (if flipped { -v } else { v }, w) + (if flipped { -v } else { v }, one_minus_w) } } fn trunc(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { - let (v, w) = v_w(mu, sigma, margin, tie); + // `v_w` returns `1 - w` rather than `w`: forming the difference here is + // what destroyed the truncated variance in the far tail. + let (v, one_minus_w) = v_w(mu, sigma, margin, tie); let mu_trunc = mu + sigma * v; - let sigma_trunc = sigma * (1.0 - w).sqrt(); + let sigma_trunc = sigma * one_minus_w.sqrt(); (mu_trunc, sigma_trunc) } @@ -762,6 +849,80 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 { #[cfg(test)] mod tests { + /// The truncated variance must stay a variance across every branch, and + /// the branches must agree where they meet. + /// + /// `v_w` now has three regimes for a tie — half-line, narrow-window, and + /// the direct subtraction — and a misplaced crossover between them is the + /// failure mode this guards. A jump at a boundary is visible here even + /// though the absolute values are not pinned. + #[test] + fn truncated_variance_is_continuous_across_the_tie_branches() { + for &alpha in &[50.0, 99.0, 100.0, 101.0, 1e3, 1e5, 1e6] { + // Sweep the window width across NARROW_WINDOW_RATIO and the + // half-line threshold, which sit at different widths per alpha. + let mut previous: Option<(f64, f64)> = None; + let mut width = alpha / (NARROW_WINDOW_RATIO * 100.0); + while width < 40.0 / alpha { + // mu = 0 puts the window at [-margin, margin]; shift it out to + // `alpha` by moving the mean instead. + let margin = width * 0.5; + let mu = -(alpha + width * 0.5); + let (v, one_minus_w) = v_w(mu, 1.0, margin, true); + + assert!(v.is_finite(), "alpha {alpha}, width {width:e}: v = {v}"); + assert!( + one_minus_w.is_finite() && one_minus_w > 0.0 && one_minus_w <= 1.0, + "alpha {alpha}, width {width:e}: 1 - w = {one_minus_w:e} is not a variance" + ); + + if let Some((pv, pw)) = previous { + // Consecutive widths differ by 2x, so the moments may not + // differ by more than a small multiple of that. + assert!( + one_minus_w / pw < 32.0 && pw / one_minus_w < 32.0, + "alpha {alpha}: 1 - w jumped from {pw:e} to {one_minus_w:e} \ + at width {width:e} — a branch boundary is misplaced" + ); + assert!( + (v - pv).abs() <= 8.0 * width.max(1e-12) + 1e-9 * v.abs(), + "alpha {alpha}: v jumped from {pv} to {v} at width {width:e}" + ); + } + previous = Some((v, one_minus_w)); + width *= 2.0; + } + } + } + + /// The narrow-window form against high-precision quadrature. + /// + /// These are the inputs where the direct `v^2 - u` subtraction had four + /// significant digits left and returned a negative variance. + #[test] + fn narrow_window_truncation_matches_quadrature() { + for &(alpha, d, expect_v, expect_w) in &[ + (1e6, 2e-6, 1_000_000.000_000_687, 2.759_383_390_335_666e-13), + (1e4, 1e-6, 10_000.000_000_499_167, 8.333_291_666_831_727e-14), + ( + 1e3, + 1e-5, + 1_000.000_004_991_666_6, + 8.333_291_666_803_818e-12, + ), + ] { + let (v, one_minus_w) = narrow_window_truncation(alpha, d); + assert!( + ((v - expect_v) / expect_v).abs() < 1e-12, + "alpha {alpha:e}: v = {v}, want {expect_v}" + ); + assert!( + ((one_minus_w - expect_w) / expect_w).abs() < 1e-8, + "alpha {alpha:e}: 1 - w = {one_minus_w:e}, want {expect_w:e}" + ); + } + } + /// A NaN must survive the fold from ANY position, not only the last. /// /// The fold runs over a `HashMap`, so "last" is per-process hash order. The diff --git a/tests/non_finite_results.rs b/tests/non_finite_results.rs index 2434172..b7db0f0 100644 --- a/tests/non_finite_results.rs +++ b/tests/non_finite_results.rs @@ -10,7 +10,9 @@ //! `!tuple_gt(..)`. These tests pin the guard from outside. use smallvec::smallvec; -use trueskill_tt::{Event, Gaussian, History, InferenceError, Member, Outcome, Team}; +use trueskill_tt::{ + ConstantDrift, Event, Gaussian, History, InferenceError, Member, Outcome, Team, +}; fn scored_fit( sigma: f64, @@ -163,3 +165,53 @@ fn a_nan_competitor_is_not_masked_by_a_healthy_one() { "{err:?}" ); } + +/// A tie observed with a narrow draw margin between far-apart competitors must +/// produce a fit, not NaN skills. +/// +/// The tie branch forms the truncated variance from `v^2 - u`, and both grow as +/// `alpha^2` while their difference stays `O(1)`. Deep enough into the tail +/// that subtraction had four digits left: measured, it returned `1 - w` +/// negative and `sqrt` of it was NaN. The half-line escape hatch did not cover +/// it, because that keys on how many window-widths from the mean the window +/// sits and a narrow window fails that however deep it is. +/// +/// These parameters are ordinary for a precise-scoring domain, and the +/// neighbouring wider-margin case always worked — so this was a cliff, not +/// "extreme inputs break". +#[test] +fn a_narrow_draw_margin_far_into_the_tail_still_fits() { + for (beta, p_draw, sd, gap) in [ + (1e-2, 1e-8, 1e-2, 10.0), + (1e-3, 1e-9, 1e-3, 1.0), + (1e-4, 1e-12, 1e-4, 1.0), + ] { + let mut h = History::builder() + .mu(0.0) + .sigma(sd) + .beta(beta) + .p_draw(p_draw) + .drift(ConstantDrift(0.0)) + .build(); + h.add_events(vec![Event { + time: 1i64, + teams: smallvec![ + Team::with_members([Member::new("a").with_prior(Gaussian::from_ms(0.0, sd))]), + Team::with_members([Member::new("b").with_prior(Gaussian::from_ms(gap, sd))]), + ], + outcome: Outcome::draw(2), + }]) + .unwrap(); + + let report = h + .converge() + .unwrap_or_else(|e| panic!("beta {beta:e}, p_draw {p_draw:e}: {e:?}")); + assert!(report.converged); + + let skill = h.current_skill(&"a").unwrap(); + assert!( + skill.mu().is_finite() && skill.sigma().is_finite() && skill.sigma() > 0.0, + "beta {beta:e}, p_draw {p_draw:e}: {skill:?}" + ); + } +}