From c8f2b16b4c61a2030607ef502ec57852bbbbb945 Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Sat, 15 Aug 2026 21:46:22 +0200 Subject: [PATCH] fix: prevent probability underflow in pairwise cost matrix Replaces premature exponentiation-based row normalization with log-sum-exp arithmetic to prevent tiny probabilities from collapsing to exactly zero. This eliminates spurious infinite costs for valid but rare transitions while preserving correct IEEE 754 semantics for genuinely unobserved pairs. Adds explicit guards against NaN in degenerate rows and includes a regression test verifying finite costs for probabilities as low as 1e-200. --- .../theory/evolutionary_distances/index.html | 46 ++++++++-- DevDocMD/theory/evolutionary_distances.md | 47 +++++++++-- src/obikphylo/src/cardcomp.rs | 84 ++++++++++++++----- 3 files changed, 145 insertions(+), 32 deletions(-) diff --git a/DevDoc/theory/evolutionary_distances/index.html b/DevDoc/theory/evolutionary_distances/index.html index c56c4b9e..71d10aa8 100644 --- a/DevDoc/theory/evolutionary_distances/index.html +++ b/DevDoc/theory/evolutionary_distances/index.html @@ -4100,12 +4100,46 @@ precision (not the model file's truncated 6 decimals). Written alongside _iqtree.model/_iqtree.fasta from the same CompactAlphabet both already use, so there is no risk of the three files disagreeing. Covered by states_csv_maps_compact_symbols_back_to_canonical_ones.

-

The zero-exchangeability pattern itself (state 0 in the user's report, -frequency 3.26%, R=0 with every other state) is not yet explained — -plausibly a genuinely unobserved transition in the calibration -(cardinality_transitions/composition_transitions count 0 for every -pair involving it), which is a legitimate, if numerically extreme, result -of -ln(0), not necessarily a bug — not investigated further.

+

Root cause of the zero-exchangeability pattern found and fixed +(2026-08-15): premature exp() in pairwise_cost_matrix underflowed +merely-tiny probabilities to exactly 0.0. The user also reported +iqtree3 emitting "Numerical underflow for lh-derivative" warnings on +the same run — a real signal, traced to obikphylo/src/cardcomp.rs's +pairwise_cost_matrix, not to the frequency computation (which is a +plain, safe f64 division, never close to underflow at any realistic +scale). The function already accumulated log_p in log-space (correct), +but then row-normalised by exponentiating each cell first +(raw[a][b] = log_p.exp()) and summing the results — f64::exp hard +underflows to exactly 0.0 for any input below roughly -709, which a +sum of several individually-small-but-nonzero probability factors +(composition/cardinality terms, best_pairing_cost's pairing terms) can +reach easily on real, skewed calibration data. Once raw[a][b] was +exactly 0.0, normalisation and -ln turned a merely tiny probability +into a +∞ cost indistinguishable from a literally unobserved +transition (p == 0.0 exactly, e.g. p_comp[i][j] never once tallied) — +conflating two different things: "never observed" (should be +∞, a +correct MLE result) and "observed, but the joint probability of this +multi-step transition is extremely small" (should be a large finite +cost).

+

Fix: row-normalise via log-sum-exp instead of exponentiating first — +row_max = max_b(log_p[a][b]), log_sum = row_max + ln(Σ_b +exp(log_p[a][b] - row_max)) (every shifted term is in (0,1], so this +never underflows for a finite log_p[a][b]), then +cost[a][b] = log_sum - log_p[a][b] directly — no intermediate +probability is ever materialised. This falls out of IEEE 754 arithmetic +without a special case: a genuinely-unobserved factor (log_p[a][b] == +-∞, from the existing if p > 0.0 {...} else { NEG_INFINITY } guards +already in the log-accumulation loop) still yields cost = +∞ exactly +(finite − (−∞) = +∞), preserving the correct semantics for that case, +while every merely-tiny-but-nonzero transition now gets a large but +finite cost. A degenerate all--∞ row (a state with literally zero +probability of transitioning to anything, row_max == -∞) is guarded +explicitly to avoid a -∞ − (-∞) = NaN in the log-sum-exp itself. +Covered by cardcomp::tests::underflow_prone_transition_gets_finite_cost_not_infinite +(all off-diagonal composition probabilities set to 1e-200, well past +where the old exp()-first code would have underflowed to 0.0, cost +asserted finite). Every pre-existing cardcomp test still passes +unchanged (numerically identical results when no underflow occurs).

References

The Mash mutation-rate model this discussion contrasts with: (Fan et al. 2015; Marbl Lab 2026)1 2.

diff --git a/DevDocMD/theory/evolutionary_distances.md b/DevDocMD/theory/evolutionary_distances.md index 23d1805d..d2d227bf 100644 --- a/DevDocMD/theory/evolutionary_distances.md +++ b/DevDocMD/theory/evolutionary_distances.md @@ -2090,12 +2090,47 @@ precision (not the model file's truncated 6 decimals). Written alongside already use, so there is no risk of the three files disagreeing. Covered by `states_csv_maps_compact_symbols_back_to_canonical_ones`. -The zero-exchangeability pattern itself (state 0 in the user's report, -frequency 3.26%, `R=0` with every other state) is not yet explained — -plausibly a genuinely unobserved transition in the calibration -(`cardinality_transitions`/`composition_transitions` count `0` for every -pair involving it), which is a legitimate, if numerically extreme, result -of `-ln(0)`, not necessarily a bug — not investigated further. +**Root cause of the zero-exchangeability pattern found and fixed +(2026-08-15): premature `exp()` in `pairwise_cost_matrix` underflowed +merely-tiny probabilities to exactly `0.0`.** The user also reported +`iqtree3` emitting "Numerical underflow for lh-derivative" warnings on +the same run — a real signal, traced to `obikphylo/src/cardcomp.rs`'s +`pairwise_cost_matrix`, not to the frequency computation (which is a +plain, safe `f64` division, never close to underflow at any realistic +scale). The function already accumulated `log_p` in log-space (correct), +but then row-normalised by exponentiating each cell *first* +(`raw[a][b] = log_p.exp()`) and summing the results — `f64::exp` hard +underflows to exactly `0.0` for any input below roughly `-709`, which a +sum of several individually-small-but-nonzero probability factors +(composition/cardinality terms, `best_pairing_cost`'s pairing terms) can +reach easily on real, skewed calibration data. Once `raw[a][b]` was +exactly `0.0`, normalisation and `-ln` turned a merely tiny probability +into a `+∞` cost indistinguishable from a *literally* unobserved +transition (`p == 0.0` exactly, e.g. `p_comp[i][j]` never once tallied) — +conflating two different things: "never observed" (should be `+∞`, a +correct MLE result) and "observed, but the joint probability of this +multi-step transition is extremely small" (should be a large *finite* +cost). + +**Fix**: row-normalise via log-sum-exp instead of exponentiating first — +`row_max = max_b(log_p[a][b])`, `log_sum = row_max + ln(Σ_b +exp(log_p[a][b] - row_max))` (every shifted term is in `(0,1]`, so this +never underflows for a finite `log_p[a][b]`), then +`cost[a][b] = log_sum - log_p[a][b]` directly — no intermediate +probability is ever materialised. This falls out of IEEE 754 arithmetic +without a special case: a genuinely-unobserved factor (`log_p[a][b] == +-∞`, from the existing `if p > 0.0 {...} else { NEG_INFINITY }` guards +already in the log-accumulation loop) still yields `cost = +∞` exactly +(`finite − (−∞) = +∞`), preserving the correct semantics for that case, +while every merely-tiny-but-nonzero transition now gets a large but +*finite* cost. A degenerate all-`-∞` row (a state with literally zero +probability of transitioning to anything, `row_max == -∞`) is guarded +explicitly to avoid a `-∞ − (-∞) = NaN` in the log-sum-exp itself. +Covered by `cardcomp::tests::underflow_prone_transition_gets_finite_cost_not_infinite` +(all off-diagonal composition probabilities set to `1e-200`, well past +where the old `exp()`-first code would have underflowed to `0.0`, cost +asserted finite). Every pre-existing `cardcomp` test still passes +unchanged (numerically identical results when no underflow occurs). ## References diff --git a/src/obikphylo/src/cardcomp.rs b/src/obikphylo/src/cardcomp.rs index 9a6db64b..e89d7b27 100644 --- a/src/obikphylo/src/cardcomp.rs +++ b/src/obikphylo/src/cardcomp.rs @@ -146,52 +146,66 @@ fn best_pairing_cost(lost: &[u8], gained: &[u8], p_comp: &[[f64; 4]; 4]) -> f64 /// only by composition matching (shared-base retention and paired /// substitutions), never by a state pair's cardinality difference alone. pub fn pairwise_cost_matrix(p_card: &[[f64; 5]; 5], p_comp: &[[f64; 4]; 4], free_loss: bool) -> [[f64; 16]; 16] { - let mut raw = [[0.0f64; 16]; 16]; + let mut log_p = [[0.0f64; 16]; 16]; // ln(P), *before* row-normalisation for a in 0u8..16 { for b in 0u8..16 { let shared = a & b; let lost: Vec = (0..4).filter(|&i| a & (1 << i) != 0 && b & (1 << i) == 0).collect(); let gained: Vec = (0..4).filter(|&i| b & (1 << i) != 0 && a & (1 << i) == 0).collect(); - let mut log_p = 0.0; // accumulate ln(P), so 0.0 = probability 1 + let mut lp = 0.0; // accumulate ln(P), so 0.0 = probability 1 if !free_loss { let card_a = a.count_ones() as usize; let card_b = b.count_ones() as usize; let p_c = p_card[card_a][card_b]; - log_p += if p_c > 0.0 { p_c.ln() } else { f64::NEG_INFINITY }; + lp += if p_c > 0.0 { p_c.ln() } else { f64::NEG_INFINITY }; } for i in 0..4u8 { if shared & (1 << i) != 0 { let p = p_comp[i as usize][i as usize]; - log_p += if p > 0.0 { p.ln() } else { f64::NEG_INFINITY }; + lp += if p > 0.0 { p.ln() } else { f64::NEG_INFINITY }; } } - log_p -= best_pairing_cost(&lost, &gained, p_comp); + lp -= best_pairing_cost(&lost, &gained, p_comp); - raw[a as usize][b as usize] = log_p.exp(); + log_p[a as usize][b as usize] = lp; } } - // Row-normalise to a proper transition probability matrix. - let mut p = [[0.0f64; 16]; 16]; - for a in 0..16 { - let row_sum: f64 = raw[a].iter().sum(); - if row_sum > 0.0 { - for b in 0..16 { - p[a][b] = raw[a][b] / row_sum; - } - } - } - - // Cost = -ln(P), then symmetrise (see doc comment: parsimony on an - // unrooted tree requires a symmetric cost matrix). + // Row-normalise and convert to cost entirely in log-space + // (log-sum-exp), never exponentiating a raw `log_p` value directly — + // for a transition reachable only through several low-probability + // steps, `log_p.exp()` can underflow to exactly `0.0` (anything below + // roughly `-709` does, in `f64`), silently turning a real, if small, + // probability into a hard `+∞` cost. Observed in practice: a + // calibrated matrix with many such exact-zero entries is numerically + // unstable for IQ-TREE's own likelihood-derivative computation + // ("Numerical underflow for lh-derivative"). `log_sum_exp` never + // exponentiates anything above `0` (every term is shifted by the + // row's own max first), so it stays accurate across the full range + // `f64` can represent, not just what survives a direct `exp()`. let mut cost = [[0.0f64; 16]; 16]; for a in 0..16 { + let row_max = log_p[a].iter().cloned().fold(f64::NEG_INFINITY, f64::max); + if row_max == f64::NEG_INFINITY { + // Every transition out of this state has probability 0 in the + // calibration data — genuinely unreachable, not underflow. + cost[a] = [f64::INFINITY; 16]; + continue; + } + let log_sum = row_max + log_p[a].iter().map(|&lp| (lp - row_max).exp()).sum::().ln(); for b in 0..16 { - cost[a][b] = if p[a][b] > 0.0 { -p[a][b].ln() } else { f64::INFINITY }; + // `log_sum - log_p[a][b]` is `+∞` automatically when + // `log_p[a][b] == -∞` (IEEE 754: finite − (−∞) = +∞) — no + // separate branch needed for a genuinely zero-probability + // transition. + cost[a][b] = log_sum - log_p[a][b]; } } + + // Symmetrise (see doc comment: parsimony on an unrooted tree requires + // a symmetric cost matrix). let mut sym = [[0.0f64; 16]; 16]; for a in 0..16 { for b in 0..16 { @@ -289,4 +303,34 @@ mod tests { } } } + + #[test] + fn underflow_prone_transition_gets_finite_cost_not_infinite() { + // Every off-diagonal composition probability is tiny (1e-200) but + // not exactly zero — small enough that the naive `log_p.exp()` + // (exponentiate before row-normalising) hard-underflows to `0.0` + // in `f64` well before the probability is truly zero, silently + // turning a real, if minuscule, probability into `+∞` cost. The + // log-sum-exp normalisation in `pairwise_cost_matrix` must keep + // this finite instead. + let p_card = [[0.2f64; 5]; 5]; + let tiny = 1e-200; + let mut p_comp = [[tiny; 4]; 4]; + for i in 0..4 { + p_comp[i][i] = 1.0 - 3.0 * tiny; + } + + let cost = pairwise_cost_matrix(&p_card, &p_comp, false); + + // {A,C} -> {G,T}: no shared bases, both members substituted — the + // shape most prone to underflow (several tiny-probability factors + // multiplied together). + let a = 0b0011u8; // A, C + let b = 0b1100u8; // G, T + assert!( + cost[a as usize][b as usize].is_finite(), + "cost must stay finite for a merely tiny (not exactly zero) probability, got {}", + cost[a as usize][b as usize] + ); + } }