WIP: piv implementation with tests

emiddleton · emiddleton · commit 970232c68e59 · 2021-06-01T17:44:59.000+09:00
diff --git 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
@@ -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!()
+        }
+    }
+}
diff --git a/ndarray-linalg/src/rank.rs 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<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)
+    }
+}
diff --git a/ndarray-linalg/tests/pinv.rs 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<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);
+}
diff --git a/ndarray-linalg/tests/rank.rs 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<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());
+}
