qualia_core_db/solvers/learning/splines/
smoothing.rs

1//! Penalized smoothing spline (ISL ch 7.5) — least squares with a roughness
2//! penalty that shrinks the wiggly (knot) part of the fit.
3//!
4//! Minimise `‖y − Bβ‖² + λ·βᵀPβ`, where `B` is the truncated-power spline basis
5//! (shared with [`super::RegressionSpline`]) and `P` penalizes only the
6//! truncated-power (knot) coefficients — the part that controls smoothness. `λ = 0`
7//! reproduces the (interpolating) regression spline; large `λ` shrinks the knot
8//! terms toward a global polynomial (a smooth fit). The penalized normal equations
9//! `(BᵀB + λP)β = Bᵀy` are solved with `linear_algebra::cholesky` (no new solver).
10//! Kernel-class `DenseLinear`.
11
12use crate::solvers::learning::splines::{basis_len, basis_row};
13use crate::solvers::learning::LearningError;
14use crate::solvers::linear_algebra::cholesky::{cholesky_factor, cholesky_solve};
15use crate::solvers::linear_algebra::gemm::{gemm, matvec, Transpose};
16
17/// A fitted smoothing spline.
18#[derive(Debug, Clone)]
19pub struct SmoothingSpline {
20    pub degree: usize,
21    pub knots: Vec<f64>,
22    pub coefficients: Vec<f64>,
23    pub lambda: f64,
24}
25
26impl SmoothingSpline {
27    /// Fit with smoothing parameter `lambda ≥ 0`. Fails closed on shape mismatch /
28    /// `degree == 0`.
29    pub fn fit(
30        x: &[f64],
31        y: &[f64],
32        n: usize,
33        degree: usize,
34        knots: &[f64],
35        lambda: f64,
36    ) -> Result<Self, LearningError> {
37        if n == 0 || x.len() != n || y.len() != n || degree == 0 || lambda < 0.0 {
38            return Err(LearningError::InvalidDimension);
39        }
40        let m = basis_len(degree, knots.len());
41        // Design matrix B (n × m).
42        let mut b = vec![0.0; n * m];
43        for i in 0..n {
44            basis_row(x[i], degree, knots, &mut b[i * m..(i + 1) * m]);
45        }
46        // A = BᵀB + λP, where P = diag(0,…,0,1,…,1) penalizes only the knot terms.
47        let mut a = vec![0.0; m * m];
48        gemm(
49            Transpose::Yes,
50            Transpose::No,
51            m,
52            m,
53            n,
54            1.0,
55            &b,
56            &b,
57            0.0,
58            &mut a,
59        )?;
60        for j in (degree + 1)..m {
61            a[j * m + j] += lambda;
62        }
63        // rhs = Bᵀy.
64        let mut rhs = vec![0.0; m];
65        matvec(Transpose::Yes, m, n, &b, y, &mut rhs)?;
66        // Solve.
67        let mut l = vec![0.0; m * m];
68        cholesky_factor(m, &a, &mut l).map_err(|_| LearningError::Singular)?;
69        let mut coefficients = vec![0.0; m];
70        cholesky_solve(m, &l, &rhs, &mut coefficients)?;
71        Ok(Self {
72            degree,
73            knots: knots.to_vec(),
74            coefficients,
75            lambda,
76        })
77    }
78
79    pub fn predict_one(&self, x: f64) -> f64 {
80        let m = basis_len(self.degree, self.knots.len());
81        let mut row = vec![0.0; m];
82        basis_row(x, self.degree, &self.knots, &mut row);
83        row.iter().zip(&self.coefficients).map(|(b, c)| b * c).sum()
84    }
85
86    pub fn predict(&self, x: &[f64]) -> Vec<f64> {
87        x.iter().map(|&xi| self.predict_one(xi)).collect()
88    }
89}
90
91#[cfg(test)]
92mod tests {
93    use super::*;
94    use crate::solvers::learning::splines::RegressionSpline;
95
96    /// Roughness proxy: sum of squared second differences of the fitted values.
97    fn roughness(vals: &[f64]) -> f64 {
98        let mut r = 0.0;
99        for i in 1..vals.len() - 1 {
100            let d2 = vals[i + 1] - 2.0 * vals[i] + vals[i - 1];
101            r += d2 * d2;
102        }
103        r
104    }
105
106    #[test]
107    fn lambda_zero_matches_regression_spline() {
108        let x: Vec<f64> = (0..20).map(|i| i as f64 * 0.3).collect();
109        let y: Vec<f64> = x.iter().map(|&xi| xi.sin()).collect();
110        let knots = vec![1.5, 3.0, 4.5];
111        let ss = SmoothingSpline::fit(&x, &y, 20, 3, &knots, 0.0).unwrap();
112        let rs = RegressionSpline::fit(&x, &y, 20, 3, &knots).unwrap();
113        for i in 0..20 {
114            assert!((ss.predict_one(x[i]) - rs.predict_one(x[i])).abs() < 1e-6);
115        }
116    }
117
118    #[test]
119    fn larger_lambda_is_smoother() {
120        // Noisy data; a larger penalty yields a smoother (less wiggly) fit.
121        let n = 40;
122        let x: Vec<f64> = (0..n).map(|i| i as f64 * 0.2).collect();
123        let y: Vec<f64> = x
124            .iter()
125            .enumerate()
126            .map(|(i, &xi)| xi.sin() + ((i % 2) as f64 - 0.5) * 0.6)
127            .collect();
128        let knots: Vec<f64> = (1..8).map(|k| k as f64).collect();
129        let light = SmoothingSpline::fit(&x, &y, n, 3, &knots, 0.01).unwrap();
130        let heavy = SmoothingSpline::fit(&x, &y, n, 3, &knots, 100.0).unwrap();
131        let r_light = roughness(&light.predict(&x));
132        let r_heavy = roughness(&heavy.predict(&x));
133        assert!(
134            r_heavy < r_light,
135            "heavier penalty must be smoother: {r_heavy} !< {r_light}"
136        );
137    }
138
139    #[test]
140    fn guards() {
141        assert_eq!(
142            SmoothingSpline::fit(&[1.0, 2.0], &[1.0], 2, 3, &[], 1.0).unwrap_err(),
143            LearningError::InvalidDimension
144        );
145    }
146}
qualia_core_db/solvers/learning/splines/smoothing.rs

qualia_core_db/solvers/learning/splines/
smoothing.rs
