diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs new file mode 100644 index 00000000..391b18c4 --- /dev/null +++ b/src/krylov/householder.rs @@ -0,0 +1,200 @@ +//! Householder reflection +//! +//! - [Householder transformation - Wikipedia](https://en.wikipedia.org/wiki/Householder_transformation) +//! + +use super::*; +use crate::{inner::*, norm::*}; +use num_traits::One; + +/// Calc a reflactor `w` from a vector `x` +pub fn calc_reflector(x: &mut ArrayBase) -> A +where + A: Scalar + Lapack, + S: DataMut, +{ + let norm = x.norm_l2(); + let alpha = -x[0].mul_real(norm / x[0].abs()); + x[0] -= alpha; + let inv_rev_norm = A::Real::one() / x.norm_l2(); + azip!(mut a(x) in { *a = a.mul_real(inv_rev_norm)}); + alpha +} + +/// Take a reflection `P = I - 2ww^T` +/// +/// Panic +/// ------ +/// - if the size of `w` and `a` mismaches +pub fn reflect(w: &ArrayBase, a: &mut ArrayBase) +where + A: Scalar + Lapack, + S1: Data, + S2: DataMut, +{ + assert_eq!(w.len(), a.len()); + let n = a.len(); + let c = A::from(2.0).unwrap() * w.inner(&a); + for l in 0..n { + a[l] -= c * w[l]; + } +} + +/// Iterative orthogonalizer using Householder reflection +#[derive(Debug, Clone)] +pub struct Householder { + /// Dimension of orthogonalizer + dim: usize, + + /// Store Householder reflector. + /// + /// The coefficient is copied into another array, and this does not contain + v: Vec>, +} + +impl Householder { + /// Create a new orthogonalizer + pub fn new(dim: usize) -> Self { + Householder { dim, v: Vec::new() } + } + + /// Take a Reflection `P = I - 2ww^T` + fn fundamental_reflection(&self, k: usize, a: &mut ArrayBase) + where + S: DataMut, + { + assert!(k < self.v.len()); + assert_eq!(a.len(), self.dim, "Input array size mismaches to the dimension"); + reflect(&self.v[k].slice(s![k..]), &mut a.slice_mut(s![k..])); + } + + /// Take forward reflection `P = P_l ... P_1` + pub fn forward_reflection(&self, a: &mut ArrayBase) + where + S: DataMut, + { + assert!(a.len() == self.dim); + let l = self.v.len(); + for k in 0..l { + self.fundamental_reflection(k, a); + } + } + + /// Take backward reflection `P = P_1 ... P_l` + pub fn backward_reflection(&self, a: &mut ArrayBase) + where + S: DataMut, + { + assert!(a.len() == self.dim); + let l = self.v.len(); + for k in (0..l).rev() { + self.fundamental_reflection(k, a); + } + } + + fn eval_residual(&self, a: &ArrayBase) -> A::Real + where + S: Data, + { + let l = self.v.len(); + a.slice(s![l..]).norm_l2() + } +} + +impl Orthogonalizer for Householder { + type Elem = A; + + fn dim(&self) -> usize { + self.dim + } + + fn len(&self) -> usize { + self.v.len() + } + + fn coeff(&self, a: ArrayBase) -> Array1 + where + S: Data, + { + let mut a = a.into_owned(); + self.forward_reflection(&mut a); + let res = self.eval_residual(&a); + let k = self.len(); + let mut c = Array1::zeros(k + 1); + azip!(mut c(c.slice_mut(s![..k])), a(a.slice(s![..k])) in { *c = a }); + c[k] = A::from_real(res); + c + } + + fn append(&mut self, mut a: ArrayBase, rtol: A::Real) -> Result, Array1> + where + S: DataMut, + { + assert_eq!(a.len(), self.dim); + let k = self.len(); + + self.forward_reflection(&mut a); + let mut coef = Array::zeros(k + 1); + for i in 0..k { + coef[i] = a[i]; + } + if self.is_full() { + return Err(coef); // coef[k] must be zero in this case + } + + let alpha = calc_reflector(&mut a.slice_mut(s![k..])); + coef[k] = alpha; + + if alpha.abs() < rtol { + // linearly dependent + return Err(coef); + } + self.v.push(a.into_owned()); + Ok(coef) + } + + fn get_q(&self) -> Q { + assert!(self.len() > 0); + let mut a = Array::zeros((self.dim(), self.len())); + for (i, mut col) in a.axis_iter_mut(Axis(1)).enumerate() { + col[i] = A::one(); + self.backward_reflection(&mut col); + } + a + } +} + +/// Online QR decomposition using Householder reflection +pub fn householder( + iter: impl Iterator>, + dim: usize, + rtol: A::Real, + strategy: Strategy, +) -> (Q, R) +where + A: Scalar + Lapack, + S: Data, +{ + let h = Householder::new(dim); + qr(iter, h, rtol, strategy) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::assert::*; + use num_traits::Zero; + + #[test] + fn check_reflector() { + let mut a = array![c64::new(1.0, 1.0), c64::new(1.0, 0.0), c64::new(0.0, 1.0)]; + let mut w = a.clone(); + calc_reflector(&mut w); + reflect(&w, &mut a); + close_l2( + &a, + &array![-c64::new(2.0.sqrt(), 2.0.sqrt()), c64::zero(), c64::zero()], + 1e-9, + ); + } +} diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs index 87e57c9b..9caff111 100644 --- a/src/krylov/mgs.rs +++ b/src/krylov/mgs.rs @@ -4,25 +4,6 @@ use super::*; use crate::{generate::*, inner::*, norm::Norm}; /// Iterative orthogonalizer using modified Gram-Schmit procedure -/// -/// ```rust -/// # use ndarray::*; -/// # 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); -/// -/// 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 @@ -31,7 +12,7 @@ pub struct MGS { q: Vec>, } -impl MGS { +impl MGS { /// Create an empty orthogonalizer pub fn new(dimension: usize) -> Self { Self { @@ -39,6 +20,28 @@ impl MGS { q: Vec::new(), } } + + /// Orthogonalize given vector against to the current basis + /// + /// - Returned array is coefficients and residual norm + /// - `a` will contain the residual vector + /// + pub fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 + where + S: DataMut, + { + assert_eq!(a.len(), self.dim()); + let mut coef = Array1::zeros(self.len() + 1); + for i in 0..self.len() { + let q = &self.q[i]; + let c = q.inner(&a); + 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); + coef + } } impl Orthogonalizer for MGS { @@ -52,22 +55,13 @@ impl Orthogonalizer for MGS { self.q.len() } - fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 + fn coeff(&self, a: ArrayBase) -> Array1 where A: Lapack, - S: DataMut, + S: Data, { - assert_eq!(a.len(), self.dim()); - let mut coef = Array1::zeros(self.len() + 1); - for i in 0..self.len() { - let q = &self.q[i]; - let c = q.inner(&a); - 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); - coef + let mut a = a.into_owned(); + self.orthogonalize(&mut a) } fn append(&mut self, a: ArrayBase, rtol: A::Real) -> Result, Array1> @@ -92,7 +86,7 @@ impl Orthogonalizer for MGS { } } -/// Online QR decomposition of vectors using modified Gram-Schmit algorithm +/// Online QR decomposition using modified Gram-Schmit algorithm pub fn mgs( iter: impl Iterator>, dim: usize, diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs index 677b1601..b53df0c7 100644 --- a/src/krylov/mod.rs +++ b/src/krylov/mod.rs @@ -3,8 +3,10 @@ use crate::types::*; use ndarray::*; -mod mgs; +pub mod householder; +pub mod mgs; +pub use householder::{householder, Householder}; pub use mgs::{mgs, MGS}; /// Q-matrix @@ -22,6 +24,28 @@ pub type Q = Array2; pub type R = Array2; /// Trait for creating orthogonal basis from iterator of arrays +/// +/// Example +/// ------- +/// +/// ```rust +/// # use ndarray::*; +/// # 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); +/// +/// 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); +/// } +/// ``` pub trait Orthogonalizer { type Elem: Scalar; @@ -40,15 +64,18 @@ pub trait Orthogonalizer { self.len() == 0 } - /// Orthogonalize given vector using current basis + /// Calculate the coefficient to the given basis and residual norm + /// + /// - The length of the returned array must be `self.len() + 1` + /// - Last component is the residual norm /// /// Panic /// ------- /// - if the size of the input array mismatches to the dimension /// - fn orthogonalize(&self, a: &mut ArrayBase) -> Array1 + fn coeff(&self, a: ArrayBase) -> Array1 where - S: DataMut; + S: Data; /// Add new vector if the residual is larger than relative tolerance /// diff --git a/tests/householder.rs b/tests/householder.rs new file mode 100644 index 00000000..adc4d9b2 --- /dev/null +++ b/tests/householder.rs @@ -0,0 +1,96 @@ +use ndarray::*; +use ndarray_linalg::{krylov::*, *}; + +fn over(rtol: A::Real) { + const N: usize = 4; + let a: Array2 = random((N, N * 2)); + + // Terminate + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let a_sub = a.slice(s![.., 0..N]); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I"); + assert_close_l2!(&q.dot(&r), &a_sub, rtol; "Check A = QR"); + + // Skip + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip); + let a_sub = a.slice(s![.., 0..N]); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + assert_close_l2!(&q.dot(&r), &a_sub, rtol); + + // Full + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Full); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol); + assert_close_l2!(&q.dot(&r), &a, rtol); +} + +#[test] +fn over_f32() { + over::(1e-5); +} +#[test] +fn over_f64() { + over::(1e-9); +} +#[test] +fn over_c32() { + over::(1e-5); +} +#[test] +fn over_c64() { + over::(1e-9); +} + +fn full(rtol: A::Real) { + const N: usize = 5; + let a: Array2 = random((N, N)); + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I"); + assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR"); +} + +#[test] +fn full_f32() { + full::(1e-5); +} +#[test] +fn full_f64() { + full::(1e-9); +} +#[test] +fn full_c32() { + full::(1e-5); +} +#[test] +fn full_c64() { + full::(1e-9); +} + +fn half(rtol: A::Real) { + const N: usize = 4; + let a: Array2 = random((N, N / 2)); + let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate); + let qc: Array2 = conjugate(&q); + assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol; "Check Q^H Q = I"); + assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR"); +} + +#[test] +fn half_f32() { + half::(1e-5); +} +#[test] +fn half_f64() { + half::(1e-9); +} +#[test] +fn half_c32() { + half::(1e-5); +} +#[test] +fn half_c64() { + half::(1e-9); +}