From 970232c68e59ac569080fb9529bf9c22719c6ae6 Mon Sep 17 00:00:00 2001 From: Edward Middleton Date: Tue, 1 Jun 2021 17:44:59 +0900 Subject: [PATCH 1/3] WIP: piv implementation with tests --- ndarray-linalg/src/lib.rs | 4 + ndarray-linalg/src/pinv.rs | 116 ++++++++++++++++++++++++ ndarray-linalg/src/rank.rs | 27 ++++++ ndarray-linalg/tests/pinv.rs | 167 +++++++++++++++++++++++++++++++++++ ndarray-linalg/tests/rank.rs | 38 ++++++++ 5 files changed, 352 insertions(+) create mode 100644 ndarray-linalg/src/pinv.rs create mode 100644 ndarray-linalg/src/rank.rs create mode 100644 ndarray-linalg/tests/pinv.rs create mode 100644 ndarray-linalg/tests/rank.rs diff --git a/ndarray-linalg/src/lib.rs b/ndarray-linalg/src/lib.rs index ba9d10c5..aef7f3cb 100644 --- a/ndarray-linalg/src/lib.rs +++ b/ndarray-linalg/src/lib.rs @@ -64,7 +64,9 @@ pub mod lobpcg; pub mod norm; pub mod operator; pub mod opnorm; +pub mod pinv; pub mod qr; +pub mod rank; pub mod solve; pub mod solveh; pub mod svd; @@ -88,7 +90,9 @@ pub use lobpcg::{TruncatedEig, TruncatedOrder, TruncatedSvd}; pub use norm::*; pub use operator::*; pub use opnorm::*; +pub use pinv::*; pub use qr::*; +pub use rank::*; pub use solve::*; pub use solveh::*; pub use svd::*; diff --git a/ndarray-linalg/src/pinv.rs b/ndarray-linalg/src/pinv.rs new file mode 100644 index 00000000..9f5a36a2 --- /dev/null +++ b/ndarray-linalg/src/pinv.rs @@ -0,0 +1,116 @@ +//! Moore-Penrose pseudo-inverse of a Matrices +//! +//! [](https://hadrienj.github.io/posts/Deep-Learning-Book-Series-2.9-The-Moore-Penrose-Pseudoinverse/) + +use crate::{error::*, svd::SVDInplace, types::*}; +use ndarray::*; +use num_traits::Float; + +/// pseudo-inverse of a matrix reference +pub trait Pinv { + type E; + type C; + fn pinv(&self, threshold: Option) -> Result; +} + +/// pseudo-inverse +pub trait PInvInto { + type E; + type C; + fn pinv_into(self, rcond: Option) -> Result; +} + +/// pseudo-inverse for a mutable reference of a matrix +pub trait PInvInplace { + type E; + type C; + fn pinv_inplace(&mut self, rcond: Option) -> Result; +} + +impl PInvInto for ArrayBase +where + A: Scalar + Lapack, + S: DataMut, +{ + type E = A::Real; + type C = Array2; + + fn pinv_into(mut self, rcond: Option) -> Result { + self.pinv_inplace(rcond) + } +} + +impl Pinv for ArrayBase +where + A: Scalar + Lapack, + S: Data, +{ + type E = A::Real; + type C = Array2; + + fn pinv(&self, rcond: Option) -> Result { + let a = self.to_owned(); + a.pinv_into(rcond) + } +} + +impl PInvInplace for ArrayBase +where + A: Scalar + Lapack, + S: DataMut, +{ + type E = A::Real; + type C = Array2; + + fn pinv_inplace(&mut self, rcond: Option) -> Result { + if let (Some(u), s, Some(v_h)) = self.svd_inplace(true, true)? { + // threshold = ε⋅max(m, n)⋅max(Σ) + // NumPy defaults rcond to 1e-15 which is about 10 * f64 machine epsilon + let rcond = rcond.unwrap_or_else(|| { + let (n, m) = self.dim(); + Self::E::epsilon() * Self::E::real(n.max(m)) + }); + let threshold = rcond * s[0]; + + // Determine how many singular values to keep and compute the + // values of `V Σ⁺` (up to `num_keep` columns). + let (num_keep, v_s_inv) = { + let mut v_h_t = v_h.reversed_axes(); + let mut num_keep = 0; + for (&sing_val, mut v_h_t_col) in s.iter().zip(v_h_t.columns_mut()) { + if sing_val > threshold { + let sing_val_recip = sing_val.recip(); + v_h_t_col.map_inplace(|v_h_t| { + *v_h_t = A::from_real(sing_val_recip) * v_h_t.conj() + }); + num_keep += 1; + } else { + /* + if sing_val != Self::E::real(0.0) { + panic!( + "for {:#?} singular value {:?} smaller then threshold {:?}", + &self, &sing_val, &threshold + ); + } + */ + break; + } + } + v_h_t.slice_axis_inplace(Axis(1), Slice::from(..num_keep)); + (num_keep, v_h_t) + }; + + // Compute `U^H` (up to `num_keep` rows). + let u_h = { + let mut u_t = u.reversed_axes(); + u_t.slice_axis_inplace(Axis(0), Slice::from(..num_keep)); + u_t.map_inplace(|x| *x = x.conj()); + u_t + }; + + Ok(v_s_inv.dot(&u_h)) + } else { + unreachable!() + } + } +} diff --git a/ndarray-linalg/src/rank.rs b/ndarray-linalg/src/rank.rs new file mode 100644 index 00000000..759960de --- /dev/null +++ b/ndarray-linalg/src/rank.rs @@ -0,0 +1,27 @@ +///! Computes the rank of a matrix using single value decomposition +use ndarray::*; + +use super::error::*; +use super::svd::SVD; +use super::types::*; +use num_traits::Float; + +pub trait Rank { + fn rank(&self) -> Result; +} + +impl Rank for ArrayBase +where + A: Scalar + Lapack, + S: Data, +{ + fn rank(&self) -> Result { + let (_, sv, _) = self.svd(false, false)?; + + let (n, m) = self.dim(); + let tol = A::Real::epsilon() * A::Real::real(n.max(m)) * sv[0]; + + let output = sv.iter().take_while(|v| v > &&tol).count(); + Ok(output) + } +} diff --git a/ndarray-linalg/tests/pinv.rs b/ndarray-linalg/tests/pinv.rs new file mode 100644 index 00000000..43088e02 --- /dev/null +++ b/ndarray-linalg/tests/pinv.rs @@ -0,0 +1,167 @@ +use ndarray::arr2; +use ndarray::*; +use ndarray_linalg::rank::Rank; +use ndarray_linalg::*; +use rand::{seq::SliceRandom, thread_rng}; + +/// creates a zero matrix which always has rank zero +pub fn zero_rank(sh: Sh) -> ArrayBase +where + A: Scalar, + S: DataOwned, + D: Dimension, + Sh: ShapeBuilder, +{ + ArrayBase::zeros(sh) +} + +/// creates a random matrix and repeatedly creates a linear dependency between rows until the +/// rank drops. +pub fn partial_rank(sh: Sh) -> Array2 +where + A: Scalar + Lapack, + Sh: ShapeBuilder, +{ + let mut rng = thread_rng(); + let mut result: Array2 = random(sh); + println!("before: {:?}", result); + + let (n, m) = result.dim(); + println!("(n, m) => ({:?},{:?})", n, m); + + // create randomized row iterator + let min_dim = n.min(m); + let mut row_indexes = (0..min_dim).into_iter().collect::>(); + row_indexes.as_mut_slice().shuffle(&mut rng); + let mut row_index_iter = row_indexes.iter().cycle(); + + for count in 1..=10 { + println!("count: {}", count); + let (&x, &y) = ( + row_index_iter.next().unwrap(), + row_index_iter.next().unwrap(), + ); + let (from_row_index, to_row_index) = if x < y { (x, y) } else { (y, x) }; + println!("(r_f, r_t) => ({:?},{:?})", from_row_index, to_row_index); + + let mut it = result.outer_iter_mut(); + let from_row = it.nth(from_row_index).unwrap(); + let mut to_row = it.nth(to_row_index - (from_row_index + 1)).unwrap(); + + // set the to_row with the value of the from_row multiplied by rand_multiple + let rand_multiple = A::rand(&mut rng); + println!("rand_multiple: {:?}", rand_multiple); + Zip::from(&mut to_row) + .and(&from_row) + .for_each(|r1, r2| *r1 = *r2 * rand_multiple); + + if let Ok(rank) = result.rank() { + println!("result: {:?}", result); + println!("rank: {:?}", rank); + if rank > 0 && rank < min_dim { + return result; + } + } + } + unreachable!("unable to generate random partial rank matrix after making 10 mutations") +} + +/// creates a random matrix and insures it is full rank. +pub fn full_rank(sh: Sh) -> Array2 +where + A: Scalar + Lapack, + Sh: ShapeBuilder + Clone, +{ + for _ in 0..10 { + let r: Array2 = random(sh.clone()); + let (n, m) = r.dim(); + let n = n.min(m); + if let Ok(rank) = r.rank() { + println!("result: {:?}", r); + println!("rank: {:?}", rank); + if rank == n { + return r; + } + } + } + unreachable!("unable to generate random full rank matrix in 10 tries") +} + +fn test(a: &Array2, tolerance: T::Real) { + println!("a = \n{:?}", &a); + let a_plus: Array2<_> = a.pinv(None).unwrap(); + println!("a_plus = \n{:?}", &a_plus); + let ident = a.dot(&a_plus); + assert_close_l2!(&ident.dot(a), &a, tolerance); + assert_close_l2!(&a_plus.dot(&ident), &a_plus, tolerance); +} + +macro_rules! test_both_impl { + ($type:ty, $test:tt, $n:expr, $m:expr, $t:expr) => { + paste::item! { + #[test] + fn []() { + let a: Array2<$type> = $test(($n, $m)); + test::<$type>(&a, $t); + } + + #[test] + fn []() { + let a = $test(($n, $m).f()); + test::<$type>(&a, $t); + } + } + }; +} + +macro_rules! test_pinv_impl { + ($type:ty, $n:expr, $m:expr, $a:expr) => { + test_both_impl!($type, zero_rank, $n, $m, $a); + test_both_impl!($type, partial_rank, $n, $m, $a); + test_both_impl!($type, full_rank, $n, $m, $a); + }; +} + +test_pinv_impl!(f32, 3, 3, 1e-4); +test_pinv_impl!(f32, 4, 3, 1e-4); +test_pinv_impl!(f32, 3, 4, 1e-4); + +test_pinv_impl!(c32, 3, 3, 1e-4); +test_pinv_impl!(c32, 4, 3, 1e-4); +test_pinv_impl!(c32, 3, 4, 1e-4); + +test_pinv_impl!(f64, 3, 3, 1e-12); +test_pinv_impl!(f64, 4, 3, 1e-12); +test_pinv_impl!(f64, 3, 4, 1e-12); + +test_pinv_impl!(c64, 3, 3, 1e-12); +test_pinv_impl!(c64, 4, 3, 1e-12); +test_pinv_impl!(c64, 3, 4, 1e-12); + +// +// This matrix was taken from 7.1.1 Test1 in +// "On Moore-Penrose Pseudoinverse Computation for Stiffness Matrices Resulting +// from Higher Order Approximation" by Marek Klimczak +// https://doi.org/10.1155/2019/5060397 +// +#[test] +fn pinv_test_single_value_less_then_threshold_3x3() { + #[rustfmt::skip] + let a: Array2 = arr2(&[ + [ 1., -1., 0.], + [-1., 2., -1.], + [ 0., -1., 1.] + ], + ); + #[rustfmt::skip] + let a_plus_actual: Array2 = arr2(&[ + [ 5. / 9., -1. / 9., -4. / 9.], + [-1. / 9., 2. / 9., -1. / 9.], + [-4. / 9., -1. / 9., 5. / 9.], + ], + ); + let a_plus: Array2<_> = a.pinv(None).unwrap(); + println!("a_plus -> {:?}", &a_plus); + println!("a_plus_actual -> {:?}", &a_plus); + assert_close_l2!(&a_plus, &a_plus_actual, 1e-15); +} diff --git a/ndarray-linalg/tests/rank.rs b/ndarray-linalg/tests/rank.rs new file mode 100644 index 00000000..3cb93d85 --- /dev/null +++ b/ndarray-linalg/tests/rank.rs @@ -0,0 +1,38 @@ +use ndarray::*; +use ndarray_linalg::*; + +#[test] +fn rank_test_zero_3x3() { + #[rustfmt::skip] + let a: Array2 = arr2(&[ + [0., 0., 0.], + [0., 0., 0.], + [0., 0., 0.], + ], + ); + assert_eq!(0, a.rank().unwrap()); +} + +#[test] +fn rank_test_partial_3x3() { + #[rustfmt::skip] + let a: Array2 = arr2(&[ + [1., 2., 3.], + [4., 5., 6.], + [7., 8., 9.], + ], + ); + assert_eq!(2, a.rank().unwrap()); +} + +#[test] +fn rank_test_full_3x3() { + #[rustfmt::skip] + let a: Array2 = arr2(&[ + [1., 0., 2.], + [2., 1., 0.], + [3., 2., 1.], + ], + ); + assert_eq!(3, a.rank().unwrap()); +} From ff89e2b982b8ad042e8bd1b1b749668afebed93e Mon Sep 17 00:00:00 2001 From: Edward Middleton Date: Wed, 2 Jun 2021 13:40:22 +0900 Subject: [PATCH 2/3] added generate::random_with_rank function and updated test to use it --- ndarray-linalg/src/generate.rs | 50 +++++++++++++++++++++ ndarray-linalg/tests/pinv.rs | 82 +++++++--------------------------- 2 files changed, 66 insertions(+), 66 deletions(-) diff --git a/ndarray-linalg/src/generate.rs b/ndarray-linalg/src/generate.rs index 67def14a..6bffa5d9 100644 --- a/ndarray-linalg/src/generate.rs +++ b/ndarray-linalg/src/generate.rs @@ -1,12 +1,15 @@ //! Generator functions for matrices +use ndarray::linalg::general_mat_mul; use ndarray::*; use rand::prelude::*; use super::convert::*; use super::error::*; use super::qr::*; +use super::rank::Rank; use super::types::*; +use super::Scalar; /// Hermite conjugate matrix pub fn conjugate(a: &ArrayBase) -> ArrayBase @@ -34,6 +37,53 @@ where ArrayBase::from_shape_fn(sh, |_| A::rand(&mut rng)) } +/// Generate random array with a given rank +/// +/// The rank must be less then or equal to the smallest dimension of array +pub fn random_with_rank(shape: Sh, rank: usize) -> Array2 +where + A: Scalar + Lapack, + Sh: ShapeBuilder + Clone, +{ + // handle zero-rank case + if rank == 0 { + return Array2::zeros(shape); + } + + let (n, m) = shape.clone().into_shape().raw_dim().into_pattern(); + let min_dim = usize::min(n, m); + assert!(rank <= min_dim); + + for _ in 0..10 { + // handle full-rank case + let out = if rank == min_dim { + random(shape.clone()) + + // handle partial-rank case + } else { + // multiplying two full-rank arrays with dimensions `m × r` and `r × n` will + // produce `an m × n` array with rank `r` + // https://en.wikipedia.org/wiki/Rank_(linear_algebra)#Properties + let mut out = Array2::zeros(shape.clone()); + let left: Array2 = random([out.nrows(), rank]); + let right: Array2 = random([rank, out.ncols()]); + general_mat_mul(A::one(), &left, &right, A::zero(), &mut out); + out + }; + + // check rank + if let Ok(out_rank) = out.rank() { + if out_rank == rank { + return out; + } + } + } + + unreachable!( + "Failed to generate random matrix of desired rank within 10 tries. This is very unlikely." + ); +} + /// Generate random unitary matrix using QR decomposition /// /// Be sure that this it **NOT** a uniform distribution. Use it only for test purpose. diff --git a/ndarray-linalg/tests/pinv.rs b/ndarray-linalg/tests/pinv.rs index 43088e02..4829b137 100644 --- a/ndarray-linalg/tests/pinv.rs +++ b/ndarray-linalg/tests/pinv.rs @@ -1,90 +1,40 @@ use ndarray::arr2; use ndarray::*; -use ndarray_linalg::rank::Rank; use ndarray_linalg::*; -use rand::{seq::SliceRandom, thread_rng}; +use rand::{thread_rng, Rng}; -/// creates a zero matrix which always has rank zero -pub fn zero_rank(sh: Sh) -> ArrayBase +/// create a zero rank array +pub fn zero_rank(sh: Sh) -> Array2 where - A: Scalar, - S: DataOwned, - D: Dimension, - Sh: ShapeBuilder, + A: Scalar + Lapack, + Sh: ShapeBuilder + Clone, { - ArrayBase::zeros(sh) + random_with_rank(sh, 0) } -/// creates a random matrix and repeatedly creates a linear dependency between rows until the -/// rank drops. +/// create a random matrix with a random partial rank. pub fn partial_rank(sh: Sh) -> Array2 where A: Scalar + Lapack, - Sh: ShapeBuilder, + Sh: ShapeBuilder + Clone, { let mut rng = thread_rng(); - let mut result: Array2 = random(sh); - println!("before: {:?}", result); - - let (n, m) = result.dim(); - println!("(n, m) => ({:?},{:?})", n, m); - - // create randomized row iterator + let (m, n) = sh.clone().into_shape().raw_dim().into_pattern(); let min_dim = n.min(m); - let mut row_indexes = (0..min_dim).into_iter().collect::>(); - row_indexes.as_mut_slice().shuffle(&mut rng); - let mut row_index_iter = row_indexes.iter().cycle(); - - for count in 1..=10 { - println!("count: {}", count); - let (&x, &y) = ( - row_index_iter.next().unwrap(), - row_index_iter.next().unwrap(), - ); - let (from_row_index, to_row_index) = if x < y { (x, y) } else { (y, x) }; - println!("(r_f, r_t) => ({:?},{:?})", from_row_index, to_row_index); - - let mut it = result.outer_iter_mut(); - let from_row = it.nth(from_row_index).unwrap(); - let mut to_row = it.nth(to_row_index - (from_row_index + 1)).unwrap(); - - // set the to_row with the value of the from_row multiplied by rand_multiple - let rand_multiple = A::rand(&mut rng); - println!("rand_multiple: {:?}", rand_multiple); - Zip::from(&mut to_row) - .and(&from_row) - .for_each(|r1, r2| *r1 = *r2 * rand_multiple); - - if let Ok(rank) = result.rank() { - println!("result: {:?}", result); - println!("rank: {:?}", rank); - if rank > 0 && rank < min_dim { - return result; - } - } - } - unreachable!("unable to generate random partial rank matrix after making 10 mutations") + let rank = rng.gen_range(1..min_dim); + println!("desired rank = {}", rank); + random_with_rank(sh, rank) } -/// creates a random matrix and insures it is full rank. +/// create a random matrix and ensures it is full rank. pub fn full_rank(sh: Sh) -> Array2 where A: Scalar + Lapack, Sh: ShapeBuilder + Clone, { - for _ in 0..10 { - let r: Array2 = random(sh.clone()); - let (n, m) = r.dim(); - let n = n.min(m); - if let Ok(rank) = r.rank() { - println!("result: {:?}", r); - println!("rank: {:?}", rank); - if rank == n { - return r; - } - } - } - unreachable!("unable to generate random full rank matrix in 10 tries") + let (m, n) = sh.clone().into_shape().raw_dim().into_pattern(); + let min_dim = n.min(m); + random_with_rank(sh, min_dim) } fn test(a: &Array2, tolerance: T::Real) { From 942b7d741c530099e896d47dea86fade7e17f968 Mon Sep 17 00:00:00 2001 From: Edward Middleton Date: Wed, 2 Jun 2021 14:40:39 +0900 Subject: [PATCH 3/3] avoid reallocating arrays in loop --- ndarray-linalg/src/generate.rs | 54 ++++++++++++++++++++++------------ 1 file changed, 35 insertions(+), 19 deletions(-) diff --git a/ndarray-linalg/src/generate.rs b/ndarray-linalg/src/generate.rs index 6bffa5d9..0ee2a07a 100644 --- a/ndarray-linalg/src/generate.rs +++ b/ndarray-linalg/src/generate.rs @@ -54,28 +54,44 @@ where let min_dim = usize::min(n, m); assert!(rank <= min_dim); - for _ in 0..10 { - // handle full-rank case - let out = if rank == min_dim { - random(shape.clone()) - - // handle partial-rank case - } else { - // multiplying two full-rank arrays with dimensions `m × r` and `r × n` will - // produce `an m × n` array with rank `r` - // https://en.wikipedia.org/wiki/Rank_(linear_algebra)#Properties - let mut out = Array2::zeros(shape.clone()); - let left: Array2 = random([out.nrows(), rank]); - let right: Array2 = random([rank, out.ncols()]); + let mut rng = thread_rng(); + + // handle full-rank case + if rank == min_dim { + let mut out = random(shape); + for _ in 0..10 { + // check rank + if let Ok(out_rank) = out.rank() { + if out_rank == rank { + return out; + } + } + + out.mapv_inplace(|_| A::rand(&mut rng)); + } + + // handle partial-rank case + // + // multiplying two full-rank arrays with dimensions `m × r` and `r × n` will + // produce `an m × n` array with rank `r` + // https://en.wikipedia.org/wiki/Rank_(linear_algebra)#Properties + } else { + let mut out = Array2::zeros(shape); + let mut left: Array2 = random([out.nrows(), rank]); + let mut right: Array2 = random([rank, out.ncols()]); + + for _ in 0..10 { general_mat_mul(A::one(), &left, &right, A::zero(), &mut out); - out - }; - // check rank - if let Ok(out_rank) = out.rank() { - if out_rank == rank { - return out; + // check rank + if let Ok(out_rank) = out.rank() { + if out_rank == rank { + return out; + } } + + left.mapv_inplace(|_| A::rand(&mut rng)); + right.mapv_inplace(|_| A::rand(&mut rng)); } }