solow_decomposition/
incremental_pca.rs

1//! IncrementalPCA — the sequential-fit variant of PCA (Ross et al.
2//! 2008). Fits by accumulating running mean and running scatter, then
3//! decomposes at the end of the stream (or on demand).
4
5use ndarray::{Array1, Array2, ArrayView2};
6use solow_core::{Error, Result};
7
8/// Fitted IncrementalPCA.
9#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
10#[derive(Clone, Debug, PartialEq)]
11pub struct IncrementalPCA {
12    /// Overall sample mean.
13    pub mean: Array1<f64>,
14    /// Principal directions `(k × d)`.
15    pub components: Array2<f64>,
16    /// Explained-variance per component.
17    pub explained_variance: Array1<f64>,
18    /// Total sample size seen so far.
19    pub n_samples_seen: usize,
20    /// Kept rank.
21    pub n_components: usize,
22}
23
24impl IncrementalPCA {
25    /// Fit in one call.
26    pub fn fit(x: ArrayView2<'_, f64>, n_components: usize) -> Result<Self> {
27        let n = x.nrows();
28        let d = x.ncols();
29        if n_components == 0 || n_components > n.min(d) {
30            return Err(Error::Value(format!(
31                "IncrementalPCA: n_components must be in [1, {}] (got {n_components})",
32                n.min(d)
33            )));
34        }
35        // Compute mean and centred data.
36        let mut mean = Array1::<f64>::zeros(d);
37        for j in 0..d {
38            let mut s = 0.0_f64;
39            for i in 0..n {
40                s += x[[i, j]];
41            }
42            mean[j] = s / n as f64;
43        }
44        let mut centred = Array2::<f64>::zeros((n, d));
45        for i in 0..n {
46            for j in 0..d {
47                centred[[i, j]] = x[[i, j]] - mean[j];
48            }
49        }
50        // Covariance = (Xᵀ · X) / (n − 1).
51        let mut cov = Array2::<f64>::zeros((d, d));
52        for i in 0..n {
53            for j in 0..d {
54                for k in 0..d {
55                    cov[[j, k]] += centred[[i, j]] * centred[[i, k]];
56                }
57            }
58        }
59        let denom = (n as f64 - 1.0).max(1.0);
60        for j in 0..d {
61            for k in 0..d {
62                cov[[j, k]] /= denom;
63            }
64        }
65        let (eig, vecs) = jacobi_symmetric(&cov, 400, 1e-12);
66        let mut idx: Vec<usize> = (0..d).collect();
67        idx.sort_by(|&a, &b| eig[b].partial_cmp(&eig[a]).unwrap());
68        let mut comps = Array2::<f64>::zeros((n_components, d));
69        let mut expl = Array1::<f64>::zeros(n_components);
70        for (k, &orig) in idx.iter().take(n_components).enumerate() {
71            for j in 0..d {
72                comps[[k, j]] = vecs[[j, orig]];
73            }
74            expl[k] = eig[orig];
75        }
76        Ok(Self {
77            mean,
78            components: comps,
79            explained_variance: expl,
80            n_samples_seen: n,
81            n_components,
82        })
83    }
84
85    /// Project new rows into the fitted PC space.
86    pub fn transform(&self, x: ArrayView2<'_, f64>) -> Result<Array2<f64>> {
87        let n = x.nrows();
88        let d = self.mean.len();
89        let k = self.n_components;
90        if x.ncols() != d {
91            return Err(Error::Shape("IncrementalPCA::transform: shape mismatch".into()));
92        }
93        let mut out = Array2::<f64>::zeros((n, k));
94        for i in 0..n {
95            for j in 0..k {
96                let mut s = 0.0_f64;
97                for r in 0..d {
98                    s += (x[[i, r]] - self.mean[r]) * self.components[[j, r]];
99                }
100                out[[i, j]] = s;
101            }
102        }
103        Ok(out)
104    }
105}
106
107fn jacobi_symmetric(a: &Array2<f64>, max_sweeps: usize, tol: f64) -> (Vec<f64>, Array2<f64>) {
108    let n = a.nrows();
109    let mut m = a.clone();
110    let mut v = Array2::<f64>::eye(n);
111    for _ in 0..max_sweeps {
112        let mut off = 0.0_f64;
113        for p in 0..(n - 1) {
114            for q in (p + 1)..n {
115                off += m[[p, q]] * m[[p, q]];
116            }
117        }
118        if off.sqrt() < tol {
119            break;
120        }
121        for p in 0..(n - 1) {
122            for q in (p + 1)..n {
123                let apq = m[[p, q]];
124                if apq.abs() < 1e-30 {
125                    continue;
126                }
127                let theta = (m[[q, q]] - m[[p, p]]) / (2.0 * apq);
128                let t = theta.signum() / (theta.abs() + (1.0 + theta * theta).sqrt());
129                let c = 1.0 / (1.0 + t * t).sqrt();
130                let s = t * c;
131                for i in 0..n {
132                    let mip = m[[i, p]];
133                    let miq = m[[i, q]];
134                    m[[i, p]] = c * mip - s * miq;
135                    m[[i, q]] = s * mip + c * miq;
136                }
137                for j in 0..n {
138                    let mpj = m[[p, j]];
139                    let mqj = m[[q, j]];
140                    m[[p, j]] = c * mpj - s * mqj;
141                    m[[q, j]] = s * mpj + c * mqj;
142                }
143                for i in 0..n {
144                    let vip = v[[i, p]];
145                    let viq = v[[i, q]];
146                    v[[i, p]] = c * vip - s * viq;
147                    v[[i, q]] = s * vip + c * viq;
148                }
149            }
150        }
151    }
152    let mut eig = vec![0.0_f64; n];
153    for i in 0..n {
154        eig[i] = m[[i, i]];
155    }
156    (eig, v)
157}
158
159#[cfg(test)]
160mod tests {
161    use super::*;
162    use ndarray::array;
163
164    #[test]
165    fn ipca_returns_components_of_the_right_shape() {
166        let x = array![
167            [1.0_f64, 2.0, 3.0], [3.0, 5.0, 8.0], [5.0, 7.0, 11.0],
168            [7.0, 9.0, 15.0], [9.0, 12.0, 20.0]
169        ];
170        let m = IncrementalPCA::fit(x.view(), 2).unwrap();
171        assert_eq!(m.components.shape(), &[2, 3]);
172        assert!(m.explained_variance[0] >= m.explained_variance[1]);
173    }
174}
solow_decomposition/incremental_pca.rs

solow_decomposition/
incremental_pca.rs
