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                    // Diagonal element: L[i][i] = sqrt(A[i][i] - sum(L[i][k]^2 for k < i))
77                    let mut sum = Expression::integer(0);
78                    for k in 0..i {
79                        let l_ik = l_elements[i][k].clone();
80                        sum = Expression::add(vec![
81                            sum,
82                            Expression::pow(l_ik, Expression::integer(2)),
83                        ])
84                        .simplify();
85                    }
86
87                    let diagonal_val = Expression::add(vec![
88                        self.get_element(i, i),
89                        Expression::mul(vec![Expression::integer(-1), sum]),
90                    ])
91                    .simplify();
92
93                    l_elements[i][i] = Expression::pow(diagonal_val, Expression::rational(1, 2));
94                } else {
95                    // Lower triangular element: L[i][j] = (A[i][j] - sum(L[i][k]*L[j][k] for k < j)) / L[j][j]
96                    let mut sum = Expression::integer(0);
97                    for k in 0..j {
98                        let l_ik = l_elements[i][k].clone();
99                        let l_jk = l_elements[j][k].clone();
100                        sum = Expression::add(vec![sum, Expression::mul(vec![l_ik, l_jk])])
101                            .simplify();
102                    }
103
104                    let numerator = Expression::add(vec![
105                        self.get_element(i, j),
106                        Expression::mul(vec![Expression::integer(-1), sum]),
107                    ])
108                    .simplify();
109
110                    // Use canonical form for division: a / b = a * b^(-1)
111                    l_elements[i][j] = Expression::mul(vec![
112                        numerator,
113                        Expression::pow(l_elements[j][j].clone(), Expression::integer(-1)),
114                    ])
115                    .simplify();
116                }
117            }
118        }
119
120        Some(CholeskyDecomposition {
121            l: Matrix::dense(l_elements),
122        })
123    }
124
125    /// Check if matrix is positive definite using Cholesky test
126    ///
127    /// # Examples
128    ///
129    /// ```
130    /// use mathhook_core::matrices::Matrix;
131    /// use mathhook_core::Expression;
132    ///
133    /// let identity = Matrix::identity(3);
134    /// assert!(identity.is_positive_definite_cholesky());
135    ///
136    /// let scalar = Matrix::scalar(2, Expression::integer(5));
137    /// assert!(scalar.is_positive_definite_cholesky());
138    /// ```
139    pub fn is_positive_definite_cholesky(&self) -> bool {
140        if !self.is_symmetric() {
141            return false;
142        }
143
144        match self {
145            Matrix::Identity(_) => true,
146            Matrix::Scalar(data) => {
147                // Check if scalar_value > 0 (simplified check)
148                !data.scalar_value.is_zero() && data.scalar_value != Expression::integer(-1)
149            }
150            Matrix::Diagonal(data) => {
151                // Check if all diagonal elements > 0 (simplified check)
152                data.diagonal_elements.iter().all(|elem| !elem.is_zero())
153            }
154            _ => {
155                // Use Cholesky decomposition test
156                self.cholesky_decomposition().is_some()
157            }
158        }
159    }
160}
mathhook_core/matrices/decomposition/cholesky.rs

mathhook_core/matrices/decomposition/
cholesky.rs
