WIP: piv implementation with tests

Author: emiddleton

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::*;

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<Self::E>) -> Result<Self::C>;
+}
+
+/// pseudo-inverse
+pub trait PInvInto {
+ type E;
+ type C;
+ fn pinv_into(self, rcond: Option<Self::E>) -> Result<Self::C>;
+}
+
+/// pseudo-inverse for a mutable reference of a matrix
+pub trait PInvInplace {
+ type E;
+ type C;
+ fn pinv_inplace(&mut self, rcond: Option<Self::E>) -> Result<Self::C>;
+}
+
+impl<A, S> PInvInto for ArrayBase<S, Ix2>
+where
+ A: Scalar + Lapack,
+ S: DataMut<Elem = A>,
+{
+ type E = A::Real;
+ type C = Array2<A>;
+
+ fn pinv_into(mut self, rcond: Option<Self::E>) -> Result<Self::C> {
+ self.pinv_inplace(rcond)
+ }
+}
+
+impl<A, S> Pinv for ArrayBase<S, Ix2>
+where
+ A: Scalar + Lapack,
+ S: Data<Elem = A>,
+{
+ type E = A::Real;
+ type C = Array2<A>;
+
+ fn pinv(&self, rcond: Option<Self::E>) -> Result<Self::C> {
+ let a = self.to_owned();
+ a.pinv_into(rcond)
+ }
+}
+
+impl<A, S> PInvInplace for ArrayBase<S, Ix2>
+where
+ A: Scalar + Lapack,
+ S: DataMut<Elem = A>,
+{
+ type E = A::Real;
+ type C = Array2<A>;
+
+ fn pinv_inplace(&mut self, rcond: Option<Self::E>) -> Result<Self::C> {
+ 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!()
+ }
+ }
+}

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<Ix>;
+}
+
+impl<A, S> Rank for ArrayBase<S, Ix2>
+where
+ A: Scalar + Lapack,
+ S: Data<Elem = A>,
+{
+ fn rank(&self) -> Result<Ix> {
+ 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)
+ }
+}

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<A, S, Sh, D>(sh: Sh) -> ArrayBase<S, D>
+where
+ A: Scalar,
+ S: DataOwned<Elem = A>,
+ D: Dimension,
+ Sh: ShapeBuilder<Dim = D>,
+{
+ ArrayBase::zeros(sh)
+}
+
+/// creates a random matrix and repeatedly creates a linear dependency between rows until the
+/// rank drops.
+pub fn partial_rank<A, Sh>(sh: Sh) -> Array2<A>
+where
+ A: Scalar + Lapack,
+ Sh: ShapeBuilder<Dim = Ix2>,
+{
+ let mut rng = thread_rng();
+ let mut result: Array2<A> = 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::<Vec<usize>>();
+ 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<A, Sh>(sh: Sh) -> Array2<A>
+where
+ A: Scalar + Lapack,
+ Sh: ShapeBuilder<Dim = Ix2> + Clone,
+{
+ for _ in 0..10 {
+ let r: Array2<A> = 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<T: Scalar + Lapack>(a: &Array2<T>, 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 [<pinv_test_ $type _ $test _ $n x $m _r>]() {
+ let a: Array2<$type> = $test(($n, $m));
+ test::<$type>(&a, $t);
+ }
+
+ #[test]
+ fn [<pinv_test_ $type _ $test _ $n x $m _c>]() {
+ 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<f64> = arr2(&[
+ [ 1., -1., 0.],
+ [-1., 2., -1.],
+ [ 0., -1., 1.]
+ ],
+ );
+ #[rustfmt::skip]
+ let a_plus_actual: Array2<f64> = 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);
+}

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<f64> = 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<f64> = 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<f64> = arr2(&[
+ [1., 0., 2.],
+ [2., 1., 0.],
+ [3., 2., 1.],
+ ],
+ );
+ assert_eq!(3, a.rank().unwrap());
+}