From 4f004d7ba9be4709ae98e19a71f80e3d07de2c7a Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 14 May 2019 00:41:23 +0900 Subject: [PATCH 1/6] Move mgs -> krylov/mgs --- src/{ => krylov}/mgs.rs | 0 src/krylov/mod.rs | 1 + src/lib.rs | 4 +++- 3 files changed, 4 insertions(+), 1 deletion(-) rename src/{ => krylov}/mgs.rs (100%) create mode 100644 src/krylov/mod.rs diff --git a/src/mgs.rs b/src/krylov/mgs.rs similarity index 100% rename from src/mgs.rs rename to src/krylov/mgs.rs diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs new file mode 100644 index 00000000..e0af8799 --- /dev/null +++ b/src/krylov/mod.rs @@ -0,0 +1 @@ +pub mod mgs; diff --git a/src/lib.rs b/src/lib.rs index 429183dd..36a22611 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -44,9 +44,9 @@ pub mod eigh; pub mod error; pub mod generate; pub mod inner; +pub mod krylov; pub mod lapack; pub mod layout; -pub mod mgs; pub mod norm; pub mod operator; pub mod opnorm; @@ -76,3 +76,5 @@ pub use svd::*; pub use trace::*; pub use triangular::*; pub use types::*; + +pub use krylov::mgs; From 252b70ebdd12fedeed18406fe403af82004b9a9f Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 14 May 2019 00:47:35 +0900 Subject: [PATCH 2/6] Move definitions to mod.rs --- src/krylov/mgs.rs | 39 ++------------------------------------- src/krylov/mod.rs | 42 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 37 deletions(-) diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs index 9e25385c..370acde2 100644 --- a/src/krylov/mgs.rs +++ b/src/krylov/mgs.rs @@ -1,7 +1,7 @@ //! Modified Gram-Schmit orthogonalizer -use crate::{generate::*, inner::*, norm::Norm, types::*}; -use ndarray::*; +use super::*; +use crate::{generate::*, inner::*, norm::Norm}; /// Iterative orthogonalizer using modified Gram-Schmit procedure #[derive(Debug, Clone)] @@ -12,20 +12,6 @@ pub struct MGS { q: Vec>, } -/// Q-matrix -/// -/// - Maybe **NOT** square -/// - Unitary for existing columns -/// -pub type Q = Array2; - -/// R-matrix -/// -/// - Maybe **NOT** square -/// - Upper triangle -/// -pub type R = Array2; - impl MGS { /// Create an empty orthogonalizer pub fn new(dimension: usize) -> Self { @@ -129,27 +115,6 @@ impl MGS { } } -/// 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, - - /// 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 - /// ``` - Full, -} - /// Online QR decomposition of vectors using modified Gram-Schmit algorithm pub fn mgs( iter: impl Iterator>, diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs index e0af8799..eae12bf2 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -1 +1,43 @@ +//! Krylov subspace + +use crate::types::*; +use ndarray::*; + pub mod mgs; + +pub use mgs::MGS; + +/// Q-matrix +/// +/// - Maybe **NOT** square +/// - Unitary for existing columns +/// +pub type Q = Array2; + +/// R-matrix +/// +/// - Maybe **NOT** square +/// - Upper triangle +/// +pub type R = Array2; + +/// 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, + + /// 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 + /// ``` + Full, +} From e79d4b9a06570e52c9ec2be03ed5d415608c5acc Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 14 May 2019 00:54:47 +0900 Subject: [PATCH 3/6] Add Orthogonalizer trait and general impl of online QR decomposition --- src/krylov/mod.rs | 90 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs index eae12bf2..d4de9b33 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -21,6 +21,59 @@ pub type Q = Array2; /// pub type R = Array2; +pub trait Orthogonalizer { + type Elem: Scalar; + + /// Dimension of input array + fn dim(&self) -> usize; + + /// Number of cached basis + fn len(&self) -> usize; + + /// check if the basis spans entire space + fn is_full(&self) -> bool { + self.len() == self.dim() + } + + fn is_empty(&self) -> bool { + self.len() == 0 + } + + /// Orthogonalize given vector using current basis + /// + /// Panic + /// ------- + /// - if the size of the input array mismatches to the dimension + /// + fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 + where + S: DataMut; + + /// Add new vector if the residual is larger than relative tolerance + /// + /// Returns + /// -------- + /// Coefficients to the `i`-th Q-vector + /// + /// - The size of array must be `self.len() + 1` + /// - The last element is the residual norm of input vector + /// + /// Panic + /// ------- + /// - if the size of the input array mismatches to the dimension + /// + fn append( + &mut self, + a: ArrayBase, + rtol: ::Real, + ) -> Result, Array1> + where + S: DataMut; + + /// Get Q-matrix of generated basis + fn get_q(&self) -> Q; +} + /// Strategy for linearly dependent vectors appearing in iterative QR decomposition #[derive(Clone, Copy, Debug, Eq, PartialEq)] pub enum Strategy { @@ -41,3 +94,40 @@ pub enum Strategy { /// ``` Full, } + +/// Online QR decomposition using arbitrary orthogonalizer +pub fn qr( + iter: impl Iterator>, + mut ortho: impl Orthogonalizer, + rtol: A::Real, + strategy: Strategy, +) -> (Q, R) +where + A: Scalar + Lapack, + S: Data, +{ + assert_eq!(ortho.len(), 0); + + let mut coefs = Vec::new(); + for a in iter { + match ortho.append(a.into_owned(), rtol) { + Ok(coef) => coefs.push(coef), + Err(coef) => match strategy { + Strategy::Terminate => break, + Strategy::Skip => continue, + Strategy::Full => coefs.push(coef), + }, + } + } + let n = ortho.len(); + let m = coefs.len(); + 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) +} From 3e912da4fd47ce832356331a3a38d846917c2b0e Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 14 May 2019 01:00:27 +0900 Subject: [PATCH 4/6] impl Orthogonalizer for MGS --- src/krylov/mgs.rs | 104 +++++++++++++--------------------------------- tests/mgs.rs | 2 +- 2 files changed, 31 insertions(+), 75 deletions(-) diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs index 370acde2..ba86bb34 100644 --- a/src/krylov/mgs.rs +++ b/src/krylov/mgs.rs @@ -4,6 +4,25 @@ use super::*; use crate::{generate::*, inner::*, norm::Norm}; /// Iterative orthogonalizer using modified Gram-Schmit procedure +/// +/// ```rust +/// # use ndarray::*; +/// # use ndarray_linalg::{mgs::*, krylov::*, *}; +/// 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); +/// +/// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); +/// close_l2(&coef, &array![1.0, 1.0], 1e-9); +/// +/// // 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); +/// } +/// ``` #[derive(Debug, Clone)] pub struct MGS { /// Dimension of base space @@ -20,36 +39,20 @@ impl MGS { q: Vec::new(), } } +} + +impl Orthogonalizer for MGS { + type Elem = A; - /// Dimension of input array - pub fn dim(&self) -> usize { + 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 { + fn len(&self) -> usize { self.q.len() } - /// Orthogonalize given vector using current basis - /// - /// Panic - /// ------- - /// - if the size of the input array mismatches to the dimension - /// - pub fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 + fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 where A: Lapack, S: DataMut, @@ -67,32 +70,7 @@ impl MGS { coef } - /// Add new vector if the residual is larger than relative tolerance - /// - /// ```rust - /// # use ndarray::*; - /// # 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); - /// - /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap(); - /// close_l2(&coef, &array![1.0, 1.0], 1e-9); - /// - /// // 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); - /// } - /// ``` - /// - /// Panic - /// ------- - /// - if the size of the input array mismatches to the dimension - /// - pub fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Result, Array1> + fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Result, Array1> where A: Lapack, S: Data, @@ -109,8 +87,7 @@ impl MGS { Ok(coef) } - /// Get orthogonal basis as Q matrix - pub fn get_q(&self) -> Q { + fn get_q(&self) -> Q { hstack(&self.q).unwrap() } } @@ -126,27 +103,6 @@ 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) { - Ok(coef) => coefs.push(coef), - Err(coef) => match strategy { - Strategy::Terminate => break, - Strategy::Skip => continue, - Strategy::Full => coefs.push(coef), - }, - } - } - let n = ortho.len(); - let m = coefs.len(); - 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) + let mgs = MGS::new(dim); + qr(iter, mgs, rtol, strategy) } diff --git a/tests/mgs.rs b/tests/mgs.rs index 2212d98f..0e43ea68 100644 --- a/tests/mgs.rs +++ b/tests/mgs.rs @@ -1,5 +1,5 @@ use ndarray::*; -use ndarray_linalg::{mgs::*, *}; +use ndarray_linalg::{krylov::*, mgs::*, *}; fn qr_full() { const N: usize = 5; From 98bc788f763074034dbe8bbb2eea33c625fd5a52 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 14 May 2019 01:04:01 +0900 Subject: [PATCH 5/6] Revise namesapce --- src/krylov/mgs.rs | 2 +- src/krylov/mod.rs | 4 ++-- src/lib.rs | 2 -- tests/mgs.rs | 2 +- 4 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs index ba86bb34..87e57c9b 100644 --- a/src/krylov/mgs.rs +++ b/src/krylov/mgs.rs @@ -7,7 +7,7 @@ use crate::{generate::*, inner::*, norm::Norm}; /// /// ```rust /// # use ndarray::*; -/// # use ndarray_linalg::{mgs::*, krylov::*, *}; +/// # use ndarray_linalg::{krylov::*, *}; /// 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/src/krylov/mod.rs b/src/krylov/mod.rs index d4de9b33..1991929e 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -3,9 +3,9 @@ use crate::types::*; use ndarray::*; -pub mod mgs; +mod mgs; -pub use mgs::MGS; +pub use mgs::{mgs, MGS}; /// Q-matrix /// diff --git a/src/lib.rs b/src/lib.rs index 36a22611..90f7c003 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -76,5 +76,3 @@ pub use svd::*; pub use trace::*; pub use triangular::*; pub use types::*; - -pub use krylov::mgs; diff --git a/tests/mgs.rs b/tests/mgs.rs index 0e43ea68..35c860de 100644 --- a/tests/mgs.rs +++ b/tests/mgs.rs @@ -1,5 +1,5 @@ use ndarray::*; -use ndarray_linalg::{krylov::*, mgs::*, *}; +use ndarray_linalg::{krylov::*, *}; fn qr_full() { const N: usize = 5; From 1db5de497688369ef718318ebb5e78c0d37ee127 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Tue, 14 May 2019 01:06:41 +0900 Subject: [PATCH 6/6] Fill document --- src/krylov/mod.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs index 1991929e..677b1601 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -21,6 +21,7 @@ pub type Q = Array2; /// pub type R = Array2; +/// Trait for creating orthogonal basis from iterator of arrays pub trait Orthogonalizer { type Elem: Scalar;