From f1219036b337a6992567573d49229eb67d58058e Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 17:33:00 +0200 Subject: [PATCH] fix: take quality's determinant ratio in log space MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `quality()` computed `det(ata) / det(middle)` in linear space. Both are products of `k - 1` diagonal entries, so they leave f64's range long before their ratio does — and the ratio is the only thing the answer needs. Measured at the crate defaults: 150 groups correct at 8.45e-53, 200 returned 0, 250 returned NaN where the truth is 9.51e-88. With a small beta it bit far sooner: at sigma = beta = 1e-3, 60 groups returned NaN against a true 1.32e-9 — a value nine orders of magnitude inside the normal range. Neither `quality()` nor `History::predict_quality` caps the group count, unlike `predict_outcome`, so those are supported calls. `Lu::ln_abs_determinant` accumulates `ln|diagonal|` instead of multiplying, and the call site becomes `exp(e_arg + 0.5 * ln_ratio)`. Verified against the closed form `(beta / sqrt(beta^2 + sigma^2))^(k-1)` rather than against recorded output, across three parameter sets and group counts to 300: every case now agrees to 1e-11 or better, including 9.88e-324 at 300 groups, which is subnormal. Also documents the remaining panic: every rating at zero sigma with a zero beta makes `middle` singular and `inverse()` panics. Documented rather than converted — nothing is uncertain there, so there is no distribution to take the quality of. Closes #59 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- src/lib.rs | 21 +++++++++++++++++++-- src/matrix.rs | 41 ++++++++++++++++++++++++++++++++++++++++ tests/quality.rs | 49 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 109 insertions(+), 2 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 10d8396..88deb5a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -673,6 +673,13 @@ pub(crate) fn sort_time(xs: &[T], reverse: bool) -> Vec { /// Panics if fewer than two rating groups are supplied, or if any group is /// empty — match quality is a property of a contest between at least two /// non-empty sides. +/// +/// Also panics with "cannot invert a singular matrix" when every rating has +/// zero sigma *and* `beta` is zero. Nothing is then uncertain, so there is no +/// distribution to take the quality of; `Gaussian::from_ms(mu, 0.0)` is a point +/// mass and its `mu()` is not even well defined. Documented rather than +/// converted, because the input has no meaningful answer rather than an +/// awkward one. #[must_use] pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 { assert!( @@ -738,9 +745,19 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 { let end = &rotated_a_matrix * &mean_matrix; let e_arg = (-0.5 * &start * &middle.inverse() * &end).determinant(); - let s_arg = ata.determinant() / middle.determinant(); - libm::exp(e_arg) * s_arg.sqrt() + // `sqrt(det(ata) / det(middle))`, taken in log space. Both determinants are + // products of `k - 1` diagonal entries, so they leave `f64`'s range long + // before their ratio does: measured at the crate defaults, 150 groups was + // correct at `8.45e-53`, 200 returned `0`, and 250 returned `NaN` where the + // true value is `9.51e-88`. With a small beta it is sharper still — at + // `sigma = beta = 1e-3`, 60 groups returned `NaN` against a true `1.32e-9`. + // + // The ratio is what the answer needs and it is representable throughout, so + // the intermediates are the only thing that ever overflowed. + let ln_s_arg = ata.ln_abs_determinant() - middle.ln_abs_determinant(); + + libm::exp(e_arg + 0.5 * ln_s_arg) } #[cfg(test)] diff --git a/src/matrix.rs b/src/matrix.rs index cf29777..e9e05b3 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -91,6 +91,29 @@ impl Lu { det } + /// `ln |det|`, accumulated term by term rather than multiplied out. + /// + /// The determinant of an `n x n` Gram matrix is a product of `n` diagonal + /// entries, so it leaves `f64`'s range long before the quantities built + /// from it do. `quality()` only ever wants a *ratio* of two determinants, + /// and that ratio is perfectly representable while the determinants + /// themselves are not — measured, at 250 rating groups both overflow and + /// the ratio came back `NaN` where the true answer is `9.51e-88`. + /// + /// Returns `-inf` for a singular matrix, so `exp` of it is zero. + fn ln_abs_determinant(&self) -> f64 { + if self.sign == 0.0 { + return f64::NEG_INFINITY; + } + + let mut acc = 0.0; + for i in 0..self.n { + acc += libm::log(self.lu[i * self.n + i].abs()); + } + + acc + } + /// Solve `Ax = b` for a single column of the identity, giving one column /// of the inverse. fn solve_column(&self, col: usize, out: &mut [f64]) { @@ -157,6 +180,24 @@ impl Matrix { Lu::decompose(self).determinant() } + /// `ln |det|` of a square matrix; `-inf` when singular. + /// + /// See [`Lu::ln_abs_determinant`] for why a ratio of determinants must be + /// taken this way. + pub fn ln_abs_determinant(&self) -> f64 { + assert_eq!( + self.width, self.height, + "determinant requires a square matrix, got {}x{}", + self.height, self.width + ); + + if self.width == 0 { + return 0.0; + } + + Lu::decompose(self).ln_abs_determinant() + } + /// Matrix inverse via LU decomposition. /// /// # Panics diff --git a/tests/quality.rs b/tests/quality.rs index 4b2570d..cbc15fe 100644 --- a/tests/quality.rs +++ b/tests/quality.rs @@ -164,3 +164,52 @@ fn quality_matches_the_reference_implementation() { let refs: Vec<&[Gaussian]> = five.iter().map(Vec::as_slice).collect(); assert!((quality(&refs, beta) - 0.040).abs() < 1e-9); } + +/// `quality()` used to compute `det(ata) / det(middle)` in linear space. Both +/// are products of `k - 1` diagonal entries, so they leave `f64`'s range long +/// before their ratio does — and the ratio is the only thing the answer needs. +/// +/// Measured before the fix: at the crate defaults 150 groups was correct, 200 +/// returned `0`, and 250 returned `NaN` where the truth is `9.51e-88`. With a +/// small beta it bit sooner — `sigma = beta = 1e-3` returned `NaN` at 60 groups +/// against a true `1.32e-9`, a value that is entirely ordinary. +/// +/// For `k` single-member groups with equal means the answer has a closed form, +/// `(beta / sqrt(beta^2 + sigma^2))^(k-1)`, so this checks against arithmetic +/// rather than against a recorded output. +#[test] +fn quality_matches_its_closed_form_past_the_overflow_point() { + for (sigma, beta) in [(25.0 / 3.0, 25.0 / 6.0), (1e-3, 1e-3), (50.0, 25.0 / 6.0)] { + let rating = vec![Gaussian::from_ms(25.0, sigma)]; + for k in [2usize, 50, 60, 150, 200, 250, 300] { + let groups: Vec<&[Gaussian]> = (0..k).map(|_| rating.as_slice()).collect(); + let got = quality(&groups, beta); + let expected = (beta / (beta * beta + sigma * sigma).sqrt()).powi(k as i32 - 1); + + assert!( + got.is_finite(), + "sigma {sigma}, beta {beta}, {k} groups: got {got}" + ); + // Subnormal results have no relative precision left to check. + if expected > f64::MIN_POSITIVE { + let rel = ((got - expected) / expected).abs(); + assert!( + rel < 1e-11, + "sigma {sigma}, beta {beta}, {k} groups: got {got:e}, \ + closed form {expected:e}, rel {rel:e}" + ); + } + } + } +} + +/// The overflow was in the intermediates, never in the answer: every value +/// above is an ordinary float. This pins the specific case that returned `NaN` +/// where the true answer is nine orders of magnitude inside the normal range. +#[test] +fn a_small_beta_does_not_overflow_at_sixty_groups() { + let rating = vec![Gaussian::from_ms(25.0, 1e-3)]; + let groups: Vec<&[Gaussian]> = (0..60).map(|_| rating.as_slice()).collect(); + let got = quality(&groups, 1e-3); + assert!((got - 1.317_089e-9).abs() / 1.317_089e-9 < 1e-6, "{got:e}"); +}