diff --git a/Cargo.toml b/Cargo.toml index 8b545fb6e..5fd427b06 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -75,6 +75,7 @@ members = [ "algorithms/linfa-ica", "algorithms/linfa-bayes", "algorithms/linfa-elasticnet", + "algorithms/linfa-pls", "datasets", ] diff --git a/README.md b/README.md index 56babd58e..396e848a8 100644 --- a/README.md +++ b/README.md @@ -41,6 +41,7 @@ Where does `linfa` stand right now? [Are we learning yet?](http://www.arewelearn | [hierarchical](algorithms/linfa-hierarchical/) | Agglomerative hierarchical clustering | Tested | Unsupervised learning | Cluster and build hierarchy of clusters | | [bayes](algorithms/linfa-bayes/) | Naive Bayes | Tested | Supervised learning | Contains Gaussian Naive Bayes | | [ica](algorithms/linfa-ica/) | Independent component analysis | Tested | Unsupervised learning | Contains FastICA implementation | +| [pls](algorithms/linfa-pls/) | Partial Least Squares | Tested | Supervised learning | Contains PLS estimators for dimensionality reduction and regression | We believe that only a significant community effort can nurture, build, and sustain a machine learning ecosystem in Rust - there is no other way forward. diff --git a/algorithms/linfa-pls/Cargo.toml b/algorithms/linfa-pls/Cargo.toml new file mode 100644 index 000000000..ea6fa517c --- /dev/null +++ b/algorithms/linfa-pls/Cargo.toml @@ -0,0 +1,40 @@ +[package] +name = "linfa-pls" +version = "0.3.1" +edition = "2018" +authors = ["relf "] +description = "Partial Least Squares family methods" +license = "MIT/Apache-2.0" + +repository = "https://github.com/rust-ml/linfa" +readme = "README.md" + +keywords = ["pls", "machine-learning", "linfa", "supervised"] +categories = ["algorithms", "mathematics", "science"] + +[features] +default = [] +serde = ["serde_crate", "ndarray/serde"] + +[dependencies.serde_crate] +package = "serde" +optional = true +version = "1.0" +default-features = false +features = ["std", "derive"] + +[dependencies] +ndarray = { version = "0.13", default-features=false } +ndarray-linalg = "0.12" +ndarray-stats = "0.3" +ndarray-rand = "0.11" +rand_isaac = "0.2" +num-traits = "0.2" +paste = "1.0" + +linfa = { version = "0.3.1", path = "../.." } + +[dev-dependencies] +linfa-datasets = { version = "0.3.1", path = "../../datasets", features = ["linnerud"] } +rand_isaac = "0.2" +approx = "0.3" diff --git a/algorithms/linfa-pls/examples/pls_regression.rs b/algorithms/linfa-pls/examples/pls_regression.rs new file mode 100644 index 000000000..2cea60e48 --- /dev/null +++ b/algorithms/linfa-pls/examples/pls_regression.rs @@ -0,0 +1,40 @@ +use linfa::prelude::*; +use linfa_pls::{PlsRegression, Result}; +use ndarray::{Array, Array1, Array2}; +use ndarray_rand::rand::SeedableRng; +use ndarray_rand::rand_distr::StandardNormal; +use ndarray_rand::RandomExt; +use rand_isaac::Isaac64Rng; + +#[allow(clippy::many_single_char_names)] +fn main() -> Result<()> { + let n = 1000; + let q = 3; + let p = 10; + let mut rng = Isaac64Rng::seed_from_u64(42); + + // X shape (n, p) random + let x: Array2 = Array::random_using((n, p), StandardNormal, &mut rng); + + // B shape (p, q) such that B[0, ..] = 1, B[1, ..] = 2; otherwise zero + let mut b: Array2 = Array2::zeros((p, q)); + b.row_mut(0).assign(&Array1::ones(q)); + b.row_mut(1).assign(&Array1::from_elem(q, 2.)); + + // Y shape (n, q) such that yj = 1*x1 + 2*x2 + noise(Normal(5, 1)) + let y = x.dot(&b) + Array::random_using((n, q), StandardNormal, &mut rng).mapv(|v: f64| v + 5.); + + let ds = Dataset::new(x, y); + let pls = PlsRegression::params(3) + .scale(true) + .max_iterations(200) + .fit(&ds)?; + + println!("True B (such that: Y = XB + noise)"); + println!("{:?}", b); + + // PLS regression coefficients is an estimation of B + println!("Estimated B"); + println!("{:1.1}", pls.coefficients()); + Ok(()) +} diff --git a/algorithms/linfa-pls/src/errors.rs b/algorithms/linfa-pls/src/errors.rs new file mode 100644 index 000000000..a39c568bb --- /dev/null +++ b/algorithms/linfa-pls/src/errors.rs @@ -0,0 +1,36 @@ +use ndarray_linalg::error::LinalgError; +use std::error::Error; +use std::fmt::{self, Display}; + +pub type Result = std::result::Result; + +#[derive(Debug)] +pub enum PlsError { + NotEnoughSamplesError(String), + BadComponentNumberError(String), + PowerMethodNotConvergedError(String), + LinalgError(String), +} + +impl Display for PlsError { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + match self { + Self::NotEnoughSamplesError(message) => write!(f, "Not enough samples: {}", message), + Self::BadComponentNumberError(message) => { + write!(f, "Bad component number: {}", message) + } + Self::PowerMethodNotConvergedError(message) => { + write!(f, "Power method not converged: {}", message) + } + Self::LinalgError(message) => write!(f, "Linear algebra error: {}", message), + } + } +} + +impl Error for PlsError {} + +impl From for PlsError { + fn from(error: LinalgError) -> PlsError { + PlsError::LinalgError(error.to_string()) + } +} diff --git a/algorithms/linfa-pls/src/lib.rs b/algorithms/linfa-pls/src/lib.rs new file mode 100644 index 000000000..557795112 --- /dev/null +++ b/algorithms/linfa-pls/src/lib.rs @@ -0,0 +1,242 @@ +//! # Partial Least Squares +//! +//! `linfa-pls` provides an implementation of methods in the PLS (Partial Least Squares) family. +//! The PLS method is a statistical method that finds a linear relationship between +//! input variables and output variables by projecting them onto a new subspace formed +//! by newly chosen variables (aka latent variables), which are linear +//! combinations of the input variables. The subspace is choosen to maximize the +//! covariance between responses and independant variables. +//! +//! This approach is particularly useful when the original data are characterized by +//! a large number of highly collinear variables measured on a small number of samples. +//! +//! The implementation is a port of the scikit-learn 0.24 cross-decomposition module. +//! +//! ## References +//! +//! * [A survey of Partial Least Squares (PLS) methods, with emphasis on the two-block case JA Wegelin](https://stat.uw.edu/sites/default/files/files/reports/2000/tr371.pdf) +//! * [Scikit-Learn User Guide](https://scikit-learn.org/stable/modules/cross_decomposition.html#cross-decomposition) +//! +//! ## Example +//! +//! ```rust, ignore +//! use linfa::prelude::*; +//! use linfa_pls::{errors::Result, PlsRegression}; +//! use ndarray::array; +//! +//! // Load linnerud datase 20 samples, 3 input features, 3 output features +//! let ds = linnerud(); +//! +//! // Fit PLS2 method using 2 principal components (latent variables) +//! let pls = PlsRegression::params(2).fit(&ds)?; +//! +//! // We can either apply the dimension reduction to a dataset +//! let reduced_ds = pls.transform(ds); +//! +//! // ... or predict outputs given a new input sample. +//! let exercices = array![[14., 146., 61.], [6., 80., 60.]]; +//! let physio_measures = pls.predict(exercices); +//! ``` +mod errors; +mod pls_generic; +mod pls_svd; +mod utils; + +use crate::pls_generic::*; + +use linfa::{traits::Fit, traits::PredictRef, traits::Transformer, DatasetBase, Float}; +use ndarray::{Array2, ArrayBase, Data, Ix2}; +use ndarray_linalg::{Lapack, Scalar}; + +pub use errors::*; +pub use pls_svd::*; + +macro_rules! pls_algo { ($name:ident) => { + paste::item! { + + pub struct [](PlsParams); + + impl [] { + /// Set the maximum number of iterations of the power method when algorithm='Nipals'. Ignored otherwise. + pub fn max_iterations(mut self, max_iter: usize) -> Self { + self.0.max_iter = max_iter; + self + } + + /// Set the tolerance used as convergence criteria in the power method: the algorithm + /// stops whenever the squared norm of u_i - u_{i-1} is less than tol, where u corresponds + /// to the left singular vector. + pub fn tolerance(mut self, tolerance: F) -> Self { + self.0.tolerance = tolerance; + self + } + + /// Set whether to scale the dataset + pub fn scale(mut self, scale: bool) -> Self { + self.0.scale = scale; + self + } + + /// Set the algorithm used to estimate the first singular vectors of the cross-covariance matrix. + /// `Nipals` uses the power method while `Svd` will compute the whole SVD. + pub fn algorithm(mut self, algorithm: Algorithm) -> Self { + self.0.algorithm = algorithm; + self + } + } + + pub struct [](Pls); + impl [] { + + pub fn params(n_components: usize) -> [] { + [](Pls::[<$name:lower>](n_components)) + } + + /// Singular vectors of the cross-covariance matrices + pub fn weights(&self) -> (&Array2, &Array2) { + self.0.weights() + } + + /// Loadings of records and targets + pub fn loadings(&self) -> (&Array2, &Array2) { + self.0.loadings() + } + + /// Projection matrices used to transform records and targets + pub fn rotations(&self) -> (&Array2, &Array2) { + self.0.rotations() + } + + /// The coefficients of the linear model such that Y is approximated as Y = X.coefficients + pub fn coefficients(&self) -> &Array2 { + self.0.coefficients() + } + + /// Transform the given dataset in the projected space back to the original space. + pub fn inverse_transform( + &self, + dataset: DatasetBase< + ArrayBase, Ix2>, + ArrayBase, Ix2>, + >, + ) -> DatasetBase, Array2> { + self.0.inverse_transform(dataset) + } + } + + impl> Fit<'_, ArrayBase, ArrayBase> + for [] + { + type Object = Result<[]>; + fn fit( + &self, + dataset: &DatasetBase, ArrayBase>, + ) -> Result<[]> { + let pls = self.0.fit(dataset)?; + Ok([](pls)) + } + } + + impl> Transformer< + DatasetBase, ArrayBase>, + DatasetBase, Array2>, + > for [] + { + /// Apply dimension reduction to the given dataset + fn transform( + &self, + dataset: DatasetBase, ArrayBase>, + ) -> DatasetBase, Array2> { + self.0.transform(dataset) + } + } + + impl> PredictRef, Array2> for [] { + /// Given an input matrix `X`, with shape `(n_samples, n_features)`, + /// `predict` returns the target variable according to [] method + /// learned from the training data distribution. + fn predict_ref<'a>(&'a self, x: &ArrayBase) -> Array2 { + self.0.predict_ref(x) + } + } + } +}} + +pls_algo!(Regression); +pls_algo!(Canonical); +pls_algo!(Cca); + +#[cfg(test)] +mod test { + use super::*; + use approx::assert_abs_diff_eq; + use linfa::{traits::Fit, traits::Predict, traits::Transformer}; + use linfa_datasets::linnerud; + use ndarray::array; + + macro_rules! test_pls_algo { + (Svd) => { + paste::item! { + #[test] + fn []() -> Result<()> { + let ds = linnerud(); + let pls = PlsSvd::::params(3).fit(&ds)?; + let _ds1 = pls.transform(ds); + Ok(()) + } + } + }; + + ($name:ident, $expected:expr) => { + paste::item! { + #[test] + fn []() -> Result<()> { + let ds = linnerud(); + let pls = []::::params(2).fit(&ds)?; + let _ds1 = pls.transform(ds); + let exercices = array![[14., 146., 61.], [6., 80., 60.]]; + let physios = pls.predict(exercices); + assert_abs_diff_eq!($expected, physios.targets(), epsilon=1e-2); + Ok(()) + } + } + }; + } + + // Prediction values were checked against scikit-learn 0.24.1 + test_pls_algo!( + Canonical, + array![ + [180.56979423, 33.29543984, 56.90850758], + [190.854022, 38.91963398, 53.26914489] + ] + ); + test_pls_algo!( + Regression, + array![ + [172.39580643, 34.11919145, 57.15430526], + [192.11167813, 38.05058858, 53.99844922] + ] + ); + test_pls_algo!( + Cca, + array![ + [181.56238421, 34.42502589, 57.31447865], + [205.11767414, 40.23445194, 52.26494323] + ] + ); + test_pls_algo!(Svd); + + #[test] + fn test_one_component_equivalence() -> Result<()> { + // PlsRegression, PlsSvd and PLSCanonical should all be equivalent when n_components is 1 + let ds = linnerud(); + let regression = PlsRegression::params(1).fit(&ds)?.transform(linnerud()); + let canonical = PlsCanonical::params(1).fit(&ds)?.transform(linnerud()); + let svd = PlsSvd::::params(1).fit(&ds)?.transform(linnerud()); + + assert_abs_diff_eq!(regression.records(), canonical.records(), epsilon = 1e-5); + assert_abs_diff_eq!(svd.records(), canonical.records(), epsilon = 1e-5); + Ok(()) + } +} diff --git a/algorithms/linfa-pls/src/pls_generic.rs b/algorithms/linfa-pls/src/pls_generic.rs new file mode 100644 index 000000000..8172f087f --- /dev/null +++ b/algorithms/linfa-pls/src/pls_generic.rs @@ -0,0 +1,834 @@ +use crate::errors::{PlsError, Result}; +use crate::utils; +use linfa::{ + dataset::Records, traits::Fit, traits::PredictRef, traits::Transformer, Dataset, DatasetBase, + Float, +}; +use ndarray::{Array1, Array2, ArrayBase, Data, Ix2}; +use ndarray_linalg::{svd::*, Lapack, Scalar}; +use ndarray_stats::QuantileExt; + +#[cfg_attr( + feature = "serde", + derive(Serialize, Deserialize), + serde(crate = "serde_crate") +)] +#[derive(Debug, Clone)] +pub(crate) struct Pls { + x_mean: Array1, + x_std: Array1, + y_mean: Array1, + y_std: Array1, + x_weights: Array2, // U + y_weights: Array2, // V + x_scores: Array2, // xi + y_scores: Array2, // Omega + x_loadings: Array2, // Gamma + y_loadings: Array2, // Delta + x_rotations: Array2, + y_rotations: Array2, + coefficients: Array2, + n_iters: Array1, +} + +#[derive(PartialEq, Debug, Clone, Copy)] +pub enum Algorithm { + Nipals, + Svd, +} + +#[derive(PartialEq, Debug, Clone, Copy)] +pub(crate) enum DeflationMode { + Regression, + Canonical, +} + +#[derive(PartialEq, Debug, Clone, Copy)] +pub(crate) enum Mode { + A, + B, +} + +/// Generic PLS algorithm. +/// Main ref: Wegelin, a survey of Partial Least Squares (PLS) methods, +/// with emphasis on the two-block case +/// https://www.stat.washington.edu/research/reports/2000/tr371.pdf +impl Pls { + // Constructor for PlsRegression method + pub fn regression(n_components: usize) -> PlsParams { + PlsParams::new(n_components) + } + + // Constructor for PlsCanonical method + pub fn canonical(n_components: usize) -> PlsParams { + PlsParams::new(n_components).deflation_mode(DeflationMode::Canonical) + } + + // Constructor for PlsCca method + pub fn cca(n_components: usize) -> PlsParams { + PlsParams::new(n_components) + .deflation_mode(DeflationMode::Canonical) + .mode(Mode::B) + } + + pub fn weights(&self) -> (&Array2, &Array2) { + (&self.x_weights, &self.y_weights) + } + + #[cfg(test)] + pub fn scores(&self) -> (&Array2, &Array2) { + (&self.x_scores, &self.y_scores) + } + + pub fn loadings(&self) -> (&Array2, &Array2) { + (&self.x_loadings, &self.y_loadings) + } + + pub fn rotations(&self) -> (&Array2, &Array2) { + (&self.x_rotations, &self.y_rotations) + } + + pub fn coefficients(&self) -> &Array2 { + &self.coefficients + } + + pub fn inverse_transform( + &self, + dataset: DatasetBase< + ArrayBase, Ix2>, + ArrayBase, Ix2>, + >, + ) -> DatasetBase, Array2> { + let mut x_orig = dataset.records().dot(&self.x_loadings.t()); + x_orig = &x_orig * &self.x_std; + x_orig = &x_orig + &self.x_mean; + let mut y_orig = dataset.targets().dot(&self.y_loadings.t()); + y_orig = &y_orig * &self.y_std; + y_orig = &y_orig + &self.y_mean; + Dataset::new(x_orig, y_orig) + } +} + +impl> + Transformer< + DatasetBase, ArrayBase>, + DatasetBase, Array2>, + > for Pls +{ + fn transform( + &self, + dataset: DatasetBase, ArrayBase>, + ) -> DatasetBase, Array2> { + let mut x_norm = dataset.records() - &self.x_mean; + x_norm /= &self.x_std; + let mut y_norm = dataset.targets() - &self.y_mean; + y_norm /= &self.y_std; + // Apply rotations + let x_proj = x_norm.dot(&self.x_rotations); + let y_proj = y_norm.dot(&self.y_rotations); + Dataset::new(x_proj, y_proj) + } +} + +impl> PredictRef, Array2> for Pls { + fn predict_ref(&self, x: &ArrayBase) -> Array2 { + let mut x = x - &self.x_mean; + x /= &self.x_std; + x.dot(&self.coefficients) + &self.y_mean + } +} + +#[derive(Debug, Clone)] +pub(crate) struct PlsParams { + pub(crate) n_components: usize, + pub(crate) max_iter: usize, + pub(crate) tolerance: F, + pub(crate) scale: bool, + pub(crate) algorithm: Algorithm, + deflation_mode: DeflationMode, + mode: Mode, +} + +impl PlsParams { + pub fn new(n_components: usize) -> PlsParams { + PlsParams { + n_components, + max_iter: 500, + tolerance: F::from(1e-6).unwrap(), + scale: true, + algorithm: Algorithm::Nipals, + deflation_mode: DeflationMode::Regression, + mode: Mode::A, + } + } + + #[cfg(test)] + pub fn max_iterations(mut self, max_iter: usize) -> Self { + self.max_iter = max_iter; + self + } + + #[cfg(test)] + pub fn tolerance(mut self, tolerance: F) -> Self { + self.tolerance = tolerance; + self + } + + #[cfg(test)] + pub fn scale(mut self, scale: bool) -> Self { + self.scale = scale; + self + } + + #[cfg(test)] + pub fn algorithm(mut self, algorithm: Algorithm) -> Self { + self.algorithm = algorithm; + self + } + + pub fn deflation_mode(mut self, deflation_mode: DeflationMode) -> Self { + self.deflation_mode = deflation_mode; + self + } + + pub fn mode(mut self, mode: Mode) -> Self { + self.mode = mode; + self + } +} + +impl> Fit<'_, ArrayBase, ArrayBase> + for PlsParams +{ + type Object = Result>; + + fn fit(&self, dataset: &DatasetBase, ArrayBase>) -> Result> { + let records = dataset.records(); + let targets = dataset.targets(); + + let n = records.nrows(); + let p = records.ncols(); + let q = targets.ncols(); + + if n < 2 { + return Err(PlsError::NotEnoughSamplesError(format!( + "should be greater than 1, got {}", + dataset.records().nsamples() + ))); + } + + let n_components = self.n_components; + let rank_upper_bound = match self.deflation_mode { + DeflationMode::Regression => { + // With PLSRegression n_components is bounded by the rank of (x.T x) + // see Wegelin page 25 + p + } + DeflationMode::Canonical => { + // With CCA and PLSCanonical, n_components is bounded by the rank of + // X and the rank of Y: see Wegelin page 12 + n.min(p.min(q)) + } + }; + + if 1 > n_components || n_components > rank_upper_bound { + return Err(PlsError::BadComponentNumberError(format!( + "n_components should be in [1, {}], got {}", + rank_upper_bound, n_components + ))); + } + let norm_y_weights = self.deflation_mode == DeflationMode::Canonical; + let (mut xk, mut yk, x_mean, y_mean, x_std, y_std) = + utils::center_scale_dataset(&dataset, self.scale); + + let mut x_weights = Array2::::zeros((p, n_components)); // U + let mut y_weights = Array2::::zeros((q, n_components)); // V + let mut x_scores = Array2::::zeros((n, n_components)); // xi + let mut y_scores = Array2::::zeros((n, n_components)); // Omega + let mut x_loadings = Array2::::zeros((p, n_components)); // Gamma + let mut y_loadings = Array2::::zeros((q, n_components)); // Delta + let mut n_iters = Array1::zeros(n_components); + + // This whole thing corresponds to the algorithm in section 4.1 of the + // review from Wegelin. See above for a notation mapping from code to + // paper. + let eps = F::epsilon(); + for k in 0..n_components { + // Find first left and right singular vectors of the x.T.dot(Y) + // cross-covariance matrix. + + let (mut x_weights_k, mut y_weights_k) = match self.algorithm { + Algorithm::Nipals => { + // Replace columns that are all close to zero with zeros + for mut yj in yk.gencolumns_mut() { + if *(yj.mapv(|y| y.abs()).max().unwrap()) < F::from(10.).unwrap() * eps { + yj.assign(&Array1::zeros(yj.len())); + } + } + + let (x_weights_k, y_weights_k, n_iter) = + self.get_first_singular_vectors_power_method(&xk, &yk, norm_y_weights)?; + n_iters[k] = n_iter; + (x_weights_k, y_weights_k) + } + Algorithm::Svd => self.get_first_singular_vectors_svd(&xk, &yk)?, + }; + utils::svd_flip_1d(&mut x_weights_k, &mut y_weights_k); + + // compute scores, i.e. the projections of x and Y + let x_scores_k = xk.dot(&x_weights_k); + let y_ss = if norm_y_weights { + F::from(1.).unwrap() + } else { + y_weights_k.dot(&y_weights_k) + }; + let y_scores_k = yk.dot(&y_weights_k) / y_ss; + + // Deflation: subtract rank-one approx to obtain xk+1 and yk+1 + let x_loadings_k = x_scores_k.dot(&xk) / x_scores_k.dot(&x_scores_k); + xk = xk - utils::outer(&x_scores_k, &x_loadings_k); // outer product + + let y_loadings_k = match self.deflation_mode { + DeflationMode::Canonical => { + // regress yk on y_score + let y_loadings_k = y_scores_k.dot(&yk) / y_scores_k.dot(&y_scores_k); + yk = yk - utils::outer(&y_scores_k, &y_loadings_k); // outer product + y_loadings_k + } + DeflationMode::Regression => { + // regress yk on x_score + let y_loadings_k = x_scores_k.dot(&yk) / x_scores_k.dot(&x_scores_k); + yk = yk - utils::outer(&x_scores_k, &y_loadings_k); // outer product + y_loadings_k + } + }; + + x_weights.column_mut(k).assign(&x_weights_k); + y_weights.column_mut(k).assign(&y_weights_k); + x_scores.column_mut(k).assign(&x_scores_k); + y_scores.column_mut(k).assign(&y_scores_k); + x_loadings.column_mut(k).assign(&x_loadings_k); + y_loadings.column_mut(k).assign(&y_loadings_k); + } + // x was approximated as xi . Gamma.T + x_(R+1) xi . Gamma.T is a sum + // of n_components rank-1 matrices. x_(R+1) is whatever is left + // to fully reconstruct x, and can be 0 if x is of rank n_components. + // Similiarly, Y was approximated as Omega . Delta.T + Y_(R+1) + + // Compute transformation matrices (rotations_). See User Guide. + let x_rotations = x_weights.dot(&utils::pinv2(&x_loadings.t().dot(&x_weights), None)); + let y_rotations = y_weights.dot(&utils::pinv2(&y_loadings.t().dot(&y_weights), None)); + + let mut coefficients = x_rotations.dot(&y_loadings.t()); + coefficients *= &y_std; + + Ok(Pls { + x_mean, + x_std, + y_mean, + y_std, + x_weights, + y_weights, + x_scores, + y_scores, + x_loadings, + y_loadings, + x_rotations, + y_rotations, + coefficients, + n_iters, + }) + } +} + +impl PlsParams { + /// Return the first left and right singular vectors of x'Y. + /// Provides an alternative to the svd(x'Y) and uses the power method instead. + fn get_first_singular_vectors_power_method( + &self, + x: &ArrayBase, Ix2>, + y: &ArrayBase, Ix2>, + norm_y_weights: bool, + ) -> Result<(Array1, Array1, usize)> { + let eps = F::epsilon(); + + let mut y_score = Array1::ones(y.ncols()); + for col in y.t().genrows() { + if *col.mapv(|v| v.abs()).max().unwrap() > eps { + y_score = col.to_owned(); + break; + } + } + + let mut x_pinv = None; + let mut y_pinv = None; + if self.mode == Mode::B { + x_pinv = Some(utils::pinv2(&x, Some(F::from(10.).unwrap() * eps))); + y_pinv = Some(utils::pinv2(&y, Some(F::from(10.).unwrap() * eps))); + } + + // init to big value for first convergence check + let mut x_weights_old = Array1::::from_elem(x.ncols(), F::from(100.).unwrap()); + + let mut n_iter = 1; + let mut x_weights = Array1::::ones(x.ncols()); + let mut y_weights = Array1::::ones(y.ncols()); + let mut converged = false; + while n_iter < self.max_iter { + x_weights = match self.mode { + Mode::A => x.t().dot(&y_score) / y_score.dot(&y_score), + Mode::B => x_pinv.to_owned().unwrap().dot(&y_score), + }; + x_weights /= (x_weights.dot(&x_weights)).sqrt() + eps; + let x_score = x.dot(&x_weights); + + y_weights = match self.mode { + Mode::A => y.t().dot(&x_score) / x_score.dot(&x_score), + Mode::B => y_pinv.to_owned().unwrap().dot(&x_score), + }; + + if norm_y_weights { + y_weights /= (y_weights.dot(&y_weights)).sqrt() + eps + } + + let ya = y.dot(&y_weights); + let yb = y_weights.dot(&y_weights) + eps; + y_score = ya.mapv(|v| v / yb); + + let x_weights_diff = &x_weights - &x_weights_old; + if x_weights_diff.dot(&x_weights_diff) < self.tolerance || y.ncols() == 1 { + converged = true; + break; + } else { + x_weights_old = x_weights.to_owned(); + n_iter += 1; + } + } + if n_iter == self.max_iter && !converged { + Err(PlsError::PowerMethodNotConvergedError(format!( + "Singular vector computation power method: max iterations ({}) reached", + self.max_iter + ))) + } else { + Ok((x_weights, y_weights, n_iter)) + } + } + + fn get_first_singular_vectors_svd( + &self, + x: &ArrayBase, Ix2>, + y: &ArrayBase, Ix2>, + ) -> Result<(Array1, Array1)> { + let c = x.t().dot(y); + let (u, _, vt) = c.svd(true, true)?; + let u = u.unwrap().column(0).to_owned(); + let vt = vt.unwrap().row(0).to_owned(); + Ok((u, vt)) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_abs_diff_eq; + use linfa::dataset::Records; + use linfa::traits::Predict; + use linfa_datasets::linnerud; + use ndarray::{array, stack, Array, Axis}; + use ndarray_rand::rand::SeedableRng; + use ndarray_rand::rand_distr::StandardNormal; + use ndarray_rand::RandomExt; + use rand_isaac::Isaac64Rng; + + fn assert_matrix_orthonormal(m: &Array2) { + assert_abs_diff_eq!(&m.t().dot(m), &Array::eye(m.ncols()), epsilon = 1e-7); + } + + fn assert_matrix_orthogonal(m: &Array2) { + let k = m.t().dot(m); + assert_abs_diff_eq!(&k, &Array::from_diag(&k.diag()), epsilon = 1e-7); + } + + #[test] + fn test_pls_canonical_basics() -> Result<()> { + // Basic checks for PLSCanonical + let dataset = linnerud(); + let records = dataset.records(); + + let pls = Pls::canonical(records.ncols()).fit(&dataset)?; + + let (x_weights, y_weights) = pls.weights(); + assert_matrix_orthonormal(x_weights); + assert_matrix_orthonormal(y_weights); + + let (x_scores, y_scores) = pls.scores(); + assert_matrix_orthogonal(x_scores); + assert_matrix_orthogonal(y_scores); + + // Check X = TP' and Y = UQ' + let (p, q) = pls.loadings(); + let t = x_scores; + let u = y_scores; + + // Need to scale first + let (xc, yc, ..) = utils::center_scale_dataset(&dataset, true); + assert_abs_diff_eq!(&xc, &t.dot(&p.t()), epsilon = 1e-7); + assert_abs_diff_eq!(&yc, &u.dot(&q.t()), epsilon = 1e-7); + + // Check that rotations on training data lead to scores + let ds = pls.transform(dataset); + assert_abs_diff_eq!(ds.records(), x_scores, epsilon = 1e-7); + assert_abs_diff_eq!(ds.targets(), y_scores, epsilon = 1e-7); + + Ok(()) + } + + #[test] + fn test_sanity_check_pls_regression() { + let dataset = linnerud(); + let pls = Pls::regression(3) + .fit(&dataset) + .expect("PLS fitting failed"); + + // The results were checked against scikit-learn 0.24 PlsRegression + let expected_x_weights = array![ + [0.61330704, -0.00443647, 0.78983213], + [0.74697144, -0.32172099, -0.58183269], + [0.25668686, 0.94682413, -0.19399983] + ]; + + let expected_x_loadings = array![ + [0.61470416, -0.24574278, 0.78983213], + [0.65625755, -0.14396183, -0.58183269], + [0.51733059, 1.00609417, -0.19399983] + ]; + + let expected_y_weights = array![ + [-0.32456184, 0.29892183, 0.20316322], + [-0.42439636, 0.61970543, 0.19320542], + [0.13143144, -0.26348971, -0.17092916] + ]; + + let expected_y_loadings = array![ + [-0.32456184, 0.29892183, 0.20316322], + [-0.42439636, 0.61970543, 0.19320542], + [0.13143144, -0.26348971, -0.17092916] + ]; + assert_abs_diff_eq!(pls.x_weights, expected_x_weights, epsilon = 1e-6); + assert_abs_diff_eq!(pls.x_loadings, expected_x_loadings, epsilon = 1e-6); + assert_abs_diff_eq!(pls.y_weights, expected_y_weights, epsilon = 1e-6); + assert_abs_diff_eq!(pls.y_loadings, expected_y_loadings, epsilon = 1e-6); + } + + #[test] + fn test_sanity_check_pls_regression_constant_column_y() { + let mut dataset = linnerud(); + let nrows = dataset.targets.nrows(); + dataset.targets.column_mut(0).assign(&Array1::ones(nrows)); + let pls = Pls::regression(3) + .fit(&dataset) + .expect("PLS fitting failed"); + + // The results were checked against scikit-learn 0.24 PlsRegression + let expected_x_weights = array![ + [0.6273573, 0.007081799, 0.7786994], + [0.7493417, -0.277612681, -0.6011807], + [0.2119194, 0.960666981, -0.1794690] + ]; + + let expected_x_loadings = array![ + [0.6273512, -0.22464538, 0.7786994], + [0.6643156, -0.09871193, -0.6011807], + [0.5125877, 1.01407380, -0.1794690] + ]; + + let expected_y_loadings = array![ + [0.0000000, 0.0000000, 0.0000000], + [-0.4357300, 0.5828479, 0.2174802], + [0.1353739, -0.2486423, -0.1810386] + ]; + assert_abs_diff_eq!(pls.x_weights, expected_x_weights, epsilon = 1e-6); + assert_abs_diff_eq!(pls.x_loadings, expected_x_loadings, epsilon = 1e-6); + // For the PLSRegression with default parameters, y_loadings == y_weights + assert_abs_diff_eq!(pls.y_loadings, expected_y_loadings, epsilon = 1e-6); + assert_abs_diff_eq!(pls.y_weights, expected_y_loadings, epsilon = 1e-6); + } + + #[test] + fn test_sanity_check_pls_canonical() -> Result<()> { + // Sanity check for PLSCanonical + // The results were checked against the R-package plspm + let dataset = linnerud(); + let pls = Pls::canonical(dataset.records().ncols()).fit(&dataset)?; + + let expected_x_weights = array![ + [-0.61330704, 0.25616119, -0.74715187], + [-0.74697144, 0.11930791, 0.65406368], + [-0.25668686, -0.95924297, -0.11817271] + ]; + + let expected_x_rotations = array![ + [-0.61330704, 0.41591889, -0.62297525], + [-0.74697144, 0.31388326, 0.77368233], + [-0.25668686, -0.89237972, -0.24121788] + ]; + + let expected_y_weights = array![ + [0.58989127, 0.7890047, 0.1717553], + [0.77134053, -0.61351791, 0.16920272], + [-0.23887670, -0.03267062, 0.97050016] + ]; + + let expected_y_rotations = array![ + [0.58989127, 0.7168115, 0.30665872], + [0.77134053, -0.70791757, 0.19786539], + [-0.23887670, -0.00343595, 0.94162826] + ]; + + let (x_weights, y_weights) = pls.weights(); + let (x_rotations, y_rotations) = pls.rotations(); + assert_abs_diff_eq!( + expected_x_rotations.mapv(|v| v.abs()), + x_rotations.mapv(|v| v.abs()), + epsilon = 1e-7 + ); + assert_abs_diff_eq!( + expected_x_weights.mapv(|v| v.abs()), + x_weights.mapv(|v| v.abs()), + epsilon = 1e-7 + ); + assert_abs_diff_eq!( + expected_y_rotations.mapv(|v| v.abs()), + y_rotations.mapv(|v| v.abs()), + epsilon = 1e-7 + ); + assert_abs_diff_eq!( + expected_y_weights.mapv(|v| v.abs()), + y_weights.mapv(|v| v.abs()), + epsilon = 1e-7 + ); + + let x_rotations_sign_flip = (x_rotations / &expected_x_rotations).mapv(|v| v.signum()); + let x_weights_sign_flip = (x_weights / &expected_x_weights).mapv(|v| v.signum()); + let y_rotations_sign_flip = (y_rotations / &expected_y_rotations).mapv(|v| v.signum()); + let y_weights_sign_flip = (y_weights / &expected_y_weights).mapv(|v| v.signum()); + assert_abs_diff_eq!(x_rotations_sign_flip, x_weights_sign_flip); + assert_abs_diff_eq!(y_rotations_sign_flip, y_weights_sign_flip); + + assert_matrix_orthonormal(x_weights); + assert_matrix_orthonormal(y_weights); + + let (x_scores, y_scores) = pls.scores(); + assert_matrix_orthogonal(x_scores); + assert_matrix_orthogonal(y_scores); + Ok(()) + } + + #[test] + fn test_sanity_check_pls_canonical_random() { + // Sanity check for PLSCanonical on random data + // The results were checked against the R-package plspm + let n = 500; + let p_noise = 10; + let q_noise = 5; + + // 2 latents vars: + let mut rng = Isaac64Rng::seed_from_u64(100); + let l1: Array1 = Array1::random_using(n, StandardNormal, &mut rng); + let l2: Array1 = Array1::random_using(n, StandardNormal, &mut rng); + let mut latents = Array::zeros((4, n)); + latents.row_mut(0).assign(&l1); + latents.row_mut(0).assign(&l1); + latents.row_mut(0).assign(&l2); + latents.row_mut(0).assign(&l2); + latents = latents.reversed_axes(); + + let mut x = &latents + &Array2::::random_using((n, 4), StandardNormal, &mut rng); + let mut y = latents + &Array2::::random_using((n, 4), StandardNormal, &mut rng); + + x = stack( + Axis(1), + &[ + x.view(), + Array2::random_using((n, p_noise), StandardNormal, &mut rng).view(), + ], + ) + .unwrap(); + y = stack( + Axis(1), + &[ + y.view(), + Array2::random_using((n, q_noise), StandardNormal, &mut rng).view(), + ], + ) + .unwrap(); + + let ds = Dataset::new(x, y); + let pls = Pls::canonical(3) + .fit(&ds) + .expect("PLS canonical fitting failed"); + + let (x_weights, y_weights) = pls.weights(); + assert_matrix_orthonormal(x_weights); + assert_matrix_orthonormal(y_weights); + + let (x_scores, y_scores) = pls.scores(); + assert_matrix_orthogonal(x_scores); + assert_matrix_orthogonal(y_scores); + } + + #[test] + fn test_scale_and_stability() -> Result<()> { + // scale=True is equivalent to scale=False on centered/scaled data + // This allows to check numerical stability over platforms as well + + let ds = linnerud(); + let (x_s, y_s, ..) = utils::center_scale_dataset(&ds, true); + let ds_s = Dataset::new(x_s, y_s); + + let ds_score = Pls::regression(2) + .scale(true) + .tolerance(1e-3) + .fit(&ds)? + .transform(ds.to_owned()); + let ds_s_score = Pls::regression(2) + .scale(false) + .tolerance(1e-3) + .fit(&ds_s)? + .transform(ds_s.to_owned()); + + assert_abs_diff_eq!(ds_s_score.records(), ds_score.records(), epsilon = 1e-4); + assert_abs_diff_eq!(ds_s_score.targets(), ds_score.targets(), epsilon = 1e-4); + Ok(()) + } + + #[test] + fn test_one_component_equivalence() -> Result<()> { + // PlsRegression, PlsSvd and PLSCanonical should all be equivalent when n_components is 1 + let ds = linnerud(); + let ds2 = linnerud(); + let regression = Pls::regression(1).fit(&ds)?.transform(ds); + let canonical = Pls::canonical(1).fit(&ds2)?.transform(ds2); + + assert_abs_diff_eq!(regression.records(), canonical.records(), epsilon = 1e-7); + Ok(()) + } + + #[test] + fn test_convergence_fail() { + let ds = linnerud(); + assert!( + Pls::canonical(ds.records().nfeatures()) + .max_iterations(2) + .fit(&ds) + .is_err(), + "PLS power method should not converge, hence raise an error" + ); + } + + #[test] + fn test_bad_component_number() { + let ds = linnerud(); + assert!( + Pls::cca(ds.records().nfeatures() + 1).fit(&ds).is_err(), + "n_components too large should raise an error" + ); + assert!( + Pls::canonical(0).fit(&ds).is_err(), + "n_components=0 should raise an error" + ); + } + + #[test] + fn test_singular_value_helpers() -> Result<()> { + // Make sure SVD and power method give approximately the same results + let ds = linnerud(); + + let (mut u1, mut v1, _) = PlsParams::new(2).get_first_singular_vectors_power_method( + ds.records(), + ds.targets(), + true, + )?; + let (mut u2, mut v2) = + PlsParams::new(2).get_first_singular_vectors_svd(ds.records(), ds.targets())?; + + utils::svd_flip_1d(&mut u1, &mut v1); + utils::svd_flip_1d(&mut u2, &mut v2); + + let rtol = 1e-1; + assert_abs_diff_eq!(u1, u2, epsilon = rtol); + assert_abs_diff_eq!(v1, v2, epsilon = rtol); + Ok(()) + } + + macro_rules! test_pls_algo_nipals_svd { + ($($name:ident, )*) => { + paste::item! { + $( + #[test] + fn []() -> Result<()> { + let ds = linnerud(); + let pls = Pls::[<$name>](3).fit(&ds)?; + let ds1 = pls.transform(ds.to_owned()); + let ds2 = Pls::[<$name>](3).algorithm(Algorithm::Svd).fit(&ds)?.transform(ds); + assert_abs_diff_eq!(ds1.records(), ds2.records(), epsilon=1e-2); + let exercices = array![[14., 146., 61.], [6., 80., 60.]]; + let physios = pls.predict(exercices); + println!("Physiologicals = {:?}", physios.targets()); + Ok(()) + } + )* + } + }; + } + + test_pls_algo_nipals_svd! { + canonical, regression, + } + + #[test] + fn test_cca() -> Result<()> { + // values checked against scikit-learn 0.24.1 CCA + let ds = linnerud(); + let cca = Pls::cca(3).fit(&ds)?; + let ds = cca.transform(ds); + let expected_x = array![ + [0.09597886, 0.13862931, -1.0311966], + [-0.7170194, 0.25195026, -0.83049671], + [-0.76492193, 0.37601463, 1.20714686], + [-0.03734329, -0.9746487, 0.79363542], + [0.42809962, -0.50053551, 0.40089685], + [-0.54141144, -0.29403268, -0.47221389], + [-0.29901672, -0.67023009, 0.17945745], + [-0.11425233, -0.43360723, -0.47235823], + [1.29212153, -0.9373391, 0.02572464], + [-0.17770025, 3.4785377, 0.8486413], + [0.39344638, -1.28718499, 1.43816035], + [0.52667844, 0.82080301, -0.02624471], + [0.74616393, 0.54578854, 0.01825073], + [-1.42623443, -0.00884605, -0.24019883], + [-0.72026991, -0.73588273, 0.2241694], + [0.4237932, 0.99977428, -0.1667137], + [-0.88437821, -0.73784626, -0.01073894], + [1.05159992, 0.26381077, -0.83138216], + [1.26196754, -0.18618728, -0.12863494], + [-0.53730151, -0.10896789, -0.92590428] + ]; + assert_abs_diff_eq!(expected_x, ds.records(), epsilon = 1e-2); + Ok(()) + } + + #[test] + fn test_transform_and_inverse() -> Result<()> { + let ds = linnerud(); + let pls = Pls::canonical(3).fit(&ds)?; + + let ds_proj = pls.transform(ds); + let ds_orig = pls.inverse_transform(ds_proj); + + let ds = linnerud(); + assert_abs_diff_eq!(ds.records(), ds_orig.records(), epsilon = 1e-6); + assert_abs_diff_eq!(ds.targets(), ds_orig.targets(), epsilon = 1e-6); + Ok(()) + } +} diff --git a/algorithms/linfa-pls/src/pls_svd.rs b/algorithms/linfa-pls/src/pls_svd.rs new file mode 100644 index 000000000..c1adf558b --- /dev/null +++ b/algorithms/linfa-pls/src/pls_svd.rs @@ -0,0 +1,182 @@ +use crate::errors::{PlsError, Result}; +use crate::utils; +use linfa::{dataset::Records, traits::Fit, traits::Transformer, DatasetBase, Float}; +use ndarray::{s, Array1, Array2, ArrayBase, Data, Ix2}; +use ndarray_linalg::{svd::*, Lapack, Scalar}; + +#[cfg_attr( + feature = "serde", + derive(Serialize, Deserialize), + serde(crate = "serde_crate") +)] +#[derive(Debug, Clone)] +pub struct PlsSvdParams { + n_components: usize, + scale: bool, +} + +impl PlsSvdParams { + pub fn new(n_components: usize) -> PlsSvdParams { + PlsSvdParams { + n_components, + scale: true, + } + } + + pub fn scale(mut self, scale: bool) -> Self { + self.scale = scale; + self + } +} + +impl Default for PlsSvdParams { + fn default() -> Self { + Self::new(2) + } +} + +#[allow(clippy::many_single_char_names)] +impl> Fit<'_, ArrayBase, ArrayBase> + for PlsSvdParams +{ + type Object = Result>; + + fn fit( + &self, + dataset: &DatasetBase, ArrayBase>, + ) -> Result> { + if dataset.nsamples() < 2 { + return Err(PlsError::NotEnoughSamplesError(format!( + "should be greater than 1, got {}", + dataset.records().nsamples() + ))); + } + // we'll compute the SVD of the cross-covariance matrix = X.T.dot(Y) + // This matrix rank is at most min(n_samples, n_features, n_targets) so + // n_components cannot be bigger than that. + + let rank_upper_bound = dataset + .nsamples() + .min(dataset.nfeatures()) + .min(dataset.targets().ncols()); + if 1 > self.n_components || self.n_components > rank_upper_bound { + return Err(PlsError::BadComponentNumberError(format!( + "n_components should be in [1, {}], got {}", + rank_upper_bound, self.n_components + ))); + } + let (x, y, x_mean, y_mean, x_std, y_std) = + utils::center_scale_dataset(&dataset, self.scale); + + // Compute SVD of cross-covariance matrix + let c = x.t().dot(&y); + let (u, _, vt) = c.svd(true, true).unwrap(); + let u = u.unwrap().slice(s![.., ..self.n_components]).to_owned(); + let vt = vt.unwrap().slice(s![..self.n_components, ..]).to_owned(); + let (u, vt) = utils::svd_flip(&u, &vt); + let v = vt.reversed_axes(); + + let x_weights = u; + let y_weights = v; + + Ok(PlsSvd { + x_mean, + x_std, + y_mean, + y_std, + x_weights, + y_weights, + }) + } +} + +pub struct PlsSvd { + x_mean: Array1, + x_std: Array1, + y_mean: Array1, + y_std: Array1, + x_weights: Array2, + y_weights: Array2, +} + +impl PlsSvd { + pub fn params(n_components: usize) -> PlsSvdParams { + PlsSvdParams { + n_components, + scale: true, + } + } + + pub(crate) fn means(&self) -> (&Array1, &Array1) { + (&self.x_mean, &self.y_mean) + } + + pub(crate) fn stds(&self) -> (&Array1, &Array1) { + (&self.x_std, &self.y_std) + } + + pub fn weights(&self) -> (&Array2, &Array2) { + (&self.x_weights, &self.y_weights) + } +} + +impl> + Transformer< + DatasetBase, ArrayBase>, + DatasetBase, Array2>, + > for PlsSvd +{ + fn transform( + &self, + dataset: DatasetBase, ArrayBase>, + ) -> DatasetBase, Array2> { + let (x_mean, y_mean) = &self.means(); + let (x_std, y_std) = &self.stds(); + let (x_weights, y_weights) = &self.weights(); + let xr = (dataset.records() - *x_mean) / *x_std; + let x_scores = xr.dot(*x_weights); + let yr = (dataset.targets() - *y_mean) / *y_std; + let y_scores = yr.dot(*y_weights); + DatasetBase::new(x_scores, y_scores) + } +} + +#[cfg(test)] +mod test { + use super::*; + use approx::assert_abs_diff_eq; + use linfa_datasets::linnerud; + use ndarray::array; + + #[test] + fn test_svd() -> Result<()> { + // values checked against scikit-learn 0.24.1 PlsSVD + let ds = linnerud(); + let pls = PlsSvd::::params(3).fit(&ds)?; + let ds = pls.transform(ds); + let expected_x = array![ + [-0.37144954, -0.0544441, -0.82290137], + [-1.34032497, 0.19638169, -0.71715313], + [-0.08234873, 0.58492788, 0.86557407], + [-0.35496515, -0.62863268, 0.74383396], + [0.46311708, -0.39856773, 0.39748814], + [-1.30584148, -0.20072641, -0.3591439], + [-0.86178968, -0.43791399, 0.2111225], + [-0.79728366, -0.3790222, -0.32195725], + [1.14229739, -0.93000533, 0.19761764], + [3.03443501, 2.81149299, 0.22224139], + [0.40921689, -0.84959246, 1.30923934], + [1.40508381, 0.53658054, -0.09910248], + [1.53073864, 0.29558804, -0.01949986], + [-2.2227316, 0.19806308, -0.2536748], + [-1.49897159, -0.4114628, 0.23494514], + [1.3140941, 0.67110308, -0.2366431], + [-1.88043225, -0.41844445, 0.04307104], + [1.23661961, -0.09041449, -0.63734812], + [1.60595982, -0.37158339, -0.01919568], + [-1.42542371, -0.12332727, -0.73851355] + ]; + assert_abs_diff_eq!(expected_x, ds.records(), epsilon = 1e-6); + Ok(()) + } +} diff --git a/algorithms/linfa-pls/src/utils.rs b/algorithms/linfa-pls/src/utils.rs new file mode 100644 index 000000000..7fd95f53f --- /dev/null +++ b/algorithms/linfa-pls/src/utils.rs @@ -0,0 +1,125 @@ +use linfa::{DatasetBase, Float}; +use ndarray::{s, Array1, Array2, ArrayBase, Axis, Data, DataMut, Ix1, Ix2, Zip}; +use ndarray_linalg::{svd::*, Lapack, Scalar}; +use ndarray_stats::QuantileExt; + +pub fn outer( + a: &ArrayBase, Ix1>, + b: &ArrayBase, Ix1>, +) -> Array2 { + let mut outer = Array2::zeros((a.len(), b.len())); + Zip::from(outer.genrows_mut()).and(a).apply(|mut out, ai| { + out.assign(&b.mapv(|v| *ai * v)); + }); + outer +} + +pub fn pinv2>( + x: &ArrayBase, + cond: Option, +) -> Array2 { + let (opt_u, s, opt_vh) = x.svd(true, true).unwrap(); + let u = opt_u.unwrap(); + let vh = opt_vh.unwrap(); + + let cond = cond.unwrap_or( + F::from(*s.max().unwrap()).unwrap() + * F::from(x.nrows().max(x.ncols())).unwrap() + * F::epsilon(), + ); + + let rank = s.fold(0, |mut acc, v| { + if F::from(*v).unwrap() > cond { + acc += 1 + }; + acc + }); + + let mut ucut = u.slice(s![.., ..rank]).to_owned(); + ucut /= &s.slice(s![..rank]).mapv(|v| F::from(v).unwrap()); + ucut.dot(&vh.slice(s![..rank, ..])) + .mapv(|v| v.conj()) + .t() + .to_owned() +} + +#[allow(clippy::type_complexity)] +pub fn center_scale_dataset>( + dataset: &DatasetBase, ArrayBase>, + scale: bool, +) -> ( + Array2, + Array2, + Array1, + Array1, + Array1, + Array1, +) { + let (xnorm, x_mean, x_std) = center_scale(&dataset.records(), scale); + let (ynorm, y_mean, y_std) = center_scale(&dataset.targets(), scale); + (xnorm, ynorm, x_mean, y_mean, x_std, y_std) +} + +fn center_scale( + x: &ArrayBase, Ix2>, + scale: bool, +) -> (Array2, Array1, Array1) { + let x_mean = x.mean_axis(Axis(0)).unwrap(); + let (xnorm, x_std) = if scale { + let mut x_std = x.std_axis(Axis(0), F::one()); + x_std.mapv_inplace(|v| if v == F::zero() { F::one() } else { v }); + ((x - &x_mean) / &x_std, x_std) + } else { + ((x - &x_mean), Array1::ones(x.ncols())) + }; + + (xnorm, x_mean, x_std) +} + +pub fn svd_flip_1d( + x_weights: &mut ArrayBase, Ix1>, + y_weights: &mut ArrayBase, Ix1>, +) { + let biggest_abs_val_idx = x_weights.mapv(|v| v.abs()).argmax().unwrap(); + let sign: F = x_weights[biggest_abs_val_idx].signum(); + x_weights.map_inplace(|v| *v *= sign); + y_weights.map_inplace(|v| *v *= sign); +} + +pub fn svd_flip( + u: &ArrayBase, Ix2>, + v: &ArrayBase, Ix2>, +) -> (Array2, Array2) { + // columns of u, rows of v + let abs_u = u.mapv(|v| v.abs()); + let max_abs_val_indices = abs_u.map_axis(Axis(0), |col| col.argmax().unwrap()); + let mut signs = Array1::::zeros(u.ncols()); + let range: Vec = (0..u.ncols()).collect(); + Zip::from(&mut signs) + .and(&max_abs_val_indices) + .and(&range) + .apply(|s, &i, &j| *s = u[[i, j]].signum()); + (u * &signs, v * &signs.insert_axis(Axis(1))) +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_abs_diff_eq; + use ndarray::array; + + #[test] + fn test_outer() { + let a = array![1., 2., 3.]; + let b = array![2., 3.]; + let expected = array![[2., 3.], [4., 6.], [6., 9.]]; + assert_abs_diff_eq!(expected, outer(&a, &b)); + } + + #[test] + fn test_pinv2() { + let a = array![[1., 2., 3.], [4., 5., 6.], [7., 8., 10.]]; + let a_pinv2 = pinv2(&a, None); + assert_abs_diff_eq!(a.dot(&a_pinv2), Array2::eye(3), epsilon = 1e-6) + } +} diff --git a/docs/website/content/about.md b/docs/website/content/about.md index 62917a5e3..c89490a17 100644 --- a/docs/website/content/about.md +++ b/docs/website/content/about.md @@ -29,6 +29,7 @@ Where does `linfa` stand right now? [Are we learning yet?](http://www.arewelearn | [hierarchical](https://github.com/rust-ml/linfa/tree/master/algorithms/linfa-hierarchical) | Agglomerative hierarchical clustering | Tested | Unsupervised learning | Cluster and build hierarchy of clusters | | [bayes](https://github.com/rust-ml/linfa/tree/master/algorithms/linfa-bayes) | Naive Bayes | Tested | Supervised learning | Contains Gaussian Naive Bayes | | [ica](https://github.com/rust-ml/linfa/tree/master/algorithms/linfa-ica) | Independent component analysis | Tested | Unsupervised learning | Contains FastICA implementation | +| [pls](https://github.com/rust-ml/linfa/tree/master/algorithms/linfa-pls) | Partial Least Squares | Tested | Supervised learning | Contains PLS estimators for dimensionality reduction and regression |