Numerical-Stability Report: `feature_engineering/numerical.py` Moments¶

Date: 2026-04-24 Scope: Audit + benchmark of moment-computation kernels in mlframe.feature_engineering.numerical for catastrophic cancellation and numerical precision issues. Companion to the recent multi-output landing's Welford ensembling refactor.

Note: the numba moment kernels referenced as numerical.py:NNN below have since moved into the sibling mlframe/feature_engineering/_numerical_numba.py (re-exported from numerical.py). Current homes: compute_numerical_aggregates_numba at _numerical_numba.py:~48, compute_moments_slope_mi at _numerical_numba.py:~770 (including the weighted-moments bug-fix site, now at _numerical_numba.py:~646-657). compute_simple_stats_numba remains in numerical.py:~277.

TL;DR¶

Kahan-Babuška-Neumaier compensated 2-pass variance is the clear winner for mean / var — beats naive on every input distribution tested, numerically exact (relative error 0) on 6 of 7 distributions.
Welford-Pébay generalised online moments wins for skew / kurt on the hard catastrophic cases (170-315× better) but slightly worse on easy cases (long arrays of normal data).
One real bug found and fixed: numerical.py:740 — weighted_skew += w_summand * w_d * next_weight was double-accumulating skew and never accumulating kurt. Affects every caller using compute_moments_slope_mi with weights.
Runtime cost ~2× for Welford / Kahan vs naive — acceptable for the precision gains on hot-path features.

Audit findings (10 total, 7 HIGH severity)¶

#	Function	File:line	Severity	Issue
1	`compute_simple_stats_numba`	numerical.py:147-156	HIGH	Two-pass var, but uncompensated inner sum
2	`compute_moments_slope_mi`	numerical.py:724-741	HIGH	Uncompensated skew/kurt sums + cancellation in `std^3` divider
3	`compute_moments_slope_mi`	numerical.py:705-709	HIGH	Uncompensated `slope_over` / `slope_under` / `r_sum`
4	`compute_numerical_aggregates_numba`	numerical.py:281-291, 325-328	HIGH	Uncompensated quadratic / cubic / harmonic mean accumulators
5	`compute_numerical_aggregates_numba`	numerical.py:340-354	MED	Heuristic geomean log-mode threshold, no log-sum-exp guard
6	`compute_moments_slope_mi`	numerical.py:796	MED	Intercept = `mean_value - slope * xvals_mean` cancellation
7	`compute_moments_slope_mi` (weighted)	numerical.py:726-740	HIGH	Uncompensated weighted moments + bug at line 740
8	`rolling_moving_average`	numerical.py:1065	MED	Initial-window `np.sum` uncompensated; Kahan only on rolling
9	`compute_moments_slope_mi` (MAD)	numerical.py:717, 749	LOW-MED	Uncompensated abs-deviation sum
10	`compute_numerical_aggregates_numba` (drawdown)	numerical.py:442-497	MED	Recursive call inherits all upstream instabilities

Bug found and fixed on the spot¶

# numerical.py:738-740 (BEFORE)
kurt += summand * d
if weights is not None:
    weighted_skew += w_summand * w_d * next_weight  # WRONG accumulator!

# numerical.py:738-740 (AFTER, fixed 2026-04-24)
kurt += summand * d
if weights is not None:
    weighted_kurt += w_summand * w_d * next_weight  # accumulator name corrected

weighted_skew was being doubly accumulated (line 736 + line 740); weighted_kurt was never accumulated, so the downstream normaliser at line 771 always divided 0 by factor and emitted -3.0 for every weighted call. This silently produced wrong feature values for every caller passing weights to compute_moments_slope_mi.

Stable-kernel implementations¶

Added mlframe/feature_engineering/_numerical_stable.py with:

welford_mean_var_seq(arr) — single-pass Welford for mean+var
welford_moments_seq(arr) — Pébay generalised online for mean/var/skew/kurt
kahan_sum_seq(arr) — Neumaier-improved Kahan compensated sum
kahan_dot_seq(a, b) — Kahan-compensated dot product (slope/correlation numerator)
kahan_two_pass_var_seq(arr) — best-of-both: exact Kahan mean (pass 1) + Kahan sum-of-squared-deviations (pass 2)

Plus naive baselines (matching the patterns in numerical.py) for direct comparison.

