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] + ); + } }