qualia_core_db/solvers/learning/classification/
svm.rs

1//! Support Vector Machine (ISL ch 9, PRML ch 7) — soft-margin binary classifier
2//! trained by simplified Sequential Minimal Optimization (SMO) on the dual, with a
3//! linear or RBF (Gaussian) kernel. Kernel SVM separates classes a linear boundary
4//! cannot. Labels are boolean (true = +1, false = −1). Kernel-class `DenseLinear`
5//! (the kernel matrix) + `Divergent` (the SMO working-set loop) → CPU here.
6
7use crate::solvers::learning::LearningError;
8
9/// The kernel `K(a, b)`.
10#[derive(Debug, Clone, Copy, PartialEq)]
11pub enum Kernel {
12    /// `⟨a, b⟩`.
13    Linear,
14    /// `exp(−γ·‖a − b‖²)`.
15    Rbf { gamma: f64 },
16}
17
18impl Kernel {
19    fn eval(self, a: &[f64], b: &[f64]) -> f64 {
20        match self {
21            Kernel::Linear => a.iter().zip(b).map(|(x, y)| x * y).sum(),
22            Kernel::Rbf { gamma } => {
23                let d2: f64 = a.iter().zip(b).map(|(x, y)| (x - y) * (x - y)).sum();
24                (-gamma * d2).exp()
25            }
26        }
27    }
28}
29
30/// A fitted SVM: the support vectors (the training points with non-zero `α`) plus
31/// the bias and kernel.
32#[derive(Debug, Clone)]
33pub struct Svm {
34    sv_x: Vec<f64>,       // n_sv × p
35    sv_alpha_y: Vec<f64>, // αᵢ·yᵢ per support vector
36    b: f64,
37    kernel: Kernel,
38    p: usize,
39}
40
41struct Lcg(u64);
42impl Lcg {
43    fn below(&mut self, bound: usize) -> usize {
44        self.0 = self
45            .0
46            .wrapping_mul(6364136223846793005)
47            .wrapping_add(1442695040888963407);
48        ((self.0 >> 33) as usize) % bound.max(1)
49    }
50}
51
52/// Fill the `n×n` kernel matrix on the CPU — the always-correct floor: one
53/// `kernel.eval` per upper-triangle pair, mirrored into a symmetric matrix. This is
54/// byte-identical to the original inline loop and is used off-accelerator, below
55/// threshold, or if a forge offload errors.
56fn fill_kernel_cpu(k: &mut [f64], x: &[f64], n: usize, p: usize, kernel: Kernel) {
57    for i in 0..n {
58        for j in i..n {
59            let v = kernel.eval(&x[i * p..(i + 1) * p], &x[j * p..(j + 1) * p]);
60            k[i * n + j] = v;
61            k[j * n + i] = v;
62        }
63    }
64}
65
66/// Fit a soft-margin SVM by simplified SMO. `c` is the regularization (box) bound,
67/// `max_passes` the number of consecutive no-change sweeps to declare convergence.
68/// Fails closed on shape mismatch / a single-class target.
69pub fn fit(
70    x: &[f64],
71    y: &[bool],
72    n: usize,
73    p: usize,
74    c: f64,
75    kernel: Kernel,
76    max_passes: usize,
77    tol: f64,
78) -> Result<Svm, LearningError> {
79    if n == 0 || p == 0 || x.len() != n * p || y.len() != n {
80        return Err(LearningError::InvalidDimension);
81    }
82    if !(c > 0.0) {
83        return Err(LearningError::InsufficientData);
84    }
85    let yi: Vec<f64> = y.iter().map(|&b| if b { 1.0 } else { -1.0 }).collect();
86    let n_pos = y.iter().filter(|&&b| b).count();
87    if n_pos == 0 || n_pos == n {
88        return Err(LearningError::InsufficientData); // need both classes
89    }
90
91    // Precompute the kernel matrix (the `DenseLinear` kernel class) — best path on this
92    // machine, with the CPU floor [`fill_kernel_cpu`]. Above `GEMM_GPU_THRESHOLD` and with
93    // an accelerator present, the whole `n×n` matrix is built from one dense pass:
94    //   • `Linear`: the Gram matrix `X·Xᵀ` straight through the engine GEMM;
95    //   • `Rbf`:    `exp(−γ·‖xᵢ−xⱼ‖²)` over `dispatch::pairwise_sq_dist_f64`'s best-path
96    //     squared-distance matrix, with the elementwise `exp` on the CPU.
97    // Off accelerator or sub-threshold, the exact symmetric CPU loop runs — byte-identical
98    // to before. A forge error never propagates: it falls back to the CPU floor.
99    let mut k = vec![0.0; n * n];
100    // GPU best-path kernel-matrix build via the forge — only when it's compiled in
101    // (native + wgsl-forge). On wasm32 the exact symmetric CPU kernel loop runs.
102    #[cfg(all(not(target_arch = "wasm32"), feature = "wgsl-forge"))]
103    {
104        let work = n.saturating_mul(n).saturating_mul(p);
105        let caps = crate::wgsl_forge::dispatch::caps();
106        let accelerated =
107            (caps.cuda || caps.wgpu) && work >= crate::wgsl_forge::dispatch::GEMM_GPU_THRESHOLD;
108        if accelerated {
109            match kernel {
110                Kernel::Linear => {
111                    use crate::solvers::linear_algebra::gemm::{gemm, Transpose};
112                    // Gram = X·Xᵀ ([n×n]): op(A)=X (n×p), op(B)=Xᵀ (p×n).
113                    if gemm(
114                        Transpose::No,
115                        Transpose::Yes,
116                        n,
117                        n,
118                        p,
119                        1.0,
120                        x,
121                        x,
122                        0.0,
123                        &mut k,
124                    )
125                    .is_err()
126                    {
127                        fill_kernel_cpu(&mut k, x, n, p, kernel);
128                    }
129                }
130                Kernel::Rbf { gamma } => {
131                    let d = crate::wgsl_forge::dispatch::pairwise_sq_dist_f64(x, x, n, n, p);
132                    for (kij, &dij) in k.iter_mut().zip(d.iter()) {
133                        *kij = (-gamma * dij).exp();
134                    }
135                }
136            }
137        } else {
138            fill_kernel_cpu(&mut k, x, n, p, kernel);
139        }
140    }
141    #[cfg(not(all(not(target_arch = "wasm32"), feature = "wgsl-forge")))]
142    fill_kernel_cpu(&mut k, x, n, p, kernel);
143
144    let mut alpha = vec![0.0; n];
145    let mut b = 0.0;
146    let mut rng = Lcg(0x9E3779B97F4A7C15);
147
148    // f(xᵢ) = Σ_m α_m y_m K(m,i) + b.
149    let f = |alpha: &[f64], b: f64, i: usize, k: &[f64]| -> f64 {
150        let mut s = b;
151        for m in 0..n {
152            if alpha[m] != 0.0 {
153                s += alpha[m] * yi[m] * k[m * n + i];
154            }
155        }
156        s
157    };
158
159    let mut passes = 0;
160    let max_iter = max_passes.max(1);
161    let hard_cap = 10_000; // total outer sweeps guard
162    let mut sweeps = 0;
163    while passes < max_iter && sweeps < hard_cap {
164        sweeps += 1;
165        let mut num_changed = 0;
166        for i in 0..n {
167            let ei = f(&alpha, b, i, &k) - yi[i];
168            if (yi[i] * ei < -tol && alpha[i] < c) || (yi[i] * ei > tol && alpha[i] > 0.0) {
169                // Pick j ≠ i.
170                let mut j = rng.below(n);
171                if j == i {
172                    j = (j + 1) % n;
173                }
174                let ej = f(&alpha, b, j, &k) - yi[j];
175                let (ai_old, aj_old) = (alpha[i], alpha[j]);
176                // Bounds on α_j.
177                let (lo, hi) = if yi[i] != yi[j] {
178                    ((aj_old - ai_old).max(0.0), c + (aj_old - ai_old).min(0.0))
179                } else {
180                    ((ai_old + aj_old - c).max(0.0), (ai_old + aj_old).min(c))
181                };
182                if (hi - lo).abs() < 1e-12 {
183                    continue;
184                }
185                let eta = 2.0 * k[i * n + j] - k[i * n + i] - k[j * n + j];
186                if eta >= 0.0 {
187                    continue;
188                }
189                let mut aj = aj_old - yi[j] * (ei - ej) / eta;
190                aj = aj.clamp(lo, hi);
191                if (aj - aj_old).abs() < 1e-9 {
192                    continue;
193                }
194                let ai = ai_old + yi[i] * yi[j] * (aj_old - aj);
195                // Bias update.
196                let b1 = b
197                    - ei
198                    - yi[i] * (ai - ai_old) * k[i * n + i]
199                    - yi[j] * (aj - aj_old) * k[i * n + j];
200                let b2 = b
201                    - ej
202                    - yi[i] * (ai - ai_old) * k[i * n + j]
203                    - yi[j] * (aj - aj_old) * k[j * n + j];
204                alpha[i] = ai;
205                alpha[j] = aj;
206                b = if ai > 0.0 && ai < c {
207                    b1
208                } else if aj > 0.0 && aj < c {
209                    b2
210                } else {
211                    0.5 * (b1 + b2)
212                };
213                num_changed += 1;
214            }
215        }
216        if num_changed == 0 {
217            passes += 1;
218        } else {
219            passes = 0;
220        }
221    }
222
223    // Keep only the support vectors (α > 0).
224    let mut sv_x = Vec::new();
225    let mut sv_alpha_y = Vec::new();
226    for i in 0..n {
227        if alpha[i] > 1e-8 {
228            sv_x.extend_from_slice(&x[i * p..(i + 1) * p]);
229            sv_alpha_y.push(alpha[i] * yi[i]);
230        }
231    }
232    Ok(Svm {
233        sv_x,
234        sv_alpha_y,
235        b,
236        kernel,
237        p,
238    })
239}
240
241impl Svm {
242    /// The signed decision value `Σ αᵢyᵢ K(svᵢ, q) + b`.
243    pub fn decision_row(&self, q: &[f64]) -> f64 {
244        let mut s = self.b;
245        for (k, &ay) in self.sv_alpha_y.iter().enumerate() {
246            s += ay
247                * self
248                    .kernel
249                    .eval(&self.sv_x[k * self.p..(k + 1) * self.p], q);
250        }
251        s
252    }
253
254    /// Predicted class (`true` = +1) = `decision ≥ 0`.
255    pub fn predict_row(&self, q: &[f64]) -> bool {
256        self.decision_row(q) >= 0.0
257    }
258
259    pub fn n_support_vectors(&self) -> usize {
260        self.sv_alpha_y.len()
261    }
262}
263
264#[cfg(test)]
265mod tests {
266    use super::*;
267
268    #[test]
269    fn linear_svm_separates_linearly_separable() {
270        // Two clearly separated 2-D classes.
271        let x = [
272            0.0, 0.0, 1.0, 0.5, 0.5, 1.0, -0.5, 0.2, 6.0, 6.0, 5.5, 6.5, 6.5, 5.5, 5.0, 6.0,
273        ];
274        let y = [false, false, false, false, true, true, true, true];
275        let svm = fit(&x, &y, 8, 2, 1.0, Kernel::Linear, 20, 1e-3).unwrap();
276        assert!(!svm.predict_row(&[0.2, 0.2]));
277        assert!(svm.predict_row(&[6.0, 6.0]));
278        assert!(svm.n_support_vectors() >= 1);
279        // Perfect training accuracy on a separable set.
280        for i in 0..8 {
281            assert_eq!(svm.predict_row(&x[i * 2..i * 2 + 2]), y[i]);
282        }
283    }
284
285    #[test]
286    fn rbf_svm_handles_nonlinear_boundary() {
287        // Concentric-ish: inner points class 0, outer ring class 1 — not linearly
288        // separable, but an RBF kernel handles it.
289        let mut x = Vec::new();
290        let mut y = Vec::new();
291        // inner cluster (class false) near origin
292        for &(a, b) in &[(0.0, 0.0), (0.3, 0.0), (0.0, 0.3), (-0.3, 0.0), (0.0, -0.3)] {
293            x.push(a);
294            x.push(b);
295            y.push(false);
296        }
297        // outer ring (class true)
298        for &(a, b) in &[
299            (3.0, 0.0),
300            (-3.0, 0.0),
301            (0.0, 3.0),
302            (0.0, -3.0),
303            (2.1, 2.1),
304            (-2.1, -2.1),
305        ] {
306            x.push(a);
307            x.push(b);
308            y.push(true);
309        }
310        let n = 11;
311        let svm = fit(&x, &y, n, 2, 1.0, Kernel::Rbf { gamma: 0.5 }, 50, 1e-3).unwrap();
312        assert!(
313            !svm.predict_row(&[0.1, 0.1]),
314            "inner point should be class 0"
315        );
316        assert!(
317            svm.predict_row(&[3.0, 0.0]),
318            "outer point should be class 1"
319        );
320        assert!(svm.predict_row(&[0.0, -3.0]));
321    }
322
323    #[test]
324    fn guards() {
325        assert_eq!(
326            fit(
327                &[1.0, 2.0],
328                &[true, true][..1],
329                1,
330                2,
331                1.0,
332                Kernel::Linear,
333                5,
334                1e-3
335            )
336            .unwrap_err(),
337            LearningError::InsufficientData
338        );
339        let x = [0.0, 0.0, 1.0, 1.0];
340        assert_eq!(
341            fit(&x, &[true, true], 2, 2, 1.0, Kernel::Linear, 5, 1e-3).unwrap_err(),
342            LearningError::InsufficientData
343        );
344    }
345}
qualia_core_db/solvers/learning/classification/svm.rs

qualia_core_db/solvers/learning/classification/
svm.rs