Benchmark results¶

Inputs (sized to expose cancellation; numpy float64):

Name	Distribution	N	Notes
`normal_N1k`	N(0, 1)	1,000	easy baseline
`normal_N100k`	N(0, 1)	100,000	drift accumulation
`normal_N1M`	N(0, 1)	1,000,000	longer drift
`large_mean_small_var`	1e6 + N(0, 0.1)	100,000	mild cancellation
`small_variance_1e9`	1e9 + N(0, 1e-3)	100,000	EXTREME cancellation
`lognormal`	exp(N(0, 1))	100,000	wide range
`high_kurtosis`	90% N(0,1) + 10% N(0,10)	100,000	heavy tails

MEAN / VARIANCE — relative error vs `np.var(ddof=0)`¶

distribution	naive 2-pass	Welford	Kahan 2-pass	best
normal_N1k	2.3e-16	5.7e-16	0.0	kahan
normal_N100k	4.2e-15	2.2e-16	2.2e-16	kahan
normal_N1M	4.2e-14	1.1e-15	0.0	kahan
large_mean_small_var	6.8e-15	1.0e-10	0.0	kahan
small_variance_1e9	3.6e-5	5.5e-6	0.0	kahan
lognormal	1.8e-15	1.6e-15	0.0	kahan
high_kurtosis	9.8e-16	2.1e-15	3.3e-16	kahan

Key observation: Kahan-2pass gives literal zero error (numerical exactness at float64) on 6 of 7 distributions. The catastrophic case small_variance_1e9 — where naive loses 5 digits of precision (3.6e-5 relative error means only 4-5 correct digits out of 16 possible) — Kahan recovers full precision.

MOMENTS (skew + kurt) — relative error vs `scipy.stats.skew/kurtosis(bias=True)`¶

distribution	naive_skew	welf_skew	naive_kurt	welf_kurt
normal_N1k	2.5e-15	6.4e-16	8.8e-14	7.3e-14
normal_N100k	4.4e-15	1.0e-14	2.9e-12	2.2e-13
normal_N1M	1.4e-13	2.9e-13	1.3e-11	7.4e-13
large_mean_small_var	3.9e-5	2.3e-7	1.4e-6	3.3e-7
small_variance_1e9	2.78 ⚠️	8.8e-3	8.8e-3	2.2e-3
lognormal	5.5e-15	6.8e-15	2.4e-15	1.9e-14
high_kurtosis	0.0	2.5e-15	1.6e-15	1.1e-14

Critical finding — on small_variance_1e9: - Naive returns skew = 2.78 while reference is ~0.011 — total garbage (sign and magnitude both wrong) - Welford-Pébay returns 8.8e-3 — recovers 2-3 correct digits - This is exactly the regime where data has high baseline + small noise (sensor data, financial prices, etc.) — common in feature-engineering

On easy cases (long arrays of normal data) Welford is slightly worse — per-element rounding accumulates over 100k-1M iterations. This is a legitimate trade-off: prefer Welford when input MIGHT be ill-conditioned; prefer naive 2-pass when input is reliably well-conditioned. Hybrid: detect via cheap pre-check (e.g., var/mean^2 < 1e-12 → switch to Kahan).

KAHAN SUM — relative error vs `np.sum(arr)` (numpy pairwise)¶

distribution	naive_sum	kahan_sum	improvement
normal_N1k	8.6e-16	3.7e-16	2.3×
normal_N100k	3.0e-15	0.0	∞ (exact)
normal_N1M	3.9e-15	1.5e-16	26×
large_mean_small_var	1.3e-14	0.0	∞
small_variance_1e9	5.9e-15	0.0	∞
lognormal	8.8e-15	0.0	∞
high_kurtosis	2.6e-14	3.4e-16	75×

