Improve documentation and add MP law test

Author: bytesnake

ndarray-linalg/Cargo.toml

@@ -50,6 +50,7 @@ paste = "1.0"
 criterion = "0.3"
 # Keep the same version as ndarray's dependency!
 approx = { version = "0.4", features = ["num-complex"] }
+rand_xoshiro = "0.6"
 
 [[bench]]
 name = "truncated_eig"

ndarray-linalg/src/lobpcg/eig.rs

@@ -1,9 +1,10 @@
+//! Truncated eigenvalue decomposition
+//!
+
 use super::lobpcg::{lobpcg, LobpcgResult, Order};
 use crate::{generate, Scalar};
 use lax::Lapack;
 
-///! Implements truncated eigenvalue decomposition
-///
 use ndarray::prelude::*;
 use ndarray::stack;
 use ndarray::ScalarOperand;
@@ -15,6 +16,20 @@ use num_traits::{Float, NumCast};
 /// parameter like maximal iteration, precision and constraint matrix. Furthermore it allows
 /// conversion into a iterative solver where each iteration step yields a new eigenvalue/vector
 /// pair.
+///
+/// # Example
+///
+/// ```rust
+/// let diag = arr1(&[1., 2., 3., 4., 5.]);
+/// let a = Array2::from_diag(&diag);
+///
+/// let eig = TruncatedEig::new(a, Order::Largest)
+///    .precision(1e-5)
+///    .maxiter(500);
+///
+/// let res = eig.decompose();
+/// ```
+
 pub struct TruncatedEig<A: Scalar> {
     order: Order,
     problem: Array2<A>,
@@ -25,6 +40,11 @@ pub struct TruncatedEig<A: Scalar> {
 }
 
 impl<A: Float + Scalar + ScalarOperand + Lapack + PartialOrd + Default> TruncatedEig<A> {
+    /// Create a new truncated eigenproblem solver
+    ///
+    /// # Properties
+    /// * `problem`: problem matrix
+    /// * `order`: ordering of the eigenvalues with [TruncatedOrder](crate::TruncatedOrder)
     pub fn new(problem: Array2<A>, order: Order) -> TruncatedEig<A> {
         TruncatedEig {
             precision: 1e-5,
@@ -36,31 +56,68 @@ impl<A: Float + Scalar + ScalarOperand + Lapack + PartialOrd + Default> Truncate
         }
     }
 
+    /// Set desired precision
+    ///
+    /// This argument specifies the desired precision, which is passed to the LOBPCG solver. It
+    /// controls at which point the opimization of each eigenvalue is stopped. The precision is
+    /// global and applied to all eigenvalues with respect to their L2 norm.
+    ///
+    /// If the precision can't be reached and the maximum number of iteration is reached, then an
+    /// error is returned in [LobpcgResult](crate::lobpcg::LobpcgResult).
     pub fn precision(mut self, precision: f32) -> Self {
         self.precision = precision;
 
         self
     }
 
+    /// Set the maximal number of iterations
+    ///
+    /// The LOBPCG is an iterative approach to eigenproblems and stops when this maximum 
+    /// number of iterations are reached.
     pub fn maxiter(mut self, maxiter: usize) -> Self {
         self.maxiter = maxiter;
 
         self
     }
 
+    /// Construct a solution, which is orthogonal to this
+    ///
+    /// If a number of eigenvectors are already known, then this function can be used to construct
+    /// a orthogonal subspace. Also used with an iterative approach.
     pub fn orthogonal_to(mut self, constraints: Array2<A>) -> Self {
         self.constraints = Some(constraints);
 
         self
     }
 
+    /// Apply a preconditioner
+    ///
+    /// A preconditioning matrix can speed up the solving process by improving the spectral
+    /// distribution of the eigenvalues. It requires prior knowledge of the problem.
     pub fn precondition_with(mut self, preconditioner: Array2<A>) -> Self {
         self.preconditioner = Some(preconditioner);
 
         self
     }
 
-    // calculate the eigenvalues decompose
+    /// Calculate the eigenvalue decomposition
+    ///
+    /// # Parameters
+    ///
+    ///  * `num`: number of eigenvalues ordered by magnitude
+    ///
+    /// # Example
+    ///
+    /// ```rust
+    /// let diag = arr1(&[1., 2., 3., 4., 5.]);
+    /// let a = Array2::from_diag(&diag);
+    ///
+    /// let eig = TruncatedEig::new(a, Order::Largest)
+    ///    .precision(1e-5)
+    ///    .maxiter(500);
+    ///
+    /// let res = eig.decompose();
+    /// ```
     pub fn decompose(&self, num: usize) -> LobpcgResult<A> {
         let x: Array2<f64> = generate::random((self.problem.len_of(Axis(0)), num));
         let x = x.mapv(|x| NumCast::from(x).unwrap());
@@ -104,10 +161,24 @@ impl<A: Float + Scalar + ScalarOperand + Lapack + PartialOrd + Default> IntoIter
     }
 }
 
-/// Truncate eigenproblem iterator
+/// Truncated eigenproblem iterator
 ///
 /// This wraps a truncated eigenproblem and provides an iterator where each step yields a new
 /// eigenvalue/vector pair. Useful for generating pairs until a certain condition is met.
+///
+/// # Example
+///
+/// ```rust
+/// let teig = TruncatedEig::new(a, Order::Largest)
+///     .precision(1e-5)
+///     .maxiter(500);
+/// 
+/// // solve eigenproblem until eigenvalues get smaller than 0.5
+/// let res = teig.into_iter()
+///     .take_while(|x| x.0[0] > 0.5)
+///     .flat_map(|x| x.0)
+///     .collect();
+/// ```
 pub struct TruncatedEigIterator<A: Scalar> {
     step_size: usize,
     remaining: usize,

ndarray-linalg/src/lobpcg/mod.rs

@@ -2,6 +2,6 @@ mod eig;
 mod lobpcg;
 mod svd;
 
-pub use eig::TruncatedEig;
+pub use eig::{TruncatedEig, TruncatedEigIterator};
 pub use lobpcg::{lobpcg, LobpcgResult, Order as TruncatedOrder};
-pub use svd::TruncatedSvd;
+pub use svd::{TruncatedSvd, MagnitudeCorrection};

ndarray-linalg/src/lobpcg/svd.rs

@@ -173,6 +173,11 @@ impl<A: Float + Scalar + ScalarOperand + Lapack + PartialOrd + Default> Truncate
     }
 }
 
