From f74d95542fb8f1e1b6c87c3a0a5957ed20c46e56 Mon Sep 17 00:00:00 2001 From: vangel76 Date: Tue, 3 Nov 2026 23:38:55 +0100 Subject: [PATCH 1/5] LINEST: solve by Householder QR instead of the normal equations (#3322) The normal equations square the condition number, so a degree-5 polynomial on x = 1..12 came back with errors up to 2.6e-6 and x, x^2 near 1e6 gave slope 2.99976 and intercept 251 for y = 3x + 7. LINEST, LOGEST, TREND and GROWTH now factor the centred, column-scaled data with Householder QR. Collinear columns are still dropped under the same relative tolerance. Co-Authored-By: Claude Sonnet 5.5 --- crates/functions/src/stat/regression.rs | 222 ++++++++++++++---------- crates/functions/src/stat/tests.rs | 148 ++++++++++++++++ 2 files changed, 278 insertions(+), 92 deletions(-) diff --git a/crates/functions/src/stat/regression.rs b/crates/functions/src/stat/regression.rs index 55bf0aec7..3da3ee644 100644 --- a/crates/functions/src/stat/regression.rs +++ b/crates/functions/src/stat/regression.rs @@ -55,12 +55,12 @@ fn strict_num(v: &Value) -> R { } } -/// Most predictor columns LINEST, LOGEST, TREND and GROWTH accept. The normal equations are a -/// k×k matrix inverted in O(k³), so an unbounded k (one `SEQUENCE(1,16384)` argument) allocated +/// Most predictor columns LINEST, LOGEST, TREND and GROWTH accept. The fit works on a k×k +/// triangle in O(k³), so an unbounded k (one `SEQUENCE(1,16384)` argument) allocated /// gigabytes and could abort the process. pub(super) const MAX_REGRESSORS: usize = 255; -/// Most multiply-adds (rows × k²) a regression may spend building its normal equations. +/// Most multiply-adds (rows × k²) a regression may spend on its factorisation. const MAX_REGRESSION_WORK: usize = 1_000_000_000; /// Rejects regressions too large to fit: more than [`MAX_REGRESSORS`] predictors is `#VALUE!`, @@ -121,59 +121,130 @@ struct Fit { ssresid: f64, } -/// Inverts a symmetric positive semi-definite matrix by Gauss–Jordan on the diagonal, dropping -/// (zeroing) collinear variables. Returns the inverse and the dropped flags. -fn sym_inverse(m: &[Vec]) -> (Vec>, Vec) { - let p = m.len(); - let mut a: Vec> = m.to_vec(); - let mut inv: Vec> = (0..p).map(|i| (0..p).map(|j| if i == j { 1.0 } else { 0.0 }).collect()).collect(); - let mut dropped = vec![false; p]; - let diag: Vec = (0..p).map(|i| m.get(i).and_then(|r| r.get(i)).copied().unwrap_or(0.0).abs()).collect(); - for j in 0..p { - let piv = a.get(j).and_then(|r| r.get(j)).copied().unwrap_or(0.0); - let scale = diag.get(j).copied().unwrap_or(0.0); - if piv.abs() <= 1e-12 * scale || piv.abs() < 1e-300 { - if let Some(d) = dropped.get_mut(j) { - *d = true; - } - for row in a.iter_mut() { - if let Some(v) = row.get_mut(j) { - *v = 0.0; - } - } - if let Some(r) = a.get_mut(j) { - r.iter_mut().for_each(|v| *v = 0.0); - } - if let Some(r) = inv.get_mut(j) { - r.iter_mut().for_each(|v| *v = 0.0); - } +/// Euclidean norm without intermediate overflow or underflow. +fn norm2(v: &[f64]) -> f64 { + let m = v.iter().fold(0.0f64, |m, x| m.max(x.abs())); + if m == 0.0 || !m.is_finite() { + return m; + } + m * v.iter().map(|x| (x / m) * (x / m)).sum::().sqrt() +} + +/// 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 { + coef: Vec, + dropped: Vec, + /// Residual sum of squares of the fit. + ssresid: f64, + /// Diagonal of `(X'X)^-1` over the kept columns (0 for dropped ones). + diag: Vec, + /// `xm' (X'X)^-1 xm` over the kept columns. + quad: f64, +} + +/// Least squares by Householder QR on unit-norm columns, instead of the normal equations, whose +/// condition number is the square of the data's. Columns are taken left to right; one whose part +/// 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 { + let n = y.len(); + let k = cols.len(); + let mut scale = vec![1.0; k]; + let mut dropped = vec![false; k]; + // Householder vectors (unit norm, acting on rows `at..n`) and the columns of R. + let mut refl: Vec<(usize, Vec)> = Vec::new(); + let mut r_cols: Vec> = Vec::new(); + 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()) { + dropped[j] = true; continue; } - let arow: Vec = a.get(j).map(|r| r.iter().map(|v| v / piv).collect()).unwrap_or_default(); - let irow: Vec = inv.get(j).map(|r| r.iter().map(|v| v / piv).collect()).unwrap_or_default(); - for i in 0..p { - if i == j { - continue; - } - let f = a.get(i).and_then(|r| r.get(j)).copied().unwrap_or(0.0); - if f == 0.0 { - continue; - } - if let Some(r) = a.get_mut(i) { - r.iter_mut().zip(&arow).for_each(|(v, w)| *v -= f * w); - } - if let Some(r) = inv.get_mut(i) { - r.iter_mut().zip(&irow).for_each(|(v, w)| *v -= f * w); - } + col.iter_mut().for_each(|v| *v /= len); + scale[j] = len; + for (at, v) in &refl { + apply_reflector(*at, v, col); } - if let Some(r) = a.get_mut(j) { - *r = arow; + let at = refl.len(); + let tail = col.get(at..).unwrap_or(&[]); + let alpha = norm2(tail); + if alpha * alpha <= 1e-12 { + dropped[j] = true; + continue; } - if let Some(r) = inv.get_mut(j) { - *r = irow; + let sign = if tail.first().is_some_and(|t| *t < 0.0) { -1.0 } else { 1.0 }; + let mut v: Vec = tail.to_vec(); + if let Some(h) = v.first_mut() { + *h += sign * alpha; + } + let vn = norm2(&v); + v.iter_mut().for_each(|x| *x /= vn); + let mut rc: Vec = col.get(..at).unwrap_or(&[]).to_vec(); + rc.push(-sign * alpha); + r_cols.push(rc); + refl.push((at, v)); + kept.push(j); + } + for (at, v) in &refl { + apply_reflector(*at, v, &mut y); + } + let r = kept.len(); + // Back substitution R c = (Q'y)[..r], then undo the column scaling. + let mut c = vec![0.0; r]; + for i in (0..r).rev() { + let mut t = y.get(i).copied().unwrap_or(0.0); + for jj in i + 1..r { + t -= r_cols.get(jj).and_then(|rc| rc.get(i)).copied().unwrap_or(0.0) * c.get(jj).copied().unwrap_or(0.0); + } + if let Some(ci) = c.get_mut(i) { + *ci = t / r_cols.get(i).and_then(|rc| rc.get(i)).copied().unwrap_or(f64::NAN); } } - (inv, dropped) + let ssresid = y.get(r..n).unwrap_or(&[]).iter().map(|v| v * v).sum::(); + // R^-1 (upper triangular), row i of the inverse at rinv[i]. + let rij = |i: usize, j: usize| r_cols.get(j).and_then(|rc| rc.get(i)).copied().unwrap_or(0.0); + let mut rinv = vec![vec![0.0; r]; r]; + for t in 0..r { + for i in (0..=t).rev() { + let mut v = if i == t { 1.0 } else { 0.0 }; + for jj in i + 1..=t { + v -= rij(i, jj) * rinv.get(jj).and_then(|row| row.get(t)).copied().unwrap_or(0.0); + } + if let Some(cell) = rinv.get_mut(i).and_then(|row| row.get_mut(t)) { + *cell = v / rij(i, i); + } + } + } + let mut coef = vec![0.0; k]; + let mut diag = vec![0.0; k]; + for (i, &j) in kept.iter().enumerate() { + let s = scale.get(j).copied().unwrap_or(1.0); + if let Some(v) = coef.get_mut(j) { + *v = c.get(i).copied().unwrap_or(0.0) / s; + } + let row = rinv.get(i); + if let Some(v) = diag.get_mut(j) { + *v = row.map_or(0.0, |rw| rw.iter().map(|x| x * x).sum::()) / (s * s); + } + } + // xm' (S R'R S)^-1 xm = |R^-T S^-1 xm|². + let w: Vec = kept.iter().map(|&j| xm.get(j).copied().unwrap_or(0.0) / scale.get(j).copied().unwrap_or(1.0)).collect(); + let quad = (0..r) + .map(|t| { + (0..=t).map(|i| rinv.get(i).and_then(|row| row.get(t)).copied().unwrap_or(0.0) * w.get(i).copied().unwrap_or(0.0)).sum::().powi(2) + }) + .sum::(); + QrFit { coef, dropped, ssresid, diag, quad } +} + +/// Applies the Householder reflection `I - 2vv'` on rows `at..` of `x`. +fn apply_reflector(at: usize, v: &[f64], x: &mut [f64]) { + let Some(tail) = x.get_mut(at..) else { return }; + let dot: f64 = v.iter().zip(tail.iter()).map(|(a, b)| a * b).sum(); + tail.iter_mut().zip(v).for_each(|(t, a)| *t -= 2.0 * dot * a); } fn regress(y: &[f64], x: &[Vec], k: usize, constant: bool) -> R { @@ -189,35 +260,13 @@ fn regress(y: &[f64], x: &[Vec], k: usize, constant: bool) -> R { } else { (vec![0.0; k], 0.0) }; - // Normal equations on (centred) data. - let mut xtx = vec![vec![0.0; k]; k]; - let mut xty = vec![0.0; k]; - for (row, yv) in x.iter().zip(y) { - let d: Vec = (0..k).map(|j| row.get(j).copied().unwrap_or(0.0) - xm.get(j).copied().unwrap_or(0.0)).collect(); - let yd = yv - ym; - for (i, di) in d.iter().enumerate() { - if let Some(t) = xty.get_mut(i) { - *t += di * yd; - } - if let Some(r) = xtx.get_mut(i) { - for (j, dj) in d.iter().enumerate() { - if let Some(v) = r.get_mut(j) { - *v += di * dj; - } - } - } - } - } - let (inv, dropped) = sym_inverse(&xtx); - let coef: Vec = (0..k).map(|i| inv.get(i).map(|r| r.iter().zip(&xty).map(|(a, b)| a * b).sum()).unwrap_or(0.0)).collect(); + // (Centred) predictor columns and response. + let cols: Vec> = + (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 b = if constant { ym - coef.iter().zip(&xm).map(|(c, m)| c * m).sum::() } else { 0.0 }; - let mut ssresid = 0.0; - let mut sstot = 0.0; - for (row, yv) in x.iter().zip(y) { - let pred = b + coef.iter().zip(row).map(|(c, v)| c * v).sum::(); - ssresid += (yv - pred).powi(2); - sstot += if constant { (yv - ym).powi(2) } else { yv * yv }; - } let used = dropped.iter().filter(|d| !**d).count() + usize::from(constant); let df = nf - used as f64; let ssreg = sstot - ssresid; @@ -230,20 +279,9 @@ fn regress(y: &[f64], x: &[Vec], k: usize, constant: bool) -> R { _ => Err(CellError::Num), }; let s2 = sey2.unwrap_or(f64::NAN); - let se: Vec = (0..k) - .map(|i| if dropped.get(i).copied().unwrap_or(true) { 0.0 } else { (inv.get(i).and_then(|r| r.get(i)).copied().unwrap_or(0.0) * s2).sqrt() }) - .collect(); - let se_b = if constant { - let mut q = 0.0; - for (i, mi) in xm.iter().enumerate() { - for (j, mj) in xm.iter().enumerate() { - q += mi * mj * inv.get(i).and_then(|r| r.get(j)).copied().unwrap_or(0.0); - } - } - Some((s2 * (1.0 / nf + q)).sqrt()) - } else { - None - }; + let se: Vec = + (0..k).map(|i| if dropped.get(i).copied().unwrap_or(true) { 0.0 } else { (diag.get(i).copied().unwrap_or(0.0) * s2).sqrt() }).collect(); + let se_b = if constant { Some((s2 * (1.0 / nf + quad)).sqrt()) } else { None }; Ok(Fit { coef, b, se, se_b, r2, sey, f, df, ssreg, ssresid }) } diff --git a/crates/functions/src/stat/tests.rs b/crates/functions/src/stat/tests.rs index ba59116cd..e4c3082fe 100644 --- a/crates/functions/src/stat/tests.rs +++ b/crates/functions/src/stat/tests.rs @@ -311,6 +311,154 @@ fn linest_family() { is_err(ev("LINEST", vec![c(&[1.0, 2.0]), c(&[1.0, 2.0, 3.0])]), CellError::Ref); } +/// LINEST on columns of predictors, `stats` on, as rows of numbers (`None` for an error cell). +fn linest_rows(y: &[f64], xs: &[Vec], constant: bool) -> Vec>> { + let x = arr((0..y.len()).map(|i| xs.iter().map(|col| nv(col[i])).collect()).collect()); + let r = rows_of(&ev("LINEST", vec![c(y), av(x), b(constant), b(true)])); + r.iter().map(|row| row.iter().map(Value::as_f64).collect()).collect() +} + +fn rel_close(got: Option, want: f64, tol: f64, what: &str) { + let got = got.unwrap_or(f64::NAN); + assert!((got - want).abs() <= tol * want.abs().max(1e-300), "{what}: got {got:e}, want {want:e}, rel tol {tol:e}"); +} + +// B2: LINEST solved the normal equations, whose condition number is the square of the data's. +// Expected values below come from exact rational arithmetic: the normal equations solved with +// Python `fractions.Fraction` (no rounding at all), rounded once to f64 (script in the PR text). +#[test] +fn linest_ill_conditioned_polynomial_is_exact() { + // y = 1 + x + x² + x³ + x⁴ + x⁵ on x = 1..12: every coefficient is exactly 1. The normal + // equations were off by up to 2.6e-6. + let xs: Vec = (1..=12).map(f64::from).collect(); + let ys: Vec = xs.iter().map(|x| (0..=5).map(|j| x.powi(j)).sum()).collect(); + let cols: Vec> = (1..=5).map(|j| xs.iter().map(|x| x.powi(j)).collect()).collect(); + let r = linest_rows(&ys, &cols, true); + for (i, v) in r[0].iter().enumerate() { + rel_close(*v, 1.0, 1e-9, &format!("coefficient {i}")); + } +} + +#[test] +fn linest_large_offset_regressors_recover_the_line() { + // y = 3x + 7 exactly with x = 1e6 + i and regressors [x, x²]: the true fit has x² coefficient + // 0, x coefficient 3, intercept 7. The normal equations returned 2.99976 and 251.14. + let xs: Vec = (0..12).map(|i| 1_000_000.0 + f64::from(i)).collect(); + let ys: Vec = xs.iter().map(|x| 3.0 * x + 7.0).collect(); + let r = linest_rows(&ys, &[xs.clone(), xs.iter().map(|x| x * x).collect()], true); + assert!(r[0][0].unwrap().abs() < 1e-12, "x² coefficient {:?}", r[0][0]); + rel_close(r[0][1], 3.0, 1e-9, "x coefficient"); + // The intercept subtracts 3e6 from 3e6 + 7 and carries the x² coefficient times 1e12, so + // it is good to about 1e-4 in absolute terms, not better. + assert!((r[0][2].unwrap() - 7.0).abs() < 1e-3, "intercept {:?}", r[0][2]); + // One regressor: exact. + let one = linest_rows(&ys, std::slice::from_ref(&xs), true); + rel_close(one[0][0], 3.0, 1e-12, "slope"); + rel_close(one[0][1], 7.0, 1e-9, "intercept"); +} + +#[test] +fn linest_noisy_fits_match_exact_rational_least_squares() { + // Two predictors with noise, constant fitted. Order of each row: x2, x1, intercept. + let x1 = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]; + let x2 = [2.0, 1.0, 4.0, 3.0, 7.0, 5.0, 9.0, 6.0]; + let y = [5.0, 6.0, 11.0, 12.0, 20.0, 19.0, 30.0, 28.0]; + let r = linest_rows(&y, &[x1.to_vec(), x2.to_vec()], true); + let (coef, se) = ([1.4270152505446623, 2.452069716775599, -1.2592592592592593], [0.3376442444771937, 0.3679392851360571, 1.0847489958915801]); + for i in 0..3 { + rel_close(r[0][i], coef[i], 1e-12, "coef"); + rel_close(r[1][i], se[i], 1e-12, "se"); + } + rel_close(r[2][0], 0.985070139277778, 1e-12, "r2"); + rel_close(r[2][1], 1.367057173604725, 1e-12, "sey"); + rel_close(r[3][0], 164.94965318256004, 1e-11, "F"); + assert_eq!(r[3][1], Some(5.0)); + rel_close(r[4][0], 616.5307734204792, 1e-12, "ssreg"); + rel_close(r[4][1], 9.344226579520697, 1e-11, "ssresid"); + assert_eq!(r[2][2], None); + // const = FALSE: through the origin, R² and F from the uncentred sums. + let r = linest_rows(&y, &[x1.to_vec(), x2.to_vec()], false); + rel_close(r[0][0], 1.3529411764705883, 1e-12, "coef"); + rel_close(r[0][1], 2.303921568627451, 1e-12, "coef"); + assert_eq!(r[0][2], Some(0.0)); + rel_close(r[1][0], 0.34102994374364876, 1e-12, "se"); + rel_close(r[1][1], 0.35495521934550084, 1e-12, "se"); + assert_eq!(r[1][2], None); + rel_close(r[2][0], 0.995718966041848, 1e-12, "r2"); + rel_close(r[2][1], 1.4061024795535124, 1e-12, "sey"); + rel_close(r[3][0], 697.7652892561983, 1e-11, "F"); + assert_eq!(r[3][1], Some(6.0)); + rel_close(r[4][0], 2759.1372549019607, 1e-12, "ssreg"); + rel_close(r[4][1], 11.862745098039216, 1e-11, "ssresid"); + // Cubic and quintic with noise on x = 1..12; order x^k .. x, intercept. + let xs: Vec = (1..=12).map(f64::from).collect(); + let noise = [3.0, -2.0, 4.0, -1.0, 2.0, -3.0, 1.0, -4.0, 2.0, 3.0, -2.0, 1.0]; + let cub: Vec = xs.iter().zip(noise).map(|(x, e)| 1.0 + x + x * x + x.powi(3) + e).collect(); + let r = linest_rows(&cub, &(1..=3).map(|j| xs.iter().map(|x| x.powi(j)).collect()).collect::>>(), true); + let coef = [0.9955969955969955, 1.1452991452991452, -0.33514633514633513, 4.373737373737374]; + let se = [0.02755463331301788, 0.5434124199147191, 3.102549978874372, 4.851885058051787]; + for i in 0..4 { + rel_close(r[0][i], coef[i], 1e-11, "cubic coef"); + rel_close(r[1][i], se[i], 1e-10, "cubic se"); + } + rel_close(r[2][0], 0.9999837292113533, 1e-12, "cubic r2"); + rel_close(r[4][1], 70.35586635586635, 1e-9, "cubic ssresid"); + let quint: Vec = xs.iter().zip(noise).map(|(x, e)| (0..=5).map(|j| x.powi(j)).sum::() + 50.0 * e).collect(); + let r = linest_rows(&quint, &(1..=5).map(|j| xs.iter().map(|x| x.powi(j)).collect()).collect::>>(), true); + let coef = [0.9330693815987934, 3.022273070067188, -20.97398190045249, 108.20489681886741, -250.24940010969422, 271.45454545454544]; + let se = [0.20046539191103213, 6.538253746916523, 78.39635356777087, 421.2608646012634, 976.3504976093745, 736.5866411764556]; + for i in 0..6 { + rel_close(r[0][i], coef[i], 1e-8, "quintic coef"); + rel_close(r[1][i], se[i], 1e-8, "quintic se"); + } + rel_close(r[2][1], 168.58174051770763, 1e-9, "quintic sey"); + rel_close(r[4][1], 170518.81941587824, 1e-7, "quintic ssresid"); +} + +#[test] +fn linest_drops_collinear_columns() { + // x2 = 2·x1: the later column is dropped (coefficient and standard error 0), the fit is the + // one on x1 alone, and the degrees of freedom count only the kept column. + let x1 = [1.0, 2.0, 3.0, 4.0, 5.0]; + let y = [2.0, 4.0, 5.0, 4.0, 5.0]; + let x2: Vec = x1.iter().map(|v| 2.0 * v).collect(); + let r = linest_rows(&y, &[x1.to_vec(), x2], true); + // Exact: slope 0.6, intercept 2.2, ssresid 2.4, ssreg 3.6, df 3 (Fraction arithmetic). + assert_eq!(r[0][0], Some(0.0)); + rel_close(r[0][1], 0.6, 1e-12, "slope"); + rel_close(r[0][2], 2.2, 1e-12, "intercept"); + assert_eq!(r[1][0], Some(0.0)); + rel_close(r[1][1], (0.8f64 / 10.0).sqrt(), 1e-12, "se slope"); + assert_eq!(r[3][1], Some(3.0)); + rel_close(r[4][0], 3.6, 1e-12, "ssreg"); + rel_close(r[4][1], 2.4, 1e-12, "ssresid"); + // A constant column is collinear with the intercept and dropped as well. + let r = linest_rows(&y, &[vec![7.0; 5], x1.to_vec()], true); + assert_eq!(r[0][1], Some(0.0)); + 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"); + // 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); +} + +#[test] +fn logest_is_linest_of_ln_y() { + // y = 5·2^x: m = 2, b = 5. + let x = [0.0, 1.0, 2.0, 3.0, 4.0]; + let y: Vec = x.iter().map(|v| 5.0 * 2f64.powf(*v)).collect(); + let l = rows_of(&ev("LOGEST", vec![c(&y), c(&x)])); + close_tol(l[0][0].clone(), 2.0, 1e-12); + close_tol(l[0][1].clone(), 5.0, 1e-12); + // The statistics are those of the fit on ln y, untransformed. + let ln: Vec = y.iter().map(|v| v.ln()).collect(); + let ls = rows_of(&ev("LOGEST", vec![c(&y), c(&x), b(true), b(true)])); + let ns = rows_of(&ev("LINEST", vec![c(&ln), c(&x), b(true), b(true)])); + assert_eq!(ls[1..], ns[1..]); + close_tol(ls[0][0].clone(), ns[0][0].as_f64().unwrap().exp(), 1e-14); +} + // #2177: one predictor row of 16,384 columns built three 16,384² matrices (6 GB) and an O(k³) // inversion; a wider TOROW aborted the process. It is now refused before allocating. #[test] From d0b40c3b3029adae0a600aa3f874c154d5e0eaee Mon Sep 17 00:00:00 2001 From: vangel76 Date: Wed, 4 Nov 2026 00:05:38 +0100 Subject: [PATCH 2/5] 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); From a5c74edfecd7152add1748c5285e15a92bd859cf Mon Sep 17 00:00:00 2001 From: vangel76 Date: Wed, 4 Nov 2026 00:18:00 +0100 Subject: [PATCH 3/5] linest test: use f64::from directly (clippy) Co-Authored-By: Claude Sonnet 5.5 --- crates/functions/src/stat/tests.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/functions/src/stat/tests.rs b/crates/functions/src/stat/tests.rs index ccbbf8a3a..1dd208111 100644 --- a/crates/functions/src/stat/tests.rs +++ b/crates/functions/src/stat/tests.rs @@ -443,7 +443,7 @@ fn linest_drops_collinear_columns() { // 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 xs: Vec = (1..=n).map(f64::from).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}"); From 583c8a4819bd8f1a94431d041203d3e1c53e23b9 Mon Sep 17 00:00:00 2001 From: vangel76 Date: Sat, 7 Nov 2026 17:01:20 +0100 Subject: [PATCH 4/5] 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); From 86ed3d93b49a75fbd89ac4bb38f8c5d13528b7f8 Mon Sep 17 00:00:00 2001 From: vangel76 Date: Sun, 8 Nov 2026 12:19:34 +0100 Subject: [PATCH 5/5] LINEST: treat only identical raw values as a constant predictor The 64-ulp spread threshold still classified distinct representable predictors such as [1e15, 1e15 + 1, 1e15 + 2] as constant and dropped them. A column is now constant only when every raw value is the same double; the 1e15 case is added to the offset regression test. Co-Authored-By: Claude Opus 5.5 --- crates/functions/src/stat/regression.rs | 11 ++++------- crates/functions/src/stat/tests.rs | 9 +++++++++ 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/crates/functions/src/stat/regression.rs b/crates/functions/src/stat/regression.rs index fe9515860..622022ee9 100644 --- a/crates/functions/src/stat/regression.rs +++ b/crates/functions/src/stat/regression.rs @@ -130,10 +130,6 @@ fn norm2(v: &[f64]) -> f64 { m * v.iter().map(|x| (x / m) * (x / m)).sum::().sqrt() } -/// 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. struct QrFit { @@ -276,11 +272,12 @@ fn regress(y: &[f64], x: &[Vec], k: usize, constant: bool) -> R { 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() }; // Only centring leaves rounding residue behind; without an intercept a constant column is data. + // Constant means every raw value is the same double: distinct values, however close, vary. 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 + let mut col = x.iter().map(|r| r.get(j).copied().unwrap_or(0.0)); + let first = col.next(); + constant && first.is_some_and(|f| col.all(|v| v == f)) }) .collect(); let QrFit { coef, dropped, ssresid, diag, quad } = qr_solve(cols, yc, &xm, &constant_cols); diff --git a/crates/functions/src/stat/tests.rs b/crates/functions/src/stat/tests.rs index 1ce6ae1ce..10fd00e09 100644 --- a/crates/functions/src/stat/tests.rs +++ b/crates/functions/src/stat/tests.rs @@ -582,6 +582,15 @@ fn linest_drops_collinear_columns() { 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"); + // Same at 1e15, where the spread is only a few ulps of the offset: still distinct doubles. + let x15 = [1e15, 1e15 + 1.0, 1e15 + 2.0]; + let r = linest_rows(&[0.0, 1.0, 2.0], &[x15.to_vec()], true); + rel_close(r[0][0], 1.0, 1e-9, "slope on a predictor with a 1e15 offset"); + rel_close(r[0][1], -1e15, 1e-9, "intercept on a predictor with a 1e15 offset"); + assert_eq!(r[3][1], Some(1.0), "df: the 1e15 predictor is kept"); + assert!(r[4][1].is_some_and(|v| v.abs() < 1e-9), "ssresid: the 1e15 line is exact, got {:?}", r[4][1]); + let fitted = rows_of(&ev("TREND", vec![c(&[0.0, 1.0, 2.0]), c(&x15)])); + close(fitted[2][0].clone(), 2.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);