pr_ml/svm/

binary.rs

//! Binary classification support vector machine.

use super::{FittedSVMDataPoint, Kernel, RowVector, SVMParams};

impl<const D: usize, K> SVMParams<D, K>
where
    K: Kernel<D>,
{
    /// Fits a [`BinarySVM`] to the given data.
    #[allow(clippy::similar_names, reason = "necessary for clarity")]
    #[allow(clippy::too_many_lines, reason = "SMO algorithm implementation")]
    pub fn fit_binary<I>(self, data: I) -> BinarySVM<D, K>
    where
        I: IntoIterator<Item = (RowVector<D>, bool)>,
    {
        // Collect data into vectors
        let mut data_points: Vec<_> = data
            .into_iter()
            .map(|(x, y)| FittedSVMDataPoint { x, y, alpha: 0.0 })
            .collect();

        let mut bias = 0.0;
        let kernel = self.kernel;
        let n = data_points.len();

        if n == 0 {
            return BinarySVM {
                kernel,
                bias,
                data_points,
            };
        }

        // SMO algorithm - no upfront kernel cache to save memory
        for _ in 0..self.max_iter {
            let mut num_changed_alphas = 0;

            for i in 0..n {
                // Calculate E_i = f(x_i) - y_i
                let mut f_xi = bias;
                for j in 0..n {
                    if data_points[j].alpha > 0.0 {
                        let y_j = if data_points[j].y { 1.0 } else { -1.0 };
                        f_xi += data_points[j].alpha
                            * y_j
                            * kernel.compute(&data_points[j].x, &data_points[i].x);
                    }
                }
                let y_i = if data_points[i].y { 1.0 } else { -1.0 };
                let e_i = f_xi - y_i;

                // Check KKT conditions
                let r_i = e_i * y_i;
                if (r_i < -self.tol && data_points[i].alpha < self.c)
                    || (r_i > self.tol && data_points[i].alpha > 0.0)
                {
                    // Select j != i
                    let mut j = i;
                    while j == i {
                        j = (j + 1) % n;
                    }

                    // Calculate E_j
                    let mut f_x_j = bias;
                    for k in 0..n {
                        if data_points[k].alpha > 0.0 {
                            let y_k = if data_points[k].y { 1.0 } else { -1.0 };
                            f_x_j += data_points[k].alpha
                                * y_k
                                * kernel.compute(&data_points[k].x, &data_points[j].x);
                        }
                    }
                    let y_j = if data_points[j].y { 1.0 } else { -1.0 };
                    let e_j = f_x_j - y_j;

                    // Save old alphas
                    let alpha_i_old = data_points[i].alpha;
                    let alpha_j_old = data_points[j].alpha;

                    // Compute bounds L and H
                    let (l, h) = if data_points[i].y == data_points[j].y {
                        let l = f32::max(0.0, alpha_i_old + alpha_j_old - self.c);
                        let h = f32::min(self.c, alpha_i_old + alpha_j_old);
                        (l, h)
                    } else {
                        let l = f32::max(0.0, alpha_j_old - alpha_i_old);
                        let h = f32::min(self.c, self.c + alpha_j_old - alpha_i_old);
                        (l, h)
                    };
                    if (h - l).abs() < 1e-8 {
                        continue;
                    }

                    // Compute eta
                    let k_ii = kernel.compute(&data_points[i].x, &data_points[i].x);
                    let k_jj = kernel.compute(&data_points[j].x, &data_points[j].x);
                    let k_ij = kernel.compute(&data_points[i].x, &data_points[j].x);
                    let eta = 2.0f32.mul_add(k_ij, -k_ii) - k_jj;
                    if eta >= 0.0 {
                        continue;
                    }

                    // Update alpha_j
                    let alpha_j_new = alpha_j_old - y_j * (e_i - e_j) / eta;
                    let alpha_j_new = f32::min(h, f32::max(l, alpha_j_new));
                    if (alpha_j_new - alpha_j_old).abs() < 1e-5 {
                        continue;
                    }

                    // Update alpha_i & bias
                    let alpha_i_new = (y_i * y_j).mul_add(alpha_j_old - alpha_j_new, alpha_i_old);

                    let b1 = (y_j * (alpha_j_new - alpha_j_old)).mul_add(
                        -k_ij,
                        (y_i * (alpha_i_new - alpha_i_old)).mul_add(-k_ii, bias - e_i),
                    );

                    let b2 = (y_j * (alpha_j_new - alpha_j_old)).mul_add(
                        -k_jj,
                        (y_i * (alpha_i_new - alpha_i_old)).mul_add(-k_ij, bias - e_j),
                    );

                    bias = if alpha_i_new > 0.0 && alpha_i_new < self.c {
                        b1
                    } else if alpha_j_new > 0.0 && alpha_j_new < self.c {
                        b2
                    } else {
                        f32::midpoint(b1, b2)
                    };

                    // Update alphas
                    data_points[i].alpha = alpha_i_new;
                    data_points[j].alpha = alpha_j_new;
                    num_changed_alphas += 1;
                }
            }

            if num_changed_alphas == 0 {
                break;
            }
        }

        BinarySVM {
            kernel,
            bias,
            data_points,
        }
    }
}

/// A fitted binary classification support vector machine.
///
/// # Type Parameters
///
/// - `D` - The dimension or number of features.
/// - `K` - The type of the kernel function.
#[derive(Debug, Clone, PartialEq)]
pub struct BinarySVM<const D: usize, K>
where
    K: Kernel<D>,
{
    /// The kernel function.
    kernel: K,
    /// Bias term.
    bias: f32,
    /// Data points used in training.
    data_points: Vec<FittedSVMDataPoint<D>>,
}

impl<const D: usize, K> BinarySVM<D, K>
where
    K: Kernel<D>,
{
    /// Creates a new [`SVMParams`] for fitting a [`BinarySVM`].
    #[must_use]
    pub fn params() -> SVMParams<D, K> {
        SVMParams::new()
    }

    /// Predicts the class label for a given input vector.
    pub fn predict(&self, x: &RowVector<D>) -> bool {
        self.decision_function(x) >= 0.0
    }

    /// Computes the decision value for a given input vector.
    ///
    /// This is useful for getting a confidence score. A larger positive value indicates higher confidence in the positive class, while a larger negative value indicates higher confidence in the negative class.
    pub fn decision_function(&self, x: &RowVector<D>) -> f32 {
        // Compute decision function: f(x) = sum(alpha_i * y_i * K(x_i, x)) + bias
        let mut decision_value = self.bias;

        for point in &self.data_points {
            if point.alpha > 0.0 {
                let y = if point.y { 1.0 } else { -1.0 };
                decision_value += point.alpha * y * self.kernel.compute(&point.x, x);
            }
        }

        decision_value
    }

    /// Returns an iterator over the support vectors and their corresponding alpha values.
    pub fn support_vectors(&self) -> impl Iterator<Item = &FittedSVMDataPoint<D>> {
        self.data_points.iter().filter(|p| p.alpha > 0.0)
    }

    /// Returns the number of support vectors.
    ///
    /// Support vectors are data points with non-zero alpha values.
    pub fn num_support_vectors(&self) -> usize {
        self.support_vectors().count()
    }

    /// Returns the bias term of the SVM.
    pub const fn bias(&self) -> f32 {
        self.bias
    }
}