Skip to content

feat(linear): add native multi-output matrix support to LinearRegression #433

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mi7plus wants to merge 3 commits into
base: development
Choose a base branch
from
Open

feat(linear): add native multi-output matrix support to LinearRegression #433

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d3402b4
feat(linear): add native multi-output matrix support to LinearRegression
mi7plus Aug 20, 2026
a41057e
style: apply rustfmt formatting to linear regression tests
mi7plus Aug 20, 2026
2635d3c
fix(linear): resolve reviewer feedback for multi-output regression
mi7plus Aug 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
262 changes: 218 additions & 44 deletions src/linear/linear_regression.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,13 +16,14 @@
//! The QR decomposition is more computationally efficient and more numerically stable than calculating the normal equation directly,
//! but does not work for all data matrices. Unlike the QR decomposition, all matrices have an SVD decomposition.
//!
//! Example:
//! ### Examples
//!
//! #### Single-Output Regression
//! ```
//! use smartcore::linalg::basic::matrix::DenseMatrix;
//! use smartcore::linear::linear_regression::*;
//!
//! // Longley dataset (https://www.statsmodels.org/stable/datasets/generated/longley.html)
//! // Longley dataset ([https://www.statsmodels.org/stable/datasets/generated/longley.html](https://www.statsmodels.org/stable/datasets/generated/longley.html))
//! let x = DenseMatrix::from_2d_array(&[
//! &[234.289, 235.6, 159.0, 107.608, 1947., 60.323],
//! &[259.426, 232.5, 145.6, 108.632, 1948., 61.122],
Expand Down Expand Up @@ -52,6 +53,23 @@
//! let y_hat = lr.predict(&x).unwrap();
//! ```
//!
//! #### Multi-Output Regression
//! ```
//! use smartcore::linalg::basic::matrix::DenseMatrix;
//! use smartcore::linear::linear_regression::*;
//!
//! let x = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[2.0, 1.0], &[3.0, 4.0], &[4.0, 3.0]]).unwrap();
//! let y = DenseMatrix::from_2d_array(&[&[4.0, 1.5], &[5.0, 3.0], &[10.0, 0.5], &[11.0, 2.0]]).unwrap();
//!
//! let model = LinearRegression::<f64, f64, DenseMatrix<f64>, Vec<f64>>::fit_matrix(
//! &x,
//! &y,
//! Default::default(),
//! ).unwrap();
//!
//! let y_hat = model.predict_matrix(&x).unwrap();
//! ```
//!
//! ## References:
//!
//! * ["Pattern Recognition and Machine Learning", C.M. Bishop, Linear Models for Regression](https://www.microsoft.com/en-us/research/uploads/prod/2006/01/Bishop-Pattern-Recognition-and-Machine-Learning-2006.pdf)
Expand All @@ -60,6 +78,7 @@
//!
//! <script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
//! <script id="MathJax-script" async src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>

use std::fmt::Debug;
use std::marker::PhantomData;

Expand Down Expand Up @@ -109,10 +128,11 @@ pub struct LinearRegression<
TX: Number + RealNumber,
TY: Number,
X: Array2<TX> + QRDecomposable<TX> + SVDDecomposable<TX>,
Y: Array1<TY>,
Y,
> {
coefficients: Option<X>,
intercept: Option<TX>,
intercept: Option<X>,
_phantom_tx: PhantomData<TX>,
_phantom_ty: PhantomData<TY>,
_phantom_y: PhantomData<Y>,
}
Expand Down Expand Up @@ -184,17 +204,33 @@ impl<
TX: Number + RealNumber,
TY: Number,
X: Array2<TX> + QRDecomposable<TX> + SVDDecomposable<TX>,
Y: Array1<TY>,
Y,
> PartialEq for LinearRegression<TX, TY, X, Y>
{
fn eq(&self, other: &Self) -> bool {
self.intercept == other.intercept
&& self.coefficients().shape() == other.coefficients().shape()
&& self
.coefficients()
.iterator(0)
.zip(other.coefficients().iterator(0))
.all(|(&a, &b)| (a - b).abs() <= TX::epsilon())
let intercepts_match = match (&self.intercept, &other.intercept) {
(Some(a), Some(b)) => {
a.shape() == b.shape()
&& a.iterator(0)
.zip(b.iterator(0))
.all(|(&x, &y)| (x - y).abs() <= TX::epsilon())
}
(None, None) => true,
_ => false,
};

let coefficients_match = match (&self.coefficients, &other.coefficients) {
(Some(a), Some(b)) => {
a.shape() == b.shape()
&& a.iterator(0)
.zip(b.iterator(0))
.all(|(&x, &y)| (x - y).abs() <= TX::epsilon())
}
(None, None) => true,
_ => false,
};

intercepts_match && coefficients_match
}
}

