A tiny non-zero drift returns a posterior variance 12 000x too small, as Ok #57

New Issue

2026-09-09T14:51:05Z

logaritmisk commented

2026-09-09 14:51:05 +00:00

src/history.rs:972 (collapse rule), :995-998 (assembly), src/joint.rs:52 (pivot rejection).

time_expanded_joint collapses consecutive appearances into one variable only when drift <= 0.0 exactly. For any smaller-but-positive drift it writes 1.0 / drift into the diagonal and off-diagonal. At drift = 1e-16 that is 1e16, and 1e16 + 0.28 rounds to 1e16 — the prior and contrast contributions are annihilated in the stored f64, so the matrix the solver sees is no longer the model's matrix.

Measured

8 competitors x 15 slices, sigma=6 beta=1 score_sigma=2 gamma=0.5, varying only Member::with_drift_scale(s), then History::joint() -> posterior_of(a - b):

`drift_scale`	drift var	f64 result	exact (mpmath, 200 dps)	error
0 (collapsed)	0	0.373056994818652843	0.373056994818653	—
1e-4	2.5e-9	0.373057041109483012		1.2e-7
1e-5	2.5e-11	0.373053396579676677		9.7e-6
1e-6	2.5e-13	0.372555558924591246		1.3e-3, silent
1e-7 … 1e-9	2.5e-15 … 2.5e-19	`Err(JointUnavailable)`	0.3730569948186…	misleading error
1e-10	2.5e-21	3.05172676872997e-5	0.373056994818653	12 200x too small, no error

At drift_scale = 1e-10 the caller is handed sigma = 0.0055 where the truth is 0.6108 — a 111x overconfident credible interval, returned as Ok.

This is not inherent conditioning

The same matrix assembled and solved in mpmath at 200 dps converges monotonically and smoothly to the collapsed value, and is flat from 1e-16 down to 1e-40:

drift_var = 1e-8   EXACT var = 0.373057083328942514
drift_var = 1e-12  EXACT var = 0.373056994827503879
drift_var = 1e-16  EXACT var = 0.373056994818653735
drift_var = 1e-40  EXACT var = 0.37305699481865285   <- the collapsed value

The quantity is perfectly well conditioned. Only the chosen representation — an explicit 1/drift precision — is not.

The obvious fix is wrong, and it was measured

Symmetric diagonal (Jacobi) equilibration does not help. Factoring D^-1/2 A D^-1/2 and whitening D^-1/2 c is equal or worse:

drift var	cond2	plain f64	Jacobi f64	exact
1e-6	6.5e8	0.986303749377213	0.986303748712236	0.986303749293778
1e-10	6.5e12	0.986300339105363	0.986295901649363	0.986301370100957
1e-12	6.4e14	0.986116684054720	0.985988196997571	0.986301369865393
1e-14	1.4e17	`None`	0.973624975751773	0.986301369863037

30x worse at 1e-10. The reason it cannot help: the damage happens at assembly, in the += at history.rs:995-998, before the factorisation ever runs. Scaling a matrix that has already lost its small entries recovers nothing.

Fix

The collapse rule needs a relative threshold rather than an exact-zero one: collapse when 1/drift swamps the other precisions in that row — i.e. when adding the prior and contrast terms to 1/drift would not change it. That is the condition that actually matters, and it degrades gracefully into the existing drift == 0 case.

Secondary defect

When Cholesky::factor returns None here, joint() reports:

the precision matrix is not positive-definite, which means a competitor has neither a proper prior nor any evidence

Every competitor in that fixture has a proper prior and plenty of evidence. The diagnosis is wrong and would send someone hunting the wrong bug for a long time.

Reachability

Fully public: Member::with_drift_scale (added in 617bc07) or History::builder().drift(ConstantDrift(1e-8)). Anyone modelling a "nearly static" competitor — a course layout, a reference bot — reaches for a small scale rather than exactly 0.0, and the docs for with_drift_scale actively suggest pinning as a use case.

Found by a floating-point audit, 2026-09-09.

