//! Ridge regression (L2-regularized linear regression).
//!
//! This module provides [`Ridge`], which fits a linear model with L2
//! regularization using the closed-form solution:
//!
//! ```text
//! w = (X^T X + alpha * I)^{-1} X^T y
//! ```
//!
//! The regularization parameter `alpha` controls the strength of the
//! L2 penalty, shrinking coefficients toward zero.
//!
//! ## REQ status (per `.design/linear/ridge.md`, mirrors `sklearn/linear_model/_ridge.py` @ 1.5.2)
//!
//! Mirrors `sklearn.linear_model.Ridge` (`_ridge.py:1016`), default dense path
//! `solver='auto'`→`'cholesky'` with `fit_intercept` via centering (intercept unpenalized).
//! coef_/intercept_ match the live sklearn oracle to 1e-8 across alpha∈{0.1,1,10,100}.
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-1 (L2 cholesky fit, intercept unpenalized) | SHIPPED | `Fit for Ridge` (centering + `linalg::solve_ridge`). Consumer: `RsRidge` in `ferrolearn-python/src/regressors.rs`. |
//! | REQ-2 (predict = X·coef + intercept) | SHIPPED | `Predict for FittedRidge`. |
//! | REQ-3 (fit_intercept incl. false) | SHIPPED | `with_fit_intercept`. |
//! | REQ-4 (HasCoefficients introspection) | SHIPPED | `HasCoefficients for FittedRidge`. |
//! | REQ-5 (alpha≥0 validation; alpha=0 → OLS incl. rank-deficient min-norm) | SHIPPED | negative-alpha → `InvalidParameter`; alpha=0 singular falls back `solve_ridge` → `solve_lstsq` (ferray min-norm), mirroring sklearn cholesky→SVD (`_ridge.py:752-756`). Closed #392; test `divergence_ridge_alpha_zero_rank_deficient_min_norm`. |
//! | REQ-6 (multi-output 2-D Y → 2-D coef_) | SHIPPED | `FittedRidgeMulti<F>` + `Fit<Array2<F>, Array2<F>> for Ridge` (per-target `solve_ridge`, `coefficients()` `(n_features, n_targets)` / `intercepts()` `(n_targets,)`), mirroring sklearn `Ridge` 2-D `coef_` `(n_targets, n_features)` / `intercept_` `(n_targets,)` (`_ridge.py:543`/`:550`). Non-test consumer (R-HONEST-4 — prior NOT-STARTED note predates the multi-output Python binding #29): `RsRidgeMultiOutput in ferrolearn-python/src/regressors.rs` (field `Option<FittedRidgeMulti<f64>>`, transposes `coefficients()` to `(n_targets, n_features)`); `_regressors.py::Ridge.fit` routes the `y.ndim == 2 && y.shape[1] > 1` path to `_RsRidgeMultiOutput`. Verified: `tests/divergence_regressors.py::test_ridge_multioutput_matches_sklearn` (predict/coef ≤ 1e-8), `ferrolearn.Ridge(alpha=a).fit(X, Y_2d)` matches sklearn ≤ 1e-16. |
//! | REQ-7 (per-target alpha array) | SHIPPED | `Ridge<F>` adds `pub alpha_per_target: Option<Array1<F>>` (default `None`) + `with_alpha_per_target(Array1<F>)` builder. On the multi-output `Fit<Array2, Array2>` path, when `Some(alphas)` it validates `alphas.len() == n_targets` (else `ShapeMismatch`) and each `alphas[k] >= 0` (else `InvalidParameter`), then solves each target column `k` independently with its own penalty via `linalg::solve_ridge` on the SAME centered (fit_intercept) / raw design the scalar path uses — mathematically identical to an independent scalar-`alpha` Ridge per column, mirroring sklearn's array-valued `alpha` (`_ridge.py:701-712`). `None` is byte-identical to the historic scalar `solve_ridge_multi` path. Oracle tests: `ridge_per_target_alpha_matches_sklearn` (alpha `[0.5, 2.0]`, coef col0 `[0.79891892, 1.43891892]`/col1 `[0.78, 0.355]`, intercepts `[-0.75351351, -0.065]`), `ridge_per_target_alpha_equals_independent_scalar_fits`, `ridge_per_target_alpha_length_mismatch_errors`, `ridge_multi_scalar_alpha_unchanged` (regression guard). Closes #385. |
//! | REQ-8a (dense solver variants auto/cholesky/svd + solver_) | SHIPPED | `pub enum RidgeSolver { #[default] Auto, Cholesky, Svd }` + `Ridge<F>` gains `pub solver: RidgeSolver` (default `Auto`) + `with_solver`; `FittedRidge<F>` gains `solver_` + `pub fn solver()`. The single-output `fit_with_sample_weight` resolves `Auto`→`Cholesky` (`resolve_solver`, `_ridge.py:830`) and dispatches the unconstrained dense solve via `solve_unconstrained`: `Svd` → `linalg::solve_ridge_svd` (`coef = V·diag(sᵢ/(sᵢ²+alpha))·Uᵀy` from the thin SVD via `ferray::linalg::svd`, the analog of sklearn `_solve_svd` `_ridge.py:200-216`); `Cholesky`/`Auto` → the unchanged `linalg::solve_ridge` (byte-identical). All dense solvers yield the IDENTICAL unique solution (strictly convex). `solver_` stores the resolved value. Governs only the single-output unconstrained dense path. Consumer: `Fit<Array2, Array1>::fit for Ridge` (production: `RsRidge::fit` in `ferrolearn-python/src/regressors.rs`; `ridge_cv.rs`). Oracle tests `ridge_solver_svd_matches_sklearn_and_cholesky` (coef `[0.8228070175, 1.3561403509]`, intercept `-0.5768421053`), `ridge_solver_resolution`, `ridge_solver_default_cholesky_unchanged`, `ridge_solver_svd_no_intercept`. Closes #386 (8a). |
//! | REQ-8b (iterative solver variants lsqr/sparse_cg/sag/saga/lbfgs) | NOT-STARTED | #386 — needs iterative/SGD substrate + RNG. Not represented as `RidgeSolver` variants until the substrate lands. |
//! | REQ-9 (positive=True) | SHIPPED | `Ridge<F>` adds `pub positive: bool` (default `false`, `_ridge.py:902`/`:911`) + `with_positive(bool)` builder. When `self.positive`, `fit_with_sample_weight` routes the coefficient solve through `solve_nonneg_ridge`, which first computes the UNCONSTRAINED Cholesky optimum (`crate::linalg::solve_ridge`); if that solve succeeds AND all components are `≥ 0` the positivity constraint is INACTIVE and (by strict convexity + KKT) that feasible unconstrained minimizer IS the constrained minimizer, so it is returned EXACTLY as `(w_unc, 0)` (the exact closed-form Cholesky optimum, `n_iter_ = Some(0)`). Otherwise (solve errored, or any component `< 0` → constraint active) it falls through to projected coordinate descent minimizing `0.5·‖A·w−b‖² + 0.5·alpha·‖w‖²` s.t. `w ≥ 0` (`new = max(0, (A[:,j]ᵀr + col_sq[j]·old)/(col_sq[j] + alpha))`, incremental residual update, `max_iter=self.max_iter.unwrap_or(1000)`/`self.tol`). Either way it runs on the SAME centered (and `√w`-rescaled) design `solve_ridge` uses, then recovers `intercept = y_off − x_off·coef` (fit_intercept) / `0` identically. R-DEV-6: sklearn's default positive-Ridge L-BFGS-B (`_solve_lbfgs`, `_ridge.py:300`, objective `0.5·‖Xw−y‖²+0.5·alpha·‖w‖²`, bounds `[(0,inf)]`, dispatched at `_ridge.py:923-928`) runs at `tol=1e-4` and UNDER-CONVERGES on interior optima by ~1e-4; at `tol=1e-12` it reaches the true optimum = ferrolearn's exact Cholesky optimum (to ~1e-11). ferrolearn therefore ships the TRUE non-negative-ridge optimum, matching a CONVERGED sklearn exactly and being at least as precise as the default. `positive=false` (default) is byte-identical to the unconstrained Cholesky path. Oracle tests: `ridge_positive_matches_sklearn` (active constraint, alpha=1, fit_intercept coef `[1.19891304, 0.0]`, intercept `-6.17744565`, all ≥ 0, differs from unconstrained `[0.95708502, -1.85401484]`), `ridge_positive_false_unchanged` (byte-identical guard), `ridge_positive_all_nonneg_equals_unconstrained` (inactive-constraint sanity); divergence suite `divergence_ridge_positive_interior` (interior optimum == converged sklearn tol=1e-12 == exact unconstrained Cholesky). Closes #387, #2131. |
//! | REQ-10 (max_iter/tol + n_iter_) | SHIPPED | `Ridge<F>` adds `pub max_iter: Option<usize>` (default `None`) and `pub tol: F` (default `1e-4`) with `with_max_iter`/`with_tol` builders. `FittedRidge<F>` adds `n_iter_: Option<usize>` (always `None` for the direct Cholesky solver) with `pub fn n_iter(&self) -> Option<usize>`. Mirrors sklearn ctor `max_iter=None, tol=1e-4` (`_ridge.py:899-900`) and `n_iter_` set at `_ridge.py:994`; `max_iter`/`tol` are no-ops for the direct solver (closed-form normal equations, no iteration) — matching sklearn's direct `cholesky`/`svd` paths which also yield `n_iter_=None`. Test: `ridge_max_iter_tol_niter_defaults_and_builders`. Closes #388. |
//! | REQ-11 (sample_weight) | SHIPPED | `Ridge::fit_with_sample_weight(x, y, sample_weight: Option<&Array1<F>>)` solves WEIGHTED ridge `min Σᵢ wᵢ(yᵢ−xᵢ·coef)² + alpha·‖coef‖²`: weighted offsets `x_off[j]=Σwᵢx[i,j]/Σwᵢ`, `y_off=Σwᵢyᵢ/Σwᵢ` (fit_intercept), centering, then `√wᵢ` row-rescaling (`_rescale_data`, `_ridge.py:682-688`), `linalg::solve_ridge(&Xs, &ys, alpha)` with the penalty `alpha` UNSCALED (since `Xsᵀ·Xs == Xᵀ·W·X`), `intercept = y_off − x_off·coef`; `fit_intercept=false` skips centering (raw `√w`-rescale, intercept 0). `Fit::fit` delegates `fit_with_sample_weight(x, y, None)` (None byte-identical to the historic centering + `solve_ridge` body; alpha=0 OLS min-norm fallback preserved). Oracle tests `ridge_fit_sample_weight_with_intercept_matches_sklearn` (alpha=1 coef `[0.9233502538, 1.39678511]`, intercept `-0.8033840948`, differs from unweighted `[0.8228070175, 1.3561403509]`), `ridge_fit_sample_weight_no_intercept_matches_sklearn` (alpha=2 coef `[0.7273779983, 1.3737799835]`, intercept 0), `ridge_fit_none_sample_weight_equals_unweighted` (byte-identical guard). Closes #389. |
//! | REQ-12 (copy_X/random_state) | SHIPPED | `Ridge<F>` adds `pub copy_x: bool` (default `true`) and `pub random_state: Option<u64>` (default `None`) fields with `with_copy_x`/`with_random_state` builders. `copy_x` ABI-only (fit never mutates `x`); `random_state` stored-but-no-op for the deterministic Cholesky solver (only `sag`/`saga` use it, `_ridge.py:898`/`:903`). Test: `ridge_copy_x_random_state_defaults_and_builders`. Closes #390. |
//! | REQ-13 (ferray substrate) | NOT-STARTED | #391 (alpha=0 fallback already on ferray::linalg::lstsq; coef return tied to #359). |
//! | REQ-14 (non-finite input rejected) | SHIPPED | Both fit entries reject any NaN/+/-inf in X, y, or `sample_weight` BEFORE centering/solve with `FerroError::InvalidParameter`, mirroring sklearn's `_validate_data(force_all_finite=True)` (`_ridge.py:1242`) + `_check_sample_weight` (default `force_all_finite=True`) → `ValueError("Input X contains NaN.")` / `"... contains infinity ..."`. The single-output `fit_with_sample_weight` (shared entry `Fit::fit` delegates to) checks X/y/sample_weight; the SEPARATE multi-output arm `Fit<Array2, Array2>::fit` checks X/y independently. `.iter().any(|v| !v.is_finite())` catches both NaN and Inf; the finite path is byte-identical (the guard never fires on finite input). Verified vs the live sklearn 1.5.2 oracle (R-CHAR-3): `Ridge().fit` raises `ValueError` for NaN/+inf/-inf in X, NaN/inf in y, and NaN/inf in sample_weight (`tests/divergence_linear_nonfinite_batch2.rs::ridge_*`). Non-test consumer: the existing `Fit::fit` / `RsRidge` consumers. (#2259) |
//!
//! acto-critic: core L2 numerics (coef/intercept, alpha scaling, fit_intercept, f32) match the
//! live oracle; one divergence (#392, alpha=0 rank-deficient min-norm) found and fixed.
//! Two states only per goal.md R-DEFER-2.
//!
//! # Examples
//!
//! ```
//! use ferrolearn_linear::Ridge;
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array1, Array2};
//!
//! let model = Ridge::<f64>::new();
//! let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
//! let y = array![2.0, 4.0, 6.0, 8.0];
//!
//! let fitted = model.fit(&x, &y).unwrap();
//! let preds = fitted.predict(&x).unwrap();
//! ```
use ferray::linalg::LinalgFloat;
use ferrolearn_core::error::FerroError;
use ferrolearn_core::introspection::HasCoefficients;
use ferrolearn_core::pipeline::{FittedPipelineEstimator, PipelineEstimator};
use ferrolearn_core::traits::{Fit, Predict};
use ndarray::{Array1, Array2, Axis, ScalarOperand};
use num_traits::{Float, FromPrimitive};
use crate::linalg;
/// Solver used for the dense single-output Ridge coefficient solve.
///
/// Mirrors a (dense) subset of scikit-learn's `solver` parameter
/// (`sklearn/linear_model/_ridge.py:830` `resolve_solver`). Only the
/// closed-form dense solvers are implemented here; the iterative solvers
/// (`lsqr`/`sparse_cg`/`sag`/`saga`/`lbfgs`) are tracked separately as
/// REQ-8b (blocker #386) and are intentionally NOT represented as variants
/// until their iterative/SGD substrate lands.
///
/// Every dense solver returns the IDENTICAL ridge solution — the problem is
/// strictly convex (`alpha > 0`, or `alpha = 0` with full-rank `X`), so the
/// minimizer is unique. The variant therefore only governs *how* the unique
/// solution is computed, not *what* it is.
///
/// This governs only the single-output (`Fit<Array2, Array1>`) dense
/// unconstrained solve. The multi-output (`Fit<Array2, Array2>`) path and the
/// `positive=true` constrained path are unaffected by this setting.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum RidgeSolver {
/// `'auto'` — resolve automatically. For the dense path this resolves to
/// [`Cholesky`](RidgeSolver::Cholesky) (sklearn `resolve_solver`,
/// `_ridge.py:830`: dense, non-positive data → `'cholesky'`).
#[default]
Auto,
/// `'cholesky'` — solve the normal equations `(XᵀX + alpha·I)·w = Xᵀy`
/// via a Cholesky factorization (`_solve_cholesky`, `_ridge.py:741`).
Cholesky,
/// `'svd'` — solve via the singular value decomposition of `X`
/// (`_solve_svd`, `_ridge.py:200`). Numerically identical to
/// [`Cholesky`](RidgeSolver::Cholesky) on the strictly-convex ridge
/// problem (unique minimizer).
Svd,
}
impl RidgeSolver {
/// Resolve `Auto` to the concrete dense solver scikit-learn picks for the
/// dense, non-positive path: `'cholesky'` (sklearn `resolve_solver`,
/// `_ridge.py:830`). `Cholesky`/`Svd` resolve to themselves.
#[must_use]
fn resolve(self) -> RidgeSolver {
match self {
RidgeSolver::Auto => RidgeSolver::Cholesky,
other => other,
}
}
}
/// Ridge regression (L2-regularized least squares).
///
/// Adds an L2 penalty to the ordinary least squares objective, which
/// shrinks coefficients toward zero and can help with multicollinearity.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
#[derive(Debug, Clone)]
pub struct Ridge<F> {
/// Regularization strength. Larger values specify stronger
/// regularization.
pub alpha: F,
/// Optional per-target L2 penalties for the multi-output
/// `Fit<Array2, Array2>` path (sklearn accepts `alpha` as an array of
/// shape `(n_targets,)`, validated at
/// `sklearn/linear_model/_ridge.py:701-712`). When `Some`, each target
/// column `k` is fitted with its own penalty `alpha[k]`, overriding the
/// scalar [`alpha`](Self::alpha); this is mathematically identical to
/// fitting each target column with an independent scalar-`alpha` Ridge.
/// Default `None` (the scalar `alpha` applies to every target).
pub alpha_per_target: Option<Array1<F>>,
/// Whether to fit an intercept (bias) term.
pub fit_intercept: bool,
/// Whether `X` may be overwritten during fit (sklearn `copy_X`,
/// `_ridge.py:898`). ferrolearn's `fit` never mutates `x` (it reads
/// via `.mean_axis()`/`.outer_iter()`), so the observable
/// non-mutation contract holds regardless; the field is exposed for
/// ABI parity with sklearn. Default `true`, matching sklearn's
/// `copy_X=True` default (`_ridge.py:898`).
pub copy_x: bool,
/// Random seed for the `sag`/`saga` solvers (sklearn `random_state`,
/// `_ridge.py:903`). ferrolearn's default solver is deterministic
/// Cholesky (`solver='auto'`→`'cholesky'`), so this field is stored
/// for ABI parity and has no effect on the computed coefficients.
/// Default `None`, matching sklearn's `random_state=None` (`_ridge.py:903`).
pub random_state: Option<u64>,
/// Maximum number of iterations for iterative solvers (sklearn `max_iter`,
/// `_ridge.py:899`). Exposed for sklearn ABI parity; the implemented
/// direct Cholesky solver solves the normal equations in closed form with
/// no iteration, so this field is stored but has no effect on the computed
/// result. When an iterative solver is added (future REQ-8 #386), this
/// will control convergence. Default `None`, matching sklearn's default
/// (`_ridge.py:899`).
pub max_iter: Option<usize>,
/// Tolerance for iterative solvers (sklearn `tol`, `_ridge.py:900`).
/// Exposed for sklearn ABI parity; the implemented direct Cholesky solver
/// solves the normal equations in closed form with no iteration, so this
/// field is stored but has no effect on the computed result. When an
/// iterative solver is added (future REQ-8 #386), this will control
/// convergence. Default `1e-4`, matching sklearn's default (`_ridge.py:900`).
pub tol: F,
/// When `true`, constrain the fitted coefficients to be non-negative
/// (sklearn `positive`, `_ridge.py:902`/`:911`). sklearn solves the
/// non-negative ridge QP `min 0.5·‖X·w − y‖² + 0.5·alpha·‖w‖²` subject to
/// `w ≥ 0` via the L-BFGS-B solver (`_solve_lbfgs`, `_ridge.py:300`,
/// dispatched at `:923-928`); ferrolearn solves the same unique optimum
/// with projected coordinate descent. Default `false`, matching sklearn's
/// `positive=False` (`_ridge.py:902`); `positive=false` is byte-identical
/// to the unconstrained Cholesky path.
pub positive: bool,
/// Closed-form dense solver for the single-output unconstrained solve
/// (sklearn `solver`, `_ridge.py:830` `resolve_solver`). One of
/// [`RidgeSolver::Auto`] (default, resolves to `Cholesky` for the dense
/// path), [`RidgeSolver::Cholesky`], or [`RidgeSolver::Svd`]. All dense
/// solvers yield the IDENTICAL (unique) ridge solution; this only selects
/// the factorization used. Governs ONLY the single-output
/// (`Fit<Array2, Array1>`) unconstrained path — the multi-output and
/// `positive=true` paths ignore it. The iterative solvers
/// (`lsqr`/`sparse_cg`/`sag`/`saga`/`lbfgs`) are tracked as REQ-8b (#386).
/// Default [`RidgeSolver::Auto`], matching sklearn's `solver='auto'`
/// (`_ridge.py:903`).
pub solver: RidgeSolver,
}
impl<F: Float> Ridge<F> {
/// Create a new `Ridge` with default settings.
///
/// Defaults: `alpha = 1.0`, `fit_intercept = true`, `copy_x = true`,
/// `random_state = None`, `max_iter = None`, `tol = 1e-4`,
/// `positive = false`, `solver = Auto` — mirroring sklearn's ctor defaults
/// (`sklearn/linear_model/_ridge.py:895-903`).
#[must_use]
pub fn new() -> Self {
Self {
alpha: F::one(),
alpha_per_target: None,
fit_intercept: true,
copy_x: true,
random_state: None,
max_iter: None,
tol: F::from(1e-4).unwrap_or_else(F::epsilon),
positive: false,
solver: RidgeSolver::Auto,
}
}
/// Set the regularization strength.
#[must_use]
pub fn with_alpha(mut self, alpha: F) -> Self {
self.alpha = alpha;
self
}
/// Set per-target L2 penalties for the multi-output
/// `Fit<Array2, Array2>` path (sklearn array-valued `alpha`,
/// `sklearn/linear_model/_ridge.py:701-712`).
///
/// When set, each target column `k` of `Y` is fitted with its own penalty
/// `alphas[k]`, overriding the scalar [`alpha`](Self::alpha) on the
/// multi-output fit. This is mathematically identical to fitting each
/// target column with an independent scalar-`alpha` Ridge. The array
/// length must equal the number of target columns and every entry must be
/// non-negative (validated at fit time).
#[must_use]
pub fn with_alpha_per_target(mut self, alphas: Array1<F>) -> Self {
self.alpha_per_target = Some(alphas);
self
}
/// Set whether to fit an intercept term.
#[must_use]
pub fn with_fit_intercept(mut self, fit_intercept: bool) -> Self {
self.fit_intercept = fit_intercept;
self
}
/// Set the `copy_X` flag (sklearn `copy_X`, `_ridge.py:898`).
///
/// ferrolearn's fit never mutates `x`, so this is exposed for ABI
/// parity with sklearn and does not change the computed result.
#[must_use]
pub fn with_copy_x(mut self, copy_x: bool) -> Self {
self.copy_x = copy_x;
self
}
/// Set the `random_state` seed (sklearn `random_state`, `_ridge.py:903`).
///
/// Only used by sklearn's `sag`/`saga` solvers. ferrolearn's default
/// solver is deterministic Cholesky, so this is stored for ABI parity
/// and has no effect on the computed coefficients.
#[must_use]
pub fn with_random_state(mut self, random_state: Option<u64>) -> Self {
self.random_state = random_state;
self
}
/// Set the maximum number of iterations for iterative solvers (sklearn
/// `max_iter`, `_ridge.py:899`).
///
/// ferrolearn's direct Cholesky solver solves the normal equations in
/// closed form with no iteration, so this is stored for sklearn ABI
/// parity and does not affect the computed result. When an iterative
/// solver is added (future REQ-8 #386), this will take effect.
#[must_use]
pub fn with_max_iter(mut self, max_iter: Option<usize>) -> Self {
self.max_iter = max_iter;
self
}
/// Set the convergence tolerance for iterative solvers (sklearn `tol`,
/// `_ridge.py:900`).
///
/// ferrolearn's direct Cholesky solver solves the normal equations in
/// closed form with no iteration, so this is stored for sklearn ABI
/// parity and does not affect the computed result. When an iterative
/// solver is added (future REQ-8 #386), this will take effect.
#[must_use]
pub fn with_tol(mut self, tol: F) -> Self {
self.tol = tol;
self
}
/// Set whether to constrain the fitted coefficients to be non-negative
/// (sklearn `positive`, `_ridge.py:902`/`:911`).
///
/// When `true`, ferrolearn solves the non-negative ridge QP
/// `min 0.5·‖X·w − y‖² + 0.5·alpha·‖w‖²` subject to `w ≥ 0` via projected
/// coordinate descent — the same unique optimum sklearn reaches with its
/// L-BFGS-B solver (`_solve_lbfgs`, `_ridge.py:300`). `false` (default)
/// uses the unconstrained Cholesky path, byte-identical to today.
#[must_use]
pub fn with_positive(mut self, positive: bool) -> Self {
self.positive = positive;
self
}
/// Set the dense closed-form solver for the single-output unconstrained
/// solve (sklearn `solver`, `_ridge.py:830`).
///
/// [`RidgeSolver::Auto`] (default) resolves to [`RidgeSolver::Cholesky`]
/// for the dense path; [`RidgeSolver::Svd`] solves the same unique ridge
/// problem via the SVD of `X`. All dense solvers yield the identical
/// coefficients (the minimizer is unique); only the factorization differs.
/// This governs ONLY the single-output (`Fit<Array2, Array1>`)
/// unconstrained path.
#[must_use]
pub fn with_solver(mut self, solver: RidgeSolver) -> Self {
self.solver = solver;
self
}
}
impl<F: Float + Send + Sync + ScalarOperand + FromPrimitive + LinalgFloat + 'static> Ridge<F> {
/// Fit the Ridge regression model with optional per-sample weights.
///
/// Mirrors scikit-learn's `Ridge.fit(X, y, sample_weight=None)`
/// (`sklearn/linear_model/_ridge.py`). When `sample_weight` is `Some(w)`,
/// this solves the WEIGHTED ridge problem
/// `min Σᵢ wᵢ (yᵢ − xᵢ·coef)² + alpha·‖coef‖²`:
///
/// - `fit_intercept=true`: offsets are the WEIGHTED means
/// `x_off[j] = Σᵢ wᵢ·x[i,j] / Σwᵢ`, `y_off = Σᵢ wᵢ·yᵢ / Σwᵢ`. `X` and `y`
/// are centered by those offsets, each row is then rescaled by `√wᵢ`
/// (sklearn `_rescale_data`, `_ridge.py:682-688`), the cholesky ridge solve
/// runs on the rescaled design with the penalty `alpha` UNSCALED (because
/// `Xsᵀ·Xs == Xᵀ·W·X`), and `intercept = y_off − x_off·coef`.
/// - `fit_intercept=false`: no centering; each row is rescaled by `√wᵢ`, the
/// ridge solve runs on the rescaled `X`/`y`, and `intercept = 0`.
///
/// `sample_weight=None` is BYTE-IDENTICAL to [`Fit::fit`] (the unweighted
/// centering + `linalg::solve_ridge` path), which delegates here. The
/// `alpha=0` → OLS min-norm fallback is preserved (handled inside
/// `linalg::solve_ridge`).
///
/// # Errors
///
/// Returns [`FerroError::ShapeMismatch`] if the number of samples in `x` and
/// `y` (or, when provided, `sample_weight`) differ.
/// Returns [`FerroError::InvalidParameter`] if `alpha` is negative.
/// Returns [`FerroError::InsufficientSamples`] if there are no samples.
/// Returns [`FerroError::NumericalInstability`] if the weighted-offset
/// denominator (`Σwᵢ`) cannot be formed or a mean cannot be computed.
/// Solve the non-negative ridge problem
/// `min 0.5·‖A·w − b‖² + 0.5·alpha·‖w‖²` subject to `w ≥ 0` via projected
/// coordinate descent.
///
/// Mirrors the unique optimum sklearn reaches with its L-BFGS-B solver
/// (`_solve_lbfgs`, `sklearn/linear_model/_ridge.py:300`, objective
/// `0.5·‖Xw−y‖² + 0.5·alpha·‖w‖²` with bounds `[(0, inf)]`). For a smooth
/// strongly-convex QP with box constraints, coordinate descent with
/// per-coordinate projection converges to that exact optimum.
///
/// Returns `(coef, n_iter)`. A column with `‖A[:,j]‖² + alpha == 0` keeps
/// its coordinate at zero (no division by zero).
///
/// Inactive-constraint short-circuit (R-DEV-6): we first compute the
/// UNCONSTRAINED ridge optimum via the exact Cholesky solve
/// (`crate::linalg::solve_ridge`). If that succeeds AND every component is
/// `>= 0`, the positivity constraint is inactive, and by strict convexity +
/// the KKT conditions a feasible unconstrained minimizer IS the constrained
/// minimizer — so we return it as `(w_unc, 0)` (zero iterations, exact). This
/// matches a CONVERGED sklearn exactly: sklearn's default positive-Ridge
/// L-BFGS-B (`_solve_lbfgs`, `_ridge.py:300`) runs at `tol=1e-4` and
/// under-converges on interior optima by ~1e-4, whereas at `tol=1e-12` it
/// reaches the same true optimum ferrolearn computes here in closed form.
/// Otherwise (the Cholesky solve errored, or any component is `< 0` so the
/// constraint is active) we fall through to projected coordinate descent.
fn solve_nonneg_ridge(&self, a: &Array2<F>, b: &Array1<F>) -> (Array1<F>, usize) {
// Inactive-constraint short-circuit (R-DEV-6): the exact unconstrained
// Cholesky optimum, when feasible (all components >= 0), IS the
// constrained optimum (strict convexity + KKT). Return it exactly.
if let Ok(w_unc) = crate::linalg::solve_ridge(a, b, self.alpha)
&& w_unc.iter().all(|&v| v >= <F as num_traits::Zero>::zero())
{
return (w_unc, 0);
}
// Active constraint (or the Cholesky solve failed) — delegate to the
// shared projected-coordinate-descent kernel
// (`crate::linalg::nonneg_ridge_cd`) so `Ridge` and `RidgeClassifier`
// solve the non-negative ridge problem identically.
crate::linalg::nonneg_ridge_cd(a, b, self.alpha, self.max_iter.unwrap_or(1000), self.tol)
}
/// Solve the unconstrained dense ridge system on the given (already
/// centered / `√w`-rescaled) design, dispatching on the RESOLVED solver.
///
/// - [`RidgeSolver::Svd`] → `linalg::solve_ridge_svd` (`_solve_svd`,
/// `_ridge.py:200`).
/// - [`RidgeSolver::Cholesky`] / [`RidgeSolver::Auto`] →
/// `linalg::solve_ridge` (Cholesky normal equations, `_ridge.py:741`) —
/// byte-identical to the historic path.
///
/// `resolved` is `self.solver.resolve()` (precomputed by the caller). Every
/// dense solver returns the same unique minimizer; the choice only affects
/// the factorization.
fn solve_unconstrained(
&self,
resolved: RidgeSolver,
x: &Array2<F>,
y: &Array1<F>,
) -> Result<Array1<F>, FerroError> {
match resolved {
RidgeSolver::Svd => linalg::solve_ridge_svd(x, y, self.alpha),
// `Auto` is resolved to `Cholesky` before this point; matched here
// for exhaustiveness with byte-identical behaviour.
RidgeSolver::Cholesky | RidgeSolver::Auto => linalg::solve_ridge(x, y, self.alpha),
}
}
pub fn fit_with_sample_weight(
&self,
x: &Array2<F>,
y: &Array1<F>,
sample_weight: Option<&Array1<F>>,
) -> Result<FittedRidge<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(),
});
}
// `<F as num_traits::Zero>::zero()`: the `LinalgFloat` bound pulls
// `ferray::Element` (which also defines a `zero`) into scope, so a bare
// `F::zero()` is ambiguous between `Element` and `num_traits::Zero`.
if self.alpha < <F as num_traits::Zero>::zero() {
return Err(FerroError::InvalidParameter {
name: "alpha".into(),
reason: "must be non-negative".into(),
});
}
if n_samples == 0 {
return Err(FerroError::InsufficientSamples {
required: 1,
actual: 0,
context: "Ridge requires at least one sample".into(),
});
}
if let Some(w) = sample_weight
&& w.len() != n_samples
{
return Err(FerroError::ShapeMismatch {
expected: vec![n_samples],
actual: vec![w.len()],
context: "sample_weight length must match number of samples in X".into(),
});
}
// Non-finite input validation (#2259). sklearn `Ridge.fit` ->
// `self._validate_data(X, y, ...)` (`_ridge.py:1242`) keeps the default
// `force_all_finite=True`, so `check_array` rejects any NaN or +/-inf in
// X OR y with a `ValueError` BEFORE the solve. sklearn also validates
// `sample_weight` via `_check_sample_weight` (default `force_all_finite=
// True`), raising on a non-finite weight. `.iter().any(|v| !v.is_finite())`
// rejects both NaN and Inf (bounds-safe, no panic, R-CODE-2), matching
// the crate idiom (`linear_regression.rs`/`lasso.rs`). The finite path is
// byte-identical (the guard never fires on finite input). This is the
// shared single-output entry; `Fit::fit` delegates here with `None`.
if x.iter().any(|v| !v.is_finite()) {
return Err(FerroError::InvalidParameter {
name: "X".into(),
reason: "Input X contains NaN or infinity.".into(),
});
}
if y.iter().any(|v| !v.is_finite()) {
return Err(FerroError::InvalidParameter {
name: "y".into(),
reason: "Input y contains NaN or infinity.".into(),
});
}
if let Some(w) = sample_weight
&& w.iter().any(|v| !v.is_finite())
{
return Err(FerroError::InvalidParameter {
name: "sample_weight".into(),
reason: "Input sample_weight contains NaN or infinity.".into(),
});
}
// Resolve the dense solver once (sklearn `resolve_solver`,
// `_ridge.py:830`): `Auto` → `Cholesky` for the dense path. This is the
// value stored as the fitted `solver_`. It governs only the
// unconstrained dense solve; the `positive` path reports the same
// resolved dense value (the `solver` field does not select the
// constrained CD solver — see the `solver` field doc).
let resolved = self.solver.resolve();
match sample_weight {
None => {
// Unweighted path — identical to the original `Fit::fit` body.
if self.fit_intercept {
// Center the data to handle the intercept.
let x_mean =
x.mean_axis(Axis(0))
.ok_or_else(|| FerroError::NumericalInstability {
message: "failed to compute column means".into(),
})?;
let y_mean = y.mean().ok_or_else(|| FerroError::NumericalInstability {
message: "failed to compute target mean".into(),
})?;
let x_centered = x - &x_mean;
let y_centered = y - y_mean;
let (w, n_iter_) = if self.positive {
let (coef, iters) = self.solve_nonneg_ridge(&x_centered, &y_centered);
(coef, Some(iters))
} else {
(
self.solve_unconstrained(resolved, &x_centered, &y_centered)?,
None,
)
};
let intercept = y_mean - x_mean.dot(&w);
Ok(FittedRidge {
coefficients: w,
intercept,
n_iter_,
solver_: resolved,
})
} else {
let (w, n_iter_) = if self.positive {
let (coef, iters) = self.solve_nonneg_ridge(x, y);
(coef, Some(iters))
} else {
(self.solve_unconstrained(resolved, x, y)?, None)
};
Ok(FittedRidge {
coefficients: w,
// Disambiguate `Element::zero` vs `num_traits::Zero::zero`
// (both in scope under the `LinalgFloat` bound).
intercept: <F as num_traits::Zero>::zero(),
n_iter_,
solver_: resolved,
})
}
}
Some(w) => {
// Per-row √w factor (sklearn `_rescale_data`, `_ridge.py:682-688`).
let w_sqrt = w.mapv(<F as Float>::sqrt);
if self.fit_intercept {
// WEIGHTED centering: offsets are the weighted means
// x_off[j] = Σ wᵢ x[i,j] / Σ wᵢ, y_off = Σ wᵢ yᵢ / Σ wᵢ
// (sklearn `_preprocess_data` weighted `_average`).
let w_sum = w.sum();
if w_sum <= <F as num_traits::Zero>::zero() {
return Err(FerroError::NumericalInstability {
message: "sum of sample_weight must be positive to center".into(),
});
}
let mut x_off = Array1::<F>::zeros(n_features);
for (i, row) in x.outer_iter().enumerate() {
let wi = w[i];
x_off = &x_off + &row.mapv(|v| v * wi);
}
x_off.mapv_inplace(|v| v / w_sum);
let y_off = y
.iter()
.zip(w.iter())
.fold(<F as num_traits::Zero>::zero(), |acc, (&yi, &wi)| {
acc + wi * yi
})
/ w_sum;
// Center, then row-rescale by √w. Penalty `alpha` is UNSCALED:
// (√w·Xc)ᵀ(√w·Xc) == Xcᵀ·W·Xc, so the closed form
// (Xsᵀ·Xs + alpha·I)·coef = Xsᵀ·ys IS the weighted ridge solve.
let x_centered = x - &x_off;
let y_centered = y - y_off;
let x_scaled = &x_centered * &w_sqrt.view().insert_axis(Axis(1));
let y_scaled = &y_centered * &w_sqrt;
let (coef, n_iter_) = if self.positive {
let (c, iters) = self.solve_nonneg_ridge(&x_scaled, &y_scaled);
(c, Some(iters))
} else {
(
self.solve_unconstrained(resolved, &x_scaled, &y_scaled)?,
None,
)
};
let intercept = y_off - x_off.dot(&coef);
Ok(FittedRidge {
coefficients: coef,
intercept,
n_iter_,
solver_: resolved,
})
} else {
// No centering; just √w row-rescaling, intercept 0.
let x_scaled = x * &w_sqrt.view().insert_axis(Axis(1));
let y_scaled = y * &w_sqrt;
let (coef, n_iter_) = if self.positive {
let (c, iters) = self.solve_nonneg_ridge(&x_scaled, &y_scaled);
(c, Some(iters))
} else {
(
self.solve_unconstrained(resolved, &x_scaled, &y_scaled)?,
None,
)
};
Ok(FittedRidge {
coefficients: coef,
intercept: <F as num_traits::Zero>::zero(),
n_iter_,
solver_: resolved,
})
}
}
}
}
}
impl<F: Float> Default for Ridge<F> {
fn default() -> Self {
Self::new()
}
}
/// Fitted Ridge regression model.
///
/// Stores the learned coefficients and intercept. Implements [`Predict`]
/// to generate predictions and [`HasCoefficients`] for introspection.
#[derive(Debug, Clone)]
pub struct FittedRidge<F> {
/// Learned coefficient vector (one per feature).
coefficients: Array1<F>,
/// Learned intercept (bias) term.
intercept: F,
/// Number of iterations run by an iterative solver, or `None` for direct
/// solvers (sklearn `n_iter_`, `_ridge.py:994`). ferrolearn's Cholesky
/// solver is direct (no iteration), so this is always `None` — matching
/// sklearn's behaviour when `solver='cholesky'` or `solver='svd'` resolve
/// the normal equations in closed form.
n_iter_: Option<usize>,
/// The RESOLVED dense solver actually used for the coefficient solve
/// (sklearn `solver_`, the resolution of `solver='auto'`,
/// `_ridge.py:830`/`:994`). `Auto` resolves to [`RidgeSolver::Cholesky`]
/// for the dense path; an explicit `Cholesky`/`Svd` resolves to itself.
solver_: RidgeSolver,
}
impl<F> FittedRidge<F> {
/// Return the number of iterations run by an iterative solver, or `None`
/// for the direct Cholesky solver (sklearn `n_iter_`, `_ridge.py:994`).
///
/// ferrolearn implements only the direct Cholesky path, so this is always
/// `None`. When an iterative solver is added (future REQ-8 #386), it will
/// return `Some(n)`.
#[must_use]
pub fn n_iter(&self) -> Option<usize> {
self.n_iter_
}
/// Return the RESOLVED dense solver used for the coefficient solve
/// (sklearn `solver_`, `_ridge.py:830`/`:994`).
///
/// [`RidgeSolver::Auto`] resolves to [`RidgeSolver::Cholesky`] for the
/// dense path, so a default-`solver` fit reports `Cholesky`; an explicit
/// `Svd` reports `Svd`.
#[must_use]
pub fn solver(&self) -> RidgeSolver {
self.solver_
}
}
impl<F: Float + Send + Sync + ScalarOperand + FromPrimitive + LinalgFloat + 'static>
Fit<Array2<F>, Array1<F>> for Ridge<F>
{
type Fitted = FittedRidge<F>;
type Error = FerroError;
/// Fit the Ridge regression model using Cholesky decomposition.
///
/// Solves `(X^T X + alpha * I)^{-1} X^T y`.
///
/// # Errors
///
/// Returns [`FerroError::ShapeMismatch`] if the number of samples in
/// `x` and `y` differ.
/// Returns [`FerroError::InvalidParameter`] if `alpha` is negative.
fn fit(&self, x: &Array2<F>, y: &Array1<F>) -> Result<FittedRidge<F>, FerroError> {
// Unweighted ridge is the `sample_weight=None` arm of the weighted fit;
// delegating keeps the None path byte-identical to the historic body
// (centering + `solve_ridge`), mirroring sklearn's single `fit` entry
// (`_ridge.py`, `sample_weight=None` default).
self.fit_with_sample_weight(x, y, None)
}
}
impl<F: Float + Send + Sync + ScalarOperand + 'static> Predict<Array2<F>> for FittedRidge<F> {
type Output = Array1<F>;
type Error = FerroError;
/// Predict target values for the given feature matrix.
///
/// Computes `X @ coefficients + intercept`.
///
/// # Errors
///
/// Returns [`FerroError::ShapeMismatch`] if the number of features
/// does not match the fitted model.
fn predict(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
let n_features = x.ncols();
if n_features != self.coefficients.len() {
return Err(FerroError::ShapeMismatch {
expected: vec![self.coefficients.len()],
actual: vec![n_features],
context: "number of features must match fitted model".into(),
});
}
let preds = x.dot(&self.coefficients) + self.intercept;
Ok(preds)
}
}
impl<F: Float + Send + Sync + ScalarOperand + 'static> HasCoefficients<F> for FittedRidge<F> {
fn coefficients(&self) -> &Array1<F> {
&self.coefficients
}
fn intercept(&self) -> F {
self.intercept
}
}
/// Fitted multi-output Ridge regression model.
///
/// Companion to [`FittedRidge`] for the case where `Y` has multiple
/// target columns. Stores a `(n_features, n_targets)` coefficient matrix
/// and a per-target intercept vector. The Cholesky factor of
/// `X^T X + alpha * I` is computed once during [`Ridge::fit`] and shared
/// across all targets, so multi-output fitting costs the same `O(p^3)`
/// factorization as the single-output path.
#[derive(Debug, Clone)]
pub struct FittedRidgeMulti<F> {
/// Learned coefficients, shape `(n_features, n_targets)`.
coefficients: Array2<F>,
/// Per-target intercept vector, length `n_targets`. Filled with
/// zeros when `fit_intercept = false`.
intercepts: Array1<F>,
}
impl<F: Float> FittedRidgeMulti<F> {
/// Borrow the learned coefficient matrix `(n_features, n_targets)`.
#[must_use]
pub fn coefficients(&self) -> &Array2<F> {
&self.coefficients
}
/// Borrow the per-target intercept vector `(n_targets,)`.
#[must_use]
pub fn intercepts(&self) -> &Array1<F> {
&self.intercepts
}
}
impl<F: Float + Send + Sync + ScalarOperand + FromPrimitive + LinalgFloat + 'static>
Fit<Array2<F>, Array2<F>> for Ridge<F>
{
type Fitted = FittedRidgeMulti<F>;
type Error = FerroError;
/// Fit the multi-output Ridge regression model using a single
/// shared Cholesky factorization across all `Y` columns.
///
/// Solves `(X^T X + alpha * I)^{-1} X^T Y` where `Y` is
/// `(n_samples, n_targets)`.
///
/// # Errors
///
/// Returns [`FerroError::ShapeMismatch`] if the number of samples in
/// `x` and `y` differ.
/// Returns [`FerroError::InvalidParameter`] if `alpha` is negative.
fn fit(&self, x: &Array2<F>, y: &Array2<F>) -> Result<FittedRidgeMulti<F>, FerroError> {
let (n_samples, _n_features) = x.dim();
if n_samples != y.nrows() {
return Err(FerroError::ShapeMismatch {
expected: vec![n_samples],
actual: vec![y.nrows()],
context: "y rows must match number of samples in X".into(),
});
}
if n_samples == 0 {
return Err(FerroError::InsufficientSamples {
required: 1,
actual: 0,
context: "Ridge requires at least one sample".into(),
});
}
// Non-finite input validation (#2259) — SEPARATE multi-output arm.
// `Fit<Array2, Array2>::fit` does NOT delegate to
// `fit_with_sample_weight`, so the same `_validate_data(force_all_finite=
// True)` reject-at-fit contract (`_ridge.py:1242`, `multi_output=True`)
// is enforced here independently before centering/solve.
if x.iter().any(|v| !v.is_finite()) {
return Err(FerroError::InvalidParameter {
name: "X".into(),
reason: "Input X contains NaN or infinity.".into(),
});
}
if y.iter().any(|v| !v.is_finite()) {
return Err(FerroError::InvalidParameter {
name: "y".into(),
reason: "Input y contains NaN or infinity.".into(),
});
}
let n_targets = y.ncols();
// Per-target alpha array (sklearn array-valued `alpha`,
// `sklearn/linear_model/_ridge.py:701-712`) — validate the array and
// solve each target column independently with its own penalty. When
// `None`, fall through to the byte-identical scalar-alpha path below.
if let Some(alphas) = &self.alpha_per_target {
if alphas.len() != n_targets {
return Err(FerroError::ShapeMismatch {
expected: vec![n_targets],
actual: vec![alphas.len()],
context: "alpha array length must equal number of targets".into(),
});
}
for &a in alphas.iter() {
if a < <F as num_traits::Zero>::zero() {
return Err(FerroError::InvalidParameter {
name: "alpha".into(),
reason: "must be non-negative".into(),
});
}
}
let n_features = x.ncols();
let mut coefficients = Array2::<F>::zeros((n_features, n_targets));
let mut intercepts = Array1::<F>::zeros(n_targets);
if self.fit_intercept {
// Center once (identical to the scalar path), then solve each
// centered target column with its own alpha.
let x_mean =
x.mean_axis(Axis(0))
.ok_or_else(|| FerroError::NumericalInstability {
message: "failed to compute column means of X".into(),
})?;
let y_mean =
y.mean_axis(Axis(0))
.ok_or_else(|| FerroError::NumericalInstability {
message: "failed to compute column means of Y".into(),
})?;
let x_centered = x - &x_mean;
let y_centered = y - &y_mean;
for k in 0..n_targets {
let y_col = y_centered.column(k).to_owned();
let w_col = linalg::solve_ridge(&x_centered, &y_col, alphas[k])?;
intercepts[k] = y_mean[k] - x_mean.dot(&w_col);
coefficients.column_mut(k).assign(&w_col);
}
} else {
for k in 0..n_targets {
let y_col = y.column(k).to_owned();
let w_col = linalg::solve_ridge(x, &y_col, alphas[k])?;
coefficients.column_mut(k).assign(&w_col);
}
}
return Ok(FittedRidgeMulti {
coefficients,
intercepts,
});
}
if self.alpha < <F as num_traits::Zero>::zero() {
return Err(FerroError::InvalidParameter {
name: "alpha".into(),
reason: "must be non-negative".into(),
});
}
if self.fit_intercept {
// Center the data to handle the intercept (per-target).
let x_mean = x
.mean_axis(Axis(0))
.ok_or_else(|| FerroError::NumericalInstability {
message: "failed to compute column means of X".into(),
})?;
let y_mean = y
.mean_axis(Axis(0))
.ok_or_else(|| FerroError::NumericalInstability {
message: "failed to compute column means of Y".into(),
})?;
let x_centered = x - &x_mean;
let y_centered = y - &y_mean;
let w = linalg::solve_ridge_multi(&x_centered, &y_centered, self.alpha)?;
// intercept[k] = y_mean[k] - x_mean · w[:, k]
let mut intercepts = Array1::<F>::zeros(n_targets);
for k in 0..n_targets {
let col = w.column(k);
let dot = x_mean.dot(&col);
intercepts[k] = y_mean[k] - dot;
}
Ok(FittedRidgeMulti {
coefficients: w,
intercepts,
})
} else {
let w = linalg::solve_ridge_multi(x, y, self.alpha)?;
Ok(FittedRidgeMulti {
coefficients: w,
intercepts: Array1::<F>::zeros(n_targets),
})
}
}
}
impl<F: Float + Send + Sync + ScalarOperand + 'static> Predict<Array2<F>> for FittedRidgeMulti<F> {
type Output = Array2<F>;
type Error = FerroError;
/// Predict target values for the given feature matrix.
///
/// Computes `X @ coefficients + intercepts` and returns an
/// `(n_samples, n_targets)` array.
///
/// # Errors
///
/// Returns [`FerroError::ShapeMismatch`] if the number of features
/// does not match the fitted model.
fn predict(&self, x: &Array2<F>) -> Result<Array2<F>, FerroError> {
let n_features = x.ncols();
if n_features != self.coefficients.nrows() {
return Err(FerroError::ShapeMismatch {
expected: vec![self.coefficients.nrows()],
actual: vec![n_features],
context: "number of features must match fitted model".into(),
});
}
let mut preds = x.dot(&self.coefficients);
// Broadcast-add per-target intercepts.
for (k, &b) in self.intercepts.iter().enumerate() {
let mut col = preds.column_mut(k);
col.mapv_inplace(|v| v + b);
}
Ok(preds)
}
}
// Pipeline integration.
impl<F> PipelineEstimator<F> for Ridge<F>
where
F: Float + FromPrimitive + ScalarOperand + LinalgFloat + Send + Sync + 'static,
{
fn fit_pipeline(
&self,
x: &Array2<F>,
y: &Array1<F>,
) -> Result<Box<dyn FittedPipelineEstimator<F>>, FerroError> {
let fitted = self.fit(x, y)?;
Ok(Box::new(fitted))
}
}
impl<F> FittedPipelineEstimator<F> for FittedRidge<F>
where
F: Float + ScalarOperand + Send + Sync + 'static,
{
fn predict_pipeline(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
self.predict(x)
}
}
#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;
use ndarray::array;
#[test]
fn test_ridge_no_regularization() {
// With alpha=0, Ridge should behave like OLS.
let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
let y = array![3.0, 5.0, 7.0, 9.0, 11.0];
let model = Ridge::<f64>::new().with_alpha(0.0);
let fitted = model.fit(&x, &y).unwrap();
assert_relative_eq!(fitted.coefficients()[0], 2.0, epsilon = 1e-8);
assert_relative_eq!(fitted.intercept(), 1.0, epsilon = 1e-8);
}
#[test]
fn test_ridge_shrinks_coefficients() {
let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
let y = array![3.0, 5.0, 7.0, 9.0, 11.0];
let model_low = Ridge::<f64>::new().with_alpha(0.01);
let model_high = Ridge::<f64>::new().with_alpha(100.0);
let fitted_low = model_low.fit(&x, &y).unwrap();
let fitted_high = model_high.fit(&x, &y).unwrap();
// Higher alpha should shrink coefficients more.
assert!(fitted_high.coefficients()[0].abs() < fitted_low.coefficients()[0].abs());
}
#[test]
fn test_ridge_no_intercept() {
let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
let y = array![2.0, 4.0, 6.0, 8.0];
let model = Ridge::<f64>::new()
.with_alpha(0.0)
.with_fit_intercept(false);
let fitted = model.fit(&x, &y).unwrap();
assert_relative_eq!(fitted.coefficients()[0], 2.0, epsilon = 1e-10);
assert_relative_eq!(fitted.intercept(), 0.0, epsilon = 1e-10);
}
#[test]
fn test_ridge_negative_alpha() {
let x = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
let y = array![1.0, 2.0, 3.0];
let model = Ridge::<f64>::new().with_alpha(-1.0);
let result = model.fit(&x, &y);
assert!(result.is_err());
}
#[test]
fn test_ridge_shape_mismatch() {
let x = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
let y = array![1.0, 2.0];
let model = Ridge::<f64>::new();
let result = model.fit(&x, &y);
assert!(result.is_err());
}
#[test]
fn test_ridge_predict() {
let x =
Array2::from_shape_vec((4, 2), vec![1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0]).unwrap();
let y = array![1.0, 2.0, 3.0, 6.0];
let model = Ridge::<f64>::new().with_alpha(0.01);
let fitted = model.fit(&x, &y).unwrap();
let preds = fitted.predict(&x).unwrap();
assert_eq!(preds.len(), 4);
}
#[test]
fn test_ridge_pipeline_integration() {
let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
let y = array![3.0, 5.0, 7.0, 9.0];
let model = Ridge::<f64>::new();
let fitted = model.fit_pipeline(&x, &y).unwrap();
let preds = fitted.predict_pipeline(&x).unwrap();
assert_eq!(preds.len(), 4);
}
#[test]
fn test_ridge_has_coefficients() {
let x = Array2::from_shape_vec((3, 2), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
let y = array![1.0, 2.0, 3.0];
let model = Ridge::<f64>::new();
let fitted = model.fit(&x, &y).unwrap();
assert_eq!(fitted.coefficients().len(), 2);
}
#[test]
fn ridge_fit_sample_weight_with_intercept_matches_sklearn() -> Result<(), FerroError> {
// Live sklearn 1.5.2 oracle (WEIGHTED ridge, alpha=1, fit_intercept=True):
// cd /tmp && python3 -c "import numpy as np; \
// from sklearn.linear_model import Ridge; \
// X=np.array([[1.,2.],[2.,1.],[3.,4.],[4.,3.],[5.,5.]]); \
// y=np.array([3.0,2.5,7.1,6.0,11.2]); w=np.array([1.,4.,1.,1.,3.]); \
// m=Ridge(alpha=1.0).fit(X,y,sample_weight=w); \
// print([round(c,10) for c in m.coef_], round(m.intercept_,10))"
// -> [0.9233502538, 1.39678511] -0.8033840948
let x = Array2::from_shape_vec((5, 2), vec![1., 2., 2., 1., 3., 4., 4., 3., 5., 5.])
.map_err(|e| FerroError::ShapeMismatch {
expected: vec![5, 2],
actual: vec![],
context: e.to_string(),
})?;
let y = array![3.0, 2.5, 7.1, 6.0, 11.2];
let w = array![1.0, 4.0, 1.0, 1.0, 3.0];
let model = Ridge::<f64>::new().with_alpha(1.0);
let fitted = model.fit_with_sample_weight(&x, &y, Some(&w))?;
assert_relative_eq!(fitted.coefficients()[0], 0.923_350_253_8, epsilon = 1e-7);
assert_relative_eq!(fitted.coefficients()[1], 1.396_785_11, epsilon = 1e-7);
assert_relative_eq!(fitted.intercept(), -0.803_384_094_8, epsilon = 1e-7);
// Non-tautological: the weighted result MUST differ from the unweighted
// fit (oracle unweighted coef_ [0.8228070175, 1.3561403509]).
let unweighted = model.fit(&x, &y)?;
assert_relative_eq!(
unweighted.coefficients()[0],
0.822_807_017_5,
epsilon = 1e-7
);
assert_relative_eq!(
unweighted.coefficients()[1],
1.356_140_350_9,
epsilon = 1e-7
);
assert!((fitted.coefficients()[0] - unweighted.coefficients()[0]).abs() > 1e-3);
assert!((fitted.intercept() - unweighted.intercept()).abs() > 1e-3);
Ok(())
}
#[test]
fn ridge_fit_sample_weight_no_intercept_matches_sklearn() -> Result<(), FerroError> {
// Live sklearn 1.5.2 oracle (WEIGHTED ridge, alpha=2, fit_intercept=False):
// cd /tmp && python3 -c "import numpy as np; \
// from sklearn.linear_model import Ridge; \
// X=np.array([[1.,2.],[2.,1.],[3.,4.],[4.,3.],[5.,5.]]); \
// y=np.array([3.0,2.5,7.1,6.0,11.2]); w=np.array([1.,4.,1.,1.,3.]); \
// m=Ridge(alpha=2.0,fit_intercept=False).fit(X,y,sample_weight=w); \
// print([round(c,10) for c in m.coef_])"
// -> [0.7273779983, 1.3737799835]
let x = Array2::from_shape_vec((5, 2), vec![1., 2., 2., 1., 3., 4., 4., 3., 5., 5.])
.map_err(|e| FerroError::ShapeMismatch {
expected: vec![5, 2],
actual: vec![],
context: e.to_string(),
})?;
let y = array![3.0, 2.5, 7.1, 6.0, 11.2];
let w = array![1.0, 4.0, 1.0, 1.0, 3.0];
let model = Ridge::<f64>::new()
.with_alpha(2.0)
.with_fit_intercept(false);
let fitted = model.fit_with_sample_weight(&x, &y, Some(&w))?;
assert_relative_eq!(fitted.coefficients()[0], 0.727_377_998_3, epsilon = 1e-7);
assert_relative_eq!(fitted.coefficients()[1], 1.373_779_983_5, epsilon = 1e-7);
assert_eq!(fitted.intercept(), 0.0);
Ok(())
}
#[test]
fn ridge_fit_none_sample_weight_equals_unweighted() -> Result<(), FerroError> {
// Regression guard: the `None` path is BYTE-IDENTICAL to `fit`.
let x = Array2::from_shape_vec((5, 2), vec![1., 2., 2., 1., 3., 4., 4., 3., 5., 5.])
.map_err(|e| FerroError::ShapeMismatch {
expected: vec![5, 2],
actual: vec![],
context: e.to_string(),
})?;
let y = array![3.0, 2.5, 7.1, 6.0, 11.2];
let model = Ridge::<f64>::new().with_alpha(1.0);
let via_fit = model.fit(&x, &y)?;
let via_none = model.fit_with_sample_weight(&x, &y, None)?;
assert_eq!(
via_fit.coefficients()[0].to_bits(),
via_none.coefficients()[0].to_bits()
);
assert_eq!(
via_fit.coefficients()[1].to_bits(),
via_none.coefficients()[1].to_bits()
);
assert_eq!(
via_fit.intercept().to_bits(),
via_none.intercept().to_bits()
);
// Same for fit_intercept=false.
let model_ni = Ridge::<f64>::new()
.with_alpha(1.0)
.with_fit_intercept(false);
let via_fit_ni = model_ni.fit(&x, &y)?;
let via_none_ni = model_ni.fit_with_sample_weight(&x, &y, None)?;
assert_eq!(
via_fit_ni.coefficients()[0].to_bits(),
via_none_ni.coefficients()[0].to_bits()
);
assert_eq!(
via_fit_ni.intercept().to_bits(),
via_none_ni.intercept().to_bits()
);
Ok(())
}
// -- copy_x / random_state ABI-parity ----------------------------------
#[test]
fn ridge_copy_x_random_state_defaults_and_builders() -> Result<(), FerroError> {
// Defaults mirror sklearn Ridge(copy_X=True, random_state=None)
// (`sklearn/linear_model/_ridge.py:898`/`:903`).
assert!(Ridge::<f64>::new().copy_x, "copy_x default must be true");
assert_eq!(
Ridge::<f64>::new().random_state,
None,
"random_state default must be None"
);
// Builders store the supplied value.
assert!(!Ridge::<f64>::new().with_copy_x(false).copy_x);
assert!(Ridge::<f64>::new().with_copy_x(true).copy_x);
assert_eq!(
Ridge::<f64>::new().with_random_state(Some(42)).random_state,
Some(42)
);
assert_eq!(
Ridge::<f64>::new().with_random_state(None).random_state,
None
);
// No behavior change: fit produces byte-identical coef_/intercept_
// regardless of copy_x or random_state (deterministic Cholesky
// solver is unaffected by either param).
//
// Live sklearn oracle (alpha=1.0, fit_intercept=true):
// python3 -c "import numpy as np; from sklearn.linear_model import Ridge;
// X=np.array([[1.,2.],[3.,4.],[5.,6.],[7.,8.],[9.,10.]]);
// y=np.array([1.,2.,3.,4.,5.]);
// m=Ridge(alpha=1.0).fit(X,y); print(m.coef_.tolist(), m.intercept_)"
// -> [0.07692307692307693, 0.4230769230769231] 0.0
let x = Array2::from_shape_vec((5, 2), vec![1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])
.map_err(|e| FerroError::ShapeMismatch {
expected: vec![5, 2],
actual: vec![],
context: e.to_string(),
})?;
let y = array![1., 2., 3., 4., 5.];
let ref_fitted = Ridge::<f64>::new().with_alpha(1.0).fit(&x, &y)?;
let copy_false = Ridge::<f64>::new()
.with_alpha(1.0)
.with_copy_x(false)
.fit(&x, &y)?;
let rs_some = Ridge::<f64>::new()
.with_alpha(1.0)
.with_random_state(Some(42))
.fit(&x, &y)?;
// Byte-identical (deterministic Cholesky, no mutation).
assert_eq!(
ref_fitted.coefficients()[0].to_bits(),
copy_false.coefficients()[0].to_bits()
);
assert_eq!(
ref_fitted.coefficients()[1].to_bits(),
copy_false.coefficients()[1].to_bits()
);
assert_eq!(
ref_fitted.intercept().to_bits(),
copy_false.intercept().to_bits()
);
assert_eq!(
ref_fitted.coefficients()[0].to_bits(),
rs_some.coefficients()[0].to_bits()
);
assert_eq!(
ref_fitted.coefficients()[1].to_bits(),
rs_some.coefficients()[1].to_bits()
);
assert_eq!(
ref_fitted.intercept().to_bits(),
rs_some.intercept().to_bits()
);
Ok(())
}
// -- Multi-output Ridge -------------------------------------------------
#[test]
fn test_ridge_multi_recovers_two_targets_with_zero_alpha() {
// Two synthetic targets sharing the same features:
// y1 = 2*x1 - 3*x2 + 5
// y2 = 3*x1 + x2 - 1
let n = 50;
let x_data: Vec<f64> = (0..n)
.flat_map(|i| {
let i = i as f64;
[i / 10.0, (i / 7.0).sin()]
})
.collect();
let x = Array2::from_shape_vec((n, 2), x_data).unwrap();
let y_data: Vec<f64> = (0..n)
.flat_map(|i| {
let x1 = i as f64 / 10.0;
let x2 = (i as f64 / 7.0).sin();
[2.0 * x1 - 3.0 * x2 + 5.0, 3.0 * x1 + x2 - 1.0]
})
.collect();
let y = Array2::from_shape_vec((n, 2), y_data).unwrap();
let model = Ridge::<f64>::new().with_alpha(1e-8);
let fitted: FittedRidgeMulti<f64> = model.fit(&x, &y).unwrap();
let coef = fitted.coefficients();
assert_eq!(coef.shape(), &[2, 2]);
// Target 0
assert_relative_eq!(coef[[0, 0]], 2.0, epsilon = 1e-4);
assert_relative_eq!(coef[[1, 0]], -3.0, epsilon = 1e-4);
// Target 1
assert_relative_eq!(coef[[0, 1]], 3.0, epsilon = 1e-4);
assert_relative_eq!(coef[[1, 1]], 1.0, epsilon = 1e-4);
// Intercepts
assert_relative_eq!(fitted.intercepts()[0], 5.0, epsilon = 1e-4);
assert_relative_eq!(fitted.intercepts()[1], -1.0, epsilon = 1e-4);
}
#[test]
fn test_ridge_multi_no_intercept() {
// y = X @ B with no bias; verify intercepts come out zero and
// coefficients match the OLS solve.
let x = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
let y =
Array2::from_shape_vec((4, 2), vec![2.0, 4.0, 4.0, 8.0, 6.0, 12.0, 8.0, 16.0]).unwrap();
let model = Ridge::<f64>::new()
.with_alpha(0.0)
.with_fit_intercept(false);
let fitted = model.fit(&x, &y).unwrap();
// y[:, 0] = 2*x, y[:, 1] = 4*x
assert_relative_eq!(fitted.coefficients()[[0, 0]], 2.0, epsilon = 1e-8);
assert_relative_eq!(fitted.coefficients()[[0, 1]], 4.0, epsilon = 1e-8);
assert_relative_eq!(fitted.intercepts()[0], 0.0, epsilon = 1e-12);
assert_relative_eq!(fitted.intercepts()[1], 0.0, epsilon = 1e-12);
}
#[test]
fn test_ridge_multi_shrinks_with_large_alpha() {
// Heavy regularization should pull all targets toward zero
// coefficients (intercepts may still recover the means).
let n = 20;
let x = Array2::from_shape_vec((n, 1), (0..n).map(|i| i as f64).collect()).unwrap();
let y_data: Vec<f64> = (0..n)
.flat_map(|i| [(i as f64) * 10.0, (i as f64) * 5.0])
.collect();
let y = Array2::from_shape_vec((n, 2), y_data).unwrap();
let model = Ridge::<f64>::new()
.with_alpha(1e6)
.with_fit_intercept(false);
let fitted = model.fit(&x, &y).unwrap();
assert!(fitted.coefficients()[[0, 0]].abs() < 1.0);
assert!(fitted.coefficients()[[0, 1]].abs() < 1.0);
}
#[test]
fn test_ridge_multi_predict_round_trips_training_data() {
// Verify Fit→Predict round-trip on the training data: the model
// should reproduce y up to the regularization-induced bias.
let n = 40;
let x_data: Vec<f64> = (0..n)
.flat_map(|i| [(i as f64) / 10.0, ((i as f64) / 3.0).cos()])
.collect();
let x = Array2::from_shape_vec((n, 2), x_data).unwrap();
let y_data: Vec<f64> = (0..n)
.flat_map(|i| {
let x1 = (i as f64) / 10.0;
let x2 = ((i as f64) / 3.0).cos();
[1.5 * x1 + 0.7 * x2 + 0.2, -0.5 * x1 + 2.0 * x2 - 1.0]
})
.collect();
let y = Array2::from_shape_vec((n, 2), y_data).unwrap();
let model = Ridge::<f64>::new().with_alpha(1e-6);
let fitted = model.fit(&x, &y).unwrap();
let y_hat = fitted.predict(&x).unwrap();
assert_eq!(y_hat.shape(), y.shape());
// Maximum element-wise error across both targets.
let max_err = y_hat
.iter()
.zip(y.iter())
.map(|(a, b)| (a - b).abs())
.fold(0.0_f64, f64::max);
assert!(max_err < 1e-3, "max_err = {max_err}");
}
#[test]
fn test_ridge_multi_shape_mismatch() {
// 5-sample X but 3-sample Y — should error.
let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
let y = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
let model = Ridge::<f64>::new();
let err = <Ridge<f64> as Fit<Array2<f64>, Array2<f64>>>::fit(&model, &x, &y).unwrap_err();
matches!(err, FerroError::ShapeMismatch { .. });
}
#[test]
fn test_ridge_multi_single_target_matches_single_output_path() {
// The same data routed through the single-output and the
// multi-output Fit impls should produce coefficients that agree
// to within numerical noise. This pins the parallel paths from
// drifting apart over time.
let n = 30;
let x = Array2::from_shape_vec(
(n, 3),
(0..n)
.flat_map(|i| {
let i = i as f64;
[i / 10.0, (i / 5.0).sin(), (i / 11.0).cos()]
})
.collect(),
)
.unwrap();
let y_1d: Array1<f64> = x
.rows()
.into_iter()
.map(|r| 2.0 * r[0] - r[1] + 0.5 * r[2] + 3.0)
.collect();
let y_2d = Array2::from_shape_vec((n, 1), y_1d.to_vec()).unwrap();
let model = Ridge::<f64>::new().with_alpha(0.1);
let single: FittedRidge<f64> = model.fit(&x, &y_1d).unwrap();
let multi: FittedRidgeMulti<f64> = model.fit(&x, &y_2d).unwrap();
// Coefficients agree element-wise.
for j in 0..3 {
assert_relative_eq!(
single.coefficients()[j],
multi.coefficients()[[j, 0]],
epsilon = 1e-10
);
}
// Intercepts agree.
assert_relative_eq!(single.intercept(), multi.intercepts()[0], epsilon = 1e-10);
}
// -- max_iter / tol / n_iter_ ABI-parity (REQ-10) ----------------------
#[test]
fn ridge_max_iter_tol_niter_defaults_and_builders() -> Result<(), FerroError> {
// Defaults mirror sklearn Ridge(max_iter=None, tol=1e-4)
// (`sklearn/linear_model/_ridge.py:899-900`).
assert_eq!(
Ridge::<f64>::new().max_iter,
None,
"max_iter default must be None"
);
assert!(
(Ridge::<f64>::new().tol - 1e-4_f64).abs() < 1e-12,
"tol default must be 1e-4"
);
// Builders store the supplied value.
assert_eq!(
Ridge::<f64>::new().with_max_iter(Some(500)).max_iter,
Some(500)
);
assert_eq!(Ridge::<f64>::new().with_max_iter(None).max_iter, None);
assert!(
(Ridge::<f64>::new().with_tol(1e-6_f64).tol - 1e-6_f64).abs() < 1e-18,
"with_tol must store the supplied value"
);
// n_iter_ is always None for the direct Cholesky solver, matching
// sklearn: Ridge(alpha=1.0).fit(X,y).n_iter_ is None when
// solver='auto' resolves to 'cholesky' (`_ridge.py:994`).
let x = Array2::from_shape_vec((4, 2), vec![1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 2.0])
.map_err(|e| FerroError::ShapeMismatch {
expected: vec![4, 2],
actual: vec![],
context: e.to_string(),
})?;
let y = array![1.0, 2.0, 3.0, 6.0];
let fitted = Ridge::<f64>::new().fit(&x, &y)?;
assert_eq!(
fitted.n_iter(),
None,
"n_iter_ must be None for direct Cholesky"
);
// No behavior change: fit produces byte-identical coef_/intercept_
// regardless of max_iter or tol (direct Cholesky is unaffected).
//
// Live sklearn oracle (alpha=1.0, fit_intercept=true):
// python3 -c "import numpy as np; from sklearn.linear_model import Ridge;
// X=np.array([[1.,0.],[0.,1.],[1.,1.],[2.,2.]]);
// y=np.array([1.,2.,3.,6.]);
// m=Ridge(alpha=1.0).fit(X,y); print(m.coef_.tolist(), m.intercept_)"
// -> [0.875, 1.375] 0.75
let ref_fitted = Ridge::<f64>::new().with_alpha(1.0).fit(&x, &y)?;
let with_max_iter = Ridge::<f64>::new()
.with_alpha(1.0)
.with_max_iter(Some(500))
.fit(&x, &y)?;
let with_tol = Ridge::<f64>::new()
.with_alpha(1.0)
.with_tol(1e-6_f64)
.fit(&x, &y)?;
// Byte-identical (deterministic Cholesky, max_iter/tol are no-ops).
for j in 0..2 {
assert_eq!(
ref_fitted.coefficients()[j].to_bits(),
with_max_iter.coefficients()[j].to_bits(),
"coef_[{j}] must be byte-identical regardless of max_iter"
);
assert_eq!(
ref_fitted.coefficients()[j].to_bits(),
with_tol.coefficients()[j].to_bits(),
"coef_[{j}] must be byte-identical regardless of tol"
);
}
assert_eq!(
ref_fitted.intercept().to_bits(),
with_max_iter.intercept().to_bits()
);
assert_eq!(
ref_fitted.intercept().to_bits(),
with_tol.intercept().to_bits()
);
assert_eq!(with_max_iter.n_iter(), None);
assert_eq!(with_tol.n_iter(), None);
Ok(())
}
// -- positive=True non-negative ridge (REQ-9) --------------------------
#[test]
fn ridge_positive_matches_sklearn() -> Result<(), FerroError> {
// Live sklearn 1.5.2 oracle (non-negative ridge, alpha=1, positive=True):
// cd /tmp && python3 -c "import numpy as np; \
// from sklearn.linear_model import Ridge; \
// X=np.array([[1.,3.],[2.,1.],[3.,4.],[4.,2.],[5.,5.],[6.,1.],[2.,4.],[5.,2.]]); \
// y=1.0*X[:,0]-2.0*X[:,1]+np.array([0.1,-0.2,0.15,0.0,-0.1,0.05,0.2,-0.05]); \
// m=Ridge(alpha=1.0,positive=True).fit(X,y); \
// print([round(c,8) for c in m.coef_], round(m.intercept_,8))"
// -> [1.19891304, 0.0] -6.17744565
let x = Array2::from_shape_vec(
(8, 2),
vec![
1., 3., 2., 1., 3., 4., 4., 2., 5., 5., 6., 1., 2., 4., 5., 2.,
],
)
.map_err(|e| FerroError::ShapeMismatch {
expected: vec![8, 2],
actual: vec![],
context: e.to_string(),
})?;
let raw = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 2.0, 5.0];
let f1 = [3.0, 1.0, 4.0, 2.0, 5.0, 1.0, 4.0, 2.0];
let noise = [0.1, -0.2, 0.15, 0.0, -0.1, 0.05, 0.2, -0.05];
let y: Array1<f64> = (0..8).map(|i| raw[i] - 2.0 * f1[i] + noise[i]).collect();
let model = Ridge::<f64>::new().with_alpha(1.0).with_positive(true);
let fitted = model.fit(&x, &y)?;
assert_relative_eq!(fitted.coefficients()[0], 1.198_913_04, epsilon = 1e-5);
assert_relative_eq!(fitted.coefficients()[1], 0.0, epsilon = 1e-5);
assert_relative_eq!(fitted.intercept(), -6.177_445_65, epsilon = 1e-4);
// All coefficients are non-negative.
for &c in fitted.coefficients().iter() {
assert!(
c >= 0.0,
"coef {c} must be non-negative under positive=True"
);
}
// Non-tautological: the constrained fit MUST differ from the
// unconstrained ridge (oracle unconstrained coef_
// [0.95708502, -1.85401484], intercept -0.23250675 — feature 1
// is strongly negative and clamps to 0).
let unconstrained = Ridge::<f64>::new().with_alpha(1.0).fit(&x, &y)?;
assert_relative_eq!(
unconstrained.coefficients()[0],
0.957_085_02,
epsilon = 1e-5
);
assert_relative_eq!(
unconstrained.coefficients()[1],
-1.854_014_84,
epsilon = 1e-5
);
assert!(
(fitted.coefficients()[1] - unconstrained.coefficients()[1]).abs() > 1e-2,
"positive fit must differ from unconstrained on the clamped coordinate"
);
Ok(())
}
#[test]
fn ridge_positive_false_unchanged() -> Result<(), FerroError> {
// Regression guard: positive=false (default) is BYTE-IDENTICAL to the
// current unconstrained Cholesky `fit`.
let x = Array2::from_shape_vec((5, 2), vec![1., 2., 2., 1., 3., 4., 4., 3., 5., 5.])
.map_err(|e| FerroError::ShapeMismatch {
expected: vec![5, 2],
actual: vec![],
context: e.to_string(),
})?;
let y = array![3.0, 2.5, 7.1, 6.0, 11.2];
let model = Ridge::<f64>::new().with_alpha(1.0);
let baseline = model.fit(&x, &y)?;
let explicit_false = Ridge::<f64>::new()
.with_alpha(1.0)
.with_positive(false)
.fit(&x, &y)?;
for j in 0..2 {
assert_eq!(
baseline.coefficients()[j].to_bits(),
explicit_false.coefficients()[j].to_bits(),
"coef_[{j}] must be byte-identical with positive=false"
);
}
assert_eq!(
baseline.intercept().to_bits(),
explicit_false.intercept().to_bits()
);
// positive=false keeps the direct-solver n_iter_ = None contract.
assert_eq!(explicit_false.n_iter(), None);
Ok(())
}
#[test]
fn ridge_positive_all_nonneg_equals_unconstrained() -> Result<(), FerroError> {
// Sanity: data whose unconstrained ridge is already non-negative —
// the box constraint is inactive, so positive=True must reproduce the
// unconstrained coefficients.
//
// Live sklearn 1.5.2 oracle:
// cd /tmp && python3 -c "import numpy as np; \
// from sklearn.linear_model import Ridge; \
// X=np.array([[1.,2.],[2.,1.],[3.,4.],[4.,3.],[5.,5.]]); \
// y=np.array([3.0,2.5,7.1,6.0,11.2]); \
// u=Ridge(alpha=1.0).fit(X,y); p=Ridge(alpha=1.0,positive=True).fit(X,y); \
// print([round(c,8) for c in u.coef_], [round(c,8) for c in p.coef_])"
// -> [0.82280702, 1.35614035] [0.82280702, 1.35614035]
let x = Array2::from_shape_vec((5, 2), vec![1., 2., 2., 1., 3., 4., 4., 3., 5., 5.])
.map_err(|e| FerroError::ShapeMismatch {
expected: vec![5, 2],
actual: vec![],
context: e.to_string(),
})?;
let y = array![3.0, 2.5, 7.1, 6.0, 11.2];
let unconstrained = Ridge::<f64>::new().with_alpha(1.0).fit(&x, &y)?;
let positive = Ridge::<f64>::new()
.with_alpha(1.0)
.with_positive(true)
.fit(&x, &y)?;
for j in 0..2 {
assert!(unconstrained.coefficients()[j] >= 0.0);
// CD converges to the unconstrained optimum at the default
// tol=1e-4; sklearn's own L-BFGS-B agrees only to ~1e-6 on this
// fixture (oracle: 0.82280707 vs 0.82280702).
assert_relative_eq!(
positive.coefficients()[j],
unconstrained.coefficients()[j],
epsilon = 1e-4
);
}
assert_relative_eq!(
positive.intercept(),
unconstrained.intercept(),
epsilon = 1e-4
);
Ok(())
}
// -- per-target alpha array (REQ-7) ------------------------------------
/// Build the shared 5x2 fixture `X` and 5x2 `Y` used by the per-target
/// alpha tests.
fn per_target_fixture() -> Result<(Array2<f64>, Array2<f64>), FerroError> {
let shape_err = |e: ndarray::ShapeError| FerroError::ShapeMismatch {
expected: vec![5, 2],
actual: vec![],
context: e.to_string(),
};
let x = Array2::from_shape_vec((5, 2), vec![1., 2., 2., 1., 3., 4., 4., 3., 5., 5.])
.map_err(shape_err)?;
let y = Array2::from_shape_vec(
(5, 2),
vec![3.0, 1.0, 2.5, 2.0, 7.1, 3.5, 6.0, 4.2, 11.2, 6.0],
)
.map_err(shape_err)?;
Ok((x, y))
}
#[test]
fn ridge_per_target_alpha_matches_sklearn() -> Result<(), FerroError> {
// Live sklearn 1.5.2 oracle (per-target alpha [0.5, 2.0]):
// cd /tmp && python3 -c "import numpy as np; \
// from sklearn.linear_model import Ridge; \
// X=np.array([[1.,2.],[2.,1.],[3.,4.],[4.,3.],[5.,5.]]); \
// Y=np.array([[3.,1.],[2.5,2.],[7.1,3.5],[6.,4.2],[11.2,6.]]); \
// m=Ridge(alpha=np.array([0.5,2.0])).fit(X,Y); \
// print(m.coef_.tolist(), m.intercept_.tolist())"
// -> coef_ (n_targets, n_features) = [[0.79891892, 1.43891892],
// [0.78, 0.355]]
// intercept_ = [-0.75351351, -0.065]
// ferrolearn stores coefficients as (n_features, n_targets) = the
// TRANSPOSE: column 0 = target 0, column 1 = target 1.
let (x, y) = per_target_fixture()?;
let model = Ridge::<f64>::new().with_alpha_per_target(array![0.5, 2.0]);
let fitted = model.fit(&x, &y)?;
let c0 = fitted.coefficients().column(0).to_owned();
let c1 = fitted.coefficients().column(1).to_owned();
assert_relative_eq!(c0[0], 0.798_918_92, epsilon = 1e-6);
assert_relative_eq!(c0[1], 1.438_918_92, epsilon = 1e-6);
assert_relative_eq!(c1[0], 0.78, epsilon = 1e-6);
assert_relative_eq!(c1[1], 0.355, epsilon = 1e-6);
assert_relative_eq!(fitted.intercepts()[0], -0.753_513_51, epsilon = 1e-6);
assert_relative_eq!(fitted.intercepts()[1], -0.065, epsilon = 1e-6);
Ok(())
}
#[test]
fn ridge_per_target_alpha_equals_independent_scalar_fits() -> Result<(), FerroError> {
// The per-target array path must reproduce, column by column, what an
// independent scalar-alpha Ridge produces on each target column via the
// single-output (1-D) `fit` path.
let (x, y) = per_target_fixture()?;
let alphas = array![0.5, 2.0];
let multi = Ridge::<f64>::new()
.with_alpha_per_target(alphas.clone())
.fit(&x, &y)?;
for k in 0..2 {
let y_col = y.column(k).to_owned();
let scalar = Ridge::<f64>::new().with_alpha(alphas[k]).fit(&x, &y_col)?;
for j in 0..2 {
assert_relative_eq!(
multi.coefficients()[[j, k]],
scalar.coefficients()[j],
epsilon = 1e-9
);
}
assert_relative_eq!(multi.intercepts()[k], scalar.intercept(), epsilon = 1e-9);
}
Ok(())
}
#[test]
fn ridge_per_target_alpha_length_mismatch_errors() -> Result<(), FerroError> {
// alpha array of length 1 with 2 target columns must error.
let (x, y) = per_target_fixture()?;
let model = Ridge::<f64>::new().with_alpha_per_target(array![0.5]);
let result = <Ridge<f64> as Fit<Array2<f64>, Array2<f64>>>::fit(&model, &x, &y);
assert!(matches!(result, Err(FerroError::ShapeMismatch { .. })));
Ok(())
}
#[test]
fn ridge_multi_scalar_alpha_unchanged() -> Result<(), FerroError> {
// Regression guard: with NO per-target array, the multi-output scalar
// path is byte-identical to the historic `solve_ridge_multi` behavior.
//
// Live sklearn 1.5.2 oracle (scalar alpha=1.0):
// cd /tmp && python3 -c "import numpy as np; \
// from sklearn.linear_model import Ridge; \
// X=np.array([[1.,2.],[2.,1.],[3.,4.],[4.,3.],[5.,5.]]); \
// Y=np.array([[3.,1.],[2.5,2.],[7.1,3.5],[6.,4.2],[11.2,6.]]); \
// m=Ridge(alpha=1.0).fit(X,Y); \
// print(m.coef_.tolist(), m.intercept_.tolist())"
let (x, y) = per_target_fixture()?;
let baseline = Ridge::<f64>::new().with_alpha(1.0).fit(&x, &y)?;
let explicit_none = Ridge::<f64> {
alpha_per_target: None,
..Ridge::<f64>::new().with_alpha(1.0)
}
.fit(&x, &y)?;
for j in 0..2 {
for k in 0..2 {
assert_eq!(
baseline.coefficients()[[j, k]].to_bits(),
explicit_none.coefficients()[[j, k]].to_bits(),
"coef_[{j},{k}] must be byte-identical without per-target alpha"
);
}
}
for k in 0..2 {
assert_eq!(
baseline.intercepts()[k].to_bits(),
explicit_none.intercepts()[k].to_bits(),
"intercept_[{k}] must be byte-identical without per-target alpha"
);
}
Ok(())
}
// -- solver variants (auto/cholesky/svd) + solver_ (REQ-8a) ------------
/// Build the shared 5x2 fixture `X` and 1-D `y` used by the solver tests.
fn solver_fixture() -> Result<(Array2<f64>, Array1<f64>), FerroError> {
let x = Array2::from_shape_vec((5, 2), vec![1., 2., 2., 1., 3., 4., 4., 3., 5., 5.])
.map_err(|e| FerroError::ShapeMismatch {
expected: vec![5, 2],
actual: vec![],
context: e.to_string(),
})?;
let y = array![3.0, 2.5, 7.1, 6.0, 11.2];
Ok((x, y))
}
#[test]
fn ridge_solver_svd_matches_sklearn_and_cholesky() -> Result<(), FerroError> {
// Live sklearn 1.5.2 oracle (alpha=1, fit_intercept=True, solver='svd'):
// cd /tmp && python3 -c "import numpy as np; \
// from sklearn.linear_model import Ridge; \
// X=np.array([[1.,2.],[2.,1.],[3.,4.],[4.,3.],[5.,5.]]); \
// y=np.array([3.0,2.5,7.1,6.0,11.2]); \
// m=Ridge(alpha=1.0,solver='svd').fit(X,y); \
// print([round(c,10) for c in m.coef_], round(m.intercept_,10))"
// -> [0.8228070175, 1.3561403509] -0.5768421053
// Every dense solver returns this same unique ridge solution.
let (x, y) = solver_fixture()?;
let svd_fit = Ridge::<f64>::new()
.with_alpha(1.0)
.with_solver(RidgeSolver::Svd)
.fit(&x, &y)?;
assert_relative_eq!(svd_fit.coefficients()[0], 0.822_807_017_5, epsilon = 1e-7);
assert_relative_eq!(svd_fit.coefficients()[1], 1.356_140_350_9, epsilon = 1e-7);
assert_relative_eq!(svd_fit.intercept(), -0.576_842_105_3, epsilon = 1e-7);
// SVD and Cholesky solve the same strictly-convex (unique) problem, so
// they must agree to ~1e-9.
let chol_fit = Ridge::<f64>::new()
.with_alpha(1.0)
.with_solver(RidgeSolver::Cholesky)
.fit(&x, &y)?;
for j in 0..2 {
assert_relative_eq!(
svd_fit.coefficients()[j],
chol_fit.coefficients()[j],
epsilon = 1e-9
);
}
assert_relative_eq!(svd_fit.intercept(), chol_fit.intercept(), epsilon = 1e-9);
Ok(())
}
#[test]
fn ridge_solver_resolution() -> Result<(), FerroError> {
// Default `solver` is Auto; the fitted `solver_` resolves Auto→Cholesky
// for the dense path (sklearn `resolve_solver`, `_ridge.py:830`).
assert_eq!(Ridge::<f64>::new().solver, RidgeSolver::Auto);
let (x, y) = solver_fixture()?;
let auto_fit = Ridge::<f64>::new().with_alpha(1.0).fit(&x, &y)?;
assert_eq!(
auto_fit.solver(),
RidgeSolver::Cholesky,
"Auto must resolve to Cholesky for the dense path"
);
let svd_fit = Ridge::<f64>::new()
.with_alpha(1.0)
.with_solver(RidgeSolver::Svd)
.fit(&x, &y)?;
assert_eq!(svd_fit.solver(), RidgeSolver::Svd);
let chol_fit = Ridge::<f64>::new()
.with_alpha(1.0)
.with_solver(RidgeSolver::Cholesky)
.fit(&x, &y)?;
assert_eq!(chol_fit.solver(), RidgeSolver::Cholesky);
Ok(())
}
#[test]
fn ridge_solver_default_cholesky_unchanged() -> Result<(), FerroError> {
// Regression guard: the default (Auto→Cholesky) coef_/intercept_ are
// BYTE-IDENTICAL to a pre-existing direct `fit` (no `with_solver`).
let (x, y) = solver_fixture()?;
let direct = Ridge::<f64>::new().with_alpha(1.0).fit(&x, &y)?;
let explicit_chol = Ridge::<f64>::new()
.with_alpha(1.0)
.with_solver(RidgeSolver::Cholesky)
.fit(&x, &y)?;
let explicit_auto = Ridge::<f64>::new()
.with_alpha(1.0)
.with_solver(RidgeSolver::Auto)
.fit(&x, &y)?;
for j in 0..2 {
assert_eq!(
direct.coefficients()[j].to_bits(),
explicit_chol.coefficients()[j].to_bits(),
"coef_[{j}] must be byte-identical for explicit Cholesky"
);
assert_eq!(
direct.coefficients()[j].to_bits(),
explicit_auto.coefficients()[j].to_bits(),
"coef_[{j}] must be byte-identical for explicit Auto"
);
}
assert_eq!(
direct.intercept().to_bits(),
explicit_chol.intercept().to_bits()
);
assert_eq!(
direct.intercept().to_bits(),
explicit_auto.intercept().to_bits()
);
Ok(())
}
#[test]
fn ridge_solver_svd_no_intercept() -> Result<(), FerroError> {
// With fit_intercept=false, the SVD solver must match the Cholesky
// no-intercept fit (same unique solution) to ~1e-9. sklearn:
// Ridge(alpha=1.0, fit_intercept=False, solver='svd') ==
// Ridge(alpha=1.0, fit_intercept=False, solver='cholesky').
let (x, y) = solver_fixture()?;
let svd_fit = Ridge::<f64>::new()
.with_alpha(1.0)
.with_fit_intercept(false)
.with_solver(RidgeSolver::Svd)
.fit(&x, &y)?;
let chol_fit = Ridge::<f64>::new()
.with_alpha(1.0)
.with_fit_intercept(false)
.with_solver(RidgeSolver::Cholesky)
.fit(&x, &y)?;
for j in 0..2 {
assert_relative_eq!(
svd_fit.coefficients()[j],
chol_fit.coefficients()[j],
epsilon = 1e-9
);
}
assert_eq!(svd_fit.intercept(), 0.0);
assert_eq!(chol_fit.intercept(), 0.0);
assert_eq!(svd_fit.solver(), RidgeSolver::Svd);
Ok(())
}
}