mathhook_core/matrices/decomposition/
cholesky.rs

1//! Cholesky decomposition algorithms
2//!
3//! This module provides Cholesky decomposition for symmetric positive definite
4//! matrices, useful for solving linear systems and optimization problems.
5
6use crate::core::Expression;
7use crate::matrices::types::*;
8use crate::matrices::unified::Matrix;
9use crate::simplify::Simplify;
10
11/// Cholesky decomposition implementation
12impl Matrix {
13    /// Perform Cholesky decomposition for positive definite matrices
14    ///
15    /// Decomposes symmetric positive definite matrix A into A = LL^T where:
16    /// - L is lower triangular with positive diagonal elements
17    ///
18    /// # Examples
19    ///
20    /// ```
21    /// use mathhook_core::matrices::Matrix;
22    /// use mathhook_core::Expression;
23    ///
24    /// let matrix = Matrix::from_arrays([
25    ///     [4, 2],
26    ///     [2, 3]
27    /// ]);
28    ///
29    /// if let Some(chol) = matrix.cholesky_decomposition() {
30    ///     let (l_rows, l_cols) = chol.l.dimensions();
31    ///     assert_eq!(l_rows, 2);
32    ///     assert_eq!(l_cols, 2);
33    /// }
34    /// ```
35    pub fn cholesky_decomposition(&self) -> Option<CholeskyDecomposition> {
36        let (rows, cols) = self.dimensions();
37        if rows != cols || !self.is_symmetric() {
38            return None;
39        }
40
41        match self {
42            Matrix::Identity(data) => Some(CholeskyDecomposition {
43                l: Matrix::identity(data.size),
44            }),
45            Matrix::Scalar(data) => {
46                let sqrt_c = Expression::pow(data.scalar_value.clone(), Expression::rational(1, 2));
47                Some(CholeskyDecomposition {
48                    l: Matrix::scalar(data.size, sqrt_c),
49                })
50            }
51            Matrix::Diagonal(data) => {
52                let sqrt_elements: Vec<Expression> = data
53                    .diagonal_elements
54                    .iter()
55                    .map(|elem| Expression::pow(elem.clone(), Expression::rational(1, 2)))
56                    .collect();
57                Some(CholeskyDecomposition {
58                    l: Matrix::diagonal(sqrt_elements),
59                })
60            }
61            _ => {
62                // General Cholesky decomposition
63                self.general_cholesky()
64            }
65        }
66    }
67
68    /// General Cholesky decomposition implementation
69    fn general_cholesky(&self) -> Option<CholeskyDecomposition> {
70        let (n, _) = self.dimensions();
71        let mut l_elements = vec![vec![Expression::integer(0); n]; n];
72
73        for i in 0..n {
74            for j in 0..=i {
75                if i == j {
76                    // L[i][i] = sqrt(A[i][i] - sum(L[i][k]^2 for k < i))
77                    let sum =
78                        l_elements[i]
79                            .iter()
80                            .take(i)
81                            .fold(Expression::integer(0), |acc, l_ik| {
82                                Expression::add(vec![
83                                    acc,
84                                    Expression::pow(l_ik.clone(), Expression::integer(2)),
85                                ])
86                                .simplify()
87                            });
88
89                    let diagonal_val = Expression::add(vec![
90                        self.get_element(i, i),
91                        Expression::mul(vec![Expression::integer(-1), sum]),
92                    ])
93                    .simplify();
94
95                    l_elements[i][i] = Expression::pow(diagonal_val, Expression::rational(1, 2));
96                } else {
97                    // L[i][j] = (A[i][j] - sum(L[i][k]*L[j][k] for k < j)) / L[j][j]
98                    let sum = l_elements[i].iter().zip(l_elements[j].iter()).take(j).fold(
99                        Expression::integer(0),
100                        |acc, (l_ik, l_jk)| {
101                            Expression::add(vec![
102                                acc,
103                                Expression::mul(vec![l_ik.clone(), l_jk.clone()]),
104                            ])
105                            .simplify()
106                        },
107                    );
108
109                    let numerator = Expression::add(vec![
110                        self.get_element(i, j),
111                        Expression::mul(vec![Expression::integer(-1), sum]),
112                    ])
113                    .simplify();
114
115                    // Use canonical form for division: a / b = a * b^(-1)
116                    l_elements[i][j] = Expression::mul(vec![
117                        numerator,
118                        Expression::pow(l_elements[j][j].clone(), Expression::integer(-1)),
119                    ])
120                    .simplify();
121                }
122            }
123        }
124
125        Some(CholeskyDecomposition {
126            l: Matrix::dense(l_elements),
127        })
128    }
129
130    /// Check if matrix is positive definite using Cholesky test
131    ///
132    /// # Examples
133    ///
134    /// ```
135    /// use mathhook_core::matrices::Matrix;
136    /// use mathhook_core::Expression;
137    ///
138    /// let identity = Matrix::identity(3);
139    /// assert!(identity.is_positive_definite_cholesky());
140    ///
141    /// let scalar = Matrix::scalar(2, Expression::integer(5));
142    /// assert!(scalar.is_positive_definite_cholesky());
143    /// ```
144    pub fn is_positive_definite_cholesky(&self) -> bool {
145        if !self.is_symmetric() {
146            return false;
147        }
148
149        match self {
150            Matrix::Identity(_) => true,
151            Matrix::Scalar(data) => {
152                // Check if scalar_value > 0 (simplified check)
153                !data.scalar_value.is_zero() && data.scalar_value != Expression::integer(-1)
154            }
155            Matrix::Diagonal(data) => {
156                // Check if all diagonal elements > 0 (simplified check)
157                data.diagonal_elements.iter().all(|elem| !elem.is_zero())
158            }
159            _ => {
160                // Use Cholesky decomposition test
161                self.cholesky_decomposition().is_some()
162            }
163        }
164    }
165}
mathhook_core/matrices/decomposition/cholesky.rs

mathhook_core/matrices/decomposition/
cholesky.rs
