solow_decomposition/
kernel_pca.rs

1//! Kernel PCA (Schölkopf, Smola & Müller 1998).
2
3use ndarray::{Array1, Array2, ArrayView2};
4use solow_core::{Error, Result};
5use solow_manifold::isomap::jacobi_symmetric;
6
7/// Kernel family used by [`KernelPca`].
8#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
9#[derive(Copy, Clone, Debug, PartialEq)]
10pub enum KernelKind {
11    /// Linear kernel `⟨x, y⟩`.
12    Linear,
13    /// RBF (Gaussian) kernel `exp(−γ · ‖x − y‖²)`.
14    Rbf {
15        /// Bandwidth γ.
16        gamma: f64,
17    },
18    /// Polynomial kernel `(⟨x, y⟩ · γ + coef0)^degree`.
19    Polynomial {
20        /// Polynomial degree.
21        degree: u32,
22        /// Scale factor γ on the dot product.
23        gamma: f64,
24        /// Additive offset.
25        coef0: f64,
26    },
27}
28
29impl KernelKind {
30    fn apply(&self, xi: &[f64], xj: &[f64]) -> f64 {
31        match self {
32            KernelKind::Linear => dot(xi, xj),
33            KernelKind::Rbf { gamma } => {
34                let mut s = 0.0_f64;
35                for k in 0..xi.len() {
36                    let d = xi[k] - xj[k];
37                    s += d * d;
38                }
39                (-gamma * s).exp()
40            }
41            KernelKind::Polynomial {
42                degree,
43                gamma,
44                coef0,
45            } => {
46                let base = gamma * dot(xi, xj) + coef0;
47                let mut acc = 1.0_f64;
48                for _ in 0..*degree {
49                    acc *= base;
50                }
51                acc
52            }
53        }
54    }
55}
56
57fn dot(a: &[f64], b: &[f64]) -> f64 {
58    let mut s = 0.0_f64;
59    for k in 0..a.len() {
60        s += a[k] * b[k];
61    }
62    s
63}
64
65/// Fitted kernel-PCA embedding.
66#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
67#[derive(Clone, Debug, PartialEq)]
68pub struct KernelPca {
69    /// Low-dimensional embedding, `n × n_components`.
70    pub embedding: Array2<f64>,
71    /// Retained eigenvalues.
72    pub eigenvalues: Array1<f64>,
73    /// Kernel used at fit time.
74    pub kernel: KernelKind,
75}
76
77impl KernelPca {
78    /// Fit onto `x` with the given kernel and target dimension.
79    pub fn fit(x: ArrayView2<'_, f64>, kernel: KernelKind, n_components: usize) -> Result<Self> {
80        if x.nrows() == 0 || x.ncols() == 0 {
81            return Err(Error::Value("KernelPca::fit: x must be non-empty".into()));
82        }
83        if n_components < 1 || n_components > x.nrows() {
84            return Err(Error::Value(format!(
85                "KernelPca::fit: n_components in [1, n] (got {n_components})"
86            )));
87        }
88        let n = x.nrows();
89        // Build kernel matrix.
90        let mut k = Array2::<f64>::zeros((n, n));
91        for i in 0..n {
92            let xi: Vec<f64> = x.row(i).to_vec();
93            for j in 0..=i {
94                let xj: Vec<f64> = x.row(j).to_vec();
95                let v = kernel.apply(&xi, &xj);
96                k[[i, j]] = v;
97                k[[j, i]] = v;
98            }
99        }
100        // Centre the kernel matrix: K_c = K − 1_N·K − K·1_N + 1_N·K·1_N.
101        let mut row_means = vec![0.0_f64; n];
102        let mut col_means = vec![0.0_f64; n];
103        let mut grand = 0.0_f64;
104        for i in 0..n {
105            for j in 0..n {
106                row_means[i] += k[[i, j]];
107                col_means[j] += k[[i, j]];
108                grand += k[[i, j]];
109            }
110        }
111        for v in row_means.iter_mut() {
112            *v /= n as f64;
113        }
114        for v in col_means.iter_mut() {
115            *v /= n as f64;
116        }
117        grand /= (n * n) as f64;
118        let mut kc = k.clone();
119        for i in 0..n {
120            for j in 0..n {
121                kc[[i, j]] += grand - row_means[i] - col_means[j];
122            }
123        }
124        // Eigendecompose (kc is symmetric).
125        let (eigvals, eigvecs) = jacobi_symmetric(&kc, 300, 1e-12);
126        let mut order: Vec<usize> = (0..n).collect();
127        order.sort_by(|&a, &b| eigvals[b].partial_cmp(&eigvals[a]).unwrap());
128        let mut embedding = Array2::<f64>::zeros((n, n_components));
129        let mut ev_out = Array1::<f64>::zeros(n_components);
130        for c in 0..n_components {
131            let idx = order[c];
132            let ev = eigvals[idx].max(0.0);
133            ev_out[c] = ev;
134            let scale = ev.sqrt();
135            for i in 0..n {
136                embedding[[i, c]] = eigvecs[[i, idx]] * scale;
137            }
138        }
139        Ok(Self {
140            embedding,
141            eigenvalues: ev_out,
142            kernel,
143        })
144    }
145}
146
147#[cfg(test)]
148mod tests {
149    use super::*;
150    use ndarray::array;
151
152    #[test]
153    fn kernel_pca_linear_recovers_pca_directions_up_to_sign() {
154        // With a linear kernel KernelPCA is standard PCA.
155        let x = array![
156            [1.0, 0.0],
157            [0.0, 1.0],
158            [-1.0, 0.0],
159            [0.0, -1.0],
160            [0.5, 0.5],
161            [-0.5, -0.5]
162        ];
163        let kpca = KernelPca::fit(x.view(), KernelKind::Linear, 2).unwrap();
164        assert_eq!(kpca.embedding.dim(), (6, 2));
165        // Eigenvalues descending and non-negative.
166        assert!(kpca.eigenvalues[0] + 1e-9 >= kpca.eigenvalues[1]);
167        assert!(kpca.eigenvalues[0] >= 0.0);
168    }
169}
solow_decomposition/kernel_pca.rs

solow_decomposition/
kernel_pca.rs
