feanor_math/algorithms/
qr.rs

1use crate::algorithms::matmul::ComputeInnerProduct;
2use crate::field::{Field, FieldStore};
3use crate::integer::*;
4use crate::matrix::*;
5use crate::ring::*;
6use crate::rings::approx_real::{ApproxRealField, SqrtRing};
7use crate::rings::rational::*;
8
9use std::cmp::min;
10
11#[stability::unstable(feature = "enable")]
12pub trait QRDecompositionField: Field {
13
14    ///
15    /// Given a matrix `A`, computes an orthogonal matrix `Q` and an upper triangular
16    /// matrix `R` with `A = Q R`. The function writes `Q diag(x_1, ..., x_n)` to `q` and
17    /// `diag(1/x_1, ..., 1/x_n) R` to `matrix`, and returns `x_1^2, ..., x_n^2`, where
18    /// `x_1, ..., x_n` are the elements on the diagonal of `R`.
19    /// 
20    /// Returning the values as given above instead of just `Q` and `R` is done
21    /// to avoid the computation of square-roots, which may not be supported by the
22    /// underlying ring. If it is supported, you can use [`QRDecompositionField::qr_decomposition()`]
23    /// instead. Note that this means that `diag(x_1^2, ..., x_n^2)` and `R`
24    /// are the LDL-decomposition of `A^T A`.
25    /// 
26    /// # Rank-deficient matrices
27    /// 
28    /// Do not use this for matrices that do not have full rank. If the underlying ring
29    /// is exact, this will panic. For approximate rings (in particular floating-point numbers),
30    /// matrices that don't have full rank, or are very badly conditioned, will give inaccurate
31    /// results.
32    /// 
33    /// Clearly, rank-deficient matrices cannot be supported, since for those the value
34    /// `diag(1/x_1, ..., 1/x_n)` is not defined.
35    /// 
36    fn scaled_qr_decomposition<V1, V2>(&self, matrix: SubmatrixMut<V1, Self::Element>, q: SubmatrixMut<V2, Self::Element>) -> Vec<Self::Element>
37        where V1: AsPointerToSlice<Self::Element>, V2: AsPointerToSlice<Self::Element>;
38
39    ///
40    /// Given a square symmetric matrix `A`, computes a strict lower triangular matrix `L` and
41    /// a diagonal matrix `D` such that `A = L D L^T`. The function writes `L` to `matrix`
42    /// and returns the diagonal elements of `D`.
43    /// 
44    /// # Singular matrices
45    /// 
46    /// Do not use this for matrices that are singular. If the underlying ring is exact, 
47    /// this will panic. For approximate rings (in particular floating-point numbers),
48    /// matrices that don't have full rank, or are very badly conditioned, will give inaccurate
49    /// results. Note however that the matrix is not required to be positive definite, it may
50    /// have both positive and negative eigenvalues (but no zero eigenvalues).
51    /// 
52    /// Why don't we support singular matrices? Because many singular matrices don't have
53    /// an LDL decomposition. For example, the matrix `[[ 0, 1 ], [ 1, 1 ]]` doesn't.
54    /// 
55    fn ldl_decomposition<V>(&self, matrix: SubmatrixMut<V, Self::Element>) -> Vec<Self::Element>
56        where V: AsPointerToSlice<Self::Element>
57    {
58        ldl_decomposition_impl(RingRef::new(self), matrix)
59    }
60       
61    ///
62    /// Given a matrix `A`, computes an orthogonal matrix `Q` and an upper triangular
63    /// matrix `R` with `A = Q R`. These are returned in `matrix` and `q`, respectively.
64    /// 
65    /// Note that if the ring is not a [`SqrtRing`], you can still use [`QRDecompositionField::scaled_qr_decomposition()`].
66    /// 
67    /// This function supports non-full-rank matrices as well.
68    /// 
69    fn qr_decomposition<V1, V2>(&self, mut matrix: SubmatrixMut<V1, Self::Element>, mut q: SubmatrixMut<V2, Self::Element>)
70        where V1: AsPointerToSlice<Self::Element>, V2: AsPointerToSlice<Self::Element>, Self: SqrtRing
71    {
72        let d = self.scaled_qr_decomposition(matrix.reborrow(), q.reborrow());
73        for (i, scale_sqr) in d.into_iter().enumerate() {
74            let scale = self.sqrt(scale_sqr);
75            let scale_inv = self.div(&self.one(), &scale);
76            for j in 0..matrix.col_count() {
77                self.mul_assign_ref(matrix.at_mut(i, j), &scale);
78            }
79            for k in 0..q.row_count() {
80                self.mul_assign_ref(q.at_mut(k, i), &scale_inv);
81            }
82        }
83    }
84}
85
86impl<I> QRDecompositionField for RationalFieldBase<I>
87    where I: RingStore,
88        I::Type: IntegerRing
89{
90    fn scaled_qr_decomposition<V1, V2>(&self, mut matrix: SubmatrixMut<V1, Self::Element>, mut q: SubmatrixMut<V2, Self::Element>) -> Vec<Self::Element>
91        where V1: AsPointerToSlice<Self::Element>, V2: AsPointerToSlice<Self::Element>
92    {
93        // since there is no issue with numerical stability, we can do Gram-Schmidt
94        let ring = RingValue::from_ref(self);
95        let m = matrix.row_count();
96        let n = matrix.col_count();
97        assert_eq!(m, q.row_count());
98        assert_eq!(m, q.col_count());
99
100        let mut result = Vec::with_capacity(n);
101        let mut mus = Vec::with_capacity(n);
102        for i in 0..n {
103            mus.clear();
104            for j in 0..i {
105                mus.push(self.div(
106                    &<_ as ComputeInnerProduct>::inner_product_ref(self, (0..m).map(|k| (matrix.at(k, i), q.at(k, j)))),
107                    &result[j]
108                ));
109            }
110            let (mut target, orthogonalized) = q.reborrow().split_cols(i..(i + 1), 0..i);
111            for k in 0..m {
112                *target.at_mut(k, 0) = self.sub_ref_fst(
113                    matrix.at(k, i),
114                    <_ as ComputeInnerProduct>::inner_product_ref(self, (0..i).map(|j| (&mus[j], orthogonalized.at(k, j))))
115                );
116            }
117            result.push(<_ as RingStore>::sum(ring, (0..m).map(|k| ring.pow(ring.clone_el(target.at(k, 0)), 2))));
118            for (k, c) in mus.drain(..).enumerate() {
119                *matrix.at_mut(k, i) = c;
120            }
121            *matrix.at_mut(i, i) = self.one();
122            for k in (i + 1)..m {
123                *matrix.at_mut(k, i) = self.zero();
124            }
125        }
126
127        return result;
128    }
129}
130
131fn ldl_decomposition_impl<R, V>(ring: R, mut matrix: SubmatrixMut<V, El<R>>) -> Vec<El<R>>
132    where R: RingStore, 
133        R::Type: Field,
134        V: AsPointerToSlice<El<R>>
135{
136    assert_eq!(matrix.row_count(), matrix.col_count());
137    let n = matrix.row_count();
138    let mut result = Vec::with_capacity(n);
139    for i in 0..n {
140        let pivot = ring.clone_el(matrix.at(i, i));
141        if !ring.get_ring().is_approximate() && ring.is_zero(&pivot) {
142            panic!("matrix is singular")
143        }
144        let pivot_inv = ring.div(&ring.one(), matrix.at(i, i));
145        for j in i..n {
146            ring.mul_assign_ref(matrix.at_mut(j, i), &pivot_inv);
147        }
148        for k in (i + 1)..n {
149            for l in k..n {
150                let subtract = ring.mul_ref_snd(ring.mul_ref(matrix.as_const().at(k, i), matrix.as_const().at(l, i)), &pivot);
151                ring.sub_assign(matrix.at_mut(l, k), subtract);
152            }
153        }
154        result.push(pivot);
155    }
156    for i in 0..n {
157        for j in (i + 1)..n {
158            *matrix.at_mut(i, j) = ring.zero();
159        }
160    }
161    return result;
162}
163
164impl<R: ApproxRealField + SqrtRing> QRDecompositionField for R {
165
166    default fn scaled_qr_decomposition<V1, V2>(&self, mut matrix: SubmatrixMut<V1, Self::Element>, mut q: SubmatrixMut<V2, Self::Element>) -> Vec<Self::Element>
167        where V1: AsPointerToSlice<Self::Element>, V2: AsPointerToSlice<Self::Element>
168    {
169        self.qr_decomposition(matrix.reborrow(), q.reborrow());
170        let mut result = Vec::with_capacity(matrix.row_count());
171        for i in 0..matrix.row_count() {
172            let mut scale = self.clone_el(matrix.at(i, i));
173            let scale_inv = self.div(&self.one(), &scale);
174            for j in i..matrix.col_count() {
175                self.mul_assign_ref(matrix.at_mut(i, j), &scale_inv);
176            }
177            for j in 0..q.row_count() {
178                self.mul_assign_ref(q.at_mut(j, i), &scale);
179            }
180            self.square(&mut scale);
181            result.push(scale);
182        }
183        return result;
184    }
185
186    default fn ldl_decomposition<V>(&self, matrix: SubmatrixMut<V, Self::Element>) -> Vec<Self::Element>
187        where V: AsPointerToSlice<Self::Element>
188    {
189        ldl_decomposition_impl(RingRef::new(self), matrix)
190    }
191
192    default fn qr_decomposition<V1, V2>(&self, mut matrix: SubmatrixMut<V1, Self::Element>, mut q: SubmatrixMut<V2, Self::Element>)
193        where V1: AsPointerToSlice<Self::Element>, V2: AsPointerToSlice<Self::Element>
194    {
195        let ring = RingRef::new(self);
196        let m = matrix.row_count();
197        let n = matrix.col_count();
198        assert_eq!(m, q.row_count());
199        assert_eq!(m, q.col_count());
200        for i in 0..m {
201            for j in 0..m {
202                *q.at_mut(i, j) = if i == j { self.one() } else { self.zero() };
203            }
204        }
205
206        let mut householder_vector = Vec::with_capacity(m);
207        for i in 0..min(n, m) {
208            let norm_sqr = <_ as RingStore>::sum(&ring, (i..m).map(|k| ring.pow(ring.clone_el(matrix.at(k, i)), 2)));
209            let norm = self.sqrt(self.clone_el(&norm_sqr));
210            let alpha = if self.is_neg(matrix.at(i, i)) {
211                self.clone_el(&norm)
212            } else {
213                self.negate(self.clone_el(&norm))
214            };
215            // | x - alpha * e1 | / sqrt(2)
216            let scale = self.sqrt(self.sub(norm_sqr, self.mul_ref(&alpha, matrix.at(i, i))));
217            householder_vector.clear();
218            householder_vector.extend((i..m).map(|k| ring.clone_el(matrix.at(k, i))));
219            ring.sub_assign_ref(&mut householder_vector[0], &alpha);
220            for x in &mut householder_vector {
221                *x = self.div(x, &scale);
222            }
223
224            // update matrix
225            let mut rest = matrix.reborrow().submatrix(i..m, (i + 1)..n);
226            for j in 0..(n - i - 1) {
227                let inner_product = <_ as ComputeInnerProduct>::inner_product_ref(self, (0..(m - i)).map(|k| (&householder_vector[k], rest.at(k, j))));
228                for k in 0..(m - i) {
229                    ring.sub_assign(rest.at_mut(k, j), ring.mul_ref(&inner_product, &householder_vector[k]));
230                }
231            }
232
233            // update q
234            let mut rest = q.reborrow().restrict_cols(i..m);
235            for j in 0..m {
236                let inner_product = <_ as ComputeInnerProduct>::inner_product_ref(self, (0..(m - i)).map(|k| (&householder_vector[k], rest.at(j, k))));
237                for k in 0..(m - i) {
238                    ring.sub_assign(rest.at_mut(j, k), ring.mul_ref(&inner_product, &householder_vector[k]));
239                }
240            }
241
242            // update pivot
243            let mut pivot_col = matrix.reborrow().submatrix(i..m, i..(i + 1));
244            for k in 1..(m - i) {
245                *pivot_col.at_mut(k, 0) = self.zero();
246            }
247            *pivot_col.at_mut(0, 0) = alpha;
248        }
249    }
250}
251
252#[cfg(test)]
253use crate::algorithms::matmul::STANDARD_MATMUL;
254#[cfg(test)]
255use crate::matrix::{TransposableSubmatrix, TransposableSubmatrixMut};
256#[cfg(test)]
257use crate::algorithms::matmul::MatmulAlgorithm;
258#[cfg(test)]
259use crate::matrix::format_matrix;
260#[cfg(test)]
261use crate::rings::approx_real::float::Real64;
262#[cfg(test)]
263use crate::homomorphism::Homomorphism;
264#[cfg(test)]
265use crate::assert_matrix_eq;
266#[cfg(test)]
267use crate::rings::fraction::FractionFieldStore;
268#[cfg(test)]
269use crate::primitive_int::StaticRing;
270
271#[cfg(test)]
272fn assert_is_correct_qr<V1, V2, V3>(original: Submatrix<V1, f64>, q: Submatrix<V2, f64>, r: Submatrix<V3, f64>)
273    where V1: AsPointerToSlice<f64>, V2: AsPointerToSlice<f64>, V3: AsPointerToSlice<f64>
274{
275    let m = q.row_count();
276    let n = r.col_count();
277    assert_eq!(m, original.row_count());
278    assert_eq!(n, original.col_count());
279    assert_eq!(m, r.row_count());
280    let mut product = OwnedMatrix::zero(m, n, Real64::RING);
281    STANDARD_MATMUL.matmul(
282        TransposableSubmatrix::from(q),
283        TransposableSubmatrix::from(r),
284        TransposableSubmatrixMut::from(product.data_mut()),
285        Real64::RING
286    );
287    for i in 0..m {
288        for j in 0..n {
289            if !(Real64::RING.get_ring().is_approx_eq(*original.at(i, j), *product.at(i, j), 100)) {
290                println!("product does not match; Q, R are");
291                println!("{}", format_matrix(m, m, |i, j| q.at(i, j), Real64::RING));
292                println!("and");
293                println!("{}", format_matrix(m, n, |i, j| r.at(i, j), Real64::RING));
294                println!("the product is");
295                println!("{}", format_matrix(m, n, |i, j| product.at(i, j), Real64::RING));
296                panic!();
297            }
298        }
299    }
300    let mut product = OwnedMatrix::zero(m, m, Real64::RING);
301    STANDARD_MATMUL.matmul(
302        TransposableSubmatrix::from(q).transpose(),
303        TransposableSubmatrix::from(q),
304        TransposableSubmatrixMut::from(product.data_mut()),
305        Real64::RING
306    );
307    for i in 0..m {
308        for j in 0..m {
309            let expected = if i == j { 1. } else { 0. };
310            if !(Real64::RING.get_ring().is_approx_eq(expected, *product.at(i, j), 100)) {
311                println!("Q is not orthogonal");
312                println!("{}", format_matrix(m, m, |i, j| q.at(i, j), Real64::RING));
313                panic!();
314            }
315        }
316    }
317
318    for j in 0..n {
319        for i in (j + 1)..m {
320            if !(Real64::RING.get_ring().is_approx_eq(0., *r.at(i, j), 100)) {
321                println!("R is not upper triangular");
322                println!("{}", format_matrix(m, n, |i, j| r.at(i, j), Real64::RING));
323                panic!();
324            }
325        }
326    }
327}
328
329#[cfg(test)]
330fn assert_is_correct_ldl<V1, V2>(original: Submatrix<V1, f64>, l: Submatrix<V2, f64>, d: &[f64])
331    where V1: AsPointerToSlice<f64>, V2: AsPointerToSlice<f64>
332{
333    let n = l.col_count();
334    assert_eq!(n, l.row_count());
335    assert_eq!(n, original.col_count());
336    assert_eq!(n, original.row_count());
337    let l_scaled = OwnedMatrix::from_fn(n, n, |i, j| *l.at(i, j) * d[j]);
338    let mut product = OwnedMatrix::zero(n, n, Real64::RING);
339    STANDARD_MATMUL.matmul(
340        TransposableSubmatrix::from(l_scaled.data()),
341        TransposableSubmatrix::from(l).transpose(),
342        TransposableSubmatrixMut::from(product.data_mut()),
343        Real64::RING
344    );
345    for i in 0..n {
346        for j in 0..n {
347            if !(Real64::RING.get_ring().is_approx_eq(*original.at(i, j), *product.at(i, j), 100)) {
348                println!("product does not match; L is");
349                println!("{}", format_matrix(n, n, |i, j| l.at(i, j), Real64::RING));
350                println!("D is diag{:?} and the product LDL^T is", d);
351                println!("{}", format_matrix(n, n, |i, j| product.at(i, j), Real64::RING));
352                panic!();
353            }
354        }
355    }
356    for i in 0..n {
357        for j in (i + 1)..n {
358            if !(Real64::RING.get_ring().is_approx_eq(0., *l.at(i, j), 100)) {
359                println!("L is not lower triangular");
360                println!("{}", format_matrix(n, n, |i, j| l.at(i, j), Real64::RING));
361                panic!();
362            }
363        }
364    }
365}
366
367#[test]
368fn test_float_qr() {
369    let RR = Real64::RING;
370    let a = OwnedMatrix::new_with_shape(vec![0., 1., 1., 0.], 2, 2);
371    let mut r = a.clone_matrix(RR);
372    let mut q = OwnedMatrix::zero(2, 2, RR);
373    RR.get_ring().qr_decomposition(r.data_mut(), q.data_mut());
374    assert_is_correct_qr(a.data(), q.data(), r.data());
375
376    let a = OwnedMatrix::new_with_shape(vec![1., 2., 3., 4., 5., 6.], 3, 2);
377    let mut r = a.clone_matrix(RR);
378    let mut q = OwnedMatrix::zero(3, 3, RR);
379    RR.get_ring().qr_decomposition(r.data_mut(), q.data_mut());
380    assert_is_correct_qr(a.data(), q.data(), r.data());
381
382    let a = OwnedMatrix::new_with_shape(vec![1., 2., 3., 4., 5., 6.], 2, 3);
383    let mut r = a.clone_matrix(RR);
384    let mut q = OwnedMatrix::zero(2, 2, RR);
385    RR.get_ring().qr_decomposition(r.data_mut(), q.data_mut());
386    assert_is_correct_qr(a.data(), q.data(), r.data());
387
388    let a = OwnedMatrix::new_with_shape(vec![1., 1., 1., 2., 2., 3., 0., 0., 1.], 3, 3);
389    let mut r = a.clone_matrix(RR);
390    let mut q = OwnedMatrix::zero(3, 3, RR);
391    RR.get_ring().qr_decomposition(r.data_mut(), q.data_mut());
392    assert_is_correct_qr(a.data(), q.data(), r.data());
393
394    let a = OwnedMatrix::new_with_shape((1..31).map(|x| x as f64 * if x % 2 == 0 { -1.0 } else { 1.0 }).collect::<Vec<_>>(), 6, 5);
395    let mut r = a.clone_matrix(RR);
396    let mut q = OwnedMatrix::zero(6, 6, RR);
397    RR.get_ring().qr_decomposition(r.data_mut(), q.data_mut());
398    assert_is_correct_qr(a.data(), q.data(), r.data());
399}
400
401#[test]
402fn test_float_qdr() {
403    let RR = Real64::RING;
404    let a = OwnedMatrix::new_with_shape((1..10).map(|c| c as f64).collect(), 3, 3);
405    let mut r = a.clone_matrix(RR);
406    let mut q = OwnedMatrix::zero(3, 3, RR);
407    let diags = RR.get_ring().scaled_qr_decomposition(r.data_mut(), q.data_mut());
408    for i in 0..3 {
409        for j in 0..3 {
410            if i == j {
411                assert!(RR.get_ring().is_approx_eq(1., *r.at(i, j), 100));
412            }
413            RR.mul_assign(r.at_mut(i, j), diags[i].sqrt());
414            RR.mul_assign(q.at_mut(i, j), 1. / diags[j].sqrt());
415        }
416    }
417    assert_is_correct_qr(a.data(), q.data(), r.data());
418}
419
420#[test]
421fn test_float_ldl() {
422    let RR = Real64::RING;
423    let a = OwnedMatrix::new_with_shape(vec![5., 1., 1., 5.], 2, 2);
424    let mut l = a.clone_matrix(RR);
425    let d = RR.get_ring().ldl_decomposition(l.data_mut());
426    assert_is_correct_ldl(a.data(), l.data(), &d);
427
428    let a = OwnedMatrix::new_with_shape(vec![1., 2., 3., 2., 6., 5., 3., 5., 20.], 3, 3);
429    let mut l = a.clone_matrix(RR);
430    let d = RR.get_ring().ldl_decomposition(l.data_mut());
431    assert_is_correct_ldl(a.data(), l.data(), &d);
432    
433    let mut a = OwnedMatrix::zero(5, 5, RR);
434    let factor = OwnedMatrix::new((0..25).map(|c| (c as f64).powi(2)).collect(), 5);
435    STANDARD_MATMUL.matmul(
436        TransposableSubmatrix::from(factor.data()),
437        TransposableSubmatrix::from(factor.data()).transpose(),
438        TransposableSubmatrixMut::from(a.data_mut()),
439        RR
440    );
441    let mut l = a.clone_matrix(RR);
442    let d = RR.get_ring().ldl_decomposition(l.data_mut());
443    assert_is_correct_ldl(a.data(), l.data(), &d);
444
445    let a = OwnedMatrix::new_with_shape(vec![1., 2., 3., 2., 6., 5., 3., 5., -20.], 3, 3);
446    let mut l = a.clone_matrix(RR);
447    let d = RR.get_ring().ldl_decomposition(l.data_mut());
448    assert_is_correct_ldl(a.data(), l.data(), &d);
449}
450
451#[test]
452fn test_rational_qdr() {
453    let QQ = RationalField::new(StaticRing::<i64>::RING);
454    let mut actual_r = OwnedMatrix::new_with_shape((1..10).map(|x| QQ.pow(QQ.int_hom().map(x), 2)).collect(), 3, 3);
455    let mut actual_q = OwnedMatrix::zero(3, 3, &QQ);
456    let diags = QQ.get_ring().scaled_qr_decomposition(actual_r.data_mut(), actual_q.data_mut());
457    assert_el_eq!(&QQ, QQ.from_fraction(2658, 1), &diags[0]);
458    assert_el_eq!(&QQ, QQ.from_fraction(9891, 443), &diags[1]);
459    assert_el_eq!(&QQ, QQ.from_fraction(864, 1099), &diags[2]);
460
461    let mut expected_r = OwnedMatrix::identity(3, 3, &QQ);
462    *expected_r.at_mut(0, 1) = QQ.from_fraction(590, 443);
463    *expected_r.at_mut(0, 2) = QQ.from_fraction(759, 443);
464    *expected_r.at_mut(1, 2) = QQ.from_fraction(2700, 1099);
465    assert_matrix_eq!(&QQ, expected_r, actual_r);
466
467    let expected_q_num = [[486857, 1299018, 356172], [7789712, 1796865, -233904], [23855993, -613242, 69108]];
468    let expected_q_den = 443 * 1099;
469    let expected_q = OwnedMatrix::from_fn(3, 3, |i, j| QQ.from_fraction(expected_q_num[i][j], expected_q_den));
470    assert_matrix_eq!(&QQ, expected_q, actual_q);
471}
feanor_math/algorithms/qr.rs

feanor_math/algorithms/
qr.rs
