From d0b40c3b3029adae0a600aa3f874c154d5e0eaee Mon Sep 17 00:00:00 2001 From: vangel76 Date: Wed, 4 Nov 2026 00:05:38 +0100 Subject: [PATCH] linest: drop a constant column that centres to rounding noise (review point) CodeRabbit pointed out that the drop test used the centred column length, so a constant column such as [0.1; 3] (centred to about -1.4e-17) was scaled to a unit vector and kept, giving a meaningless coefficient and wrong df. A centred column shorter than 1e-12 of its uncentred length is now treated as collinear with the intercept. The suggested 1e-6 threshold was not used: it would also drop columns that sit far from zero, such as 1e9 + i. Co-Authored-By: Claude Sonnet 5.5 --- crates/functions/src/stat/regression.rs | 15 ++++++++++++--- crates/functions/src/stat/tests.rs | 18 ++++++++++++++++++ 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/crates/functions/src/stat/regression.rs b/crates/functions/src/stat/regression.rs index 3da3ee644..f5eccfff9 100644 --- a/crates/functions/src/stat/regression.rs +++ b/crates/functions/src/stat/regression.rs @@ -130,6 +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; + /// Result of [`qr_solve`]: the least-squares solution on the kept columns, in the caller's /// column scale, and what the standard errors need. struct QrFit { @@ -148,7 +151,11 @@ struct QrFit { /// outside the span of the kept columns is under 1e-6 of its own length (the 1e-12 relative /// pivot of the earlier normal-equation solver, on squares) is collinear: it is dropped, with /// coefficient and standard error 0, as Excel does. -fn qr_solve(mut cols: Vec>, mut y: Vec, xm: &[f64]) -> QrFit { +/// +/// `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 { let n = y.len(); let k = cols.len(); let mut scale = vec![1.0; k]; @@ -159,7 +166,8 @@ fn qr_solve(mut cols: Vec>, mut y: Vec, xm: &[f64]) -> QrFit { let mut kept: Vec = Vec::new(); for (j, col) in cols.iter_mut().enumerate() { let len = norm2(col); - if !(len > 1e-300 && len.is_finite()) { + let raw = raw_len.get(j).copied().unwrap_or(len); + if !(len > 1e-300 && len.is_finite() && len > NOISE * raw) { dropped[j] = true; continue; } @@ -265,7 +273,8 @@ 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 QrFit { coef, dropped, ssresid, diag, quad } = qr_solve(cols, yc, &xm); + 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); 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 e4c3082fe..ccbbf8a3a 100644 --- a/crates/functions/src/stat/tests.rs +++ b/crates/functions/src/stat/tests.rs @@ -438,6 +438,24 @@ fn linest_drops_collinear_columns() { assert_eq!(r[1][1], Some(0.0)); rel_close(r[0][0], 0.6, 1e-12, "slope"); rel_close(r[0][2], 2.2, 1e-12, "intercept"); + // A constant column whose mean does not centre to exactly 0 (the mean of three 0.1 is + // 0.10000000000000002) is just as collinear with the intercept: dropped, not scaled up to a + // unit vector. Some length hits a rounding mean for each constant. + for c in [0.1, 1.1, 1e6 + 0.3] { + for n in 3..=12 { + let xs: Vec = (1..=n).map(|i| f64::from(i)).collect(); + let ys: Vec = (1..=n).map(|i| f64::from(i * i % 7) + 0.5 * f64::from(i)).collect(); + let r = linest_rows(&ys, &[vec![c; n as usize], xs], true); + assert_eq!(r[0][1], Some(0.0), "coefficient of the constant {c}, n={n}"); + assert_eq!(r[1][1], Some(0.0), "standard error of the constant {c}, n={n}"); + assert_eq!(r[3][1], Some(f64::from(n) - 2.0), "df with the constant {c}, n={n}"); + } + } + // A column that merely sits far from zero is not noise and is kept. + let off: Vec = x1.iter().map(|v| 1e9 + v).collect(); + 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)); // 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);