`src/history.rs:972` (collapse rule), `:995-998` (assembly), `src/joint.rs:52` (pivot rejection). `time_expanded_joint` collapses consecutive appearances into one variable only when `drift <= 0.0` **exactly**. For any smaller-but-positive drift it writes `1.0 / drift` into the diagonal and off-diagonal. At `drift = 1e-16` that is `1e16`, and `1e16 + 0.28` rounds to `1e16` — the prior and contrast contributions are annihilated in the stored `f64`, so the matrix the solver sees is no longer the model's matrix. ## Measured 8 competitors x 15 slices, `sigma=6 beta=1 score_sigma=2 gamma=0.5`, varying only `Member::with_drift_scale(s)`, then `History::joint()` -> `posterior_of(a - b)`: | `drift_scale` | drift var | f64 result | exact (mpmath, 200 dps) | error | |---|---|---|---|---| | 0 (collapsed) | 0 | 0.373056994818652843 | 0.373056994818653 | — | | 1e-4 | 2.5e-9 | 0.373057041109483012 | | 1.2e-7 | | 1e-5 | 2.5e-11 | 0.373053396579676677 | | 9.7e-6 | | 1e-6 | 2.5e-13 | 0.372555558924591246 | | **1.3e-3, silent** | | 1e-7 … 1e-9 | 2.5e-15 … 2.5e-19 | **`Err(JointUnavailable)`** | 0.3730569948186… | misleading error | | **1e-10** | 2.5e-21 | **3.05172676872997e-5** | **0.373056994818653** | **12 200x too small, no error** | At `drift_scale = 1e-10` the caller is handed `sigma = 0.0055` where the truth is `0.6108` — a **111x overconfident credible interval**, returned as `Ok`. ## This is not inherent conditioning The same matrix assembled and solved in mpmath at 200 dps converges *monotonically and smoothly* to the collapsed value, and is flat from 1e-16 down to 1e-40: ``` drift_var = 1e-8 EXACT var = 0.373057083328942514 drift_var = 1e-12 EXACT var = 0.373056994827503879 drift_var = 1e-16 EXACT var = 0.373056994818653735 drift_var = 1e-40 EXACT var = 0.37305699481865285 <- the collapsed value ``` The quantity is perfectly well conditioned. Only the chosen *representation* — an explicit `1/drift` precision — is not. ## The obvious fix is wrong, and it was measured Symmetric diagonal (Jacobi) equilibration does **not** help. Factoring `D^-1/2 A D^-1/2` and whitening `D^-1/2 c` is equal or worse: | drift var | cond2 | plain f64 | Jacobi f64 | exact | |---|---|---|---|---| | 1e-6 | 6.5e8 | 0.986303749377213 | 0.986303748712236 | 0.986303749293778 | | 1e-10 | 6.5e12 | 0.986300339105363 | 0.986295901649363 | 0.986301370100957 | | 1e-12 | 6.4e14 | 0.986116684054720 | 0.985988196997571 | 0.986301369865393 | | 1e-14 | 1.4e17 | `None` | 0.973624975751773 | 0.986301369863037 | 30x worse at 1e-10. The reason it cannot help: the damage happens at **assembly**, in the `+=` at `history.rs:995-998`, before the factorisation ever runs. Scaling a matrix that has already lost its small entries recovers nothing. ## Fix The collapse rule needs a **relative** threshold rather than an exact-zero one: collapse when `1/drift` swamps the other precisions in that row — i.e. when adding the prior and contrast terms to `1/drift` would not change it. That is the condition that actually matters, and it degrades gracefully into the existing `drift == 0` case. ## Secondary defect When `Cholesky::factor` returns `None` here, `joint()` reports: > the precision matrix is not positive-definite, which means a competitor has neither a proper prior nor any evidence Every competitor in that fixture has a proper prior **and** plenty of evidence. The diagnosis is wrong and would send someone hunting the wrong bug for a long time. ## Reachability Fully public: `Member::with_drift_scale` (added in `617bc07`) or `History::builder().drift(ConstantDrift(1e-8))`. Anyone modelling a "nearly static" competitor — a course layout, a reference bot — reaches for a small scale rather than exactly `0.0`, and the docs for `with_drift_scale` actively suggest pinning as a use case. Found by a floating-point audit, 2026-09-09.

logaritmisk added the bug numerics labels 2026-09-09 14:53:47 +00:00

logaritmisk referenced this issue from a commit

2026-09-09 16:10:51 +00:00

fix!: collapse a drift too small to represent, on a relative threshold

logaritmisk closed this issue

2026-09-09 16:10:51 +00:00

logaritmisk referenced this issue from a commit

2026-09-09 16:10:52 +00:00

Merge branch 'fix/numerics-critical'

Sign in to join this conversation.