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

ndarray_linalg/krylov/
householder.rs
