Truncate svd with float eps; Take preconditioner as a closure

bytesnake · bytesnake · commit 68d8ec7c67d4 · 2020-03-19T10:08:27.000+01:00
diff --git a/examples/truncated_eig.rs b/examples/truncated_eig.rs
@@ -7,6 +7,7 @@ use ndarray_linalg::*;
 fn main() {
     let n = 10;
     let v = random_unitary(n);
+
     // set eigenvalues in decreasing order
     let t = Array1::linspace(n as f64, -(n as f64), n);
 
diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs
@@ -64,15 +64,27 @@ impl<A: Scalar + Lapack + PartialOrd + Default> TruncatedEig<A> {
         let x = Array2::random((self.problem.len_of(Axis(0)), num), Uniform::new(0.0, 1.0))
             .mapv(|x| NumCast::from(x).unwrap());
 
-        lobpcg(
-            |y| self.problem.dot(&y),
-            x,
-            self.preconditioner.clone(),
-            self.constraints.clone(),
-            self.precision,
-            self.maxiter,
-            self.order.clone(),
-        )
+        if let Some(ref preconditioner) = self.preconditioner {
+            lobpcg(
+                |y| self.problem.dot(&y),
+                x,
+                |mut y| y.assign(&preconditioner.dot(&y)),
+                self.constraints.clone(),
+                self.precision,
+                self.maxiter,
+                self.order.clone(),
+            )
+        } else {
+            lobpcg(
+                |y| self.problem.dot(&y),
+                x,
+                |_| {},
+                self.constraints.clone(),
+                self.precision,
+                self.maxiter,
+                self.order.clone(),
+            )
+        }
     }
 }
 
diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs
@@ -147,10 +147,10 @@ fn orthonormalize<T: Scalar + Lapack>(v: Array2<T>) -> Result<(Array2<T>, Array2
 /// for it. All iterations are tracked and the optimal solution returned. In case of an error a
 /// special variant `EigResult::NotConverged` additionally carries the error. This can happen when
 /// the precision of the matrix is too low (switch from `f32` to `f64` for example).
-pub fn lobpcg<A: Scalar + Lapack + PartialOrd + Default, F: Fn(ArrayView2<A>) -> Array2<A>>(
+pub fn lobpcg<A: Scalar + Lapack + PartialOrd + Default, F: Fn(ArrayView2<A>) -> Array2<A>, G: Fn(ArrayViewMut2<A>)>(
     a: F,
     mut x: Array2<A>,
-    m: Option<Array2<A>>,
+    m: G,
     y: Option<Array2<A>>,
     tol: A::Real,
     maxiter: usize,
@@ -246,9 +246,9 @@ pub fn lobpcg<A: Scalar + Lapack + PartialOrd + Default, F: Fn(ArrayView2<A>) ->
 
         // select active eigenvalues, apply pre-conditioner, orthogonalize to Y and orthonormalize
         let mut active_block_r = ndarray_mask(r.view(), &activemask);
-        if let Some(ref m) = m {
-            active_block_r = m.dot(&active_block_r);
-        }
+        // apply preconditioner
+        m(active_block_r.view_mut());
+
         if let (Some(ref y), Some(ref fact_yy)) = (&y, &fact_yy) {
             apply_constraints(active_block_r.view_mut(), fact_yy, y.view());
         }
@@ -453,7 +453,7 @@ mod tests {
         let n = a.len_of(Axis(0));
         let x: Array2<f64> = Array2::random((n, num), Uniform::new(0.0, 1.0));
 
-        let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, 2 * n, order);
+        let result = lobpcg(|y| a.dot(&y), x, |_| {}, None, 1e-10, 2 * n, order);
         match result {
             EigResult::Ok(vals, _, r_norms) | EigResult::Err(vals, _, r_norms, _) => {
                 // check convergence
@@ -510,7 +510,7 @@ mod tests {
         let x: Array2<f64> = Array2::random((10, 1), Uniform::new(0.0, 1.0));
         let y: Array2<f64> = arr2(&[[1.0, 0., 0., 0., 0., 0., 0., 0., 0., 0.]]).reversed_axes();
 
-        let result = lobpcg(|y| a.dot(&y), x, None, Some(y), 1e-10, 100, Order::Smallest);
+        let result = lobpcg(|y| a.dot(&y), x, |_| {}, Some(y), 1e-10, 100, Order::Smallest);
         dbg!(&result);
         match result {
             EigResult::Ok(vals, vecs, r_norms) | EigResult::Err(vals, vecs, r_norms, _) => {
diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs
@@ -1,18 +1,19 @@
+///! Truncated singular value decomposition
+///! 
+///! This module computes the k largest/smallest singular values/vectors for a dense matrix. 
 use super::lobpcg::{lobpcg, EigResult, Order};
 use crate::error::Result;
 use crate::{Lapack, Scalar};
 use ndarray::prelude::*;
 use ndarray_rand::rand_distr::Uniform;
 use ndarray_rand::RandomExt;
 use num_traits::{Float, NumCast};
-///! Implements truncated singular value decomposition
-///
 use std::ops::DivAssign;
 
-/// The result of an eigenvalue decomposition for SVD
+/// The result of a eigenvalue decomposition, not yet transformed into singular values/vectors
 ///
 /// Provides methods for either calculating just the singular values with reduced cost or the
-/// vectors as well.
+/// vectors with additional cost of matrix multiplication.
 #[derive(Debug)]
 pub struct TruncatedSvdResult<A> {
     eigvals: Array1<A>,
@@ -21,26 +22,31 @@ pub struct TruncatedSvdResult<A> {
     ngm: bool,
 }
 
-impl<A: Float + PartialOrd + DivAssign<A> + 'static> TruncatedSvdResult<A> {
+impl<A: Float + PartialOrd + DivAssign<A> + 'static + MagnitudeCorrection> TruncatedSvdResult<A> {
     /// Returns singular values ordered by magnitude with indices.
     fn singular_values_with_indices(&self) -> (Array1<A>, Vec<usize>) {
-        // numerate and square root eigenvalues
-        let mut a = self.eigvals.iter().map(|x| x.sqrt()).enumerate().collect::<Vec<_>>();
+        // numerate eigenvalues
+        let mut a = self.eigvals.iter().enumerate().collect::<Vec<_>>();
 
         // sort by magnitude
         a.sort_by(|(_, x), (_, y)| x.partial_cmp(&y).unwrap().reverse());
 
+        // calculate cut-off magnitude (borrowed from scipy)
+        let cutoff = A::epsilon() * // float precision
+                     A::correction() * // correction term (see trait below)
+                     *a[0].1; // max eigenvalue
+
         // filter low singular values away
         let (values, indices): (Vec<A>, Vec<usize>) = a
             .into_iter()
-            .filter(|(_, x)| *x > NumCast::from(1e-5).unwrap())
-            .map(|(a, b)| (b, a))
+            .filter(|(_, x)| *x > &cutoff)
+            .map(|(a, b)| (b.sqrt(), a))
             .unzip();
 
         (Array1::from(values), indices)
     }
 
-    /// Returns singular values orderd by magnitude
+    /// Returns singular values ordered by magnitude
     pub fn values(&self) -> Array1<A> {
         let (values, _) = self.singular_values_with_indices();
 
@@ -82,10 +88,8 @@ impl<A: Float + PartialOrd + DivAssign<A> + 'static> TruncatedSvdResult<A> {
 
 /// Truncated singular value decomposition
 ///
-/// This struct wraps the LOBPCG algorithm and provides convenient builder-pattern access to
-/// 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.
+/// Wraps the LOBPCG algorithm and provides convenient builder-pattern access to
+/// parameter like maximal iteration, precision and constraint matrix.
 pub struct TruncatedSvd<A: Scalar> {
     order: Order,
     problem: Array2<A>,
@@ -117,9 +121,15 @@ impl<A: Scalar + Lapack + PartialOrd + Default> TruncatedSvd<A> {
 
     // calculate the eigenvalue decomposition
     pub fn decompose(self, num: usize) -> Result<TruncatedSvdResult<A>> {
+        if num < 1 {
+            panic!("The number of singular values to compute should be larger than zero!");
+        }
+
         let (n, m) = (self.problem.nrows(), self.problem.ncols());
 
-        let x = Array2::random((usize::min(n, m), num), Uniform::new(0.0, 1.0)).mapv(|x| NumCast::from(x).unwrap());
+        // generate initial matrix
+        let x = Array2::random((usize::min(n, m), num), Uniform::new(0.0, 1.0))
+            .mapv(|x| NumCast::from(x).unwrap());
 
         // square precision because the SVD squares the eigenvalue as well
         let precision = self.precision * self.precision;
@@ -129,7 +139,7 @@ impl<A: Scalar + Lapack + PartialOrd + Default> TruncatedSvd<A> {
             lobpcg(
                 |y| self.problem.t().dot(&self.problem.dot(&y)),
                 x,
-                None,
+                |_| {},
                 None,
                 precision,
                 self.maxiter,
@@ -139,7 +149,7 @@ impl<A: Scalar + Lapack + PartialOrd + Default> TruncatedSvd<A> {
             lobpcg(
                 |y| self.problem.dot(&self.problem.t().dot(&y)),
                 x,
-                None,
+                |_| {},
                 None,
                 precision,
                 self.maxiter,
@@ -160,6 +170,22 @@ impl<A: Scalar + Lapack + PartialOrd + Default> TruncatedSvd<A> {
     }
 }
 
+pub trait MagnitudeCorrection {
+    fn correction() -> Self;
+}
+
+impl MagnitudeCorrection for f32 {
+    fn correction() -> Self {
+        1.0e3
+    }
+}
+
+impl MagnitudeCorrection for f64 {
+    fn correction() -> Self {
+        1.0e6
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::Order;
@@ -179,7 +205,7 @@ mod tests {
             .decompose(2)
             .unwrap();
 
-        let (_, sigma, _) = res.values_vecs();
+        let (_, sigma, _) = res.values_vectors();
 
         close_l2(&sigma, &arr1(&[5.0, 3.0]), 1e-5);
     }
