qfall_math/rational/mat_q/
cholesky_decomp.rs

1// Copyright © 2024 Marvin Beckmann
2//
3// This file is part of qFALL-math.
4//
5// qFALL-math is free software: you can redistribute it and/or modify it under
6// the terms of the Mozilla Public License Version 2.0 as published by the
7// Mozilla Foundation. See <https://mozilla.org/en-US/MPL/2.0/>.
8
9//! This module implements the Cholesky Decomposition for Hermitian positive-definite
10//! matrices.
11
12use super::MatQ;
13use crate::{
14    rational::Q,
15    traits::{
16        Concatenate, MatrixDimensions, MatrixGetEntry, MatrixGetSubmatrix, MatrixSetEntry,
17        MatrixSetSubmatrix,
18    },
19};
20
21impl MatQ {
22    /// This function performs the Cholesky decomposition (the Cholesky algorithm) and
23    /// produces a matrix `L` such that `self = L * L^T`.
24    /// This function relies on the precision of [`Q::sqrt`] and will not provide
25    /// perfect results in all cases.
26    /// Furthermore, the Cholesky decomposition requires a Hermitian positive-definite
27    /// matrix.
28    ///
29    /// Returns the Cholesky decomposition of a Hermitian positive-definite matrix.
30    ///
31    /// # Examples
32    /// ```
33    /// use qfall_math::rational::MatQ;
34    /// use std::str::FromStr;
35    ///
36    /// let matrix = MatQ::from_str("[[4, 12, -16],[12,37,-43],[-16,-43,98]]").unwrap();
37    ///
38    /// let l = matrix.cholesky_decomposition();
39    /// assert_eq!(matrix, &l * l.transpose());
40    /// ```
41    ///
42    /// # Panics ...
43    /// - if `self` is not a symmetric matrix,
44    /// - if `self` has eigenvalues smaller than `0`.
45    pub fn cholesky_decomposition(&self) -> MatQ {
46        assert!(self.is_symmetric(), "The provided matrix is not symmetric.");
47        let n = self.get_num_columns();
48
49        let mut a = self.clone();
50        let mut l = MatQ::identity(n, n);
51
52        for i in 0..n {
53            // get first entry and ith column
54            let a_ii = unsafe { a.get_entry_unchecked(0, 0) };
55            assert!(a_ii > Q::ZERO, "The matrix is not positive-definite.");
56            let column_a_i = match i {
57                0 => unsafe { a.get_column_unchecked(0) },
58                _ => unsafe { MatQ::new(i, 1).get_column_unchecked(0) }
59                    .concat_vertical(&unsafe { a.get_column_unchecked(0) })
60                    .unwrap(),
61            } * (1 / (a_ii.sqrt()));
62            // in the previous line: sqrt panics if `a_ii` is negative, i.e. if an
63            // eigenvalue is negative.
64
65            // produce L matrix recursively
66            let mut l_i = MatQ::identity(n, n);
67            unsafe { l_i.set_column_unchecked(i, &column_a_i, 0) };
68            l = l * l_i;
69
70            // update matrix A recursively
71            if i < n - 1 {
72                let b = a.get_submatrix(1, -1, 1, -1).unwrap();
73                let b_minus = (1 / a_ii)
74                    * a.get_submatrix(1, -1, 0, 0).unwrap()
75                    * a.get_submatrix(0, 0, 1, -1).unwrap();
76                a = b - b_minus;
77            }
78        }
79        l
80    }
81
82    /// This function implements the Cholesky decomposition according to FLINTs
83    /// implementation. As FLINTs algorithm is not (yet) accessible through flint-sys,
84    /// this implementation follows the implementation of the algorithm from FLINT.
85    /// This, however, also means that we will work with less precision as we will work
86    /// with conversions to [`f64`] and not use [`Q`].
87    /// In turn, this makes the function much more efficient, but *not* applicable to
88    /// large numbers.
89    ///
90    /// This function relies on the precision of [`f64::sqrt`] and will not provide
91    /// perfect results in all cases.
92    /// Furthermore, the Cholesky decomposition requires a Hermitian positive-definite
93    /// matrix.
94    ///
95    /// Returns the Cholesky decomposition of a Hermitian positive-definite matrix.
96    ///
97    /// # Examples
98    /// ```
99    /// use qfall_math::rational::MatQ;
100    /// use std::str::FromStr;
101    ///
102    /// let matrix = MatQ::from_str("[[4, 12, -16],[12,37,-43],[-16,-43,98]]").unwrap();
103    ///
104    /// let l = matrix.cholesky_decomposition_flint();
105    /// assert_eq!(matrix, &l * l.transpose());
106    /// ```
107    ///
108    /// # Panics ...
109    /// - if `self` is not a symmetric matrix,
110    /// - if `self` has eigenvalues smaller than `0`.
111    #[allow(clippy::needless_range_loop)]
112    pub fn cholesky_decomposition_flint(&self) -> MatQ {
113        assert!(self.is_symmetric(), "The provided matrix is not symmetric.");
114        let mat_dimension = self.get_num_rows() as usize;
115
116        let mut out = vec![vec![0.0; mat_dimension]; mat_dimension];
117        let mat = self.collect_entries_f64();
118
119        // This code snippet originates from [flint](https://github.com/flintlib/flint/blob/main/src/fmpz_mat/chol_d.c)
120        // it is not part of [flint-sys] as it requires a specific data-type `d_mat_t`
121        for i in 0..mat_dimension {
122            for j in 0..=i {
123                let mut s: f64 = 0.0;
124                for k in 0..j {
125                    s += out[i][k] * out[j][k]
126                }
127                if i == j {
128                    // Find this requirement in https://en.wikipedia.org/wiki/Cholesky_decomposition#The_Cholesky_algorithm
129                    // a_ii > 0 as `self` needs to be positive definite
130                    assert!(
131                        mat[i][i] > s,
132                        "The provided matrix is not positive definite."
133                    );
134
135                    out[i][j] = (mat[i][i] - s).sqrt();
136                } else {
137                    out[i][j] = (mat[i][j] - s) / out[j][j];
138                }
139            }
140        }
141
142        // Assemble Cholesky decomposition as MatQ
143        let mut res = MatQ::new(mat_dimension, mat_dimension);
144        for (i, row) in out.iter().enumerate().take(mat_dimension) {
145            for (j, entry) in row.iter().enumerate().take(mat_dimension) {
146                unsafe { res.set_entry_unchecked(i as i64, j as i64, *entry) };
147            }
148        }
149
150        res
151    }
152}
153
154#[cfg(test)]
155mod test_cholesky_decomposition {
156    use crate::{
157        rational::{MatQ, Q},
158        traits::MatrixSetEntry,
159    };
160    use std::str::FromStr;
161
162    /// Ensure that a basic example (from Wikipedia) works.
163    #[test]
164    fn valid_input() {
165        let matrix = MatQ::from_str("[[4, 12, -16],[12,37,-43],[-16,-43,98]]").unwrap();
166
167        let l = MatQ::from_str("[[2, 0, 0],[6, 1, 0],[-8, 5, 3]]").unwrap();
168        assert_eq!(l, matrix.cholesky_decomposition());
169    }
170
171    /// Ensure that the function panics if a non-square matrix is provided
172    #[test]
173    #[should_panic]
174    fn non_square() {
175        let matrix = MatQ::new(3, 2);
176
177        matrix.cholesky_decomposition();
178    }
179
180    /// Ensure that the function panics if a matrix with negative eigenvalues is provided
181    #[test]
182    #[should_panic]
183    fn non_positive_definite() {
184        let matrix: MatQ = -1 * MatQ::identity(2, 2);
185
186        matrix.cholesky_decomposition();
187    }
188
189    /// Ensure that the function panics if a non-symmetric matrix is provided
190    #[test]
191    #[should_panic]
192    fn non_symmetric() {
193        let mut matrix: MatQ = MatQ::identity(2, 2);
194        matrix.set_entry(1, 0, Q::MINUS_ONE).unwrap();
195
196        matrix.cholesky_decomposition();
197    }
198
199    /// Ensure that the function works with large entries
200    #[test]
201    fn large_entries() {
202        // matrix = [[1,-2^32],[-2^{32},2^64+1]] -> L = [[1,0],[-2^32,1]]
203        let matrix: MatQ = MatQ::from_str(&format!(
204            "[[{},-{}],[-{},{}]]",
205            -1,
206            2_i64.pow(32),
207            2_i64.pow(32),
208            u64::MAX
209        ))
210        .unwrap()
211            + 2 * MatQ::identity(2, 2);
212
213        assert_eq!(
214            matrix,
215            (matrix.cholesky_decomposition() * matrix.cholesky_decomposition().transpose())
216        );
217    }
218}
219
220#[cfg(test)]
221mod test_cholesky_decomposition_flint {
222    use crate::{
223        rational::{MatQ, Q},
224        traits::MatrixSetEntry,
225    };
226    use std::str::FromStr;
227
228    /// Ensure that a basic example (from Wikipedia) works.
229    #[test]
230    fn valid_input() {
231        let matrix = MatQ::from_str("[[4, 12, -16],[12,37,-43],[-16,-43,98]]").unwrap();
232
233        let l = MatQ::from_str("[[2, 0, 0],[6, 1, 0],[-8, 5, 3]]").unwrap();
234        assert_eq!(l, matrix.cholesky_decomposition_flint());
235    }
236
237    /// Ensure that the function panics if a non-square matrix is provided
238    #[test]
239    #[should_panic]
240    fn non_square() {
241        let matrix = MatQ::new(3, 2);
242
243        matrix.cholesky_decomposition_flint();
244    }
245
246    /// Ensure that the function panics if a matrix with negative eigenvalues is provided
247    #[test]
248    #[should_panic]
249    fn non_positive_definite() {
250        let matrix: MatQ = -1 * MatQ::identity(2, 2);
251
252        matrix.cholesky_decomposition_flint();
253    }
254
255    /// Ensure that the function panics if a non-symmetric matrix is provided
256    #[test]
257    #[should_panic]
258    fn non_symmetric() {
259        let mut matrix: MatQ = MatQ::identity(2, 2);
260        matrix.set_entry(1, 0, Q::MINUS_ONE).unwrap();
261
262        matrix.cholesky_decomposition_flint();
263    }
264}
qfall_math/rational/mat_q/cholesky_decomp.rs

qfall_math/rational/mat_q/
cholesky_decomp.rs
