diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index 6ff36eb..8c55154 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -3,13 +3,20 @@ name: build on: [push] jobs: - build: - + build_debug: runs-on: ubuntu-latest - steps: - uses: actions/checkout@v1 - name: Build run: cargo build --verbose - name: Run tests run: cargo test --verbose + + build_release: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v1 + - name: Build + run: cargo build --release --verbose + - name: Run tests + run: cargo test --verbose diff --git a/CHANGELOG.md b/CHANGELOG.md index 0e59896..a21d0b5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,15 @@ ## [Unreleased] + +## [0.2.0] + +## Changed +* Used [modified Gram-Schmidt orthogonalization](https://github.com/felipeZ/eigenvalues/issues/8) +* [Reduce scaling of the search space](https://github.com/felipeZ/eigenvalues/issues/10) + +### Fixed +* Bug in DPR correction. + ## [0.1.4] ### Changed diff --git a/Cargo.toml b/Cargo.toml index dad5911..56a529d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "eigenvalues" -version = "0.1.4" +version = "0.2.0" authors = ["felipe "] description = "algorithms to compute eigenvalue/eigenvectors of symmetric matrices" readme = "README.md" diff --git a/src/algorithms/davidson.rs b/src/algorithms/davidson.rs index 910450e..b61b580 100644 --- a/src/algorithms/davidson.rs +++ b/src/algorithms/davidson.rs @@ -22,8 +22,9 @@ extern crate nalgebra as na; use super::SpectrumTarget; use crate::matrix_operations::MatrixOperations; use crate::utils; +use crate::MGS; use na::linalg::SymmetricEigen; -use na::{DMatrix, DVector}; +use na::{DMatrix, DVector, Dynamic}; use std::f64; use std::ops::Not; @@ -42,7 +43,13 @@ impl Config { /// * `dim` - dimension of the matrix to diagonalize /// * `method` - Either DPR or GJD /// * `target` Lowest, highest or somewhere in the middle portion of the spectrum - fn new(nvalues: usize, dim: usize, method: &str, target: SpectrumTarget) -> Self { + fn new( + nvalues: usize, + dim: usize, + method: &str, + target: SpectrumTarget, + tolerance: f64, + ) -> Self { let max_search_space = if nvalues * 10 < dim { nvalues * 10 } else { @@ -60,8 +67,8 @@ impl Config { Config { method: String::from(method), spectrum_target: target, - tolerance: 1e-6, - max_iters: 100, + tolerance: tolerance, + max_iters: 50, max_search_space: max_search_space, init_dim: nvalues * 2, } @@ -84,36 +91,38 @@ impl EigenDavidson { nvalues: usize, method: &str, spectrum_target: SpectrumTarget, + tolerance: f64, ) -> Result { // Initial configuration - let conf = Config::new(nvalues, h.rows(), method, spectrum_target); + let conf = Config::new(nvalues, h.rows(), method, spectrum_target, tolerance); // Initial subpace let mut dim_sub = conf.init_dim; // 1.1 Select the initial ortogonal subspace based on lowest elements - let mut basis = generate_subspace(&h.diagonal(), &conf); + let mut basis = Self::generate_subspace(&h.diagonal(), &conf); // 1.2 Select the correction to use let corrector = CorrectionMethod::::new(&h, &conf.method); + // 2. Generate subpace matrix problem by projecting into the basis + let first_subspace = basis.columns(0, dim_sub); + let mut matrix_subspace = h.matrix_matrix_prod(first_subspace); + let mut matrix_proj = first_subspace.transpose() * &matrix_subspace; + // Outer loop block Davidson schema let mut result = Err("Algorithm didn't converge!"); for i in 0..conf.max_iters { - // 2. Generate subpace matrix problem by projecting into the basis - let subspace = basis.columns(0, dim_sub); - let matrix_proj = subspace.transpose() * &h.matrix_matrix_prod(subspace); // (&h * subspace); - // 3. compute the eigenvalues and their corresponding ritz_vectors let ord_sort = match conf.spectrum_target { SpectrumTarget::Highest => false, _ => true, }; - let eig = utils::sort_eigenpairs(SymmetricEigen::new(matrix_proj), ord_sort); + let eig = utils::sort_eigenpairs(SymmetricEigen::new(matrix_proj.clone()), ord_sort); // 4. Check for convergence // 4.1 Compute the residues - let ritz_vectors = subspace * eig.eigenvectors.columns(0, dim_sub); - let residues = compute_residues(&h, &eig.eigenvalues, &ritz_vectors); + let ritz_vectors = basis.columns(0, dim_sub) * eig.eigenvectors.columns(0, dim_sub); + let residues = Self::compute_residues(&ritz_vectors, &matrix_subspace, &eig); // 4.2 Check Converge for each pair eigenvalue/eigenvector let errors = DVector::::from_iterator( @@ -123,31 +132,61 @@ impl EigenDavidson { .column_iter() .map(|col| col.norm()), ); - // 4.3 Check if all eigenvalues/eigenvectors have converged if errors.iter().all(|&x| x < conf.tolerance) { - result = Ok(create_results(&eig.eigenvalues, &ritz_vectors, nvalues)); + result = Ok(Self::create_results( + &eig.eigenvalues, + &ritz_vectors, + nvalues, + )); break; } - // 5. Update subspace basis set // 5.1 Add the correction vectors to the current basis - if 2 * dim_sub <= conf.max_search_space { + if dim_sub + conf.init_dim <= conf.max_search_space { let correction = corrector.compute_correction(residues, &eig.eigenvalues, &ritz_vectors); - update_subspace(&mut basis, correction, dim_sub, dim_sub * 2); + update_subspace(&mut basis, correction, (dim_sub, dim_sub + conf.init_dim)); // 6. Orthogonalize the subspace - basis = orthogonalize_subspace(basis); + MGS::orthonormalize(&mut basis, dim_sub); + + // Update projected matrix + matrix_subspace = { + let mut tmp = matrix_subspace.insert_columns(dim_sub, conf.init_dim, 0.0); + let new_block = h.matrix_matrix_prod(basis.columns(dim_sub, conf.init_dim)); + let mut slice = tmp.columns_mut(dim_sub, conf.init_dim); + slice.copy_from(&new_block); + tmp + }; + + matrix_proj = { + let new_dim = dim_sub + conf.init_dim; + let new_subspace = basis.columns(0, new_dim); + let mut tmp = DMatrix::::zeros(new_dim, new_dim); + let mut slice = tmp.index_mut((..dim_sub, ..dim_sub)); + slice.copy_from(&matrix_proj); + let new_block = + new_subspace.transpose() * matrix_subspace.columns(dim_sub, conf.init_dim); + let mut slice = tmp.index_mut((.., dim_sub..)); + slice.copy_from(&new_block); + let mut slice = tmp.index_mut((dim_sub.., ..)); + slice.copy_from(&new_block.transpose()); + tmp + }; + // update counter - dim_sub *= 2; + dim_sub += conf.init_dim; // 5.2 Otherwise reduce the basis of the subspace to the current // correction } else { dim_sub = conf.init_dim; basis.fill(0.0); - update_subspace(&mut basis, ritz_vectors, 0, dim_sub); + update_subspace(&mut basis, ritz_vectors, (0, dim_sub)); + // Update projected matrix + matrix_subspace = h.matrix_matrix_prod(basis.columns(0, dim_sub)); + matrix_proj = basis.columns(0, dim_sub).transpose() * &matrix_subspace; } // Check number of iterations if i > conf.max_iters { @@ -156,26 +195,62 @@ impl EigenDavidson { } result } -} -/// Extract the requested eigenvalues/eigenvectors pairs -fn create_results( - subspace_eigenvalues: &DVector, - ritz_vectors: &DMatrix, - nvalues: usize, -) -> EigenDavidson { - let eigenvectors = DMatrix::::from_iterator( - ritz_vectors.nrows(), - nvalues, - ritz_vectors.columns(0, nvalues).iter().cloned(), - ); - let eigenvalues = DVector::::from_iterator( - nvalues, - subspace_eigenvalues.rows(0, nvalues).iter().cloned(), - ); - EigenDavidson { - eigenvalues, - eigenvectors, + /// Extract the requested eigenvalues/eigenvectors pairs + fn create_results( + subspace_eigenvalues: &DVector, + ritz_vectors: &DMatrix, + nvalues: usize, + ) -> EigenDavidson { + let eigenvectors = DMatrix::::from_iterator( + ritz_vectors.nrows(), + nvalues, + ritz_vectors.columns(0, nvalues).iter().cloned(), + ); + let eigenvalues = DVector::::from_iterator( + nvalues, + subspace_eigenvalues.rows(0, nvalues).iter().cloned(), + ); + EigenDavidson { + eigenvalues, + eigenvectors, + } + } + + /// Residue vectors + fn compute_residues( + ritz_vectors: &DMatrix, + matrix_subspace: &DMatrix, + eig: &SymmetricEigen, + ) -> DMatrix { + let dim_sub = eig.eigenvalues.nrows(); + let lambda = { + let mut tmp = DMatrix::::zeros(dim_sub, dim_sub); + tmp.set_diagonal(&eig.eigenvalues); + tmp + }; + let vs = matrix_subspace * &eig.eigenvectors; + let guess = ritz_vectors * lambda; + vs - guess + } + + /// Generate initial orthonormal subspace + fn generate_subspace(diag: &DVector, conf: &Config) -> DMatrix { + if is_sorted(diag) && conf.spectrum_target == SpectrumTarget::Lowest { + DMatrix::::identity(diag.nrows(), conf.max_search_space) + } else { + let xs = diag.as_slice().to_vec(); + let mut rs = xs.clone(); + + // update the matrix according to the spectrumtarget + sort_diagonal(&mut rs, &conf); + let mut mtx = DMatrix::::zeros(diag.nrows(), conf.max_search_space); + for i in 0..conf.max_search_space { + let index = rs.iter().position(|&x| x == xs[i]).unwrap(); + mtx[(i, index)] = 1.0; + } + mtx + } } } @@ -223,10 +298,10 @@ where ) -> DMatrix { let d = self.target.diagonal(); let mut correction = DMatrix::::zeros(self.target.rows(), residues.ncols()); - for (k, r) in eigenvalues.iter().enumerate() { - let rs = DVector::::repeat(self.target.rows(), *r); - let x = residues.column(k).component_mul(&(rs - &d)); - correction.set_column(k, &x); + for (k, lambda) in eigenvalues.iter().enumerate() { + let tmp = DVector::::repeat(self.target.rows(), *lambda) - &d; + let rs = residues.column(k).component_div(&tmp); + correction.set_column(k, &rs); } correction } @@ -243,11 +318,13 @@ where let id = DMatrix::::identity(dimx, dimx); let ones = DVector::::repeat(dimx, 1.0); let mut correction = DMatrix::::zeros(dimx, dimy); + let diag = self.target.diagonal(); for (k, r) in ritz_vectors.column_iter().enumerate() { // Create the components of the linear system let t1 = &id - r * r.transpose(); let mut t2 = self.target.clone(); - t2.set_diagonal(&(eigenvalues[k] * &ones)); + let val = &diag - &(eigenvalues[k] * &ones); + t2.set_diagonal(&val); let arr = &t1 * &t2.matrix_matrix_prod(t1.rows(0, dimx)); // Solve the linear system let decomp = arr.lu(); @@ -260,67 +337,28 @@ where } /// Update the subpace with new vectors -fn update_subspace(basis: &mut DMatrix, vectors: DMatrix, start: usize, end: usize) { - let mut i = 0; // indices for the new vector to add - for k in start..end { - basis.set_column(k, &vectors.column(i)); - i += 1; - } -} - -/// Orthogonalize the subpsace using the QR method -fn orthogonalize_subspace(basis: DMatrix) -> DMatrix { - let qr = na::linalg::QR::new(basis); - qr.q() +fn update_subspace(basis: &mut DMatrix, vectors: DMatrix, range: (usize, usize)) { + let (start, end) = range; + let mut slice = basis.index_mut((.., start..end)); + slice.copy_from(&vectors.columns(0, end - start)); } -/// Residue vectors -fn compute_residues( - h: &M, - eigenvalues: &DVector, - ritz_vectors: &DMatrix, -) -> DMatrix { - let dim_sub = eigenvalues.nrows(); - let mut residues = DMatrix::::zeros(h.rows(), dim_sub); - for k in 0..dim_sub { - let guess = eigenvalues[k] * ritz_vectors.column(k); - let vs = h.matrix_vector_prod(ritz_vectors.column(k)); - residues.set_column(k, &(vs - guess)); - } - residues -} - -/// Generate initial orthonormal subspace -fn generate_subspace(diag: &DVector, conf: &Config) -> DMatrix { - if is_sorted(diag) { - DMatrix::::identity(diag.nrows(), conf.max_search_space) - } else { - let xs = diag.as_slice().to_vec(); - let mut rs = xs.clone(); - - match conf.spectrum_target { - SpectrumTarget::Lowest => utils::sort_vector(&mut rs, true), - SpectrumTarget::Highest => utils::sort_vector(&mut rs, false), - _ => panic!("Not implemented error!"), - } - - // update the matrix according to the spectrumtarget - let mut mtx = DMatrix::::zeros(diag.nrows(), conf.max_search_space); - for i in 0..conf.max_search_space { - let index = rs.iter().position(|&x| x == xs[i]).unwrap(); - mtx[(i, index)] = 1.0; - } - mtx +fn sort_diagonal(rs: &mut Vec, conf: &Config) { + match conf.spectrum_target { + SpectrumTarget::Lowest => utils::sort_vector(rs, true), + SpectrumTarget::Highest => utils::sort_vector(rs, false), + _ => panic!("Not implemented error!"), } } /// Check if a vector is sorted in ascending order fn is_sorted(xs: &DVector) -> bool { - let mut d: Vec = xs.iter().cloned().collect(); - d.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap()); - let vs: DVector = DVector::::from_vec(d); - let r = xs - vs; - r.norm() < f64::EPSILON + for k in 1..xs.len() { + if xs[k] < xs[k - 1] { + return false; + } + } + return true; } #[cfg(test)] @@ -332,7 +370,7 @@ mod test { fn test_update_subspace() { let mut arr = DMatrix::::repeat(3, 3, 1.); let brr = DMatrix::::zeros(3, 2); - super::update_subspace(&mut arr, brr, 0, 2); + super::update_subspace(&mut arr, brr, (0, 2)); assert_eq!(arr.column(1).sum(), 0.); assert_eq!(arr.column(2).sum(), 3.); } diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs index 1d050a8..6366eb7 100644 --- a/src/algorithms/mod.rs +++ b/src/algorithms/mod.rs @@ -7,7 +7,7 @@ pub mod davidson; /// Option to compute the lowest, highest or somewhere in the middle part of the /// spectrum -#[derive(Clone)] +#[derive(Clone,PartialEq)] pub enum SpectrumTarget { Lowest, Highest, diff --git a/src/lib.rs b/src/lib.rs index c08d355..7157ee8 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -16,15 +16,16 @@ use eigenvalues::SpectrumTarget; use na::{DMatrix, DVector}; // Generate random symmetric matrix -let brr = eigenvalues::utils::generate_diagonal_dominant(10, 0.005); +let brr = eigenvalues::utils::generate_diagonal_dominant(20, 0.005); +let tolerance = 1e-4; // Compute the first 2 lowest eigenvalues/eigenvectors using the DPR method -let eig = EigenDavidson::new (brr.clone(), 2, "DPR", SpectrumTarget::Lowest).unwrap(); +let eig = EigenDavidson::new (brr.clone(), 2, "DPR", SpectrumTarget::Lowest, tolerance).unwrap(); println!("eigenvalues:{}", eig.eigenvalues); println!("eigenvectors:{}", eig.eigenvectors); // Compute the first 2 highest eigenvalues/eigenvectors using the GJD method -let eig = EigenDavidson::new (brr, 2, "GJD", SpectrumTarget::Highest).unwrap(); +let eig = EigenDavidson::new (brr, 2, "GJD", SpectrumTarget::Highest, tolerance).unwrap(); println!("eigenvalues:{}", eig.eigenvalues); println!("eigenvectors:{}", eig.eigenvectors); ``` @@ -33,5 +34,7 @@ println!("eigenvectors:{}", eig.eigenvectors); pub mod algorithms; pub mod matrix_operations; +pub mod modified_gram_schmidt; pub mod utils; pub use algorithms::{davidson, SpectrumTarget}; +pub use modified_gram_schmidt::MGS; diff --git a/src/main.rs b/src/main.rs index 58073d4..30ddc67 100644 --- a/src/main.rs +++ b/src/main.rs @@ -5,8 +5,9 @@ use eigenvalues::davidson::EigenDavidson; use eigenvalues::SpectrumTarget; fn main() { - let brr = eigenvalues::utils::generate_diagonal_dominant(10, 0.05); - let eig = EigenDavidson::new(brr, 2, "GJD", SpectrumTarget::Lowest).unwrap(); + let brr = eigenvalues::utils::generate_diagonal_dominant(20, 0.05); + let tolerance = 1e-6; + let eig = EigenDavidson::new(brr, 2, "GJD", SpectrumTarget::Highest, tolerance).unwrap(); println!("eigenvalues:{}", eig.eigenvalues); println!("eigenvectors:{}", eig.eigenvectors); } diff --git a/src/modified_gram_schmidt.rs b/src/modified_gram_schmidt.rs new file mode 100644 index 0000000..27be78d --- /dev/null +++ b/src/modified_gram_schmidt.rs @@ -0,0 +1,65 @@ +/*! + +# Modified Gram-Schmidt (MGS) + +The Gram-Schmidt method is a method for orthonormalising a set of vectors. see: +[Gram-Schmidt process](https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process) +The MGS method improves the orthogonality loss due to the finite numerical precision +on computers. + */ + +extern crate nalgebra as na; +use na::{DMatrix, DVector, DVectorSlice}; + +pub struct MGS { + pub basis: DMatrix, +} + +impl MGS { + /// The orthonormalize static method takes three argument: + /// * `vectors` to diagonalize as columns of the matrix + /// * `start` index of the column to start orthogonalizing + /// * `end` last index of the column to diagonalize (non-inclusive) + pub fn orthonormalize(basis: &mut DMatrix, start: usize) { + let end = basis.ncols(); + for i in start..end { + for j in 0..i { + let proj = MGS::project(&basis.column(j), &basis.column(i)); + basis.set_column(i, &(basis.column(i) - proj)); + } + basis.set_column(i, &basis.column(i).normalize()); + } + } + + // Project + fn project(v1: &DVectorSlice, v2: &DVectorSlice) -> DVector { + let magnitud = v1.dot(&v2) / v1.dot(&v1); + return v1 * magnitud; + } +} + +#[cfg(test)] +mod test { + extern crate nalgebra as na; + use na::DMatrix; + + #[test] + fn test_gram_schmidt() { + let dim = 20; + let vectors = DMatrix::::new_random(dim, dim); + fun_test(vectors, 0); + } + + #[test] + fn test_start_gram_schmidt() { + let arr = DMatrix::::from_vec(3, 3, vec![1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 2.0]); + fun_test(arr, 1); + } + + fn fun_test(vectors: DMatrix, start: usize) { + let mut basis = vectors.clone(); + super::MGS::orthonormalize(&mut basis, start); + let result = basis.transpose() * &basis; + assert!(result.is_identity(1e-8)); + } +} diff --git a/tests/tests_davidson.rs b/tests/tests_davidson.rs index 07f1d94..fabb633 100644 --- a/tests/tests_davidson.rs +++ b/tests/tests_davidson.rs @@ -12,12 +12,13 @@ use eigenvalues::SpectrumTarget; fn test_davidson_lowest() { let arr = generate_diagonal_dominant(10, 0.005); let eig = sort_eigenpairs(na::linalg::SymmetricEigen::new(arr.clone()), true); - let spectrum_target = SpectrumTarget::Lowest; + let tolerance = 1.0e-4; - let dav_eig = EigenDavidson::new(arr.clone(), 2, "DPR", spectrum_target.clone()).unwrap(); + let dav_eig = + EigenDavidson::new(arr.clone(), 2, "DPR", spectrum_target.clone(), tolerance).unwrap(); test_eigenpairs(&eig, dav_eig, 2); - let dav_eig = EigenDavidson::new(arr.clone(), 2, "GJD", spectrum_target).unwrap(); + let dav_eig = EigenDavidson::new(arr.clone(), 2, "GJD", spectrum_target, tolerance).unwrap(); test_eigenpairs(&eig, dav_eig, 2); } @@ -25,25 +26,31 @@ fn test_davidson_lowest() { fn test_davidson_unsorted() { // Test the algorithm when the diagonal is unsorted let mut arr = generate_diagonal_dominant(8, 0.005); + let tolerance = 1.0e-6; let vs = na::DVector::::from_vec(vec![3.0, 2.0, 4.0, 1.0, 5.0, 6.0, 7.0, 8.0]); arr.set_diagonal(&vs); let eig = sort_eigenpairs(na::linalg::SymmetricEigen::new(arr.clone()), true); - let dav_eig = EigenDavidson::new(arr, 1, "DPR", SpectrumTarget::Lowest).unwrap(); + let dav_eig = EigenDavidson::new(arr, 1, "DPR", SpectrumTarget::Lowest, tolerance).unwrap(); test_eigenpairs(&eig, dav_eig, 1); } #[test] fn test_davidson_highest() { // Test the compution of the highest eigenvalues - let arr = generate_diagonal_dominant(10, 0.005); + let dim = 20; + let nvalues = 2; + let tolerance = 1.0e-4; + let arr = generate_diagonal_dominant(dim, 0.005); let eig = sort_eigenpairs(na::linalg::SymmetricEigen::new(arr.clone()), false); let target = SpectrumTarget::Highest; - - let dav_eig = EigenDavidson::new(arr.clone(), 2, "DPR", target.clone()).unwrap(); - test_eigenpairs(&eig, dav_eig, 2); - let dav_eig = EigenDavidson::new(arr.clone(), 2, "GJD", target).unwrap(); - test_eigenpairs(&eig, dav_eig, 2); + println!("running DPR"); + let dav_eig = + EigenDavidson::new(arr.clone(), nvalues, "DPR", target.clone(), tolerance).unwrap(); + test_eigenpairs(&eig, dav_eig, nvalues); + println!("running GJD"); + let dav_eig = EigenDavidson::new(arr.clone(), nvalues, "GJD", target, tolerance).unwrap(); + test_eigenpairs(&eig, dav_eig, nvalues); } fn test_eigenpairs( @@ -56,7 +63,7 @@ fn test_eigenpairs( assert!(relative_eq!( reference.eigenvalues[i], dav_eig.eigenvalues[i], - epsilon = 1e-8 + epsilon = 1e-6 )); // Test Eigenvectors let x = reference.eigenvectors.column(i);