From dafa8bb3ef3ab7950b23e2b164ecbf5a8db2ef62 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 28 Apr 2019 19:30:13 +0900 Subject: [PATCH 01/18] Test implementation of iterative QR decomposition --- src/arnoldi.rs | 99 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 1 + 2 files changed, 100 insertions(+) create mode 100644 src/arnoldi.rs diff --git a/src/arnoldi.rs b/src/arnoldi.rs new file mode 100644 index 00000000..76676e4f --- /dev/null +++ b/src/arnoldi.rs @@ -0,0 +1,99 @@ +use crate::{generate::*, norm::Norm, types::*}; +use ndarray::*; + +#[derive(Debug, Clone)] +pub struct MGS { + dim: usize, + q: Vec>, + r: Vec>, +} + +impl MGS { + pub fn new(dim: usize) -> Self { + Self { + dim, + q: Vec::new(), + r: Vec::new(), + } + } + + pub fn dim(&self) -> usize { + self.dim + } + + pub fn len(&self) -> usize { + self.q.len() + } + + pub fn append(&mut self, a: ArrayBase) -> A::Real + where + S: Data, + { + assert_eq!(a.len(), self.dim()); + let mut a = a.into_owned(); + let mut coef = Array1::zeros(self.len() + 1); + for i in 0..self.len() { + let q = &self.q[i]; + let c = a.dot(q); + azip!(mut a, q (q) in { *a = *a - c * q } ); + coef[i] = c; + } + let nrm = a.norm_l2(); + coef[self.len()] = A::from_real(nrm); + self.r.push(coef); + azip!(mut a in { *a = *a / A::from_real(nrm) }); + self.q.push(a); + nrm + } + + pub fn get_q(&self) -> Array2 { + hstack(&self.q).unwrap() + } + + pub fn get_r(&self) -> Array2 { + let len = self.len(); + let mut r = Array2::zeros((len, len)); + for i in 0..len { + for j in 0..=i { + r[(j, i)] = self.r[i][j]; + } + } + r + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::assert::*; + + const N: usize = 5; + + #[test] + fn new() { + let mgs: MGS = MGS::new(N); + assert_eq!(mgs.dim(), N); + assert_eq!(mgs.len(), 0); + } + + #[test] + fn append_random() { + let mut mgs: MGS = MGS::new(N); + let a: Array2 = random((N, 3)); + dbg!(&a); + for col in a.axis_iter(Axis(1)) { + let res = mgs.append(col); + dbg!(res); + } + let q = mgs.get_q(); + dbg!(&q); + let r = mgs.get_r(); + dbg!(&r); + + dbg!(q.dot(&r)); + close_l2(&q.dot(&r), &a, 1e-9).unwrap(); + + dbg!(q.t().dot(&q)); + close_l2(&q.t().dot(&q), &Array2::eye(3), 1e-9).unwrap(); + } +} diff --git a/src/lib.rs b/src/lib.rs index 0aaf6527..be40eb81 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -36,6 +36,7 @@ //! - [Random matrix generators](generate/index.html) //! - [Scalar trait](types/trait.Scalar.html) +pub mod arnoldi; pub mod assert; pub mod cholesky; pub mod convert; From 9a6a100b0deba68c4b28698676d7721b8a06f5dc Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 28 Apr 2019 19:42:18 +0900 Subject: [PATCH 02/18] Add test for complex --- src/arnoldi.rs | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 76676e4f..ecd30111 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -46,10 +46,12 @@ impl MGS { nrm } + /// Get orthogonal basis as Q matrix pub fn get_q(&self) -> Array2 { hstack(&self.q).unwrap() } + /// Get each vector norm and coefficients as R matrix pub fn get_r(&self) -> Array2 { let len = self.len(); let mut r = Array2::zeros((len, len)); @@ -66,6 +68,7 @@ impl MGS { mod tests { use super::*; use crate::assert::*; + use rand::{distributions::Standard, prelude::*}; const N: usize = 5; @@ -76,10 +79,12 @@ mod tests { assert_eq!(mgs.len(), 0); } - #[test] - fn append_random() { - let mut mgs: MGS = MGS::new(N); - let a: Array2 = random((N, 3)); + fn test(rtol: A::Real) + where + Standard: Distribution, + { + let mut mgs: MGS = MGS::new(N); + let a: Array2 = crate::generate::random((N, 3)); dbg!(&a); for col in a.axis_iter(Axis(1)) { let res = mgs.append(col); @@ -91,9 +96,19 @@ mod tests { dbg!(&r); dbg!(q.dot(&r)); - close_l2(&q.dot(&r), &a, 1e-9).unwrap(); + close_l2(&q.dot(&r), &a, rtol).unwrap(); dbg!(q.t().dot(&q)); - close_l2(&q.t().dot(&q), &Array2::eye(3), 1e-9).unwrap(); + close_l2(&q.t().dot(&q), &Array2::eye(3), rtol).unwrap(); + } + + #[test] + fn test_f32() { + test::(1e-5); + } + + #[test] + fn test_c32() { + test::(1e-5); } } From 75e033c19c77e79550fcc75aa5c8ce957ebc51a0 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 28 Apr 2019 20:26:20 +0900 Subject: [PATCH 03/18] Drop rand restriction (#143) --- src/arnoldi.rs | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index ecd30111..d1cbc2a3 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -68,7 +68,6 @@ impl MGS { mod tests { use super::*; use crate::assert::*; - use rand::{distributions::Standard, prelude::*}; const N: usize = 5; @@ -79,10 +78,7 @@ mod tests { assert_eq!(mgs.len(), 0); } - fn test(rtol: A::Real) - where - Standard: Distribution, - { + fn test(rtol: A::Real) { let mut mgs: MGS = MGS::new(N); let a: Array2 = crate::generate::random((N, 3)); dbg!(&a); From 199e2d00668e6258e1b309ea28f7df2c7407dfc2 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 28 Apr 2019 21:00:05 +0900 Subject: [PATCH 04/18] Use conjugate matrix instead of transpose for test --- src/arnoldi.rs | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index d1cbc2a3..e880198f 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -94,8 +94,9 @@ mod tests { dbg!(q.dot(&r)); close_l2(&q.dot(&r), &a, rtol).unwrap(); - dbg!(q.t().dot(&q)); - close_l2(&q.t().dot(&q), &Array2::eye(3), rtol).unwrap(); + let qt: Array2<_> = conjugate(&q); + dbg!(qt.dot(&q)); + close_l2(&qt.dot(&q), &Array2::eye(3), rtol).unwrap(); } #[test] From 1092554e22be05f92236fc13173811eab00df52f Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Mon, 29 Apr 2019 23:22:48 +0900 Subject: [PATCH 05/18] Use inner product --- src/arnoldi.rs | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index e880198f..0b7e2051 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -1,4 +1,4 @@ -use crate::{generate::*, norm::Norm, types::*}; +use crate::{generate::*, inner::*, norm::Norm, types::*}; use ndarray::*; #[derive(Debug, Clone)] @@ -34,7 +34,7 @@ impl MGS { let mut coef = Array1::zeros(self.len() + 1); for i in 0..self.len() { let q = &self.q[i]; - let c = a.dot(q); + let c = q.inner(&a); azip!(mut a, q (q) in { *a = *a - c * q } ); coef[i] = c; } @@ -108,4 +108,14 @@ mod tests { fn test_c32() { test::(1e-5); } + + #[test] + fn test_f64() { + test::(1e-9); + } + + #[test] + fn test_c64() { + test::(1e-9); + } } From 86b650b66bc6c77da4187f811b43a188db44da85 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 30 Apr 2019 00:10:00 +0900 Subject: [PATCH 06/18] append_if --- src/arnoldi.rs | 65 ++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 52 insertions(+), 13 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 0b7e2051..3638eab7 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -1,5 +1,6 @@ use crate::{generate::*, inner::*, norm::Norm, types::*}; use ndarray::*; +use num_traits::Zero; #[derive(Debug, Clone)] pub struct MGS { @@ -8,6 +9,11 @@ pub struct MGS { r: Vec>, } +pub enum Dependence { + Orthogonal(A::Real), + Coefficient(Array1), +} + impl MGS { pub fn new(dim: usize) -> Self { Self { @@ -25,7 +31,16 @@ impl MGS { self.q.len() } + /// Add new vector, return residual norm pub fn append(&mut self, a: ArrayBase) -> A::Real + where + S: Data, + { + self.append_if(a, A::Real::zero()).unwrap() + } + + /// Add new vector if the residual is larger than relative tolerance + pub fn append_if(&mut self, a: ArrayBase, rtol: A::Real) -> Option where S: Data, { @@ -39,11 +54,14 @@ impl MGS { coef[i] = c; } let nrm = a.norm_l2(); + if nrm < rtol { + return None; + } coef[self.len()] = A::from_real(nrm); self.r.push(coef); azip!(mut a in { *a = *a / A::from_real(nrm) }); self.q.push(a); - nrm + Some(nrm) } /// Get orthogonal basis as Q matrix @@ -78,7 +96,7 @@ mod tests { assert_eq!(mgs.len(), 0); } - fn test(rtol: A::Real) { + fn test_append(rtol: A::Real) { let mut mgs: MGS = MGS::new(N); let a: Array2 = crate::generate::random((N, 3)); dbg!(&a); @@ -100,22 +118,43 @@ mod tests { } #[test] - fn test_f32() { - test::(1e-5); + fn append() { + test_append::(1e-5); + test_append::(1e-5); + test_append::(1e-9); + test_append::(1e-9); } - #[test] - fn test_c32() { - test::(1e-5); - } + fn test_append_if(rtol: A::Real) { + let mut mgs: MGS = MGS::new(N); + let a: Array2 = crate::generate::random((N, 8)); + dbg!(&a); + for col in a.axis_iter(Axis(1)) { + match mgs.append_if(col, rtol) { + Some(res) => { + dbg!(res); + } + None => break, + } + } + let q = mgs.get_q(); + dbg!(&q); + let r = mgs.get_r(); + dbg!(&r); - #[test] - fn test_f64() { - test::(1e-9); + dbg!(q.dot(&r)); + close_l2(&q.dot(&r), &a.slice(s![.., 0..N]), rtol).unwrap(); + + let qt: Array2<_> = conjugate(&q); + dbg!(qt.dot(&q)); + close_l2(&qt.dot(&q), &Array2::eye(N), rtol).unwrap(); } #[test] - fn test_c64() { - test::(1e-9); + fn append_if() { + test_append_if::(1e-5); + test_append_if::(1e-5); + test_append_if::(1e-9); + test_append_if::(1e-9); } } From a7543f51be5eb3208cc558225a7aae09728f43d5 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 30 Apr 2019 18:29:13 +0900 Subject: [PATCH 07/18] Do not manage R matrix in MGS --- src/arnoldi.rs | 72 ++++++++++++++++++++++++-------------------------- 1 file changed, 34 insertions(+), 38 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 3638eab7..c0f69cba 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -1,51 +1,44 @@ use crate::{generate::*, inner::*, norm::Norm, types::*}; use ndarray::*; -use num_traits::Zero; +/// Iterative orthogonalizer using modified Gram-Schmit procedure #[derive(Debug, Clone)] pub struct MGS { - dim: usize, + /// Dimension of base space + dimension: usize, + /// Basis of spanned space q: Vec>, - r: Vec>, } -pub enum Dependence { - Orthogonal(A::Real), - Coefficient(Array1), -} +pub type Residual = ArrayBase; +pub type Coefficient = Array1; impl MGS { - pub fn new(dim: usize) -> Self { + pub fn new(dimension: usize) -> Self { Self { - dim, + dimension, q: Vec::new(), - r: Vec::new(), } } pub fn dim(&self) -> usize { - self.dim + self.dimension } pub fn len(&self) -> usize { self.q.len() } - /// Add new vector, return residual norm - pub fn append(&mut self, a: ArrayBase) -> A::Real - where - S: Data, - { - self.append_if(a, A::Real::zero()).unwrap() - } - - /// Add new vector if the residual is larger than relative tolerance - pub fn append_if(&mut self, a: ArrayBase, rtol: A::Real) -> Option + /// Orthogonalize given vector using current basis + /// + /// Panic + /// ------- + /// - if the size of the input array mismaches to the dimension + pub fn orthogonalize(&self, mut a: ArrayBase) -> (Residual, Coefficient) where - S: Data, + S: DataMut, { assert_eq!(a.len(), self.dim()); - let mut a = a.into_owned(); let mut coef = Array1::zeros(self.len() + 1); for i in 0..self.len() { let q = &self.q[i]; @@ -54,32 +47,35 @@ impl MGS { coef[i] = c; } let nrm = a.norm_l2(); + coef[self.len()] = A::from_real(nrm); + (a, coef) + } + + /// Add new vector if the residual is larger than relative tolerance + /// + /// Panic + /// ------- + /// - if the size of the input array mismaches to the dimension + pub fn append_if_independent(&mut self, a: ArrayBase, rtol: A::Real) -> Option> + where + S: Data, + { + let a = a.into_owned(); + let (mut a, coef) = self.orthogonalize(a); + let nrm = coef[coef.len()].re(); if nrm < rtol { + // Linearly dependent return None; } - coef[self.len()] = A::from_real(nrm); - self.r.push(coef); azip!(mut a in { *a = *a / A::from_real(nrm) }); self.q.push(a); - Some(nrm) + Some(coef) } /// Get orthogonal basis as Q matrix pub fn get_q(&self) -> Array2 { hstack(&self.q).unwrap() } - - /// Get each vector norm and coefficients as R matrix - pub fn get_r(&self) -> Array2 { - let len = self.len(); - let mut r = Array2::zeros((len, len)); - for i in 0..len { - for j in 0..=i { - r[(j, i)] = self.r[i][j]; - } - } - r - } } #[cfg(test)] From 755f654ec7c082d43894ffa65d138b8dcfcd6ce1 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 30 Apr 2019 19:33:20 +0900 Subject: [PATCH 08/18] Rewrite tests for MGS::append --- src/arnoldi.rs | 126 ++++++++++++++++++++----------------------------- 1 file changed, 51 insertions(+), 75 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index c0f69cba..54ae6478 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -10,10 +10,25 @@ pub struct MGS { q: Vec>, } +/// Residual vector of orthogonalization pub type Residual = ArrayBase; +/// Residual vector of orthogonalization pub type Coefficient = Array1; +/// Q-matrix (unitary matrix) +pub type Q = Array2; +/// R-matrix (upper triangle) +pub type R = Array2; -impl MGS { +impl MGS { + /// Create empty linear space + /// + /// ```rust + /// # use ndarray_linalg::*; + /// const N: usize = 5; + /// let mgs = arnoldi::MGS::::new(N); + /// assert_eq!(mgs.dim(), N); + /// assert_eq!(mgs.len(), 0); + /// ``` pub fn new(dimension: usize) -> Self { Self { dimension, @@ -36,6 +51,7 @@ impl MGS { /// - if the size of the input array mismaches to the dimension pub fn orthogonalize(&self, mut a: ArrayBase) -> (Residual, Coefficient) where + A: Lapack, S: DataMut, { assert_eq!(a.len(), self.dim()); @@ -56,13 +72,27 @@ impl MGS { /// Panic /// ------- /// - if the size of the input array mismaches to the dimension - pub fn append_if_independent(&mut self, a: ArrayBase, rtol: A::Real) -> Option> + /// + /// ```rust + /// # use ndarray::*; + /// # use ndarray_linalg::*; + /// let mut mgs = arnoldi::MGS::new(3); + /// let coef = mgs.append(array![1.0, 0.0, 0.0], 1e-9).unwrap(); + /// close_l2(&coef, &array![1.0], 1e-9).unwrap(); + /// + /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); + /// close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); + /// + /// assert!(mgs.append(array![1.0, 1.0, 0.0], 1e-9).is_none()); // Cannot append dependent vector + /// ``` + pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Option> where + A: Lapack, S: Data, { let a = a.into_owned(); let (mut a, coef) = self.orthogonalize(a); - let nrm = coef[coef.len()].re(); + let nrm = coef[coef.len() - 1].re(); if nrm < rtol { // Linearly dependent return None; @@ -73,84 +103,30 @@ impl MGS { } /// Get orthogonal basis as Q matrix - pub fn get_q(&self) -> Array2 { + pub fn get_q(&self) -> Q { hstack(&self.q).unwrap() } } -#[cfg(test)] -mod tests { - use super::*; - use crate::assert::*; - - const N: usize = 5; - - #[test] - fn new() { - let mgs: MGS = MGS::new(N); - assert_eq!(mgs.dim(), N); - assert_eq!(mgs.len(), 0); - } - - fn test_append(rtol: A::Real) { - let mut mgs: MGS = MGS::new(N); - let a: Array2 = crate::generate::random((N, 3)); - dbg!(&a); - for col in a.axis_iter(Axis(1)) { - let res = mgs.append(col); - dbg!(res); +pub fn mgs(iter: impl Iterator>, dim: usize, rtol: A::Real) -> (Q, R) +where + A: Scalar + Lapack, + S: Data, +{ + let mut ortho = MGS::new(dim); + let mut coefs = Vec::new(); + for a in iter { + match ortho.append(a, rtol) { + Some(coef) => coefs.push(coef), + None => break, } - let q = mgs.get_q(); - dbg!(&q); - let r = mgs.get_r(); - dbg!(&r); - - dbg!(q.dot(&r)); - close_l2(&q.dot(&r), &a, rtol).unwrap(); - - let qt: Array2<_> = conjugate(&q); - dbg!(qt.dot(&q)); - close_l2(&qt.dot(&q), &Array2::eye(3), rtol).unwrap(); - } - - #[test] - fn append() { - test_append::(1e-5); - test_append::(1e-5); - test_append::(1e-9); - test_append::(1e-9); } - - fn test_append_if(rtol: A::Real) { - let mut mgs: MGS = MGS::new(N); - let a: Array2 = crate::generate::random((N, 8)); - dbg!(&a); - for col in a.axis_iter(Axis(1)) { - match mgs.append_if(col, rtol) { - Some(res) => { - dbg!(res); - } - None => break, - } + let m = coefs.len(); + let mut r = Array2::zeros((m, m)); + for i in 0..m { + for j in 0..=i { + r[(j, i)] = coefs[i][j]; } - let q = mgs.get_q(); - dbg!(&q); - let r = mgs.get_r(); - dbg!(&r); - - dbg!(q.dot(&r)); - close_l2(&q.dot(&r), &a.slice(s![.., 0..N]), rtol).unwrap(); - - let qt: Array2<_> = conjugate(&q); - dbg!(qt.dot(&q)); - close_l2(&qt.dot(&q), &Array2::eye(N), rtol).unwrap(); - } - - #[test] - fn append_if() { - test_append_if::(1e-5); - test_append_if::(1e-5); - test_append_if::(1e-9); - test_append_if::(1e-9); } + (ortho.get_q(), r) } From b8d0b8df7855a9b6caf6de492e60bf97e94cd92b Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 30 Apr 2019 19:42:04 +0900 Subject: [PATCH 09/18] Fix test case --- src/arnoldi.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 54ae6478..3a7baeea 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -83,7 +83,7 @@ impl MGS { /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); /// close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); /// - /// assert!(mgs.append(array![1.0, 1.0, 0.0], 1e-9).is_none()); // Cannot append dependent vector + /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_none()); // Cannot append dependent vector /// ``` pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Option> where From 10c1fd18f5b848d27d4ca028bbd266fd23b4f06a Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Fri, 3 May 2019 18:26:11 +0900 Subject: [PATCH 10/18] orthogonalize takes &mut --- src/arnoldi.rs | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 3a7baeea..9f880021 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -10,11 +10,7 @@ pub struct MGS { q: Vec>, } -/// Residual vector of orthogonalization -pub type Residual = ArrayBase; -/// Residual vector of orthogonalization -pub type Coefficient = Array1; -/// Q-matrix (unitary matrix) +/// Q-matrix (unitary) pub type Q = Array2; /// R-matrix (upper triangle) pub type R = Array2; @@ -49,7 +45,7 @@ impl MGS { /// Panic /// ------- /// - if the size of the input array mismaches to the dimension - pub fn orthogonalize(&self, mut a: ArrayBase) -> (Residual, Coefficient) + pub fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 where A: Lapack, S: DataMut, @@ -59,12 +55,12 @@ impl MGS { for i in 0..self.len() { let q = &self.q[i]; let c = q.inner(&a); - azip!(mut a, q (q) in { *a = *a - c * q } ); + azip!(mut a (&mut *a), q (q) in { *a = *a - c * q } ); coef[i] = c; } let nrm = a.norm_l2(); coef[self.len()] = A::from_real(nrm); - (a, coef) + coef } /// Add new vector if the residual is larger than relative tolerance @@ -85,13 +81,13 @@ impl MGS { /// /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_none()); // Cannot append dependent vector /// ``` - pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Option> + pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Option> where A: Lapack, S: Data, { - let a = a.into_owned(); - let (mut a, coef) = self.orthogonalize(a); + let mut a = a.into_owned(); + let coef = self.orthogonalize(&mut a); let nrm = coef[coef.len() - 1].re(); if nrm < rtol { // Linearly dependent From 6fce47a112b6e085fed45679649612a049c3d12d Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Fri, 3 May 2019 18:47:40 +0900 Subject: [PATCH 11/18] Strategy for linearly dependent vectors --- src/arnoldi.rs | 51 ++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 43 insertions(+), 8 deletions(-) diff --git a/src/arnoldi.rs b/src/arnoldi.rs index 9f880021..d9519b17 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -73,15 +73,19 @@ impl MGS { /// # use ndarray::*; /// # use ndarray_linalg::*; /// let mut mgs = arnoldi::MGS::new(3); - /// let coef = mgs.append(array![1.0, 0.0, 0.0], 1e-9).unwrap(); + /// let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); /// close_l2(&coef, &array![1.0], 1e-9).unwrap(); /// /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); /// close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); /// - /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_none()); // Cannot append dependent vector + /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); // Fail if the vector is linearly dependend + /// + /// if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) { + /// close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); // You can get coefficients of dependent vector + /// } /// ``` - pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Option> + pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Result, Array1> where A: Lapack, S: Data, @@ -91,11 +95,11 @@ impl MGS { let nrm = coef[coef.len() - 1].re(); if nrm < rtol { // Linearly dependent - return None; + return Err(coef); } azip!(mut a in { *a = *a / A::from_real(nrm) }); self.q.push(a); - Some(coef) + Ok(coef) } /// Get orthogonal basis as Q matrix @@ -104,7 +108,34 @@ impl MGS { } } -pub fn mgs(iter: impl Iterator>, dim: usize, rtol: A::Real) -> (Q, R) +/// Strategy for linearly dependent vectors appearing in iterative QR decomposition +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub enum Strategy { + /// Terminate iteration if dependent vector comes + Terminate, + + /// Skip dependent vector + Skip, + + /// Orghotonalize dependent vector without adding to Q, + /// thus R must be non-regular like following: + /// + /// ```ignore + /// x x x x x + /// 0 x x x x + /// 0 0 0 x x + /// 0 0 0 0 x + /// 0 0 0 0 0 // 0-filled to be square matrix + /// ``` + Full, +} + +pub fn mgs( + iter: impl Iterator>, + dim: usize, + rtol: A::Real, + strategy: Strategy, +) -> (Q, R) where A: Scalar + Lapack, S: Data, @@ -113,8 +144,12 @@ where let mut coefs = Vec::new(); for a in iter { match ortho.append(a, rtol) { - Some(coef) => coefs.push(coef), - None => break, + Ok(coef) => coefs.push(coef), + Err(coef) => match strategy { + Strategy::Terminate => break, + Strategy::Skip => continue, + Strategy::Full => coefs.push(coef), + }, } } let m = coefs.len(); From 3f42447dda3bb6e0023bf322b4db6f46de090028 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 4 May 2019 00:13:17 +0900 Subject: [PATCH 12/18] Add tests of msg --- src/arnoldi.rs | 12 ++++--- tests/arnoldi.rs | 83 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 91 insertions(+), 4 deletions(-) create mode 100644 tests/arnoldi.rs diff --git a/src/arnoldi.rs b/src/arnoldi.rs index d9519b17..1067c10d 100644 --- a/src/arnoldi.rs +++ b/src/arnoldi.rs @@ -130,6 +130,7 @@ pub enum Strategy { Full, } +/// Online QR decomposition of vectors using modified Gram-Schmit algorithm pub fn mgs( iter: impl Iterator>, dim: usize, @@ -152,11 +153,14 @@ where }, } } + let n = ortho.len(); let m = coefs.len(); - let mut r = Array2::zeros((m, m)); - for i in 0..m { - for j in 0..=i { - r[(j, i)] = coefs[i][j]; + let mut r = Array2::zeros((n, m).f()); + for j in 0..m { + for i in 0..n { + if i < coefs[j].len() { + r[(i, j)] = coefs[j][i]; + } } } (ortho.get_q(), r) diff --git a/tests/arnoldi.rs b/tests/arnoldi.rs new file mode 100644 index 00000000..26040094 --- /dev/null +++ b/tests/arnoldi.rs @@ -0,0 +1,83 @@ +use ndarray::*; +use ndarray_linalg::{arnoldi::*, *}; + +fn qr_full() { + const N: usize = 5; + let rtol: A::Real = A::real(1e-9); + + let a: Array2 = random((N, N)); + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + assert_close_l2!(&q.dot(&r), &a, rtol); + + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); +} + +#[test] +fn qr_full_real() { + qr_full::(); +} + +#[test] +fn qr_full_complex() { + qr_full::(); +} + +fn qr() { + const N: usize = 4; + let rtol: A::Real = A::real(1e-9); + + let a: Array2 = random((N, N / 2)); + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + assert_close_l2!(&q.dot(&r), &a, rtol); + + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol); +} + +#[test] +fn qr_real() { + qr::(); +} + +#[test] +fn qr_complex() { + qr::(); +} + +fn qr_over() { + const N: usize = 4; + let rtol: A::Real = A::real(1e-9); + + let a: Array2 = random((N, N * 2)); + + // Terminate + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let a_sub = a.slice(s![.., 0..N]); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + + // Skip + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); + let a_sub = a.slice(s![.., 0..N]); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + + // Full + let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); + assert_close_l2!(&q.dot(&r), &a, rtol); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); +} + +#[test] +fn qr_over_real() { + qr_over::(); +} + +#[test] +fn qr_over_complex() { + qr_over::(); +} From 85738df18e6551af7f6b7a32d0ed6cfc64a5fc38 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Mon, 13 May 2019 04:14:29 +0900 Subject: [PATCH 13/18] Rename arnoldi -> mgs --- src/lib.rs | 3 ++- src/{arnoldi.rs => mgs.rs} | 4 ++-- tests/{arnoldi.rs => mgs.rs} | 2 +- 3 files changed, 5 insertions(+), 4 deletions(-) rename src/{arnoldi.rs => mgs.rs} (98%) rename tests/{arnoldi.rs => mgs.rs} (98%) diff --git a/src/lib.rs b/src/lib.rs index be40eb81..3cb1b2fd 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -36,7 +36,6 @@ //! - [Random matrix generators](generate/index.html) //! - [Scalar trait](types/trait.Scalar.html) -pub mod arnoldi; pub mod assert; pub mod cholesky; pub mod convert; @@ -47,6 +46,7 @@ pub mod generate; pub mod inner; pub mod lapack; pub mod layout; +pub mod mgs; pub mod norm; pub mod operator; pub mod opnorm; @@ -66,6 +66,7 @@ pub use eigh::*; pub use generate::*; pub use inner::*; pub use layout::*; +pub use mgs::*; pub use norm::*; pub use operator::*; pub use opnorm::*; diff --git a/src/arnoldi.rs b/src/mgs.rs similarity index 98% rename from src/arnoldi.rs rename to src/mgs.rs index 1067c10d..052e5409 100644 --- a/src/arnoldi.rs +++ b/src/mgs.rs @@ -21,7 +21,7 @@ impl MGS { /// ```rust /// # use ndarray_linalg::*; /// const N: usize = 5; - /// let mgs = arnoldi::MGS::::new(N); + /// let mgs = MGS::::new(N); /// assert_eq!(mgs.dim(), N); /// assert_eq!(mgs.len(), 0); /// ``` @@ -72,7 +72,7 @@ impl MGS { /// ```rust /// # use ndarray::*; /// # use ndarray_linalg::*; - /// let mut mgs = arnoldi::MGS::new(3); + /// let mut mgs = MGS::new(3); /// let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); /// close_l2(&coef, &array![1.0], 1e-9).unwrap(); /// diff --git a/tests/arnoldi.rs b/tests/mgs.rs similarity index 98% rename from tests/arnoldi.rs rename to tests/mgs.rs index 26040094..8de94894 100644 --- a/tests/arnoldi.rs +++ b/tests/mgs.rs @@ -1,5 +1,5 @@ use ndarray::*; -use ndarray_linalg::{arnoldi::*, *}; +use ndarray_linalg::*; fn qr_full() { const N: usize = 5; From 763b42818608e12c0fda09f091c6eadf10e33910 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Mon, 13 May 2019 04:16:00 +0900 Subject: [PATCH 14/18] Revise doc test --- src/mgs.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mgs.rs b/src/mgs.rs index 052e5409..d4f21ce9 100644 --- a/src/mgs.rs +++ b/src/mgs.rs @@ -74,15 +74,15 @@ impl MGS { /// # use ndarray_linalg::*; /// let mut mgs = MGS::new(3); /// let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); - /// close_l2(&coef, &array![1.0], 1e-9).unwrap(); + /// close_l2(&coef, &array![1.0], 1e-9); /// /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); - /// close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap(); + /// close_l2(&coef, &array![1.0, 1.0], 1e-9); /// /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); // Fail if the vector is linearly dependend /// /// if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) { - /// close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); // You can get coefficients of dependent vector + /// close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9); // You can get coefficients of dependent vector /// } /// ``` pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Result, Array1> @@ -120,7 +120,7 @@ pub enum Strategy { /// Orghotonalize dependent vector without adding to Q, /// thus R must be non-regular like following: /// - /// ```ignore + /// ```text /// x x x x x /// 0 x x x x /// 0 0 0 x x From 6c92d3ccd528273fd7d16ace2942e1bcf4f02f51 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Mon, 13 May 2019 20:55:49 +0900 Subject: [PATCH 15/18] Do not re-export mgs:: namespace --- src/lib.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/lib.rs b/src/lib.rs index 3cb1b2fd..429183dd 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -66,7 +66,6 @@ pub use eigh::*; pub use generate::*; pub use inner::*; pub use layout::*; -pub use mgs::*; pub use norm::*; pub use operator::*; pub use opnorm::*; From f25067b6ea7f8cbee9c16274f63a6956b5b338b1 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Mon, 13 May 2019 20:57:53 +0900 Subject: [PATCH 16/18] Fix documents --- src/mgs.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/mgs.rs b/src/mgs.rs index d4f21ce9..8d841215 100644 --- a/src/mgs.rs +++ b/src/mgs.rs @@ -44,7 +44,8 @@ impl MGS { /// /// Panic /// ------- - /// - if the size of the input array mismaches to the dimension + /// - if the size of the input array mismatches to the dimension + /// pub fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 where A: Lapack, @@ -67,7 +68,7 @@ impl MGS { /// /// Panic /// ------- - /// - if the size of the input array mismaches to the dimension + /// - if the size of the input array mismatches to the dimension /// /// ```rust /// # use ndarray::*; @@ -79,7 +80,7 @@ impl MGS { /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); /// close_l2(&coef, &array![1.0, 1.0], 1e-9); /// - /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); // Fail if the vector is linearly dependend + /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); // Fail if the vector is linearly dependent /// /// if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) { /// close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9); // You can get coefficients of dependent vector @@ -117,15 +118,14 @@ pub enum Strategy { /// Skip dependent vector Skip, - /// Orghotonalize dependent vector without adding to Q, - /// thus R must be non-regular like following: + /// Orthogonalize dependent vector without adding to Q, + /// i.e. R must be non-square like following: /// /// ```text /// x x x x x /// 0 x x x x /// 0 0 0 x x /// 0 0 0 0 x - /// 0 0 0 0 0 // 0-filled to be square matrix /// ``` Full, } From a5114eb4c1dcd822d568f1c2a71ef785da0e3591 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Mon, 13 May 2019 22:32:01 +0900 Subject: [PATCH 17/18] Fix test namespace --- src/mgs.rs | 4 ++-- tests/mgs.rs | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mgs.rs b/src/mgs.rs index 8d841215..29219ee7 100644 --- a/src/mgs.rs +++ b/src/mgs.rs @@ -19,7 +19,7 @@ impl MGS { /// Create empty linear space /// /// ```rust - /// # use ndarray_linalg::*; + /// # use ndarray_linalg::{mgs::*, *}; /// const N: usize = 5; /// let mgs = MGS::::new(N); /// assert_eq!(mgs.dim(), N); @@ -72,7 +72,7 @@ impl MGS { /// /// ```rust /// # use ndarray::*; - /// # use ndarray_linalg::*; + /// # use ndarray_linalg::{mgs::*, *}; /// let mut mgs = MGS::new(3); /// let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); /// close_l2(&coef, &array![1.0], 1e-9); diff --git a/tests/mgs.rs b/tests/mgs.rs index 8de94894..2212d98f 100644 --- a/tests/mgs.rs +++ b/tests/mgs.rs @@ -1,5 +1,5 @@ use ndarray::*; -use ndarray_linalg::*; +use ndarray_linalg::{mgs::*, *}; fn qr_full() { const N: usize = 5; From 8cd17de8da56714e50da6131bfb22904359a30a5 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Mon, 13 May 2019 22:47:40 +0900 Subject: [PATCH 18/18] Update documents --- src/mgs.rs | 54 +++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 37 insertions(+), 17 deletions(-) diff --git a/src/mgs.rs b/src/mgs.rs index 29219ee7..9e25385c 100644 --- a/src/mgs.rs +++ b/src/mgs.rs @@ -1,3 +1,5 @@ +//! Modified Gram-Schmit orthogonalizer + use crate::{generate::*, inner::*, norm::Norm, types::*}; use ndarray::*; @@ -10,21 +12,22 @@ pub struct MGS { q: Vec>, } -/// Q-matrix (unitary) +/// Q-matrix +/// +/// - Maybe **NOT** square +/// - Unitary for existing columns +/// pub type Q = Array2; -/// R-matrix (upper triangle) + +/// R-matrix +/// +/// - Maybe **NOT** square +/// - Upper triangle +/// pub type R = Array2; impl MGS { - /// Create empty linear space - /// - /// ```rust - /// # use ndarray_linalg::{mgs::*, *}; - /// const N: usize = 5; - /// let mgs = MGS::::new(N); - /// assert_eq!(mgs.dim(), N); - /// assert_eq!(mgs.len(), 0); - /// ``` + /// Create an empty orthogonalizer pub fn new(dimension: usize) -> Self { Self { dimension, @@ -32,10 +35,24 @@ impl MGS { } } + /// Dimension of input array pub fn dim(&self) -> usize { self.dimension } + /// Number of cached basis + /// + /// ```rust + /// # use ndarray::*; + /// # use ndarray_linalg::{mgs::*, *}; + /// const N: usize = 3; + /// let mut mgs = MGS::::new(N); + /// assert_eq!(mgs.dim(), N); + /// assert_eq!(mgs.len(), 0); + /// + /// mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap(); + /// assert_eq!(mgs.len(), 1); + /// ``` pub fn len(&self) -> usize { self.q.len() } @@ -66,10 +83,6 @@ impl MGS { /// Add new vector if the residual is larger than relative tolerance /// - /// Panic - /// ------- - /// - if the size of the input array mismatches to the dimension - /// /// ```rust /// # use ndarray::*; /// # use ndarray_linalg::{mgs::*, *}; @@ -80,12 +93,19 @@ impl MGS { /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); /// close_l2(&coef, &array![1.0, 1.0], 1e-9); /// - /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); // Fail if the vector is linearly dependent + /// // Fail if the vector is linearly dependent + /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err()); /// + /// // You can get coefficients of dependent vector /// if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) { - /// close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9); // You can get coefficients of dependent vector + /// close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9); /// } /// ``` + /// + /// Panic + /// ------- + /// - if the size of the input array mismatches to the dimension + /// pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Result, Array1> where A: Lapack,