Add eigenvalue decomposition example

Author: bytesnake

examples/truncated_eig.rs

@@ -0,0 +1,22 @@
+extern crate ndarray;
+extern crate ndarray_linalg;
+
+use ndarray::*;
+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);
+
+    println!("Generate spectrum: {:?}", &t);
+
+    let t = Array2::from_diag(&t);
+    let a = v.dot(&t.dot(&v.t()));
+
+    // calculate the truncated eigenproblem decomposition
+    for (val, _) in TruncatedEig::new(a, TruncatedOrder::Largest) {
+        println!("Found eigenvalue {}", val[0]);
+    }
+}

examples/truncated_svd.rs

@@ -0,0 +1,22 @@
+extern crate ndarray;
+extern crate ndarray_linalg;
+
+use ndarray::*;
+use ndarray_linalg::*;
+
+fn main() {
+    let a = arr2(&[[3., 2., 2.], [2., 3., -2.]]);
+
+    // calculate the truncated singular value decomposition for 2 singular values
+    let result = TruncatedSvd::new(a, TruncatedOrder::Largest).decompose(2).unwrap();
+
+    // acquire singular values, left-singular vectors and right-singular vectors
+    let (u, sigma, v_t) = result.values_vectors();
+    println!("Result of the singular value decomposition A = UΣV^T:");
+    println!(" === U ===");
+    println!("{:?}", u);
+    println!(" === Σ ===");
+    println!("{:?}", Array2::from_diag(&sigma));
+    println!(" === V^T ===");
+    println!("{:?}", v_t);
+}

src/lobpcg/eig.rs

@@ -83,6 +83,7 @@ impl<A: Float + Scalar + Lapack + PartialOrd + Default> IntoIterator for Truncat
     fn into_iter(self) -> TruncatedEigIterator<A> {
         TruncatedEigIterator {
             step_size: 1,
+            remaining: self.problem.len_of(Axis(0)),
             eig: self,
         }
     }
@@ -94,14 +95,20 @@ impl<A: Float + Scalar + Lapack + PartialOrd + Default> IntoIterator for Truncat
 /// eigenvalue/vector pair. Useful for generating pairs until a certain condition is met.
 pub struct TruncatedEigIterator<A: Scalar> {
     step_size: usize,
+    remaining: usize,
     eig: TruncatedEig<A>,
 }
 
 impl<A: Float + Scalar + Lapack + PartialOrd + Default> Iterator for TruncatedEigIterator<A> {
     type Item = (Array1<A>, Array2<A>);
 
     fn next(&mut self) -> Option<Self::Item> {
-        let res = self.eig.decompose(self.step_size);
+        if self.remaining == 0 {
+            return None;
+        }
+
+        let step_size = usize::min(self.step_size, self.remaining);
+        let res = self.eig.decompose(step_size);
 
         match res {
             EigResult::Ok(vals, vecs, norms) | EigResult::Err(vals, vecs, norms, _) => {
@@ -127,6 +134,7 @@ impl<A: Float + Scalar + Lapack + PartialOrd + Default> Iterator for TruncatedEi
                 };
 
                 self.eig.constraints = Some(new_constraints);
+                self.remaining -= step_size;
 
                 Some((vals, vecs))
             }

src/lobpcg/lobpcg.rs

@@ -87,27 +87,20 @@ fn apply_constraints<A: Scalar + Lapack>(
     y: ArrayView2<A>,
 ) {
     let gram_yv = y.t().dot(&v);
-    dbg!(&gram_yv.shape());
 
     let u = gram_yv
         .gencolumns()
         .into_iter()
         .map(|x| {
-            dbg!(&x.shape());
             let res = fact_yy.solvec(&x).unwrap();
 
-            dbg!(&res);
-
             res.to_vec()
         })
         .flatten()
         .collect::<Vec<A>>();
 
     let rows = gram_yv.len_of(Axis(0));
     let u = Array2::from_shape_vec((rows, u.len() / rows), u).unwrap();
-    dbg!(&u);
-    dbg!(y.shape());
-    dbg!(&v.shape());
 
     v -= &(y.dot(&u));
 }