Kahan recovers 1-2 orders of magnitude vs uncompensated for-loop sum (numpy's pairwise sum is already partially compensated, hence baselines ~1e-15; Python for x: s += x is much worse at ~1e-13 on N=1M).

RUNTIME (Win32 Anaconda 3.11, numba JIT warm)¶

Kernel	N	time/call	vs naive
naive_mean_var_2pass	100,000	0.40 ms	1.00×
welford_mean_var	100,000	0.72 ms	1.81×
naive_moments_2pass	100,000	0.40 ms	1.00×
welford_moments	100,000	0.83 ms	2.11×
naive_mean_var_2pass	1,000,000	3.68 ms	1.00×
welford_mean_var	1,000,000	6.95 ms	1.89×
naive_moments_2pass	1,000,000	3.21 ms	1.00×
welford_moments	1,000,000	8.16 ms	2.54×

Verdict: 1.8-2.5× wallclock penalty for Welford / Kahan kernels. Acceptable for hot-path features given the precision wins on ill-conditioned inputs. Kahan-2pass for variance: similar 2× cost, better precision than Welford in mean+var case.

Recommendations¶

Shipping status (verified against _numerical_numba.py / numerical.py). Phase 1 and Phase 2 are ✅ DONE: the weighted-moments bug fix, Kahan-compensated inner variance sum (_make_compute_simple_stats(use_kahan=...)), Welford/compensated skew-kurt (_make_compute_moments_slope_mi(use_kahan, ...), compensated= param), and Kahan slope_over/under/r_sum accumulators all landed and are covered by tests/feature_engineering/test_numerical_stability_bench.py. Still open (Phase 3): items 5 (geomean always-log-mode — code still uses the multiply-then-switch heuristic), 6 (rolling initial-window still uses uncompensated np.sum(arr[:n])), and the var/mean^2 < 1e-12 cheap auto-detect switch (routing is via the explicit compensated= flag, not the automatic heuristic). Those three are the live TODO.

Phase 1 — quick wins ✅ DONE¶

Fix line 740 bug ✅ DONE

Replace inner variance sum with Kahan compensated ✅ DONE (_make_compute_simple_stats(use_kahan=...)):

# Old: std_val = std_val + summand
# New: t = std_val + summand
#      c += (std_val - t) + summand if abs(std_val) >= abs(summand) else (summand - t) + std_val
#      std_val = t
# Final: return ..., np.sqrt((std_val + c) / size)

Per-line-overhead small, recovery 1-2 orders of magnitude.

Phase 2 — moments stability ✅ DONE¶

Replace compute_moments_slope_mi skew/kurt accumulators with Welford-Pébay (welford_moments_seq is the reference impl). Acceptable runtime cost (~2x), eliminates the catastrophic-cancellation regime where naive produces wrong-sign skewness.
Replace slope_over / slope_under / r_sum accumulators with Kahan compensated sums. Cheap (~1.2× cost), critical for regression-feature correctness on long time series.

Phase 3 — geomean / harmonic / drawdown (OPEN — not started)¶

Replace heuristic geomean log-mode (line 340-354) with always-log-mode for non-negative input + Kahan log-sum:
```
geometric_mean = exp(kahan_sum_seq(np.log(np.maximum(arr, 1e-300))))
```
Initial-window sum in rolling_moving_average (line 1065) → Kahan
Drawdown stats (lines 442-497) — fix upstream first, drawdowns inherit improvements

Hot-path heuristic¶

For features called per-row on large arrays, default to Welford (safe across all input regimes). For one-shot dataset-level statistics, Kahan-2pass is faster and more precise.

The same class of bug this report documents for numerical.py -- a raw-power-sum / textbook-binomial-expansion skew/kurt formula that catastrophically cancels once values sit at a large offset relative to their spread -- has independently turned up in three other modules, all fixed the same way (two-pass centred-moment accumulation instead of algebraic expansion of raw power sums):

_binned_numeric_agg_fe.py::_global_stats_all (global fallback stats) -- fixed.
_target_encoding_fe.py::_raw_moment_sums (per-category skew/kurt) -- fixed.
_binned_numeric_agg_fe.py::_derive_cell_stats (per-cell skew/kurt) -- still uses the unstable formula; open follow-up.

Any new skew/kurt kernel should default to the two-pass centred-moment form (or Welford-Pébay, per the hot-path heuristic above) rather than expanding E[(x-mean)^k] algebraically in terms of raw power sums.

Files¶

mlframe/feature_engineering/_numerical_stable.py — new reference impls
tests/feature_engineering/test_numerical_stability_bench.py — runnable benchmarks
mlframe/feature_engineering/_numerical_numba.py:646-657 — bug fix (moved from the original numerical.py:738-749 location)

Numerical-Stability Report: feature_engineering/numerical.py Moments¶