Expand All @@ -209,6 +245,7 @@ impl<
Self {
coefficients: Option::None,
intercept: Option::None,
_phantom_tx: PhantomData,
_phantom_ty: PhantomData,
_phantom_y: PhantomData,
}
Expand All @@ -235,77 +272,122 @@ impl<
TX: Number + RealNumber,
TY: Number,
X: Array2<TX> + QRDecomposable<TX> + SVDDecomposable<TX>,
Y: Array1<TY>,
Y,
> LinearRegression<TX, TY, X, Y>
{
/// Fits Linear Regression to your data.
/// * `x` - _NxM_ matrix with _N_ observations and _M_ features in each observation.
/// * `y` - target values
/// * `parameters` - other parameters, use `Default::default()` to set parameters to default values.
pub fn fit(
/// Fits Linear Regression when target Y is a matrix (N x K).
pub fn fit_matrix<YM: Array2<TX>>(
x: &X,
y: &Y,
y: &YM,
parameters: LinearRegressionParameters,
) -> Result<LinearRegression<TX, TY, X, Y>, Failed> {
let b = X::from_iterator(
y.iterator(0).map(|&v| TX::from(v).unwrap()),
y.shape(),
1,
0,
);
let (x_nrows, num_attributes) = x.shape();
let (y_nrows, _) = b.shape();
let (y_nrows, num_targets) = y.shape();

if x_nrows != y_nrows {
return Err(Failed::fit(
"Number of rows of X doesn\'t match number of rows of Y",
"Number of rows of X doesn't match number of rows of Y",
));
}

// Augment X with column of ones for intercept fitting
let a = x.h_stack(&X::ones(x_nrows, 1));

// Convert YM into X safely via iterator mapping
let y_x = X::from_iterator(y.iterator(0).copied(), y_nrows, num_targets, 0);

// Solve A * W = Y for W matrix of size (num_attributes + 1) x K
let w = match parameters.solver {
LinearRegressionSolverName::QR => a.qr_solve_mut(b)?,
LinearRegressionSolverName::SVD => a.svd_solve_mut(b)?,
LinearRegressionSolverName::QR => a.qr_solve_mut(y_x)?,
LinearRegressionSolverName::SVD => a.svd_solve_mut(y_x)?,
};

let weights = X::from_slice(w.slice(0..num_attributes, 0..1).as_ref());
let weights = X::from_slice(w.slice(0..num_attributes, 0..num_targets).as_ref());
let intercept = X::from_slice(
w.slice(num_attributes..num_attributes + 1, 0..num_targets)
.as_ref(),
);

Ok(LinearRegression {
intercept: Some(*w.get((num_attributes, 0))),
intercept: Some(intercept),
coefficients: Some(weights),
_phantom_tx: PhantomData,
_phantom_ty: PhantomData,
_phantom_y: PhantomData,
})
}

/// Predict target values from `x`
/// * `x` - _KxM_ data where _K_ is number of observations and _M_ is number of features.
pub fn predict(&self, x: &X) -> Result<Y, Failed> {
/// Predict matrix Y_hat given input feature matrix X
pub fn predict_matrix(&self, x: &X) -> Result<X, Failed> {
let (nrows, _) = x.shape();
let bias = X::fill(nrows, 1, *self.intercept());
let mut y_hat = x.matmul(self.coefficients());
y_hat.add_mut(&bias);
Ok(Y::from_iterator(
y_hat.iterator(0).map(|&v| TY::from(v).unwrap()),
nrows,
))

let intercept = self.intercept_matrix();
let (_, num_targets) = intercept.shape();

for i in 0..nrows {
for j in 0..num_targets {
let current = y_hat.get((i, j));
let bias = intercept.get((0, j));
y_hat.set((i, j), *current + *bias);
}
}
Ok(y_hat)
}

/// Get estimates regression coefficients
/// Get estimated regression coefficients
pub fn coefficients(&self) -> &X {
self.coefficients.as_ref().unwrap()
}

/// Get estimate of intercept
/// Get primary intercept value (for single-output regression callers)
pub fn intercept(&self) -> &TX {
self.intercept.as_ref().unwrap().get((0, 0))
}

/// Get complete intercept matrix (1 x K) for multi-output models
pub fn intercept_matrix(&self) -> &X {
self.intercept.as_ref().unwrap()
}
}

impl<
TX: Number + RealNumber,
TY: Number,
X: Array2<TX> + QRDecomposable<TX> + SVDDecomposable<TX>,
Y: Array1<TY>,
> LinearRegression<TX, TY, X, Y>
{
/// Fits Linear Regression to 1D target values.
pub fn fit(
x: &X,
y: &Y,
parameters: LinearRegressionParameters,
) -> Result<LinearRegression<TX, TY, X, Y>, Failed> {
let y_mat = X::from_iterator(
y.iterator(0).map(|&v| TX::from(v).unwrap()),
y.shape(),
1,
0,
);
Self::fit_matrix(x, &y_mat, parameters)
}

/// Predict 1D target values from `x`
pub fn predict(&self, x: &X) -> Result<Y, Failed> {
let y_hat_mat = self.predict_matrix(x)?;
let (nrows, _) = x.shape();
Ok(Y::from_iterator(
y_hat_mat.iterator(0).map(|&v| TY::from(v).unwrap()),
nrows,
))
}
}

#[cfg(test)]
mod tests {
use super::*;
use crate::linalg::basic::arrays::Array;
use crate::linalg::basic::matrix::DenseMatrix;

#[test]
Expand Down Expand Up @@ -348,16 +430,17 @@ mod tests {
83.0, 88.5, 88.2, 89.5, 96.2, 98.1, 99.0, 100.0, 101.2, 104.6, 108.4, 110.8,
];

let y_hat_qr = LinearRegression::fit(
let lr = LinearRegression::fit(
&x,
&y,
LinearRegressionParameters {
solver: LinearRegressionSolverName::QR,
},
)
.and_then(|lr| lr.predict(&x))
.unwrap();

let y_hat_qr = lr.predict(&x).unwrap();

let y_hat_svd = LinearRegression::fit(&x, &y, Default::default())
.and_then(|lr| lr.predict(&x))
.unwrap();
Expand All @@ -372,6 +455,10 @@ mod tests {
.zip(y_hat_svd.iter())
.all(|(&a, &b)| (a - b).abs() <= 5.0)
);

// Confirm single-output intercept getter works properly
let intercept: f64 = *lr.intercept();
assert!((intercept - 83.0).abs() > 0.0);
}

#[cfg_attr(
Expand Down Expand Up @@ -413,4 +500,91 @@ mod tests {

assert_eq!(lr, deserialized_lr);
}

#[test]
fn multi_output_ols_fit_predict() {
let x = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[2.0, 1.0], &[3.0, 4.0], &[4.0, 3.0]])
.unwrap();

let y = DenseMatrix::from_2d_array(&[&[4.0, 1.5], &[5.0, 3.0], &[10.0, 0.5], &[11.0, 2.0]])
.unwrap();

// Test SVD solver path
let model_svd = LinearRegression::<f64, f64, DenseMatrix<f64>, Vec<f64>>::fit_matrix(
&x,
&y,
LinearRegressionParameters::default().with_solver(LinearRegressionSolverName::SVD),
)
.unwrap();
let pred_svd = model_svd.predict_matrix(&x).unwrap();

// Test QR solver path
let model_qr = LinearRegression::<f64, f64, DenseMatrix<f64>, Vec<f64>>::fit_matrix(
&x,
&y,
LinearRegressionParameters::default().with_solver(LinearRegressionSolverName::QR),
)
.unwrap();
let pred_qr = model_qr.predict_matrix(&x).unwrap();

assert_eq!(pred_svd.shape(), (4, 2));
assert_eq!(pred_qr.shape(), (4, 2));
assert_eq!(model_svd.coefficients().shape(), (2, 2));
}

#[test]
fn multi_output_numerical_correctness() {
// Target model: Y_1 = 2*X_1 + 1*X_2 + 3, Y_2 = 1*X_1 + 3*X_2 + 2
let x = DenseMatrix::from_2d_array(&[
&[1.0, 0.0],
&[0.0, 1.0],
&[2.0, 2.0],
&[3.0, 1.0],
&[1.0, 4.0],
])
.unwrap();

let y = DenseMatrix::from_2d_array(&[
&[5.0, 3.0],
&[4.0, 5.0],
&[9.0, 10.0],
&[10.0, 8.0],
&[9.0, 15.0],
])
.unwrap();

let model = LinearRegression::<f64, f64, DenseMatrix<f64>, Vec<f64>>::fit_matrix(
&x,
&y,
Default::default(),
)
.unwrap();

let coefs = model.coefficients();
let intercept = model.intercept_matrix();

// Check slope matrix (2 x 2)
assert!((coefs.get((0, 0)) - 2.0).abs() < 1e-5);
assert!((coefs.get((1, 0)) - 1.0).abs() < 1e-5);
assert!((coefs.get((0, 1)) - 1.0).abs() < 1e-5);
assert!((coefs.get((1, 1)) - 3.0).abs() < 1e-5);

// Check intercept vector (1 x 2)
assert!((intercept.get((0, 0)) - 3.0).abs() < 1e-5);
assert!((intercept.get((0, 1)) - 2.0).abs() < 1e-5);
}

#[test]
fn fit_matrix_dimension_mismatch() {
let x = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[2.0, 1.0]]).unwrap();
let y = DenseMatrix::from_2d_array(&[&[4.0, 1.5]]).unwrap();

let result = LinearRegression::<f64, f64, DenseMatrix<f64>, Vec<f64>>::fit_matrix(
&x,
&y,
Default::default(),
);

assert!(result.is_err());
}
}
Loading