From 9f2eb2b374ffd728b91664213e9e2338946f2be2 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Thu, 12 Mar 2020 20:53:24 +0100 Subject: [PATCH 01/30] Add eigenvalue decomposition for generalized problem --- src/eigh.rs | 54 ++++++++++++++++++++++++++++++++++++++++++++++ src/lapack/eigh.rs | 19 +++++++++++----- 2 files changed, 68 insertions(+), 5 deletions(-) diff --git a/src/eigh.rs b/src/eigh.rs index 5e2162b3..1a194ca7 100644 --- a/src/eigh.rs +++ b/src/eigh.rs @@ -42,6 +42,19 @@ where } } +impl EighInto for (ArrayBase, ArrayBase) +where + A: Scalar + Lapack, + S: DataMut, +{ + type EigVal = Array1; + + fn eigh_into(mut self, uplo: UPLO) -> Result<(Self::EigVal, Self)> { + let (val, _) = self.eigh_inplace(uplo)?; + Ok((val, self)) + } +} + impl Eigh for ArrayBase where A: Scalar + Lapack, @@ -56,6 +69,20 @@ where } } +impl Eigh for (ArrayBase, ArrayBase) +where + A: Scalar + Lapack, + S: Data, +{ + type EigVal = Array1; + type EigVec = Array2; + + fn eigh(&self, uplo: UPLO) -> Result<(Self::EigVal, Self::EigVec)> { + let (a,b) = (self.0.to_owned(), self.1.to_owned()); + (a,b).eigh_into(uplo).map(|x| (x.0, (x.1).0)) + } +} + impl EighInplace for ArrayBase where A: Scalar + Lapack, @@ -75,6 +102,33 @@ where } } +impl EighInplace for (ArrayBase, ArrayBase) +where + A: Scalar + Lapack, + S: DataMut, +{ + type EigVal = Array1; + + fn eigh_inplace(&mut self, uplo: UPLO) -> Result<(Self::EigVal, &mut Self)> { + let layout = self.0.square_layout()?; + // XXX Force layout to be Fortran (see #146) + match layout { + MatrixLayout::C(_) => self.0.swap_axes(0, 1), + MatrixLayout::F(_) => {} + } + + let layout = self.1.square_layout()?; + match layout { + MatrixLayout::C(_) => self.1.swap_axes(0, 1), + MatrixLayout::F(_) => {} + } + + let s = unsafe { A::eigh_generalized(true, self.0.square_layout()?, uplo, self.0.as_allocated_mut()?, self.1.as_allocated_mut()?)? }; + + Ok((ArrayBase::from(s), self)) + } +} + /// Calculate eigenvalues without eigenvectors pub trait EigValsh { type EigVal; diff --git a/src/lapack/eigh.rs b/src/lapack/eigh.rs index 610564d2..a93e702e 100644 --- a/src/lapack/eigh.rs +++ b/src/lapack/eigh.rs @@ -12,10 +12,11 @@ use super::{into_result, UPLO}; /// Wraps `*syev` for real and `*heev` for complex pub trait Eigh_: Scalar { unsafe fn eigh(calc_eigenvec: bool, l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result>; + unsafe fn eigh_generalized(calc_eigenvec: bool, l: MatrixLayout, uplo: UPLO, a: &mut [Self], b: &mut[Self]) -> Result>; } macro_rules! impl_eigh { - ($scalar:ty, $ev:path) => { + ($scalar:ty, $ev:path, $evg:path) => { impl Eigh_ for $scalar { unsafe fn eigh(calc_v: bool, l: MatrixLayout, uplo: UPLO, mut a: &mut [Self]) -> Result> { let (n, _) = l.size(); @@ -24,11 +25,19 @@ macro_rules! impl_eigh { let info = $ev(l.lapacke_layout(), jobz, uplo as u8, n, &mut a, n, &mut w); into_result(info, w) } + + unsafe fn eigh_generalized(calc_v: bool, l: MatrixLayout, uplo: UPLO, mut a: &mut [Self], mut b: &mut [Self]) -> Result> { + let (n, _) = l.size(); + let jobz = if calc_v { b'V' } else { b'N' }; + let mut w = vec![Self::Real::zero(); n as usize]; + let info = $evg(l.lapacke_layout(), 1, jobz, uplo as u8, n, &mut a, n, &mut b, n, &mut w); + into_result(info, w) + } } }; } // impl_eigh! -impl_eigh!(f64, lapacke::dsyev); -impl_eigh!(f32, lapacke::ssyev); -impl_eigh!(c64, lapacke::zheev); -impl_eigh!(c32, lapacke::cheev); +impl_eigh!(f64, lapacke::dsyev, lapacke::dsygv); +impl_eigh!(f32, lapacke::ssyev, lapacke::ssygv); +impl_eigh!(c64, lapacke::zheev, lapacke::zhegv); +impl_eigh!(c32, lapacke::cheev, lapacke::chegv); From 8b783146ccddb33dfbdeb4f8e65d0e189f371b36 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Thu, 12 Mar 2020 20:57:53 +0100 Subject: [PATCH 02/30] Unify return types for EighInplace, EighInto, Eigh --- src/eigh.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/eigh.rs b/src/eigh.rs index 1a194ca7..c52c2abe 100644 --- a/src/eigh.rs +++ b/src/eigh.rs @@ -75,11 +75,11 @@ where S: Data, { type EigVal = Array1; - type EigVec = Array2; + type EigVec = (Array2, Array2); fn eigh(&self, uplo: UPLO) -> Result<(Self::EigVal, Self::EigVec)> { let (a,b) = (self.0.to_owned(), self.1.to_owned()); - (a,b).eigh_into(uplo).map(|x| (x.0, (x.1).0)) + (a,b).eigh_into(uplo) } } From 64da1007fd79bfe96d0628486372309c4296872e Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Thu, 12 Mar 2020 21:05:21 +0100 Subject: [PATCH 03/30] Relax requirement for mass matrix B --- src/eigh.rs | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/eigh.rs b/src/eigh.rs index c52c2abe..9066281e 100644 --- a/src/eigh.rs +++ b/src/eigh.rs @@ -42,10 +42,11 @@ where } } -impl EighInto for (ArrayBase, ArrayBase) +impl EighInto for (ArrayBase, ArrayBase) where A: Scalar + Lapack, S: DataMut, + S2: DataMut, { type EigVal = Array1; @@ -69,10 +70,11 @@ where } } -impl Eigh for (ArrayBase, ArrayBase) +impl Eigh for (ArrayBase, ArrayBase) where A: Scalar + Lapack, S: Data, + S2: Data, { type EigVal = Array1; type EigVec = (Array2, Array2); @@ -102,10 +104,11 @@ where } } -impl EighInplace for (ArrayBase, ArrayBase) +impl EighInplace for (ArrayBase, ArrayBase) where A: Scalar + Lapack, S: DataMut, + S2: DataMut, { type EigVal = Array1; From be9dcc637437b23fe489fd1291ff17370ded0be6 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Sat, 14 Mar 2020 00:12:31 +0100 Subject: [PATCH 04/30] Add LOBPCG --- Cargo.toml | 1 + examples/solveh.rs | 2 +- src/lib.rs | 3 + src/lobpcg.rs | 425 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 430 insertions(+), 1 deletion(-) create mode 100644 src/lobpcg.rs diff --git a/Cargo.toml b/Cargo.toml index 08120c51..d8c495ad 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -51,6 +51,7 @@ optional = true [dev-dependencies] paste = "0.1" criterion = "0.3" +ndarray-rand = "0.11" [[bench]] name = "eigh" diff --git a/examples/solveh.rs b/examples/solveh.rs index 8ddc73a0..30b2897b 100644 --- a/examples/solveh.rs +++ b/examples/solveh.rs @@ -10,7 +10,7 @@ fn solve() -> Result<(), error::LinalgError> { let b: Array1 = random(3); println!("b = {:?}", &b); let x = a.solveh(&b)?; - println!("Ax = {:?}", a.dot(&x));; + println!("Ax = {:?}", a.dot(&x)); Ok(()) } diff --git a/src/lib.rs b/src/lib.rs index 198449cc..0a189c06 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -36,6 +36,8 @@ //! - [Random matrix generators](generate/index.html) //! - [Scalar trait](types/trait.Scalar.html) +#[macro_use] extern crate ndarray; + extern crate blas_src; extern crate lapack_src; @@ -61,6 +63,7 @@ pub mod svddc; pub mod trace; pub mod triangular; pub mod types; +pub mod lobpcg; pub use assert::*; pub use cholesky::*; diff --git a/src/lobpcg.rs b/src/lobpcg.rs new file mode 100644 index 00000000..f1364132 --- /dev/null +++ b/src/lobpcg.rs @@ -0,0 +1,425 @@ +use num_traits::NumCast; +use ndarray::prelude::*; +use ndarray::OwnedRepr; +use crate::{cholesky::*, triangular::*, eigh::*, norm::*, close_l2}; +//use sprs::CsMat; +use crate::{Scalar, Lapack}; +use crate::error::{Result, LinalgError}; + +pub enum Order { + Largest, + Smallest +} + +#[derive(Debug)] +pub enum EigResult { + Ok(Array1, Array2, Vec), + Err(Array1, Array2, Vec, LinalgError), + NoResult(LinalgError) +} + +fn sorted_eig(a: ArrayView2, b: Option>, size: usize, order: &Order) -> Result<(Array1, Array2)> { + //close_l2(&input, &input.t(), 1e-4); + + let (vals, vecs) = match b { + Some(b) => (a, b).eigh(UPLO::Upper).map(|x| (x.0, (x.1).0))?, + _ => a.eigh(UPLO::Upper)? + }; + + let n = a.len_of(Axis(0)); + + Ok(match order { + Order::Largest => (vals.slice_move(s![n-size..; -1]).mapv(|x| Scalar::from_real(x)), vecs.slice_move(s![.., n-size..; -1])), + Order::Smallest => (vals.slice_move(s![..size]).mapv(|x| Scalar::from_real(x)), vecs.slice_move(s![.., ..size])) + }) +} + +fn ndarray_mask(matrix: ArrayView2, mask: &[bool]) -> Array2 { + let (rows, cols) = (matrix.nrows(), matrix.ncols()); + + assert_eq!(mask.len(), cols); + + let n_positive = mask.iter().filter(|x| **x).count(); + + let matrix = matrix.gencolumns().into_iter().zip(mask.iter()) + .filter(|(_,x)| **x) + .map(|(x,_)| x.to_vec()) + .flatten() + .collect::>(); + + Array2::from_shape_vec((n_positive, rows), matrix).unwrap().reversed_axes() +} + +fn apply_constraints( + mut v: ArrayViewMut, + fact_yy: &CholeskyFactorized>, + y: ArrayView2 +) { + let gram_yv = y.t().dot(&v); + + let u = gram_yv.genrows().into_iter() + .map(|x| fact_yy.solvec(&x).unwrap().to_vec()) + .flatten() + .collect::>(); + + let u = Array2::from_shape_vec((5, 5), u).unwrap(); + + v -= &(y.dot(&u)); +} + +fn orthonormalize( + v: Array2 +) -> Result<(Array2, Array2)> { + let gram_vv = v.t().dot(&v); + let gram_vv_fac = gram_vv.cholesky(UPLO::Lower)?; + + close_l2(&gram_vv, &gram_vv_fac.dot(&gram_vv_fac.t()), NumCast::from(1e-5).unwrap()); + + let v_t = v.reversed_axes(); + let u = gram_vv_fac.solve_triangular(UPLO::Lower, Diag::NonUnit, &v_t)? + .reversed_axes(); + + Ok((u, gram_vv_fac)) +} + +pub fn lobpcg) -> Array2>( + a: F, + x: Array2, + m: Option>, + y: Option>, + tol: A::Real, maxiter: usize, + order: Order +) -> EigResult { + // the target matrix should be symmetric and quadratic + //assert!(sprs::is_symmetric(&A)); + + // the initital approximation should be maximal square + // n is the dimensionality of the problem + let (n, size_x) = (x.nrows(), x.ncols()); + assert!(size_x <= n); + + let size_y = match y { + Some(ref y) => y.ncols(), + _ => 0 + }; + + if (n - size_y) < 5 * size_x { + panic!("Please use a different approach, the LOBPCG method only supports the calculation of a couple of eigenvectors!"); + } + + // cap the number of iteration + let mut iter = usize::min(n, maxiter); + + // factorize yy for later use + let fact_yy = y.as_ref().map(|x| x.t().dot(x).factorizec(UPLO::Upper).unwrap()); + + // orthonormalize the initial guess and calculate matrices AX and XAX + let (x, _) = match orthonormalize(x) { + Ok(x) => x, + Err(err) => return EigResult::NoResult(err) + }; + + let ax = a(x.view()); + let xax = x.t().dot(&ax); + + // perform eigenvalue decomposition on XAX + let (mut lambda, eig_block) = match sorted_eig(xax.view(), None, size_x, &order) { + Ok(x) => x, + Err(err) => return EigResult::NoResult(err) + }; + + //dbg!(&lambda, &eig_block); + + // initiate X and AX with eigenvector + let mut x = x.dot(&eig_block); + let mut ax = ax.dot(&eig_block); + + //dbg!(&X, &AX); + let mut activemask = vec![true; size_x]; + let mut residual_norms = Vec::new(); + let mut results = vec![(lambda.clone(), x.clone())]; + + let mut previous_block_size = size_x; + + let mut ident: Array2 = Array2::eye(size_x); + let ident0: Array2 = Array2::eye(size_x); + + let mut ap: Option<(Array2, Array2)> = None; + + let final_norm = loop { + // calculate residual + let lambda_tmp = lambda.clone().insert_axis(Axis(0)); + let tmp = &x * &lambda_tmp; + + let r = &ax - &tmp; + + // calculate L2 norm of error for every eigenvalue + let tmp = r.gencolumns().into_iter().map(|x| x.norm()).collect::>(); + residual_norms.push(tmp.clone()); + + // disable eigenvalues which are below the tolerance threshold + activemask = tmp.iter().zip(activemask.iter()).map(|(x, a)| *x > tol && *a).collect(); + + // resize identity block if necessary + let current_block_size = activemask.iter().filter(|x| **x).count(); + if current_block_size != previous_block_size { + previous_block_size = current_block_size; + ident = Array2::eye(current_block_size); + } + + // if we are below the threshold for all eigenvalue or exceeded the number of iteration, + // abort + if current_block_size == 0 || iter == 0 { + break Ok(tmp); + } + + // 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); + } + if let (Some(ref y), Some(ref fact_yy)) = (&y, &fact_yy) { + apply_constraints(active_block_r.view_mut(), fact_yy, y.view()); + } + + let (r,_) = match orthonormalize(active_block_r) { + Ok(x) => x, + Err(err) => break Err(err) + }; + + let ar = a(r.view()); + + // perform the Rayleigh Ritz procedure + let xaw = x.t().dot(&ar); + let waw = r.t().dot(&ar); + let xw = x.t().dot(&r); + + // compute symmetric gram matrices + let (gram_a, gram_b, active_p, active_ap) = if let Some((ref p, ref ap)) = ap { + let active_p = ndarray_mask(p.view(), &activemask); + let active_ap = ndarray_mask(ap.view(), &activemask); + + let (active_p, p_r) = orthonormalize(active_p).unwrap(); + //dbg!(&active_P, &P_R); + let active_ap = match p_r.solve_triangular(UPLO::Lower, Diag::NonUnit, &active_ap.reversed_axes()) { + Ok(x) => x, + Err(err) => break Err(err) + }; + + let active_ap = active_ap.reversed_axes(); + + //dbg!(&active_AP); + //dbg!(&R); + + let xap = x.t().dot(&active_ap); + let wap = r.t().dot(&active_ap); + let pap = active_p.t().dot(&active_ap); + let xp = x.t().dot(&active_p); + let wp = r.t().dot(&active_p); + + ( + stack![Axis(0), + stack![Axis(1), Array2::from_diag(&lambda), xaw, xap], + stack![Axis(1), xaw.t(), waw, wap], + stack![Axis(1), xap.t(), wap.t(), pap] + ], + + stack![Axis(0), + stack![Axis(1), ident0, xw, xp], + stack![Axis(1), xw.t(), ident, wp], + stack![Axis(1), xp.t(), wp.t(), ident] + ], + Some(active_p), + Some(active_ap) + ) + } else { + ( + stack![Axis(0), + stack![Axis(1), Array2::from_diag(&lambda), xaw], + stack![Axis(1), xaw.t(), waw] + ], + stack![Axis(0), + stack![Axis(1), ident0, xw], + stack![Axis(1), xw.t(), ident] + ], + None, + None + ) + }; + + //assert!(is_symmetric(gramA.view())); + //assert!(is_symmetric(gramB.view())); + + //dbg!(&gramA, &gramB); + let (new_lambda, eig_vecs) = match sorted_eig(gram_a.view(), Some(gram_b.view()), size_x, &order) { + Ok(x) => x, + Err(err) => break Err(err) + }; + lambda = new_lambda; + + //dbg!(&lambda, &eig_vecs); + let (pp, app, eig_x) = if let (Some(_), (Some(ref active_p), Some(ref active_ap))) = (ap, (active_p, active_ap)) { + + let eig_x = eig_vecs.slice(s![..size_x, ..]); + let eig_r = eig_vecs.slice(s![size_x..size_x+current_block_size, ..]); + let eig_p = eig_vecs.slice(s![size_x+current_block_size.., ..]); + + //dbg!(&eig_X); + //dbg!(&eig_R); + //dbg!(&eig_P); + + //dbg!(&R, &AR, &active_P, &active_AP); + + let pp = r.dot(&eig_r) + active_p.dot(&eig_p); + let app = ar.dot(&eig_r) + active_ap.dot(&eig_p); + + //dbg!(&pp); + //dbg!(&app); + + (pp, app, eig_x) + } else { + let eig_x = eig_vecs.slice(s![..size_x, ..]); + let eig_r = eig_vecs.slice(s![size_x.., ..]); + + let pp = r.dot(&eig_r); + let app = ar.dot(&eig_r); + + (pp, app, eig_x) + }; + + x = x.dot(&eig_x) + &pp; + ax = ax.dot(&eig_x) + &app; + + results.push((lambda.clone(), x.clone())); + + //dbg!(&X); + //dbg!(&AX); + + ap = Some((pp, app)); + + //dbg!(&ap); + + iter -= 1; + }; + + let best_idx = residual_norms.iter() + .enumerate() + .min_by(|&(_, item1): &(usize, &Vec), &(_, item2): &(usize, &Vec)| { + let norm1: A::Real = item1.iter().map(|x| (*x)*(*x)).sum(); + let norm2: A::Real = item2.iter().map(|x| (*x)*(*x)).sum(); + norm1.partial_cmp(&norm2).unwrap() + }); + + match best_idx { + Some((idx, norms)) => { + let (vals, vecs) = results[idx].clone(); + let norms = norms.iter().map(|x| Scalar::from_real(*x)).collect(); + + match final_norm { + Ok(_) => EigResult::Ok(vals, vecs, norms), + Err(err) => EigResult::Err(vals, vecs, norms, err) + } + }, + None => { + match final_norm { + Ok(_) => panic!("Not error available!"), + Err(err) => EigResult::NoResult(err) + } + } + } +} + +#[cfg(test)] +mod tests { + use super::sorted_eig; + use super::orthonormalize; + use super::ndarray_mask; + use super::Order; + use super::lobpcg; + use super::EigResult; + use crate::close_l2; + use ndarray::prelude::*; + use crate::qr::*; + use ndarray_rand::RandomExt; + use ndarray_rand::rand_distr::Uniform; + + #[test] + fn test_sorted_eigen() { + let matrix = Array2::random((10, 10), Uniform::new(0., 10.)); + let matrix = matrix.t().dot(&matrix); + + // return all eigenvectors with largest first + let (vals, vecs) = sorted_eig(matrix.view(), None, 10, &Order::Largest).unwrap(); + + // calculate V * A * V' and compare to original matrix + let diag = Array2::from_diag(&vals); + let rec = (vecs.dot(&diag)).dot(&vecs.t()); + + close_l2(&matrix, &rec, 1e-5); + } + + #[test] + fn test_masking() { + let matrix = Array2::random((10, 5), Uniform::new(0., 10.)); + let masked_matrix = ndarray_mask(matrix.view(), &[true, true, false, true, false]); + close_l2(&masked_matrix.slice(s![.., 2]), &matrix.slice(s![.., 3]), 1e-12); + } + + #[test] + fn test_orthonormalize() { + let matrix = Array2::random((10, 10), Uniform::new(-10., 10.)); + + let (n, l) = orthonormalize(matrix.clone()).unwrap(); + + // check for orthogonality + let identity = n.dot(&n.t()); + close_l2(&identity, &Array2::eye(10), 1e-2); + + // compare returned factorization with QR decomposition + let (_, r) = matrix.qr().unwrap(); + close_l2(&r.mapv(|x: f32| x.abs()) , &l.t().mapv(|x| x.abs()), 1e-2); + } + + #[test] + fn test_eigsolver_diag() { + let diag = arr1(&[1.,2.,3.,4.,5.,6.,7.,8.,9.,10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20.]); + let a = Array2::from_diag(&diag); + let x: Array2 = Array2::random((20, 3), Uniform::new(0.0, 1.0)); + + let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, 20, Order::Smallest); + match result { + EigResult::Ok(vals, _, _) => close_l2(&vals, &arr1(&[1.0, 2.0, 3.0]), 1e-5), + EigResult::Err(vals, _,_,_) => close_l2(&vals, &arr1(&[1.0, 2.0, 3.0]), 1e-5), + EigResult::NoResult(err) => panic!("Did not converge: {:?}", err) + } + + let x: Array2 = Array2::random((20, 3), Uniform::new(0.0, 1.0)); + let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, 20, Order::Largest); + match result { + EigResult::Ok(vals, _, _) => close_l2(&vals, &arr1(&[20.0, 19.0, 18.0]), 1e-5), + EigResult::Err(vals, _,_,_) => close_l2(&vals, &arr1(&[20.0, 19.0, 18.0]), 1e-5), + EigResult::NoResult(err) => panic!("Did not converge: {:?}", err) + } + } + + #[test] + fn test_eigsolver() { + let n = 50; + let tmp = Array2::random((n, n), Uniform::new(0.0, 1.0)); + let (v, _) = orthonormalize(tmp).unwrap(); + + // set eigenvalues in decreasing order + let t = Array2::from_diag(&Array1::linspace(n as f64, 1.0, n)); + let a = v.dot(&t.dot(&v.t())); + + let x: Array2 = Array2::random((n, 5), Uniform::new(0.0, 1.0)); + + let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, 20, Order::Largest); + match result { + EigResult::Ok(vals, _, _) | EigResult::Err(vals, _, _, _) => { + close_l2(&vals, &arr1(&[50.0, 49.0, 48.0, 47.0, 46.0]), 1e-5); + }, + EigResult::NoResult(err) => panic!("Did not converge: {:?}", err) + } + } +} From 80e6a617e92f53103c8a016147ba50cbe8e70d1f Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Sat, 14 Mar 2020 00:36:55 +0100 Subject: [PATCH 05/30] Add some comments on the LOBPCG and run `cargo fmt` --- src/eigh.rs | 14 ++- src/lapack/eigh.rs | 29 +++++- src/lapack/svddc.rs | 16 +-- src/lib.rs | 5 +- src/lobpcg.rs | 233 +++++++++++++++++++++++++++----------------- tests/svddc.rs | 4 +- 6 files changed, 191 insertions(+), 110 deletions(-) diff --git a/src/eigh.rs b/src/eigh.rs index 9066281e..1fc7a031 100644 --- a/src/eigh.rs +++ b/src/eigh.rs @@ -80,8 +80,8 @@ where type EigVec = (Array2, Array2); fn eigh(&self, uplo: UPLO) -> Result<(Self::EigVal, Self::EigVec)> { - let (a,b) = (self.0.to_owned(), self.1.to_owned()); - (a,b).eigh_into(uplo) + let (a, b) = (self.0.to_owned(), self.1.to_owned()); + (a, b).eigh_into(uplo) } } @@ -126,7 +126,15 @@ where MatrixLayout::F(_) => {} } - let s = unsafe { A::eigh_generalized(true, self.0.square_layout()?, uplo, self.0.as_allocated_mut()?, self.1.as_allocated_mut()?)? }; + let s = unsafe { + A::eigh_generalized( + true, + self.0.square_layout()?, + uplo, + self.0.as_allocated_mut()?, + self.1.as_allocated_mut()?, + )? + }; Ok((ArrayBase::from(s), self)) } diff --git a/src/lapack/eigh.rs b/src/lapack/eigh.rs index a93e702e..11b40bb0 100644 --- a/src/lapack/eigh.rs +++ b/src/lapack/eigh.rs @@ -12,7 +12,13 @@ use super::{into_result, UPLO}; /// Wraps `*syev` for real and `*heev` for complex pub trait Eigh_: Scalar { unsafe fn eigh(calc_eigenvec: bool, l: MatrixLayout, uplo: UPLO, a: &mut [Self]) -> Result>; - unsafe fn eigh_generalized(calc_eigenvec: bool, l: MatrixLayout, uplo: UPLO, a: &mut [Self], b: &mut[Self]) -> Result>; + unsafe fn eigh_generalized( + calc_eigenvec: bool, + l: MatrixLayout, + uplo: UPLO, + a: &mut [Self], + b: &mut [Self], + ) -> Result>; } macro_rules! impl_eigh { @@ -26,11 +32,28 @@ macro_rules! impl_eigh { into_result(info, w) } - unsafe fn eigh_generalized(calc_v: bool, l: MatrixLayout, uplo: UPLO, mut a: &mut [Self], mut b: &mut [Self]) -> Result> { + unsafe fn eigh_generalized( + calc_v: bool, + l: MatrixLayout, + uplo: UPLO, + mut a: &mut [Self], + mut b: &mut [Self], + ) -> Result> { let (n, _) = l.size(); let jobz = if calc_v { b'V' } else { b'N' }; let mut w = vec![Self::Real::zero(); n as usize]; - let info = $evg(l.lapacke_layout(), 1, jobz, uplo as u8, n, &mut a, n, &mut b, n, &mut w); + let info = $evg( + l.lapacke_layout(), + 1, + jobz, + uplo as u8, + n, + &mut a, + n, + &mut b, + n, + &mut w, + ); into_result(info, w) } } diff --git a/src/lapack/svddc.rs b/src/lapack/svddc.rs index 9c59b7aa..22f770e3 100644 --- a/src/lapack/svddc.rs +++ b/src/lapack/svddc.rs @@ -3,10 +3,10 @@ use num_traits::Zero; use crate::error::*; use crate::layout::MatrixLayout; -use crate::types::*; use crate::svddc::UVTFlag; +use crate::types::*; -use super::{SVDOutput, into_result}; +use super::{into_result, SVDOutput}; pub trait SVDDC_: Scalar { unsafe fn svddc(l: MatrixLayout, jobz: UVTFlag, a: &mut [Self]) -> Result>; @@ -15,11 +15,7 @@ pub trait SVDDC_: Scalar { macro_rules! impl_svdd { ($scalar:ty, $gesdd:path) => { impl SVDDC_ for $scalar { - unsafe fn svddc( - l: MatrixLayout, - jobz: UVTFlag, - mut a: &mut [Self], - ) -> Result> { + unsafe fn svddc(l: MatrixLayout, jobz: UVTFlag, mut a: &mut [Self]) -> Result> { let (m, n) = l.size(); let k = m.min(n); let lda = l.lda(); @@ -51,11 +47,7 @@ macro_rules! impl_svdd { SVDOutput { s: s, u: if jobz == UVTFlag::None { None } else { Some(u) }, - vt: if jobz == UVTFlag::None { - None - } else { - Some(vt) - }, + vt: if jobz == UVTFlag::None { None } else { Some(vt) }, }, ) } diff --git a/src/lib.rs b/src/lib.rs index 0a189c06..8eac92ee 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -36,7 +36,8 @@ //! - [Random matrix generators](generate/index.html) //! - [Scalar trait](types/trait.Scalar.html) -#[macro_use] extern crate ndarray; +#[macro_use] +extern crate ndarray; extern crate blas_src; extern crate lapack_src; @@ -52,6 +53,7 @@ pub mod inner; pub mod krylov; pub mod lapack; pub mod layout; +pub mod lobpcg; pub mod norm; pub mod operator; pub mod opnorm; @@ -63,7 +65,6 @@ pub mod svddc; pub mod trace; pub mod triangular; pub mod types; -pub mod lobpcg; pub use assert::*; pub use cholesky::*; diff --git a/src/lobpcg.rs b/src/lobpcg.rs index f1364132..933f2229 100644 --- a/src/lobpcg.rs +++ b/src/lobpcg.rs @@ -1,39 +1,62 @@ -use num_traits::NumCast; +///! Locally Optimal Block Preconditioned Conjugated +///! +///! This module implements the Locally Optimal Block Preconditioned Conjugated (LOBPCG) algorithm, +///which can be used as a solver for large symmetric positive definite eigenproblems. + +use crate::error::{LinalgError, Result}; +use crate::{cholesky::*, close_l2, eigh::*, norm::*, triangular::*}; +use crate::{Lapack, Scalar}; use ndarray::prelude::*; use ndarray::OwnedRepr; -use crate::{cholesky::*, triangular::*, eigh::*, norm::*, close_l2}; -//use sprs::CsMat; -use crate::{Scalar, Lapack}; -use crate::error::{Result, LinalgError}; +use num_traits::NumCast; +/// Find largest or smallest eigenvalues pub enum Order { Largest, - Smallest + Smallest, } +/// The result of the eigensolver +/// +/// In the best case the eigensolver has converged with a result better than the given threshold, +/// then a `EigResult::Ok` gives the eigenvalues, eigenvectors and norms. If an error ocurred +/// during the process, it is returned in `EigResult::Err`, but the best result is still returned, +/// as it could be usable. If there is no result at all, then `EigResult::NoResult` is returned. +/// This happens if the algorithm fails in an early stage, for example if the matrix `A` is not SPD #[derive(Debug)] pub enum EigResult { Ok(Array1, Array2, Vec), Err(Array1, Array2, Vec, LinalgError), - NoResult(LinalgError) + NoResult(LinalgError), } -fn sorted_eig(a: ArrayView2, b: Option>, size: usize, order: &Order) -> Result<(Array1, Array2)> { - //close_l2(&input, &input.t(), 1e-4); - +/// Solve full eigenvalue problem, sort by `order` and truncate to `size` +fn sorted_eig( + a: ArrayView2, + b: Option>, + size: usize, + order: &Order, +) -> Result<(Array1, Array2)> { let (vals, vecs) = match b { Some(b) => (a, b).eigh(UPLO::Upper).map(|x| (x.0, (x.1).0))?, - _ => a.eigh(UPLO::Upper)? + _ => a.eigh(UPLO::Upper)?, }; let n = a.len_of(Axis(0)); Ok(match order { - Order::Largest => (vals.slice_move(s![n-size..; -1]).mapv(|x| Scalar::from_real(x)), vecs.slice_move(s![.., n-size..; -1])), - Order::Smallest => (vals.slice_move(s![..size]).mapv(|x| Scalar::from_real(x)), vecs.slice_move(s![.., ..size])) + Order::Largest => ( + vals.slice_move(s![n-size..; -1]).mapv(|x| Scalar::from_real(x)), + vecs.slice_move(s![.., n-size..; -1]), + ), + Order::Smallest => ( + vals.slice_move(s![..size]).mapv(|x| Scalar::from_real(x)), + vecs.slice_move(s![.., ..size]), + ), }) } +/// Masks a matrix with the given `matrix` fn ndarray_mask(matrix: ArrayView2, mask: &[bool]) -> Array2 { let (rows, cols) = (matrix.nrows(), matrix.ncols()); @@ -41,23 +64,33 @@ fn ndarray_mask(matrix: ArrayView2, mask: &[bool]) -> Array2 { let n_positive = mask.iter().filter(|x| **x).count(); - let matrix = matrix.gencolumns().into_iter().zip(mask.iter()) - .filter(|(_,x)| **x) - .map(|(x,_)| x.to_vec()) + let matrix = matrix + .gencolumns() + .into_iter() + .zip(mask.iter()) + .filter(|(_, x)| **x) + .map(|(x, _)| x.to_vec()) .flatten() .collect::>(); - Array2::from_shape_vec((n_positive, rows), matrix).unwrap().reversed_axes() + Array2::from_shape_vec((n_positive, rows), matrix) + .unwrap() + .reversed_axes() } +/// Applies constraints ensuring that a matrix is orthogonal to it +/// +/// This functions takes a matrix `v` and constraint matrix `y` and orthogonalize the `v` to `y`. fn apply_constraints( mut v: ArrayViewMut, fact_yy: &CholeskyFactorized>, - y: ArrayView2 + y: ArrayView2, ) { let gram_yv = y.t().dot(&v); - let u = gram_yv.genrows().into_iter() + let u = gram_yv + .genrows() + .into_iter() .map(|x| fact_yy.solvec(&x).unwrap().to_vec()) .flatten() .collect::>(); @@ -67,32 +100,57 @@ fn apply_constraints( v -= &(y.dot(&u)); } -fn orthonormalize( - v: Array2 -) -> Result<(Array2, Array2)> { +/// Orthonormalize `V` with Cholesky factorization +/// +/// This also returns the matrix `R` of the `QR` problem +fn orthonormalize(v: Array2) -> Result<(Array2, Array2)> { let gram_vv = v.t().dot(&v); let gram_vv_fac = gram_vv.cholesky(UPLO::Lower)?; - close_l2(&gram_vv, &gram_vv_fac.dot(&gram_vv_fac.t()), NumCast::from(1e-5).unwrap()); + close_l2( + &gram_vv, + &gram_vv_fac.dot(&gram_vv_fac.t()), + NumCast::from(1e-5).unwrap(), + ); let v_t = v.reversed_axes(); - let u = gram_vv_fac.solve_triangular(UPLO::Lower, Diag::NonUnit, &v_t)? + let u = gram_vv_fac + .solve_triangular(UPLO::Lower, Diag::NonUnit, &v_t)? .reversed_axes(); Ok((u, gram_vv_fac)) } +/// Eigenvalue solver for large symmetric positive definite (SPD) eigenproblems +/// +/// # Arguments +/// * `a` - An operator defining the problem, usually a sparse (sometimes also dense) matrix +/// multiplication. Also called the "Stiffness matrix". +/// * `x` - Initial approximation to the k eigenvectors. If `a` has shape=(n,n), then `x` should +/// have shape=(n,k). +/// * `m` - Preconditioner to `a`, by default the identity matrix. In the optimal case `m` +/// approximates the inverse of `a`. +/// * `y` - Constraints of (n,size_y), iterations are performed in the orthogonal complement of the +/// column-space of `y`. It must be full rank. +/// * `tol` - The tolerance values defines at which point the solver stops the optimization. The l2-norm +/// of the residual is compared to this value and the eigenvalue approximation returned if below +/// the threshold. +/// * `maxiter` - The maximal number of iterations +/// * `order` - Whether to solve for the largest or lowest eigenvalues +/// +/// The function returns an `EigResult` with the eigenvalue/eigenvector and achieved residual norm +/// 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) -> Array2>( a: F, x: Array2, m: Option>, y: Option>, - tol: A::Real, maxiter: usize, - order: Order + tol: A::Real, + maxiter: usize, + order: Order, ) -> EigResult { - // the target matrix should be symmetric and quadratic - //assert!(sprs::is_symmetric(&A)); - // the initital approximation should be maximal square // n is the dimensionality of the problem let (n, size_x) = (x.nrows(), x.ncols()); @@ -100,7 +158,7 @@ pub fn lobpcg) -> let size_y = match y { Some(ref y) => y.ncols(), - _ => 0 + _ => 0, }; if (n - size_y) < 5 * size_x { @@ -116,7 +174,7 @@ pub fn lobpcg) -> // orthonormalize the initial guess and calculate matrices AX and XAX let (x, _) = match orthonormalize(x) { Ok(x) => x, - Err(err) => return EigResult::NoResult(err) + Err(err) => return EigResult::NoResult(err), }; let ax = a(x.view()); @@ -125,7 +183,7 @@ pub fn lobpcg) -> // perform eigenvalue decomposition on XAX let (mut lambda, eig_block) = match sorted_eig(xax.view(), None, size_x, &order) { Ok(x) => x, - Err(err) => return EigResult::NoResult(err) + Err(err) => return EigResult::NoResult(err), }; //dbg!(&lambda, &eig_block); @@ -182,13 +240,13 @@ pub fn lobpcg) -> apply_constraints(active_block_r.view_mut(), fact_yy, y.view()); } - let (r,_) = match orthonormalize(active_block_r) { + let (r, _) = match orthonormalize(active_block_r) { Ok(x) => x, - Err(err) => break Err(err) + Err(err) => break Err(err), }; let ar = a(r.view()); - + // perform the Rayleigh Ritz procedure let xaw = x.t().dot(&ar); let waw = r.t().dot(&ar); @@ -203,7 +261,7 @@ pub fn lobpcg) -> //dbg!(&active_P, &P_R); let active_ap = match p_r.solve_triangular(UPLO::Lower, Diag::NonUnit, &active_ap.reversed_axes()) { Ok(x) => x, - Err(err) => break Err(err) + Err(err) => break Err(err), }; let active_ap = active_ap.reversed_axes(); @@ -218,51 +276,46 @@ pub fn lobpcg) -> let wp = r.t().dot(&active_p); ( - stack![Axis(0), + stack![ + Axis(0), stack![Axis(1), Array2::from_diag(&lambda), xaw, xap], stack![Axis(1), xaw.t(), waw, wap], stack![Axis(1), xap.t(), wap.t(), pap] ], - - stack![Axis(0), + stack![ + Axis(0), stack![Axis(1), ident0, xw, xp], stack![Axis(1), xw.t(), ident, wp], stack![Axis(1), xp.t(), wp.t(), ident] ], Some(active_p), - Some(active_ap) + Some(active_ap), ) } else { ( - stack![Axis(0), + stack![ + Axis(0), stack![Axis(1), Array2::from_diag(&lambda), xaw], stack![Axis(1), xaw.t(), waw] ], - stack![Axis(0), - stack![Axis(1), ident0, xw], - stack![Axis(1), xw.t(), ident] - ], + stack![Axis(0), stack![Axis(1), ident0, xw], stack![Axis(1), xw.t(), ident]], + None, None, - None ) }; - //assert!(is_symmetric(gramA.view())); - //assert!(is_symmetric(gramB.view())); - - //dbg!(&gramA, &gramB); let (new_lambda, eig_vecs) = match sorted_eig(gram_a.view(), Some(gram_b.view()), size_x, &order) { Ok(x) => x, - Err(err) => break Err(err) + Err(err) => break Err(err), }; lambda = new_lambda; //dbg!(&lambda, &eig_vecs); - let (pp, app, eig_x) = if let (Some(_), (Some(ref active_p), Some(ref active_ap))) = (ap, (active_p, active_ap)) { - + let (pp, app, eig_x) = if let (Some(_), (Some(ref active_p), Some(ref active_ap))) = (ap, (active_p, active_ap)) + { let eig_x = eig_vecs.slice(s![..size_x, ..]); - let eig_r = eig_vecs.slice(s![size_x..size_x+current_block_size, ..]); - let eig_p = eig_vecs.slice(s![size_x+current_block_size.., ..]); + let eig_r = eig_vecs.slice(s![size_x..size_x + current_block_size, ..]); + let eig_p = eig_vecs.slice(s![size_x + current_block_size.., ..]); //dbg!(&eig_X); //dbg!(&eig_R); @@ -292,23 +345,18 @@ pub fn lobpcg) -> results.push((lambda.clone(), x.clone())); - //dbg!(&X); - //dbg!(&AX); - ap = Some((pp, app)); - //dbg!(&ap); - iter -= 1; }; - let best_idx = residual_norms.iter() - .enumerate() - .min_by(|&(_, item1): &(usize, &Vec), &(_, item2): &(usize, &Vec)| { - let norm1: A::Real = item1.iter().map(|x| (*x)*(*x)).sum(); - let norm2: A::Real = item2.iter().map(|x| (*x)*(*x)).sum(); + let best_idx = residual_norms.iter().enumerate().min_by( + |&(_, item1): &(usize, &Vec), &(_, item2): &(usize, &Vec)| { + let norm1: A::Real = item1.iter().map(|x| (*x) * (*x)).sum(); + let norm2: A::Real = item2.iter().map(|x| (*x) * (*x)).sum(); norm1.partial_cmp(&norm2).unwrap() - }); + }, + ); match best_idx { Some((idx, norms)) => { @@ -317,32 +365,31 @@ pub fn lobpcg) -> match final_norm { Ok(_) => EigResult::Ok(vals, vecs, norms), - Err(err) => EigResult::Err(vals, vecs, norms, err) - } - }, - None => { - match final_norm { - Ok(_) => panic!("Not error available!"), - Err(err) => EigResult::NoResult(err) + Err(err) => EigResult::Err(vals, vecs, norms, err), } } + None => match final_norm { + Ok(_) => panic!("Not error available!"), + Err(err) => EigResult::NoResult(err), + }, } } #[cfg(test)] mod tests { - use super::sorted_eig; - use super::orthonormalize; - use super::ndarray_mask; - use super::Order; use super::lobpcg; + use super::ndarray_mask; + use super::orthonormalize; + use super::sorted_eig; use super::EigResult; + use super::Order; use crate::close_l2; - use ndarray::prelude::*; use crate::qr::*; - use ndarray_rand::RandomExt; + use ndarray::prelude::*; use ndarray_rand::rand_distr::Uniform; - + use ndarray_rand::RandomExt; + + /// Test the `sorted_eigen` function #[test] fn test_sorted_eigen() { let matrix = Array2::random((10, 10), Uniform::new(0., 10.)); @@ -358,6 +405,7 @@ mod tests { close_l2(&matrix, &rec, 1e-5); } + /// Test the masking function #[test] fn test_masking() { let matrix = Array2::random((10, 5), Uniform::new(0., 10.)); @@ -365,6 +413,7 @@ mod tests { close_l2(&masked_matrix.slice(s![.., 2]), &matrix.slice(s![.., 3]), 1e-12); } + /// Test orthonormalization of a random matrix #[test] fn test_orthonormalize() { let matrix = Array2::random((10, 10), Uniform::new(-10., 10.)); @@ -377,31 +426,37 @@ mod tests { // compare returned factorization with QR decomposition let (_, r) = matrix.qr().unwrap(); - close_l2(&r.mapv(|x: f32| x.abs()) , &l.t().mapv(|x| x.abs()), 1e-2); + close_l2(&r.mapv(|x: f32| x.abs()), &l.t().mapv(|x| x.abs()), 1e-2); } + /// Test the eigensolver with a identity matrix problem and a random initial solution #[test] fn test_eigsolver_diag() { - let diag = arr1(&[1.,2.,3.,4.,5.,6.,7.,8.,9.,10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20.]); + let diag = arr1(&[ + 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., + ]); let a = Array2::from_diag(&diag); let x: Array2 = Array2::random((20, 3), Uniform::new(0.0, 1.0)); + // find smallest eigenvalues let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, 20, Order::Smallest); match result { EigResult::Ok(vals, _, _) => close_l2(&vals, &arr1(&[1.0, 2.0, 3.0]), 1e-5), - EigResult::Err(vals, _,_,_) => close_l2(&vals, &arr1(&[1.0, 2.0, 3.0]), 1e-5), - EigResult::NoResult(err) => panic!("Did not converge: {:?}", err) + EigResult::Err(vals, _, _, _) => close_l2(&vals, &arr1(&[1.0, 2.0, 3.0]), 1e-5), + EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), } + // find largest eigenvalues let x: Array2 = Array2::random((20, 3), Uniform::new(0.0, 1.0)); let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, 20, Order::Largest); match result { EigResult::Ok(vals, _, _) => close_l2(&vals, &arr1(&[20.0, 19.0, 18.0]), 1e-5), - EigResult::Err(vals, _,_,_) => close_l2(&vals, &arr1(&[20.0, 19.0, 18.0]), 1e-5), - EigResult::NoResult(err) => panic!("Did not converge: {:?}", err) + EigResult::Err(vals, _, _, _) => close_l2(&vals, &arr1(&[20.0, 19.0, 18.0]), 1e-5), + EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), } } + /// Test the eigensolver with matrix of constructed eigenvalues #[test] fn test_eigsolver() { let n = 50; @@ -412,14 +467,16 @@ mod tests { let t = Array2::from_diag(&Array1::linspace(n as f64, 1.0, n)); let a = v.dot(&t.dot(&v.t())); + // initial random solution let x: Array2 = Array2::random((n, 5), Uniform::new(0.0, 1.0)); + // find five largest eigenvalues let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, 20, Order::Largest); match result { EigResult::Ok(vals, _, _) | EigResult::Err(vals, _, _, _) => { close_l2(&vals, &arr1(&[50.0, 49.0, 48.0, 47.0, 46.0]), 1e-5); - }, - EigResult::NoResult(err) => panic!("Did not converge: {:?}", err) + } + EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), } } } diff --git a/tests/svddc.rs b/tests/svddc.rs index 1dbbf32b..2c9204c8 100644 --- a/tests/svddc.rs +++ b/tests/svddc.rs @@ -14,7 +14,7 @@ fn test(a: &Array2, flag: UVTFlag) { assert!(u.is_none()); assert!(vt.is_none()); return; - }, + } }; let u: Array2<_> = u.unwrap(); let vt: Array2<_> = vt.unwrap(); @@ -53,7 +53,7 @@ macro_rules! test_svd_impl { let a = random(($n, $m).f()); test(&a, UVTFlag::Full); } - + #[test] fn []() { let a = random(($n, $m).f()); From b21b8463491cdb7167b38cdb48f2fe41e6805953 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Mon, 16 Mar 2020 09:21:59 +0100 Subject: [PATCH 06/30] Improve tests and add test for convergence --- src/lobpcg.rs | 76 +++++++++++++++++++++++++++++++-------------------- 1 file changed, 47 insertions(+), 29 deletions(-) diff --git a/src/lobpcg.rs b/src/lobpcg.rs index 933f2229..c528fe81 100644 --- a/src/lobpcg.rs +++ b/src/lobpcg.rs @@ -350,6 +350,7 @@ pub fn lobpcg) -> iter -= 1; }; + dbg!(&residual_norms); let best_idx = residual_norms.iter().enumerate().min_by( |&(_, item1): &(usize, &Vec), &(_, item2): &(usize, &Vec)| { let norm1: A::Real = item1.iter().map(|x| (*x) * (*x)).sum(); @@ -429,6 +430,37 @@ mod tests { close_l2(&r.mapv(|x: f32| x.abs()), &l.t().mapv(|x| x.abs()), 1e-2); } + fn assert_symmetric(a: &Array2) { + close_l2(a, &a.t(), 1e-5); + } + + fn check_eigenvalues(a: &Array2, order: Order, num: usize, ground_truth_eigvals: &[f64]) { + assert_symmetric(a); + + let n = a.len_of(Axis(0)); + let x: Array2 = Array2::random((n, num), Uniform::new(0.0, 1.0)); + + let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, n, order); + match result { + EigResult::Ok(vals, _, r_norms) | EigResult::Err(vals, _, r_norms, _) => { + // check convergence + for (i, norm) in r_norms.into_iter().enumerate() { + if norm > 0.01 { + println!("==== Assertion Failed ===="); + println!("The {} eigenvalue estimation did not converge!", i); + panic!("Too large deviation of residual norm: {} > 0.01", norm); + } + } + + // check correct order of eigenvalues + if ground_truth_eigvals.len() == num { + close_l2(&Array1::from(ground_truth_eigvals.to_vec()), &vals, 5e-2) + } + }, + EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), + } + } + /// Test the eigensolver with a identity matrix problem and a random initial solution #[test] fn test_eigsolver_diag() { @@ -436,47 +468,33 @@ mod tests { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., ]); let a = Array2::from_diag(&diag); - let x: Array2 = Array2::random((20, 3), Uniform::new(0.0, 1.0)); - - // find smallest eigenvalues - let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, 20, Order::Smallest); - match result { - EigResult::Ok(vals, _, _) => close_l2(&vals, &arr1(&[1.0, 2.0, 3.0]), 1e-5), - EigResult::Err(vals, _, _, _) => close_l2(&vals, &arr1(&[1.0, 2.0, 3.0]), 1e-5), - EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), - } - // find largest eigenvalues - let x: Array2 = Array2::random((20, 3), Uniform::new(0.0, 1.0)); - let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, 20, Order::Largest); - match result { - EigResult::Ok(vals, _, _) => close_l2(&vals, &arr1(&[20.0, 19.0, 18.0]), 1e-5), - EigResult::Err(vals, _, _, _) => close_l2(&vals, &arr1(&[20.0, 19.0, 18.0]), 1e-5), - EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), - } + check_eigenvalues(&a, Order::Largest, 3, &[20.,19.,18.]); + check_eigenvalues(&a, Order::Smallest, 3, &[1., 2., 3.]); } /// Test the eigensolver with matrix of constructed eigenvalues #[test] - fn test_eigsolver() { + fn test_eigsolver_constructed() { let n = 50; let tmp = Array2::random((n, n), Uniform::new(0.0, 1.0)); + //let (v, _) = tmp.qr_square().unwrap(); let (v, _) = orthonormalize(tmp).unwrap(); // set eigenvalues in decreasing order - let t = Array2::from_diag(&Array1::linspace(n as f64, 1.0, n)); + let t = Array2::from_diag(&Array1::linspace(n as f64, -(n as f64), n)); let a = v.dot(&t.dot(&v.t())); - // initial random solution - let x: Array2 = Array2::random((n, 5), Uniform::new(0.0, 1.0)); - // find five largest eigenvalues - let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, 20, Order::Largest); - match result { - EigResult::Ok(vals, _, _) | EigResult::Err(vals, _, _, _) => { - close_l2(&vals, &arr1(&[50.0, 49.0, 48.0, 47.0, 46.0]), 1e-5); - } - EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), - } + check_eigenvalues(&a, Order::Largest, 5, &[50.0, 48.0, 46.0, 44.0, 42.0]); + check_eigenvalues(&a, Order::Smallest, 5, &[-50.0, -48.0, -46.0, -44.0, -42.0]); + } + + #[test] + fn test_eigsolver_convergence() { + let tmp = Array2::random((50, 50), Uniform::new(0.0, 1.0)); + let a = tmp.t().dot(&tmp); + + check_eigenvalues(&a, Order::Largest, 5, &[]); } } From 4bb6a53346a2932251b5fd53291a4734552db643 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Mon, 16 Mar 2020 09:50:37 +0100 Subject: [PATCH 07/30] Move algorithm to seperate folder and add test for constraints --- src/{ => lobpcg}/lobpcg.rs | 49 ++++++++++++++++++++++++++++++++++---- 1 file changed, 45 insertions(+), 4 deletions(-) rename src/{ => lobpcg}/lobpcg.rs (90%) diff --git a/src/lobpcg.rs b/src/lobpcg/lobpcg.rs similarity index 90% rename from src/lobpcg.rs rename to src/lobpcg/lobpcg.rs index c528fe81..61f0a4cb 100644 --- a/src/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -95,7 +95,8 @@ fn apply_constraints( .flatten() .collect::>(); - let u = Array2::from_shape_vec((5, 5), u).unwrap(); + let rows = gram_yv.len_of(Axis(0)); + let u = Array2::from_shape_vec((rows, u.len() / rows), u).unwrap(); v -= &(y.dot(&u)); } @@ -144,7 +145,7 @@ fn orthonormalize(v: Array2) -> Result<(Array2, Array2 /// the precision of the matrix is too low (switch from `f32` to `f64` for example). pub fn lobpcg) -> Array2>( a: F, - x: Array2, + mut x: Array2, m: Option>, y: Option>, tol: A::Real, @@ -166,10 +167,20 @@ pub fn lobpcg) -> } // cap the number of iteration - let mut iter = usize::min(n, maxiter); + let mut iter = usize::min(n * 10, maxiter); // factorize yy for later use - let fact_yy = y.as_ref().map(|x| x.t().dot(x).factorizec(UPLO::Upper).unwrap()); + let fact_yy = match y { + Some(ref y) => { + let fact_yy = y.t().dot(y).factorizec(UPLO::Upper).unwrap(); + + apply_constraints(x.view_mut(), &fact_yy, y.view()); + Some(fact_yy) + }, + None => None + }; + + // orthonormalize the initial guess and calculate matrices AX and XAX let (x, _) = match orthonormalize(x) { @@ -497,4 +508,34 @@ mod tests { check_eigenvalues(&a, Order::Largest, 5, &[]); } + + #[test] + fn test_eigsolver_constrainted() { + let diag = arr1(&[ + 1., 2., 3., 4., 5., 6., 7., 8., 9., 10. + ]); + let a = Array2::from_diag(&diag); + let x: Array2 = Array2::random((10, 1), Uniform::new(0.0, 1.0)); + let y: Array2 = 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); + dbg!(&result); + match result { + EigResult::Ok(vals, vecs, r_norms) | EigResult::Err(vals, vecs, r_norms, _) => { + // check convergence + for (i, norm) in r_norms.into_iter().enumerate() { + if norm > 0.01 { + println!("==== Assertion Failed ===="); + println!("The {} eigenvalue estimation did not converge!", i); + panic!("Too large deviation of residual norm: {} > 0.01", norm); + } + } + + // should be the second eigenvalue + close_l2(&vals, &Array1::from(vec![2.0]), 1e-2); + close_l2(&vecs.column(0).mapv(|x| x.abs()), &arr1(&[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), 1e-5); + }, + EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), + } + } } From 90c2c35f65990f1d822a7fe0df34a8e61210d72c Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Mon, 16 Mar 2020 11:04:24 +0100 Subject: [PATCH 08/30] Add truncated eigenvalue module --- Cargo.toml | 2 +- src/lobpcg/lobpcg.rs | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index d8c495ad..eeb2ddf9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -28,6 +28,7 @@ num-traits = "0.2" cauchy = "0.2.1" num-complex = "0.2.1" rand = "0.5" +ndarray-rand = "0.11" [dependencies.ndarray] version = "0.13" @@ -51,7 +52,6 @@ optional = true [dev-dependencies] paste = "0.1" criterion = "0.3" -ndarray-rand = "0.11" [[bench]] name = "eigh" diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 61f0a4cb..ec8724fe 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -11,6 +11,7 @@ use ndarray::OwnedRepr; use num_traits::NumCast; /// Find largest or smallest eigenvalues +#[derive(Debug, Clone)] pub enum Order { Largest, Smallest, From 9c21d78c718e322f294540c6c486c2c250ff5906 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Mon, 16 Mar 2020 12:39:08 +0100 Subject: [PATCH 09/30] Implement iterator for TruncatedEig --- src/lobpcg/lobpcg.rs | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index ec8724fe..19de4388 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -88,16 +88,27 @@ fn apply_constraints( y: ArrayView2, ) { let gram_yv = y.t().dot(&v); + dbg!(&gram_yv.shape()); let u = gram_yv - .genrows() + .gencolumns() .into_iter() - .map(|x| fact_yy.solvec(&x).unwrap().to_vec()) + .map(|x| { + dbg!(&x.shape()); + let res = fact_yy.solvec(&x).unwrap(); + + dbg!(&res); + + res.to_vec() + }) .flatten() .collect::>(); 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)); } @@ -173,7 +184,7 @@ pub fn lobpcg) -> // factorize yy for later use let fact_yy = match y { Some(ref y) => { - let fact_yy = y.t().dot(y).factorizec(UPLO::Upper).unwrap(); + let fact_yy = y.t().dot(y).factorizec(UPLO::Lower).unwrap(); apply_constraints(x.view_mut(), &fact_yy, y.view()); Some(fact_yy) @@ -181,8 +192,6 @@ pub fn lobpcg) -> None => None }; - - // orthonormalize the initial guess and calculate matrices AX and XAX let (x, _) = match orthonormalize(x) { Ok(x) => x, @@ -362,7 +371,7 @@ pub fn lobpcg) -> iter -= 1; }; - dbg!(&residual_norms); + //dbg!(&residual_norms); let best_idx = residual_norms.iter().enumerate().min_by( |&(_, item1): &(usize, &Vec), &(_, item2): &(usize, &Vec)| { let norm1: A::Real = item1.iter().map(|x| (*x) * (*x)).sum(); @@ -452,7 +461,7 @@ mod tests { let n = a.len_of(Axis(0)); let x: Array2 = Array2::random((n, num), Uniform::new(0.0, 1.0)); - let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, n, order); + let result = lobpcg(|y| a.dot(&y), x, None, None, 1e-10, 2 * n, order); match result { EigResult::Ok(vals, _, r_norms) | EigResult::Err(vals, _, r_norms, _) => { // check convergence From 8188649f6f2afc831d7d8ab3407824042017b955 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Mon, 16 Mar 2020 14:03:31 +0100 Subject: [PATCH 10/30] Implement IntoIterator for TruncatedEig --- src/lobpcg/eig.rs | 134 ++++++++++++++++++++++++++++++++++++++++++++++ src/lobpcg/mod.rs | 5 ++ 2 files changed, 139 insertions(+) create mode 100644 src/lobpcg/eig.rs create mode 100644 src/lobpcg/mod.rs diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs new file mode 100644 index 00000000..d169a045 --- /dev/null +++ b/src/lobpcg/eig.rs @@ -0,0 +1,134 @@ +///! Implements truncated eigenvalue decomposition +/// + +use ndarray::prelude::*; +use ndarray::stack; +use ndarray_rand::rand_distr::Uniform; +use ndarray_rand::RandomExt; +use num_traits::{Float, NumCast}; +use crate::{Scalar, Lapack}; +use super::lobpcg::{lobpcg, EigResult, Order}; + +pub struct TruncatedEig { + order: Order, + problem: Array2, + pub constraints: Option>, + precision: A::Real, + maxiter: usize +} + +impl TruncatedEig { + pub fn new(problem: Array2, order: Order) -> TruncatedEig { + TruncatedEig { + precision: NumCast::from(1e-5).unwrap(), + maxiter: problem.len_of(Axis(0)) * 2, + constraints: None, + order, + problem + } + } + + pub fn precision(mut self, precision: A::Real) -> Self { + self.precision = precision; + + self + } + + pub fn maxiter(mut self, maxiter: usize) -> Self { + self.maxiter = maxiter; + + self + + } + + pub fn constraints(mut self, constraints: Array2) -> Self { + self.constraints = Some(constraints); + + self + } + + pub fn once(&self, num: usize) -> EigResult { + 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, None, self.constraints.clone(), self.precision, self.maxiter, self.order.clone()) + } +} + +impl IntoIterator for TruncatedEig { + type Item = (Array1, Array2); + type IntoIter = TruncatedEigIterator; + + fn into_iter(self) -> TruncatedEigIterator{ + TruncatedEigIterator { + step_size: 1, + eig: self + } + } +} + +pub struct TruncatedEigIterator { + step_size: usize, + eig: TruncatedEig +} + +impl Iterator for TruncatedEigIterator { + type Item = (Array1, Array2); + + fn next(&mut self) -> Option { + let res = self.eig.once(self.step_size); + dbg!(&res); + + match res { + EigResult::Ok(vals, vecs, norms) | EigResult::Err(vals, vecs, norms, _) => { + // abort if any eigenproblem did not converge + for r_norm in norms { + if r_norm > NumCast::from(0.1).unwrap() { + return None; + } + } + + let new_constraints = if let Some(ref constraints) = self.eig.constraints { + let eigvecs_arr = constraints.gencolumns().into_iter() + .chain(vecs.gencolumns().into_iter()) + .map(|x| x.insert_axis(Axis(1))) + .collect::>(); + + stack(Axis(1), &eigvecs_arr).unwrap() + } else { + vecs.clone() + }; + + dbg!(&new_constraints); + + self.eig.constraints = Some(new_constraints); + + Some((vals, vecs)) + }, + EigResult::NoResult(_) => None + } + } +} + +#[cfg(test)] +mod tests { + use super::TruncatedEig; + use super::Order; + use ndarray::Array2; + use ndarray_rand::rand_distr::Uniform; + use ndarray_rand::RandomExt; + + #[test] + fn test_truncated_eig() { + let a = Array2::random((50, 50), Uniform::new(0., 1.0)); + let a = a.t().dot(&a); + + let teig = TruncatedEig::new(a, Order::Largest) + .precision(1e-5) + .maxiter(500); + + let res = teig.into_iter().take(3).flat_map(|x| x.0.to_vec()).collect::>(); + dbg!(&res); + panic!(""); + } +} diff --git a/src/lobpcg/mod.rs b/src/lobpcg/mod.rs new file mode 100644 index 00000000..a82919c4 --- /dev/null +++ b/src/lobpcg/mod.rs @@ -0,0 +1,5 @@ +mod lobpcg; +mod eig; + +pub use lobpcg::{lobpcg, EigResult, Order}; +pub use eig::TruncatedEig; From 705d37abaf1e17b6befd4d4e08b27e6edf186a7c Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Mon, 16 Mar 2020 18:03:51 +0100 Subject: [PATCH 11/30] Add truncated singular value decomposition module --- src/lobpcg/eig.rs | 42 ++++++++---- src/lobpcg/lobpcg.rs | 8 --- src/lobpcg/mod.rs | 2 + src/lobpcg/svd.rs | 154 +++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 186 insertions(+), 20 deletions(-) create mode 100644 src/lobpcg/svd.rs diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs index d169a045..96ff44ea 100644 --- a/src/lobpcg/eig.rs +++ b/src/lobpcg/eig.rs @@ -9,10 +9,17 @@ use num_traits::{Float, NumCast}; use crate::{Scalar, Lapack}; use super::lobpcg::{lobpcg, EigResult, Order}; +/// Truncated eigenproblem solver +/// +/// 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. pub struct TruncatedEig { order: Order, problem: Array2, pub constraints: Option>, + preconditioner: Option>, precision: A::Real, maxiter: usize } @@ -22,6 +29,7 @@ impl TruncatedEig { TruncatedEig { precision: NumCast::from(1e-5).unwrap(), maxiter: problem.len_of(Axis(0)) * 2, + preconditioner: None, constraints: None, order, problem @@ -41,17 +49,24 @@ impl TruncatedEig { } - pub fn constraints(mut self, constraints: Array2) -> Self { + pub fn orthogonal_to(mut self, constraints: Array2) -> Self { self.constraints = Some(constraints); self } + pub fn precondition_with(mut self, preconditioner: Array2) -> Self { + self.preconditioner = Some(preconditioner); + + self + } + + // calculate the eigenvalues once pub fn once(&self, num: usize) -> EigResult { 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, None, self.constraints.clone(), self.precision, self.maxiter, self.order.clone()) + lobpcg(|y| self.problem.dot(&y), x, self.preconditioner.clone(), self.constraints.clone(), self.precision, self.maxiter, self.order.clone()) } } @@ -67,6 +82,10 @@ impl IntoIterator for Truncat } } +/// Truncate 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. pub struct TruncatedEigIterator { step_size: usize, eig: TruncatedEig @@ -77,7 +96,6 @@ impl Iterator for TruncatedEi fn next(&mut self) -> Option { let res = self.eig.once(self.step_size); - dbg!(&res); match res { EigResult::Ok(vals, vecs, norms) | EigResult::Err(vals, vecs, norms, _) => { @@ -88,6 +106,7 @@ impl Iterator for TruncatedEi } } + // add the new eigenvector to the internal constrain matrix let new_constraints = if let Some(ref constraints) = self.eig.constraints { let eigvecs_arr = constraints.gencolumns().into_iter() .chain(vecs.gencolumns().into_iter()) @@ -99,8 +118,6 @@ impl Iterator for TruncatedEi vecs.clone() }; - dbg!(&new_constraints); - self.eig.constraints = Some(new_constraints); Some((vals, vecs)) @@ -114,21 +131,22 @@ impl Iterator for TruncatedEi mod tests { use super::TruncatedEig; use super::Order; - use ndarray::Array2; - use ndarray_rand::rand_distr::Uniform; - use ndarray_rand::RandomExt; + use ndarray::{arr1, Array2}; #[test] fn test_truncated_eig() { - let a = Array2::random((50, 50), Uniform::new(0., 1.0)); - let a = a.t().dot(&a); + let diag = arr1(&[ + 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., + ]); + let a = Array2::from_diag(&diag); let teig = TruncatedEig::new(a, Order::Largest) .precision(1e-5) .maxiter(500); let res = teig.into_iter().take(3).flat_map(|x| x.0.to_vec()).collect::>(); - dbg!(&res); - panic!(""); + let ground_truth = vec![20., 19., 18.]; + + assert!(ground_truth.into_iter().zip(res.into_iter()).map(|(x,y)| (x-y)*(x-y)).sum::() < 0.01); } } diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 19de4388..599777fe 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -511,14 +511,6 @@ mod tests { check_eigenvalues(&a, Order::Smallest, 5, &[-50.0, -48.0, -46.0, -44.0, -42.0]); } - #[test] - fn test_eigsolver_convergence() { - let tmp = Array2::random((50, 50), Uniform::new(0.0, 1.0)); - let a = tmp.t().dot(&tmp); - - check_eigenvalues(&a, Order::Largest, 5, &[]); - } - #[test] fn test_eigsolver_constrainted() { let diag = arr1(&[ diff --git a/src/lobpcg/mod.rs b/src/lobpcg/mod.rs index a82919c4..854d8f5f 100644 --- a/src/lobpcg/mod.rs +++ b/src/lobpcg/mod.rs @@ -1,5 +1,7 @@ mod lobpcg; mod eig; +mod svd; pub use lobpcg::{lobpcg, EigResult, Order}; pub use eig::TruncatedEig; +pub use svd::TruncatedSvd; diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs new file mode 100644 index 00000000..00b4707f --- /dev/null +++ b/src/lobpcg/svd.rs @@ -0,0 +1,154 @@ +///! Implements truncated singular value decomposition +/// + +use std::ops::DivAssign; +use ndarray::prelude::*; +use ndarray::stack; +use ndarray_rand::rand_distr::Uniform; +use ndarray_rand::RandomExt; +use num_traits::{Float, NumCast}; +use crate::{Scalar, Lapack}; +use super::lobpcg::{lobpcg, EigResult, Order}; +use crate::error::Result; + +#[derive(Debug)] +pub struct TruncatedSvdResult { + eigvals: Array1, + eigvecs: Array2, + problem: Array2, + ngm: bool +} + +impl + 'static> TruncatedSvdResult { + fn singular_values_with_indices(&self) -> (Vec, Vec) { + let mut a = self.eigvals.iter() + .map(|x| if *x < NumCast::from(1e-5).unwrap() { NumCast::from(0.0).unwrap() } else { *x }) + .map(|x| x.sqrt()) + .enumerate() + .collect::>(); + + a.sort_by(|(_,x), (_, y)| x.partial_cmp(&y).unwrap().reverse()); + + a.into_iter().map(|(a,b)| (b,a)).unzip() + } + + pub fn values(&self) -> Vec { + let (values, indices) = self.singular_values_with_indices(); + + values + } + + pub fn values_vecs(&self) -> (Array2, Vec, Array2) { + let (values, indices) = self.singular_values_with_indices(); + let n_values = values.iter().filter(|x| **x > NumCast::from(0.0).unwrap()).count(); + + if self.ngm { + let vlarge = self.eigvecs.select(Axis(1), &indices); + let mut ularge = self.problem.dot(&vlarge); + + ularge.gencolumns_mut().into_iter() + .zip(values.iter()) + .for_each(|(mut a,b)| a.mapv_inplace(|x| x / *b)); + + let vhlarge = vlarge.reversed_axes(); + + (vhlarge, values, ularge) + } else { + let ularge = self.eigvecs.select(Axis(1), &indices); + + let mut vlarge = ularge.dot(&self.problem); + vlarge.gencolumns_mut().into_iter() + .zip(values.iter()) + .for_each(|(mut a,b)| a.mapv_inplace(|x| x / *b)); + let vhlarge = vlarge.reversed_axes(); + + (vhlarge, values, ularge) + } + } +} + +/// 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. +pub struct TruncatedSvd { + order: Order, + problem: Array2, + precision: A::Real, + maxiter: usize +} + +impl TruncatedSvd { + pub fn new(problem: Array2, order: Order) -> TruncatedSvd { + TruncatedSvd { + precision: NumCast::from(1e-5).unwrap(), + maxiter: problem.len_of(Axis(0)) * 2, + order, + problem + } + } + + pub fn precision(mut self, precision: A::Real) -> Self { + self.precision = precision; + + self + } + + pub fn maxiter(mut self, maxiter: usize) -> Self { + self.maxiter = maxiter; + + self + + } + + // calculate the eigenvalues once + pub fn once(&self, num: usize) -> Result> { + let (n,m) = (self.problem.rows(), self.problem.ncols()); + + let x = Array2::random((usize::min(n,m), num), Uniform::new(0.0, 1.0)) + .mapv(|x| NumCast::from(x).unwrap()); + + let res = if n > m { + lobpcg(|y| self.problem.t().dot(&self.problem.dot(&y)), x, None, None, self.precision, self.maxiter, self.order.clone()) + } else { + lobpcg(|y| self.problem.dot(&self.problem.t().dot(&y)), x, None, None, self.precision, self.maxiter, self.order.clone()) + }; + + match res { + EigResult::Ok(vals, vecs, _) | EigResult::Err(vals, vecs, _, _) => { + Ok(TruncatedSvdResult { + problem: self.problem.clone(), + eigvals: vals, + eigvecs: vecs, + ngm: n > m + }) + }, + _ => panic!("") + } + } +} + +#[cfg(test)] +mod tests { + use super::TruncatedSvd; + use super::Order; + use ndarray::{arr1, Array2}; + + #[test] + fn test_truncated_svd() { + let diag = arr1(&[ + 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., + ]); + let a = Array2::from_diag(&diag); + + let res = TruncatedSvd::new(a, Order::Largest) + .precision(1e-5) + .maxiter(500) + .once(3) + .unwrap(); + + dbg!(&res.values()); + } +} From b0e1a40d44a72b8d760ce4b20bd507b802f15f44 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Tue, 17 Mar 2020 12:51:19 +0100 Subject: [PATCH 12/30] Improve test for truncated SVD --- src/lobpcg/lobpcg.rs | 4 ++-- src/lobpcg/svd.rs | 52 +++++++++++++++++++++++--------------------- 2 files changed, 29 insertions(+), 27 deletions(-) diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 599777fe..48e3d84d 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -169,14 +169,14 @@ pub fn lobpcg) -> let (n, size_x) = (x.nrows(), x.ncols()); assert!(size_x <= n); - let size_y = match y { + /*let size_y = match y { Some(ref y) => y.ncols(), _ => 0, }; if (n - size_y) < 5 * size_x { panic!("Please use a different approach, the LOBPCG method only supports the calculation of a couple of eigenvectors!"); - } + }*/ // cap the number of iteration let mut iter = usize::min(n * 10, maxiter); diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index 00b4707f..97d6d2c6 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -3,7 +3,6 @@ use std::ops::DivAssign; use ndarray::prelude::*; -use ndarray::stack; use ndarray_rand::rand_distr::Uniform; use ndarray_rand::RandomExt; use num_traits::{Float, NumCast}; @@ -20,29 +19,32 @@ pub struct TruncatedSvdResult { } impl + 'static> TruncatedSvdResult { - fn singular_values_with_indices(&self) -> (Vec, Vec) { + fn singular_values_with_indices(&self) -> (Array1, Vec) { let mut a = self.eigvals.iter() - .map(|x| if *x < NumCast::from(1e-5).unwrap() { NumCast::from(0.0).unwrap() } else { *x }) .map(|x| x.sqrt()) .enumerate() .collect::>(); a.sort_by(|(_,x), (_, y)| x.partial_cmp(&y).unwrap().reverse()); - a.into_iter().map(|(a,b)| (b,a)).unzip() + let (values, indices): (Vec, Vec) = a.into_iter() + .filter(|(_,x)| *x > NumCast::from(1e-5).unwrap()) + .map(|(a,b)| (b,a)) + .unzip(); + + (Array1::from(values), indices) } - pub fn values(&self) -> Vec { - let (values, indices) = self.singular_values_with_indices(); + pub fn values(&self) -> Array1 { + let (values, _) = self.singular_values_with_indices(); values } - pub fn values_vecs(&self) -> (Array2, Vec, Array2) { + pub fn values_vecs(&self) -> (Array2, Array1, Array2) { let (values, indices) = self.singular_values_with_indices(); - let n_values = values.iter().filter(|x| **x > NumCast::from(0.0).unwrap()).count(); - if self.ngm { + let (u, v) = if self.ngm { let vlarge = self.eigvecs.select(Axis(1), &indices); let mut ularge = self.problem.dot(&vlarge); @@ -50,20 +52,19 @@ impl + 'static> TruncatedSvdResult { .zip(values.iter()) .for_each(|(mut a,b)| a.mapv_inplace(|x| x / *b)); - let vhlarge = vlarge.reversed_axes(); - - (vhlarge, values, ularge) + (ularge, vlarge) } else { let ularge = self.eigvecs.select(Axis(1), &indices); - let mut vlarge = ularge.dot(&self.problem); + let mut vlarge = self.problem.t().dot(&ularge); vlarge.gencolumns_mut().into_iter() .zip(values.iter()) .for_each(|(mut a,b)| a.mapv_inplace(|x| x / *b)); - let vhlarge = vlarge.reversed_axes(); - (vhlarge, values, ularge) - } + (ularge, vlarge) + }; + + (u, values, v.reversed_axes()) } } @@ -105,7 +106,7 @@ impl TruncatedSvd { // calculate the eigenvalues once pub fn once(&self, num: usize) -> Result> { - let (n,m) = (self.problem.rows(), self.problem.ncols()); + 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()); @@ -125,30 +126,31 @@ impl TruncatedSvd { ngm: n > m }) }, - _ => panic!("") + EigResult::NoResult(err) => Err(err) } } } #[cfg(test)] mod tests { + use crate::close_l2; use super::TruncatedSvd; use super::Order; - use ndarray::{arr1, Array2}; + use ndarray::{arr1, arr2}; #[test] fn test_truncated_svd() { - let diag = arr1(&[ - 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., - ]); - let a = Array2::from_diag(&diag); + let a = arr2(&[[3., 2., 2.], + [2., 3., -2.]]); let res = TruncatedSvd::new(a, Order::Largest) .precision(1e-5) .maxiter(500) - .once(3) + .once(2) .unwrap(); - dbg!(&res.values()); + let (u, sigma, v_t) = res.values_vecs(); + + close_l2(&sigma, &arr1(&[5.0, 3.0]), 1e-5); } } From 33340341879345fdcb1564aef591ece8fb1d5f20 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Tue, 17 Mar 2020 14:19:16 +0100 Subject: [PATCH 13/30] Add test for random matrix reconstruction with SVD --- src/lobpcg/svd.rs | 46 ++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 40 insertions(+), 6 deletions(-) diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index 97d6d2c6..6d65d277 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -10,6 +10,10 @@ use crate::{Scalar, Lapack}; use super::lobpcg::{lobpcg, EigResult, Order}; use crate::error::Result; +/// The result of an eigenvalue decomposition for SVD +/// +/// Provides methods for either calculating just the singular values with reduced cost or the +/// vectors as well. #[derive(Debug)] pub struct TruncatedSvdResult { eigvals: Array1, @@ -19,14 +23,18 @@ pub struct TruncatedSvdResult { } impl + 'static> TruncatedSvdResult { + /// Returns singular values ordered by magnitude with indices. fn singular_values_with_indices(&self) -> (Array1, Vec) { + // numerate and square root eigenvalues let mut a = self.eigvals.iter() .map(|x| x.sqrt()) .enumerate() .collect::>(); + // sort by magnitude a.sort_by(|(_,x), (_, y)| x.partial_cmp(&y).unwrap().reverse()); + // filter low singular values away let (values, indices): (Vec, Vec) = a.into_iter() .filter(|(_,x)| *x > NumCast::from(1e-5).unwrap()) .map(|(a,b)| (b,a)) @@ -35,15 +43,18 @@ impl + 'static> TruncatedSvdResult { (Array1::from(values), indices) } + /// Returns singular values orderd by magnitude pub fn values(&self) -> Array1 { let (values, _) = self.singular_values_with_indices(); values } + /// Returns singular values, left-singular vectors and right-singular vectors pub fn values_vecs(&self) -> (Array2, Array1, Array2) { let (values, indices) = self.singular_values_with_indices(); + // branch n > m (for A is [n x m]) let (u, v) = if self.ngm { let vlarge = self.eigvecs.select(Axis(1), &indices); let mut ularge = self.problem.dot(&vlarge); @@ -105,18 +116,23 @@ impl TruncatedSvd { } // calculate the eigenvalues once - pub fn once(&self, num: usize) -> Result> { + pub fn once(self, num: usize) -> Result> { 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()); + // square precision because the SVD squares the eigenvalue as well + let precision = self.precision * self.precision; + + // use problem definition with less operations required let res = if n > m { - lobpcg(|y| self.problem.t().dot(&self.problem.dot(&y)), x, None, None, self.precision, self.maxiter, self.order.clone()) + lobpcg(|y| self.problem.t().dot(&self.problem.dot(&y)), x, None, None, precision, self.maxiter, self.order.clone()) } else { - lobpcg(|y| self.problem.dot(&self.problem.t().dot(&y)), x, None, None, self.precision, self.maxiter, self.order.clone()) + lobpcg(|y| self.problem.dot(&self.problem.t().dot(&y)), x, None, None, precision, self.maxiter, self.order.clone()) }; + // convert into TruncatedSvdResult match res { EigResult::Ok(vals, vecs, _) | EigResult::Err(vals, vecs, _, _) => { Ok(TruncatedSvdResult { @@ -136,7 +152,9 @@ mod tests { use crate::close_l2; use super::TruncatedSvd; use super::Order; - use ndarray::{arr1, arr2}; + use ndarray::{arr1, arr2, Array2}; + use ndarray_rand::rand_distr::Uniform; + use ndarray_rand::RandomExt; #[test] fn test_truncated_svd() { @@ -145,12 +163,28 @@ mod tests { let res = TruncatedSvd::new(a, Order::Largest) .precision(1e-5) - .maxiter(500) + .maxiter(10) .once(2) .unwrap(); - let (u, sigma, v_t) = res.values_vecs(); + let (_, sigma, _) = res.values_vecs(); close_l2(&sigma, &arr1(&[5.0, 3.0]), 1e-5); } + + #[test] + fn test_truncated_svd_random() { + let a: Array2 = Array2::random((50, 10), Uniform::new(0.0, 1.0)); + + let res = TruncatedSvd::new(a.clone(), Order::Largest) + .precision(1e-5) + .maxiter(10) + .once(10) + .unwrap(); + + let (u, sigma, v_t) = res.values_vecs(); + let reconstructed = u.dot(&Array2::from_diag(&sigma).dot(&v_t)); + + close_l2(&a, &reconstructed, 1e-5); + } } From 4acb85ca7199135b593fd2b78b57429430487132 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Tue, 17 Mar 2020 14:38:34 +0100 Subject: [PATCH 14/30] Add example for truncated SVD --- examples/eigh.rs | 1 + src/lib.rs | 1 + src/lobpcg/eig.rs | 6 +++--- src/lobpcg/mod.rs | 2 +- src/lobpcg/svd.rs | 12 ++++++------ 5 files changed, 12 insertions(+), 10 deletions(-) diff --git a/examples/eigh.rs b/examples/eigh.rs index c9bcc941..267ab303 100644 --- a/examples/eigh.rs +++ b/examples/eigh.rs @@ -12,3 +12,4 @@ fn main() { let av = a.dot(&vecs); println!("AV = \n{:?}", av); } + diff --git a/src/lib.rs b/src/lib.rs index 8eac92ee..b5cf1031 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -74,6 +74,7 @@ pub use eigh::*; pub use generate::*; pub use inner::*; pub use layout::*; +pub use lobpcg::{TruncatedEig, TruncatedSvd, TruncatedOrder}; pub use norm::*; pub use operator::*; pub use opnorm::*; diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs index 96ff44ea..038c817e 100644 --- a/src/lobpcg/eig.rs +++ b/src/lobpcg/eig.rs @@ -61,8 +61,8 @@ impl TruncatedEig { self } - // calculate the eigenvalues once - pub fn once(&self, num: usize) -> EigResult { + // calculate the eigenvalues decompose + pub fn decompose(&self, num: usize) -> EigResult { let x = Array2::random((self.problem.len_of(Axis(0)), num), Uniform::new(0.0, 1.0)) .mapv(|x| NumCast::from(x).unwrap()); @@ -95,7 +95,7 @@ impl Iterator for TruncatedEi type Item = (Array1, Array2); fn next(&mut self) -> Option { - let res = self.eig.once(self.step_size); + let res = self.eig.decompose(self.step_size); match res { EigResult::Ok(vals, vecs, norms) | EigResult::Err(vals, vecs, norms, _) => { diff --git a/src/lobpcg/mod.rs b/src/lobpcg/mod.rs index 854d8f5f..013c70f2 100644 --- a/src/lobpcg/mod.rs +++ b/src/lobpcg/mod.rs @@ -2,6 +2,6 @@ mod lobpcg; mod eig; mod svd; -pub use lobpcg::{lobpcg, EigResult, Order}; +pub use lobpcg::{lobpcg, EigResult, Order as TruncatedOrder}; pub use eig::TruncatedEig; pub use svd::TruncatedSvd; diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index 6d65d277..daa74c0d 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -51,7 +51,7 @@ impl + 'static> TruncatedSvdResult { } /// Returns singular values, left-singular vectors and right-singular vectors - pub fn values_vecs(&self) -> (Array2, Array1, Array2) { + pub fn values_vectors(&self) -> (Array2, Array1, Array2) { let (values, indices) = self.singular_values_with_indices(); // branch n > m (for A is [n x m]) @@ -115,8 +115,8 @@ impl TruncatedSvd { } - // calculate the eigenvalues once - pub fn once(self, num: usize) -> Result> { + // calculate the eigenvalue decomposition + pub fn decompose(self, num: usize) -> Result> { let (n,m) = (self.problem.nrows(), self.problem.ncols()); let x = Array2::random((usize::min(n,m), num), Uniform::new(0.0, 1.0)) @@ -164,7 +164,7 @@ mod tests { let res = TruncatedSvd::new(a, Order::Largest) .precision(1e-5) .maxiter(10) - .once(2) + .decompose(2) .unwrap(); let (_, sigma, _) = res.values_vecs(); @@ -179,10 +179,10 @@ mod tests { let res = TruncatedSvd::new(a.clone(), Order::Largest) .precision(1e-5) .maxiter(10) - .once(10) + .decompose(10) .unwrap(); - let (u, sigma, v_t) = res.values_vecs(); + let (u, sigma, v_t) = res.values_vectors(); let reconstructed = u.dot(&Array2::from_diag(&sigma).dot(&v_t)); close_l2(&a, &reconstructed, 1e-5); From 13cf2b3b8d877f94083e32fd97c5f4324b95f4a9 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Tue, 17 Mar 2020 14:40:05 +0100 Subject: [PATCH 15/30] Run cargo fmt again --- examples/eigh.rs | 1 - src/lib.rs | 2 +- src/lobpcg/eig.rs | 53 +++++++++++++-------- src/lobpcg/lobpcg.rs | 27 ++++++----- src/lobpcg/mod.rs | 4 +- src/lobpcg/svd.rs | 110 ++++++++++++++++++++++++------------------- 6 files changed, 111 insertions(+), 86 deletions(-) diff --git a/examples/eigh.rs b/examples/eigh.rs index 267ab303..c9bcc941 100644 --- a/examples/eigh.rs +++ b/examples/eigh.rs @@ -12,4 +12,3 @@ fn main() { let av = a.dot(&vecs); println!("AV = \n{:?}", av); } - diff --git a/src/lib.rs b/src/lib.rs index b5cf1031..472cc348 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -74,7 +74,7 @@ pub use eigh::*; pub use generate::*; pub use inner::*; pub use layout::*; -pub use lobpcg::{TruncatedEig, TruncatedSvd, TruncatedOrder}; +pub use lobpcg::{TruncatedEig, TruncatedOrder, TruncatedSvd}; pub use norm::*; pub use operator::*; pub use opnorm::*; diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs index 038c817e..8b3701eb 100644 --- a/src/lobpcg/eig.rs +++ b/src/lobpcg/eig.rs @@ -1,13 +1,12 @@ +use super::lobpcg::{lobpcg, EigResult, Order}; +use crate::{Lapack, Scalar}; ///! Implements truncated eigenvalue decomposition /// - use ndarray::prelude::*; use ndarray::stack; use ndarray_rand::rand_distr::Uniform; use ndarray_rand::RandomExt; use num_traits::{Float, NumCast}; -use crate::{Scalar, Lapack}; -use super::lobpcg::{lobpcg, EigResult, Order}; /// Truncated eigenproblem solver /// @@ -21,7 +20,7 @@ pub struct TruncatedEig { pub constraints: Option>, preconditioner: Option>, precision: A::Real, - maxiter: usize + maxiter: usize, } impl TruncatedEig { @@ -31,8 +30,8 @@ impl TruncatedEig { maxiter: problem.len_of(Axis(0)) * 2, preconditioner: None, constraints: None, - order, - problem + order, + problem, } } @@ -46,7 +45,6 @@ impl TruncatedEig { self.maxiter = maxiter; self - } pub fn orthogonal_to(mut self, constraints: Array2) -> Self { @@ -66,7 +64,15 @@ impl TruncatedEig { 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()) + lobpcg( + |y| self.problem.dot(&y), + x, + self.preconditioner.clone(), + self.constraints.clone(), + self.precision, + self.maxiter, + self.order.clone(), + ) } } @@ -74,10 +80,10 @@ impl IntoIterator for Truncat type Item = (Array1, Array2); type IntoIter = TruncatedEigIterator; - fn into_iter(self) -> TruncatedEigIterator{ + fn into_iter(self) -> TruncatedEigIterator { TruncatedEigIterator { step_size: 1, - eig: self + eig: self, } } } @@ -88,7 +94,7 @@ impl IntoIterator for Truncat /// eigenvalue/vector pair. Useful for generating pairs until a certain condition is met. pub struct TruncatedEigIterator { step_size: usize, - eig: TruncatedEig + eig: TruncatedEig, } impl Iterator for TruncatedEigIterator { @@ -108,7 +114,9 @@ impl Iterator for TruncatedEi // add the new eigenvector to the internal constrain matrix let new_constraints = if let Some(ref constraints) = self.eig.constraints { - let eigvecs_arr = constraints.gencolumns().into_iter() + let eigvecs_arr = constraints + .gencolumns() + .into_iter() .chain(vecs.gencolumns().into_iter()) .map(|x| x.insert_axis(Axis(1))) .collect::>(); @@ -121,16 +129,16 @@ impl Iterator for TruncatedEi self.eig.constraints = Some(new_constraints); Some((vals, vecs)) - }, - EigResult::NoResult(_) => None + } + EigResult::NoResult(_) => None, } } } #[cfg(test)] mod tests { - use super::TruncatedEig; use super::Order; + use super::TruncatedEig; use ndarray::{arr1, Array2}; #[test] @@ -140,13 +148,18 @@ mod tests { ]); let a = Array2::from_diag(&diag); - let teig = TruncatedEig::new(a, Order::Largest) - .precision(1e-5) - .maxiter(500); - + let teig = TruncatedEig::new(a, Order::Largest).precision(1e-5).maxiter(500); + let res = teig.into_iter().take(3).flat_map(|x| x.0.to_vec()).collect::>(); let ground_truth = vec![20., 19., 18.]; - assert!(ground_truth.into_iter().zip(res.into_iter()).map(|(x,y)| (x-y)*(x-y)).sum::() < 0.01); + assert!( + ground_truth + .into_iter() + .zip(res.into_iter()) + .map(|(x, y)| (x - y) * (x - y)) + .sum::() + < 0.01 + ); } } diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 48e3d84d..bbf7d713 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -1,8 +1,7 @@ ///! Locally Optimal Block Preconditioned Conjugated ///! ///! This module implements the Locally Optimal Block Preconditioned Conjugated (LOBPCG) algorithm, -///which can be used as a solver for large symmetric positive definite eigenproblems. - +///which can be used as a solver for large symmetric positive definite eigenproblems. use crate::error::{LinalgError, Result}; use crate::{cholesky::*, close_l2, eigh::*, norm::*, triangular::*}; use crate::{Lapack, Scalar}; @@ -145,7 +144,7 @@ fn orthonormalize(v: Array2) -> Result<(Array2, Array2 /// approximates the inverse of `a`. /// * `y` - Constraints of (n,size_y), iterations are performed in the orthogonal complement of the /// column-space of `y`. It must be full rank. -/// * `tol` - The tolerance values defines at which point the solver stops the optimization. The l2-norm +/// * `tol` - The tolerance values defines at which point the solver stops the optimization. The l2-norm /// of the residual is compared to this value and the eigenvalue approximation returned if below /// the threshold. /// * `maxiter` - The maximal number of iterations @@ -188,8 +187,8 @@ pub fn lobpcg) -> apply_constraints(x.view_mut(), &fact_yy, y.view()); Some(fact_yy) - }, - None => None + } + None => None, }; // orthonormalize the initial guess and calculate matrices AX and XAX @@ -477,7 +476,7 @@ mod tests { if ground_truth_eigvals.len() == num { close_l2(&Array1::from(ground_truth_eigvals.to_vec()), &vals, 5e-2) } - }, + } EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), } } @@ -490,7 +489,7 @@ mod tests { ]); let a = Array2::from_diag(&diag); - check_eigenvalues(&a, Order::Largest, 3, &[20.,19.,18.]); + check_eigenvalues(&a, Order::Largest, 3, &[20., 19., 18.]); check_eigenvalues(&a, Order::Smallest, 3, &[1., 2., 3.]); } @@ -513,12 +512,10 @@ mod tests { #[test] fn test_eigsolver_constrainted() { - let diag = arr1(&[ - 1., 2., 3., 4., 5., 6., 7., 8., 9., 10. - ]); + let diag = arr1(&[1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]); let a = Array2::from_diag(&diag); let x: Array2 = Array2::random((10, 1), Uniform::new(0.0, 1.0)); - let y: Array2 = arr2(&[[1.0,0.,0.,0.,0.,0.,0.,0.,0.,0.]]).reversed_axes(); + let y: Array2 = 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); dbg!(&result); @@ -535,8 +532,12 @@ mod tests { // should be the second eigenvalue close_l2(&vals, &Array1::from(vec![2.0]), 1e-2); - close_l2(&vecs.column(0).mapv(|x| x.abs()), &arr1(&[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), 1e-5); - }, + close_l2( + &vecs.column(0).mapv(|x| x.abs()), + &arr1(&[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), + 1e-5, + ); + } EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), } } diff --git a/src/lobpcg/mod.rs b/src/lobpcg/mod.rs index 013c70f2..717cbc76 100644 --- a/src/lobpcg/mod.rs +++ b/src/lobpcg/mod.rs @@ -1,7 +1,7 @@ -mod lobpcg; mod eig; +mod lobpcg; mod svd; -pub use lobpcg::{lobpcg, EigResult, Order as TruncatedOrder}; pub use eig::TruncatedEig; +pub use lobpcg::{lobpcg, EigResult, Order as TruncatedOrder}; pub use svd::TruncatedSvd; diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index daa74c0d..6d1a758f 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -1,43 +1,40 @@ -///! Implements truncated singular value decomposition -/// - -use std::ops::DivAssign; +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}; -use crate::{Scalar, Lapack}; -use super::lobpcg::{lobpcg, EigResult, Order}; -use crate::error::Result; +///! Implements truncated singular value decomposition +/// +use std::ops::DivAssign; /// The result of an eigenvalue decomposition for SVD /// /// Provides methods for either calculating just the singular values with reduced cost or the -/// vectors as well. +/// vectors as well. #[derive(Debug)] pub struct TruncatedSvdResult { eigvals: Array1, eigvecs: Array2, problem: Array2, - ngm: bool + ngm: bool, } impl + 'static> TruncatedSvdResult { /// Returns singular values ordered by magnitude with indices. fn singular_values_with_indices(&self) -> (Array1, Vec) { // numerate and square root eigenvalues - let mut a = self.eigvals.iter() - .map(|x| x.sqrt()) - .enumerate() - .collect::>(); + let mut a = self.eigvals.iter().map(|x| x.sqrt()).enumerate().collect::>(); // sort by magnitude - a.sort_by(|(_,x), (_, y)| x.partial_cmp(&y).unwrap().reverse()); - + a.sort_by(|(_, x), (_, y)| x.partial_cmp(&y).unwrap().reverse()); + // filter low singular values away - let (values, indices): (Vec, Vec) = a.into_iter() - .filter(|(_,x)| *x > NumCast::from(1e-5).unwrap()) - .map(|(a,b)| (b,a)) + let (values, indices): (Vec, Vec) = a + .into_iter() + .filter(|(_, x)| *x > NumCast::from(1e-5).unwrap()) + .map(|(a, b)| (b, a)) .unzip(); (Array1::from(values), indices) @@ -58,19 +55,23 @@ impl + 'static> TruncatedSvdResult { let (u, v) = if self.ngm { let vlarge = self.eigvecs.select(Axis(1), &indices); let mut ularge = self.problem.dot(&vlarge); - - ularge.gencolumns_mut().into_iter() - .zip(values.iter()) - .for_each(|(mut a,b)| a.mapv_inplace(|x| x / *b)); + + ularge + .gencolumns_mut() + .into_iter() + .zip(values.iter()) + .for_each(|(mut a, b)| a.mapv_inplace(|x| x / *b)); (ularge, vlarge) } else { let ularge = self.eigvecs.select(Axis(1), &indices); let mut vlarge = self.problem.t().dot(&ularge); - vlarge.gencolumns_mut().into_iter() - .zip(values.iter()) - .for_each(|(mut a,b)| a.mapv_inplace(|x| x / *b)); + vlarge + .gencolumns_mut() + .into_iter() + .zip(values.iter()) + .for_each(|(mut a, b)| a.mapv_inplace(|x| x / *b)); (ularge, vlarge) }; @@ -89,7 +90,7 @@ pub struct TruncatedSvd { order: Order, problem: Array2, precision: A::Real, - maxiter: usize + maxiter: usize, } impl TruncatedSvd { @@ -97,8 +98,8 @@ impl TruncatedSvd { TruncatedSvd { precision: NumCast::from(1e-5).unwrap(), maxiter: problem.len_of(Axis(0)) * 2, - order, - problem + order, + problem, } } @@ -112,61 +113,72 @@ impl TruncatedSvd { self.maxiter = maxiter; self - } // calculate the eigenvalue decomposition pub fn decompose(self, num: usize) -> Result> { - let (n,m) = (self.problem.nrows(), self.problem.ncols()); + 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()); + 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 + // square precision because the SVD squares the eigenvalue as well let precision = self.precision * self.precision; // use problem definition with less operations required let res = if n > m { - lobpcg(|y| self.problem.t().dot(&self.problem.dot(&y)), x, None, None, precision, self.maxiter, self.order.clone()) + lobpcg( + |y| self.problem.t().dot(&self.problem.dot(&y)), + x, + None, + None, + precision, + self.maxiter, + self.order.clone(), + ) } else { - lobpcg(|y| self.problem.dot(&self.problem.t().dot(&y)), x, None, None, precision, self.maxiter, self.order.clone()) + lobpcg( + |y| self.problem.dot(&self.problem.t().dot(&y)), + x, + None, + None, + precision, + self.maxiter, + self.order.clone(), + ) }; // convert into TruncatedSvdResult match res { - EigResult::Ok(vals, vecs, _) | EigResult::Err(vals, vecs, _, _) => { - Ok(TruncatedSvdResult { - problem: self.problem.clone(), - eigvals: vals, - eigvecs: vecs, - ngm: n > m - }) - }, - EigResult::NoResult(err) => Err(err) + EigResult::Ok(vals, vecs, _) | EigResult::Err(vals, vecs, _, _) => Ok(TruncatedSvdResult { + problem: self.problem.clone(), + eigvals: vals, + eigvecs: vecs, + ngm: n > m, + }), + EigResult::NoResult(err) => Err(err), } } } #[cfg(test)] mod tests { - use crate::close_l2; - use super::TruncatedSvd; use super::Order; + use super::TruncatedSvd; + use crate::close_l2; use ndarray::{arr1, arr2, Array2}; use ndarray_rand::rand_distr::Uniform; use ndarray_rand::RandomExt; #[test] fn test_truncated_svd() { - let a = arr2(&[[3., 2., 2.], - [2., 3., -2.]]); + let a = arr2(&[[3., 2., 2.], [2., 3., -2.]]); let res = TruncatedSvd::new(a, Order::Largest) .precision(1e-5) .maxiter(10) .decompose(2) .unwrap(); - + let (_, sigma, _) = res.values_vecs(); close_l2(&sigma, &arr1(&[5.0, 3.0]), 1e-5); From f2909f0f023d5301a495f45302051484ce159eca Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Tue, 17 Mar 2020 14:51:07 +0100 Subject: [PATCH 16/30] Add eigenvalue decomposition example --- examples/truncated_eig.rs | 22 ++++++++++++++++++++++ examples/truncated_svd.rs | 22 ++++++++++++++++++++++ src/lobpcg/eig.rs | 10 +++++++++- src/lobpcg/lobpcg.rs | 7 ------- 4 files changed, 53 insertions(+), 8 deletions(-) create mode 100644 examples/truncated_eig.rs create mode 100644 examples/truncated_svd.rs diff --git a/examples/truncated_eig.rs b/examples/truncated_eig.rs new file mode 100644 index 00000000..765c62c7 --- /dev/null +++ b/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]); + } +} diff --git a/examples/truncated_svd.rs b/examples/truncated_svd.rs new file mode 100644 index 00000000..8cc164eb --- /dev/null +++ b/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); +} diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs index 8b3701eb..37d5a05e 100644 --- a/src/lobpcg/eig.rs +++ b/src/lobpcg/eig.rs @@ -83,6 +83,7 @@ impl IntoIterator for Truncat fn into_iter(self) -> TruncatedEigIterator { TruncatedEigIterator { step_size: 1, + remaining: self.problem.len_of(Axis(0)), eig: self, } } @@ -94,6 +95,7 @@ impl IntoIterator for Truncat /// eigenvalue/vector pair. Useful for generating pairs until a certain condition is met. pub struct TruncatedEigIterator { step_size: usize, + remaining: usize, eig: TruncatedEig, } @@ -101,7 +103,12 @@ impl Iterator for TruncatedEi type Item = (Array1, Array2); fn next(&mut self) -> Option { - 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 Iterator for TruncatedEi }; self.eig.constraints = Some(new_constraints); + self.remaining -= step_size; Some((vals, vecs)) } diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index bbf7d713..3d9122c6 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -87,17 +87,13 @@ fn apply_constraints( y: ArrayView2, ) { 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() @@ -105,9 +101,6 @@ fn apply_constraints( 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)); } From 68d8ec7c67d49ed9a264fd3417ba30891f7b1d0c Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Thu, 19 Mar 2020 10:08:27 +0100 Subject: [PATCH 17/30] Truncate svd with float eps; Take preconditioner as a closure --- examples/truncated_eig.rs | 1 + src/lobpcg/eig.rs | 30 +++++++++++++------ src/lobpcg/lobpcg.rs | 14 ++++----- src/lobpcg/svd.rs | 62 +++++++++++++++++++++++++++------------ 4 files changed, 73 insertions(+), 34 deletions(-) diff --git a/examples/truncated_eig.rs b/examples/truncated_eig.rs index 765c62c7..58172ec6 100644 --- 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 index 37d5a05e..edd53147 100644 --- a/src/lobpcg/eig.rs +++ b/src/lobpcg/eig.rs @@ -64,15 +64,27 @@ impl TruncatedEig { 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 index 3d9122c6..50ed232e 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -147,10 +147,10 @@ fn orthonormalize(v: Array2) -> Result<(Array2, 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) -> Array2>( +pub fn lobpcg) -> Array2, G: Fn(ArrayViewMut2)>( a: F, mut x: Array2, - m: Option>, + m: G, y: Option>, tol: A::Real, maxiter: usize, @@ -246,9 +246,9 @@ pub fn lobpcg) -> // 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 = 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 = Array2::random((10, 1), Uniform::new(0.0, 1.0)); let y: Array2 = 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 index 6d1a758f..717329f8 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -1,3 +1,6 @@ +///! 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}; @@ -5,14 +8,12 @@ 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 { eigvals: Array1, @@ -21,26 +22,31 @@ pub struct TruncatedSvdResult { ngm: bool, } -impl + 'static> TruncatedSvdResult { +impl + 'static + MagnitudeCorrection> TruncatedSvdResult { /// Returns singular values ordered by magnitude with indices. fn singular_values_with_indices(&self) -> (Array1, Vec) { - // numerate and square root eigenvalues - let mut a = self.eigvals.iter().map(|x| x.sqrt()).enumerate().collect::>(); + // numerate eigenvalues + let mut a = self.eigvals.iter().enumerate().collect::>(); // 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, Vec) = 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 { let (values, _) = self.singular_values_with_indices(); @@ -82,10 +88,8 @@ impl + 'static> TruncatedSvdResult { /// 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 { order: Order, problem: Array2, @@ -117,9 +121,15 @@ impl TruncatedSvd { // calculate the eigenvalue decomposition pub fn decompose(self, num: usize) -> Result> { + 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 TruncatedSvd { lobpcg( |y| self.problem.t().dot(&self.problem.dot(&y)), x, - None, + |_| {}, None, precision, self.maxiter, @@ -139,7 +149,7 @@ impl TruncatedSvd { lobpcg( |y| self.problem.dot(&self.problem.t().dot(&y)), x, - None, + |_| {}, None, precision, self.maxiter, @@ -160,6 +170,22 @@ impl TruncatedSvd { } } +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); } From a7ce0e323c9121bc20fabefdc020f6150423b364 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Thu, 19 Mar 2020 11:06:12 +0100 Subject: [PATCH 18/30] Use select for `ndarray_mask` --- src/lobpcg/lobpcg.rs | 55 +++++++++++--------------------------------- 1 file changed, 14 insertions(+), 41 deletions(-) diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 50ed232e..546a11a7 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -58,24 +58,13 @@ fn sorted_eig( /// Masks a matrix with the given `matrix` fn ndarray_mask(matrix: ArrayView2, mask: &[bool]) -> Array2 { - let (rows, cols) = (matrix.nrows(), matrix.ncols()); + assert_eq!(mask.len(), matrix.ncols()); - assert_eq!(mask.len(), cols); + let indices = (0..mask.len()).zip(mask.into_iter()) + .filter(|(_,b)| **b).map(|(a,_)| a) + .collect::>(); - let n_positive = mask.iter().filter(|x| **x).count(); - - let matrix = matrix - .gencolumns() - .into_iter() - .zip(mask.iter()) - .filter(|(_, x)| **x) - .map(|(x, _)| x.to_vec()) - .flatten() - .collect::>(); - - Array2::from_shape_vec((n_positive, rows), matrix) - .unwrap() - .reversed_axes() + matrix.select(Axis(1), &indices) } /// Applies constraints ensuring that a matrix is orthogonal to it @@ -193,19 +182,16 @@ pub fn lobpcg) -> let ax = a(x.view()); let xax = x.t().dot(&ax); - // perform eigenvalue decomposition on XAX + // perform eigenvalue decomposition of XAX as initialization let (mut lambda, eig_block) = match sorted_eig(xax.view(), None, size_x, &order) { Ok(x) => x, Err(err) => return EigResult::NoResult(err), }; - //dbg!(&lambda, &eig_block); - // initiate X and AX with eigenvector let mut x = x.dot(&eig_block); let mut ax = ax.dot(&eig_block); - //dbg!(&X, &AX); let mut activemask = vec![true; size_x]; let mut residual_norms = Vec::new(); let mut results = vec![(lambda.clone(), x.clone())]; @@ -219,10 +205,11 @@ pub fn lobpcg) -> let final_norm = loop { // calculate residual - let lambda_tmp = lambda.clone().insert_axis(Axis(0)); - let tmp = &x * &lambda_tmp; + let lambda_diag = Array2::from_diag(&lambda); + let lambda_x = x.dot(&lambda_diag); - let r = &ax - &tmp; + // calculate residual AX - lambdaX + let r = &ax - &lambda_x; // calculate L2 norm of error for every eigenvalue let tmp = r.gencolumns().into_iter().map(|x| x.norm()).collect::>(); @@ -248,7 +235,7 @@ pub fn lobpcg) -> let mut active_block_r = ndarray_mask(r.view(), &activemask); // apply preconditioner m(active_block_r.view_mut()); - + // apply constraints if let (Some(ref y), Some(ref fact_yy)) = (&y, &fact_yy) { apply_constraints(active_block_r.view_mut(), fact_yy, y.view()); } @@ -271,7 +258,7 @@ pub fn lobpcg) -> let active_ap = ndarray_mask(ap.view(), &activemask); let (active_p, p_r) = orthonormalize(active_p).unwrap(); - //dbg!(&active_P, &P_R); + let active_ap = match p_r.solve_triangular(UPLO::Lower, Diag::NonUnit, &active_ap.reversed_axes()) { Ok(x) => x, Err(err) => break Err(err), @@ -279,9 +266,6 @@ pub fn lobpcg) -> let active_ap = active_ap.reversed_axes(); - //dbg!(&active_AP); - //dbg!(&R); - let xap = x.t().dot(&active_ap); let wap = r.t().dot(&active_ap); let pap = active_p.t().dot(&active_ap); @@ -291,7 +275,7 @@ pub fn lobpcg) -> ( stack![ Axis(0), - stack![Axis(1), Array2::from_diag(&lambda), xaw, xap], + stack![Axis(1), lambda_diag, xaw, xap], stack![Axis(1), xaw.t(), waw, wap], stack![Axis(1), xap.t(), wap.t(), pap] ], @@ -308,7 +292,7 @@ pub fn lobpcg) -> ( stack![ Axis(0), - stack![Axis(1), Array2::from_diag(&lambda), xaw], + stack![Axis(1), lambda_diag, xaw], stack![Axis(1), xaw.t(), waw] ], stack![Axis(0), stack![Axis(1), ident0, xw], stack![Axis(1), xw.t(), ident]], @@ -323,25 +307,15 @@ pub fn lobpcg) -> }; lambda = new_lambda; - //dbg!(&lambda, &eig_vecs); let (pp, app, eig_x) = if let (Some(_), (Some(ref active_p), Some(ref active_ap))) = (ap, (active_p, active_ap)) { let eig_x = eig_vecs.slice(s![..size_x, ..]); let eig_r = eig_vecs.slice(s![size_x..size_x + current_block_size, ..]); let eig_p = eig_vecs.slice(s![size_x + current_block_size.., ..]); - //dbg!(&eig_X); - //dbg!(&eig_R); - //dbg!(&eig_P); - - //dbg!(&R, &AR, &active_P, &active_AP); - let pp = r.dot(&eig_r) + active_p.dot(&eig_p); let app = ar.dot(&eig_r) + active_ap.dot(&eig_p); - //dbg!(&pp); - //dbg!(&app); - (pp, app, eig_x) } else { let eig_x = eig_vecs.slice(s![..size_x, ..]); @@ -363,7 +337,6 @@ pub fn lobpcg) -> iter -= 1; }; - //dbg!(&residual_norms); let best_idx = residual_norms.iter().enumerate().min_by( |&(_, item1): &(usize, &Vec), &(_, item2): &(usize, &Vec)| { let norm1: A::Real = item1.iter().map(|x| (*x) * (*x)).sum(); From f25a177d60ad871a8748a7634e361465f32a11f0 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Thu, 19 Mar 2020 11:27:53 +0100 Subject: [PATCH 19/30] Only save best result for LOBPCG --- src/lobpcg/lobpcg.rs | 36 ++++++++++++------------------------ 1 file changed, 12 insertions(+), 24 deletions(-) diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 546a11a7..7ebae90c 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -194,7 +194,7 @@ pub fn lobpcg) -> let mut activemask = vec![true; size_x]; let mut residual_norms = Vec::new(); - let mut results = vec![(lambda.clone(), x.clone())]; + let mut best_result = None; let mut previous_block_size = size_x; @@ -215,6 +215,12 @@ pub fn lobpcg) -> let tmp = r.gencolumns().into_iter().map(|x| x.norm()).collect::>(); residual_norms.push(tmp.clone()); + // compare best result and update if we improved + let sum_rnorm: A::Real = tmp.iter().cloned().sum(); + if best_result.as_ref().map(|x: &(_,_,Vec)| x.2.iter().cloned().sum::() > sum_rnorm).unwrap_or(true) { + best_result = Some((lambda.clone(), x.clone(), tmp.clone())); + } + // disable eigenvalues which are below the tolerance threshold activemask = tmp.iter().zip(activemask.iter()).map(|(x, a)| *x > tol && *a).collect(); @@ -330,35 +336,17 @@ pub fn lobpcg) -> x = x.dot(&eig_x) + &pp; ax = ax.dot(&eig_x) + &app; - results.push((lambda.clone(), x.clone())); - ap = Some((pp, app)); iter -= 1; }; - let best_idx = residual_norms.iter().enumerate().min_by( - |&(_, item1): &(usize, &Vec), &(_, item2): &(usize, &Vec)| { - let norm1: A::Real = item1.iter().map(|x| (*x) * (*x)).sum(); - let norm2: A::Real = item2.iter().map(|x| (*x) * (*x)).sum(); - norm1.partial_cmp(&norm2).unwrap() - }, - ); + let (vals, vecs, rnorm) = best_result.unwrap(); + let rnorm = rnorm.into_iter().map(|x| Scalar::from_real(x)).collect(); - match best_idx { - Some((idx, norms)) => { - let (vals, vecs) = results[idx].clone(); - let norms = norms.iter().map(|x| Scalar::from_real(*x)).collect(); - - match final_norm { - Ok(_) => EigResult::Ok(vals, vecs, norms), - Err(err) => EigResult::Err(vals, vecs, norms, err), - } - } - None => match final_norm { - Ok(_) => panic!("Not error available!"), - Err(err) => EigResult::NoResult(err), - }, + match final_norm { + Ok(_) => EigResult::Ok(vals, vecs, rnorm), + Err(err) => EigResult::Err(vals, vecs, rnorm, err) } } From cea3c0eda56e59ca541bfa38495dea77adf7fa47 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Thu, 19 Mar 2020 11:46:45 +0100 Subject: [PATCH 20/30] Add benchmarks for truncated eigenvalue decomposition --- Cargo.toml | 5 +++++ benches/truncated_eig.rs | 41 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+) create mode 100644 benches/truncated_eig.rs diff --git a/Cargo.toml b/Cargo.toml index eeb2ddf9..5bd0aac7 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -53,6 +53,11 @@ optional = true paste = "0.1" criterion = "0.3" +[[bench]] +name = "truncated_eig" +harness = false + [[bench]] name = "eigh" harness = false + diff --git a/benches/truncated_eig.rs b/benches/truncated_eig.rs new file mode 100644 index 00000000..b24aac7c --- /dev/null +++ b/benches/truncated_eig.rs @@ -0,0 +1,41 @@ +#[macro_use] +extern crate criterion; + +use criterion::Criterion; +use ndarray::*; +use ndarray_linalg::*; + +macro_rules! impl_teig { + ($n:expr) => { + paste::item! { + fn [](c: &mut Criterion) { + c.bench_function(&format!("truncated_eig{}", $n), |b| { + let a: Array2 = random(($n, $n)); + let a = a.t().dot(&a); + + b.iter(move || { + let _result = TruncatedEig::new(a.clone(), TruncatedOrder::Largest).decompose(1); + }) + }); + c.bench_function(&format!("truncated_eig{}_t", $n), |b| { + let a: Array2 = random(($n, $n).f()); + let a = a.t().dot(&a); + + b.iter(|| { + let _result = TruncatedEig::new(a.clone(), TruncatedOrder::Largest).decompose(1); + }) + }); + } + } + }; +} + +impl_teig!(4); +impl_teig!(8); +impl_teig!(16); +impl_teig!(32); +impl_teig!(64); +impl_teig!(128); + +criterion_group!(teig, teig4, teig8, teig16, teig32, teig64, teig128); +criterion_main!(teig); From 6c94817aca0963420a41db92a3b3c3ed96669a8e Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Thu, 19 Mar 2020 13:07:46 +0100 Subject: [PATCH 21/30] Add restarting and improve performance with explicit gram flag --- src/lobpcg/eig.rs | 7 ++- src/lobpcg/lobpcg.rs | 129 +++++++++++++++++++++++++++++-------------- src/lobpcg/svd.rs | 3 +- 3 files changed, 94 insertions(+), 45 deletions(-) diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs index edd53147..187f7055 100644 --- a/src/lobpcg/eig.rs +++ b/src/lobpcg/eig.rs @@ -4,6 +4,7 @@ use crate::{Lapack, Scalar}; /// use ndarray::prelude::*; use ndarray::stack; +use ndarray::ScalarOperand; use ndarray_rand::rand_distr::Uniform; use ndarray_rand::RandomExt; use num_traits::{Float, NumCast}; @@ -23,7 +24,7 @@ pub struct TruncatedEig { maxiter: usize, } -impl TruncatedEig { +impl TruncatedEig { pub fn new(problem: Array2, order: Order) -> TruncatedEig { TruncatedEig { precision: NumCast::from(1e-5).unwrap(), @@ -88,7 +89,7 @@ impl TruncatedEig { } } -impl IntoIterator for TruncatedEig { +impl IntoIterator for TruncatedEig { type Item = (Array1, Array2); type IntoIter = TruncatedEigIterator; @@ -111,7 +112,7 @@ pub struct TruncatedEigIterator { eig: TruncatedEig, } -impl Iterator for TruncatedEigIterator { +impl Iterator for TruncatedEigIterator { type Item = (Array1, Array2); fn next(&mut self) -> Option { diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 7ebae90c..638c7cf1 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -7,7 +7,8 @@ use crate::{cholesky::*, close_l2, eigh::*, norm::*, triangular::*}; use crate::{Lapack, Scalar}; use ndarray::prelude::*; use ndarray::OwnedRepr; -use num_traits::NumCast; +use ndarray::ScalarOperand; +use num_traits::{NumCast, Float}; /// Find largest or smallest eigenvalues #[derive(Debug, Clone)] @@ -136,7 +137,7 @@ fn orthonormalize(v: Array2) -> Result<(Array2, 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) -> Array2, G: Fn(ArrayViewMut2)>( +pub fn lobpcg) -> Array2, G: Fn(ArrayViewMut2)>( a: F, mut x: Array2, m: G, @@ -193,15 +194,17 @@ pub fn lobpcg) -> let mut ax = ax.dot(&eig_block); let mut activemask = vec![true; size_x]; - let mut residual_norms = Vec::new(); + let mut residual_norms_history = Vec::new(); let mut best_result = None; let mut previous_block_size = size_x; let mut ident: Array2 = Array2::eye(size_x); let ident0: Array2 = Array2::eye(size_x); + let two: A = NumCast::from(2.0).unwrap(); let mut ap: Option<(Array2, Array2)> = None; + let mut explicit_gram_flag = false; let final_norm = loop { // calculate residual @@ -212,17 +215,17 @@ pub fn lobpcg) -> let r = &ax - &lambda_x; // calculate L2 norm of error for every eigenvalue - let tmp = r.gencolumns().into_iter().map(|x| x.norm()).collect::>(); - residual_norms.push(tmp.clone()); + let residual_norms = r.gencolumns().into_iter().map(|x| x.norm()).collect::>(); + residual_norms_history.push(residual_norms.clone()); // compare best result and update if we improved - let sum_rnorm: A::Real = tmp.iter().cloned().sum(); + let sum_rnorm: A::Real = residual_norms.iter().cloned().sum(); if best_result.as_ref().map(|x: &(_,_,Vec)| x.2.iter().cloned().sum::() > sum_rnorm).unwrap_or(true) { - best_result = Some((lambda.clone(), x.clone(), tmp.clone())); + best_result = Some((lambda.clone(), x.clone(), residual_norms.clone())); } // disable eigenvalues which are below the tolerance threshold - activemask = tmp.iter().zip(activemask.iter()).map(|(x, a)| *x > tol && *a).collect(); + activemask = residual_norms.iter().zip(activemask.iter()).map(|(x, a)| *x > tol && *a).collect(); // resize identity block if necessary let current_block_size = activemask.iter().filter(|x| **x).count(); @@ -234,17 +237,19 @@ pub fn lobpcg) -> // if we are below the threshold for all eigenvalue or exceeded the number of iteration, // abort if current_block_size == 0 || iter == 0 { - break Ok(tmp); + break Ok(residual_norms); } // select active eigenvalues, apply pre-conditioner, orthogonalize to Y and orthonormalize let mut active_block_r = ndarray_mask(r.view(), &activemask); // apply preconditioner m(active_block_r.view_mut()); - // apply constraints + // apply constraints to the preconditioned residuals if let (Some(ref y), Some(ref fact_yy)) = (&y, &fact_yy) { apply_constraints(active_block_r.view_mut(), fact_yy, y.view()); } + // orthogonalize the preconditioned residual to x + active_block_r -= &x.dot(&x.t().dot(&active_block_r)); let (r, _) = match orthonormalize(active_block_r) { Ok(x) => x, @@ -253,57 +258,99 @@ pub fn lobpcg) -> let ar = a(r.view()); + let max_rnorm = if A::epsilon() > NumCast::from(1e-8).unwrap() { + NumCast::from(1.0).unwrap() + } else { + NumCast::from(1.0e-8).unwrap() + }; + + // if we once below the max_rnorm enable explicit gram flag + let max_norm = residual_norms.into_iter().fold(A::Real::neg_infinity(), A::Real::max); + + explicit_gram_flag = max_norm <= max_rnorm || explicit_gram_flag; + // perform the Rayleigh Ritz procedure - let xaw = x.t().dot(&ar); - let waw = r.t().dot(&ar); - let xw = x.t().dot(&r); + let xar = x.t().dot(&ar); + let mut rar = r.t().dot(&ar); - // compute symmetric gram matrices - let (gram_a, gram_b, active_p, active_ap) = if let Some((ref p, ref ap)) = ap { + let (xax, xx, rr, xr) = if explicit_gram_flag { + rar = (&rar + &rar.t()) / two; + let xax = x.t().dot(&ax); + + ( + (&xax + &xax.t()) / two, + x.t().dot(&x), + r.t().dot(&r), + x.t().dot(&r) + ) + } else { + ( + lambda_diag, + ident0.clone(), + ident.clone(), + Array2::zeros((size_x, current_block_size)) + ) + }; + + let p_ap: Option<(_,_)> = ap.as_ref().and_then(|(p, ap)| { let active_p = ndarray_mask(p.view(), &activemask); let active_ap = ndarray_mask(ap.view(), &activemask); - let (active_p, p_r) = orthonormalize(active_p).unwrap(); + if let Ok((active_p, p_r)) = orthonormalize(active_p) { + if let Ok(active_ap) = p_r.solve_triangular(UPLO::Lower, Diag::NonUnit, &active_ap.reversed_axes()) { + let active_ap = active_ap.reversed_axes(); + + Some((active_p, active_ap)) + } else { + None + } + } else { + None + } + }); - let active_ap = match p_r.solve_triangular(UPLO::Lower, Diag::NonUnit, &active_ap.reversed_axes()) { - Ok(x) => x, - Err(err) => break Err(err), + // compute symmetric gram matrices + let (gram_a, gram_b) = if let Some((active_p, active_ap)) = &p_ap { + let xap = x.t().dot(active_ap); + let rap = r.t().dot(active_ap); + let pap = active_p.t().dot(active_ap); + let xp = x.t().dot(active_p); + let rp = r.t().dot(active_p); + let (pap, pp) = if explicit_gram_flag { + ( + (&pap + &pap.t()) / two, + active_p.t().dot(active_p) + ) + } else { + (pap, ident.clone()) }; - let active_ap = active_ap.reversed_axes(); - - let xap = x.t().dot(&active_ap); - let wap = r.t().dot(&active_ap); - let pap = active_p.t().dot(&active_ap); - let xp = x.t().dot(&active_p); - let wp = r.t().dot(&active_p); - ( stack![ Axis(0), - stack![Axis(1), lambda_diag, xaw, xap], - stack![Axis(1), xaw.t(), waw, wap], - stack![Axis(1), xap.t(), wap.t(), pap] + stack![Axis(1), xax, xar, xap], + stack![Axis(1), xar.t(), rar, rap], + stack![Axis(1), xap.t(), rap.t(), pap] ], stack![ Axis(0), - stack![Axis(1), ident0, xw, xp], - stack![Axis(1), xw.t(), ident, wp], - stack![Axis(1), xp.t(), wp.t(), ident] + stack![Axis(1), xx, xr, xp], + stack![Axis(1), xr.t(), rr, rp], + stack![Axis(1), xp.t(), rp.t(), pp] ], - Some(active_p), - Some(active_ap), ) } else { ( stack![ Axis(0), - stack![Axis(1), lambda_diag, xaw], - stack![Axis(1), xaw.t(), waw] + stack![Axis(1), xax, xar], + stack![Axis(1), xar.t(), rar] + ], + stack![ + Axis(0), + stack![Axis(1), xx, xr], + stack![Axis(1), xr.t(), rr] ], - stack![Axis(0), stack![Axis(1), ident0, xw], stack![Axis(1), xw.t(), ident]], - None, - None, ) }; @@ -313,7 +360,7 @@ pub fn lobpcg) -> }; lambda = new_lambda; - let (pp, app, eig_x) = if let (Some(_), (Some(ref active_p), Some(ref active_ap))) = (ap, (active_p, active_ap)) + let (pp, app, eig_x) = if let (Some(_), Some((active_p, active_ap))) = (ap, p_ap) { let eig_x = eig_vecs.slice(s![..size_x, ..]); let eig_r = eig_vecs.slice(s![size_x..size_x + current_block_size, ..]); diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index 717329f8..feb55c12 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -5,6 +5,7 @@ use super::lobpcg::{lobpcg, EigResult, Order}; use crate::error::Result; use crate::{Lapack, Scalar}; use ndarray::prelude::*; +use ndarray::ScalarOperand; use ndarray_rand::rand_distr::Uniform; use ndarray_rand::RandomExt; use num_traits::{Float, NumCast}; @@ -97,7 +98,7 @@ pub struct TruncatedSvd { maxiter: usize, } -impl TruncatedSvd { +impl TruncatedSvd { pub fn new(problem: Array2, order: Order) -> TruncatedSvd { TruncatedSvd { precision: NumCast::from(1e-5).unwrap(), From 233fc8a0b8736723b7264b66ac597ff00139007b Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Thu, 19 Mar 2020 13:15:42 +0100 Subject: [PATCH 22/30] Restart the eigenvalue decomposition as well --- src/lobpcg/lobpcg.rs | 39 ++++++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 638c7cf1..26e93d5c 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -264,7 +264,7 @@ pub fn lobpcg = ap.as_ref().and_then(|(p, ap)| { - let active_p = ndarray_mask(p.view(), &activemask); - let active_ap = ndarray_mask(ap.view(), &activemask); + let p_ap = ap.as_ref() + .and_then(|(p, ap)| { + let active_p = ndarray_mask(p.view(), &activemask); + let active_ap = ndarray_mask(ap.view(), &activemask); - if let Ok((active_p, p_r)) = orthonormalize(active_p) { - if let Ok(active_ap) = p_r.solve_triangular(UPLO::Lower, Diag::NonUnit, &active_ap.reversed_axes()) { - let active_ap = active_ap.reversed_axes(); - - Some((active_p, active_ap)) - } else { - None - } - } else { - None - } - }); + orthonormalize(active_p).map(|x| (active_ap, x)).ok() + }) + .and_then(|(active_ap, (active_p, p_r))| { + let active_ap = active_ap.reversed_axes(); + p_r.solve_triangular(UPLO::Lower, Diag::NonUnit, &active_ap) + .map(|active_ap| (active_p, active_ap.reversed_axes())) + .ok() + }); // compute symmetric gram matrices let (gram_a, gram_b) = if let Some((active_p, active_ap)) = &p_ap { @@ -356,7 +353,15 @@ pub fn lobpcg x, - Err(err) => break Err(err), + Err(err) => { + // restart if the eigproblem decomposition failed + if ap.is_some() { + ap = None; + continue; + } else { + break Err(err); + } + } }; lambda = new_lambda; From 93f8b44af64ad592694e8c01c61d3346f5778ff0 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Thu, 19 Mar 2020 13:36:16 +0100 Subject: [PATCH 23/30] Remove unnecessary match blocks --- src/lobpcg/lobpcg.rs | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 26e93d5c..8d65764c 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -204,7 +204,7 @@ pub fn lobpcg, Array2)> = None; - let mut explicit_gram_flag = false; + let mut explicit_gram_flag = true; let final_norm = loop { // calculate residual @@ -258,7 +258,8 @@ pub fn lobpcg NumCast::from(1e-8).unwrap() { + // check whether `A` is of type `f32` or `f64` + let max_rnorm_float = if A::epsilon() > NumCast::from(1e-8).unwrap() { NumCast::from(1.0).unwrap() } else { NumCast::from(1.0e-8).unwrap() @@ -266,8 +267,7 @@ pub fn lobpcg EigResult::Ok(vals, vecs, rnorm), Err(err) => EigResult::Err(vals, vecs, rnorm, err) @@ -466,7 +468,7 @@ mod tests { let n = a.len_of(Axis(0)); let x: Array2 = Array2::random((n, num), Uniform::new(0.0, 1.0)); - let result = lobpcg(|y| a.dot(&y), x, |_| {}, None, 1e-10, 2 * n, order); + let result = lobpcg(|y| a.dot(&y), x, |_| {}, None, 1e-5, n, order); match result { EigResult::Ok(vals, _, r_norms) | EigResult::Err(vals, _, r_norms, _) => { // check convergence From e9f9f6a5bdb7185dd0bf9b341aa398f2f4dd66fb Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Thu, 19 Mar 2020 13:52:21 +0100 Subject: [PATCH 24/30] Increase precision in orthonormalization test --- src/lobpcg/lobpcg.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 8d65764c..a695b561 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -396,7 +396,7 @@ pub fn lobpcg EigResult::Ok(vals, vecs, rnorm), @@ -445,7 +445,7 @@ mod tests { /// Test orthonormalization of a random matrix #[test] fn test_orthonormalize() { - let matrix = Array2::random((10, 10), Uniform::new(-10., 10.)); + let matrix: Array2 = Array2::random((10, 10), Uniform::new(-10., 10.)); let (n, l) = orthonormalize(matrix.clone()).unwrap(); @@ -455,7 +455,7 @@ mod tests { // compare returned factorization with QR decomposition let (_, r) = matrix.qr().unwrap(); - close_l2(&r.mapv(|x: f32| x.abs()), &l.t().mapv(|x| x.abs()), 1e-2); + close_l2(&r.mapv(|x| x.abs()), &l.t().mapv(|x| x.abs()), 1e-2); } fn assert_symmetric(a: &Array2) { From 564bb77efaf84b7697270f0381115f923bb86526 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Sun, 22 Mar 2020 11:48:42 +0100 Subject: [PATCH 25/30] Improve comments on LOBPCG --- src/lobpcg/lobpcg.rs | 135 ++++++++++++++++++++++++------------------- 1 file changed, 75 insertions(+), 60 deletions(-) diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index a695b561..a25afe90 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -73,7 +73,7 @@ fn ndarray_mask(matrix: ArrayView2, mask: &[bool]) -> Array2 { /// This functions takes a matrix `v` and constraint matrix `y` and orthogonalize the `v` to `y`. fn apply_constraints( mut v: ArrayViewMut, - fact_yy: &CholeskyFactorized>, + cholesky_yy: &CholeskyFactorized>, y: ArrayView2, ) { let gram_yv = y.t().dot(&v); @@ -82,7 +82,7 @@ fn apply_constraints( .gencolumns() .into_iter() .map(|x| { - let res = fact_yy.solvec(&x).unwrap(); + let res = cholesky_yy.solvec(&x).unwrap(); res.to_vec() }) @@ -120,23 +120,22 @@ fn orthonormalize(v: Array2) -> Result<(Array2, Array2 /// /// # Arguments /// * `a` - An operator defining the problem, usually a sparse (sometimes also dense) matrix -/// multiplication. Also called the "Stiffness matrix". -/// * `x` - Initial approximation to the k eigenvectors. If `a` has shape=(n,n), then `x` should +/// multiplication. Also called the "stiffness matrix". +/// * `x` - Initial approximation of the k eigenvectors. If `a` has shape=(n,n), then `x` should /// have shape=(n,k). -/// * `m` - Preconditioner to `a`, by default the identity matrix. In the optimal case `m` -/// approximates the inverse of `a`. +/// * `m` - Preconditioner to `a`, by default the identity matrix. Should approximate the inverse +/// of `a`. /// * `y` - Constraints of (n,size_y), iterations are performed in the orthogonal complement of the /// column-space of `y`. It must be full rank. -/// * `tol` - The tolerance values defines at which point the solver stops the optimization. The l2-norm -/// of the residual is compared to this value and the eigenvalue approximation returned if below -/// the threshold. +/// * `tol` - The tolerance values defines at which point the solver stops the optimization. The approximation +/// of a eigenvalue stops when then l2-norm of the residual is below this threshold. /// * `maxiter` - The maximal number of iterations /// * `order` - Whether to solve for the largest or lowest eigenvalues /// /// The function returns an `EigResult` with the eigenvalue/eigenvector and achieved residual norm /// 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). +/// the precision of the matrix is too low (switch then from `f32` to `f64` for example). pub fn lobpcg) -> Array2, G: Fn(ArrayViewMut2)>( a: F, mut x: Array2, @@ -163,37 +162,37 @@ pub fn lobpcg { - let fact_yy = y.t().dot(y).factorizec(UPLO::Lower).unwrap(); + // calculate cholesky factorization of YY' and apply constraints to initial guess + let cholesky_yy = y.as_ref().map(|y| { + let cholesky_yy = y.t().dot(y).factorizec(UPLO::Lower).unwrap(); + apply_constraints(x.view_mut(), &cholesky_yy, y.view()); + cholesky_yy + }); - apply_constraints(x.view_mut(), &fact_yy, y.view()); - Some(fact_yy) - } - None => None, - }; - - // orthonormalize the initial guess and calculate matrices AX and XAX + // orthonormalize the initial guess let (x, _) = match orthonormalize(x) { Ok(x) => x, Err(err) => return EigResult::NoResult(err), }; + // calculate AX and XAX for Rayleigh quotient let ax = a(x.view()); let xax = x.t().dot(&ax); - // perform eigenvalue decomposition of XAX as initialization + // perform eigenvalue decomposition of XAX let (mut lambda, eig_block) = match sorted_eig(xax.view(), None, size_x, &order) { Ok(x) => x, Err(err) => return EigResult::NoResult(err), }; - // initiate X and AX with eigenvector + // initiate approximation of the eigenvector let mut x = x.dot(&eig_block); let mut ax = ax.dot(&eig_block); + // track residual below threshold let mut activemask = vec![true; size_x]; + + // track residuals and best result let mut residual_norms_history = Vec::new(); let mut best_result = None; @@ -203,7 +202,7 @@ pub fn lobpcg = Array2::eye(size_x); let two: A = NumCast::from(2.0).unwrap(); - let mut ap: Option<(Array2, Array2)> = None; + let mut previous_p_ap: Option<(Array2, Array2)> = None; let mut explicit_gram_flag = true; let final_norm = loop { @@ -245,8 +244,8 @@ pub fn lobpcg x, Err(err) => { // restart if the eigproblem decomposition failed - if ap.is_some() { - ap = None; + if previous_p_ap.is_some() { + previous_p_ap = None; continue; - } else { + } else { // or break if restart is not possible break Err(err); } } }; lambda = new_lambda; - let (pp, app, eig_x) = if let Some((active_p, active_ap)) = p_ap + // approximate eigenvector X and conjugate vectors P with solution of eigenproblem + let (p, ap, tau) = if let Some((active_p, active_ap)) = p_ap { - let eig_x = eig_vecs.slice(s![..size_x, ..]); - let eig_r = eig_vecs.slice(s![size_x..size_x + current_block_size, ..]); - let eig_p = eig_vecs.slice(s![size_x + current_block_size.., ..]); - - let pp = r.dot(&eig_r) + active_p.dot(&eig_p); - let app = ar.dot(&eig_r) + active_ap.dot(&eig_p); - - (pp, app, eig_x) + // tau are eigenvalues to basis of X + let tau = eig_vecs.slice(s![..size_x, ..]); + // alpha are eigenvalues to basis of R + let alpha = eig_vecs.slice(s![size_x..size_x + current_block_size, ..]); + // gamma are eigenvalues to basis of P + let gamma = eig_vecs.slice(s![size_x + current_block_size.., ..]); + + // update AP and P in span{R, P} as linear combination + let updated_p = r.dot(&alpha) + active_p.dot(&gamma); + let updated_ap = ar.dot(&alpha) + active_ap.dot(&gamma); + + (updated_p, updated_ap, tau) } else { - let eig_x = eig_vecs.slice(s![..size_x, ..]); - let eig_r = eig_vecs.slice(s![size_x.., ..]); + // tau are eigenvalues to basis of X + let tau = eig_vecs.slice(s![..size_x, ..]); + // alpha are eigenvalues to basis of R + let alpha = eig_vecs.slice(s![size_x.., ..]); - let pp = r.dot(&eig_r); - let app = ar.dot(&eig_r); + // update AP and P as linear combination of the residual matrix R + let updated_p = r.dot(&alpha); + let updated_ap = ar.dot(&alpha); - (pp, app, eig_x) + (updated_p, updated_ap, tau) }; - x = x.dot(&eig_x) + &pp; - ax = ax.dot(&eig_x) + &app; + // update approximation of X as linear combinations of span{X, P, R} + x = x.dot(&tau) + &p; + ax = ax.dot(&tau) + ≈ - ap = Some((pp, app)); + previous_p_ap = Some((p, ap)); iter -= 1; }; + // retrieve best result and convert norm into `A` let (vals, vecs, rnorm) = best_result.unwrap(); let rnorm = rnorm.into_iter().map(|x| Scalar::from_real(x)).collect(); - //dbg!(&residual_norms_history); + dbg!(&residual_norms_history); match final_norm { Ok(_) => EigResult::Ok(vals, vecs, rnorm), @@ -468,21 +483,22 @@ mod tests { let n = a.len_of(Axis(0)); let x: Array2 = Array2::random((n, num), Uniform::new(0.0, 1.0)); - let result = lobpcg(|y| a.dot(&y), x, |_| {}, None, 1e-5, n, order); + let result = lobpcg(|y| a.dot(&y), x, |_| {}, None, 1e-5, n * 2, order); + dbg!(&result); match result { EigResult::Ok(vals, _, r_norms) | EigResult::Err(vals, _, r_norms, _) => { // check convergence for (i, norm) in r_norms.into_iter().enumerate() { - if norm > 0.01 { + if norm > 1e-5 { println!("==== Assertion Failed ===="); - println!("The {} eigenvalue estimation did not converge!", i); + println!("The {}th eigenvalue estimation did not converge!", i); panic!("Too large deviation of residual norm: {} > 0.01", norm); } } // check correct order of eigenvalues if ground_truth_eigvals.len() == num { - close_l2(&Array1::from(ground_truth_eigvals.to_vec()), &vals, 5e-2) + close_l2(&Array1::from(ground_truth_eigvals.to_vec()), &vals, num as f64 * 5e-4) } } EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), @@ -519,30 +535,29 @@ mod tests { } #[test] - fn test_eigsolver_constrainted() { + fn test_eigsolver_constrained() { let diag = arr1(&[1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]); let a = Array2::from_diag(&diag); let x: Array2 = Array2::random((10, 1), Uniform::new(0.0, 1.0)); - let y: Array2 = arr2(&[[1.0, 0., 0., 0., 0., 0., 0., 0., 0., 0.]]).reversed_axes(); + let y: Array2 = arr2(&[[1.0, 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 1.0, 0., 0., 0., 0., 0., 0., 0., 0.]]).reversed_axes(); - let result = lobpcg(|y| a.dot(&y), x, |_| {}, Some(y), 1e-10, 100, Order::Smallest); - dbg!(&result); + let result = lobpcg(|y| a.dot(&y), x, |_| {}, Some(y), 1e-10, 50, Order::Smallest); match result { EigResult::Ok(vals, vecs, r_norms) | EigResult::Err(vals, vecs, r_norms, _) => { // check convergence for (i, norm) in r_norms.into_iter().enumerate() { if norm > 0.01 { println!("==== Assertion Failed ===="); - println!("The {} eigenvalue estimation did not converge!", i); + println!("The {}th eigenvalue estimation did not converge!", i); panic!("Too large deviation of residual norm: {} > 0.01", norm); } } - // should be the second eigenvalue - close_l2(&vals, &Array1::from(vec![2.0]), 1e-2); + // should be the third eigenvalue + close_l2(&vals, &Array1::from(vec![3.0]), 1e-10); close_l2( &vecs.column(0).mapv(|x| x.abs()), - &arr1(&[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), + &arr1(&[0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), 1e-5, ); } From a6ac4b3ab53264beadfe45f8994d76b87c3f26db Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Mon, 23 Mar 2020 17:20:43 +0100 Subject: [PATCH 26/30] Rename `EigResult` to `LobpcgResult` --- src/lobpcg/eig.rs | 8 ++-- src/lobpcg/lobpcg.rs | 101 ++++++++++++++++++++++--------------------- src/lobpcg/mod.rs | 2 +- src/lobpcg/svd.rs | 13 +++--- 4 files changed, 63 insertions(+), 61 deletions(-) diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs index 187f7055..e69734e7 100644 --- a/src/lobpcg/eig.rs +++ b/src/lobpcg/eig.rs @@ -1,4 +1,4 @@ -use super::lobpcg::{lobpcg, EigResult, Order}; +use super::lobpcg::{lobpcg, LobpcgResult, Order}; use crate::{Lapack, Scalar}; ///! Implements truncated eigenvalue decomposition /// @@ -61,7 +61,7 @@ impl Truncate } // calculate the eigenvalues decompose - pub fn decompose(&self, num: usize) -> EigResult { + pub fn decompose(&self, num: usize) -> LobpcgResult { let x = Array2::random((self.problem.len_of(Axis(0)), num), Uniform::new(0.0, 1.0)) .mapv(|x| NumCast::from(x).unwrap()); @@ -124,7 +124,7 @@ impl Iterator let res = self.eig.decompose(step_size); match res { - EigResult::Ok(vals, vecs, norms) | EigResult::Err(vals, vecs, norms, _) => { + LobpcgResult::Ok(vals, vecs, norms) | LobpcgResult::Err(vals, vecs, norms, _) => { // abort if any eigenproblem did not converge for r_norm in norms { if r_norm > NumCast::from(0.1).unwrap() { @@ -151,7 +151,7 @@ impl Iterator Some((vals, vecs)) } - EigResult::NoResult(_) => None, + LobpcgResult::NoResult(_) => None, } } } diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index a25afe90..40df72c5 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -8,7 +8,7 @@ use crate::{Lapack, Scalar}; use ndarray::prelude::*; use ndarray::OwnedRepr; use ndarray::ScalarOperand; -use num_traits::{NumCast, Float}; +use num_traits::{Float, NumCast}; /// Find largest or smallest eigenvalues #[derive(Debug, Clone)] @@ -20,12 +20,12 @@ pub enum Order { /// The result of the eigensolver /// /// In the best case the eigensolver has converged with a result better than the given threshold, -/// then a `EigResult::Ok` gives the eigenvalues, eigenvectors and norms. If an error ocurred -/// during the process, it is returned in `EigResult::Err`, but the best result is still returned, -/// as it could be usable. If there is no result at all, then `EigResult::NoResult` is returned. +/// then a `LobpcgResult::Ok` gives the eigenvalues, eigenvectors and norms. If an error ocurred +/// during the process, it is returned in `LobpcgResult::Err`, but the best result is still returned, +/// as it could be usable. If there is no result at all, then `LobpcgResult::NoResult` is returned. /// This happens if the algorithm fails in an early stage, for example if the matrix `A` is not SPD #[derive(Debug)] -pub enum EigResult { +pub enum LobpcgResult { Ok(Array1, Array2, Vec), Err(Array1, Array2, Vec, LinalgError), NoResult(LinalgError), @@ -61,8 +61,10 @@ fn sorted_eig( fn ndarray_mask(matrix: ArrayView2, mask: &[bool]) -> Array2 { assert_eq!(mask.len(), matrix.ncols()); - let indices = (0..mask.len()).zip(mask.into_iter()) - .filter(|(_,b)| **b).map(|(a,_)| a) + let indices = (0..mask.len()) + .zip(mask.into_iter()) + .filter(|(_, b)| **b) + .map(|(a, _)| a) .collect::>(); matrix.select(Axis(1), &indices) @@ -70,7 +72,7 @@ fn ndarray_mask(matrix: ArrayView2, mask: &[bool]) -> Array2 { /// Applies constraints ensuring that a matrix is orthogonal to it /// -/// This functions takes a matrix `v` and constraint matrix `y` and orthogonalize the `v` to `y`. +/// This functions takes a matrix `v` and constraint-matrix `y` and orthogonalize `v` to `y`. fn apply_constraints( mut v: ArrayViewMut, cholesky_yy: &CholeskyFactorized>, @@ -132,11 +134,15 @@ fn orthonormalize(v: Array2) -> Result<(Array2, Array2 /// * `maxiter` - The maximal number of iterations /// * `order` - Whether to solve for the largest or lowest eigenvalues /// -/// The function returns an `EigResult` with the eigenvalue/eigenvector and achieved residual norm +/// The function returns an `LobpcgResult` with the eigenvalue/eigenvector and achieved residual norm /// 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 +/// special variant `LobpcgResult::NotConverged` additionally carries the error. This can happen when /// the precision of the matrix is too low (switch then from `f32` to `f64` for example). -pub fn lobpcg) -> Array2, G: Fn(ArrayViewMut2)>( +pub fn lobpcg< + A: Float + Scalar + Lapack + ScalarOperand + PartialOrd + Default, + F: Fn(ArrayView2) -> Array2, + G: Fn(ArrayViewMut2), +>( a: F, mut x: Array2, m: G, @@ -144,7 +150,7 @@ pub fn lobpcg EigResult { +) -> LobpcgResult { // the initital approximation should be maximal square // n is the dimensionality of the problem let (n, size_x) = (x.nrows(), x.ncols()); @@ -172,7 +178,7 @@ pub fn lobpcg x, - Err(err) => return EigResult::NoResult(err), + Err(err) => return LobpcgResult::NoResult(err), }; // calculate AX and XAX for Rayleigh quotient @@ -182,7 +188,7 @@ pub fn lobpcg x, - Err(err) => return EigResult::NoResult(err), + Err(err) => return LobpcgResult::NoResult(err), }; // initiate approximation of the eigenvector @@ -219,12 +225,20 @@ pub fn lobpcg)| x.2.iter().cloned().sum::() > sum_rnorm).unwrap_or(true) { + if best_result + .as_ref() + .map(|x: &(_, _, Vec)| x.2.iter().cloned().sum::() > sum_rnorm) + .unwrap_or(true) + { best_result = Some((lambda.clone(), x.clone(), residual_norms.clone())); } // disable eigenvalues which are below the tolerance threshold - activemask = residual_norms.iter().zip(activemask.iter()).map(|(x, a)| *x > tol && *a).collect(); + activemask = residual_norms + .iter() + .zip(activemask.iter()) + .map(|(x, a)| *x > tol && *a) + .collect(); // resize identity block if necessary let current_block_size = activemask.iter().filter(|x| **x).count(); @@ -279,23 +293,19 @@ pub fn lobpcg EigResult::Ok(vals, vecs, rnorm), - Err(err) => EigResult::Err(vals, vecs, rnorm, err) + Ok(_) => LobpcgResult::Ok(vals, vecs, rnorm), + Err(err) => LobpcgResult::Err(vals, vecs, rnorm, err), } } @@ -425,7 +424,7 @@ mod tests { use super::ndarray_mask; use super::orthonormalize; use super::sorted_eig; - use super::EigResult; + use super::LobpcgResult; use super::Order; use crate::close_l2; use crate::qr::*; @@ -486,7 +485,7 @@ mod tests { let result = lobpcg(|y| a.dot(&y), x, |_| {}, None, 1e-5, n * 2, order); dbg!(&result); match result { - EigResult::Ok(vals, _, r_norms) | EigResult::Err(vals, _, r_norms, _) => { + LobpcgResult::Ok(vals, _, r_norms) | LobpcgResult::Err(vals, _, r_norms, _) => { // check convergence for (i, norm) in r_norms.into_iter().enumerate() { if norm > 1e-5 { @@ -501,7 +500,7 @@ mod tests { close_l2(&Array1::from(ground_truth_eigvals.to_vec()), &vals, num as f64 * 5e-4) } } - EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), + LobpcgResult::NoResult(err) => panic!("Did not converge: {:?}", err), } } @@ -539,11 +538,15 @@ mod tests { let diag = arr1(&[1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]); let a = Array2::from_diag(&diag); let x: Array2 = Array2::random((10, 1), Uniform::new(0.0, 1.0)); - let y: Array2 = arr2(&[[1.0, 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 1.0, 0., 0., 0., 0., 0., 0., 0., 0.]]).reversed_axes(); + let y: Array2 = arr2(&[ + [1.0, 0., 0., 0., 0., 0., 0., 0., 0., 0.], + [0., 1.0, 0., 0., 0., 0., 0., 0., 0., 0.], + ]) + .reversed_axes(); let result = lobpcg(|y| a.dot(&y), x, |_| {}, Some(y), 1e-10, 50, Order::Smallest); match result { - EigResult::Ok(vals, vecs, r_norms) | EigResult::Err(vals, vecs, r_norms, _) => { + LobpcgResult::Ok(vals, vecs, r_norms) | LobpcgResult::Err(vals, vecs, r_norms, _) => { // check convergence for (i, norm) in r_norms.into_iter().enumerate() { if norm > 0.01 { @@ -561,7 +564,7 @@ mod tests { 1e-5, ); } - EigResult::NoResult(err) => panic!("Did not converge: {:?}", err), + LobpcgResult::NoResult(err) => panic!("Did not converge: {:?}", err), } } } diff --git a/src/lobpcg/mod.rs b/src/lobpcg/mod.rs index 717cbc76..0c0bce47 100644 --- a/src/lobpcg/mod.rs +++ b/src/lobpcg/mod.rs @@ -3,5 +3,5 @@ mod lobpcg; mod svd; pub use eig::TruncatedEig; -pub use lobpcg::{lobpcg, EigResult, Order as TruncatedOrder}; +pub use lobpcg::{lobpcg, LobpcgResult, Order as TruncatedOrder}; pub use svd::TruncatedSvd; diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index feb55c12..8222000f 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -1,7 +1,7 @@ ///! Truncated singular value decomposition -///! -///! This module computes the k largest/smallest singular values/vectors for a dense matrix. -use super::lobpcg::{lobpcg, EigResult, Order}; +///! +///! This module computes the k largest/smallest singular values/vectors for a dense matrix. +use super::lobpcg::{lobpcg, LobpcgResult, Order}; use crate::error::Result; use crate::{Lapack, Scalar}; use ndarray::prelude::*; @@ -129,8 +129,7 @@ impl Truncate let (n, m) = (self.problem.nrows(), self.problem.ncols()); // generate initial matrix - let x = Array2::random((usize::min(n, m), num), Uniform::new(0.0, 1.0)) - .mapv(|x| NumCast::from(x).unwrap()); + 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; @@ -160,13 +159,13 @@ impl Truncate // convert into TruncatedSvdResult match res { - EigResult::Ok(vals, vecs, _) | EigResult::Err(vals, vecs, _, _) => Ok(TruncatedSvdResult { + LobpcgResult::Ok(vals, vecs, _) | LobpcgResult::Err(vals, vecs, _, _) => Ok(TruncatedSvdResult { problem: self.problem.clone(), eigvals: vals, eigvecs: vecs, ngm: n > m, }), - EigResult::NoResult(err) => Err(err), + LobpcgResult::NoResult(err) => Err(err), } } } From 6a12b2a7ee066c49df60455e8acb28cbcfc4edf4 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Wed, 25 Mar 2020 10:48:10 +0100 Subject: [PATCH 27/30] Continue if full eigendecomposition fails --- src/lobpcg/lobpcg.rs | 112 ++++++++++++++++++++++--------------------- 1 file changed, 58 insertions(+), 54 deletions(-) diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 40df72c5..9fd0fb65 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -6,8 +6,7 @@ use crate::error::{LinalgError, Result}; use crate::{cholesky::*, close_l2, eigh::*, norm::*, triangular::*}; use crate::{Lapack, Scalar}; use ndarray::prelude::*; -use ndarray::OwnedRepr; -use ndarray::ScalarOperand; +use ndarray::{OwnedRepr, ScalarOperand, Data}; use num_traits::{Float, NumCast}; /// Find largest or smallest eigenvalues @@ -32,19 +31,19 @@ pub enum LobpcgResult { } /// Solve full eigenvalue problem, sort by `order` and truncate to `size` -fn sorted_eig( - a: ArrayView2, - b: Option>, +fn sorted_eig, A: Scalar + Lapack>( + a: ArrayBase, + b: Option>, size: usize, order: &Order, ) -> Result<(Array1, Array2)> { + let n = a.len_of(Axis(0)); + let (vals, vecs) = match b { Some(b) => (a, b).eigh(UPLO::Upper).map(|x| (x.0, (x.1).0))?, _ => a.eigh(UPLO::Upper)?, }; - let n = a.len_of(Axis(0)); - Ok(match order { Order::Largest => ( vals.slice_move(s![n-size..; -1]).mapv(|x| Scalar::from_real(x)), @@ -320,55 +319,60 @@ pub fn lobpcg< .ok() }); - // compute symmetric gram matrices - let (gram_a, gram_b) = if let Some((active_p, active_ap)) = &p_ap { - let xap = x.t().dot(active_ap); - let rap = r.t().dot(active_ap); - let pap = active_p.t().dot(active_ap); - let xp = x.t().dot(active_p); - let rp = r.t().dot(active_p); - let (pap, pp) = if explicit_gram_flag { - ((&pap + &pap.t()) / two, active_p.t().dot(active_p)) - } else { - (pap, ident.clone()) - }; + // compute symmetric gram matrices and calculate solution of eigenproblem + // + // first try to compute the eigenvalue decomposition of the span{R, X, P}, + // if this fails (or the algorithm was restarted), then just use span{R, X} + let result = p_ap.as_ref() + .ok_or(LinalgError::Lapack { return_code: 1 }) + .and_then(|(active_p, active_ap)| { + let xap = x.t().dot(active_ap); + let rap = r.t().dot(active_ap); + let pap = active_p.t().dot(active_ap); + let xp = x.t().dot(active_p); + let rp = r.t().dot(active_p); + let (pap, pp) = if explicit_gram_flag { + ((&pap + &pap.t()) / two, active_p.t().dot(active_p)) + } else { + (pap, ident.clone()) + }; + + sorted_eig( + stack![ + Axis(0), + stack![Axis(1), xax, xar, xap], + stack![Axis(1), xar.t(), rar, rap], + stack![Axis(1), xap.t(), rap.t(), pap] + ], + Some(stack![ + Axis(0), + stack![Axis(1), xx, xr, xp], + stack![Axis(1), xr.t(), rr, rp], + stack![Axis(1), xp.t(), rp.t(), pp] + ]), + size_x, + &order + ) + }) + .or_else(|_| { + sorted_eig( + stack![Axis(0), stack![Axis(1), xax, xar], stack![Axis(1), xar.t(), rar]], + Some(stack![Axis(0), stack![Axis(1), xx, xr], stack![Axis(1), xr.t(), rr]]), + size_x, + &order + ) + }); - ( - stack![ - Axis(0), - stack![Axis(1), xax, xar, xap], - stack![Axis(1), xar.t(), rar, rap], - stack![Axis(1), xap.t(), rap.t(), pap] - ], - stack![ - Axis(0), - stack![Axis(1), xx, xr, xp], - stack![Axis(1), xr.t(), rr, rp], - stack![Axis(1), xp.t(), rp.t(), pp] - ], - ) - } else { - ( - stack![Axis(0), stack![Axis(1), xax, xar], stack![Axis(1), xar.t(), rar]], - stack![Axis(0), stack![Axis(1), xx, xr], stack![Axis(1), xr.t(), rr]], - ) - }; - // apply Rayleigh-Ritz method for (A - lambda) to calculate optimal expansion coefficients - let (new_lambda, eig_vecs) = match sorted_eig(gram_a.view(), Some(gram_b.view()), size_x, &order) { - Ok(x) => x, - Err(err) => { - // restart if the eigproblem decomposition failed - if previous_p_ap.is_some() { - previous_p_ap = None; - continue; - } else { - // or break if restart is not possible - break Err(err); - } - } - }; - lambda = new_lambda; + // update eigenvalues and eigenvectors (lambda is also used in the next iteration) + let eig_vecs; + match result { + Ok((x, y)) => { + lambda = x; + eig_vecs = y; + }, + Err(x) => break Err(x) + } // approximate eigenvector X and conjugate vectors P with solution of eigenproblem let (p, ap, tau) = if let Some((active_p, active_ap)) = p_ap { From 35564fa5a3b2abff2897e80afd907351aa863414 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Wed, 25 Mar 2020 11:15:43 +0100 Subject: [PATCH 28/30] Remove debugging output --- src/lobpcg/lobpcg.rs | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 9fd0fb65..101eb98b 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -414,8 +414,6 @@ pub fn lobpcg< let (vals, vecs, rnorm) = best_result.unwrap(); let rnorm = rnorm.into_iter().map(|x| Scalar::from_real(x)).collect(); - dbg!(&residual_norms_history); - match final_norm { Ok(_) => LobpcgResult::Ok(vals, vecs, rnorm), Err(err) => LobpcgResult::Err(vals, vecs, rnorm, err), From f8ad553d1e3d92b799b020bbbdef9d368a7304e6 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Wed, 1 Apr 2020 13:23:17 +0200 Subject: [PATCH 29/30] Remove generic type from precision for simpler use --- src/lobpcg/eig.rs | 6 +++--- src/lobpcg/lobpcg.rs | 3 ++- src/lobpcg/svd.rs | 6 +++--- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs index e69734e7..673c31ca 100644 --- a/src/lobpcg/eig.rs +++ b/src/lobpcg/eig.rs @@ -20,14 +20,14 @@ pub struct TruncatedEig { problem: Array2, pub constraints: Option>, preconditioner: Option>, - precision: A::Real, + precision: f32, maxiter: usize, } impl TruncatedEig { pub fn new(problem: Array2, order: Order) -> TruncatedEig { TruncatedEig { - precision: NumCast::from(1e-5).unwrap(), + precision: 1e-5, maxiter: problem.len_of(Axis(0)) * 2, preconditioner: None, constraints: None, @@ -36,7 +36,7 @@ impl Truncate } } - pub fn precision(mut self, precision: A::Real) -> Self { + pub fn precision(mut self, precision: f32) -> Self { self.precision = precision; self diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index 101eb98b..fd1ea37e 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -146,7 +146,7 @@ pub fn lobpcg< mut x: Array2, m: G, y: Option>, - tol: A::Real, + tol: f32, maxiter: usize, order: Order, ) -> LobpcgResult { @@ -166,6 +166,7 @@ pub fn lobpcg< // cap the number of iteration let mut iter = usize::min(n * 10, maxiter); + let tol = NumCast::from(tol).unwrap(); // calculate cholesky factorization of YY' and apply constraints to initial guess let cholesky_yy = y.as_ref().map(|y| { diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index 8222000f..0bdc861d 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -94,21 +94,21 @@ impl + 'static + MagnitudeCorrection> Trunc pub struct TruncatedSvd { order: Order, problem: Array2, - precision: A::Real, + precision: f32, maxiter: usize, } impl TruncatedSvd { pub fn new(problem: Array2, order: Order) -> TruncatedSvd { TruncatedSvd { - precision: NumCast::from(1e-5).unwrap(), + precision: 1e-5, maxiter: problem.len_of(Axis(0)) * 2, order, problem, } } - pub fn precision(mut self, precision: A::Real) -> Self { + pub fn precision(mut self, precision: f32) -> Self { self.precision = precision; self From ae2ce6a72c1bb473f699c07bcc0fe45eee4b2d6f Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Tue, 28 Apr 2020 11:41:02 +0200 Subject: [PATCH 30/30] Remove dependency to `ndarray-rand`; Fix restarting issue --- Cargo.toml | 1 - src/lobpcg/eig.rs | 8 +++----- src/lobpcg/lobpcg.rs | 20 ++++++++++---------- src/lobpcg/svd.rs | 14 ++++++-------- 4 files changed, 19 insertions(+), 24 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 5bd0aac7..a2e642b0 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -28,7 +28,6 @@ num-traits = "0.2" cauchy = "0.2.1" num-complex = "0.2.1" rand = "0.5" -ndarray-rand = "0.11" [dependencies.ndarray] version = "0.13" diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs index 673c31ca..6f4d1ac4 100644 --- a/src/lobpcg/eig.rs +++ b/src/lobpcg/eig.rs @@ -1,12 +1,10 @@ use super::lobpcg::{lobpcg, LobpcgResult, Order}; -use crate::{Lapack, Scalar}; +use crate::{Lapack, Scalar, generate}; ///! Implements truncated eigenvalue decomposition /// use ndarray::prelude::*; use ndarray::stack; use ndarray::ScalarOperand; -use ndarray_rand::rand_distr::Uniform; -use ndarray_rand::RandomExt; use num_traits::{Float, NumCast}; /// Truncated eigenproblem solver @@ -62,8 +60,8 @@ impl Truncate // calculate the eigenvalues decompose pub fn decompose(&self, num: usize) -> LobpcgResult { - let x = Array2::random((self.problem.len_of(Axis(0)), num), Uniform::new(0.0, 1.0)) - .mapv(|x| NumCast::from(x).unwrap()); + let x: Array2 = generate::random((self.problem.len_of(Axis(0)), num)); + let x = x.mapv(|x| NumCast::from(x).unwrap()); if let Some(ref preconditioner) = self.preconditioner { lobpcg( diff --git a/src/lobpcg/lobpcg.rs b/src/lobpcg/lobpcg.rs index fd1ea37e..b9f3e52c 100644 --- a/src/lobpcg/lobpcg.rs +++ b/src/lobpcg/lobpcg.rs @@ -304,7 +304,7 @@ pub fn lobpcg< }; // mask and orthonormalize P and AP - let p_ap = previous_p_ap + let mut p_ap = previous_p_ap .as_ref() .and_then(|(p, ap)| { let active_p = ndarray_mask(p.view(), &activemask); @@ -356,6 +356,8 @@ pub fn lobpcg< ) }) .or_else(|_| { + p_ap = None; + sorted_eig( stack![Axis(0), stack![Axis(1), xax, xar], stack![Axis(1), xar.t(), rar]], Some(stack![Axis(0), stack![Axis(1), xx, xr], stack![Axis(1), xr.t(), rr]]), @@ -431,14 +433,13 @@ mod tests { use super::Order; use crate::close_l2; use crate::qr::*; + use crate::generate; use ndarray::prelude::*; - use ndarray_rand::rand_distr::Uniform; - use ndarray_rand::RandomExt; /// Test the `sorted_eigen` function #[test] fn test_sorted_eigen() { - let matrix = Array2::random((10, 10), Uniform::new(0., 10.)); + let matrix: Array2 = generate::random((10, 10)) * 10.0; let matrix = matrix.t().dot(&matrix); // return all eigenvectors with largest first @@ -454,7 +455,7 @@ mod tests { /// Test the masking function #[test] fn test_masking() { - let matrix = Array2::random((10, 5), Uniform::new(0., 10.)); + let matrix: Array2 = generate::random((10, 5)) * 10.0; let masked_matrix = ndarray_mask(matrix.view(), &[true, true, false, true, false]); close_l2(&masked_matrix.slice(s![.., 2]), &matrix.slice(s![.., 3]), 1e-12); } @@ -462,7 +463,7 @@ mod tests { /// Test orthonormalization of a random matrix #[test] fn test_orthonormalize() { - let matrix: Array2 = Array2::random((10, 10), Uniform::new(-10., 10.)); + let matrix: Array2 = generate::random((10, 10)) * 10.0; let (n, l) = orthonormalize(matrix.clone()).unwrap(); @@ -483,10 +484,9 @@ mod tests { assert_symmetric(a); let n = a.len_of(Axis(0)); - let x: Array2 = Array2::random((n, num), Uniform::new(0.0, 1.0)); + let x: Array2 = generate::random((n, num)); let result = lobpcg(|y| a.dot(&y), x, |_| {}, None, 1e-5, n * 2, order); - dbg!(&result); match result { LobpcgResult::Ok(vals, _, r_norms) | LobpcgResult::Err(vals, _, r_norms, _) => { // check convergence @@ -523,7 +523,7 @@ mod tests { #[test] fn test_eigsolver_constructed() { let n = 50; - let tmp = Array2::random((n, n), Uniform::new(0.0, 1.0)); + let tmp = generate::random((n, n)); //let (v, _) = tmp.qr_square().unwrap(); let (v, _) = orthonormalize(tmp).unwrap(); @@ -540,7 +540,7 @@ mod tests { fn test_eigsolver_constrained() { let diag = arr1(&[1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]); let a = Array2::from_diag(&diag); - let x: Array2 = Array2::random((10, 1), Uniform::new(0.0, 1.0)); + let x: Array2 = generate::random((10, 1)); let y: Array2 = arr2(&[ [1.0, 0., 0., 0., 0., 0., 0., 0., 0., 0.], [0., 1.0, 0., 0., 0., 0., 0., 0., 0., 0.], diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index 0bdc861d..bc85f06f 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -3,11 +3,9 @@ ///! This module computes the k largest/smallest singular values/vectors for a dense matrix. use super::lobpcg::{lobpcg, LobpcgResult, Order}; use crate::error::Result; -use crate::{Lapack, Scalar}; +use crate::{Lapack, Scalar, generate}; use ndarray::prelude::*; use ndarray::ScalarOperand; -use ndarray_rand::rand_distr::Uniform; -use ndarray_rand::RandomExt; use num_traits::{Float, NumCast}; use std::ops::DivAssign; @@ -129,7 +127,8 @@ impl Truncate let (n, m) = (self.problem.nrows(), self.problem.ncols()); // generate initial matrix - let x = Array2::random((usize::min(n, m), num), Uniform::new(0.0, 1.0)).mapv(|x| NumCast::from(x).unwrap()); + let x: Array2 = generate::random((usize::min(n, m), num)); + let x = x.mapv(|x| NumCast::from(x).unwrap()); // square precision because the SVD squares the eigenvalue as well let precision = self.precision * self.precision; @@ -190,10 +189,9 @@ impl MagnitudeCorrection for f64 { mod tests { use super::Order; use super::TruncatedSvd; - use crate::close_l2; + use crate::{close_l2, generate}; + use ndarray::{arr1, arr2, Array2}; - use ndarray_rand::rand_distr::Uniform; - use ndarray_rand::RandomExt; #[test] fn test_truncated_svd() { @@ -212,7 +210,7 @@ mod tests { #[test] fn test_truncated_svd_random() { - let a: Array2 = Array2::random((50, 10), Uniform::new(0.0, 1.0)); + let a: Array2 = generate::random((50, 10)); let res = TruncatedSvd::new(a.clone(), Order::Largest) .precision(1e-5)