Implement Arnoldi method

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/inner.rs

@@ -0,0 +1,21 @@
+use crate::types::*;
+use ndarray::*;
+
+pub trait Inner<S: Data> {
+    /// Inner product `(self.conjugate, rhs)`
+    fn inner(&self, rhs: &ArrayBase<S, Ix1>) -> S::Elem;
+}
+
+impl<A, S1, S2> Inner<S1> for ArrayBase<S2, Ix1>
+where
+    A: Scalar,
+    S1: Data<Elem = A>,
+    S2: Data<Elem = A>,
+{
+    fn inner(&self, rhs: &ArrayBase<S1, Ix1>) -> A {
+        Zip::from(self)
+            .and(rhs)
+            .fold_while(A::zero(), |acc, s, r| FoldWhile::Continue(acc + s.conj() * *r))
+            .into_inner()
+    }
+}

src/krylov/householder.rs

@@ -0,0 +1,107 @@
+use super::*;
+use crate::{inner::*, norm::*};
+use num_traits::Zero;
+
+/// Iterative orthogonalizer using Householder reflection
+#[derive(Debug, Clone)]
+pub struct Householder<A: Scalar> {
+    /// Dimension of orthogonalizer
+    dim: usize,
+
+    /// Store Householder reflector.
+    ///
+    /// The coefficient is copied into another array, and this does not contain
+    v: Vec<Array1<A>>,
+}
+
+impl<A: Scalar> Householder<A> {
+    /// Create a new orthogonalizer
+    pub fn new(dim: usize) -> Self {
+        Householder { dim, v: Vec::new() }
+    }
+
+    /// 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, "Input array size mismaches to the dimension");
+
+        let w = self.v[k].slice(s![k..]);
+        let mut a_slice = a.slice_mut(s![k..]);
+        let c = A::from(2.0).unwrap() * w.inner(&a_slice);
+        for l in 0..self.dim - k {
+            a_slice[l] -= c * w[l];
+        }
+    }
+}
+
+impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
+    type Elem = A;
+
+    fn dim(&self) -> usize {
+        self.dim
+    }
+
+    fn len(&self) -> usize {
+        self.v.len()
+    }
+
+    fn orthogonalize<S>(&self, a: &mut ArrayBase<S, Ix1>) -> A::Real
+    where
+        S: DataMut<Elem = A>,
+    {
+        for k in 0..self.len() {
+            self.reflect(k, a);
+        }
+        if self.is_full() {
+            return Zero::zero();
+        }
+        // residual norm
+        a.slice(s![self.len()..]).norm_l2()
+    }
+
+    fn append<S>(&mut self, mut a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
+    where
+        S: DataMut<Elem = A>,
+    {
+        assert_eq!(a.len(), self.dim);
+        let k = self.len();
+        let alpha = self.orthogonalize(&mut a);
+        let mut coef = Array::zeros(k + 1);
+        for i in 0..k {
+            coef[i] = a[i];
+        }
+        if alpha < rtol {
+            // linearly dependent
+            coef[k] = A::from_real(alpha);
+            return Err(coef);
+        }
+
+        // Add reflector
+        assert!(k < a.len()); // this must hold because `alpha == 0` if k >= a.len()
+        let alpha = if a[k].abs() > Zero::zero() {
+            a[k].mul_real(alpha / a[k].abs())
+        } else {
+            A::from_real(alpha)
+        };
+        coef[k] = alpha;
+
+        a[k] -= alpha;
+        let norm = a.slice(s![k..]).norm_l2();
+        azip!(mut a (a.slice_mut(s![..k])) in { *a = Zero::zero() }); // this can be omitted
+        azip!(mut a (a.slice_mut(s![k..])) in { *a = a.div_real(norm) });
+        self.v.push(a.into_owned());
+        Ok(coef)
+    }
+
+    fn get_q(&self) -> Q<A> {
+        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();
+            for l in (0..self.len()).rev() {
+                self.reflect(l, &mut col);
+            }
+        }
+        a
+    }
+}

src/krylov/mgs.rs

@@ -0,0 +1,101 @@
+use super::*;
+use crate::{generate::*, inner::*, norm::Norm};
+
+/// Iterative orthogonalizer using modified Gram-Schmit procedure
+#[derive(Debug, Clone)]
+pub struct MGS<A> {
+    /// Dimension of base space
+    dim: usize,
+    /// Basis of spanned space
+    q: Vec<Array1<A>>,
+}
+
+impl<A: Scalar> MGS<A> {
+    pub fn new(dim: usize) -> Self {
+        Self { dim, q: Vec::new() }
+    }
+
+    fn ortho<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
+    where
+        A: Lapack,
+        S: DataMut<Elem = A>,
+    {
+        assert_eq!(a.len(), self.dim);
+        let mut coef = Array1::zeros(self.q.len() + 1);
+        for i in 0..self.q.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.q.len()] = A::from(nrm).unwrap();
+        coef
+    }
+}
+
+impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
+    type Elem = A;
+
+    fn dim(&self) -> usize {
+        self.dim
+    }
+
+    fn len(&self) -> usize {
+        self.q.len()
+    }
+
+    fn orthogonalize<S>(&self, a: &mut ArrayBase<S, Ix1>) -> A::Real
+    where
+        S: DataMut<Elem = A>,
+    {
+        let coef = self.ortho(a);
+        // Write coefficients into `a`
+        azip!(mut a (a.slice_mut(s![0..self.len()])), coef in { *a = coef });
+        // 0-fill for remaining
+        azip!(mut a (a.slice_mut(s![self.len()..])) in { *a = A::zero() });
+        coef[self.len()].re()
+    }
+
+    fn append<S>(&mut self, mut a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
+    where
+        S: DataMut<Elem = A>,
+    {
+        let coef = self.ortho(&mut a);
+        let nrm = coef[self.len()].re();
+        if nrm < rtol {
+            // Linearly dependent
+            return Err(coef);
+        }
+        azip!(mut a in { *a = *a / A::from_real(nrm) });
+        self.q.push(a.into_owned());
+        Ok(coef)
+    }
+
+    fn get_q(&self) -> Q<A> {
+        hstack(&self.q).unwrap()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::assert::*;
+
+    #[test]
+    fn mgs_append() {
+        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);
+
+        assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err());
+
+        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);
+        }
+    }
+
+}