array_matrix/matrix/
qr_householder.rs

1use num_complex::Complex;
2use num_traits::{Float, Zero, One};
3
4use crate::{Matrix, SquareMatrix, matrix_init, MMul, Transpose};
5
6pub trait QRHouseholder: Matrix
7{
8    type OutputQ;
9    type OutputR;
10
11    /// Returns the Householder QR-decomposition of the given matrix
12    /// 
13    /// A = QR
14    /// 
15    /// # Examples
16    /// 
17    /// ```rust
18    /// let a = [
19    ///     [1.0, 2.0],
20    ///     [3.0, 4.0]
21    /// ];
22    /// let (q, r) = a.qr_householder();
23    /// ```
24    fn qr_householder(&self) -> (Self::OutputQ, Self::OutputR);
25}
26
27impl<F: Float, const L: usize, const H: usize> QRHouseholder for [[Complex<F>; L]; H]
28where
29    Self: Matrix,
30    [[Complex<F>; H]; H]: SquareMatrix
31{
32    type OutputQ = [[Complex<F>; H]; H];
33    type OutputR = [[Complex<F>; L]; H];
34
35    fn qr_householder(&self) -> (Self::OutputQ, Self::OutputR)
36    {
37        assert!(H >= L);
38        let mut a: Vec<Vec<Complex<F>>> = (0..H)
39            .map(|r| (0..L)
40                .map(|c| self[r][c])
41                .collect()
42            ).collect();
43        let mut q: [[Complex<F>; H]; H] = SquareMatrix::identity();
44        for t in 0..L.min(H - 1)
45        {
46            let x: Vec<Complex<F>> = (0..H - t)
47                .map(|r| a[r][0])
48                .collect();
49            let x_abs = x.iter()
50                .map(|xn| xn.norm_sqr())
51                .reduce(|a, b| a + b)
52                .unwrap_or(F::zero())
53                .sqrt();
54            let alpha = -Complex::cis(x[t].arg())*x_abs;
55            let mut u = x;
56            u[0] = u[0] - alpha;
57            let u_abs = u.iter()
58                .map(|un| un.norm_sqr())
59                .reduce(|a, b| a + b)
60                .unwrap_or(F::zero())
61                .sqrt();
62            let v: Vec<Complex<F>> = u.iter()
63                .map(|un| un/u_abs)
64                .collect();
65            let q_: Vec<Vec<Complex<F>>> = (0..H - t)
66                .map(|r| (0..H - t)
67                    .map(|c| if r == c {Complex::<F>::one()} else {Complex::zero()} - v[r]*v[c].conj()*F::from(2.0).unwrap())
68                    .collect()
69                ).collect();
70            let qt: [[Complex<F>; H]; H] = matrix_init(|r, c| if r >= t && c >= t
71                {
72                    q_[c - t][r - t]
73                }
74                else
75                {
76                    if r == c {Complex::one()} else {Complex::zero()}
77                }
78            );
79            q = q.mul(qt);
80            a = (1..H - t)
81                .map(|r| (1..L - t)
82                    .map(|c| (0..H - t)
83                        .map(|i| q_[r][i]*a[i][c])
84                        .reduce(|a, b| a + b)
85                        .unwrap_or(Complex::zero())
86                    ).collect()
87                ).collect()
88        }
89        let r = q.transpose().mul(self.clone());
90        return (q, r)
91    }
92}
93
94impl<const L: usize, const H: usize> QRHouseholder for [[f32; L]; H]
95where
96    Self: Matrix,
97    [[Complex<f32>; L]; H]: QRHouseholder
98{
99    type OutputQ = <[[Complex<f32>; L]; H] as QRHouseholder>::OutputQ;
100    type OutputR = <[[Complex<f32>; L]; H] as QRHouseholder>::OutputR;
101    fn qr_householder(&self) -> (Self::OutputQ, Self::OutputR)
102    {
103        self.map(|ar| ar.map(|arc| Complex::from(arc))).qr_householder()
104    }
105}
106
107impl<const L: usize, const H: usize> QRHouseholder for [[f64; L]; H]
108where
109    Self: Matrix,
110    [[Complex<f64>; L]; H]: QRHouseholder
111{
112    type OutputQ = <[[Complex<f64>; L]; H] as QRHouseholder>::OutputQ;
113    type OutputR = <[[Complex<f64>; L]; H] as QRHouseholder>::OutputR;
114    fn qr_householder(&self) -> (Self::OutputQ, Self::OutputR)
115    {
116        self.map(|ar| ar.map(|arc| Complex::from(arc))).qr_householder()
117    }
118}
array_matrix/matrix/qr_householder.rs

array_matrix/matrix/
qr_householder.rs
