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 <noreply@anthropic.com>

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::<f64>().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<Vec<f64>>, mut y: Vec<f64>, 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<Vec<f64>>, mut y: Vec<f64>, 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<Vec<f64>>, mut y: Vec<f64>, xm: &[f64]) -> QrFit {
     let mut kept: Vec<usize> = 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<f64>], k: usize, constant: bool) -> R<Fit> {
         (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<f64> = 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<f64> = (0..k).map(|j| norm2(&x.iter().map(|r| r.get(j).copied().unwrap_or(0.0)).collect::<Vec<_>>())).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::<f64>() } else { 0.0 };
     let used = dropped.iter().filter(|d| !**d).count() + usize::from(constant);
     let df = nf - used as f64;

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<f64> = (1..=n).map(|i| f64::from(i)).collect();
+            let ys: Vec<f64> = (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<f64> = 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);