Householder reflection

Author: termoshtt

src/arnoldi.rs

@@ -65,10 +65,6 @@ impl<A: Scalar> MGS<A> {
 
     /// Add new vector if the residual is larger than relative tolerance
     ///
-    /// Panic
-    /// -------
-    /// - if the size of the input array mismaches to the dimension
-    ///
     /// ```rust
     /// # use ndarray::*;
     /// # use ndarray_linalg::*;
@@ -85,6 +81,11 @@ impl<A: Scalar> MGS<A> {
     ///     close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); // You can get coefficients of dependent vector
     /// }
     /// ```
+    ///
+    /// Panic
+    /// -------
+    /// - if the size of the input array mismaches to the dimension
+    ///
     pub fn append<S>(&mut self, a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
     where
         A: Lapack,

src/householder.rs

@@ -0,0 +1,84 @@
+use crate::{inner::*, norm::*, types::*};
+use ndarray::*;
+
+/// Iterative orthogonalizer using Householder reflection
+#[derive(Debug, Clone)]
+pub struct Householder<A> {
+    dim: usize,
+    v: Vec<Array1<A>>,
+}
+
+impl<A: Scalar> Householder<A> {
+    pub fn new(dim: usize) -> Self {
+        Householder { dim, v: Vec::new() }
+    }
+
+    pub fn len(&self) -> usize {
+        self.v.len()
+    }
+
+    /// Take a Reflection `P = I - 2ww^T`
+    fn reflect<S: DataMut<Elem = A>>(&self, k: usize, a: &mut ArrayBase<S, Ix1>) {
+        assert!(k < self.v.len());
+        assert_eq!(a.len(), self.dim);
+        let w = self.v[k].slice(s![k..]);
+        let c = A::from(2.0).unwrap() * w.inner(&a.slice(s![k..]));
+        for l in k..self.dim {
+            a[l] -= c * w[l];
+        }
+    }
+
+    /// Orghotonalize input vector by reflectors
+    ///
+    /// Panic
+    /// -------
+    /// - if the size of the input array mismaches to the dimension
+    pub fn orthogonalize<S>(&self, a: &mut ArrayBase<S, Ix1>) -> A::Real
+    where
+        A: Lapack,
+        S: DataMut<Elem = A>,
+    {
+        for k in 0..self.len() {
+            self.reflect(k, a);
+        }
+        // residual norm
+        a.slice(s![self.len()..]).norm_l2()
+    }
+
+    /// Orghotonalize input vector by reflectors
+    ///
+    /// ```rust
+    /// # use ndarray::*;
+    /// # use ndarray_linalg::*;
+    /// let mut mgs = arnoldi::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();
+    ///
+    /// 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_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
+    /// }
+    /// ```
+    ///
+    /// Panic
+    /// -------
+    /// - if the size of the input array mismaches to the dimension
+    ///
+    pub fn append<S>(&mut self, mut a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
+    where
+        A: Lapack,
+        S: DataMut<Elem = A>,
+    {
+        let residual = self.orthogonalize(&mut a);
+        let coef = a.slice(s![..self.len()]).into_owned();
+        if residual < rtol {
+            return Err(coef);
+        }
+        self.v.push(a.into_owned());
+        Ok(coef)
+    }
+}

src/lib.rs

@@ -25,6 +25,7 @@ pub mod diagonal;
 pub mod eigh;
 pub mod error;
 pub mod generate;
+pub mod householder;
 pub mod inner;
 pub mod lapack;
 pub mod layout;