From d201ca13819771b7e048bd4c5bc3cd6747adc6ff Mon Sep 17 00:00:00 2001 From: YuhanLiin Date: Tue, 24 May 2022 21:40:25 -0400 Subject: [PATCH 1/8] Eliminate dependency on xoshiro --- Cargo.toml | 1 - src/lobpcg/eig.rs | 38 ++++++++++++----------------------- src/lobpcg/svd.rs | 51 ++++++++++++++++++----------------------------- 3 files changed, 32 insertions(+), 58 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 6f76611..1b7966c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -14,7 +14,6 @@ ndarray = { version = "0.15", features = ["approx"] } num-traits = "0.2.0" thiserror = "1" rand = { version = "0.8", optional=true } -rand_xoshiro = { version = "0.6", optional=true } [dev-dependencies] approx = "0.4" diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs index bae8002..1cba8f9 100644 --- a/src/lobpcg/eig.rs +++ b/src/lobpcg/eig.rs @@ -12,10 +12,6 @@ use num_traits::NumCast; use rand::Rng; use std::iter::Sum; -#[cfg(feature = "rand_xoshiro")] -use rand_xoshiro::{rand_core::SeedableRng, Xoshiro256Plus}; -//#[cfg(feature="rand_xoshiro")] - #[derive(Debug, Clone)] /// Truncated eigenproblem solver /// @@ -29,11 +25,13 @@ use rand_xoshiro::{rand_core::SeedableRng, Xoshiro256Plus}; /// ```rust /// use ndarray::{arr1, Array2}; /// use ndarray_linalg_rs::{Order, lobpcg::TruncatedEig}; +/// use rand::SeedableRng; +/// use rand_xoshiro::Xoshiro256Plus; /// /// let diag = arr1(&[1., 2., 3., 4., 5.]); /// let a = Array2::from_diag(&diag); /// -/// let mut eig = TruncatedEig::new_from_seed(a, Order::Largest, 42) +/// let mut eig = TruncatedEig::new_with_rng(a, Order::Largest, Xoshiro256Plus::seed_from_u64(42)) /// .precision(1e-5) /// .maxiter(500); /// @@ -49,22 +47,6 @@ pub struct TruncatedEig { rng: R, } -#[cfg(feature = "rand_xoshiro")] -impl TruncatedEig { - /// Create a new truncated eigenproblem solver - /// - /// # Properties - /// * `problem`: problem matrix - /// * `order`: ordering of the eigenvalues with [Order](crate::Order) - pub fn new_from_seed( - problem: Array2, - order: Order, - seed: u64, - ) -> TruncatedEig { - Self::new_with_rng(problem, order, Xoshiro256Plus::seed_from_u64(seed)) - } -} - impl TruncatedEig { /// Create a new truncated eigenproblem solver /// @@ -141,11 +123,13 @@ impl TruncatedEig { /// ```rust /// use ndarray::{arr1, Array2}; /// use ndarray_linalg_rs::{Order, lobpcg::TruncatedEig}; + /// use rand::SeedableRng; + /// use rand_xoshiro::Xoshiro256Plus; /// /// let diag = arr1(&[1., 2., 3., 4., 5.]); /// let a = Array2::from_diag(&diag); /// - /// let mut eig = TruncatedEig::new_from_seed(a, Order::Largest, 42) + /// let mut eig = TruncatedEig::new_with_rng(a, Order::Largest, Xoshiro256Plus::seed_from_u64(42)) /// .precision(1e-5) /// .maxiter(500); /// @@ -217,11 +201,13 @@ impl TruncatedEig { /// ```rust /// use ndarray::{arr1, Array2}; /// use ndarray_linalg_rs::{Order, lobpcg::TruncatedEig}; +/// use rand::SeedableRng; +/// use rand_xoshiro::Xoshiro256Plus; /// /// let diag = arr1(&[1., 2., 3., 4., 5.]); /// let a = Array2::from_diag(&diag); /// -/// let teig = TruncatedEig::new_from_seed(a, Order::Largest, 42) +/// let teig = TruncatedEig::new_with_rng(a, Order::Largest, Xoshiro256Plus::seed_from_u64(42)) /// .precision(1e-5) /// .maxiter(500); /// @@ -306,11 +292,13 @@ impl Iterator for TruncatedEigIterator { } } -#[cfg(all(test, feature = "rand_xoshiro"))] +#[cfg(test)] mod tests { use super::Order; use super::TruncatedEig; use ndarray::{arr1, Array2}; + use rand::SeedableRng; + use rand_xoshiro::Xoshiro256Plus; #[test] fn test_truncated_eig() { @@ -320,7 +308,7 @@ mod tests { ]); let a = Array2::from_diag(&diag); - let teig = TruncatedEig::new_from_seed(a, Order::Largest, 42) + let teig = TruncatedEig::new_with_rng(a, Order::Largest, Xoshiro256Plus::seed_from_u64(42)) .precision(1e-5) .maxiter(500); diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index c563f08..058c775 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -11,9 +11,6 @@ use std::iter::Sum; use rand::Rng; -#[cfg(feature = "rand_xoshiro")] -use rand_xoshiro::{rand_core::SeedableRng, Xoshiro256Plus}; - /// 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 @@ -115,23 +112,6 @@ pub struct TruncatedSvd { rng: R, } -#[cfg(feature = "rand_xoshiro")] -impl TruncatedSvd { - /// Create a new truncated SVD problem - /// - /// # Parameters - /// * `problem`: rectangular matrix which is decomposed - /// * `order`: whether to return large or small (close to zero) singular values - /// * `seed`: seed of the random number generator - pub fn new_from_seed( - problem: Array2, - order: Order, - seed: u64, - ) -> TruncatedSvd { - Self::new_with_rng(problem, order, Xoshiro256Plus::seed_from_u64(seed)) - } -} - impl TruncatedSvd { /// Create a new truncated SVD problem /// @@ -183,11 +163,13 @@ impl TruncatedSvd { /// ```rust /// use ndarray::{arr1, Array2}; /// use ndarray_linalg_rs::{Order, lobpcg::TruncatedSvd}; + /// use rand::SeedableRng; + /// use rand_xoshiro::Xoshiro256Plus; /// /// let diag = arr1(&[1., 2., 3., 4., 5.]); /// let a = Array2::from_diag(&diag); /// - /// let eig = TruncatedSvd::new_from_seed(a, Order::Largest, 42) + /// let eig = TruncatedSvd::new_with_rng(a, Order::Largest, Xoshiro256Plus::seed_from_u64(42)) /// .precision(1e-5) /// .maxiter(500); /// @@ -280,7 +262,7 @@ impl MagnitudeCorrection for f64 { } } -#[cfg(all(test, feature = "rand_xoshiro"))] +#[cfg(test)] mod tests { use super::Order; use super::TruncatedSvd; @@ -306,7 +288,7 @@ mod tests { fn test_truncated_svd() { let a = arr2(&[[3., 2., 2.], [2., 3., -2.]]); - let res = TruncatedSvd::new_from_seed(a, Order::Largest, 42) + let res = TruncatedSvd::new_with_rng(a, Order::Largest, Xoshiro256Plus::seed_from_u64(42)) .precision(1e-5) .maxiter(10) .decompose(2) @@ -321,11 +303,15 @@ mod tests { fn test_truncated_svd_random() { let a: Array2 = random((50, 10)); - let res = TruncatedSvd::new_from_seed(a.clone(), Order::Largest, 42) - .precision(1e-5) - .maxiter(10) - .decompose(10) - .unwrap(); + let res = TruncatedSvd::new_with_rng( + a.clone(), + Order::Largest, + Xoshiro256Plus::seed_from_u64(42), + ) + .precision(1e-5) + .maxiter(10) + .decompose(10) + .unwrap(); let (u, sigma, v_t) = res.values_vectors(); let reconstructed = u.dot(&Array2::from_diag(&sigma).dot(&v_t)); @@ -349,10 +335,11 @@ mod tests { // generate normal distribution random data with N >> p let data = Array2::random_using((1000, 500), StandardNormal, &mut rng) / 1000f64.sqrt(); - let res = TruncatedSvd::new_from_seed(data, Order::Largest, 42) - .precision(1e-3) - .decompose(500) - .unwrap(); + let res = + TruncatedSvd::new_with_rng(data, Order::Largest, Xoshiro256Plus::seed_from_u64(42)) + .precision(1e-3) + .decompose(500) + .unwrap(); let sv = res.values().mapv(|x: f64| x * x); From f8b139514bb1583a7ee6fdeade58faa03982f87c Mon Sep 17 00:00:00 2001 From: YuhanLiin Date: Tue, 24 May 2022 22:22:24 -0400 Subject: [PATCH 2/8] Add new lobpcg integration test and move Marchenko-Pastur into it --- Cargo.toml | 5 ++++ src/lobpcg/svd.rs | 56 ---------------------------------------- tests/lobpcg.rs | 65 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 70 insertions(+), 56 deletions(-) create mode 100644 tests/lobpcg.rs diff --git a/Cargo.toml b/Cargo.toml index 1b7966c..926706a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -25,3 +25,8 @@ rand_xoshiro = { version = "0.6" } [features] default = ["iterative"] iterative = ["rand"] + +[[test]] +name = "lobpcg" +path = "tests/lobpcg.rs" +required-features = ["iterative"] diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index 058c775..5dbf71c 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -318,60 +318,4 @@ mod tests { assert_abs_diff_eq!(&a, &reconstructed, epsilon = 1e-5); } - - /// Eigenvalue structure in high dimensions - /// - /// This test checks that the eigenvalues are following the Marchensko-Pastur law. The data is - /// standard uniformly distributed (i.e. E(x) = 0, E^2(x) = 1) and we have twice the amount of - /// data when compared to features. The probability density of the eigenvalues should then follow - /// a special densitiy function, described by the Marchenko-Pastur law. - /// - /// See also https://en.wikipedia.org/wiki/Marchenko%E2%80%93Pastur_distribution - #[test] - fn test_marchenko_pastur() { - // create random number generator - let mut rng = Xoshiro256Plus::seed_from_u64(3); - - // generate normal distribution random data with N >> p - let data = Array2::random_using((1000, 500), StandardNormal, &mut rng) / 1000f64.sqrt(); - - let res = - TruncatedSvd::new_with_rng(data, Order::Largest, Xoshiro256Plus::seed_from_u64(42)) - .precision(1e-3) - .decompose(500) - .unwrap(); - - let sv = res.values().mapv(|x: f64| x * x); - - // we have created a random spectrum and can apply the Marchenko-Pastur law - // with variance 1 and p/n = 0.5 - let (a, b) = ( - 1. * (1. - 0.5f64.sqrt()).powf(2.0), - 1. * (1. + 0.5f64.sqrt()).powf(2.0), - ); - - // check that the spectrum has correct boundaries - assert_abs_diff_eq!(b, sv[0], epsilon = 0.1); - assert_abs_diff_eq!(a, sv[sv.len() - 1], epsilon = 0.1); - - // estimate density empirical and compare with Marchenko-Pastur law - let mut i = 0; - 'outer: for th in Array1::linspace(0.1f64, 2.8, 28).slice(s![..;-1]) { - let mut count = 0; - while sv[i] >= *th { - count += 1; - i += 1; - - if i == sv.len() { - break 'outer; - } - } - - let x = th + 0.05; - let mp_law = ((b - x) * (x - a)).sqrt() / std::f64::consts::PI / x; - let empirical = count as f64 / 500. / ((2.8 - 0.1) / 28.); - - assert_abs_diff_eq!(mp_law, empirical, epsilon = 0.05); - } - } } diff --git a/tests/lobpcg.rs b/tests/lobpcg.rs new file mode 100644 index 0000000..358a439 --- /dev/null +++ b/tests/lobpcg.rs @@ -0,0 +1,65 @@ +use approx::assert_abs_diff_eq; +use ndarray::prelude::*; +use ndarray_rand::{rand_distr::StandardNormal, RandomExt}; +use proptest::prelude::*; + +use ndarray_linalg_rs::lobpcg::*; +use rand::SeedableRng; +use rand_xoshiro::Xoshiro256Plus; + +mod common; + +/// Eigenvalue structure in high dimensions +/// +/// This test checks that the eigenvalues are following the Marchensko-Pastur law. The data is +/// standard uniformly distributed (i.e. E(x) = 0, E^2(x) = 1) and we have twice the amount of +/// data when compared to features. The probability density of the eigenvalues should then follow +/// a special densitiy function, described by the Marchenko-Pastur law. +/// +/// See also https://en.wikipedia.org/wiki/Marchenko%E2%80%93Pastur_distribution +#[test] +fn test_marchenko_pastur() { + // create random number generator + let mut rng = Xoshiro256Plus::seed_from_u64(3); + + // generate normal distribution random data with N >> p + let data = Array2::random_using((1000, 500), StandardNormal, &mut rng) / 1000f64.sqrt(); + + let res = TruncatedSvd::new_with_rng(data, Order::Largest, Xoshiro256Plus::seed_from_u64(42)) + .precision(1e-3) + .decompose(500) + .unwrap(); + + let sv = res.values().mapv(|x: f64| x * x); + + // we have created a random spectrum and can apply the Marchenko-Pastur law + // with variance 1 and p/n = 0.5 + let (a, b) = ( + 1. * (1. - 0.5f64.sqrt()).powf(2.0), + 1. * (1. + 0.5f64.sqrt()).powf(2.0), + ); + + // check that the spectrum has correct boundaries + assert_abs_diff_eq!(b, sv[0], epsilon = 0.1); + assert_abs_diff_eq!(a, sv[sv.len() - 1], epsilon = 0.1); + + // estimate density empirical and compare with Marchenko-Pastur law + let mut i = 0; + 'outer: for th in Array1::linspace(0.1f64, 2.8, 28).slice(s![..;-1]) { + let mut count = 0; + while sv[i] >= *th { + count += 1; + i += 1; + + if i == sv.len() { + break 'outer; + } + } + + let x = th + 0.05; + let mp_law = ((b - x) * (x - a)).sqrt() / std::f64::consts::PI / x; + let empirical = count as f64 / 500. / ((2.8 - 0.1) / 28.); + + assert_abs_diff_eq!(mp_law, empirical, epsilon = 0.05); + } +} From 34c3b6fc421d8cbec1f9792b37ef0357f07fe747 Mon Sep 17 00:00:00 2001 From: YuhanLiin Date: Wed, 25 May 2022 00:31:41 -0400 Subject: [PATCH 3/8] Add proptest for lobpcg eig and found problematic test case --- Cargo.toml | 5 ----- src/lobpcg/mod.rs | 6 ++++-- src/lobpcg/svd.rs | 3 +-- tests/cholesky.rs | 18 +++--------------- tests/common.rs | 22 ++++++++++++++++++++++ tests/eigh.rs | 15 +++------------ tests/lobpcg.rs | 46 ++++++++++++++++++++++++++++++++++++++++++++-- 7 files changed, 77 insertions(+), 38 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 926706a..1b7966c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -25,8 +25,3 @@ rand_xoshiro = { version = "0.6" } [features] default = ["iterative"] iterative = ["rand"] - -[[test]] -name = "lobpcg" -path = "tests/lobpcg.rs" -required-features = ["iterative"] diff --git a/src/lobpcg/mod.rs b/src/lobpcg/mod.rs index fada3cd..b19f95d 100644 --- a/src/lobpcg/mod.rs +++ b/src/lobpcg/mod.rs @@ -45,8 +45,10 @@ where /// as it could be of value. If there is no result at all, then the second field is `None`. /// This happens if the algorithm fails in an early stage, for example if the matrix `A` is not SPD pub type LobpcgResult = std::result::Result, (LinalgError, Option>)>; + +#[derive(Debug, Clone, PartialEq)] pub struct Lobpcg { - eigvals: Array1, - eigvecs: Array2, + pub eigvals: Array1, + pub eigvecs: Array2, rnorm: Vec, } diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index 5dbf71c..081c58c 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -268,8 +268,7 @@ mod tests { use super::TruncatedSvd; use approx::assert_abs_diff_eq; - use ndarray::{arr1, arr2, s, Array1, Array2, NdFloat}; - use ndarray_rand::{rand_distr::StandardNormal, RandomExt}; + use ndarray::{arr1, arr2, Array2, NdFloat}; use rand::distributions::{Distribution, Standard}; use rand::SeedableRng; use rand_xoshiro::Xoshiro256Plus; diff --git a/tests/cholesky.rs b/tests/cholesky.rs index 4fc2e08..c093b90 100644 --- a/tests/cholesky.rs +++ b/tests/cholesky.rs @@ -6,18 +6,6 @@ use ndarray_linalg_rs::{cholesky::*, triangular::*}; mod common; -prop_compose! { - fn hpd_arr() - (arr in common::square_arr()) -> Array2 { - let dim = arr.nrows(); - let mut mul = arr.t().dot(&arr); - for i in 0..dim { - mul[(i, i)] += 1.0; - } - mul - } -} - fn run_cholesky_test(orig: Array2) { let chol = orig.cholesky().unwrap(); assert_abs_diff_eq!(chol.dot(&chol.t()), orig, epsilon = 1e-7); @@ -69,17 +57,17 @@ fn run_invc_test(a: Array2) { proptest! { #![proptest_config(ProptestConfig::with_cases(1000))] #[test] - fn cholesky_test(arr in hpd_arr()) { + fn cholesky_test(arr in common::hpd_arr()) { run_cholesky_test(arr) } #[test] - fn solvec_test((a, x) in common::system_of_arr(hpd_arr())) { + fn solvec_test((a, x) in common::system_of_arr(common::hpd_arr())) { run_solvec_test(a, x) } #[test] - fn invc_test(arr in hpd_arr()) { + fn invc_test(arr in common::hpd_arr()) { run_invc_test(arr) } } diff --git a/tests/common.rs b/tests/common.rs index 61b993d..12e6734 100644 --- a/tests/common.rs +++ b/tests/common.rs @@ -2,6 +2,7 @@ use std::ops::RangeInclusive; +use approx::assert_abs_diff_eq; use ndarray::prelude::*; use proptest::prelude::*; use proptest_derive::Arbitrary; @@ -49,6 +50,18 @@ prop_compose! { } } +prop_compose! { + pub fn hpd_arr() + (arr in square_arr()) -> Array2 { + let dim = arr.nrows(); + let mut mul = arr.t().dot(&arr); + for i in 0..dim { + mul[(i, i)] += 1.0; + } + mul + } +} + prop_compose! { pub fn rect_arr()(rows in DIM_RANGE, cols in DIM_RANGE) (arr in matrix(rows, cols)) -> Array2 { @@ -93,3 +106,12 @@ pub fn system_of_arr( ) }) } + +pub fn check_eigh(arr: &Array2, vals: &Array1, vecs: &Array2) { + // Original array multiplied with eigenvec should equal eigenval times eigenvec + for (i, v) in vecs.axis_iter(Axis(1)).enumerate() { + let av = arr.dot(&v); + let ev = v.mapv(|x| vals[i] * x); + assert_abs_diff_eq!(av, ev, epsilon = 1e-5); + } +} diff --git a/tests/eigh.rs b/tests/eigh.rs index 3ea115a..e88903a 100644 --- a/tests/eigh.rs +++ b/tests/eigh.rs @@ -6,15 +6,6 @@ use ndarray_linalg_rs::eigh::*; mod common; -fn check_eigh(arr: &Array2, vals: &Array1, vecs: &Array2) { - // Original array multiplied with eigenvec should equal eigenval times eigenvec - for (i, v) in vecs.axis_iter(Axis(1)).enumerate() { - let av = arr.dot(&v); - let ev = v.mapv(|x| vals[i] * x); - assert_abs_diff_eq!(av, ev, epsilon = 1e-5); - } -} - fn run_eigh_test(arr: Array2) { let n = arr.nrows(); let d = arr.eigh().unwrap(); @@ -23,7 +14,7 @@ fn run_eigh_test(arr: Array2) { // Eigenvecs should be orthogonal let s = vecs.t().dot(&vecs); assert_abs_diff_eq!(s, Array2::eye(n), epsilon = 1e-5); - check_eigh(&arr, &vals, &vecs); + common::check_eigh(&arr, &vals, &vecs); let (evals, evecs) = arr.clone().eigh_into().unwrap(); assert_abs_diff_eq!(evals, vals); @@ -33,12 +24,12 @@ fn run_eigh_test(arr: Array2) { // Check if ascending eigen is actually sorted and valid let (vals, vecs) = d.clone().sort_eig_asc(); - check_eigh(&arr, &vals, &vecs); + common::check_eigh(&arr, &vals, &vecs); assert!(vals.windows(2).into_iter().all(|w| w[0] <= w[1])); // Check if descending eigen is actually sorted and valid let (vals, vecs) = d.sort_eig_desc(); - check_eigh(&arr, &vals, &vecs); + common::check_eigh(&arr, &vals, &vecs); assert!(vals.windows(2).into_iter().all(|w| w[0] >= w[1])); } diff --git a/tests/lobpcg.rs b/tests/lobpcg.rs index 358a439..39759f6 100644 --- a/tests/lobpcg.rs +++ b/tests/lobpcg.rs @@ -2,11 +2,12 @@ use approx::assert_abs_diff_eq; use ndarray::prelude::*; use ndarray_rand::{rand_distr::StandardNormal, RandomExt}; use proptest::prelude::*; - -use ndarray_linalg_rs::lobpcg::*; use rand::SeedableRng; use rand_xoshiro::Xoshiro256Plus; +use ndarray_linalg_rs::eigh::*; +use ndarray_linalg_rs::lobpcg::*; + mod common; /// Eigenvalue structure in high dimensions @@ -63,3 +64,44 @@ fn test_marchenko_pastur() { assert_abs_diff_eq!(mp_law, empirical, epsilon = 0.05); } } + +fn run_lobpcg_eig_test(arr: Array2, num: usize, ordering: Order) { + let (eigvals, _) = arr.eigh().unwrap().sort_eig(ordering == Order::Largest); + let res = TruncatedEig::new_with_rng(arr.clone(), ordering, Xoshiro256Plus::seed_from_u64(42)) + .precision(1e-3) + .decompose(num) + .unwrap_or_else(|e| e.1.unwrap()); + + assert_abs_diff_eq!(eigvals.slice(s![..num]), res.eigvals, epsilon = 1e-5); + common::check_eigh(&arr, &res.eigvals, &res.eigvecs); +} + +fn generate_order() -> impl Strategy { + prop_oneof![Just(Order::Largest), Just(Order::Smallest)] +} + +prop_compose! { + pub fn hpd_arr_num()(arr in common::hpd_arr()) + (num in (1..arr.ncols()), arr in Just(arr)) -> (Array2, usize) { + (arr, num) + } +} + +proptest! { + #![proptest_config(ProptestConfig::with_cases(1000))] + #[test] + fn lobpcg_eig_test((arr, num) in hpd_arr_num(), ordering in generate_order()) { + run_lobpcg_eig_test(arr, num, ordering); + } +} + +#[test] +fn problematic_eig_matrix() { + let arr = array![ + [1.0, 0.0, 0.0, 0.0], + [0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 7854.796948298437, 2495.5155877621937], + [0.0, 0.0, 2495.5155877621937, 5995.696530257453] + ]; + run_lobpcg_eig_test(arr, 3, Order::Largest); +} From 659ecc144743c02bbbd427e4be739a32e1395f10 Mon Sep 17 00:00:00 2001 From: YuhanLiin Date: Wed, 25 May 2022 01:10:22 -0400 Subject: [PATCH 4/8] Add proptests for lobpcg SVD and found problematic case --- tests/lobpcg.rs | 81 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/tests/lobpcg.rs b/tests/lobpcg.rs index 39759f6..1e9e0cb 100644 --- a/tests/lobpcg.rs +++ b/tests/lobpcg.rs @@ -7,6 +7,7 @@ use rand_xoshiro::Xoshiro256Plus; use ndarray_linalg_rs::eigh::*; use ndarray_linalg_rs::lobpcg::*; +use ndarray_linalg_rs::svd::*; mod common; @@ -105,3 +106,83 @@ fn problematic_eig_matrix() { ]; run_lobpcg_eig_test(arr, 3, Order::Largest); } + +fn run_lobpcg_svd_test(arr: Array2, ordering: Order) { + let (_, s, _) = arr + .svd(false, false) + .unwrap() + .sort_svd(ordering == Order::Largest); + let (u, ts, vt) = + TruncatedSvd::new_with_rng(arr.clone(), ordering, Xoshiro256Plus::seed_from_u64(42)) + .precision(1e-3) + .maxiter(10) + .decompose(arr.ncols()) + .unwrap() + .values_vectors(); + + assert_abs_diff_eq!(s, ts, epsilon = 1e-5); + assert_abs_diff_eq!(u.dot(&Array2::from_diag(&ts)).dot(&vt), arr, epsilon = 1e-5); +} + +proptest! { + #![proptest_config(ProptestConfig::with_cases(256))] + #[test] + fn lobpcg_svd_test(arr in common::hpd_arr(), ordering in generate_order()) { + run_lobpcg_svd_test(arr, ordering); + } +} + +#[test] +fn problematic_svd_matrix() { + let arr = array![ + [ + 18703.111084031745, + 5398.592802934647, + -2798.4524863262, + 3142.0598040221316, + 10654.718971270437, + 2928.7057369452755 + ], + [ + 5398.592802934647, + 35574.82803149514, + -29613.112978401838, + -12632.782177317926, + -16546.07166801079, + -13607.176833471722 + ], + [ + -2798.4524863262, + -29613.112978401838, + 29022.408309489085, + 8718.392706824303, + 12376.7396224986, + 17995.47911319261 + ], + [ + 3142.0598040221316, + -12632.782177317926, + 8718.392706824303, + 22884.5878990548, + -598.390397885349, + -8629.726579767677 + ], + [ + 10654.718971270437, + -16546.07166801079, + 12376.7396224986, + -598.390397885349, + 27757.334483403938, + 15535.051898142627 + ], + [ + 2928.7057369452755, + -13607.176833471722, + 17995.47911319261, + -8629.726579767677, + 15535.051898142627, + 31748.677025662313 + ] + ]; + run_lobpcg_svd_test(arr, Order::Largest); +} From a0bd2053bf9bda0d4c308f2e0025a8f73386c9a3 Mon Sep 17 00:00:00 2001 From: YuhanLiin Date: Wed, 25 May 2022 01:30:47 -0400 Subject: [PATCH 5/8] Change sorting traits to use Order enum --- src/eigh.rs | 27 +++++++++++++-------------- src/lib.rs | 2 ++ src/lobpcg/algorithm.rs | 13 ++----------- src/svd.rs | 17 ++++++++--------- tests/lobpcg.rs | 7 ++----- tests/svd.rs | 2 +- 6 files changed, 28 insertions(+), 40 deletions(-) diff --git a/src/eigh.rs b/src/eigh.rs index 456676e..e0dedd3 100644 --- a/src/eigh.rs +++ b/src/eigh.rs @@ -3,7 +3,8 @@ use ndarray::{s, Array1, Array2, ArrayBase, Data, DataMut, Ix2, NdFloat}; use crate::{ - check_square, givens::GivensRotation, index::*, tridiagonal::SymmetricTridiagonal, Result, + check_square, givens::GivensRotation, index::*, tridiagonal::SymmetricTridiagonal, Order, + Result, }; fn symmetric_eig>( @@ -272,29 +273,28 @@ impl> EigValsh for ArrayBase { /// /// Will panic if shape or layout of inputs differ from eigen output, or if input contains NaN. pub trait EigSort: Sized { - fn sort_eig(self, descending: bool) -> Self; + fn sort_eig(self, order: Order) -> Self; /// Sort eigendecomposition by the eigenvalues in ascending order fn sort_eig_asc(self) -> Self { - self.sort_eig(false) + self.sort_eig(Order::Smallest) } /// Sort eigendecomposition by the eigenvalues in descending order fn sort_eig_desc(self) -> Self { - self.sort_eig(true) + self.sort_eig(Order::Largest) } } /// Implementation on output of `EigValsh` traits impl EigSort for Array1 { - fn sort_eig(mut self, descending: bool) -> Self { + fn sort_eig(mut self, order: Order) -> Self { // Panics on non-standard layouts, which is fine because our eigenvals have standard layout let slice = self.as_slice_mut().unwrap(); // Panic only happens with NaN values - if descending { - slice.sort_by(|a, b| b.partial_cmp(a).unwrap()); - } else { - slice.sort_by(|a, b| a.partial_cmp(b).unwrap()); + match order { + Order::Largest => slice.sort_by(|a, b| b.partial_cmp(a).unwrap()), + Order::Smallest => slice.sort_by(|a, b| a.partial_cmp(b).unwrap()), } self } @@ -302,14 +302,13 @@ impl EigSort for Array1 { /// Implementation on output of `Eigh` traits impl EigSort for (Array1, Array2) { - fn sort_eig(self, descending: bool) -> Self { + fn sort_eig(self, order: Order) -> Self { let (mut vals, vecs) = self; let mut value_idx: Vec<_> = vals.iter().copied().enumerate().collect(); // Panic only happens with NaN values - if descending { - value_idx.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap()); - } else { - value_idx.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap()); + match order { + Order::Largest => value_idx.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap()), + Order::Smallest => value_idx.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap()), } let mut out = Array2::zeros(vecs.dim()); diff --git a/src/lib.rs b/src/lib.rs index 86a6fac..660937a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -70,6 +70,8 @@ pub(crate) fn check_square(arr: &ArrayBase) -> Result } /// Find largest or smallest eigenvalues +/// +/// Corresponds to descending and ascending order #[derive(Debug, Copy, Clone, PartialEq, Eq)] pub enum Order { Largest, diff --git a/src/lobpcg/algorithm.rs b/src/lobpcg/algorithm.rs index 67b40ec..e07e4f1 100644 --- a/src/lobpcg/algorithm.rs +++ b/src/lobpcg/algorithm.rs @@ -31,25 +31,16 @@ fn sorted_eig( size: usize, order: Order, ) -> Result<(Array1, Array2)> { - let n = a.len_of(Axis(0)); - let res = match b { Some(b) => generalized_eig(a, b)?, _ => a.eigh_into()?, }; // sort and ensure that signs are deterministic - let (vals, vecs) = res.sort_eig(false); + let (vals, vecs) = res.sort_eig(order); let s = vecs.row(0).mapv(|x| x.signum()); let vecs = vecs * s; - - Ok(match order { - Order::Largest => ( - vals.slice_move(s![n-size..; -1]), - vecs.slice_move(s![.., n-size..; -1]), - ), - Order::Smallest => (vals.slice_move(s![..size]), vecs.slice_move(s![.., ..size])), - }) + Ok((vals.slice_move(s![..size]), vecs.slice_move(s![.., ..size]))) } /// Masks a matrix with the given `matrix` diff --git a/src/svd.rs b/src/svd.rs index 07b3643..ea64106 100644 --- a/src/svd.rs +++ b/src/svd.rs @@ -8,7 +8,7 @@ use ndarray::{s, Array1, Array2, ArrayBase, Axis, Data, DataMut, Ix2, NdFloat}; use crate::{ bidiagonal::Bidiagonal, eigh::wilkinson_shift, givens::GivensRotation, index::*, LinalgError, - Result, + Order, Result, }; fn svd>( @@ -482,29 +482,28 @@ impl> SVD for ArrayBase { /// /// Will panic if shape of inputs differs from shape of SVD output, or if input contains NaN. pub trait SvdSort: Sized { - fn sort_svd(self, descending: bool) -> Self; + fn sort_svd(self, order: Order) -> Self; /// Sort SVD decomposition by the singular values in ascending order fn sort_svd_asc(self) -> Self { - self.sort_svd(false) + self.sort_svd(Order::Smallest) } /// Sort SVD decomposition by the singular values in descending order fn sort_svd_desc(self) -> Self { - self.sort_svd(true) + self.sort_svd(Order::Largest) } } /// Implemented on the output of the `SVD` traits impl SvdSort for (Option>, Array1, Option>) { - fn sort_svd(self, descending: bool) -> Self { + fn sort_svd(self, order: Order) -> Self { let (u, mut s, vt) = self; let mut value_idx: Vec<_> = s.iter().copied().enumerate().collect(); // Panic only happens with NaN values - if descending { - value_idx.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap()); - } else { - value_idx.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap()); + match order { + Order::Largest => value_idx.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap()), + Order::Smallest => value_idx.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap()), } let apply_ordering = |arr: &Array2, ax, values_idx: &Vec<_>| { diff --git a/tests/lobpcg.rs b/tests/lobpcg.rs index 1e9e0cb..af7496c 100644 --- a/tests/lobpcg.rs +++ b/tests/lobpcg.rs @@ -67,7 +67,7 @@ fn test_marchenko_pastur() { } fn run_lobpcg_eig_test(arr: Array2, num: usize, ordering: Order) { - let (eigvals, _) = arr.eigh().unwrap().sort_eig(ordering == Order::Largest); + let (eigvals, _) = arr.eigh().unwrap().sort_eig(ordering); let res = TruncatedEig::new_with_rng(arr.clone(), ordering, Xoshiro256Plus::seed_from_u64(42)) .precision(1e-3) .decompose(num) @@ -108,10 +108,7 @@ fn problematic_eig_matrix() { } fn run_lobpcg_svd_test(arr: Array2, ordering: Order) { - let (_, s, _) = arr - .svd(false, false) - .unwrap() - .sort_svd(ordering == Order::Largest); + let (_, s, _) = arr.svd(false, false).unwrap().sort_svd(ordering); let (u, ts, vt) = TruncatedSvd::new_with_rng(arr.clone(), ordering, Xoshiro256Plus::seed_from_u64(42)) .precision(1e-3) diff --git a/tests/svd.rs b/tests/svd.rs index cf6a8c2..a7e336f 100644 --- a/tests/svd.rs +++ b/tests/svd.rs @@ -61,7 +61,7 @@ fn run_svd_test(arr: Array2) { proptest! { #![proptest_config(ProptestConfig::with_cases(1000))] #[test] - fn bidiagonal_test(arr in common::rect_arr()) { + fn svd_test(arr in common::rect_arr()) { run_svd_test(arr); } } From a780910fd5638b134ee5db5ba03198a08a07fea7 Mon Sep 17 00:00:00 2001 From: Lorenz Schmidt Date: Sun, 7 Aug 2022 21:51:42 +0200 Subject: [PATCH 6/8] Use dense eigensolver when eigenvalues exceed 1/5th of problem size --- src/lobpcg/eig.rs | 21 +++++++++++++++++++++ src/lobpcg/svd.rs | 36 +++++++++++++++++++++++++++++++----- 2 files changed, 52 insertions(+), 5 deletions(-) diff --git a/src/lobpcg/eig.rs b/src/lobpcg/eig.rs index 1cba8f9..e7106d2 100644 --- a/src/lobpcg/eig.rs +++ b/src/lobpcg/eig.rs @@ -2,6 +2,7 @@ //! use super::random; use crate::{ + eigh::{EigSort, Eigh}, lobpcg::{lobpcg, Lobpcg, LobpcgResult}, Order, Result, }; @@ -148,6 +149,26 @@ impl TruncatedEig { let x: Array2 = random((self.problem.len_of(Axis(0)), num), &mut self.rng); let x = x.mapv(|x| NumCast::from(x).unwrap()); + // use dense eigenproblem solver if more than 1/5 eigenvalues requested + if num * 5 > self.problem.nrows() { + let (eigvals, eigvecs) = self + .problem + .eigh() + .map_err(|e| (e, None))? + .sort_eig(self.order); + + let (eigvals, eigvecs) = ( + eigvals.slice_move(s![..num]), + eigvecs.slice_move(s![.., ..num]), + ); + + return Ok(Lobpcg { + eigvals, + eigvecs, + rnorm: Vec::new(), + }); + } + if let Some(ref preconditioner) = self.preconditioner { lobpcg( |y| self.problem.dot(&y), diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index 081c58c..d959590 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -2,6 +2,7 @@ ///! ///! This module computes the k largest/smallest singular values/vectors for a dense matrix. use crate::{ + eigh::{EigSort, Eigh}, lobpcg::{lobpcg, random, Lobpcg}, Order, Result, }; @@ -35,14 +36,14 @@ impl TruncatedSvdResult { 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 + //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 > &cutoff) + //.filter(|(_, x)| *x > &cutoff) .map(|(a, b)| (b.sqrt(), a)) .unzip(); @@ -188,6 +189,31 @@ impl TruncatedSvd { } let (n, m) = (self.problem.nrows(), self.problem.ncols()); + let ngm = n > m; + + // use dense eigenproblem solver if more than 1/5 eigenvalues requested + if num * 5 > n.min(m) { + let problem = if ngm { + self.problem.t().dot(&self.problem) + } else { + self.problem.dot(&self.problem.t()) + }; + + let (eigvals, eigvecs) = problem.eigh()?.sort_eig(self.order); + + let (eigvals, eigvecs) = ( + eigvals.slice_move(s![..num]), + eigvecs.slice_move(s![..num, ..]), + ); + + return Ok(TruncatedSvdResult { + eigvals, + eigvecs, + problem: self.problem, + order: self.order, + ngm, + }); + } // generate initial matrix let x: Array2 = random((usize::min(n, m), num), &mut self.rng); @@ -234,7 +260,7 @@ impl TruncatedSvd { eigvals, eigvecs, order: self.order, - ngm: n > m, + ngm, }), Err((err, None)) => Err(err), } From 81897dd42b55f970438ca7526f695a5088a3c2e2 Mon Sep 17 00:00:00 2001 From: YuhanLiin Date: Wed, 10 Aug 2022 00:12:31 -0400 Subject: [PATCH 7/8] Add new problematic eig case, this time with only 1 eigval --- tests/common.rs | 4 +- tests/eigh.rs | 6 +-- tests/lobpcg.rs | 114 ++++++++++++++++++++++++++++++++++++++++++++---- 3 files changed, 110 insertions(+), 14 deletions(-) diff --git a/tests/common.rs b/tests/common.rs index 12e6734..456fc4a 100644 --- a/tests/common.rs +++ b/tests/common.rs @@ -107,11 +107,11 @@ pub fn system_of_arr( }) } -pub fn check_eigh(arr: &Array2, vals: &Array1, vecs: &Array2) { +pub fn check_eigh(arr: &Array2, vals: &Array1, vecs: &Array2, eps: f64) { // Original array multiplied with eigenvec should equal eigenval times eigenvec for (i, v) in vecs.axis_iter(Axis(1)).enumerate() { let av = arr.dot(&v); let ev = v.mapv(|x| vals[i] * x); - assert_abs_diff_eq!(av, ev, epsilon = 1e-5); + assert_abs_diff_eq!(av, ev, epsilon = eps); } } diff --git a/tests/eigh.rs b/tests/eigh.rs index e88903a..c153cac 100644 --- a/tests/eigh.rs +++ b/tests/eigh.rs @@ -14,7 +14,7 @@ fn run_eigh_test(arr: Array2) { // Eigenvecs should be orthogonal let s = vecs.t().dot(&vecs); assert_abs_diff_eq!(s, Array2::eye(n), epsilon = 1e-5); - common::check_eigh(&arr, &vals, &vecs); + common::check_eigh(&arr, &vals, &vecs, 1e-5); let (evals, evecs) = arr.clone().eigh_into().unwrap(); assert_abs_diff_eq!(evals, vals); @@ -24,12 +24,12 @@ fn run_eigh_test(arr: Array2) { // Check if ascending eigen is actually sorted and valid let (vals, vecs) = d.clone().sort_eig_asc(); - common::check_eigh(&arr, &vals, &vecs); + common::check_eigh(&arr, &vals, &vecs, 1e-5); assert!(vals.windows(2).into_iter().all(|w| w[0] <= w[1])); // Check if descending eigen is actually sorted and valid let (vals, vecs) = d.sort_eig_desc(); - common::check_eigh(&arr, &vals, &vecs); + common::check_eigh(&arr, &vals, &vecs, 1e-5); assert!(vals.windows(2).into_iter().all(|w| w[0] >= w[1])); } diff --git a/tests/lobpcg.rs b/tests/lobpcg.rs index af7496c..805f0aa 100644 --- a/tests/lobpcg.rs +++ b/tests/lobpcg.rs @@ -67,14 +67,15 @@ fn test_marchenko_pastur() { } fn run_lobpcg_eig_test(arr: Array2, num: usize, ordering: Order) { - let (eigvals, _) = arr.eigh().unwrap().sort_eig(ordering); + let (eigvals, eigvecs) = arr.eigh().unwrap().sort_eig(ordering); let res = TruncatedEig::new_with_rng(arr.clone(), ordering, Xoshiro256Plus::seed_from_u64(42)) .precision(1e-3) .decompose(num) .unwrap_or_else(|e| e.1.unwrap()); - assert_abs_diff_eq!(eigvals.slice(s![..num]), res.eigvals, epsilon = 1e-5); - common::check_eigh(&arr, &res.eigvals, &res.eigvecs); + assert_abs_diff_eq!(eigvals.slice(s![..num]), res.eigvals, epsilon = 1e-4); + dbg!(&res, eigvecs.t(), eigvals); + common::check_eigh(&arr, &res.eigvals, &res.eigvecs, 1e-3); } fn generate_order() -> impl Strategy { @@ -99,12 +100,107 @@ proptest! { #[test] fn problematic_eig_matrix() { let arr = array![ - [1.0, 0.0, 0.0, 0.0], - [0.0, 1.0, 0.0, 0.0], - [0.0, 0.0, 7854.796948298437, 2495.5155877621937], - [0.0, 0.0, 2495.5155877621937, 5995.696530257453] + [ + 25009.24150520947, + 7711.177084952972, + 15036.339477113419, + 12415.633983053796, + -7826.364457097828, + -5794.809136983454, + 8014.50769983203, + 8393.010531289503, + -2558.5576771065075 + ], + [ + 7711.177084952972, + 23918.11648584492, + 10709.25998973371, + 10935.552386262338, + -770.6732260414797, + 2870.7995903313977, + 8982.076596422414, + -241.77063771128996, + 14084.61163861909 + ], + [ + 15036.339477113419, + 10709.25998973371, + 15325.033678523443, + 10313.402310703299, + -7271.420592402391, + -7517.875028349669, + 9638.500971592777, + 10156.061917050705, + 2838.341689108284 + ], + [ + 12415.633983053796, + 10935.552386262338, + 10313.402310703299, + 10952.911961656562, + -4536.386138355166, + 814.3228458599201, + 6535.288606359948, + 1408.8837450771002, + 4207.472717265064 + ], + [ + -7826.364457097828, + -770.6732260414797, + -7271.420592402391, + -4536.386138355166, + 11258.72499230996, + 5527.44886857787, + -3382.5562839168374, + -6540.344015025949, + 4230.603271776603 + ], + [ + -5794.809136983454, + 2870.7995903313977, + -7517.875028349669, + 814.3228458599201, + 5527.44886857787, + 13254.282602430641, + -3862.930728271387, + -13758.627064114638, + 5570.374917816143 + ], + [ + 8014.50769983203, + 8982.076596422414, + 9638.500971592777, + 6535.288606359948, + -3382.5562839168374, + -3862.930728271387, + 6720.537128975621, + 5637.822486334431, + 4785.3669453016 + ], + [ + 8393.010531289503, + -241.77063771128996, + 10156.061917050705, + 1408.8837450771002, + -6540.344015025949, + -13758.627064114638, + 5637.822486334431, + 14787.598102375805, + -4326.736278042166 + ], + [ + -2558.5576771065075, + 14084.61163861909, + 2838.341689108284, + 4207.472717265064, + 4230.603271776603, + 5570.374917816143, + 4785.366945301336, + -4326.736278042166, + 15299.749581864759 + ] ]; - run_lobpcg_eig_test(arr, 3, Order::Largest); + run_lobpcg_eig_test(arr, 1, Order::Smallest); } fn run_lobpcg_svd_test(arr: Array2, ordering: Order) { @@ -117,7 +213,7 @@ fn run_lobpcg_svd_test(arr: Array2, ordering: Order) { .unwrap() .values_vectors(); - assert_abs_diff_eq!(s, ts, epsilon = 1e-5); + assert_abs_diff_eq!(s, ts, epsilon = 1e-3); assert_abs_diff_eq!(u.dot(&Array2::from_diag(&ts)).dot(&vt), arr, epsilon = 1e-5); } From e4d8256e0491ba96a00fcaba4ad0b2793ce19958 Mon Sep 17 00:00:00 2001 From: YuhanLiin Date: Wed, 10 Aug 2022 00:47:24 -0400 Subject: [PATCH 8/8] Restore SVD cutoff --- src/lobpcg/svd.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/lobpcg/svd.rs b/src/lobpcg/svd.rs index 7880549..406f6f6 100644 --- a/src/lobpcg/svd.rs +++ b/src/lobpcg/svd.rs @@ -36,14 +36,14 @@ impl TruncatedSvdResult { 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 + 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 > &cutoff) + .filter(|(_, x)| *x > &cutoff) .map(|(a, b)| (b.sqrt(), a)) .unzip();