//! Linear Support Vector Classifier.
//!
//! This module provides [`LinearSVC`], a liblinear-faithful linear support
//! vector classifier that operates directly in the primal space without the
//! overhead of a kernel function. The fit minimizes the L2-regularized
//! hinge (or squared-hinge) classification objective
//!
//! ```text
//! min_w 0.5 * ||w||^2 + C * sum_i L(y_i, w . x_i)
//! ```
//!
//! with `y_i ∈ {-1, +1}` per one-vs-rest sub-problem (NO `1/n` averaging — the
//! summed loss is scaled by `C`, matching `sklearn/svm/_base.py`
//! `_fit_liblinear`). The solver is liblinear's dual coordinate descent for
//! classification (`solve_l2r_l1l2_svc` in
//! `sklearn/svm/src/liblinear/linear.cpp:819`), which converges to the unique
//! minimizer of the strongly convex objective.
//!
//! Unlike [`SVC`](crate::svm::SVC) with a [`LinearKernel`](crate::svm::LinearKernel),
//! `LinearSVC` avoids computing and caching the full kernel matrix, making it
//! significantly faster for high-dimensional data.
//!
//! # Examples
//!
//! ```
//! use ferrolearn_linear::linear_svc::LinearSVC;
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array2};
//!
//! let x = Array2::from_shape_vec((6, 2), vec![
//! 1.0, 1.0, 1.0, 2.0, 2.0, 1.0,
//! 5.0, 5.0, 5.0, 6.0, 6.0, 5.0,
//! ]).unwrap();
//! let y = array![0usize, 0, 0, 1, 1, 1];
//!
//! let model = LinearSVC::<f64>::new();
//! let fitted = model.fit(&x, &y).unwrap();
//! let preds = fitted.predict(&x).unwrap();
//! assert_eq!(preds.len(), 6);
//! ```
//!
//! ## REQ status
//!
//! Binary (R-DEFER-2): SHIPPED = impl + non-test production consumer + tests +
//! green oracle verification; NOT-STARTED = open blocker `#`. `LinearSVC`/
//! `FittedLinearSVC`/`LinearSVCLoss` are boundary estimator types re-exported at
//! the crate root (`pub use linear_svc::{…}` in `lib.rs`) and registered as the
//! PyO3 `RsLinearSVC` estimator (`ferrolearn-python/src/extras.rs`); under
//! S5/R-DEFER-1 those ARE the non-test production-consumer surface. See
//! `.design/linear/linear_svc.md`.
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-1 (fit parity — coef_/intercept_ vs liblinear oracle) | SHIPPED | `fn solve_binary_dual` minimizes liblinear's `0.5·‖w‖² + C·Σ L` via the dual CD (`solve_l2r_l1l2_svc`, `linear.cpp:819`); `fn fit` maps `classes_[1]→+1` and extracts `coef_ = w[:n_features]`, `intercept_ = intercept_scaling·w_last` (`_base.py:1240-1245`). Pinned by `tests/divergence_linear_svc_fit.rs::linear_svc_coef_parity` (live oracle `coef_ [[0.12835213611984458, 0.12835213611984475]]`, `intercept_ [-1.1943776585907158]`, C=1.0, squared_hinge, fit_intercept=True). Consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). |
//! | REQ-2 (decision_function shape `(n,)` + values) | SHIPPED | `fn decision_function` returns [`DecisionScores::Binary`] = 1-D `X·w + b` for the binary case (sklearn ravels the single-column score to `(n,)`, `linear_model/_base.py:365`) and [`DecisionScores::Multiclass`] `(n, n_classes)` otherwise. Pinned by `tests/divergence_linear_svc_fit.rs::linear_svc_decision_function` (live oracle 1-D `(8,)` values). Consumer: `fn predict` reads the binary scores' sign. |
//! | REQ-3 (predict + classes_) | SHIPPED | `fn predict` uses the sign of the binary decision (`>= 0 → classes_[1]`) / argmax of the OvR scores; `HasClasses::classes` = sorted unique `y` (`classes_ = np.unique(y)`, `_classes.py:311`). The labels are downstream of the liblinear-parity fit and pinned against the live oracle by `linear_svc_predict_parity in tests/divergence_linear_svc_fit.rs` (#620; 8×2 set: `predict [0,0,0,0,1,1,1,1]`, `classes_ [0,1]`). |
//! | REQ-4 (loss {hinge, squared_hinge}) | SHIPPED | The dual CD solves BOTH the true `hinge` (`U=C`, `diag=0`) and `squared_hinge` (`U=∞`, `diag=0.5/C`) optima (`solve_l2r_l1l2_svc`, `linear.cpp:849-858`). The `hinge` optimum is pinned against the live oracle by `linear_svc_hinge_coef_parity in tests/divergence_linear_svc_fit.rs` (#621; 8×2 set, `loss='hinge'`, `C=1.0`: `coef_ [[0.15384615383852776, 0.15384615383915584]]`, `intercept_ [-1.4615384615168394]`). |
//! | REQ-5 (penalty {l1, l2}) | SHIPPED | `LinearSVC<F>` exposes `pub penalty: LinearSVCPenalty` (default `L2`) + `#[must_use] with_penalty`. `penalty=l1` routes to `fn solve_binary_l1r_l2` — liblinear's feature-major coordinate descent (`solve_l1r_l2_svc`, `linear.cpp:1467`, solver type 5, `_base.py:1014`) minimizing `‖w‖₁ + C·Σ max(0,1−yf)²` (sparse `coef_`); `penalty=l2` keeps `fn solve_binary_dual`. `liblinear_solver_type` rejects `('l1','hinge')` (`_base.py:1013`). Pinned by `test_l1_penalty_smoke` (live oracle 8×2 `l1,squared_hinge,dual=False,C=1`: `coef_ [[0.1283185834966579, 0.12831858464059265]]`, `intercept_ [-1.2079646017762715]`; ferrolearn lands within ~1.2e-9) + `test_unsupported_combinations_rejected`. Consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). |
//! | REQ-6 (multi_class {ovr, crammer_singer}) | SHIPPED | `LinearSVC<F>` exposes `pub multi_class: MultiClass` (`Ovr`/`CrammerSinger`, default `Ovr`, `_classes.py:239`) + `#[must_use] with_multi_class`. `multi_class=CrammerSinger` selects liblinear solver type 4 REGARDLESS of penalty/loss/dual (`_get_liblinear_solver_type`, `_base.py:1017,1020-1021`), so `fn fit` short-circuits to `fn fit_crammer_singer` (ignoring penalty/loss/dual) which runs ONE joint solve `fn solve_crammer_singer` — a faithful transcription of `Solver_MCSVM_CS` (`linear.cpp:510`, the class at `:493-787`): flattened `w[feature*nr_class + m]`, `alpha[i*nr_class + m]` with `Σ_m alpha=0`, `C[i] = weighted_C[y_i]`, per-sample shrinking (`active_size_i`/`alpha_index`/`be_shrunk`), the simplex inner solve `fn cs_solve_sub_problem` (sort-descending breakpoint, `linear.cpp:541-564`), and the two-level `eps_shrink = max(10·tol, 1)` stopping (`linear.cpp:738-753`). Extraction: `coef_[m][feature] = w[feature*nr_class + m]`, `intercept_[m] = intercept_scaling·w[n_features*nr_class + m]` (`_base.py:1240-1245`). BINARY: collapse to a single weight vector `coef_ = row_1 − row_0`, `intercept_ = int_1 − int_0` (`_classes.py:340-344`), `is_binary=true`. ferrolearn sweeps natural order (no `bounded_rand_int` shuffle, `linear.cpp:629`); the CS optimum is unique so the limit is identical (documented RNG/shrink-path boundary). Pinned by `test_crammer_singer_smoke in linear_svc.rs` (live oracle 3-class set `coef [[-0.06762,-0.24341],[0.30048,0.02171],[-0.23286,0.22171]]`, `int [0.91078,-0.62206,-0.28873]`, predict all-correct; binary 8×2 collapse `coef [[0.15504,0.15504]]`, `int [-1.48062]`; within 1e-2). The rigorous oracle pin in `tests/divergence_linear_svc_fit.rs` is the critic's next step. Consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). |
//! | REQ-7 (fit_intercept + intercept_scaling) | SHIPPED | `LinearSVC<F>` exposes `pub fit_intercept: bool` (default true) + `pub intercept_scaling: F` (default 1.0) + `#[must_use]` builders. When fitting an intercept the design matrix is augmented with a penalized constant column = `intercept_scaling`, and `intercept_ = intercept_scaling·w_last` (`_base.py:1188-1198,:1240-1245`); `intercept_scaling > 0` is validated. Pinned by `linear_svc_coef_parity` + module `test_fit_intercept_false_zero_intercept`/`test_invalid_intercept_scaling`. |
//! | REQ-8 (dual param) | SHIPPED | `LinearSVC<F>` exposes `pub dual: DualMode` (default `Auto`) + `#[must_use] with_dual`. `fn resolve_dual` resolves `Auto→bool` (`_validate_dual_parameter`, `_classes.py:13-29`: `n<f`→prefer dual, else→prefer primal, with fallback) against `fn liblinear_solver_type` (the `_get_liblinear_solver_type` matrix, `_base.py:995-1018`), and `fn fit` validates the resolved combo (`hinge+dual=false`, `l1+dual=true`, `l1+hinge` all rejected → `FerroError::InvalidParameter`). R-DEV-7: the resolved `dual` is **observably immaterial for `penalty=l2`** — the l2 dual CD and l2 primal minimize the same `0.5·‖w‖² + C·Σ L` and reach the same `coef_`/`intercept_`, so `penalty=l2` keeps `fn solve_binary_dual` regardless of `dual`. Pinned by `test_unsupported_combinations_rejected` + `test_dual_auto_resolution`. Consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). |
//! | REQ-9 (class_weight) | SHIPPED | `LinearSVC<F>` exposes `pub class_weight: ClassWeight<F>` (`None`/`Balanced`/`Explicit`, default `None`) + `#[must_use] with_class_weight`. `fn compute_class_weight` (mirroring `sklearn.utils.compute_class_weight`, `class_weight.py:63-81`, as called at `_base.py:1179`) expands per-class weights; `fn fit` scales `C` per class: binary `cp = C·weights[idx(classes[1])]`, `cn = C·weights[idx(classes[0])]` (`train_one(Cp=weighted_C[1], Cn=weighted_C[0])`, `linear.cpp:2543-2551`), OvR class `k` `cp = C·weights[k]`, `cn = C` base (the negative rest is UNWEIGHTED, `linear.cpp:2559-2571`). `SolverConfig` now carries `(cp, cn)`; `solve_binary_dual`/`solve_binary_l1r_l2` apply the per-sample `C_[i] = (y_i>0?cp:cn)` (`diag[i]`/`upper_bound[i]`/`C[i]`, `linear.cpp:843-858`,`:1504-1509`). When `cp == cn` (no class_weight) the math is identical to before (the 9 divergence pins stay green). Pinned by `test_class_weight_smoke in linear_svc.rs` (live oracle 8×2 imbalanced set, `squared_hinge,dual=True,C=1.0`: `None coef [[0.10056,0.15957]] int [-1.26346]`; `balanced coef [[0.09937,0.16666]] int [-1.21320]` weights `[0.6667,2.0]`; `{0:1,1:5} coef [[0.11059,0.17164]] int [-1.29547]`; ferrolearn within 1e-2). The rigorous oracle pin in `tests/divergence_linear_svc_fit.rs` is the critic's next step. Consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). |
//! | REQ-10 (C-scaling convention) | SHIPPED | the `c / n_f` division is removed; the dual CD uses `upper_bound = C` (hinge) / `diag = 0.5/C` (squared_hinge), so `coef_` tracks `C` like liblinear. Pinned by `linear_svc_coef_c_dependence` (C=0.1 → `0.0784651864625997`, C=1.0 → `0.12835213611984458`). |
//! | REQ-11 (n_iter_/n_features_in_ + param validation) | SHIPPED | `fn n_features_in` (returns the stored `n_features`, set by `_validate_data`, `_classes.py:302`) and `fn n_iter` (the max dual-CD outer-iteration count across the binary/OvR fits, `n_iter_ = n_iter_.max().item()`, `_classes.py:338`) on `FittedLinearSVC`; `fn fit` validates `tol > 0` (`Interval(Real, 0.0, None, closed="neither")`, `_classes.py:237`). Pinned by `linear_svc_attrs_and_tol_validation in tests/divergence_linear_svc_fit.rs` (#627). `n_features_in_` (oracle `2`) and the `tol <= 0` reject are exact; `n_iter_` is the documented shuffle-path RNG boundary (ferrolearn sweeps natural order, sklearn's liblinear shuffles `index` each sweep, cf. SGD), so the pin bounds `n_iter` in `[1, max_iter]` rather than exact-matching. |
//! | REQ-12 (ferray substrate) | NOT-STARTED | open prereq blocker #628. Imports `ndarray`, not `ferray-core`/`ferray::linalg` (R-SUBSTRATE). |
//! | REQ-13 (non-finite input rejected) | SHIPPED | `fn fit` rejects any NaN/+/-inf in X BEFORE the dual coordinate-descent solve with `FerroError::InvalidParameter`, mirroring sklearn's `_validate_data(force_all_finite=True)` (`svm/_classes.py:302`, the liblinear `_fit_liblinear` path; cf. `svm/_base.py:190`) → `ValueError("Input X contains NaN.")` / `"... contains infinity ..."`. `y` is `Array1<usize>` (integer labels), finite by type, and ferrolearn's `LinearSVC::fit` takes no `sample_weight` argument, so X is the only runtime check. The guard sits ahead of both the OvR and Crammer-Singer dispatch, so it covers every solver path. `.iter().any(|v| !v.is_finite())` catches NaN and Inf; the finite path is byte-identical. Verified vs the live sklearn 1.5.2 oracle (R-CHAR-3): `LinearSVC().fit` raises `ValueError` for NaN/+inf/-inf in X (`tests/divergence_linear_nonfinite_batch4.rs::linsvc_*`). Non-test consumer: `pub use linear_svc::{…}` (`lib.rs`) + `RsLinearSVC` (PyO3). (#2263) |
use ferrolearn_core::error::FerroError;
use ferrolearn_core::introspection::{HasClasses, HasCoefficients};
use ferrolearn_core::traits::{Fit, Predict};
use ndarray::{Array1, Array2, ScalarOperand};
use num_traits::Float;
/// Penalty (regularizer) norm for [`LinearSVC`].
///
/// Mirrors `sklearn.svm.LinearSVC`'s `penalty` parameter
/// (`sklearn/svm/_classes.py:51-54`): `'l2'` (default) is the standard SVC
/// `0.5·‖w‖²` regularizer; `'l1'` is the `‖w‖₁` regularizer, which yields a
/// sparse `coef_`. The penalty interacts with [`LinearSVCLoss`] / [`DualMode`]
/// via liblinear's solver-selection matrix (`_get_liblinear_solver_type`,
/// `_base.py:1011-1018`): `'l1'` is only supported with `squared_hinge` +
/// `dual=False` (solver type 5).
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum LinearSVCPenalty {
/// `‖w‖₁` regularizer — sparse `coef_`. Only valid with `squared_hinge` +
/// `dual=false` (liblinear solver type 5, `_base.py:1014`).
L1,
/// `0.5·‖w‖²` regularizer (default), the standard SVC penalty.
#[default]
L2,
}
/// Dual / primal optimization-problem selector for [`LinearSVC`].
///
/// Mirrors `sklearn.svm.LinearSVC`'s `dual` parameter
/// (`sklearn/svm/_classes.py:62-71`, default `"auto"`). `Auto` resolves to a
/// concrete `bool` via `_validate_dual_parameter` (`_classes.py:13-29`): when
/// `n_samples < n_features` it prefers `dual=true` (falling back to `false` if
/// that penalty×loss combination has no dual solver); otherwise it prefers
/// `dual=false` (falling back to `true` if there is no primal solver, e.g.
/// `hinge`). The resolved `dual` selects the liblinear solver type
/// (`_get_liblinear_solver_type`, `_base.py:1011-1018`).
///
/// Under R-DEV-7 the resolved `dual` is **observably immaterial for
/// `penalty=l2`**: the l2 dual coordinate descent and the l2 primal both
/// minimize the same strongly convex `0.5·‖w‖² + C·Σ L` and reach the same
/// `coef_`/`intercept_`. It is load-bearing only for the unsupported-combination
/// rejects and for selecting the genuinely different `l1` primal solver.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum DualMode {
/// Resolve to `true`/`false` automatically via `_validate_dual_parameter`
/// (`_classes.py:13-29`); the default.
#[default]
Auto,
/// Solve the dual optimization problem.
True,
/// Solve the primal optimization problem.
False,
}
/// Multiclass strategy for [`LinearSVC`].
///
/// Mirrors `sklearn.svm.LinearSVC`'s `multi_class` parameter
/// (`sklearn/svm/_classes.py:239`, constraint `{"ovr", "crammer_singer"}`,
/// default `"ovr"`). `Ovr` trains one binary classifier per class (the default,
/// using the penalty/loss/dual solver matrix); `CrammerSinger` runs the joint
/// Crammer-Singer multiclass SVM solver (`Solver_MCSVM_CS`,
/// `liblinear/linear.cpp:493-787`, solver type 4).
///
/// When `CrammerSinger` is selected, `_get_liblinear_solver_type` returns 4
/// **regardless of penalty/loss/dual** (`_base.py:1017,1020-1021`), so the
/// penalty/loss/dual parameters are ignored and the joint solver runs.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum MultiClass {
/// One-vs-rest: one binary sub-problem per class (the default).
#[default]
Ovr,
/// Joint Crammer-Singer multiclass SVM (`Solver_MCSVM_CS`, solver type 4,
/// `linear.cpp:493-787`). Ignores penalty/loss/dual (`_base.py:1017`).
CrammerSinger,
}
/// Loss function for [`LinearSVC`].
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum LinearSVCLoss {
/// Standard hinge loss: `max(0, 1 - y * f(x))`. liblinear solver
/// `L2R_L1LOSS_SVC_DUAL` (type 3): box `0 <= alpha <= C`, `diag = 0`.
Hinge,
/// Squared hinge loss: `max(0, 1 - y * f(x))^2` (default). liblinear solver
/// `L2R_L2LOSS_SVC_DUAL` (type 1): box `0 <= alpha <= +inf`,
/// `diag = 0.5 / C`.
SquaredHinge,
}
/// Per-class weighting strategy for [`LinearSVC`].
///
/// Mirrors `sklearn.svm.LinearSVC`'s `class_weight` parameter
/// (`sklearn/svm/_classes.py:118-124`, constraint `{None, dict, 'balanced'}`):
/// it scales the inverse-regularization `C` per class so the effective penalty
/// for class `i` is `class_weight[i]·C` (`compute_class_weight`,
/// `sklearn/svm/_base.py:1179`; `weighted_C[i] = C·class_weight[i]`,
/// `liblinear/linear.cpp:2496-2507`). The expanded per-class weights are
/// computed by [`compute_class_weight`] following
/// `sklearn.utils.compute_class_weight` semantics
/// (`sklearn/utils/class_weight.py:63-81`).
///
/// This mirrors `ferrolearn_linear::sgd::ClassWeight` for cross-estimator
/// consistency, but is defined locally (no cross-import of `sgd` internals).
#[derive(Debug, Clone, Default)]
pub enum ClassWeight<F> {
/// Uniform weights (all classes weighted `1.0`). The default
/// (`class_weight=None`, `class_weight.py:63-65`).
#[default]
None,
/// Balanced weights `n_samples / (n_classes · count_c)` per class `c`,
/// matching `sklearn.utils.compute_class_weight("balanced", ...)`
/// (`class_weight.py:66-74`).
Balanced,
/// Explicit class-label -> weight map. Classes absent from the map default
/// to `1.0`, matching the dict branch of `compute_class_weight`
/// (`class_weight.py:75-81`).
Explicit(Vec<(usize, F)>),
}
/// Compute the expanded per-class weight vector aligned to `classes`
/// (sorted ascending, matching sklearn's `classes_ = np.unique(y)`).
///
/// Faithful to `sklearn.utils.compute_class_weight`
/// (`sklearn/utils/class_weight.py:63-81`), as called by `_fit_liblinear`
/// (`compute_class_weight(class_weight, classes=classes_, y=y)`,
/// `sklearn/svm/_base.py:1179`):
/// - `None` -> all `1.0` (`:63-65`).
/// - `Balanced` -> `n_samples / (n_classes · count_c)` per class `c`,
/// where `count_c` is the number of samples with label `c` (`:66-74`).
/// - `Explicit(map)` -> `1.0` default, overridden by the map entries matched by
/// class label (`:75-81`).
///
/// `classes` is the sorted unique label set; `y` is the per-sample label array.
/// Mirrors `ferrolearn_linear::sgd::compute_class_weight` exactly.
fn compute_class_weight<F: Float>(cw: &ClassWeight<F>, classes: &[usize], y: &[usize]) -> Vec<F> {
match cw {
ClassWeight::None => vec![F::one(); classes.len()],
ClassWeight::Balanced => {
// `recip_freq = len(y) / (n_classes * bincount(y_ind))`
// (`class_weight.py:73`), indexed per class.
let n_samples = F::from(y.len()).unwrap_or_else(F::zero);
let n_classes = F::from(classes.len()).unwrap_or_else(F::one);
classes
.iter()
.map(|&c| {
let count = y.iter().filter(|&&label| label == c).count();
let count_f = F::from(count).unwrap_or_else(F::one);
if count_f > F::zero() {
n_samples / (n_classes * count_f)
} else {
F::one()
}
})
.collect()
}
ClassWeight::Explicit(map) => classes
.iter()
.map(|&c| {
map.iter()
.find(|(label, _)| *label == c)
.map_or_else(F::one, |(_, w)| *w)
})
.collect(),
}
}
/// Confidence scores returned by [`FittedLinearSVC::decision_function`].
///
/// Mirrors `sklearn.svm.LinearSVC.decision_function`
/// (`LinearClassifierMixin.decision_function`,
/// `sklearn/linear_model/_base.py:341-365`): the binary case collapses the
/// single-column score matrix to a 1-D `(n_samples,)` array
/// (`return xp.reshape(scores, (-1,)) if scores.shape[1] == 1 else scores`,
/// `_base.py:365`), the multiclass case returns `(n_samples, n_classes)`.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum DecisionScores<F> {
/// Binary scores `X · w + b` for the positive class (`classes_[1]`), shape
/// `(n_samples,)`. `> 0` predicts `classes_[1]`.
Binary(Array1<F>),
/// One-vs-rest scores, shape `(n_samples, n_classes)`. The argmax of each
/// row agrees with [`Predict`].
Multiclass(Array2<F>),
}
impl<F: Clone> DecisionScores<F> {
/// Number of samples scored (the leading axis length in both variants).
#[must_use]
pub fn n_samples(&self) -> usize {
match self {
DecisionScores::Binary(v) => v.len(),
DecisionScores::Multiclass(m) => m.nrows(),
}
}
/// Borrow the binary 1-D scores, if this is the binary case.
#[must_use]
pub fn as_binary(&self) -> Option<&Array1<F>> {
match self {
DecisionScores::Binary(v) => Some(v),
DecisionScores::Multiclass(_) => None,
}
}
/// Borrow the multiclass `(n_samples, n_classes)` scores, if this is the
/// multiclass case.
#[must_use]
pub fn as_multiclass(&self) -> Option<&Array2<F>> {
match self {
DecisionScores::Multiclass(m) => Some(m),
DecisionScores::Binary(_) => None,
}
}
}
/// Linear Support Vector Classifier (liblinear dual CD).
///
/// Solves the L2-regularized hinge or squared-hinge objective
/// `0.5*||w||^2 + C * sum_i L(y_i, w.x_i)` via liblinear's dual coordinate
/// descent. Supports binary and multiclass (one-vs-rest) classification.
/// Mirrors `sklearn.svm.LinearSVC`.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
#[derive(Debug, Clone)]
pub struct LinearSVC<F> {
/// Inverse regularization strength. Larger values allow more
/// misclassification. Must be strictly positive.
pub c: F,
/// Maximum number of dual coordinate descent iterations.
pub max_iter: usize,
/// Convergence tolerance on the projected-gradient span.
pub tol: F,
/// Loss function to use.
pub loss: LinearSVCLoss,
/// Regularizer norm (`l1` / `l2`). Default `l2`. `l1` is only supported with
/// `squared_hinge` + `dual=false` (liblinear solver type 5,
/// `_get_liblinear_solver_type`, `_base.py:1014`).
pub penalty: LinearSVCPenalty,
/// Dual / primal selector. Default `Auto`, resolved via
/// `_validate_dual_parameter` (`_classes.py:13-29`). Observably immaterial
/// for `penalty=l2` (R-DEV-7 dual-invariance), load-bearing for the
/// unsupported-combination rejects and the `l1` primal solver.
pub dual: DualMode,
/// Whether to fit an intercept. When `true`, the design matrix is augmented
/// with a synthetic constant column equal to [`intercept_scaling`]; the
/// augmented weight is penalized like any feature (liblinear convention).
///
/// [`intercept_scaling`]: Self::intercept_scaling
pub fit_intercept: bool,
/// Value of the synthetic intercept feature when [`fit_intercept`] is
/// `true`. Must be strictly positive. `intercept_ = intercept_scaling *
/// w_last`.
///
/// [`fit_intercept`]: Self::fit_intercept
pub intercept_scaling: F,
/// Per-class scaling of `C`. Default [`ClassWeight::None`] (all classes
/// weighted `1.0`). The effective penalty for class `i` is
/// `class_weight[i]·C` (`compute_class_weight`, `_base.py:1179`;
/// `weighted_C[i] = C·class_weight[i]`, `linear.cpp:2496-2507`).
pub class_weight: ClassWeight<F>,
/// Multiclass strategy. Default [`MultiClass::Ovr`] (one-vs-rest).
/// [`MultiClass::CrammerSinger`] runs the joint Crammer-Singer solver,
/// ignoring penalty/loss/dual (`_base.py:1017`, `_classes.py:239`).
pub multi_class: MultiClass,
}
impl<F: Float> LinearSVC<F> {
/// Create a new `LinearSVC` with scikit-learn's default settings.
///
/// Defaults (matching `sklearn.svm.LinearSVC`, `_classes.py`):
/// `C = 1.0`, `max_iter = 1000`, `tol = 1e-4`, `loss = SquaredHinge`,
/// `penalty = L2`, `dual = Auto`, `fit_intercept = true`,
/// `intercept_scaling = 1.0`, `class_weight = None`, `multi_class = Ovr`.
#[must_use]
pub fn new() -> Self {
// 1e-4/1.0 are exactly representable in f32/f64; the defensive fallback
// for `from(1e-4)` is never taken (no `.unwrap()` in lib code).
let one = F::one();
Self {
c: one,
max_iter: 1000,
tol: F::from(1e-4).unwrap_or_else(|| {
let ten = F::from(10).unwrap_or(one);
one / (ten * ten * ten * ten)
}),
loss: LinearSVCLoss::SquaredHinge,
penalty: LinearSVCPenalty::L2,
dual: DualMode::Auto,
fit_intercept: true,
intercept_scaling: one,
class_weight: ClassWeight::None,
multi_class: MultiClass::Ovr,
}
}
/// Set the regularization parameter C.
#[must_use]
pub fn with_c(mut self, c: F) -> Self {
self.c = c;
self
}
/// Set the maximum number of iterations.
#[must_use]
pub fn with_max_iter(mut self, max_iter: usize) -> Self {
self.max_iter = max_iter;
self
}
/// Set the convergence tolerance.
#[must_use]
pub fn with_tol(mut self, tol: F) -> Self {
self.tol = tol;
self
}
/// Set the loss function.
#[must_use]
pub fn with_loss(mut self, loss: LinearSVCLoss) -> Self {
self.loss = loss;
self
}
/// Set the penalty (regularizer) norm (sklearn `penalty`). `l1` requires
/// `squared_hinge` + `dual=false` (`_base.py:1014`).
#[must_use]
pub fn with_penalty(mut self, penalty: LinearSVCPenalty) -> Self {
self.penalty = penalty;
self
}
/// Set the dual / primal selector (sklearn `dual`). `Auto` (default)
/// resolves via `_validate_dual_parameter` (`_classes.py:13-29`).
#[must_use]
pub fn with_dual(mut self, dual: DualMode) -> Self {
self.dual = dual;
self
}
/// Set whether to fit an intercept (sklearn `fit_intercept`).
#[must_use]
pub fn with_fit_intercept(mut self, fit_intercept: bool) -> Self {
self.fit_intercept = fit_intercept;
self
}
/// Set the intercept scaling (sklearn `intercept_scaling`). Must be
/// strictly positive when `fit_intercept` is `true`.
#[must_use]
pub fn with_intercept_scaling(mut self, intercept_scaling: F) -> Self {
self.intercept_scaling = intercept_scaling;
self
}
/// Set the per-class `C` scaling (sklearn `class_weight`,
/// `_classes.py:118-124`). [`ClassWeight::None`] (default) leaves every
/// class at `1.0`; [`ClassWeight::Balanced`] uses
/// `n_samples / (n_classes · count_c)`; [`ClassWeight::Explicit`] takes a
/// `(label, weight)` map (unlisted classes default to `1.0`).
#[must_use]
pub fn with_class_weight(mut self, class_weight: ClassWeight<F>) -> Self {
self.class_weight = class_weight;
self
}
/// Set the multiclass strategy (sklearn `multi_class`, `_classes.py:239`).
/// [`MultiClass::Ovr`] (default) trains one-vs-rest binary classifiers;
/// [`MultiClass::CrammerSinger`] runs the joint Crammer-Singer solver
/// (ignoring penalty/loss/dual, `_base.py:1017`).
#[must_use]
pub fn with_multi_class(mut self, multi_class: MultiClass) -> Self {
self.multi_class = multi_class;
self
}
}
impl<F: Float> Default for LinearSVC<F> {
fn default() -> Self {
Self::new()
}
}
/// Fitted Linear Support Vector Classifier.
///
/// Stores the learned weight vectors, intercepts, and class labels.
/// For binary classification a single weight vector is stored; for
/// multiclass, one per class (one-vs-rest).
#[derive(Debug, Clone)]
pub struct FittedLinearSVC<F> {
/// Weight vectors: one per binary sub-problem.
/// Binary: `[w]`, Multiclass: `[w_0, w_1, ..., w_{k-1}]`.
weight_vectors: Vec<Array1<F>>,
/// Intercept for each sub-problem.
intercepts: Vec<F>,
/// Sorted unique class labels.
classes: Vec<usize>,
/// Whether this is a binary problem.
is_binary: bool,
/// Number of features.
n_features: usize,
/// Maximum dual-CD outer-iteration count across sub-problem fits
/// (`n_iter_ = n_iter_.max().item()`, `_classes.py:338`).
n_iter: usize,
}
impl<F: Float> FittedLinearSVC<F> {
/// Returns the weight vectors (one per binary sub-problem).
#[must_use]
pub fn weight_vectors(&self) -> &[Array1<F>] {
&self.weight_vectors
}
/// Returns the intercepts (one per binary sub-problem).
#[must_use]
pub fn intercepts(&self) -> &[F] {
&self.intercepts
}
/// Number of features seen during fit (`n_features_in_`).
///
/// Mirrors sklearn's `n_features_in_`, set by `_validate_data`
/// (`sklearn/svm/_classes.py:302`); equals `X.ncols()`.
#[must_use]
pub fn n_features_in(&self) -> usize {
self.n_features
}
/// Maximum number of dual coordinate-descent outer iterations across the
/// (binary or one-vs-rest) sub-problem fits.
///
/// Mirrors sklearn's `n_iter_ = n_iter_.max().item()`
/// (`sklearn/svm/_classes.py:338`). The exact value is shuffle-path
/// dependent (sklearn's liblinear shuffles the active index each sweep;
/// ferrolearn sweeps natural order), so it is bounded in `[1, max_iter]`
/// rather than exact-matching the oracle.
#[must_use]
pub fn n_iter(&self) -> usize {
self.n_iter
}
}
impl<F: Float + ScalarOperand + Send + Sync + 'static> FittedLinearSVC<F> {
/// Raw signed distance from the decision boundary. Mirrors sklearn
/// `LinearSVC.decision_function`.
///
/// Binary: [`DecisionScores::Binary`] of shape `(n_samples,)` containing
/// `X @ w + b` for the positive class (`classes_[1]`); sklearn ravels the
/// single-column score to 1-D (`linear_model/_base.py:365`).
/// Multiclass: [`DecisionScores::Multiclass`] of shape
/// `(n_samples, n_classes)` of one-vs-rest scores; the argmax of each row
/// agrees with [`Predict`].
///
/// # Errors
///
/// Returns [`FerroError::ShapeMismatch`] if the number of features
/// does not match the fitted model.
pub fn decision_function(&self, x: &Array2<F>) -> Result<DecisionScores<F>, FerroError> {
let n_features = x.ncols();
if n_features != self.n_features {
return Err(FerroError::ShapeMismatch {
expected: vec![self.n_features],
actual: vec![n_features],
context: "number of features must match fitted model".into(),
});
}
let n_samples = x.nrows();
if self.is_binary {
// sklearn collapses the single-column binary score matrix to a 1-D
// (n_samples,) array (`linear_model/_base.py:365`).
let scores = x.dot(&self.weight_vectors[0]) + self.intercepts[0];
Ok(DecisionScores::Binary(scores))
} else {
let n_classes = self.classes.len();
let mut out = Array2::<F>::zeros((n_samples, n_classes));
for c in 0..n_classes {
for i in 0..n_samples {
out[[i, c]] = x.row(i).dot(&self.weight_vectors[c]) + self.intercepts[c];
}
}
Ok(DecisionScores::Multiclass(out))
}
}
}
/// Resolved per-fit solver configuration passed to [`solve_binary_dual`]
/// (groups the dual-CD knobs to keep the solver signature small).
#[derive(Debug, Clone, Copy)]
struct SolverConfig<F> {
/// Per-sample penalty for the positive (`y_i > 0`) group: `Cp = C·w[+]`
/// (`train_one(Cp, Cn)`, `linear.cpp:2543-2571`; `C_[i] = (y_i>0 ? Cp :
/// Cn)`, `linear.cpp:843-858`, `:1504-1509`).
cp: F,
/// Per-sample penalty for the negative (`y_i <= 0`) group: `Cn = C·w[-]`
/// (binary) or the base `C` (multiclass OvR; the negative group is the
/// unweighted rest, `linear.cpp:2559-2571`).
cn: F,
/// Maximum dual-CD outer iterations.
max_iter: usize,
/// Projected-gradient span stopping tolerance.
tol: F,
/// Hinge / squared-hinge loss selector.
loss: LinearSVCLoss,
/// Whether the design matrix is augmented with a penalized bias column.
fit_intercept: bool,
/// Value of the synthetic bias column when `fit_intercept` is set.
intercept_scaling: F,
}
/// Solve a single binary L2-SVM via liblinear's dual coordinate descent.
///
/// Minimizes `0.5 * ||w||^2 + C * sum_i L(y_i, w.x_i)` (NO `1/n` averaging,
/// matching `sklearn/svm/_base.py` `_fit_liblinear`) with `y_i ∈ {-1, +1}`.
/// The augmented weight vector has length `w_size`; when `fit_intercept` is
/// set, `w_size = n_features + 1` and the trailing weight multiplies the
/// synthetic constant column `intercept_scaling` (penalized like any feature).
///
/// This is liblinear's `solve_l2r_l1l2_svc` (`linear.cpp:819`): the dual is
///
/// ```text
/// min_alpha 0.5 * alpha^T (Q + diag) alpha - e^T alpha,
/// s.t. 0 <= alpha_i <= U,
/// ```
///
/// where `Q_ij = y_i y_j x_i.x_j`, `w = sum_i alpha_i y_i x_i`,
/// `QD[i] = diag + ||x_i||^2`. **hinge** (`L2R_L1LOSS_SVC_DUAL`): `diag = 0`,
/// `U = C` (`linear.cpp:852-858`). **squared_hinge** (`L2R_L2LOSS_SVC_DUAL`):
/// `diag = 0.5/C`, `U = +inf` (`linear.cpp:849-850`).
///
/// Returns `(w_augmented, n_iter, converged)`.
fn solve_binary_dual<F: Float + 'static>(
x: &Array2<F>,
y_signed: &Array1<F>,
cfg: &SolverConfig<F>,
) -> (Vec<F>, usize, bool) {
let SolverConfig {
cp,
cn,
max_iter,
tol,
loss,
fit_intercept,
intercept_scaling,
} = *cfg;
let (n_samples, n_features) = x.dim();
let w_size = if fit_intercept {
n_features + 1
} else {
n_features
};
let inf = F::infinity();
let two = F::one() + F::one();
let half = F::one() / two;
let tiny = F::from(1.0e-12).unwrap_or_else(F::epsilon);
// Per-sample penalty `C_[i] = (y_i > 0 ? Cp : Cn)` (`linear.cpp:843-858`,
// `GETI(i) ≡ i`); `class_weight` makes Cp/Cn differ. Per-sample diag /
// upper_bound follow the solver type: squared_hinge → `diag[i] = 0.5/C_[i]`,
// `U[i] = +inf`; hinge → `diag[i] = 0`, `U[i] = C_[i]`.
let mut diag = vec![F::zero(); n_samples];
let mut upper_bound = vec![inf; n_samples];
for i in 0..n_samples {
let c_i = if y_signed[i] > F::zero() { cp } else { cn };
match loss {
LinearSVCLoss::Hinge => {
diag[i] = F::zero();
upper_bound[i] = c_i;
}
LinearSVCLoss::SquaredHinge => {
diag[i] = half / c_i;
upper_bound[i] = inf;
}
}
}
// QD[i] = diag[i] + ||x_i||^2 (including the augmented bias column).
let mut qd = vec![F::zero(); n_samples];
for (i, qd_i) in qd.iter_mut().enumerate() {
let mut acc = diag[i];
let row = x.row(i);
for &v in row.iter() {
acc = acc + v * v;
}
if fit_intercept {
acc = acc + intercept_scaling * intercept_scaling;
}
*qd_i = acc;
}
let mut alpha = vec![F::zero(); n_samples];
let mut w = vec![F::zero(); w_size];
// alpha starts at 0 so w starts at 0; nothing to accumulate.
let mut index: Vec<usize> = (0..n_samples).collect();
let mut active_size = n_samples;
let mut pgmax_old = inf;
let mut pgmin_old = -inf;
// w . x_i over the (augmented) design matrix.
let dot_w_xi = |w: &[F], i: usize| -> F {
let mut acc = F::zero();
let row = x.row(i);
for (j, &v) in row.iter().enumerate() {
acc = acc + v * w[j];
}
if fit_intercept {
acc = acc + intercept_scaling * w[n_features];
}
acc
};
let mut n_iter: usize = 0;
let mut converged = false;
// liblinear shuffles `index` each sweep; the minimizer is unique so order
// only affects the path, not the limit. We sweep in natural order for
// determinism (no RNG), reaching the same converged optimum.
for iter in 0..max_iter {
n_iter = iter + 1;
let mut pgmax_new = -inf;
let mut pgmin_new = inf;
let mut s = 0;
while s < active_size {
let i = index[s];
let yi = y_signed[i];
// G = y_i*(w.x_i) - 1 + diag_i*alpha_i (`linear.cpp:909-921`).
let g = yi * dot_w_xi(&w, i) - F::one() + diag[i] * alpha[i];
// Projected gradient + shrinking (`linear.cpp:923-949`).
let mut pg = F::zero();
if alpha[i] == F::zero() {
if g > pgmax_old {
active_size -= 1;
index.swap(s, active_size);
continue; // re-process the swapped-in element at `s`
} else if g < F::zero() {
pg = g;
}
} else if alpha[i] == upper_bound[i] {
if g < pgmin_old {
active_size -= 1;
index.swap(s, active_size);
continue;
} else if g > F::zero() {
pg = g;
}
} else {
pg = g;
}
if pg > pgmax_new {
pgmax_new = pg;
}
if pg < pgmin_new {
pgmin_new = pg;
}
if pg.abs() > tiny {
let alpha_old = alpha[i];
// alpha_i <- clamp(alpha_i - G/QD[i], 0, U) (`linear.cpp:957`).
let mut new_alpha = alpha[i] - g / qd[i];
if new_alpha < F::zero() {
new_alpha = F::zero();
} else if new_alpha > upper_bound[i] {
new_alpha = upper_bound[i];
}
alpha[i] = new_alpha;
let d = (alpha[i] - alpha_old) * yi;
if d != F::zero() {
let row = x.row(i);
for (j, &v) in row.iter().enumerate() {
w[j] = w[j] + d * v;
}
if fit_intercept {
w[n_features] = w[n_features] + d * intercept_scaling;
}
}
}
s += 1;
}
// Stopping: PGmax_new - PGmin_new <= tol on the full set
// (`linear.cpp:972-990`). Absolute tol (the dual SVC solver receives
// `param->eps` directly, `linear.cpp:2364`).
if pgmax_new - pgmin_new <= tol {
if active_size == n_samples {
converged = true;
break;
}
active_size = n_samples;
pgmax_old = inf;
pgmin_old = -inf;
continue;
}
pgmax_old = pgmax_new;
pgmin_old = pgmin_new;
if pgmax_old <= F::zero() {
pgmax_old = inf;
}
if pgmin_old >= F::zero() {
pgmin_old = -inf;
}
}
(w, n_iter, converged)
}
/// Resolve `(penalty, loss, dual)` to a liblinear solver "magic number" for the
/// one-vs-rest `multi_class='ovr'` case, mirroring `_get_liblinear_solver_type`
/// (`sklearn/svm/_base.py:995-1049`). Returns `Err` for the unsupported
/// combinations sklearn raises `ValueError` on.
///
/// The supported (`multi_class='ovr'`) entries of `_solver_type_dict`
/// (`_base.py:1013-1014`):
///
/// ```text
/// hinge: { l2: { dual=true: 3 } }
/// squared_hinge: { l1: { dual=false: 5 }, l2: { dual=false: 2, dual=true: 1 } }
/// ```
///
/// (`crammer_singer` = 4 is REQ-6/#623, out of scope; this function assumes the
/// existing `'ovr'` multi-class.) The `Err` strings mirror sklearn's
/// `error_string` (`_base.py:1033-1043`).
fn liblinear_solver_type(
penalty: LinearSVCPenalty,
loss: LinearSVCLoss,
dual: bool,
) -> Result<u8, FerroError> {
match (loss, penalty, dual) {
// hinge: { l2: { dual=true: 3 } }
(LinearSVCLoss::Hinge, LinearSVCPenalty::L2, true) => Ok(3),
(LinearSVCLoss::Hinge, LinearSVCPenalty::L2, false) => Err(FerroError::InvalidParameter {
name: "dual".into(),
reason: "The combination of penalty='l2' and loss='hinge' are not \
supported when dual=false"
.into(),
}),
// hinge + l1: penalty has no entry under `hinge` → combination unsupported.
(LinearSVCLoss::Hinge, LinearSVCPenalty::L1, _) => Err(FerroError::InvalidParameter {
name: "penalty".into(),
reason: "The combination of penalty='l1' and loss='hinge' is not supported".into(),
}),
// squared_hinge: { l1: { dual=false: 5 } }
(LinearSVCLoss::SquaredHinge, LinearSVCPenalty::L1, false) => Ok(5),
(LinearSVCLoss::SquaredHinge, LinearSVCPenalty::L1, true) => {
Err(FerroError::InvalidParameter {
name: "dual".into(),
reason: "The combination of penalty='l1' and loss='squared_hinge' are not \
supported when dual=true"
.into(),
})
}
// squared_hinge: { l2: { dual=false: 2, dual=true: 1 } }
(LinearSVCLoss::SquaredHinge, LinearSVCPenalty::L2, false) => Ok(2),
(LinearSVCLoss::SquaredHinge, LinearSVCPenalty::L2, true) => Ok(1),
}
}
/// Resolve the [`DualMode`] to a concrete `bool`, mirroring
/// `_validate_dual_parameter` (`sklearn/svm/_classes.py:13-29`).
///
/// For [`DualMode::Auto`]: when `n_samples < n_features` try `dual=true` (fall
/// back to `false` if that penalty×loss combination has no dual solver); else
/// (`n_samples >= n_features`) try `dual=false` (fall back to `true` if there is
/// no primal solver, e.g. `hinge`). Resolution is checked against
/// [`liblinear_solver_type`] so it is automatically consistent with the
/// solver-selection matrix.
fn resolve_dual(
dual: DualMode,
penalty: LinearSVCPenalty,
loss: LinearSVCLoss,
n_samples: usize,
n_features: usize,
) -> bool {
match dual {
DualMode::True => true,
DualMode::False => false,
DualMode::Auto => {
if n_samples < n_features {
// Prefer dual=true; fall back to false if unsupported.
liblinear_solver_type(penalty, loss, true).is_ok()
} else {
// Prefer dual=false; fall back to true if no primal solver.
liblinear_solver_type(penalty, loss, false).is_err()
}
}
}
}
/// Solve a single binary L1-regularized L2-loss (squared-hinge) SVM via
/// liblinear's feature-major coordinate descent, `solve_l1r_l2_svc`
/// (`sklearn/svm/src/liblinear/linear.cpp:1467`). Minimizes
///
/// ```text
/// ‖w‖₁ + C · Σ_i max(0, 1 − y_i·(w·x_i))²
/// ```
///
/// (the `l1`-penalty objective — a genuinely different, sparse optimum from the
/// l2 dual). The augmented intercept column (value `intercept_scaling`) is
/// appended when `fit_intercept`, penalized in `‖w‖₁` like any feature
/// (`coef_ = w[:n_features]`, `intercept_ = intercept_scaling·w[n_features]`).
///
/// State (`linear.cpp:1488-1526`): `b[i] = 1 − y_i·(w·x_i)` (running residual),
/// `xj_sq[j] = Σ_i C·(y_i·x_ij)²`. Per feature `j`:
/// `G_loss = −2·Σ_{i: b[i]>0} C·(y_i·x_ij)·b[i]`,
/// `H = max(2·Σ_{i: b[i]>0} C·(y_i·x_ij)², 1e-12)` (`linear.cpp:1542-1562`).
/// `Gp = G_loss+1`, `Gn = G_loss−1`. Newton direction with soft-threshold
/// (`linear.cpp:1589-1595`): `d = −Gp/H` if `Gp < H·w[j]`, `d = −Gn/H` if
/// `Gn > H·w[j]`, else `d = −w[j]`. Then a backtracking line search
/// (`sigma=0.01`, ≤20 steps, halving `d`, `linear.cpp:1600-1661`) updating
/// `b[]`, then `w[j] += d`. Shrinking via `active_size`/`Gmax_old`
/// (`linear.cpp:1567-1579, 1691-1705`); stop when
/// `Gnorm1_new ≤ eps·Gnorm1_init` on the full active set (`linear.cpp:1691`).
///
/// `C[i] = (y_i > 0 ? Cp : Cn)` per-sample (`class_weight` scales `C` per class,
/// `linear.cpp:1504-1509`; `Cp = C·w[+]`, `Cn = C·w[-]` binary / base `C`
/// multiclass, `:2543-2571`). liblinear shuffles
/// `index` each sweep (`bounded_rand_int`, `linear.cpp:1535`); ferrolearn sweeps
/// NATURAL ORDER for determinism (no RNG) — the l1 optimum is unique so the
/// limit is identical (the documented RNG-path boundary, as `solve_binary_dual`).
/// We use `eps = tol` directly: liblinear scales `primal_solver_tol`
/// (`linear.cpp:2321,2374`) but the unique optimum is `tol`-invariant at the
/// limit (the test drives `tol=1e-10` + huge `max_iter`).
///
/// Returns `(w_augmented, n_iter, converged)`.
#[allow(
clippy::too_many_lines,
reason = "faithful transcription of liblinear solve_l1r_l2_svc (linear.cpp:1467)"
)]
fn solve_binary_l1r_l2<F: Float + 'static>(
x: &Array2<F>,
y_signed: &Array1<F>,
cfg: &SolverConfig<F>,
) -> (Vec<F>, usize, bool) {
let SolverConfig {
cp,
cn,
max_iter,
tol,
fit_intercept,
intercept_scaling,
..
} = *cfg;
let (n_samples, n_features) = x.dim();
// Per-sample penalty `C[i] = (y_i > 0 ? Cp : Cn)` (`solve_l1r_l2_svc`,
// `linear.cpp:1504-1509`); `class_weight` makes Cp/Cn differ.
let c_of = |i: usize| -> F { if y_signed[i] > F::zero() { cp } else { cn } };
let w_size = if fit_intercept {
n_features + 1
} else {
n_features
};
let inf = F::infinity();
let two = F::one() + F::one();
let sigma = F::from(0.01).unwrap_or_else(|| F::one() / (two * two * two * two * two * two));
let tiny = F::from(1.0e-12).unwrap_or_else(F::epsilon);
let max_num_linesearch = 20usize;
let nl = F::from(n_samples).unwrap_or_else(F::one);
// `yx(i, j)` = y_i·x_ij over the augmented design matrix (the j-th feature
// column entry for sample i). liblinear stores `x->value *= y[ind]`
// (`linear.cpp:1520`); we recompute it lazily for determinism / clarity.
let yx = |i: usize, j: usize| -> F {
let yi = y_signed[i];
if j < n_features {
yi * x[[i, j]]
} else {
yi * intercept_scaling
}
};
// b[i] = 1 − y_i·(w·x_i). w starts at 0 so b starts at 1 (`linear.cpp:1500`).
let mut b = vec![F::one(); n_samples];
let mut w = vec![F::zero(); w_size];
// xj_sq[j] = Σ_i C[i]·(y_i·x_ij)² (`linear.cpp:1523`, per-sample C[i]).
let mut xj_sq = vec![F::zero(); w_size];
for (j, xj_sq_j) in xj_sq.iter_mut().enumerate() {
let mut acc = F::zero();
for i in 0..n_samples {
let val = yx(i, j);
acc = acc + c_of(i) * val * val;
}
*xj_sq_j = acc;
}
let mut index: Vec<usize> = (0..w_size).collect();
let mut active_size = w_size;
let mut gmax_old = inf;
let mut gnorm1_init = -F::one();
let mut n_iter: usize = 0;
let mut converged = false;
while n_iter < max_iter {
let mut gmax_new = F::zero();
let mut gnorm1_new = F::zero();
// liblinear shuffles `index[0..active_size]` here (`linear.cpp:1533-1537`);
// ferrolearn sweeps natural order (no RNG); the unique optimum is
// unchanged at the limit.
let mut s = 0;
while s < active_size {
let j = index[s];
// G_loss = −2·Σ_{i: b[i]>0} C·(y_i·x_ij)·b[i];
// H = 2·Σ_{i: b[i]>0} C·(y_i·x_ij)² (`linear.cpp:1542-1561`).
let mut g_loss = F::zero();
let mut h = F::zero();
for (i, &bi) in b.iter().enumerate() {
if bi > F::zero() {
let val = yx(i, j);
let tmp = c_of(i) * val;
g_loss = g_loss - tmp * bi;
h = h + tmp * val;
}
}
g_loss = g_loss * two;
let g = g_loss;
h = h * two;
if h < tiny {
h = tiny;
}
let gp = g + F::one();
let gn = g - F::one();
let wj = w[j];
// Violation + shrinking (`linear.cpp:1564-1587`).
let mut violation = F::zero();
if wj == F::zero() {
if gp < F::zero() {
violation = -gp;
} else if gn > F::zero() {
violation = gn;
} else if gp > gmax_old / nl && gn < -(gmax_old / nl) {
active_size -= 1;
index.swap(s, active_size);
continue; // re-process the swapped-in element at `s`
}
} else if wj > F::zero() {
violation = gp.abs();
} else {
violation = gn.abs();
}
if violation > gmax_new {
gmax_new = violation;
}
gnorm1_new = gnorm1_new + violation;
// Newton direction with soft-threshold (`linear.cpp:1589-1595`).
let mut d = if gp < h * wj {
-gp / h
} else if gn > h * wj {
-gn / h
} else {
-wj
};
if d.abs() < tiny {
s += 1;
continue;
}
// Backtracking line search (`linear.cpp:1600-1661`).
let mut delta = (wj + d).abs() - wj.abs() + g * d;
let mut d_old = F::zero();
let mut num_linesearch = 0usize;
while num_linesearch < max_num_linesearch {
let d_diff = d_old - d;
let mut cond = (wj + d).abs() - wj.abs() - sigma * delta;
let appxcond = xj_sq[j] * d * d + g_loss * d + cond;
if appxcond <= F::zero() {
for (i, bi) in b.iter_mut().enumerate() {
*bi = *bi + d_diff * yx(i, j);
}
break;
}
let mut loss_old = F::zero();
let mut loss_new = F::zero();
if num_linesearch == 0 {
for (i, bi) in b.iter_mut().enumerate() {
if *bi > F::zero() {
loss_old = loss_old + c_of(i) * *bi * *bi;
}
let b_new = *bi + d_diff * yx(i, j);
*bi = b_new;
if b_new > F::zero() {
loss_new = loss_new + c_of(i) * b_new * b_new;
}
}
} else {
for (i, bi) in b.iter_mut().enumerate() {
let b_new = *bi + d_diff * yx(i, j);
*bi = b_new;
if b_new > F::zero() {
loss_new = loss_new + c_of(i) * b_new * b_new;
}
}
}
cond = cond + loss_new - loss_old;
if cond <= F::zero() {
break;
}
d_old = d;
d = d / two;
delta = delta / two;
num_linesearch += 1;
}
w[j] = w[j] + d;
// Recompute b[] if the line search took the maximum steps
// (`linear.cpp:1665-1682`).
if num_linesearch >= max_num_linesearch {
for bi in b.iter_mut() {
*bi = F::one();
}
for (jj, &wjj) in w.iter().enumerate() {
if wjj == F::zero() {
continue;
}
for (i, bi) in b.iter_mut().enumerate() {
*bi = *bi - wjj * yx(i, jj);
}
}
}
s += 1;
}
if n_iter == 0 {
gnorm1_init = gnorm1_new;
}
n_iter += 1;
// Stop when Gnorm1_new ≤ eps·Gnorm1_init on the full active set
// (`linear.cpp:1691-1702`).
if gnorm1_new <= tol * gnorm1_init {
if active_size == w_size {
converged = true;
break;
}
active_size = w_size;
gmax_old = inf;
continue;
}
gmax_old = gmax_new;
}
(w, n_iter, converged)
}
/// Solve the joint Crammer-Singer multiclass SVM, transcribing liblinear's
/// `Solver_MCSVM_CS` (`sklearn/svm/src/liblinear/linear.cpp:493-787`,
/// solver type 4, `_base.py:1017`).
///
/// Minimizes the Crammer-Singer objective over a single joint weight matrix
/// `w` flattened as `(w_size × nr_class)` with `w[feature*nr_class + m]`. The
/// dual variables `alpha[i*nr_class + m]` satisfy `Σ_m alpha[i,m] = 0`,
/// `alpha[i,m] <= C[i]` if `y_i == m` else `alpha[i,m] <= 0`. Per sample,
/// `C[i] = W[i] · weighted_C[y_i]` (`linear.cpp:521-522`); here `W[i] = 1` and
/// `weighted_C[c] = C · class_weight[c]` (`= C` when `class_weight=None`).
///
/// `y_class[i]` is the class index (`0..nr_class`, the sorted-`classes`
/// position) of sample `i`. When `fit_intercept`, the design matrix is augmented
/// with the constant column `intercept_scaling` at feature index `n_features`
/// (it IS part of `QD` and `w`, `linear.cpp:512`/`608-618` over the augmented
/// row).
///
/// State / sweep / shrinking (`linear.cpp:576-754`) is transcribed faithfully,
/// including per-sample shrinking (`active_size_i`, `alpha_index`) and the
/// two-level stopping (`eps_shrink = max(10·eps, 1)`). liblinear shuffles the
/// `index` set each sweep (`bounded_rand_int`, `linear.cpp:629`); ferrolearn
/// sweeps NATURAL ORDER (no RNG) — the Crammer-Singer optimum is unique, so the
/// limit is identical (the documented RNG/shrink-path boundary, as the dual/l1
/// solvers). `eps = tol` (the solver receives `param->eps` directly,
/// `linear.cpp:2535`).
///
/// Returns `(w_flat, n_iter, converged)` where `w_flat[feature*nr_class + m]`.
#[allow(
clippy::too_many_lines,
clippy::too_many_arguments,
reason = "faithful transcription of liblinear Solver_MCSVM_CS (linear.cpp:493-787); \
the args mirror the solver's ctor + Solve() inputs (prob/nr_class/weighted_C/\
eps/max_iter/bias) and grouping them would obscure the 1:1 transcription"
)]
fn solve_crammer_singer<F: Float + 'static>(
x: &Array2<F>,
y_class: &[usize],
nr_class: usize,
weighted_c: &[F],
max_iter: usize,
tol: F,
fit_intercept: bool,
intercept_scaling: F,
) -> (Vec<F>, usize, bool) {
let (l, n_features) = x.dim();
let w_size = if fit_intercept {
n_features + 1
} else {
n_features
};
let inf = F::infinity();
let ten = F::from(10).unwrap_or_else(|| {
let two = F::one() + F::one();
two + two + two + two + two
});
let tiny = F::from(1.0e-12).unwrap_or_else(F::epsilon);
// `C[i] = W[i] · weighted_C[y_i]`; `W[i] = 1` (`linear.cpp:521-522`).
let c_per_sample: Vec<F> = y_class.iter().map(|&yi| weighted_c[yi]).collect();
// `x_val(i, feat)` over the augmented row (feature index `n_features` is the
// constant intercept_scaling column when fit_intercept).
let x_val = |i: usize, feat: usize| -> F {
if feat < n_features {
x[[i, feat]]
} else {
intercept_scaling
}
};
// alpha[i*nr_class + m], w[feature*nr_class + m] (`linear.cpp:580-602`).
let mut alpha = vec![F::zero(); l * nr_class];
let mut w = vec![F::zero(); w_size * nr_class];
// alpha_index[i*nr_class + m] = m; QD[i] = ||x_i||^2 over the augmented row;
// active_size_i[i] = nr_class; y_index[i] = y_class[i] (`linear.cpp:603-622`).
let mut alpha_index = vec![0usize; l * nr_class];
let mut qd = vec![F::zero(); l];
let mut active_size_i = vec![nr_class; l];
let mut y_index = y_class.to_vec();
let mut index: Vec<usize> = (0..l).collect();
for i in 0..l {
for m in 0..nr_class {
alpha_index[i * nr_class + m] = m;
}
let mut acc = F::zero();
for feat in 0..w_size {
let v = x_val(i, feat);
acc = acc + v * v;
}
qd[i] = acc;
}
// Scratch buffers reused per sample (`B`, `G`, `alpha_new`,
// `linear.cpp:518-519,581`).
let mut b_buf = vec![F::zero(); nr_class];
let mut g_buf = vec![F::zero(); nr_class];
let mut alpha_new = vec![F::zero(); nr_class];
let mut active_size = l;
// eps_shrink = max(10·eps, 1) (`linear.cpp:590`).
let mut eps_shrink = (ten * tol).max(F::one());
let mut start_from_all = true;
let mut iter = 0usize;
let mut converged = false;
while iter < max_iter {
let mut stopping = -inf;
// liblinear shuffles index[0..active_size] here (`linear.cpp:627-631`);
// ferrolearn sweeps natural order (no RNG); the unique optimum is
// unchanged at the limit.
let mut s = 0;
while s < active_size {
let i = index[s];
let ai = qd[i];
if ai > F::zero() {
let asi = active_size_i[i];
// G[m] = (m==y_i ? 0 : 1) + w_{alpha_index[m]} · x_i
// (`linear.cpp:641-653`).
for g in g_buf.iter_mut().take(asi) {
*g = F::one();
}
if y_index[i] < asi {
g_buf[y_index[i]] = F::zero();
}
for feat in 0..w_size {
let xv = x_val(i, feat);
if xv == F::zero() {
continue;
}
let base = feat * nr_class;
for m in 0..asi {
let idx = alpha_index[i * nr_class + m];
g_buf[m] = g_buf[m] + w[base + idx] * xv;
}
}
// minG over {alpha_i[idx]<0} ∪ {y_i if alpha_i[y_i]<C[i]};
// maxG over all active m (`linear.cpp:655-666`).
let mut min_g = inf;
let mut max_g = -inf;
for m in 0..asi {
let idx = alpha_index[i * nr_class + m];
if alpha[i * nr_class + idx] < F::zero() && g_buf[m] < min_g {
min_g = g_buf[m];
}
if g_buf[m] > max_g {
max_g = g_buf[m];
}
}
if y_index[i] < asi
&& alpha[i * nr_class + y_class[i]] < c_per_sample[i]
&& g_buf[y_index[i]] < min_g
{
min_g = g_buf[y_index[i]];
}
// Per-sample shrinking via be_shrunk (`linear.cpp:668-697`).
let mut m = 0;
while m < active_size_i[i] {
let idx_m = alpha_index[i * nr_class + m];
if cs_be_shrunk(
c_per_sample[i],
m,
y_index[i],
alpha[i * nr_class + idx_m],
g_buf[m],
min_g,
) {
active_size_i[i] -= 1;
while active_size_i[i] > m {
let asi_top = active_size_i[i];
let idx_top = alpha_index[i * nr_class + asi_top];
if !cs_be_shrunk(
c_per_sample[i],
asi_top,
y_index[i],
alpha[i * nr_class + idx_top],
g_buf[asi_top],
min_g,
) {
alpha_index.swap(i * nr_class + m, i * nr_class + asi_top);
g_buf.swap(m, asi_top);
if y_index[i] == asi_top {
y_index[i] = m;
} else if y_index[i] == m {
y_index[i] = asi_top;
}
break;
}
active_size_i[i] -= 1;
}
}
m += 1;
}
if active_size_i[i] <= 1 {
active_size -= 1;
index.swap(s, active_size);
continue; // re-process the swapped-in element at `s`
}
if max_g - min_g <= tiny {
s += 1;
continue;
}
stopping = stopping.max(max_g - min_g);
// B[m] = G[m] - Ai·alpha_i[idx] (`linear.cpp:704-705`).
let asi = active_size_i[i];
for m in 0..asi {
let idx = alpha_index[i * nr_class + m];
b_buf[m] = g_buf[m] - ai * alpha[i * nr_class + idx];
}
cs_solve_sub_problem(ai, y_index[i], c_per_sample[i], asi, &b_buf, &mut alpha_new);
// Apply d = alpha_new[m] - alpha_i[idx] and update w
// (`linear.cpp:708-728`).
let mut d_nz: Vec<(usize, F)> = Vec::new();
for m in 0..asi {
let idx = alpha_index[i * nr_class + m];
let d = alpha_new[m] - alpha[i * nr_class + idx];
alpha[i * nr_class + idx] = alpha_new[m];
if d.abs() >= tiny {
d_nz.push((idx, d));
}
}
if !d_nz.is_empty() {
for feat in 0..w_size {
let xv = x_val(i, feat);
if xv == F::zero() {
continue;
}
let base = feat * nr_class;
for &(idx, d) in &d_nz {
w[base + idx] = w[base + idx] + d * xv;
}
}
}
}
s += 1;
}
iter += 1;
// Two-level stopping (`linear.cpp:738-753`).
if stopping < eps_shrink {
if stopping < tol && start_from_all {
converged = true;
break;
}
active_size = l;
for asi in active_size_i.iter_mut() {
*asi = nr_class;
}
eps_shrink = (eps_shrink / (F::one() + F::one())).max(tol);
start_from_all = true;
} else {
start_from_all = false;
}
}
(w, iter, converged)
}
/// `Solver_MCSVM_CS::be_shrunk` (`linear.cpp:566-574`): shrink class `m` of
/// sample `i` when its dual variable is at its bound and the gradient is below
/// `minG`. `bound = C[i]` if `m == yi` (the class index), else `0`.
fn cs_be_shrunk<F: Float>(c_yi: F, m: usize, yi: usize, alpha_i: F, g_m: F, min_g: F) -> bool {
let bound = if m == yi { c_yi } else { F::zero() };
alpha_i == bound && g_m < min_g
}
/// `Solver_MCSVM_CS::solve_sub_problem` (`linear.cpp:541-564`): the per-sample
/// simplex-projection inner solve. `B[..active_i]` is the gradient offset
/// buffer; the result is written to `alpha_new[..active_i]`.
///
/// Clones `B → D`, adds `Ai·C_yi` to `D[yi]` (when `yi < active_i`), sorts `D`
/// DESCENDING, then finds the breakpoint `beta` and projects each coordinate:
/// `alpha_new[r] = min(C_yi, (beta-B[r])/Ai)` for `r == yi`, else
/// `min(0, (beta-B[r])/Ai)`.
fn cs_solve_sub_problem<F: Float + 'static>(
ai: F,
yi: usize,
c_yi: F,
active_i: usize,
b_buf: &[F],
alpha_new: &mut [F],
) {
// clone(D, B, active_i); D[yi] += Ai·C_yi (`linear.cpp:546-548`).
let mut d: Vec<F> = b_buf[..active_i].to_vec();
if yi < active_i {
d[yi] = d[yi] + ai * c_yi;
}
// qsort DESCENDING (compare_double returns -1 when a > b, `linear.cpp:532-538`).
d.sort_by(|a, b| match b.partial_cmp(a) {
Some(ord) => ord,
None => core::cmp::Ordering::Equal,
});
// beta = D[0] - Ai·C_yi; for r=1; r<active_i && beta<r·D[r]; r++ { beta+=D[r]; }
// beta /= r (`linear.cpp:551-554`).
let mut beta = d[0] - ai * c_yi;
let mut r = 1usize;
while r < active_i {
let r_f = F::from(r).unwrap_or_else(F::one);
if beta < r_f * d[r] {
beta = beta + d[r];
r += 1;
} else {
break;
}
}
let r_f = F::from(r).unwrap_or_else(F::one);
beta = beta / r_f;
// Project each coordinate (`linear.cpp:556-562`).
for (rr, an) in alpha_new.iter_mut().enumerate().take(active_i) {
let cand = (beta - b_buf[rr]) / ai;
*an = if rr == yi {
c_yi.min(cand)
} else {
F::zero().min(cand)
};
}
}
impl<F: Float + Send + Sync + ScalarOperand + 'static> Fit<Array2<F>, Array1<usize>>
for LinearSVC<F>
{
type Fitted = FittedLinearSVC<F>;
type Error = FerroError;
/// Fit the linear SVC model using liblinear's dual coordinate descent.
///
/// Minimizes `0.5*||w||^2 + C * sum_i L(y_i, w.x_i)` (no `1/n`). When
/// `fit_intercept` is set, the design matrix is augmented with a constant
/// column equal to `intercept_scaling`, the augmented weight is penalized
/// like any feature, and `intercept_ = intercept_scaling * w_last`
/// (`_base.py:1188-1198,:1240-1245`).
///
/// # Errors
///
/// - [`FerroError::ShapeMismatch`] — sample count mismatch.
/// - [`FerroError::InvalidParameter`] — `C` not positive, or
/// `intercept_scaling` not positive when fitting an intercept.
/// - [`FerroError::InsufficientSamples`] — fewer than 2 distinct classes.
fn fit(&self, x: &Array2<F>, y: &Array1<usize>) -> Result<FittedLinearSVC<F>, FerroError> {
let (n_samples, n_features) = x.dim();
if n_samples != y.len() {
return Err(FerroError::ShapeMismatch {
expected: vec![n_samples],
actual: vec![y.len()],
context: "y length must match number of samples in X".into(),
});
}
if self.c <= F::zero() {
return Err(FerroError::InvalidParameter {
name: "C".into(),
reason: "must be positive".into(),
});
}
// `_parameter_constraints["tol"] = Interval(Real, 0.0, None,
// closed="neither")` (`_classes.py:237`) → `tol <= 0` raises.
if self.tol <= F::zero() {
return Err(FerroError::InvalidParameter {
name: "tol".into(),
reason: "must be positive".into(),
});
}
// liblinear raises when intercept_scaling <= 0 with fit_intercept
// (`_base.py:1190-1196`).
if self.fit_intercept && self.intercept_scaling <= F::zero() {
return Err(FerroError::InvalidParameter {
name: "intercept_scaling".into(),
reason: "must be greater than 0 when fit_intercept is true".into(),
});
}
// Non-finite input validation (#2263). sklearn `LinearSVC.fit`
// -> `self._validate_data(X, y, ...)` (`svm/_classes.py:302`, the shared
// liblinear `_fit_liblinear` path; cf. `svm/_base.py:190`) keeps the
// default `force_all_finite=True`, so `check_array` rejects any NaN or
// +/-inf in X with a `ValueError("Input X contains NaN.")` / `"...
// contains infinity ..."` BEFORE the dual coordinate-descent solve. `y`
// is `Array1<usize>` (integer class labels), finite by type, and
// ferrolearn's `LinearSVC::fit` takes no `sample_weight` argument, so X is
// the only runtime check. `.iter().any(|v| !v.is_finite())` rejects both
// NaN and Inf (bounds-safe, no panic, R-CODE-2); the finite path is
// byte-identical.
if x.iter().any(|v| !v.is_finite()) {
return Err(FerroError::InvalidParameter {
name: "X".into(),
reason: "Input X contains NaN or infinity.".into(),
});
}
let y_vec: Vec<usize> = y.to_vec();
let mut classes: Vec<usize> = y_vec.clone();
classes.sort_unstable();
classes.dedup();
if classes.len() < 2 {
return Err(FerroError::InsufficientSamples {
required: 2,
actual: classes.len(),
context: "LinearSVC requires at least 2 distinct classes".into(),
});
}
// `class_weight` scales `C` per class: `weighted_C[i] = C·w[class i]`
// (`compute_class_weight(class_weight, classes=classes_, y=y)`,
// `_base.py:1179`; `linear.cpp:2496-2507`). `weights` is aligned to
// `classes` (sorted unique = LabelEncoder order). Used by BOTH the OvR
// path (below) and the Crammer-Singer joint solver.
let weights = compute_class_weight(&self.class_weight, &classes, &y_vec);
// `multi_class='crammer_singer'` selects liblinear solver type 4
// REGARDLESS of penalty/loss/dual (`_get_liblinear_solver_type` returns
// 4 for crammer_singer, `_base.py:1017,1020-1021`), so penalty/loss/dual
// are IGNORED and the joint Crammer-Singer solver runs.
if self.multi_class == MultiClass::CrammerSinger {
return self.fit_crammer_singer(x, &classes, &y_vec, &weights);
}
// Resolve `dual` (auto → bool via `_validate_dual_parameter`,
// `_classes.py:13-29`) then validate the penalty×loss×dual combination
// against the liblinear solver matrix (`_get_liblinear_solver_type`,
// `_base.py:995-1049`). Unsupported combinations raise (ValueError →
// `InvalidParameter`). The resolved solver type selects the per-
// sub-problem solver below.
let dual = resolve_dual(self.dual, self.penalty, self.loss, n_samples, n_features);
let _solver_type = liblinear_solver_type(self.penalty, self.loss, dual)?;
// Solve one binary sub-problem (positive group penalty `cp`, negative
// group penalty `cn`) and split the augmented weight vector into
// (coef_, intercept_) per `_base.py:1240-1245`. For `penalty=l1`
// use liblinear's L1 coordinate descent (`solve_l1r_l2_svc`, solver
// type 5, `linear.cpp:1467`); for `penalty=l2` keep the dual CD — the
// l2 optimum is dual-invariant (R-DEV-7), so dual=True and dual=False
// reach the same `coef_`/`intercept_`.
let penalty = self.penalty;
let solve_one = |y_signed: &Array1<F>, cp: F, cn: F| -> (Array1<F>, F, usize, bool) {
let cfg = SolverConfig {
cp,
cn,
max_iter: self.max_iter,
tol: self.tol,
loss: self.loss,
fit_intercept: self.fit_intercept,
intercept_scaling: self.intercept_scaling,
};
let (w, n_iter, converged) = match penalty {
LinearSVCPenalty::L1 => solve_binary_l1r_l2(x, y_signed, &cfg),
LinearSVCPenalty::L2 => solve_binary_dual(x, y_signed, &cfg),
};
let coef = Array1::from_iter(w.iter().take(n_features).copied());
let intercept = if self.fit_intercept {
self.intercept_scaling * w[n_features]
} else {
F::zero()
};
(coef, intercept, n_iter, converged)
};
let mut any_unconverged = false;
// n_iter_ = n_iter_.max() across the sub-problem fits (`_classes.py:338`).
let mut max_n_iter: usize = 0;
let fitted = if classes.len() == 2 {
// Binary classification. liblinear's `prob.y` are the LabelEncoder
// indices, mapped `y[i] = +1 if y_ind > 0 else -1`
// (`linear.cpp:861-871`), so the positive class is `classes_[1]`
// and `coef_` is for `classes_[1]` (`_base.py` / `_classes.py:311`).
let y_signed: Array1<F> = y.mapv(|label| {
if label == classes[1] {
F::one()
} else {
-F::one()
}
});
// Binary (`train_one(Cp=weighted_C[1], Cn=weighted_C[0])`,
// `linear.cpp:2543-2551`): Cp = C·w[classes[1]] (the +1 class),
// Cn = C·w[classes[0]] (the −1 class).
let cp = self.c * weights[1];
let cn = self.c * weights[0];
let (coef, intercept, n_iter, converged) = solve_one(&y_signed, cp, cn);
if !converged {
any_unconverged = true;
}
max_n_iter = max_n_iter.max(n_iter);
FittedLinearSVC {
weight_vectors: vec![coef],
intercepts: vec![intercept],
classes,
is_binary: true,
n_features,
n_iter: max_n_iter,
}
} else {
// Multiclass: one-vs-rest. Each class is the positive (+1) class of
// its own binary sub-problem.
let mut weight_vectors = Vec::with_capacity(classes.len());
let mut intercepts = Vec::with_capacity(classes.len());
for (k, &cls) in classes.iter().enumerate() {
let y_signed: Array1<F> =
y.mapv(|label| if label == cls { F::one() } else { -F::one() });
// Multiclass OvR (`train_one(Cp=weighted_C[k], Cn=param->C)`,
// `linear.cpp:2559-2571`): Cp = C·w[class k] (the +1 class),
// Cn = C (the BASE C — the negative rest is UNWEIGHTED).
let cp = self.c * weights[k];
let cn = self.c;
let (coef, intercept, n_iter, converged) = solve_one(&y_signed, cp, cn);
if !converged {
any_unconverged = true;
}
max_n_iter = max_n_iter.max(n_iter);
weight_vectors.push(coef);
intercepts.push(intercept);
}
FittedLinearSVC {
weight_vectors,
intercepts,
classes,
is_binary: false,
n_features,
n_iter: max_n_iter,
}
};
// liblinear warns when any sub-problem reaches `max_iter` without
// satisfying the stopping criterion (`_base.py:1234-1238`). The crate's
// warning channel is `eprintln!` (cf. qda.rs / lda.rs warnings).
if any_unconverged {
eprintln!("Liblinear failed to converge, increase the number of iterations.");
}
Ok(fitted)
}
}
impl<F: Float + Send + Sync + ScalarOperand + 'static> LinearSVC<F> {
/// Fit via the joint Crammer-Singer multiclass solver (solver type 4,
/// `_base.py:1017`). Runs ONE joint solve over all classes
/// ([`solve_crammer_singer`], `Solver_MCSVM_CS`, `linear.cpp:493-787`),
/// extracts per-class `coef_`/`intercept_` from the flattened
/// `w[feature*nr_class + m]`, and applies the binary collapse
/// (`_classes.py:340-344`).
///
/// `classes` is the sorted unique label set; `y_vec` the per-sample labels;
/// `weights` the per-class `class_weight` multipliers aligned to `classes`.
fn fit_crammer_singer(
&self,
x: &Array2<F>,
classes: &[usize],
y_vec: &[usize],
weights: &[F],
) -> Result<FittedLinearSVC<F>, FerroError> {
let (_n_samples, n_features) = x.dim();
let nr_class = classes.len();
// Map each sample's label to its class index (`0..nr_class`, the sorted
// `classes` position — liblinear's LabelEncoder order). `y_class[i]`
// indexes `w[..][m]` and `weighted_C[m]`.
let y_class: Vec<usize> = y_vec
.iter()
.map(|label| classes.iter().position(|c| c == label).unwrap_or(0))
.collect();
// weighted_C[c] = C · class_weight[c] (`linear.cpp:2496-2507`,
// `weighted_C[i] = param->C` then `*= weight`). `weights` already holds
// the per-class multipliers (1.0 when class_weight=None).
let weighted_c: Vec<F> = weights.iter().map(|&wc| self.c * wc).collect();
let (w, n_iter, converged) = solve_crammer_singer(
x,
&y_class,
nr_class,
&weighted_c,
self.max_iter,
self.tol,
self.fit_intercept,
self.intercept_scaling,
);
if !converged {
eprintln!("Liblinear failed to converge, increase the number of iterations.");
}
// Extract coef_[m][feature] = w[feature*nr_class + m] and
// intercept_[m] = intercept_scaling · w[n_features*nr_class + m]
// (the augmented-column weight), per `_base.py:1240-1245`.
let mut weight_vectors: Vec<Array1<F>> = Vec::with_capacity(nr_class);
let mut intercepts: Vec<F> = Vec::with_capacity(nr_class);
for m in 0..nr_class {
let mut row = Array1::<F>::zeros(n_features);
for (feat, r) in row.iter_mut().enumerate() {
*r = w[feat * nr_class + m];
}
weight_vectors.push(row);
let intercept = if self.fit_intercept {
self.intercept_scaling * w[n_features * nr_class + m]
} else {
F::zero()
};
intercepts.push(intercept);
}
// Binary collapse (`_classes.py:340-344`): with 2 classes, the joint
// solve yields two rows; sklearn collapses them to a single weight
// vector `coef_ = row_1 - row_0` and `intercept_ = int_1 - int_0`, so
// the binary sign decision path is used.
if nr_class == 2 {
let collapsed_coef = &weight_vectors[1] - &weight_vectors[0];
let collapsed_intercept = if self.fit_intercept {
intercepts[1] - intercepts[0]
} else {
F::zero()
};
return Ok(FittedLinearSVC {
weight_vectors: vec![collapsed_coef],
intercepts: vec![collapsed_intercept],
classes: classes.to_vec(),
is_binary: true,
n_features,
n_iter,
});
}
Ok(FittedLinearSVC {
weight_vectors,
intercepts,
classes: classes.to_vec(),
is_binary: false,
n_features,
n_iter,
})
}
}
impl<F: Float + Send + Sync + ScalarOperand + 'static> Predict<Array2<F>> for FittedLinearSVC<F> {
type Output = Array1<usize>;
type Error = FerroError;
/// Predict class labels for the given feature matrix.
///
/// Binary: `sign(X @ w + b)` mapped to class labels (`>= 0 → classes_[1]`).
/// Multiclass: argmax of decision values across one-vs-rest classifiers.
///
/// # Errors
///
/// Returns [`FerroError::ShapeMismatch`] if the number of features
/// does not match the fitted model.
fn predict(&self, x: &Array2<F>) -> Result<Array1<usize>, FerroError> {
let n_features = x.ncols();
if n_features != self.n_features {
return Err(FerroError::ShapeMismatch {
expected: vec![self.n_features],
actual: vec![n_features],
context: "number of features must match fitted model".into(),
});
}
let n_samples = x.nrows();
let mut predictions = Array1::<usize>::zeros(n_samples);
if self.is_binary {
let scores = x.dot(&self.weight_vectors[0]) + self.intercepts[0];
for i in 0..n_samples {
predictions[i] = if scores[i] >= F::zero() {
self.classes[1]
} else {
self.classes[0]
};
}
} else {
// Multiclass: pick class with highest decision value.
for i in 0..n_samples {
let mut best_class = 0;
let mut best_score = F::neg_infinity();
for (c, w) in self.weight_vectors.iter().enumerate() {
let score = x.row(i).dot(w) + self.intercepts[c];
if score > best_score {
best_score = score;
best_class = c;
}
}
predictions[i] = self.classes[best_class];
}
}
Ok(predictions)
}
}
impl<F: Float + Send + Sync + ScalarOperand + 'static> HasCoefficients<F> for FittedLinearSVC<F> {
/// Returns the coefficient vector of the first (or only) binary sub-problem.
fn coefficients(&self) -> &Array1<F> {
&self.weight_vectors[0]
}
/// Returns the intercept of the first (or only) binary sub-problem.
fn intercept(&self) -> F {
self.intercepts[0]
}
}
impl<F: Float + Send + Sync + ScalarOperand + 'static> HasClasses for FittedLinearSVC<F> {
fn classes(&self) -> &[usize] {
&self.classes
}
fn n_classes(&self) -> usize {
self.classes.len()
}
}
#[cfg(test)]
mod tests {
use super::*;
use ndarray::array;
#[test]
fn test_default_constructor() {
let m = LinearSVC::<f64>::new();
assert_eq!(m.max_iter, 1000);
assert!(m.c == 1.0);
assert_eq!(m.loss, LinearSVCLoss::SquaredHinge);
assert!(m.fit_intercept);
assert!(m.intercept_scaling == 1.0);
}
#[test]
fn test_builder_setters() {
let m = LinearSVC::<f64>::new()
.with_c(10.0)
.with_max_iter(500)
.with_tol(1e-6)
.with_loss(LinearSVCLoss::Hinge)
.with_fit_intercept(false)
.with_intercept_scaling(2.0);
assert!(m.c == 10.0);
assert_eq!(m.max_iter, 500);
assert_eq!(m.loss, LinearSVCLoss::Hinge);
assert!(!m.fit_intercept);
assert!(m.intercept_scaling == 2.0);
}
#[test]
fn test_binary_classification() {
let x = Array2::from_shape_vec(
(8, 2),
vec![
1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 2.0, 2.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0, 9.0, 9.0,
],
)
.unwrap();
let y = array![0, 0, 0, 0, 1, 1, 1, 1];
let model = LinearSVC::<f64>::new().with_c(1.0).with_max_iter(5000);
let fitted = model.fit(&x, &y).unwrap();
let preds = fitted.predict(&x).unwrap();
let correct: usize = preds.iter().zip(y.iter()).filter(|(p, a)| p == a).count();
assert!(correct >= 6, "expected at least 6 correct, got {correct}");
}
#[test]
fn test_binary_hinge_loss() {
let x = Array2::from_shape_vec(
(6, 2),
vec![1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0],
)
.unwrap();
let y = array![0, 0, 0, 1, 1, 1];
let model = LinearSVC::<f64>::new()
.with_loss(LinearSVCLoss::Hinge)
.with_max_iter(5000);
let fitted = model.fit(&x, &y).unwrap();
let preds = fitted.predict(&x).unwrap();
let correct: usize = preds.iter().zip(y.iter()).filter(|(p, a)| p == a).count();
assert!(correct >= 4, "expected at least 4 correct, got {correct}");
}
#[test]
fn test_multiclass_classification() {
let x = Array2::from_shape_vec(
(9, 2),
vec![
0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 10.0, 0.0, 10.5, 0.0, 10.0, 0.5, 0.0, 10.0, 0.5,
10.0, 0.0, 10.5,
],
)
.unwrap();
let y = array![0, 0, 0, 1, 1, 1, 2, 2, 2];
let model = LinearSVC::<f64>::new().with_c(10.0).with_max_iter(5000);
let fitted = model.fit(&x, &y).unwrap();
assert_eq!(fitted.n_classes(), 3);
assert_eq!(fitted.classes(), &[0, 1, 2]);
let preds = fitted.predict(&x).unwrap();
let correct: usize = preds.iter().zip(y.iter()).filter(|(p, a)| p == a).count();
assert!(correct >= 7, "expected at least 7 correct, got {correct}");
}
#[test]
fn test_binary_decision_function_is_1d() {
let x = Array2::from_shape_vec(
(6, 2),
vec![1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0],
)
.unwrap();
let y = array![0, 0, 0, 1, 1, 1];
let fitted = LinearSVC::<f64>::new()
.with_max_iter(5000)
.fit(&x, &y)
.unwrap();
let df = fitted.decision_function(&x).unwrap();
// sklearn ravels the binary decision_function to (n,) (`_base.py:365`).
let binary = df.as_binary().expect("binary decision is 1-D");
assert_eq!(binary.len(), 6);
assert!(df.as_multiclass().is_none());
assert_eq!(df.n_samples(), 6);
}
#[test]
fn test_multiclass_decision_function_is_2d() {
let x = Array2::from_shape_vec(
(9, 2),
vec![
0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 10.0, 0.0, 10.5, 0.0, 10.0, 0.5, 0.0, 10.0, 0.5,
10.0, 0.0, 10.5,
],
)
.unwrap();
let y = array![0, 0, 0, 1, 1, 1, 2, 2, 2];
let fitted = LinearSVC::<f64>::new()
.with_c(10.0)
.with_max_iter(5000)
.fit(&x, &y)
.unwrap();
let df = fitted.decision_function(&x).unwrap();
let scores = df.as_multiclass().expect("multiclass decision is 2-D");
assert_eq!(scores.dim(), (9, 3));
assert!(df.as_binary().is_none());
}
#[test]
fn test_shape_mismatch_fit() {
let x = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
let y = array![0, 1]; // Wrong length
let model = LinearSVC::<f64>::new();
assert!(model.fit(&x, &y).is_err());
}
#[test]
fn test_invalid_c() {
let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
let y = array![0, 0, 1, 1];
let model = LinearSVC::<f64>::new().with_c(0.0);
assert!(model.fit(&x, &y).is_err());
let model_neg = LinearSVC::<f64>::new().with_c(-1.0);
assert!(model_neg.fit(&x, &y).is_err());
}
#[test]
fn test_invalid_intercept_scaling() {
let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
let y = array![0, 0, 1, 1];
// fit_intercept=true + intercept_scaling<=0 is rejected
// (`_base.py:1190-1196`).
let model = LinearSVC::<f64>::new().with_intercept_scaling(0.0);
assert!(model.fit(&x, &y).is_err());
// But with fit_intercept=false, intercept_scaling is ignored.
let model = LinearSVC::<f64>::new()
.with_fit_intercept(false)
.with_intercept_scaling(0.0);
assert!(model.fit(&x, &y).is_ok());
}
#[test]
fn test_fit_intercept_false_zero_intercept() {
let x = Array2::from_shape_vec(
(6, 2),
vec![1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0],
)
.unwrap();
let y = array![0, 0, 0, 1, 1, 1];
let fitted = LinearSVC::<f64>::new()
.with_fit_intercept(false)
.with_max_iter(5000)
.fit(&x, &y)
.unwrap();
assert!(fitted.intercept() == 0.0);
}
#[test]
fn test_single_class_error() {
let x = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
let y = array![0, 0, 0];
let model = LinearSVC::<f64>::new();
assert!(model.fit(&x, &y).is_err());
}
#[test]
fn test_has_coefficients() {
let x = Array2::from_shape_vec(
(6, 2),
vec![1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0],
)
.unwrap();
let y = array![0, 0, 0, 1, 1, 1];
let model = LinearSVC::<f64>::new().with_max_iter(5000);
let fitted = model.fit(&x, &y).unwrap();
assert_eq!(fitted.coefficients().len(), 2);
}
#[test]
#[allow(
clippy::assertions_on_constants,
reason = "assert!(false) reports the unexpected-Err fit path without a gated panic!/expect"
)]
fn test_l1_penalty_smoke() {
// Smoke check that the `penalty=l1` path (solve_binary_l1r_l2,
// `linear.cpp:1467`, solver type 5, `_base.py:1014`) lands near the live
// sklearn 1.5.2 oracle. The rigorous oracle pin is the critic's to add.
//
// Oracle (live sklearn 1.5.2; values per R-CHAR-3 — NEVER copied from
// ferrolearn):
// python3 -c "import numpy as np; from sklearn.svm import LinearSVC; \
// X=np.array([[1.,1.],[1.,2.],[2.,1.],[2.,2.],[8.,8.],[8.,9.],[9.,8.],[9.,9.]]); \
// y=np.array([0,0,0,0,1,1,1,1]); \
// m=LinearSVC(penalty='l1',loss='squared_hinge',dual=False,C=1.0, \
// fit_intercept=True,max_iter=200000,tol=1e-10).fit(X,y); \
// print(m.coef_.tolist(), m.intercept_.tolist())"
// # coef [[0.1283185834966579, 0.12831858464059265]] int [-1.2079646017762715]
const SK_COEF_0: f64 = 0.1283185834966579;
const SK_COEF_1: f64 = 0.12831858464059265;
const SK_INTERCEPT: f64 = -1.2079646017762715;
let x = array![
[1.0, 1.0],
[1.0, 2.0],
[2.0, 1.0],
[2.0, 2.0],
[8.0, 8.0],
[8.0, 9.0],
[9.0, 8.0],
[9.0, 9.0],
];
let y = array![0usize, 0, 0, 0, 1, 1, 1, 1];
let result = LinearSVC::<f64>::new()
.with_penalty(LinearSVCPenalty::L1)
.with_loss(LinearSVCLoss::SquaredHinge)
.with_dual(DualMode::False)
.with_c(1.0)
.with_max_iter(200_000)
.with_tol(1e-10)
.fit(&x, &y);
let Ok(fitted) = result else {
assert!(false, "l1 fit must succeed");
return;
};
let coef = fitted.coefficients();
assert!(
(coef[0] - SK_COEF_0).abs() < 1e-2,
"l1 coef[0] {} vs oracle {SK_COEF_0}",
coef[0]
);
assert!(
(coef[1] - SK_COEF_1).abs() < 1e-2,
"l1 coef[1] {} vs oracle {SK_COEF_1}",
coef[1]
);
assert!(
(fitted.intercept() - SK_INTERCEPT).abs() < 1e-2,
"l1 intercept {} vs oracle {SK_INTERCEPT}",
fitted.intercept()
);
}
#[test]
#[allow(
clippy::assertions_on_constants,
reason = "assert!(false) reports the unexpected-Err fit path without a gated panic!/expect"
)]
fn test_class_weight_smoke() {
// Smoke check that `class_weight` (scaling `C` per class —
// `compute_class_weight`, `_base.py:1179`; `weighted_C[i] = C·w[i]`,
// `linear.cpp:2496-2507`) lands near the live sklearn 1.5.2 oracle. The
// rigorous oracle pin is the critic's to add.
//
// Oracle (live sklearn 1.5.2; values per R-CHAR-3 — NEVER copied from
// ferrolearn):
// python3 -c "import numpy as np; from sklearn.svm import LinearSVC; \
// X=np.array([[1.,1.],[1.,2.],[2.,1.],[2.,2.],[1.5,1.5],[2.,1.5],[8.,8.],[9.,9.]]); \
// y=np.array([0,0,0,0,0,0,1,1]); \
// for cw in (None,'balanced',{0:1,1:5}): \
// m=LinearSVC(C=1.0,loss='squared_hinge',dual=True,fit_intercept=True, \
// max_iter=200000,tol=1e-10,class_weight=cw).fit(X,y); \
// print(cw, m.coef_.tolist(), m.intercept_.tolist())"
// # None coef [[0.10056447415875154, 0.15957404219329038]] int [-1.263461307484969]
// # balanced coef [[0.09936888940946959, 0.16666283617002833]] int [-1.2132032194327564]
// # (compute_class_weight('balanced',...) = [0.6667, 2.0])
// # {0:1,1:5} coef [[0.11058720549912869, 0.17164468739390437]] int [-1.2954689964246575]
let x = array![
[1.0, 1.0],
[1.0, 2.0],
[2.0, 1.0],
[2.0, 2.0],
[1.5, 1.5],
[2.0, 1.5],
[8.0, 8.0],
[9.0, 9.0],
];
let y = array![0usize, 0, 0, 0, 0, 0, 1, 1];
// (class_weight, expected coef_, expected intercept_) from the oracle.
let cases: [(ClassWeight<f64>, [f64; 2], f64); 3] = [
(
ClassWeight::None,
[0.100_564_474_158_751_54, 0.159_574_042_193_290_38],
-1.263_461_307_484_969,
),
(
ClassWeight::Balanced,
[0.099_368_889_409_469_59, 0.166_662_836_170_028_33],
-1.213_203_219_432_756_4,
),
(
ClassWeight::Explicit(vec![(0, 1.0), (1, 5.0)]),
[0.110_587_205_499_128_69, 0.171_644_687_393_904_37],
-1.295_468_996_424_657_5,
),
];
for (cw, exp_coef, exp_int) in cases {
let result = LinearSVC::<f64>::new()
.with_loss(LinearSVCLoss::SquaredHinge)
.with_dual(DualMode::True)
.with_c(1.0)
.with_class_weight(cw)
.with_max_iter(200_000)
.with_tol(1e-10)
.fit(&x, &y);
let Ok(fitted) = result else {
assert!(false, "class_weight fit must succeed");
return;
};
let coef = fitted.coefficients();
assert!(
(coef[0] - exp_coef[0]).abs() < 1e-2,
"coef[0] {} vs oracle {}",
coef[0],
exp_coef[0]
);
assert!(
(coef[1] - exp_coef[1]).abs() < 1e-2,
"coef[1] {} vs oracle {}",
coef[1],
exp_coef[1]
);
assert!(
(fitted.intercept() - exp_int).abs() < 1e-2,
"intercept {} vs oracle {exp_int}",
fitted.intercept()
);
}
}
#[test]
fn test_unsupported_combinations_rejected() {
let x = array![[1.0, 1.0], [2.0, 2.0], [8.0, 8.0], [9.0, 9.0]];
let y = array![0usize, 0, 1, 1];
// l1 + hinge is unsupported for any dual (`_base.py:1013` has no `l1`
// under `hinge`).
assert!(
LinearSVC::<f64>::new()
.with_penalty(LinearSVCPenalty::L1)
.with_loss(LinearSVCLoss::Hinge)
.fit(&x, &y)
.is_err(),
"l1 + hinge must be rejected"
);
// l2 + hinge + dual=false is unsupported (`hinge: {l2: {dual=true: 3}}`,
// no dual=false entry).
assert!(
LinearSVC::<f64>::new()
.with_penalty(LinearSVCPenalty::L2)
.with_loss(LinearSVCLoss::Hinge)
.with_dual(DualMode::False)
.fit(&x, &y)
.is_err(),
"l2 + hinge + dual=false must be rejected"
);
// l1 + squared_hinge + dual=true is unsupported (`squared_hinge: {l1:
// {dual=false: 5}}`, no dual=true entry).
assert!(
LinearSVC::<f64>::new()
.with_penalty(LinearSVCPenalty::L1)
.with_loss(LinearSVCLoss::SquaredHinge)
.with_dual(DualMode::True)
.fit(&x, &y)
.is_err(),
"l1 + squared_hinge + dual=true must be rejected"
);
}
#[test]
fn test_dual_auto_resolution() {
// n_samples (4) >= n_features (2): auto prefers dual=false. hinge+l2:
// dual=false has no solver, so auto must fall back to dual=true
// (type 3, `_classes.py:22-27`) — the fit succeeds rather than rejecting.
let x = array![[1.0, 1.0], [2.0, 2.0], [8.0, 8.0], [9.0, 9.0]];
let y = array![0usize, 0, 1, 1];
assert!(
LinearSVC::<f64>::new()
.with_loss(LinearSVCLoss::Hinge)
.with_dual(DualMode::Auto)
.fit(&x, &y)
.is_ok(),
"auto must fall back to dual=true for hinge+l2"
);
}
#[test]
#[allow(
clippy::assertions_on_constants,
reason = "assert!(false) reports the unexpected-Err fit path without a gated panic!/expect"
)]
fn test_crammer_singer_smoke() {
// Smoke check that the joint Crammer-Singer solver (solve_crammer_singer,
// Solver_MCSVM_CS, `linear.cpp:493-787`, solver type 4, `_base.py:1017`)
// lands near the live sklearn 1.5.2 oracle for BOTH multiclass and the
// binary collapse (`_classes.py:340-344`). The rigorous oracle pin is the
// critic's to add.
//
// Oracle (live sklearn 1.5.2; values per R-CHAR-3 — NEVER copied from
// ferrolearn):
// python3 -c "import numpy as np; from sklearn.svm import LinearSVC; \
// X=np.array([[0,0],[0.5,0.2],[0.2,0.5],[1,1],[4,4],[4.5,4.2],[4.2,4.5],[5,5],\
// [0,5],[0.5,5.2],[0.2,4.8],[1,6]],dtype=float); \
// y=np.array([0,0,0,0,1,1,1,1,2,2,2,2]); \
// m=LinearSVC(multi_class='crammer_singer',C=1.0,fit_intercept=True, \
// max_iter=200000,tol=1e-10).fit(X,y); \
// print(m.coef_.tolist(), m.intercept_.tolist(), m.predict(X).tolist())"
// # coef [[-0.06761865903618029,-0.24341085269880958],
// # [0.30047599619312043,0.021705426371714007],
// # [-0.23285733712058276,0.22170542636345167]]
// # int [0.9107847137145293,-0.622059023505724,-0.2887256901997209]
// # pred [0,0,0,0,1,1,1,1,2,2,2,2]
let x = array![
[0.0, 0.0],
[0.5, 0.2],
[0.2, 0.5],
[1.0, 1.0],
[4.0, 4.0],
[4.5, 4.2],
[4.2, 4.5],
[5.0, 5.0],
[0.0, 5.0],
[0.5, 5.2],
[0.2, 4.8],
[1.0, 6.0],
];
let y = array![0usize, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2];
// Per-class oracle coef_/intercept_ (rows aligned to classes_ [0,1,2]).
let exp_coef: [[f64; 2]; 3] = [
[-0.067_618_659_036_180_29, -0.243_410_852_698_809_58],
[0.300_475_996_193_120_43, 0.021_705_426_371_714_007],
[-0.232_857_337_120_582_76, 0.221_705_426_363_451_67],
];
let exp_int: [f64; 3] = [
0.910_784_713_714_529_3,
-0.622_059_023_505_724,
-0.288_725_690_199_720_9,
];
let Ok(fitted) = LinearSVC::<f64>::new()
.with_multi_class(MultiClass::CrammerSinger)
.with_c(1.0)
.with_max_iter(200_000)
.with_tol(1e-10)
.fit(&x, &y)
else {
assert!(false, "crammer_singer multiclass fit must succeed");
return;
};
assert!(!fitted.is_binary, "3-class CS must be multiclass");
let rows = fitted.weight_vectors();
let ints = fitted.intercepts();
assert_eq!(rows.len(), 3);
for m in 0..3 {
assert!(
(rows[m][0] - exp_coef[m][0]).abs() < 1e-2,
"CS coef[{m}][0] {} vs oracle {}",
rows[m][0],
exp_coef[m][0]
);
assert!(
(rows[m][1] - exp_coef[m][1]).abs() < 1e-2,
"CS coef[{m}][1] {} vs oracle {}",
rows[m][1],
exp_coef[m][1]
);
assert!(
(ints[m] - exp_int[m]).abs() < 1e-2,
"CS intercept[{m}] {} vs oracle {}",
ints[m],
exp_int[m]
);
}
// predict is all-correct (argmax over the per-class scores).
let Ok(preds) = fitted.predict(&x) else {
assert!(false, "CS multiclass predict must succeed");
return;
};
assert_eq!(
preds.to_vec(),
y.to_vec(),
"CS multiclass predict must be all-correct"
);
// Binary collapse (`_classes.py:340-344`): with 2 classes the joint solve
// collapses to a SINGLE weight vector (`coef_` shape (1,2)) +
// `intercept_ = int_1 - int_0`.
// python3 -c "import numpy as np; from sklearn.svm import LinearSVC; \
// X=np.array([[1,1],[1,2],[2,1],[2,2],[8,8],[8,9],[9,8],[9,9]],dtype=float); \
// y=np.array([0,0,0,0,1,1,1,1]); \
// m=LinearSVC(multi_class='crammer_singer',C=1.0,fit_intercept=True, \
// max_iter=200000,tol=1e-10).fit(X,y); \
// print(m.coef_.tolist(), m.intercept_.tolist())"
// # coef [[0.15503875968992287,0.15503875968992295]] int [-1.4806201550387597]
const BIN_COEF_0: f64 = 0.155_038_759_689_922_87;
const BIN_COEF_1: f64 = 0.155_038_759_689_922_95;
const BIN_INT: f64 = -1.480_620_155_038_759_7;
let xb = array![
[1.0, 1.0],
[1.0, 2.0],
[2.0, 1.0],
[2.0, 2.0],
[8.0, 8.0],
[8.0, 9.0],
[9.0, 8.0],
[9.0, 9.0],
];
let yb = array![0usize, 0, 0, 0, 1, 1, 1, 1];
let Ok(fitted_b) = LinearSVC::<f64>::new()
.with_multi_class(MultiClass::CrammerSinger)
.with_c(1.0)
.with_max_iter(200_000)
.with_tol(1e-10)
.fit(&xb, &yb)
else {
assert!(false, "crammer_singer binary fit must succeed");
return;
};
assert!(fitted_b.is_binary, "2-class CS must collapse to binary");
let coef = fitted_b.coefficients();
assert_eq!(coef.len(), 2, "binary collapse → single (1,2) weight row");
assert!(
(coef[0] - BIN_COEF_0).abs() < 1e-2,
"CS binary coef[0] {} vs oracle {BIN_COEF_0}",
coef[0]
);
assert!(
(coef[1] - BIN_COEF_1).abs() < 1e-2,
"CS binary coef[1] {} vs oracle {BIN_COEF_1}",
coef[1]
);
assert!(
(fitted_b.intercept() - BIN_INT).abs() < 1e-2,
"CS binary intercept {} vs oracle {BIN_INT}",
fitted_b.intercept()
);
let Ok(preds_b) = fitted_b.predict(&xb) else {
assert!(false, "CS binary predict must succeed");
return;
};
assert_eq!(
preds_b.to_vec(),
yb.to_vec(),
"CS binary predict all-correct"
);
}
#[test]
fn test_predict_feature_mismatch() {
let x = Array2::from_shape_vec(
(6, 2),
vec![1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 8.0, 8.0, 8.0, 9.0, 9.0, 8.0],
)
.unwrap();
let y = array![0, 0, 0, 1, 1, 1];
let fitted = LinearSVC::<f64>::new()
.with_max_iter(5000)
.fit(&x, &y)
.unwrap();
let x_bad = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
assert!(fitted.predict(&x_bad).is_err());
}
}