From 583c8a4819bd8f1a94431d041203d3e1c53e23b9 Mon Sep 17 00:00:00 2001 From: vangel76 Date: Sat, 7 Nov 2026 17:01:20 +0100 Subject: [PATCH] LINEST: keep predictors that vary by little next to a large offset A centred column was dropped when its length fell under 1e-12 of the raw column, which also removed real variation such as [1e13, 1e13+1, 1e13+2]. Constant columns are now recognised on the raw values, within a few ulps. Co-Authored-By: Claude Sonnet 5.5 --- crates/functions/src/stat/regression.rs | 29 ++++++++++++++++--------- crates/functions/src/stat/tests.rs | 4 ++++ 2 files changed, 23 insertions(+), 10 deletions(-) diff --git a/crates/functions/src/stat/regression.rs b/crates/functions/src/stat/regression.rs index f5eccfff9..fe9515860 100644 --- a/crates/functions/src/stat/regression.rs +++ b/crates/functions/src/stat/regression.rs @@ -130,8 +130,9 @@ fn norm2(v: &[f64]) -> f64 { m * v.iter().map(|x| (x / m) * (x / m)).sum::().sqrt() } -/// Relative size below which a centred column is taken to be rounding error (about 4500 ulps). -const NOISE: f64 = 1e-12; +/// A column whose values differ by no more than this many ulps of its largest magnitude is +/// taken to be constant: the differences are rounding in how it was computed, not data. +const CONSTANT_ULPS: f64 = 64.0; /// Result of [`qr_solve`]: the least-squares solution on the kept columns, in the caller's /// column scale, and what the standard errors need. @@ -152,10 +153,12 @@ struct QrFit { /// pivot of the earlier normal-equation solver, on squares) is collinear: it is dropped, with /// coefficient and standard error 0, as Excel does. /// -/// `raw_len[j]` is the length of column `j` before centring. A centred column shorter than -/// [`NOISE`] of that is rounding residue of a constant column (`[0.1; 3]` centres to about -/// -1.4e-17, not 0): collinear with the intercept, so dropped, not scaled up to unit length. -fn qr_solve(mut cols: Vec>, mut y: Vec, xm: &[f64], raw_len: &[f64]) -> QrFit { +/// `constant[j]` says column `j` was constant before centring. Centring it leaves rounding +/// residue, not 0 (`[0.1; 3]` centres to about -1.4e-17): it is collinear with the intercept, so +/// it is dropped instead of being scaled up to unit length. Whether a column is constant is +/// decided on its raw values, not by comparing the centred length with the offset: a predictor +/// like `[1e13, 1e13 + 1, 1e13 + 2]` is far smaller than its offset and still varies. +fn qr_solve(mut cols: Vec>, mut y: Vec, xm: &[f64], constant: &[bool]) -> QrFit { let n = y.len(); let k = cols.len(); let mut scale = vec![1.0; k]; @@ -166,8 +169,7 @@ fn qr_solve(mut cols: Vec>, mut y: Vec, xm: &[f64], raw_len: &[f64 let mut kept: Vec = Vec::new(); for (j, col) in cols.iter_mut().enumerate() { let len = norm2(col); - let raw = raw_len.get(j).copied().unwrap_or(len); - if !(len > 1e-300 && len.is_finite() && len > NOISE * raw) { + if constant.get(j).copied().unwrap_or(false) || !(len > 1e-300 && len.is_finite()) { dropped[j] = true; continue; } @@ -273,8 +275,15 @@ fn regress(y: &[f64], x: &[Vec], k: usize, constant: bool) -> R { (0..k).map(|j| x.iter().map(|r| r.get(j).copied().unwrap_or(0.0) - xm.get(j).copied().unwrap_or(0.0)).collect()).collect(); let yc: Vec = y.iter().map(|v| v - ym).collect(); let sstot: f64 = if constant { yc.iter().map(|v| v * v).sum() } else { y.iter().map(|v| v * v).sum() }; - let raw_len: Vec = (0..k).map(|j| norm2(&x.iter().map(|r| r.get(j).copied().unwrap_or(0.0)).collect::>())).collect(); - let QrFit { coef, dropped, ssresid, diag, quad } = qr_solve(cols, yc, &xm, &raw_len); + // Only centring leaves rounding residue behind; without an intercept a constant column is data. + let constant_cols: Vec = (0..k) + .map(|j| { + let col = x.iter().map(|r| r.get(j).copied().unwrap_or(0.0)); + let (lo, hi, big) = col.fold((f64::INFINITY, f64::NEG_INFINITY, 0.0f64), |(lo, hi, big), v| (lo.min(v), hi.max(v), big.max(v.abs()))); + constant && hi - lo <= CONSTANT_ULPS * f64::EPSILON * big + }) + .collect(); + let QrFit { coef, dropped, ssresid, diag, quad } = qr_solve(cols, yc, &xm, &constant_cols); let b = if constant { ym - coef.iter().zip(&xm).map(|(c, m)| c * m).sum::() } else { 0.0 }; let used = dropped.iter().filter(|d| !**d).count() + usize::from(constant); let df = nf - used as f64; diff --git a/crates/functions/src/stat/tests.rs b/crates/functions/src/stat/tests.rs index 1dd208111..29ca2d987 100644 --- a/crates/functions/src/stat/tests.rs +++ b/crates/functions/src/stat/tests.rs @@ -456,6 +456,10 @@ fn linest_drops_collinear_columns() { let r = linest_rows(&y, &[off], true); rel_close(r[0][0], 0.6, 1e-6, "slope on an offset column"); assert_eq!(r[3][1], Some(3.0)); + // Distinct values far from zero, whose spread is below 1e-12 of the offset, still vary. + let r = linest_rows(&[0.0, 1.0, 2.0], &[vec![1e13, 1e13 + 1.0, 1e13 + 2.0]], true); + rel_close(r[0][0], 1.0, 1e-9, "slope on a predictor with a 1e13 offset"); + assert_eq!(r[3][1], Some(1.0), "df: the predictor is kept"); // More predictors than observations: the surplus columns drop, no panic. let r = linest_rows(&[1.0, 2.0], &[vec![1.0, 2.0], vec![2.0, 1.0], vec![0.0, 5.0]], true); assert_eq!(r[0].len(), 4);