From c827c7f455802dbe3029921ab17d64ff67c60d1f Mon Sep 17 00:00:00 2001
From: Lishen <>
Date: Sat, 4 May 2024 22:35:09 +0300
Subject: [PATCH] add pseudo inverse for symmetric eigenvalue decomposition

---
 src/linalg/symmetric_eigen.rs | 66 +++++++++++++++++++++++++++++++++++
 tests/linalg/eigen.rs         | 17 +++++++++
 2 files changed, 83 insertions(+)

diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs
index 70cad8a1..90aaf1bd 100644
--- a/src/linalg/symmetric_eigen.rs
+++ b/src/linalg/symmetric_eigen.rs
@@ -8,6 +8,7 @@ use crate::allocator::Allocator;
 use crate::base::{DefaultAllocator, Matrix2, OMatrix, OVector, SquareMatrix, Vector2};
 use crate::dimension::{Dim, DimDiff, DimSub, U1};
 use crate::storage::Storage;
+use crate::RowOVector;
 use simba::scalar::ComplexField;
 
 use crate::linalg::givens::GivensRotation;
@@ -295,6 +296,71 @@ where
         u_t.adjoint_mut();
         &self.eigenvectors * u_t
     }
+
+    /// Computes the pseudo-inverse of this matrix.
+    ///
+    /// Calculate a generalized inverse of a complex Hermitian/real symmetric
+    /// matrix using its eigenvalue decomposition and including all eigenvalues
+    /// with 'large' absolute value.
+    ///
+    /// # Arguments
+    ///
+    /// * `atol` − absolute threshold term, if `None` provided value is 0.
+    /// * `rtol` − relative threshold term, if `None` provided value is `N * eps` where
+    /// `eps` is the machine precision value of the `T::RealField`.
+    #[must_use]
+    pub fn pseudo_inverse(
+        &self,
+        atol: Option<T::RealField>,
+        rtol: Option<T::RealField>,
+    ) -> OMatrix<T, D, D>
+    where
+        DefaultAllocator: Allocator<usize, D>,
+        DefaultAllocator: Allocator<T, U1, D>,
+    {
+        let u = &self.eigenvectors;
+        let s = &self.eigenvalues;
+        let max_s = s.camax();
+        let atol = atol.unwrap_or(T::RealField::zero());
+        let rtol =
+            rtol.unwrap_or(T::RealField::default_epsilon() * crate::convert(u.ncols() as f64));
+        assert!(
+            rtol >= T::RealField::zero() && atol >= T::RealField::zero(),
+            "atol and rtol values must be positive.",
+        );
+        let val = atol + max_s * rtol;
+        let mut above_cutoff = OVector::<usize, D>::zeros_generic(u.shape_generic().0, U1);
+        let mut r_take = 0;
+        for i in 0..s.len() {
+            if s[i].clone().abs() > val {
+                above_cutoff[r_take] = i;
+                r_take += 1;
+            }
+        }
+        let psigma_diag = RowOVector::<T, D>::from_fn_generic(U1, u.shape_generic().0, |_, j| {
+            if j < r_take {
+                T::from_real(s[above_cutoff[j]].clone().recip())
+            } else {
+                T::zero()
+            }
+        });
+        let u = OMatrix::<T, D, D>::from_fn_generic(
+            u.shape_generic().0,
+            u.shape_generic().1,
+            |i, j| {
+                if j < r_take {
+                    u[(i, above_cutoff[j])].clone()
+                } else {
+                    T::zero()
+                }
+            },
+        );
+        let mut up = u.clone();
+        for i in 0..u.nrows() {
+            up.row_mut(i).component_mul_assign(&psigma_diag);
+        }
+        up * u.conjugate().transpose()
+    }
 }
 
 /// Computes the wilkinson shift, i.e., the 2x2 symmetric matrix eigenvalue to its tailing
diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs
index 7a944c44..135c166b 100644
--- a/tests/linalg/eigen.rs
+++ b/tests/linalg/eigen.rs
@@ -62,6 +62,23 @@ mod proptest_tests {
 
                         prop_assert!(relative_eq!(m.lower_triangle(), recomp.lower_triangle(), epsilon = 1.0e-5))
                     }
+
+                    #[test]
+                    fn symmetric_eigen_pseudo_inverse(m in dmatrix_($scalar)) {
+                        let eig = m.clone().symmetric_eigen();
+                        let pinv = eig.pseudo_inverse(None, None);
+                        prop_assert!(relative_eq!(
+                            m,
+                            &m*&pinv*&m,
+                            epsilon = 1.0e-5
+                        ));
+                        prop_assert!(relative_eq!(
+                            pinv,
+                            &pinv*m*&pinv,
+                            epsilon = 1.0e-5
+                        ));
+                    }
+
                 }
             }
         }