ndarray_linalg/krylov/
householder.rs

1//! Householder reflection
2//!
3//! - [Householder transformation - Wikipedia](https://en.wikipedia.org/wiki/Householder_transformation)
4//!
5
6use super::*;
7use crate::{inner::*, norm::*};
8use num_traits::One;
9
10/// Calc a reflactor `w` from a vector `x`
11pub fn calc_reflector<A>(x: &mut ArrayRef<A, Ix1>)
12where
13    A: Scalar + Lapack,
14{
15    assert!(!x.is_empty());
16    let norm = x.norm_l2();
17    let alpha = -x[0].mul_real(norm / x[0].abs());
18    x[0] -= alpha;
19    let inv_rev_norm = A::Real::one() / x.norm_l2();
20    azip!((a in x) *a = a.mul_real(inv_rev_norm));
21}
22
23/// Take a reflection `P = I - 2ww^T`
24///
25/// Panic
26/// ------
27/// - if the size of `w` and `a` mismaches
28pub fn reflect<A>(w: &ArrayRef<A, Ix1>, a: &mut ArrayRef<A, Ix1>)
29where
30    A: Scalar + Lapack,
31{
32    assert_eq!(w.len(), a.len());
33    let n = a.len();
34    let c = A::from(2.0).unwrap() * w.inner(a);
35    for l in 0..n {
36        a[l] -= c * w[l];
37    }
38}
39
40/// Iterative orthogonalizer using Householder reflection
41#[derive(Debug, Clone)]
42pub struct Householder<A: Scalar> {
43    /// Dimension of orthogonalizer
44    dim: usize,
45
46    /// Store Householder reflector.
47    ///
48    /// The coefficient is copied into another array, and this does not contain
49    v: Vec<Array1<A>>,
50
51    /// Tolerance
52    tol: A::Real,
53}
54
55impl<A: Scalar + Lapack> Householder<A> {
56    /// Create a new orthogonalizer
57    pub fn new(dim: usize, tol: A::Real) -> Self {
58        Householder {
59            dim,
60            v: Vec::new(),
61            tol,
62        }
63    }
64
65    /// Take a Reflection `P = I - 2ww^T`
66    fn fundamental_reflection(&self, k: usize, a: &mut ArrayRef<A, Ix1>) {
67        assert!(k < self.v.len());
68        assert_eq!(
69            a.len(),
70            self.dim,
71            "Input array size mismaches to the dimension"
72        );
73        reflect(&self.v[k].slice(s![k..]), &mut a.slice_mut(s![k..]));
74    }
75
76    /// Take forward reflection `P = P_l ... P_1`
77    pub fn forward_reflection(&self, a: &mut ArrayRef<A, Ix1>) {
78        assert!(a.len() == self.dim);
79        let l = self.v.len();
80        for k in 0..l {
81            self.fundamental_reflection(k, a);
82        }
83    }
84
85    /// Take backward reflection `P = P_1 ... P_l`
86    pub fn backward_reflection(&self, a: &mut ArrayRef<A, Ix1>) {
87        assert!(a.len() == self.dim);
88        let l = self.v.len();
89        for k in (0..l).rev() {
90            self.fundamental_reflection(k, a);
91        }
92    }
93
94    /// Compose coefficients array using reflected vector
95    fn compose_coefficients(&self, a: &ArrayRef<A, Ix1>) -> Coefficients<A> {
96        let k = self.len();
97        let res = a.slice(s![k..]).norm_l2();
98        let mut c = Array1::zeros(k + 1);
99        azip!((c in c.slice_mut(s![..k]), &a in a.slice(s![..k])) *c = a);
100        if k < a.len() {
101            let ak = a[k];
102            c[k] = -ak.mul_real(res / ak.abs());
103        } else {
104            c[k] = A::from_real(res);
105        }
106        c
107    }
108
109    /// Construct the residual vector from reflected vector
110    fn construct_residual(&self, a: &mut ArrayRef<A, Ix1>) {
111        let k = self.len();
112        azip!((a in a.slice_mut(s![..k])) *a = A::zero());
113        self.backward_reflection(a);
114    }
115}
116
117impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
118    type Elem = A;
119
120    fn dim(&self) -> usize {
121        self.dim
122    }
123
124    fn len(&self) -> usize {
125        self.v.len()
126    }
127
128    fn tolerance(&self) -> A::Real {
129        self.tol
130    }
131
132    fn decompose(&self, a: &mut ArrayRef<A, Ix1>) -> Array1<A> {
133        self.forward_reflection(a);
134        let coef = self.compose_coefficients(a);
135        self.construct_residual(a);
136        coef
137    }
138
139    fn coeff<S>(&self, a: ArrayBase<S, Ix1>) -> Array1<A>
140    where
141        S: Data<Elem = A>,
142    {
143        let mut a = a.into_owned();
144        self.forward_reflection(&mut a);
145        self.compose_coefficients(&a)
146    }
147
148    fn div_append(&mut self, a: &mut ArrayRef<A, Ix1>) -> AppendResult<A> {
149        assert_eq!(a.len(), self.dim);
150        let k = self.len();
151        self.forward_reflection(a);
152        let coef = self.compose_coefficients(a);
153        if coef[k].abs() < self.tol {
154            return AppendResult::Dependent(coef);
155        }
156        calc_reflector(&mut a.slice_mut(s![k..]));
157        self.v.push(a.to_owned());
158        self.construct_residual(a);
159        AppendResult::Added(coef)
160    }
161
162    fn append<S>(&mut self, a: ArrayBase<S, Ix1>) -> AppendResult<A>
163    where
164        S: Data<Elem = A>,
165    {
166        assert_eq!(a.len(), self.dim);
167        let mut a = a.into_owned();
168        let k = self.len();
169        self.forward_reflection(&mut a);
170        let coef = self.compose_coefficients(&a);
171        if coef[k].abs() < self.tol {
172            return AppendResult::Dependent(coef);
173        }
174        calc_reflector(&mut a.slice_mut(s![k..]));
175        self.v.push(a.to_owned());
176        AppendResult::Added(coef)
177    }
178
179    fn get_q(&self) -> Q<A> {
180        assert!(self.len() > 0);
181        let mut a = Array::zeros((self.dim(), self.len()));
182        for (i, mut col) in a.axis_iter_mut(Axis(1)).enumerate() {
183            col[i] = A::one();
184            self.backward_reflection(&mut col);
185        }
186        a
187    }
188}
189
190/// Online QR decomposition using Householder reflection
191pub fn householder<A, S>(
192    iter: impl Iterator<Item = ArrayBase<S, Ix1>>,
193    dim: usize,
194    rtol: A::Real,
195    strategy: Strategy,
196) -> (Q<A>, R<A>)
197where
198    A: Scalar + Lapack,
199    S: Data<Elem = A>,
200{
201    let h = Householder::new(dim, rtol);
202    qr(iter, h, strategy)
203}
204
205#[cfg(test)]
206mod tests {
207    use super::*;
208    use crate::assert::*;
209    use num_traits::Zero;
210
211    #[test]
212    fn check_reflector() {
213        let mut a = array![c64::new(1.0, 1.0), c64::new(1.0, 0.0), c64::new(0.0, 1.0)];
214        let mut w = a.clone();
215        calc_reflector(&mut w);
216        reflect(&w, &mut a);
217        close_l2(
218            &a,
219            &array![-c64::new(2.0.sqrt(), 2.0.sqrt()), c64::zero(), c64::zero()],
220            1e-9,
221        );
222    }
223}
ndarray_linalg/krylov/householder.rs

ndarray_linalg/krylov/
householder.rs
