qualia_core_db/solvers/linear_algebra/
cholesky.rs

1//! Cholesky decomposition `A = L·Lᵀ` for symmetric positive-definite matrices.
2//!
3//! Functionality parity with nalgebra's `linalg::cholesky`, implemented in the
4//! qualia idiom: **zero allocation**, operating on caller-owned row-major
5//! slices with explicit dimension. No `DMatrix`, no heap, no dependency.
6//!
7//! Cholesky is the fast, numerically-stable path for SPD systems (covariance
8//! solves, least-squares normal equations, Kalman updates, interior-point steps).
9
10use crate::solvers::SolversError;
11
12/// Compute the lower-triangular Cholesky factor `L` of the `n×n` symmetric
13/// positive-definite matrix `a` (row-major), writing `L` row-major into `l`
14/// (lower triangle filled, strictly-upper zeroed). `a` and `l` must each be
15/// length `n*n`.
16///
17/// Returns [`SolversError::SingularMatrix`] if a non-positive pivot is reached
18/// (i.e. `a` is not positive-definite) — fail closed, never a bogus factor.
19/// Only the lower triangle of `a` is read, so a symmetric `a` need not be exact
20/// in its upper half.
21pub fn cholesky_factor(n: usize, a: &[f64], l: &mut [f64]) -> Result<(), SolversError> {
22    if a.len() != n * n || l.len() != n * n {
23        return Err(SolversError::InvalidDimension);
24    }
25    for x in l.iter_mut() {
26        *x = 0.0;
27    }
28    for j in 0..n {
29        // Diagonal: L[j][j] = sqrt(A[j][j] - Σ_{k<j} L[j][k]²)
30        let mut diag = a[j * n + j];
31        for k in 0..j {
32            diag -= l[j * n + k] * l[j * n + k];
33        }
34        if !(diag > 0.0) {
35            return Err(SolversError::SingularMatrix);
36        }
37        let ljj = diag.sqrt();
38        l[j * n + j] = ljj;
39
40        // Below the diagonal: L[i][j] = (A[i][j] - Σ_{k<j} L[i][k]·L[j][k]) / L[j][j]
41        for i in (j + 1)..n {
42            let mut s = a[i * n + j];
43            for k in 0..j {
44                s -= l[i * n + k] * l[j * n + k];
45            }
46            l[i * n + j] = s / ljj;
47        }
48    }
49    Ok(())
50}
51
52/// Solve `A·x = b` for SPD `A`, given its Cholesky factor `l` (from
53/// [`cholesky_factor`]): forward-substitute `L·y = b`, then back-substitute
54/// `Lᵀ·x = y`. `l` is `n*n`; `b` and `x` are length `n`. The solution is written
55/// into `x` (which is also used as scratch for `y`).
56pub fn cholesky_solve(n: usize, l: &[f64], b: &[f64], x: &mut [f64]) -> Result<(), SolversError> {
57    if l.len() != n * n || b.len() != n || x.len() != n {
58        return Err(SolversError::InvalidDimension);
59    }
60    // Forward: L·y = b   (y accumulated in x)
61    for i in 0..n {
62        let mut s = b[i];
63        for k in 0..i {
64            s -= l[i * n + k] * x[k];
65        }
66        x[i] = s / l[i * n + i];
67    }
68    // Backward: Lᵀ·x = y
69    for i in (0..n).rev() {
70        let mut s = x[i];
71        for k in (i + 1)..n {
72            s -= l[k * n + i] * x[k];
73        }
74        x[i] = s / l[i * n + i];
75    }
76    Ok(())
77}
78
79/// Determinant of an SPD matrix from its Cholesky factor: `det(A) = Π L[i][i]²`.
80pub fn cholesky_determinant(n: usize, l: &[f64]) -> f64 {
81    let mut prod = 1.0;
82    for i in 0..n {
83        let d = l[i * n + i];
84        prod *= d * d;
85    }
86    prod
87}
88
89#[cfg(test)]
90mod tests {
91    use super::*;
92
93    const EPS: f64 = 1e-9;
94
95    // Textbook SPD matrix with a known integer factor.
96    // A = [[4,12,-16],[12,37,-43],[-16,-43,98]] = L·Lᵀ,
97    // L = [[2,0,0],[6,1,0],[-8,5,3]].
98    const A3: [f64; 9] = [4.0, 12.0, -16.0, 12.0, 37.0, -43.0, -16.0, -43.0, 98.0];
99
100    #[test]
101    fn factor_matches_known_lower_triangle() {
102        let mut l = [0.0; 9];
103        cholesky_factor(3, &A3, &mut l).unwrap();
104        let expect = [2.0, 0.0, 0.0, 6.0, 1.0, 0.0, -8.0, 5.0, 3.0];
105        for i in 0..9 {
106            assert!(
107                (l[i] - expect[i]).abs() < EPS,
108                "l[{i}]={} != {}",
109                l[i],
110                expect[i]
111            );
112        }
113    }
114
115    #[test]
116    fn reconstructs_a() {
117        let mut l = [0.0; 9];
118        cholesky_factor(3, &A3, &mut l).unwrap();
119        // L·Lᵀ == A
120        for i in 0..3 {
121            for j in 0..3 {
122                let mut s = 0.0;
123                for k in 0..3 {
124                    s += l[i * 3 + k] * l[j * 3 + k];
125                }
126                assert!((s - A3[i * 3 + j]).abs() < 1e-6);
127            }
128        }
129    }
130
131    #[test]
132    fn solves_linear_system() {
133        let mut l = [0.0; 9];
134        cholesky_factor(3, &A3, &mut l).unwrap();
135        let b = [1.0, 2.0, 3.0];
136        let mut x = [0.0; 3];
137        cholesky_solve(3, &l, &b, &mut x).unwrap();
138        // Verify A·x ≈ b
139        for i in 0..3 {
140            let mut s = 0.0;
141            for j in 0..3 {
142                s += A3[i * 3 + j] * x[j];
143            }
144            assert!((s - b[i]).abs() < 1e-6, "row {i}: {} != {}", s, b[i]);
145        }
146    }
147
148    #[test]
149    fn determinant_via_factor() {
150        let mut l = [0.0; 9];
151        cholesky_factor(3, &A3, &mut l).unwrap();
152        // det(A) = (2·1·3)² = 36
153        assert!((cholesky_determinant(3, &l) - 36.0).abs() < 1e-6);
154    }
155
156    #[test]
157    fn rejects_non_positive_definite() {
158        // Symmetric but indefinite (negative eigenvalue): [[1,2],[2,1]].
159        let a = [1.0, 2.0, 2.0, 1.0];
160        let mut l = [0.0; 4];
161        assert_eq!(
162            cholesky_factor(2, &a, &mut l),
163            Err(SolversError::SingularMatrix)
164        );
165    }
166
167    #[test]
168    fn rejects_bad_dims() {
169        let mut l = [0.0; 9];
170        assert_eq!(
171            cholesky_factor(2, &A3, &mut l),
172            Err(SolversError::InvalidDimension)
173        );
174    }
175}
qualia_core_db/solvers/linear_algebra/cholesky.rs

qualia_core_db/solvers/linear_algebra/
cholesky.rs
