diff --git a/ndarray-linalg/src/generate.rs b/ndarray-linalg/src/generate.rs
index 67def14a..0ee2a07a 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, Si, So>(a: &ArrayBase<Si, Ix2>) -> ArrayBase<So, Ix2>
@@ -34,6 +37,69 @@ 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<A, Sh>(shape: Sh, rank: usize) -> Array2<A>
+where
+ A: Scalar + Lapack,
+ Sh: ShapeBuilder<Dim = Ix2> + 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);
+
+ 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<A> = random([out.nrows(), rank]);
+ let mut right: Array2<A> = random([rank, out.ncols()]);
+
+ for _ in 0..10 {
+ general_mat_mul(A::one(), &left, &right, A::zero(), &mut 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));
+ }
+ }
+
+ 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/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<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
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<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
new file mode 100644
index 00000000..4829b137
--- /dev/null
+++ b/ndarray-linalg/tests/pinv.rs
@@ -0,0 +1,117 @@
+use ndarray::arr2;
+use ndarray::*;
+use ndarray_linalg::*;
+use rand::{thread_rng, Rng};
+
+/// create a zero rank array
+pub fn zero_rank<A, Sh>(sh: Sh) -> Array2<A>
+where
+ A: Scalar + Lapack,
+ Sh: ShapeBuilder<Dim = Ix2> + Clone,
+{
+ random_with_rank(sh, 0)
+}
+
+/// create a random matrix with a random partial rank.
+pub fn partial_rank<A, Sh>(sh: Sh) -> Array2<A>
+where
+ A: Scalar + Lapack,
+ Sh: ShapeBuilder<Dim = Ix2> + Clone,
+{
+ let mut rng = thread_rng();
+ let (m, n) = sh.clone().into_shape().raw_dim().into_pattern();
+ let min_dim = n.min(m);
+ let rank = rng.gen_range(1..min_dim);
+ println!("desired rank = {}", rank);
+ random_with_rank(sh, rank)
+}
+
+/// create a random matrix and ensures it is full rank.
+pub fn full_rank<A, Sh>(sh: Sh) -> Array2<A>
+where
+ A: Scalar + Lapack,
+ Sh: ShapeBuilder<Dim = Ix2> + Clone,
+{
+ 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<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
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<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());
+}