+/// Magnitude Correction
+///
+/// The magnitude correction changes the cut-off point at which an eigenvector belongs to the
+/// null-space and its eigenvalue is therefore zero. The correction is multiplied by the floating
+/// point epsilon and therefore dependent on the floating type.
 pub trait MagnitudeCorrection {
     fn correction() -> Self;
 }
@@ -196,6 +201,8 @@ mod tests {
     use crate::{close_l2, generate};
 
     use ndarray::{arr1, arr2, Array2};
+    use rand_xoshiro::Xoshiro256Plus;
+    use approx::assert_abs_diff_eq;
 
     #[test]
     fn test_truncated_svd() {
@@ -227,4 +234,56 @@ mod tests {
 
         close_l2(&a, &reconstructed, 1e-5);
     }
+
+    /// Eigenvalue structure in high dimensions
+    /// 
+    /// This test checks that the eigenvalues are following the Marchensko-Pastur law. The data is
+    /// standard uniformly distributed (i.e. E(x) = 0, E^2(x) = 1) and we have twice the amount of
+    /// data when compared to features. The probability density of the eigenvalues should then follow
+    /// a special densitiy function, described by the Marchenko-Pastur law.
+    /// 
+    /// See also https://en.wikipedia.org/wiki/Marchenko%E2%80%93Pastur_distribution
+    #[test]
+    fn test_marchenko_pastur() {
+        // create random number generator
+        let mut rng = SmallRng::seed_from_u64(3);
+
+        // generate normal distribution random data with N >> p
+        let data = Array2::random_using((1000, 500), StandardNormal, &mut rng);
+        let dataset = Dataset::from(data / 1000f64.sqrt());
+
+        let model = Pca::params(500).fit(&dataset);
+        let sv = model.singular_values().mapv(|x| x * x); 
+
+        // we have created a random spectrum and can apply the Marchenko-Pastur law
+        // with variance 1 and p/n = 0.5
+        let (a, b) = ( 
+            1. * (1. - 0.5f64.sqrt()).powf(2.0),
+            1. * (1. + 0.5f64.sqrt()).powf(2.0),
+        );  
+
+        // check that the spectrum has correct boundaries
+        assert_abs_diff_eq!(b, sv[0], epsilon = 0.1);
+        assert_abs_diff_eq!(a, sv[sv.len() - 1], epsilon = 0.1);
+
+        // estimate density empirical and compare with Marchenko-Pastur law
+        let mut i = 0;
+        'outer: for th in Array1::linspace(0.1, 2.8, 28).into_iter().rev() {
+            let mut count = 0;
+            while sv[i] >= *th {
+                count += 1;
+                i += 1;
+
+                if i == sv.len() {
+                    break 'outer;
+                }   
+            }   
+
+            let x = th + 0.05;
+            let mp_law = ((b - x) * (x - a)).sqrt() / std::f64::consts::PI / x;
+            let empirical = count as f64 / 500. / ((2.8 - 0.1) / 28.);
+
+            assert_abs_diff_eq!(mp_law, empirical, epsilon = 0.05);
+        }   
+    } 
 }