From d3dae39f35883b4273effd8ee808424132232b16 Mon Sep 17 00:00:00 2001 From: Johan Larsson Date: Thu, 8 Oct 2026 22:43:19 +0200 Subject: [PATCH] feat: add lazymatrix feature --- .github/workflows/ci.yml | 4 +- CHANGELOG.md | 5 + Cargo.toml | 4 + README.md | 23 + src/linear/elastic_net.rs | 323 +++++++++++- src/linear/lasso.rs | 290 ++++++++++- src/linear/lasso_optimizer.rs | 75 ++- src/linear/lazy_normalization.rs | 830 +++++++++++++++++++++++++++++++ src/linear/mod.rs | 2 + 9 files changed, 1527 insertions(+), 29 deletions(-) create mode 100644 src/linear/lazy_normalization.rs diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7715aed7..01bcc12f 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -57,6 +57,8 @@ jobs: - "--features serde" - "--features datasets" - "--features ndarray-bindings" + - "--features lazy-normalization" + - "--features lazy-normalization,ndarray-bindings" - "" steps: - uses: actions/checkout@v4 @@ -73,7 +75,7 @@ jobs: - name: Stable Build run: cargo build --no-default-features ${{ matrix.features }} - name: Tests - if: matrix.features == '--features ndarray-bindings' + if: matrix.features == '--features ndarray-bindings' || contains(matrix.features, 'lazy-normalization') run: cargo test --no-default-features ${{ matrix.features }} msrv: diff --git a/CHANGELOG.md b/CHANGELOG.md index 2c0ff5c5..e32e47a4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Added +- Optional `lazy-normalization` support for Lasso and Elastic Net. It applies normalization to a borrowed design and keeps Elastic Net's identity-block augmentation implicit, while retaining the dense Gram workspace. Enable the feature and set `with_lazy_normalization(true)` to select this path. The default remains eager. The new parameter field is available only with the feature and defaults to false during deserialization (#479). + +### Changed +- **Breaking change:** Enabling `lazy-normalization` adds a public `lazy_normalization` field to `LassoParameters`, `LassoSearchParameters`, `ElasticNetParameters`, and `ElasticNetSearchParameters`. Struct literals without `..Default::default()` must add `lazy_normalization: false` to retain the eager path, or `true` to select the lazy path. This also applies when another dependency enables the feature through Cargo feature unification (#479). ## [0.6.16] ### Fixed diff --git a/Cargo.toml b/Cargo.toml index 4ec3a4b5..3e8ab12e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -24,6 +24,7 @@ exclude = [ [dependencies] approx = "0.5.1" cfg-if = "1.0.0" +lazymatrix = { version = "0.6.0", optional = true, default-features = false } ndarray = { version = "0.17", optional = true } num-traits = "0.2.12" num = "0.4" @@ -46,6 +47,9 @@ serde = ["dep:serde", "dep:typetag"] # Optional bindings for the ndarray crate (DenseMatrix <-> ndarray). ndarray-bindings = ["dep:ndarray"] +# Apply Lasso and Elastic Net normalization without copying the design matrix. +lazy-normalization = ["dep:lazymatrix"] + # Sample datasets (iris, boston, digits, ...) and dataset generators. # Enables std_rand and serde because datasets ship serialized bundles. datasets = ["std_rand", "serde"] diff --git a/README.md b/README.md index b219b2f9..dac4ac5c 100644 --- a/README.md +++ b/README.md @@ -100,6 +100,29 @@ smartcore organizes algorithms into clear modules with consistent traits: Recent refactors emphasize reusable components in trees/forests and expanded multiclass SVM capabilities. XGBoost-style regression and single-linkage clustering have been added. See CHANGELOG for API changes and migration notes. +## Optional lazy normalization + +The optional `lazy-normalization` feature lets Lasso and Elastic Net borrow the +design matrix during normalization. Elastic Net also applies its identity-block +augmentation without storing the augmented design. Enable it with `features = +["lazy-normalization"]`, then select it with +`LassoParameters::default().with_lazy_normalization(true)` or +`ElasticNetParameters::default().with_lazy_normalization(true)`. The default +remains eager, even when the feature is enabled. The `normalize` option still +controls whether columns are centered and scaled. Parameter searches use the +fixed `lazy_normalization` setting. It works with `DenseMatrix` and +`ndarray-bindings`. + +Lazy normalization reduces memory use by computing normalized values when the +solver needs them. This adds repeated arithmetic and can increase fit time for +both estimators, including Elastic Net on large inputs. The solver still stores +a dense Gram matrix with one row and column per feature. + +Enabling the feature adds a public `lazy_normalization` field to both estimators' +parameter and search parameter structs. Code that constructs these structs +without `..Default::default()` must add the field. Missing fields default to +`false` during deserialization. + ## Data access and readers - CSV readers: Read matrices from CSV with configurable delimiter and header rows, with helpful error messages and testing utilities (including non-IO reader abstractions). diff --git a/src/linear/elastic_net.rs b/src/linear/elastic_net.rs index 30544bcd..f1045420 100644 --- a/src/linear/elastic_net.rs +++ b/src/linear/elastic_net.rs @@ -68,12 +68,20 @@ use crate::numbers::floatnum::FloatNumber; use crate::numbers::realnum::RealNumber; use crate::linear::lasso_optimizer::InteriorPointOptimizer; +#[cfg(feature = "lazy-normalization")] +use crate::linear::lazy_normalization::LazyDesign; /// Elastic net parameters #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] #[derive(Debug, Clone)] #[must_use] pub struct ElasticNetParameters { + #[cfg(feature = "lazy-normalization")] + #[cfg_attr(feature = "serde", serde(default))] + /// Borrow the design matrix instead of storing its normalized values. + /// This reduces memory use but can increase fit time. The default is false. + /// This option does not change whether normalization is applied. + pub lazy_normalization: bool, #[cfg_attr(feature = "serde", serde(default))] /// Regularization parameter. pub alpha: f64, @@ -104,6 +112,24 @@ pub struct ElasticNet, Y } impl ElasticNetParameters { + /// Select lazy matrix operations. Requires the `lazy-normalization` feature. + /// This reduces memory use but can increase fit time. + /// + /// ``` + /// # #[cfg(feature = "lazy-normalization")] + /// # { + /// use smartcore::linear::elastic_net::ElasticNetParameters; + /// + /// let parameters = ElasticNetParameters::default().with_lazy_normalization(true); + /// assert!(parameters.lazy_normalization); + /// # } + /// ``` + #[cfg(feature = "lazy-normalization")] + pub fn with_lazy_normalization(mut self, lazy_normalization: bool) -> Self { + self.lazy_normalization = lazy_normalization; + self + } + /// Regularization parameter. pub fn with_alpha(mut self, alpha: f64) -> Self { self.alpha = alpha; @@ -136,6 +162,8 @@ impl ElasticNetParameters { impl Default for ElasticNetParameters { fn default() -> Self { ElasticNetParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: false, alpha: 1.0, l1_ratio: 0.5, normalize: true, @@ -150,6 +178,13 @@ impl Default for ElasticNetParameters { #[derive(Debug, Clone)] #[must_use] pub struct ElasticNetSearchParameters { + #[cfg(feature = "lazy-normalization")] + #[cfg_attr(feature = "serde", serde(default))] + /// Borrow the design matrix instead of storing its normalized values. + /// This reduces memory use but can increase fit time. The default is false. + /// This setting applies to all parameter combinations in the search. + /// This option does not change whether normalization is applied. + pub lazy_normalization: bool, #[cfg_attr(feature = "serde", serde(default))] /// Regularization parameter. pub alpha: Vec, @@ -209,6 +244,8 @@ impl Iterator for ElasticNetSearchParametersIterator { } let next = ElasticNetParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: self.lasso_regression_search_parameters.lazy_normalization, alpha: self.lasso_regression_search_parameters.alpha[self.current_alpha], l1_ratio: self.lasso_regression_search_parameters.alpha[self.current_l1_ratio], normalize: self.lasso_regression_search_parameters.normalize[self.current_normalize], @@ -257,6 +294,8 @@ impl Default for ElasticNetSearchParameters { let default_params = ElasticNetParameters::default(); ElasticNetSearchParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: false, alpha: vec![default_params.alpha], l1_ratio: vec![default_params.l1_ratio], normalize: vec![default_params.normalize], @@ -320,7 +359,7 @@ impl, Y: Array1> y: &Y, parameters: ElasticNetParameters, ) -> Result, Failed> { - let (n, p) = x.shape(); + let (n, _) = x.shape(); if y.shape() != n { return Err(Failed::fit("Number of rows in X should = len(y)")); @@ -334,6 +373,32 @@ impl, Y: Array1> let y_mean = TX::from_f64(y.mean_by()).unwrap(); + #[cfg(feature = "lazy-normalization")] + let (w, b) = if parameters.lazy_normalization { + Self::fit_lazy_coefficients(x, y, ¶meters, l1_reg, l2_reg, y_mean)? + } else { + Self::fit_eager_coefficients(x, y, ¶meters, l1_reg, l2_reg, y_mean)? + }; + #[cfg(not(feature = "lazy-normalization"))] + let (w, b) = Self::fit_eager_coefficients(x, y, ¶meters, l1_reg, l2_reg, y_mean)?; + + Ok(ElasticNet { + intercept: Some(b), + coefficients: Some(w), + _phantom_ty: PhantomData, + _phantom_y: PhantomData, + }) + } + + fn fit_eager_coefficients( + x: &X, + y: &Y, + parameters: &ElasticNetParameters, + l1_reg: TX, + l2_reg: TX, + y_mean: TX, + ) -> Result<(X, TX), Failed> { + let (_, p) = x.shape(); let (w, b) = if parameters.normalize { let (scaled_x, col_mean, col_std) = Self::rescale_x(x)?; @@ -383,13 +448,73 @@ impl, Y: Array1> (X::from_column(&w), y_mean) }; + Ok((w, b)) + } - Ok(ElasticNet { - intercept: Some(b), - coefficients: Some(w), - _phantom_ty: PhantomData, - _phantom_y: PhantomData, - }) + #[cfg(feature = "lazy-normalization")] + fn fit_lazy_coefficients( + x: &X, + y: &Y, + parameters: &ElasticNetParameters, + l1_reg: TX, + l2_reg: TX, + y_mean: TX, + ) -> Result<(X, TX), Failed> { + let (n, p) = x.shape(); + let (w, b) = if parameters.normalize { + let x = { + let (col_mean, col_std) = Self::normalization_stats(x)?; + LazyDesign::new(x, Some(col_mean), Some(col_std)).with_l2(l2_reg) + }; + let (col_mean, col_std, gamma) = (x.centers(), x.scales(), x.gamma()); + let y = Self::augment_y(y, n + p); + + let mut optimizer = InteriorPointOptimizer::new_lazy(&x); + + let mut w = optimizer.optimize_lazy( + &x, + &y, + l1_reg * gamma, + parameters.max_iter, + TX::from_f64(parameters.tol).unwrap(), + true, + )?; + + for i in 0..p { + w.set(i, gamma * *w.get(i) / col_std[i]); + } + + let mut b = TX::zero(); + + for i in 0..p { + b += *w.get(i) * col_mean[i]; + } + + b = y_mean - b; + + (X::from_column(&w), b) + } else { + let x = LazyDesign::new(x, None, None).with_l2(l2_reg); + let (y, gamma) = (Self::augment_y(y, n + p), x.gamma()); + + let mut optimizer = InteriorPointOptimizer::new_lazy(&x); + + let mut w = optimizer.optimize_lazy( + &x, + &y, + l1_reg * gamma, + parameters.max_iter, + TX::from_f64(parameters.tol).unwrap(), + true, + )?; + + for i in 0..p { + w.set(i, gamma * *w.get(i)); + } + + (X::from_column(&w), y_mean) + }; + Ok((w, b)) } /// Predict target values from `x` @@ -422,7 +547,7 @@ impl, Y: Array1> self.intercept.as_ref().unwrap() } - fn rescale_x(x: &X) -> Result<(X, Vec, Vec), Failed> { + fn normalization_stats(x: &X) -> Result<(Vec, Vec), Failed> { let col_mean: Vec = x .mean_by(0) .iter() @@ -440,21 +565,31 @@ impl, Y: Array1> } } + Ok((col_mean, col_std)) + } + + fn rescale_x(x: &X) -> Result<(X, Vec, Vec), Failed> { + let (col_mean, col_std) = Self::normalization_stats(x)?; let mut scaled_x = x.clone(); scaled_x.scale_mut(&col_mean, &col_std, 0); Ok((scaled_x, col_mean, col_std)) } + fn augment_y(y: &Y, length: usize) -> Vec { + let mut result = Vec::::zeros(length); + for (dest, &value) in result.iter_mut().zip(y.iterator(0)) { + *dest = TX::from(value).unwrap(); + } + result + } + fn augment_x_and_y(x: &X, y: &Y, l2_reg: TX) -> (X, Vec, TX) { let (n, p) = x.shape(); let gamma = TX::one() / (TX::one() + l2_reg).sqrt(); let padding = gamma * l2_reg.sqrt(); - let mut y2 = Vec::::zeros(n + p); - for i in 0..y.shape() { - y2.set(i, TX::from(*y.get(i)).unwrap()); - } + let y2 = Self::augment_y(y, n + p); let mut x2 = X::zeros(n + p, p); @@ -476,9 +611,167 @@ mod tests { use crate::linalg::basic::matrix::DenseMatrix; use crate::metrics::mean_absolute_error; + #[cfg(feature = "lazy-normalization")] + fn compare_eager_fit>(tolerance: f64) { + for lazy_normalization in [false, true] { + use crate::linalg::basic::arrays::ArrayView1; + for (n, p) in [(8, 2), (3, 5)] { + let values: Vec = (0..n * p) + .map(|i| T::from_f64(((i * 7 + i / 3) % 13) as f64 / 3.0).unwrap()) + .collect(); + let x = X::from_iterator(values.iter().copied(), n, p, 0); + let y: Vec = (0..n) + .map(|i| { + T::from_f64(3.0).unwrap() + T::from_f64(2.0).unwrap() * values[i * p] + - values[i * p + 1] + }) + .collect(); + for normalize in [false, true] { + // The existing interior-point solver is unstable at a zero L1 penalty. + for ratio in [0.05, 0.5, 1.0] { + let parameters = ElasticNetParameters::default() + .with_alpha(0.04) + .with_lazy_normalization(lazy_normalization) + .with_l1_ratio(ratio) + .with_normalize(normalize); + let l1 = T::from_f64(0.04 * ratio * n as f64).unwrap(); + let l2 = T::from_f64(0.04 * (1.0 - ratio) * n as f64).unwrap(); + let (eager, centers, scales) = if normalize { + ElasticNet::>::rescale_x(&x).unwrap() + } else { + (x.clone(), vec![T::zero(); p], vec![T::one(); p]) + }; + let (eager, augmented_y, gamma) = + ElasticNet::>::augment_x_and_y(&eager, &y, l2); + let mut optimizer = InteriorPointOptimizer::new(&eager, p); + let expected = optimizer.optimize( + &eager, + &augmented_y, + l1 * gamma, + parameters.max_iter, + T::from_f64(parameters.tol).unwrap(), + true, + ); + let model = ElasticNet::fit(&x, &y, parameters); + let (mut expected, model) = match (expected, model) { + (Ok(expected), Ok(model)) => (expected, model), + (Err(expected), Err(actual)) => { + assert_eq!(actual.to_string(), expected.to_string()); + continue; + } + (expected, actual) => panic!( + "different fit outcomes: eager={expected:?}, lazy={actual:?}; n={n}, p={p}, ratio={ratio}, normalize={normalize}, scalar={}", + std::any::type_name::() + ), + }; + for (w, scale) in expected.iter_mut().zip(scales) { + *w = gamma * *w / scale; + } + for (&a, &b) in model.coefficients().iterator(0).zip(&expected) { + assert!( + (a - b).abs().to_f64().unwrap() < tolerance, + "{a} != {b}; n={n}, p={p}, ratio={ratio}" + ); + } + let intercept = T::from_f64(y.mean_by()).unwrap() + - expected.iter().zip(centers).map(|(&w, c)| w * c).sum::(); + assert!( + (*model.intercept() - intercept).abs().to_f64().unwrap() + < 10.0 * tolerance + ); + let prediction = model.predict(&x).unwrap(); + for (i, &actual) in prediction.iter().enumerate() { + let expected = + (0..p).map(|j| *x.get((i, j)) * expected[j]).sum::() + intercept; + assert!((actual - expected).abs().to_f64().unwrap() < 10.0 * tolerance); + } + assert!(x.iterator(0).zip(&values).all(|(a, b)| a == b)); + } + } + } + } + } + + #[cfg(feature = "lazy-normalization")] + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn lazy_fit_matches_eager_for_normalization_and_penalty_options() { + compare_eager_fit::>(1e-6); + compare_eager_fit::>(1e-2); + #[cfg(feature = "ndarray-bindings")] + { + compare_eager_fit::>(1e-6); + compare_eager_fit::>(1e-2); + } + } + + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn normalization_rejects_constant_and_near_constant_columns() { + for step in [0.0, f64::EPSILON / 16.0] { + let x = DenseMatrix::from_2d_array(&[&[0.0, 1.0], &[step, 2.0], &[2.0 * step, 4.0]]) + .unwrap(); + let error = ElasticNet::fit(&x, &vec![1., 2., 3.], ElasticNetParameters::default()) + .unwrap_err(); + assert!( + error + .to_string() + .contains("Cannot rescale constant column 0") + ); + #[cfg(feature = "lazy-normalization")] + { + let lazy_error = ElasticNet::fit( + &x, + &vec![1., 2., 3.], + ElasticNetParameters::default().with_lazy_normalization(true), + ) + .unwrap_err(); + assert_eq!(lazy_error, error); + } + } + } + + #[cfg(feature = "lazy-normalization")] + #[test] + fn lazy_selection_defaults_and_search() { + assert!(!ElasticNetParameters::default().lazy_normalization); + assert!(!ElasticNetSearchParameters::default().lazy_normalization); + for selected in [false, true] { + let search = ElasticNetSearchParameters { + lazy_normalization: selected, + alpha: vec![0.1, 0.2], + ..Default::default() + }; + let parameters: Vec<_> = search.into_iter().collect(); + assert_eq!(parameters.len(), 2); + assert!(parameters.iter().all(|p| p.lazy_normalization == selected)); + } + } + + #[cfg(all(feature = "lazy-normalization", feature = "serde"))] + #[test] + fn missing_lazy_selection_deserializes_as_eager() { + let parameters: ElasticNetParameters = serde_json::from_str("{}").unwrap(); + let search: ElasticNetSearchParameters = serde_json::from_str("{}").unwrap(); + assert!(!parameters.lazy_normalization); + assert!(!search.lazy_normalization); + let parameters = ElasticNetParameters::default().with_lazy_normalization(true); + let json = serde_json::to_string(¶meters).unwrap(); + let restored: ElasticNetParameters = serde_json::from_str(&json).unwrap(); + assert!(restored.lazy_normalization); + } + #[test] fn search_parameters() { let parameters = ElasticNetSearchParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: false, alpha: vec![0., 1.], max_iter: vec![10, 100], ..Default::default() @@ -534,6 +827,8 @@ mod tests { &x, &y, ElasticNetParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: false, alpha: 1.0, l1_ratio: 0.5, normalize: false, @@ -586,6 +881,8 @@ mod tests { &x, &y, ElasticNetParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: false, alpha: 1.0, l1_ratio: 1.0, normalize: true, @@ -599,6 +896,8 @@ mod tests { &x, &y, ElasticNetParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: false, alpha: 1.0, l1_ratio: 0.0, normalize: true, diff --git a/src/linear/lasso.rs b/src/linear/lasso.rs index 45493ea3..a57810eb 100644 --- a/src/linear/lasso.rs +++ b/src/linear/lasso.rs @@ -32,6 +32,8 @@ use crate::api::{Predictor, SupervisedEstimator}; use crate::error::Failed; use crate::linalg::basic::arrays::{Array1, Array2, ArrayView1}; use crate::linear::lasso_optimizer::InteriorPointOptimizer; +#[cfg(feature = "lazy-normalization")] +use crate::linear::lazy_normalization::LazyDesign; use crate::numbers::basenum::Number; use crate::numbers::floatnum::FloatNumber; use crate::numbers::realnum::RealNumber; @@ -41,6 +43,12 @@ use crate::numbers::realnum::RealNumber; #[derive(Debug, Clone)] #[must_use] pub struct LassoParameters { + #[cfg(feature = "lazy-normalization")] + #[cfg_attr(feature = "serde", serde(default))] + /// Borrow the design matrix instead of storing its normalized values. + /// This reduces memory use but can increase fit time. The default is false. + /// This option does not change whether normalization is applied. + pub lazy_normalization: bool, #[cfg_attr(feature = "serde", serde(default))] /// Controls the strength of the penalty to the loss function. pub alpha: f64, @@ -70,6 +78,24 @@ pub struct Lasso, Y: Arr } impl LassoParameters { + /// Select lazy matrix operations. Requires the `lazy-normalization` feature. + /// This reduces memory use but can increase fit time. + /// + /// ``` + /// # #[cfg(feature = "lazy-normalization")] + /// # { + /// use smartcore::linear::lasso::LassoParameters; + /// + /// let parameters = LassoParameters::default().with_lazy_normalization(true); + /// assert!(parameters.lazy_normalization); + /// # } + /// ``` + #[cfg(feature = "lazy-normalization")] + pub fn with_lazy_normalization(mut self, lazy_normalization: bool) -> Self { + self.lazy_normalization = lazy_normalization; + self + } + /// Regularization parameter. pub fn with_alpha(mut self, alpha: f64) -> Self { self.alpha = alpha; @@ -101,6 +127,8 @@ impl LassoParameters { impl Default for LassoParameters { fn default() -> Self { LassoParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: false, alpha: 1f64, normalize: true, tol: 1e-4, @@ -154,6 +182,13 @@ impl, Y: Array1> Pre #[derive(Debug, Clone)] #[must_use] pub struct LassoSearchParameters { + #[cfg(feature = "lazy-normalization")] + #[cfg_attr(feature = "serde", serde(default))] + /// Borrow the design matrix instead of storing its normalized values. + /// This reduces memory use but can increase fit time. The default is false. + /// This setting applies to all parameter combinations in the search. + /// This option does not change whether normalization is applied. + pub lazy_normalization: bool, #[cfg_attr(feature = "serde", serde(default))] /// Controls the strength of the penalty to the loss function. pub alpha: Vec, @@ -212,6 +247,8 @@ impl Iterator for LassoSearchParametersIterator { } let next = LassoParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: self.lasso_search_parameters.lazy_normalization, alpha: self.lasso_search_parameters.alpha[self.current_alpha], normalize: self.lasso_search_parameters.normalize[self.current_normalize], tol: self.lasso_search_parameters.tol[self.current_tol], @@ -257,6 +294,8 @@ impl Default for LassoSearchParameters { let default_params = LassoParameters::default(); LassoSearchParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: false, alpha: vec![default_params.alpha], normalize: vec![default_params.normalize], tol: vec![default_params.tol], @@ -300,6 +339,30 @@ impl, Y: Array1> Las let l1_reg = TX::from_f64(parameters.alpha * n as f64).unwrap(); + #[cfg(feature = "lazy-normalization")] + let (w, b) = if parameters.lazy_normalization { + Self::fit_lazy_coefficients(x, &y, ¶meters, l1_reg)? + } else { + Self::fit_eager_coefficients(x, &y, ¶meters, l1_reg)? + }; + #[cfg(not(feature = "lazy-normalization"))] + let (w, b) = Self::fit_eager_coefficients(x, &y, ¶meters, l1_reg)?; + + Ok(Lasso { + intercept: b, + coefficients: Some(w), + _phantom_ty: PhantomData, + _phantom_y: PhantomData, + }) + } + + fn fit_eager_coefficients( + x: &X, + y: &Vec, + parameters: &LassoParameters, + l1_reg: TX, + ) -> Result<(X, Option), Failed> { + let (_, p) = x.shape(); let (w, b) = if parameters.normalize { let (scaled_x, col_mean, col_std) = Self::rescale_x(x)?; @@ -307,7 +370,7 @@ impl, Y: Array1> Las let mut w = optimizer.optimize( &scaled_x, - &y, + y, l1_reg, parameters.max_iter, TX::from_f64(parameters.tol).unwrap(), @@ -334,7 +397,7 @@ impl, Y: Array1> Las let w = optimizer.optimize( x, - &y, + y, l1_reg, parameters.max_iter, TX::from_f64(parameters.tol).unwrap(), @@ -350,13 +413,72 @@ impl, Y: Array1> Las }, ) }; + Ok((w, b)) + } - Ok(Lasso { - intercept: b, - coefficients: Some(w), - _phantom_ty: PhantomData, - _phantom_y: PhantomData, - }) + #[cfg(feature = "lazy-normalization")] + fn fit_lazy_coefficients( + x: &X, + y: &Vec, + parameters: &LassoParameters, + l1_reg: TX, + ) -> Result<(X, Option), Failed> { + let (_, p) = x.shape(); + let (w, b) = if parameters.normalize { + let scaled_x = { + let (col_mean, col_std) = Self::normalization_stats(x)?; + LazyDesign::new(x, Some(col_mean), Some(col_std)) + }; + let (col_mean, col_std) = (scaled_x.centers(), scaled_x.scales()); + + let mut optimizer = InteriorPointOptimizer::new_lazy(&scaled_x); + + let mut w = optimizer.optimize_lazy( + &scaled_x, + y, + l1_reg, + parameters.max_iter, + TX::from_f64(parameters.tol).unwrap(), + parameters.fit_intercept, + )?; + + for (j, col_std_j) in col_std.iter().enumerate().take(p) { + w[j] /= *col_std_j; + } + + let b = if parameters.fit_intercept { + let mut xw_mean = TX::zero(); + for (i, col_mean_i) in col_mean.iter().enumerate().take(p) { + xw_mean += w[i] * *col_mean_i; + } + + Some(TX::from_f64(y.mean_by()).unwrap() - xw_mean) + } else { + None + }; + (X::from_column(&w), b) + } else { + let mut optimizer = InteriorPointOptimizer::new(x, p); + + let w = optimizer.optimize( + x, + y, + l1_reg, + parameters.max_iter, + TX::from_f64(parameters.tol).unwrap(), + parameters.fit_intercept, + )?; + + ( + X::from_column(&w), + if parameters.fit_intercept { + Some(TX::from_f64(y.mean_by()).unwrap()) + } else { + None + }, + ) + }; + Ok((w, b)) } /// Predict target values from `x` @@ -389,7 +511,7 @@ impl, Y: Array1> Las self.intercept.as_ref().unwrap() } - fn rescale_x(x: &X) -> Result<(X, Vec, Vec), Failed> { + fn normalization_stats(x: &X) -> Result<(Vec, Vec), Failed> { let col_mean: Vec = x .mean_by(0) .iter() @@ -407,6 +529,11 @@ impl, Y: Array1> Las } } + Ok((col_mean, col_std)) + } + + fn rescale_x(x: &X) -> Result<(X, Vec, Vec), Failed> { + let (col_mean, col_std) = Self::normalization_stats(x)?; let mut scaled_x = x.clone(); scaled_x.scale_mut(&col_mean, &col_std, 0); Ok((scaled_x, col_mean, col_std)) @@ -420,9 +547,150 @@ mod tests { use crate::linalg::basic::matrix::DenseMatrix; use crate::metrics::mean_absolute_error; + #[cfg(feature = "lazy-normalization")] + fn compare_eager_fit>(tolerance: f64) { + for lazy_normalization in [false, true] { + let values = [ + 1., 2., 2., 1., 3., 5., 4., 3., 5., 8., 6., 4., 7., 6., 8., 7., + ]; + let x = X::from_iterator(values.into_iter().map(|v| T::from_f64(v).unwrap()), 8, 2, 0); + let y: Vec = (0..8) + .map(|i| T::from_f64(3.0 + 2.0 * values[2 * i] - 0.4 * values[2 * i + 1]).unwrap()) + .collect(); + for normalize in [false, true] { + for fit_intercept in [false, true] { + let parameters = LassoParameters::default() + .with_alpha(0.03) + .with_lazy_normalization(lazy_normalization) + .with_normalize(normalize) + .with_fit_intercept(fit_intercept); + let (eager, centers, scales) = if normalize { + Lasso::>::rescale_x(&x).unwrap() + } else { + (x.clone(), vec![T::zero(); 2], vec![T::one(); 2]) + }; + let mut optimizer = InteriorPointOptimizer::new(&eager, 2); + let mut expected = optimizer + .optimize( + &eager, + &y, + T::from_f64(0.03 * 8.0).unwrap(), + parameters.max_iter, + T::from_f64(parameters.tol).unwrap(), + fit_intercept, + ) + .unwrap(); + for (w, scale) in expected.iter_mut().zip(scales) { + *w /= scale; + } + let model = Lasso::fit(&x, &y, parameters).unwrap(); + for (&a, &b) in model.coefficients().iterator(0).zip(&expected) { + assert!((a - b).abs().to_f64().unwrap() < tolerance, "{a} != {b}"); + } + if fit_intercept { + let intercept = T::from_f64(y.mean_by()).unwrap() + - expected.iter().zip(centers).map(|(&w, c)| w * c).sum::(); + assert!( + (*model.intercept() - intercept).abs().to_f64().unwrap() < tolerance + ); + let prediction = model.predict(&x).unwrap(); + for (i, &actual) in prediction.iter().enumerate() { + let expected = + (0..2).map(|j| *x.get((i, j)) * expected[j]).sum::() + intercept; + assert!((actual - expected).abs().to_f64().unwrap() < 10.0 * tolerance); + } + } else { + assert!(model.intercept.is_none()); + } + assert!( + x.iterator(0) + .zip(values) + .all(|(&actual, original)| actual == T::from_f64(original).unwrap()) + ); + } + } + } + } + + #[cfg(feature = "lazy-normalization")] + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn lazy_fit_matches_eager_for_normalization_and_intercept_options() { + compare_eager_fit::>(1e-8); + compare_eager_fit::>(5e-3); + #[cfg(feature = "ndarray-bindings")] + { + compare_eager_fit::>(1e-8); + compare_eager_fit::>(5e-3); + } + } + + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn normalization_rejects_constant_and_near_constant_columns() { + for step in [0.0, f64::EPSILON / 16.0] { + let x = DenseMatrix::from_2d_array(&[&[0.0, 1.0], &[step, 2.0], &[2.0 * step, 4.0]]) + .unwrap(); + let error = Lasso::fit(&x, &vec![1., 2., 3.], LassoParameters::default()).unwrap_err(); + assert!( + error + .to_string() + .contains("Cannot rescale constant column 0") + ); + #[cfg(feature = "lazy-normalization")] + { + let lazy_error = Lasso::fit( + &x, + &vec![1., 2., 3.], + LassoParameters::default().with_lazy_normalization(true), + ) + .unwrap_err(); + assert_eq!(lazy_error, error); + } + } + } + + #[cfg(feature = "lazy-normalization")] + #[test] + fn lazy_selection_defaults_and_search() { + assert!(!LassoParameters::default().lazy_normalization); + assert!(!LassoSearchParameters::default().lazy_normalization); + for selected in [false, true] { + let search = LassoSearchParameters { + lazy_normalization: selected, + alpha: vec![0.1, 0.2], + ..Default::default() + }; + let parameters: Vec<_> = search.into_iter().collect(); + assert_eq!(parameters.len(), 2); + assert!(parameters.iter().all(|p| p.lazy_normalization == selected)); + } + } + + #[cfg(all(feature = "lazy-normalization", feature = "serde"))] + #[test] + fn missing_lazy_selection_deserializes_as_eager() { + let parameters: LassoParameters = serde_json::from_str("{}").unwrap(); + let search: LassoSearchParameters = serde_json::from_str("{}").unwrap(); + assert!(!parameters.lazy_normalization); + assert!(!search.lazy_normalization); + let parameters = LassoParameters::default().with_lazy_normalization(true); + let json = serde_json::to_string(¶meters).unwrap(); + let restored: LassoParameters = serde_json::from_str(&json).unwrap(); + assert!(restored.lazy_normalization); + } + #[test] fn search_parameters() { let parameters = LassoSearchParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: false, alpha: vec![0., 1.], max_iter: vec![10, 100], fit_intercept: vec![false, true], @@ -493,6 +761,8 @@ mod tests { &x, &y, LassoParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: false, alpha: 0.1, normalize: false, tol: 1e-4, @@ -547,6 +817,8 @@ mod tests { &x, &y, LassoParameters { + #[cfg(feature = "lazy-normalization")] + lazy_normalization: false, alpha: 0.1, normalize: false, tol: 1e-8, diff --git a/src/linear/lasso_optimizer.rs b/src/linear/lasso_optimizer.rs index 5969347a..f3675d08 100644 --- a/src/linear/lasso_optimizer.rs +++ b/src/linear/lasso_optimizer.rs @@ -16,6 +16,9 @@ use crate::linalg::basic::arrays::{Array1, Array2, ArrayView1, MutArray, MutArra use crate::linear::bg_solver::BiconjugateGradientSolver; use crate::numbers::floatnum::FloatNumber; +#[cfg(feature = "lazy-normalization")] +use crate::linear::lazy_normalization::LazyDesign; + /// Interior Point Optimizer pub struct InteriorPointOptimizer> { ata: X, @@ -28,8 +31,12 @@ pub struct InteriorPointOptimizer> { impl> InteriorPointOptimizer { /// Initialize a new Interior Point Optimizer pub fn new(a: &X, n: usize) -> InteriorPointOptimizer { + Self::from_gram(a.ab(true, a, false), n) + } + + fn from_gram(ata: X, n: usize) -> Self { InteriorPointOptimizer { - ata: a.ab(true, a, false), + ata, d1: vec![T::zero(); n], d2: vec![T::zero(); n], prb: vec![T::zero(); n], @@ -37,6 +44,33 @@ impl> InteriorPointOptimizer { } } + #[cfg(feature = "lazy-normalization")] + pub(super) fn new_lazy(design: &LazyDesign<'_, T, X>) -> Self { + Self::from_gram(design.gram(), design.shape().1) + } + + #[cfg(feature = "lazy-normalization")] + pub(super) fn optimize_lazy( + &mut self, + design: &LazyDesign<'_, T, X>, + y: &Vec, + lambda: T, + max_iter: usize, + tol: T, + fit_intercept: bool, + ) -> Result, Failed> { + self.optimize_with_products( + design.shape(), + y, + lambda, + max_iter, + tol, + fit_intercept, + |v| design.matvec(v), + |v| design.transpose_matvec(v), + ) + } + /// Run the optimization pub fn optimize( &mut self, @@ -47,7 +81,34 @@ impl> InteriorPointOptimizer { tol: T, fit_intercept: bool, ) -> Result, Failed> { - let (n, p) = x.shape(); + self.optimize_with_products( + x.shape(), + y, + lambda, + max_iter, + tol, + fit_intercept, + |v| v.xa(true, x), + |v| v.xa(false, x), + ) + } + + fn optimize_with_products( + &mut self, + shape: (usize, usize), + y: &Vec, + lambda: T, + max_iter: usize, + tol: T, + fit_intercept: bool, + matvec: F, + transpose_matvec: G, + ) -> Result, Failed> + where + F: Fn(&Vec) -> Vec, + G: Fn(&Vec) -> Vec, + { + let (n, p) = shape; let p_f64 = T::from_usize(p).unwrap(); let lambda = lambda.max(T::epsilon()); @@ -96,7 +157,7 @@ impl> InteriorPointOptimizer { let lambda_f64 = lambda.to_f64().unwrap(); for ntiter in 0..max_iter { - let mut z = w.xa(true, x); + let mut z = matvec(&w); for i in 0..n { z[i] -= y[i]; @@ -104,7 +165,7 @@ impl> InteriorPointOptimizer { } // CALCULATE DUALITY GAP - let xnu = nu.xa(false, x); + let xnu = transpose_matvec(&nu); let max_xnu = xnu.norm(f64::INFINITY); if max_xnu > lambda_f64 { let lnu = T::from_f64(lambda_f64 / max_xnu).unwrap(); @@ -136,7 +197,7 @@ impl> InteriorPointOptimizer { self.d2[i] = (q1i * q1i - q2i * q2i) / t; } - let mut gradphi = z.xa(false, x); + let mut gradphi = transpose_matvec(&z); for i in 0..p { let g1 = T::two() * gradphi[i] - (q1[i] - q2[i]) / t; @@ -157,7 +218,7 @@ impl> InteriorPointOptimizer { pcgtol *= min_pcgtol; } - let error = self.solve_mut(x, &grad, &mut dxu, pcgtol, pcgmaxi)?; + let error = self.solve_mut(&self.ata, &grad, &mut dxu, pcgtol, pcgmaxi)?; if error > pcgtol { pitr = pcgmaxi; } @@ -184,7 +245,7 @@ impl> InteriorPointOptimizer { .fold(T::neg_infinity(), |max, v| v.max(max)) < T::zero() { - let mut newz = neww.xa(true, x); + let mut newz = matvec(&neww); for i in 0..n { newz[i] -= y[i]; } diff --git a/src/linear/lazy_normalization.rs b/src/linear/lazy_normalization.rs new file mode 100644 index 00000000..783f7446 --- /dev/null +++ b/src/linear/lazy_normalization.rs @@ -0,0 +1,830 @@ +//! Borrowed normalized designs for the interior-point method of Kim et al. +//! +//! See . + +use std::convert::Infallible; +use std::marker::PhantomData; + +use lazymatrix::{ + DotSlice, ElemDivAssign, LazyMatrix, MatTransposeVecInto, MatVecInto, MatrixErrorType, + MatrixShape, MatrixWrite, Scalar, ScaledSubSlice, SubScalarAssign, SumEntries, VectorView, + WeightedGramInto, WeightedGramKernel, +}; + +use crate::linalg::basic::arrays::{Array1, Array2}; +use crate::numbers::floatnum::FloatNumber; + +/// Keep the design borrowed while coefficient matrices retain the caller's storage type. +pub(super) struct LazyDesign<'a, T: FloatNumber, X: Array2> { + normalized: LazyMatrix, T>, + augmentation: Option<(T, T)>, +} + +impl<'a, T: FloatNumber, X: Array2> LazyDesign<'a, T, X> { + pub(super) fn new(x: &'a X, centers: Option>, scales: Option>) -> Self { + Self { + normalized: LazyMatrix::from_parts(MatrixAdapter(x, PhantomData), centers, scales), + augmentation: None, + } + } + + pub(super) fn with_l2(mut self, penalty: T) -> Self { + let gamma = T::one() / (T::one() + penalty).sqrt(); + self.augmentation = Some((gamma, gamma * penalty.sqrt())); + self + } + + pub(super) fn shape(&self) -> (usize, usize) { + let n = self.normalized.nrows(); + let p = self.normalized.ncols(); + (n + if self.augmentation.is_some() { p } else { 0 }, p) + } + + pub(super) fn centers(&self) -> &[T] { + self.normalized.centers().unwrap() + } + + pub(super) fn scales(&self) -> &[T] { + self.normalized.scales().unwrap() + } + + pub(super) fn gamma(&self) -> T { + self.augmentation.map_or(T::one(), |(gamma, _)| gamma) + } + + pub(super) fn matvec(&self, x: &Vec) -> Vec { + let mut out = Vec::with_capacity(self.shape().0); + out.resize(self.normalized.nrows(), T::zero()); + self.normalized.matvec_into(x, &mut out).unwrap(); + if let Some((gamma, padding)) = self.augmentation { + for value in &mut out { + *value *= gamma; + } + out.extend(x.iter().map(|&value| padding * value)); + } + out + } + + pub(super) fn transpose_matvec(&self, x: &[T]) -> Vec { + check_vector_length(x.len(), self.shape().0); + let n = self.normalized.nrows(); + let mut out = Vec::zeros(self.normalized.ncols()); + self.normalized + .mat_transpose_vec_into(&VectorRef(&x[..n]), &mut out) + .unwrap(); + if let Some((gamma, padding)) = self.augmentation { + for (value, &tail) in out.iter_mut().zip(&x[n..]) { + *value = gamma * *value + padding * tail; + } + } + out + } + + pub(super) fn gram(&self) -> X { + let p = self.normalized.ncols(); + let mut out = X::zeros(p, p); + self.normalized + .weighted_gram_into( + &UnitWeights(self.normalized.nrows()), + &mut MatrixOutput(&mut out, PhantomData), + ) + .unwrap(); + if let Some((gamma, padding)) = self.augmentation { + out.mul_scalar_mut(gamma * gamma); + for j in 0..p { + out.add_element_mut((j, j), padding * padding); + } + } + out + } +} + +struct MatrixAdapter<'a, T: FloatNumber, X: Array2>(&'a X, PhantomData); + +impl> MatrixShape for MatrixAdapter<'_, T, X> { + fn nrows(&self) -> usize { + self.0.shape().0 + } + fn ncols(&self) -> usize { + self.0.shape().1 + } +} + +impl> MatrixErrorType for MatrixAdapter<'_, T, X> { + type Error = Infallible; +} + +impl> MatVecInto> for MatrixAdapter<'_, T, X> { + fn matvec_into(&self, x: &Vec, out: &mut Vec) -> Result<(), Infallible> { + self.matvec_normalized_into::(x, None, None, out) + } + + fn matvec_normalized_into( + &self, + x: &Vec, + centers: Option<&[F]>, + scales: Option<&[F]>, + out: &mut Vec, + ) -> Result<(), Infallible> + where + Vec: ElemDivAssign + DotSlice + SubScalarAssign, + { + check_vector_length(x.len(), self.ncols()); + check_vector_length(out.len(), self.nrows()); + out.fill(T::zero()); + if self.nrows() == 0 { + return Ok(()); + } + // Column traversal reuses each view while preserving each row's summation order. + for (j, &v) in x.iter().enumerate() { + let center = centers.map_or(T::zero(), |c| T::from(c[j]).unwrap()); + let scale = scales.map_or(T::one(), |s| T::from(s[j]).unwrap()); + // Keep division when a reciprocal would overflow or lose normal precision. + let reciprocal = T::one() / scale; + let use_reciprocal = scale.is_normal() && reciprocal.is_normal(); + let column = self.0.get_col(j); + for (&raw, result) in column.iterator(0).zip(out.iter_mut()) { + // Center before accumulation to retain small variations at large offsets. + let normalized = if use_reciprocal { + (raw - center) * reciprocal + } else { + (raw - center) / scale + }; + *result += normalized * v; + } + } + Ok(()) + } +} + +struct VectorRef<'a, T>(&'a [T]); + +impl SumEntries for VectorRef<'_, T> { + fn sum_entries(&self) -> T { + self.0.iter().copied().sum() + } +} + +impl<'v, T: FloatNumber, X: Array2> MatTransposeVecInto, Vec> + for MatrixAdapter<'_, T, X> +{ + fn mat_transpose_vec_into( + &self, + x: &VectorRef<'v, T>, + out: &mut Vec, + ) -> Result<(), Infallible> { + self.mat_transpose_vec_normalized_into::(x, None, None, out) + } + + fn mat_transpose_vec_normalized_into( + &self, + x: &VectorRef<'v, T>, + centers: Option<&[F]>, + scales: Option<&[F]>, + out: &mut Vec, + ) -> Result<(), Infallible> + where + VectorRef<'v, T>: SumEntries, + Vec: ScaledSubSlice + ElemDivAssign, + { + check_vector_length(x.0.len(), self.nrows()); + check_vector_length(out.len(), self.ncols()); + if self.nrows() == 0 { + out.fill(T::zero()); + return Ok(()); + } + for (j, result) in out.iter_mut().enumerate() { + let center = centers.map_or(T::zero(), |c| T::from(c[j]).unwrap()); + let scale = scales.map_or(T::one(), |s| T::from(s[j]).unwrap()); + // Keep division when a reciprocal would overflow or lose normal precision. + let reciprocal = T::one() / scale; + let use_reciprocal = scale.is_normal() && reciprocal.is_normal(); + let column = self.0.get_col(j); + *result = column + .iterator(0) + .zip(x.0) + .fold(T::zero(), |sum, (&raw, &v)| { + let normalized = if use_reciprocal { + (raw - center) * reciprocal + } else { + (raw - center) / scale + }; + sum + normalized * v + }); + } + Ok(()) + } +} + +impl> WeightedGramKernel for MatrixAdapter<'_, T, X> { + fn weighted_gram_normalized_into( + &self, + weights: &W, + centers: Option<&[T]>, + scales: Option<&[T]>, + out: &mut O, + ) -> Result<(), Infallible> + where + W: VectorView + ?Sized, + O: MatrixWrite + ?Sized, + { + // Panels reuse normalized values across column pairs without storing the design. + const ROWS: usize = 256; + const COLS: usize = 32; + let (n, p) = self.0.shape(); + let rows = n.min(ROWS); + let cols = p.min(COLS); + let mut left = Vec::zeros(rows * cols); + let mut right = Vec::zeros(rows * cols); + let mut sums = Vec::zeros(cols * cols); + for j in (0..p).step_by(COLS) { + let nj = (p - j).min(COLS); + for k in (0..=j).step_by(COLS) { + let nk = (p - k).min(COLS); + sums.fill(T::zero()); + for row in (0..n).step_by(ROWS) { + let nr = (n - row).min(ROWS); + self.normalized_panel(row, nr, j, nj, centers, scales, &mut left); + if j == k { + right[..nr * nk].copy_from_slice(&left[..nr * nj]); + } else { + self.normalized_panel(row, nr, k, nk, centers, scales, &mut right); + } + for column in left[..nr * nj].chunks_exact_mut(nr) { + for (i, value) in column.iter_mut().enumerate() { + *value = weights.get(row + i) * *value; + } + } + for a in 0..nj { + let left_column = &left[a * nr..(a + 1) * nr]; + let width = if j == k { a + 1 } else { nk }; + for b in 0..width { + let right_column = &right[b * nr..(b + 1) * nr]; + let sum = &mut sums[a * cols + b]; + // Carry the accumulator across panels to retain row order. + for (&x, &y) in left_column.iter().zip(right_column) { + *sum += x * y; + } + } + } + } + for a in 0..nj { + let width = if j == k { a + 1 } else { nk }; + for b in 0..width { + let value = sums[a * cols + b]; + out.set(j + a, k + b, value); + if j + a != k + b { + out.set(k + b, j + a, value); + } + } + } + } + } + Ok(()) + } +} + +impl> MatrixAdapter<'_, T, X> { + fn normalized_panel( + &self, + row: usize, + nrows: usize, + column: usize, + ncols: usize, + centers: Option<&[T]>, + scales: Option<&[T]>, + out: &mut [T], + ) { + for j in 0..ncols { + let center = centers.map_or(T::zero(), |c| c[column + j]); + let scale = scales.map_or(T::one(), |s| s[column + j]); + let view = self.0.slice(row..row + nrows, column + j..column + j + 1); + for (value, &raw) in out[j * nrows..(j + 1) * nrows] + .iter_mut() + .zip(view.iterator(0)) + { + *value = (raw - center) / scale; + } + } + } +} + +struct UnitWeights(usize); + +impl VectorView for UnitWeights { + fn len(&self) -> usize { + self.0 + } + fn get(&self, index: usize) -> T { + if index >= self.0 { + std::panic::panic_any("weight index is out of bounds"); + } + T::one() + } +} + +fn check_vector_length(actual: usize, expected: usize) { + if actual != expected { + std::panic::panic_any("matrix and vector dimensions do not agree"); + } +} + +struct MatrixOutput<'a, T: FloatNumber, X: Array2>(&'a mut X, PhantomData); + +impl> MatrixShape for MatrixOutput<'_, T, X> { + fn nrows(&self) -> usize { + self.0.shape().0 + } + fn ncols(&self) -> usize { + self.0.shape().1 + } +} + +impl> MatrixWrite for MatrixOutput<'_, T, X> { + fn set(&mut self, row: usize, column: usize, value: T) { + self.0.set((row, column), value); + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::linalg::basic::arrays::{Array, Array1, Array2, MutArray, MutArrayView2}; + use crate::linalg::basic::matrix::DenseMatrix; + use crate::linear::lasso_optimizer::InteriorPointOptimizer; + use crate::numbers::floatnum::FloatNumber; + + fn check_products>(x: X, tolerance: f64) { + let value = |v| T::from_f64(v).unwrap(); + let original = x.clone(); + for centered in [false, true] { + for scaled in [false, true] { + let centers = centered.then(|| vec![value(1.0), value(2.0)]); + let scales = scaled.then(|| vec![value(2.0), value(4.0)]); + let mut eager = x.clone(); + eager.scale_mut( + centers.as_deref().unwrap_or(&[T::zero(); 2]), + scales.as_deref().unwrap_or(&[T::one(); 2]), + 0, + ); + let design = LazyDesign::new(&x, centers, scales); + let coefficients = vec![value(2.0), value(-1.0)]; + let response = vec![value(1.0), value(-2.0), value(3.0)]; + assert_close( + &design.matvec(&coefficients), + &coefficients.xa(true, &eager), + tolerance, + ); + assert_close( + &design.transpose_matvec(&response), + &response.xa(false, &eager), + tolerance, + ); + let gram = design.gram(); + assert_close( + &gram.iterator(0).copied().collect::>(), + &eager + .ab(true, &eager, false) + .iterator(0) + .copied() + .collect::>(), + tolerance, + ); + let left = design + .matvec(&coefficients) + .iter() + .zip(&response) + .map(|(&a, &b)| a * b) + .sum::(); + let right = coefficients + .iter() + .zip(design.transpose_matvec(&response)) + .map(|(&a, b)| a * b) + .sum::(); + assert!((left - right).abs().to_f64().unwrap() <= tolerance); + } + } + assert!(x.iterator(0).zip(original.iterator(0)).all(|(a, b)| a == b)); + } + + fn assert_close(actual: &[T], expected: &[T], tolerance: f64) { + assert_eq!(actual.len(), expected.len()); + for (&a, &b) in actual.iter().zip(expected) { + assert!((a - b).abs().to_f64().unwrap() <= tolerance, "{a} != {b}"); + } + } + + fn assert_relative_close(actual: &[T], expected: &[T], tolerance: f64) { + assert_eq!(actual.len(), expected.len()); + for (&a, &b) in actual.iter().zip(expected) { + let scale = a.abs().max(b.abs()).max(T::one()).to_f64().unwrap(); + assert!( + (a - b).abs().to_f64().unwrap() <= tolerance * scale, + "{a} != {b}" + ); + } + } + + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn adapter_products_overwrite_output_without_normalization() { + for (n, p) in [(3, 2), (0, 2), (3, 0)] { + let x = DenseMatrix::from_iterator((0..n * p).map(|i| i as f64 + 1.0), n, p, 0); + let adapter = MatrixAdapter(&x, PhantomData); + let coefficients = vec![2.0; p]; + let response = vec![-1.0; n]; + let mut product = vec![f64::NAN; n]; + adapter.matvec_into(&coefficients, &mut product).unwrap(); + assert_eq!(product, coefficients.xa(true, &x)); + + let mut transpose_product = vec![f64::NAN; p]; + adapter + .mat_transpose_vec_into(&VectorRef(&response), &mut transpose_product) + .unwrap(); + assert_eq!(transpose_product, response.xa(false, &x)); + assert_eq!(VectorRef(&response).sum_entries(), -(n as f64)); + + let design = LazyDesign::new(&x, None, None); + let mut gram = DenseMatrix::fill(p, p, f64::NAN); + design + .normalized + .weighted_gram_into(&UnitWeights(n), &mut MatrixOutput(&mut gram, PhantomData)) + .unwrap(); + assert_eq!(gram, x.ab(true, &x, false)); + } + } + + #[cfg(not(target_arch = "wasm32"))] + #[test] + fn adapter_products_reject_incorrect_vector_lengths() { + use std::panic::{AssertUnwindSafe, catch_unwind}; + + let x = DenseMatrix::from_2d_array(&[&[1., 2.], &[3., 4.], &[5., 6.]]).unwrap(); + let adapter = MatrixAdapter(&x, PhantomData); + for (input, output) in [(1, 3), (3, 3), (2, 2), (2, 4)] { + let mut out = vec![0.; output]; + let error = catch_unwind(AssertUnwindSafe(|| { + adapter.matvec_into(&vec![1.; input], &mut out).unwrap(); + })) + .unwrap_err(); + assert_eq!( + error.downcast_ref::<&str>(), + Some(&"matrix and vector dimensions do not agree") + ); + } + for (input, output) in [(2, 2), (4, 2), (3, 1), (3, 3)] { + let mut out = vec![0.; output]; + let error = catch_unwind(AssertUnwindSafe(|| { + adapter + .mat_transpose_vec_into(&VectorRef(&vec![1.; input]), &mut out) + .unwrap(); + })) + .unwrap_err(); + assert_eq!( + error.downcast_ref::<&str>(), + Some(&"matrix and vector dimensions do not agree") + ); + } + for augmented in [false, true] { + let design = LazyDesign::new(&x, None, None); + let design = if augmented { + design.with_l2(0.7) + } else { + design + }; + for length in [design.shape().0 - 1, design.shape().0 + 1] { + let error = catch_unwind(AssertUnwindSafe(|| { + design.transpose_matvec(&vec![1.; length]); + })) + .unwrap_err(); + assert_eq!( + error.downcast_ref::<&str>(), + Some(&"matrix and vector dimensions do not agree") + ); + } + } + } + + #[cfg(not(target_arch = "wasm32"))] + #[test] + #[should_panic(expected = "weight index is out of bounds")] + fn unit_weights_reject_out_of_bounds_index() { + let _: f64 = UnitWeights(3).get(3); + } + + fn check_extreme_product_scales() { + for scale in [ + T::min_positive_value() * T::epsilon(), + T::min_positive_value(), + ::max_value(), + -::max_value(), + T::infinity(), + T::neg_infinity(), + T::nan(), + ] { + let raw = if scale.is_finite() { scale } else { T::one() }; + let x = DenseMatrix::new(2, 1, vec![raw, -raw], true).unwrap(); + let design = LazyDesign::new(&x, None, Some(vec![scale])); + let expected = vec![raw / scale, -raw / scale]; + let actual = design.matvec(&vec![T::one()]); + for (&a, &b) in actual.iter().zip(&expected) { + if b.is_nan() { + assert!(a.is_nan()); + } else { + assert_eq!(a, b); + } + } + let actual = design.transpose_matvec(&[T::one(), -T::one()]); + let expected = expected[0] - expected[1]; + if expected.is_nan() { + assert!(actual[0].is_nan()); + } else { + assert_eq!(actual[0], expected); + } + } + } + + #[test] + fn products_preserve_extreme_scales() { + check_extreme_product_scales::(); + check_extreme_product_scales::(); + } + + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn products_match_eager_for_both_dense_layouts_and_precisions() { + for column_major in [false, true] { + let data = if column_major { + vec![1., 2., 4., 4., 6., 10.] + } else { + vec![1., 4., 2., 6., 4., 10.] + }; + check_products( + DenseMatrix::new(3, 2, data.clone(), column_major).unwrap(), + 1e-12, + ); + check_products( + DenseMatrix::new(3, 2, data.iter().map(|&v| v as f32).collect(), column_major) + .unwrap(), + 1e-5, + ); + } + } + + #[cfg(feature = "ndarray-bindings")] + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn products_match_eager_for_ndarray() { + use ndarray::ShapeBuilder; + for column_major in [false, true] { + let shape = (3, 2).set_f(column_major); + let data = if column_major { + vec![1., 2., 4., 4., 6., 10.] + } else { + vec![1., 4., 2., 6., 4., 10.] + }; + check_products( + ndarray::Array2::from_shape_vec(shape, data.clone()).unwrap(), + 1e-12, + ); + check_products( + ndarray::Array2::from_shape_vec(shape, data.iter().map(|&v| v as f32).collect()) + .unwrap(), + 1e-5, + ); + } + } + + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn centered_products_preserve_small_variation_at_large_offsets() { + let x = DenseMatrix::from_2d_array(&[ + &[1e12 + 1., 1e12 + 4.], + &[1e12 + 2., 1e12 + 6.], + &[1e12 + 4., 1e12 + 10.], + ]) + .unwrap(); + let design = LazyDesign::new(&x, Some(vec![1e12 + 1., 1e12 + 2.]), Some(vec![2., 4.])); + assert_eq!(design.matvec(&vec![2., -1.]), vec![-0.5, 0., 1.]); + assert_eq!(design.transpose_matvec(&vec![1., -2., 3.]), vec![3.5, 4.5]); + let nonbinary = LazyDesign::new(&x, Some(vec![1e12 + 1., 1e12 + 2.]), Some(vec![3., 7.])); + let mut eager = x.clone(); + eager.scale_mut(&[1e12 + 1., 1e12 + 2.], &[3., 7.], 0); + assert_relative_close( + &nonbinary.matvec(&vec![2., -1.]), + &vec![2., -1.].xa(true, &eager), + 8.0 * f64::EPSILON, + ); + assert_relative_close( + &nonbinary.transpose_matvec(&[1., -2., 3.]), + &vec![1., -2., 3.].xa(false, &eager), + 8.0 * f64::EPSILON, + ); + let gram = design.gram(); + assert_eq!(*gram.get((0, 0)), 2.5); + assert_eq!(*gram.get((0, 1)), 3.5); + assert_eq!(*gram.get((1, 1)), 5.25); + } + + fn check_panel_gram>() { + for (n, p) in [ + (0, 3), + (3, 0), + (1, 1), + (255, 31), + (256, 32), + (257, 33), + (513, 65), + ] { + for axis in [0, 1] { + let x = X::from_iterator( + (0..n * p).map(|index| { + let (i, j) = if axis == 0 { + (index / p, index % p) + } else { + (index % n, index / n) + }; + T::from_f64(100_000.0 + 2.0 * j as f64 + ((i * 7 + j) % 23) as f64 / 8.0) + .unwrap() + }), + n, + p, + axis, + ); + let weights: Vec = (0..n) + .map(|i| T::from_f64([1.0, 0.0, -0.5, 2.0][i % 4]).unwrap()) + .collect(); + for centered in [false, true] { + for scaled in [false, true] { + let centers: Vec = (0..p) + .map(|j| { + T::from_f64(if centered { + 100_000.0 + 2.0 * j as f64 + } else { + 0.0 + }) + .unwrap() + }) + .collect(); + let scales: Vec = (0..p) + .map(|j| { + T::from_f64(if scaled { [3.0, 0.5, -2.0][j % 3] } else { 1.0 }) + .unwrap() + }) + .collect(); + let mut eager = x.clone(); + eager.scale_mut(¢ers, &scales, 0); + let design = LazyDesign::new( + &x, + centered.then_some(centers), + scaled.then_some(scales), + ); + let coefficients: Vec = (0..p) + .map(|j| T::from_f64((j % 5) as f64 - 2.0).unwrap()) + .collect(); + let tolerance = 128.0 * T::epsilon().to_f64().unwrap(); + assert_relative_close( + &design.matvec(&coefficients), + &coefficients.xa(true, &eager), + tolerance, + ); + assert_relative_close( + &design.transpose_matvec(&weights), + &weights.xa(false, &eager), + tolerance, + ); + let mut actual = X::fill(p, p, T::nan()); + design + .normalized + .weighted_gram_into( + &weights, + &mut MatrixOutput(&mut actual, PhantomData), + ) + .unwrap(); + for j in 0..p { + for k in 0..=j { + let expected = (0..n).fold(T::zero(), |sum, i| { + sum + weights[i] * *eager.get((i, j)) * *eager.get((i, k)) + }); + assert_eq!( + *actual.get((j, k)), + expected, + "shape=({n},{p}), axis={axis}, pair=({j},{k})" + ); + assert_eq!(*actual.get((k, j)), expected); + } + } + } + } + } + } + } + + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn panel_gram_preserves_eager_accumulation_across_boundaries() { + check_panel_gram::>(); + check_panel_gram::>(); + #[cfg(feature = "ndarray-bindings")] + { + check_panel_gram::>(); + check_panel_gram::>(); + } + } + + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn panel_gram_preserves_nonfinite_products() { + let x = DenseMatrix::from_2d_array(&[ + &[0.0, f64::INFINITY, 2.0], + &[1.0, 3.0, f64::NAN], + &[-2.0, 0.0, 4.0], + ]) + .unwrap(); + let weights = vec![0.0, 1.0, f64::INFINITY]; + let design = LazyDesign::new(&x, None, Some(vec![1.0, -2.0, f64::INFINITY])); + let mut eager = x.clone(); + eager.scale_mut(&[0.0; 3], &[1.0, -2.0, f64::INFINITY], 0); + let mut actual = DenseMatrix::zeros(3, 3); + design + .normalized + .weighted_gram_into(&weights, &mut MatrixOutput(&mut actual, PhantomData)) + .unwrap(); + for j in 0..3 { + for k in 0..=j { + let expected = (0..3).fold(0.0, |sum, i| { + sum + weights[i] * *eager.get((i, j)) * *eager.get((i, k)) + }); + for index in [(j, k), (k, j)] { + if expected.is_nan() { + assert!(actual.get(index).is_nan()); + } else { + assert_eq!(*actual.get(index), expected); + } + } + } + } + } + + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn implicit_augmentation_matches_materialized_products_and_fit() { + let x = DenseMatrix::from_2d_array(&[&[1., 4.], &[2., 6.], &[4., 10.]]).unwrap(); + for penalty in [0.0_f64, 0.7, 100.0] { + let gamma = 1.0 / (1.0 + penalty).sqrt(); + let padding = gamma * penalty.sqrt(); + let design = + LazyDesign::new(&x, Some(vec![1., 2.]), Some(vec![2., 4.])).with_l2(penalty); + let mut eager = DenseMatrix::zeros(5, 2); + for j in 0..2 { + for i in 0..3 { + eager.set((i, j), gamma * (x.get((i, j)) - [1., 2.][j]) / [2., 4.][j]); + } + eager.set((3 + j, j), padding); + } + let w = vec![2., -1.]; + let y = vec![7., 9., 13., 0., 0.]; + assert_close(&design.matvec(&w), &w.xa(true, &eager), 1e-12); + assert_close(&design.transpose_matvec(&y), &y.xa(false, &eager), 1e-12); + assert_close( + &design.gram().iterator(0).copied().collect::>(), + &eager + .ab(true, &eager, false) + .iterator(0) + .copied() + .collect::>(), + 1e-12, + ); + let mut lazy_opt = InteriorPointOptimizer::new_lazy(&design); + let mut eager_opt = InteriorPointOptimizer::new(&eager, 2); + let actual = lazy_opt + .optimize_lazy(&design, &y, 0.2 * gamma, 1000, 1e-8, true) + .unwrap(); + let expected = eager_opt + .optimize(&eager, &y, 0.2 * gamma, 1000, 1e-8, true) + .unwrap(); + assert_close(&actual, &expected, 1e-6); + } + } +} diff --git a/src/linear/mod.rs b/src/linear/mod.rs index fb3bf8a4..d3f4f129 100644 --- a/src/linear/mod.rs +++ b/src/linear/mod.rs @@ -24,6 +24,8 @@ pub mod bg_solver; pub mod elastic_net; pub mod lasso; pub mod lasso_optimizer; +#[cfg(feature = "lazy-normalization")] +mod lazy_normalization; pub mod linear_regression; pub mod logistic_regression; pub mod ridge_regression;