From 5592c9f7e8f966d13e1c26cf947fbd9e23c1b2bd Mon Sep 17 00:00:00 2001 From: "Andreas J. Weiss" Date: Tue, 6 May 2025 13:51:35 -0400 Subject: [PATCH 1/8] basic ldl implementation, based on udu --- src/linalg/decomposition.rs | 14 ++++- src/linalg/ldl.rs | 103 ++++++++++++++++++++++++++++++++++++ src/linalg/mod.rs | 2 + tests/linalg/ldl.rs | 76 ++++++++++++++++++++++++++ tests/linalg/mod.rs | 1 + 5 files changed, 195 insertions(+), 1 deletion(-) create mode 100644 src/linalg/ldl.rs create mode 100644 tests/linalg/ldl.rs diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index 957d38ab1..716fab221 100644 --- a/src/linalg/decomposition.rs +++ b/src/linalg/decomposition.rs @@ -2,7 +2,7 @@ use crate::storage::Storage; use crate::{ Allocator, Bidiagonal, Cholesky, ColPivQR, ComplexField, DefaultAllocator, Dim, DimDiff, DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, Matrix, OMatrix, RealField, Schur, - SymmetricEigen, SymmetricTridiagonal, LU, QR, SVD, U1, UDU, + SymmetricEigen, SymmetricTridiagonal, LU, QR, SVD, U1, UDU, LDL }; /// # Rectangular matrix decomposition @@ -281,6 +281,18 @@ impl> Matrix { Hessenberg::new(self.into_owned()) } + /// Attempts to compute the LDL decomposition of this matrix. + /// + /// The input matrix `self` is assumed to be symmetric and this decomposition will only read + /// the lower-triangular part of `self`. + pub fn ldl(self) -> Option> + where + T: RealField, + DefaultAllocator: Allocator + Allocator, + { + LDL::new(self.into_owned()) + } + /// Computes the Schur decomposition of a square matrix. pub fn schur(self) -> Schur where diff --git a/src/linalg/ldl.rs b/src/linalg/ldl.rs new file mode 100644 index 000000000..38967afbd --- /dev/null +++ b/src/linalg/ldl.rs @@ -0,0 +1,103 @@ +#[cfg(feature = "serde-serialize-no-std")] +use serde::{Deserialize, Serialize}; + +use crate::allocator::Allocator; +use crate::base::{Const, DefaultAllocator, OMatrix, OVector}; +use crate::dimension::Dim; +use simba::scalar::RealField; +use std::cmp::Ordering; + +/// LDL factorization. +#[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] +#[cfg_attr( + feature = "serde-serialize-no-std", + serde(bound(serialize = "OVector: Serialize, OMatrix: Serialize")) +)] +#[cfg_attr( + feature = "serde-serialize-no-std", + serde(bound( + deserialize = "OVector: Deserialize<'de>, OMatrix: Deserialize<'de>" + )) +)] +#[derive(Clone, Debug)] +pub struct LDL +where + DefaultAllocator: Allocator + Allocator, +{ + /// The lower triangular matrix resulting from the factorization + pub l: OMatrix, + /// The diagonal matrix resulting from the factorization + pub d: OVector, +} + +impl Copy for LDL +where + DefaultAllocator: Allocator + Allocator, + OVector: Copy, + OMatrix: Copy, +{ +} + +impl LDL +where + DefaultAllocator: Allocator + Allocator, +{ + /// Computes the LDL^T factorization. + /// + /// The input matrix `p` is assumed to be symmetric and this decomposition will only read + /// the lower-triangular part of `p`. + pub fn new(p: OMatrix) -> Option { + let n = p.ncols(); + + let n_dim = p.shape_generic().1; + + let mut d = OVector::::zeros_generic(n_dim, Const::<1>); + let mut l = OMatrix::::zeros_generic(n_dim, n_dim); + + for j in 0..n { + let mut d_j = p[(j, j)].clone(); + + if j > 0 { + for k in 0..j { + d_j -= d[k].clone() * l[(j, k)].clone().powi(2); + } + } + + d[j] = d_j; + + for i in j..n { + let mut l_ij = p[(j, i)].clone(); + + for k in 0..j { + l_ij -= d[k].clone() * l[(j, k)].clone() * l[(i, k)].clone(); + } + + if matches!(d[j].partial_cmp(&T::zero())?, Ordering::Equal) { + l[(i, j)] = T::zero(); + } else { + l[(i, j)] = l_ij / d[j].clone(); + } + } + } + + Some(Self { l, d }) + } + + /// Returns the lower triangular matrix as if generated by the Cholesky decomposition. + pub fn cholesky_l(&self) -> OMatrix { + let n_dim = self.l.shape_generic().1; + + &self.l + * OMatrix::from_diagonal(&OVector::from_iterator_generic( + n_dim, + Const::<1>, + self.d.iter().map(|value| value.clone().sqrt()), + )) + } + + /// Returns the diagonal elements as a matrix + #[must_use] + pub fn d_matrix(&self) -> OMatrix { + OMatrix::from_diagonal(&self.d) + } +} diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index ac1cbb9b2..d8eac0fc1 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -17,6 +17,7 @@ pub mod givens; mod hessenberg; pub mod householder; mod inverse; +mod ldl; mod lu; mod permutation_sequence; mod pow; @@ -39,6 +40,7 @@ pub use self::cholesky::*; pub use self::col_piv_qr::*; pub use self::full_piv_lu::*; pub use self::hessenberg::*; +pub use self::ldl::*; pub use self::lu::*; pub use self::permutation_sequence::*; pub use self::qr::*; diff --git a/tests/linalg/ldl.rs b/tests/linalg/ldl.rs new file mode 100644 index 000000000..f7f50649a --- /dev/null +++ b/tests/linalg/ldl.rs @@ -0,0 +1,76 @@ +use na::Matrix3; + +#[test] +#[rustfmt::skip] +fn ldl_simple() { + let m = Matrix3::new( + 2.0, -1.0, 0.0, + -1.0, 2.0, -1.0, + 0.0, -1.0, 2.0); + + let ldl = m.ldl().unwrap(); + + // Rebuild + let p = ldl.l * ldl.d_matrix() * ldl.l.transpose(); + + assert!(relative_eq!(m, p, epsilon = 3.0e-16)); +} + +#[test] +#[should_panic] +#[rustfmt::skip] +fn ldl_non_sym_panic() { + let m = Matrix3::new( + 2.0, -1.0, 0.0, + 1.0, -2.0, 3.0, + -2.0, 1.0, 0.3); + + let ldl = m.ldl().unwrap(); + // Rebuild + let p = ldl.l * ldl.d_matrix() * ldl.l.transpose(); + + assert!(relative_eq!(m, p, epsilon = 3.0e-16)); +} + +#[cfg(feature = "proptest-support")] +mod proptest_tests { + #[allow(unused_imports)] + use crate::core::helper::{RandComplex, RandScalar}; + + macro_rules! gen_tests( + ($module: ident, $scalar: expr) => { + mod $module { + #[allow(unused_imports)] + use crate::core::helper::{RandScalar, RandComplex}; + use crate::proptest::*; + use proptest::{prop_assert, proptest}; + + proptest! { + #[test] + fn ldl(m in dmatrix_($scalar)) { + let m = &m * m.adjoint(); + + if let Some(ldl) = m.clone().ldl() { + let p = &ldl.l * &ldl.d_matrix() * &ldl.l.transpose(); + println!("m: {}, p: {}", m, p); + + prop_assert!(relative_eq!(m, p, epsilon = 1.0e-7)); + } + } + + #[test] + fn ldl_static(m in matrix4_($scalar)) { + let m = m.hermitian_part(); + + if let Some(ldl) = m.ldl() { + let p = ldl.l * ldl.d_matrix() * ldl.l.transpose(); + prop_assert!(relative_eq!(m, p, epsilon = 1.0e-7)); + } + } + } + } + } + ); + + gen_tests!(f64, PROPTEST_F64); +} diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs index d9bd6cd91..33f6668f7 100644 --- a/tests/linalg/mod.rs +++ b/tests/linalg/mod.rs @@ -8,6 +8,7 @@ mod exp; mod full_piv_lu; mod hessenberg; mod inverse; +mod ldl; mod lu; mod pow; mod qr; From b8759354eafd1bd22b745a5dff9d3fd50aa7c834 Mon Sep 17 00:00:00 2001 From: "Andreas J. Weiss" Date: Tue, 6 May 2025 13:57:44 -0400 Subject: [PATCH 2/8] test for cholesky eq --- tests/linalg/ldl.rs | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/linalg/ldl.rs b/tests/linalg/ldl.rs index f7f50649a..cbfd36df2 100644 --- a/tests/linalg/ldl.rs +++ b/tests/linalg/ldl.rs @@ -16,6 +16,20 @@ fn ldl_simple() { assert!(relative_eq!(m, p, epsilon = 3.0e-16)); } +#[test] +#[rustfmt::skip] +fn ldl_cholesky() { + let m = Matrix3::new( + 2.0, -1.0, 0.0, + -1.0, 2.0, -1.0, + 0.0, -1.0, 2.0); + + let chol= m.cholesky().unwrap(); + let ldl = m.ldl().unwrap(); + + assert!(relative_eq!(ldl.cholesky_l(), chol.l(), epsilon = 3.0e-16)); +} + #[test] #[should_panic] #[rustfmt::skip] @@ -26,6 +40,7 @@ fn ldl_non_sym_panic() { -2.0, 1.0, 0.3); let ldl = m.ldl().unwrap(); + // Rebuild let p = ldl.l * ldl.d_matrix() * ldl.l.transpose(); From 00ca97b0df453d64b2d06f9eda6c04a7c3404c3b Mon Sep 17 00:00:00 2001 From: "Andreas J. Weiss" Date: Wed, 14 May 2025 09:51:30 -0400 Subject: [PATCH 3/8] hermitian ldl, better tests --- src/linalg/decomposition.rs | 4 ++-- src/linalg/ldl.rs | 19 +++++++++---------- tests/linalg/ldl.rs | 37 +++++++++++++++++++++++++++---------- 3 files changed, 38 insertions(+), 22 deletions(-) diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index 716fab221..583588a77 100644 --- a/src/linalg/decomposition.rs +++ b/src/linalg/decomposition.rs @@ -2,7 +2,7 @@ use crate::storage::Storage; use crate::{ Allocator, Bidiagonal, Cholesky, ColPivQR, ComplexField, DefaultAllocator, Dim, DimDiff, DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, Matrix, OMatrix, RealField, Schur, - SymmetricEigen, SymmetricTridiagonal, LU, QR, SVD, U1, UDU, LDL + SymmetricEigen, SymmetricTridiagonal, LDL, LU, QR, SVD, U1, UDU, }; /// # Rectangular matrix decomposition @@ -287,7 +287,7 @@ impl> Matrix { /// the lower-triangular part of `self`. pub fn ldl(self) -> Option> where - T: RealField, + T: ComplexField, DefaultAllocator: Allocator + Allocator, { LDL::new(self.into_owned()) diff --git a/src/linalg/ldl.rs b/src/linalg/ldl.rs index 38967afbd..d029b71f4 100644 --- a/src/linalg/ldl.rs +++ b/src/linalg/ldl.rs @@ -4,8 +4,7 @@ use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; use crate::base::{Const, DefaultAllocator, OMatrix, OVector}; use crate::dimension::Dim; -use simba::scalar::RealField; -use std::cmp::Ordering; +use simba::scalar::ComplexField; /// LDL factorization. #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] @@ -20,7 +19,7 @@ use std::cmp::Ordering; )) )] #[derive(Clone, Debug)] -pub struct LDL +pub struct LDL where DefaultAllocator: Allocator + Allocator, { @@ -30,7 +29,7 @@ where pub d: OVector, } -impl Copy for LDL +impl Copy for LDL where DefaultAllocator: Allocator + Allocator, OVector: Copy, @@ -38,13 +37,13 @@ where { } -impl LDL +impl LDL where DefaultAllocator: Allocator + Allocator, { /// Computes the LDL^T factorization. /// - /// The input matrix `p` is assumed to be symmetric and this decomposition will only read + /// The input matrix `p` is assumed to be hermitian c and this decomposition will only read /// the lower-triangular part of `p`. pub fn new(p: OMatrix) -> Option { let n = p.ncols(); @@ -59,20 +58,20 @@ where if j > 0 { for k in 0..j { - d_j -= d[k].clone() * l[(j, k)].clone().powi(2); + d_j -= l[(j, k)].clone() * l[(j, k)].clone().conjugate() * d[k].clone(); } } d[j] = d_j; for i in j..n { - let mut l_ij = p[(j, i)].clone(); + let mut l_ij = p[(i, j)].clone(); for k in 0..j { - l_ij -= d[k].clone() * l[(j, k)].clone() * l[(i, k)].clone(); + l_ij -= l[(j, k)].clone().conjugate() * l[(i, k)].clone() * d[k].clone(); } - if matches!(d[j].partial_cmp(&T::zero())?, Ordering::Equal) { + if d[j] == T::zero() { l[(i, j)] = T::zero(); } else { l[(i, j)] = l_ij / d[j].clone(); diff --git a/tests/linalg/ldl.rs b/tests/linalg/ldl.rs index cbfd36df2..6f3571882 100644 --- a/tests/linalg/ldl.rs +++ b/tests/linalg/ldl.rs @@ -1,28 +1,45 @@ -use na::Matrix3; +use na::{Complex, Matrix3}; +use num::Zero; #[test] #[rustfmt::skip] fn ldl_simple() { let m = Matrix3::new( - 2.0, -1.0, 0.0, - -1.0, 2.0, -1.0, - 0.0, -1.0, 2.0); + Complex::new(2.0, 0.0), Complex::new(-1.0, 0.5), Complex::zero(), + Complex::new(-1.0, -0.5), Complex::new(2.0, 0.0), Complex::new(-1.0, 0.0), + Complex::zero(), Complex::new(-1.0, 0.0), Complex::new(2.0, 0.0)); - let ldl = m.ldl().unwrap(); + let ldl = m.lower_triangle().ldl().unwrap(); // Rebuild - let p = ldl.l * ldl.d_matrix() * ldl.l.transpose(); + let p = ldl.l * ldl.d_matrix() * ldl.l.adjoint(); - assert!(relative_eq!(m, p, epsilon = 3.0e-16)); + assert!(relative_eq!(m, p, epsilon = 3.0e-12)); +} + +#[test] +#[rustfmt::skip] +fn ldl_partial() { + let m = Matrix3::new( + Complex::new(2.0, 0.0), Complex::zero(), Complex::zero(), + Complex::zero(), Complex::zero(), Complex::zero(), + Complex::zero(), Complex::zero(), Complex::new(2.0, 0.0)); + + let ldl = m.lower_triangle().ldl().unwrap(); + + // Rebuild + let p = ldl.l * ldl.d_matrix() * ldl.l.adjoint(); + + assert!(relative_eq!(m, p, epsilon = 3.0e-12)); } #[test] #[rustfmt::skip] fn ldl_cholesky() { let m = Matrix3::new( - 2.0, -1.0, 0.0, - -1.0, 2.0, -1.0, - 0.0, -1.0, 2.0); + Complex::new(2.0, 0.0), Complex::new(-1.0, 0.5), Complex::zero(), + Complex::new(-1.0, -0.5), Complex::new(2.0, 0.0), Complex::new(-1.0, 0.0), + Complex::zero(), Complex::new(-1.0, 0.0), Complex::new(2.0, 0.0)); let chol= m.cholesky().unwrap(); let ldl = m.ldl().unwrap(); From d562fbdba82cc6858520d38f8675f0b97dc4a69e Mon Sep 17 00:00:00 2001 From: "Andreas J. Weiss" Date: Tue, 25 Nov 2025 11:29:28 -0500 Subject: [PATCH 4/8] update docs, renamed cholesky_l to lsqrt and added extra checks --- src/linalg/ldl.rs | 31 ++++++++++++++++++++++--------- tests/linalg/ldl.rs | 4 ++-- 2 files changed, 24 insertions(+), 11 deletions(-) diff --git a/src/linalg/ldl.rs b/src/linalg/ldl.rs index d029b71f4..45f4b1868 100644 --- a/src/linalg/ldl.rs +++ b/src/linalg/ldl.rs @@ -6,7 +6,8 @@ use crate::base::{Const, DefaultAllocator, OMatrix, OVector}; use crate::dimension::Dim; use simba::scalar::ComplexField; -/// LDL factorization. +/// The LDL / LDL^T factorization of a Hermitian matrix A = LDL^T where L is a +/// lower unit-triangular matrix and D is diagonal matrix. #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize-no-std", @@ -25,7 +26,7 @@ where { /// The lower triangular matrix resulting from the factorization pub l: OMatrix, - /// The diagonal matrix resulting from the factorization + /// The diagonal matrix, in vector form, resulting from the factorization pub d: OVector, } @@ -41,9 +42,9 @@ impl LDL where DefaultAllocator: Allocator + Allocator, { - /// Computes the LDL^T factorization. + /// Computes the LDL / LDL^T factorization. /// - /// The input matrix `p` is assumed to be hermitian c and this decomposition will only read + /// The input matrix `p` is assumed to be hermitian and this decomposition will only read /// the lower-triangular part of `p`. pub fn new(p: OMatrix) -> Option { let n = p.ncols(); @@ -82,19 +83,31 @@ where Some(Self { l, d }) } - /// Returns the lower triangular matrix as if generated by the Cholesky decomposition. - pub fn cholesky_l(&self) -> OMatrix { + /// Returns the matrix L * sqrt(D). + /// + /// This function can be used to generate a lower triangular matrix as if it were + /// generated by the Cholesky decomposition without the requirement of positive definiteness. + /// + /// This function returns `None` if the square root of any of the values in the diagonal matrix D is not finite. + pub fn lsqrtd(&self) -> Option> { let n_dim = self.l.shape_generic().1; - &self.l + let lsqrtd = &self.l * OMatrix::from_diagonal(&OVector::from_iterator_generic( n_dim, Const::<1>, self.d.iter().map(|value| value.clone().sqrt()), - )) + )); + + // Check for any non-finite numbers in lsqrtd and return None if necessary. + if !lsqrtd.iter().fold(true, |acc, next| acc & next.is_finite()) { + None + } else { + Some(lsqrtd) + } } - /// Returns the diagonal elements as a matrix + /// Returns the diagonal elements as a matrix. #[must_use] pub fn d_matrix(&self) -> OMatrix { OMatrix::from_diagonal(&self.d) diff --git a/tests/linalg/ldl.rs b/tests/linalg/ldl.rs index 6f3571882..271e91a1f 100644 --- a/tests/linalg/ldl.rs +++ b/tests/linalg/ldl.rs @@ -35,7 +35,7 @@ fn ldl_partial() { #[test] #[rustfmt::skip] -fn ldl_cholesky() { +fn ldl_lsqrtd() { let m = Matrix3::new( Complex::new(2.0, 0.0), Complex::new(-1.0, 0.5), Complex::zero(), Complex::new(-1.0, -0.5), Complex::new(2.0, 0.0), Complex::new(-1.0, 0.0), @@ -44,7 +44,7 @@ fn ldl_cholesky() { let chol= m.cholesky().unwrap(); let ldl = m.ldl().unwrap(); - assert!(relative_eq!(ldl.cholesky_l(), chol.l(), epsilon = 3.0e-16)); + assert!(relative_eq!(ldl.lsqrtd().unwrap(), chol.l(), epsilon = 3.0e-16)); } #[test] From 43fd28073385eba7294360bcc3ba5e3ed79e039c Mon Sep 17 00:00:00 2001 From: "Andreas J. Weiss" Date: Wed, 3 Dec 2025 11:28:21 -0500 Subject: [PATCH 5/8] in place ldl factorization --- src/linalg/ldl.rs | 128 +++++++++++++++++++++++--------------------- tests/linalg/ldl.rs | 15 ++++-- 2 files changed, 78 insertions(+), 65 deletions(-) diff --git a/src/linalg/ldl.rs b/src/linalg/ldl.rs index 45f4b1868..217171054 100644 --- a/src/linalg/ldl.rs +++ b/src/linalg/ldl.rs @@ -6,34 +6,24 @@ use crate::base::{Const, DefaultAllocator, OMatrix, OVector}; use crate::dimension::Dim; use simba::scalar::ComplexField; -/// The LDL / LDL^T factorization of a Hermitian matrix A = LDL^T where L is a -/// lower unit-triangular matrix and D is diagonal matrix. +/// The LDL / LDL^T factorization of a Hermitian matrix A = LDL^T where L is a lower unit-triangular matrix and D is diagonal matrix. #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize-no-std", - serde(bound(serialize = "OVector: Serialize, OMatrix: Serialize")) + serde(bound(serialize = "OMatrix: Serialize")) )] #[cfg_attr( feature = "serde-serialize-no-std", - serde(bound( - deserialize = "OVector: Deserialize<'de>, OMatrix: Deserialize<'de>" - )) + serde(bound(deserialize = "OMatrix: Deserialize<'de>")) )] #[derive(Clone, Debug)] -pub struct LDL +pub struct LDL(OMatrix) where - DefaultAllocator: Allocator + Allocator, -{ - /// The lower triangular matrix resulting from the factorization - pub l: OMatrix, - /// The diagonal matrix, in vector form, resulting from the factorization - pub d: OVector, -} + DefaultAllocator: Allocator + Allocator; impl Copy for LDL where DefaultAllocator: Allocator + Allocator, - OVector: Copy, OMatrix: Copy, { } @@ -42,61 +32,45 @@ impl LDL where DefaultAllocator: Allocator + Allocator, { - /// Computes the LDL / LDL^T factorization. - /// - /// The input matrix `p` is assumed to be hermitian and this decomposition will only read - /// the lower-triangular part of `p`. - pub fn new(p: OMatrix) -> Option { - let n = p.ncols(); - - let n_dim = p.shape_generic().1; - - let mut d = OVector::::zeros_generic(n_dim, Const::<1>); - let mut l = OMatrix::::zeros_generic(n_dim, n_dim); - - for j in 0..n { - let mut d_j = p[(j, j)].clone(); - - if j > 0 { - for k in 0..j { - d_j -= l[(j, k)].clone() * l[(j, k)].clone().conjugate() * d[k].clone(); - } - } + /// Returns the diagonal elements as a vector. + #[must_use] + pub fn d(&self) -> OVector { + self.0.diagonal() + } - d[j] = d_j; + /// Returns the diagonal elements as a matrix. + #[must_use] + pub fn d_matrix(&self) -> OMatrix { + OMatrix::from_diagonal(&self.0.diagonal()) + } - for i in j..n { - let mut l_ij = p[(i, j)].clone(); + /// Returns the lower triangular matrix. + #[must_use] + pub fn l_matrix(&self) -> OMatrix { + let mut l = self.0.clone(); - for k in 0..j { - l_ij -= l[(j, k)].clone().conjugate() * l[(i, k)].clone() * d[k].clone(); - } + l.column_iter_mut() + .enumerate() + .for_each(|(idx, mut column)| { + column[idx] = T::one(); + }); - if d[j] == T::zero() { - l[(i, j)] = T::zero(); - } else { - l[(i, j)] = l_ij / d[j].clone(); - } - } - } - - Some(Self { l, d }) + l } /// Returns the matrix L * sqrt(D). - /// - /// This function can be used to generate a lower triangular matrix as if it were - /// generated by the Cholesky decomposition without the requirement of positive definiteness. - /// /// This function returns `None` if the square root of any of the values in the diagonal matrix D is not finite. + /// + /// This function can be used to generate a lower triangular matrix as if it were generated by the Cholesky decomposition, without the requirement of positive definiteness. pub fn lsqrtd(&self) -> Option> { - let n_dim = self.l.shape_generic().1; + let n_dim = self.0.shape_generic().1; - let lsqrtd = &self.l + let lsqrtd: crate::Matrix>::Buffer> = &self + .l_matrix() * OMatrix::from_diagonal(&OVector::from_iterator_generic( n_dim, Const::<1>, - self.d.iter().map(|value| value.clone().sqrt()), + self.d().iter().map(|value| value.clone().sqrt()), )); // Check for any non-finite numbers in lsqrtd and return None if necessary. @@ -107,9 +81,41 @@ where } } - /// Returns the diagonal elements as a matrix. - #[must_use] - pub fn d_matrix(&self) -> OMatrix { - OMatrix::from_diagonal(&self.d) + /// Computes the LDL / LDL^T factorization. + pub fn new(mut matrix: OMatrix) -> Option { + for j in 0..matrix.ncols() { + let mut d_j: T = matrix[(j, j)].clone(); + + if j > 0 { + for k in 0..j { + d_j -= matrix[(j, k)].clone() + * matrix[(j, k)].clone().conjugate() + * matrix[(k, k)].clone(); + } + } + + matrix[(j, j)] = d_j; + + for i in (j + 1)..matrix.ncols() { + let mut l_ij = matrix[(i, j)].clone(); + + for k in 0..j { + l_ij -= matrix[(j, k)].clone().conjugate() + * matrix[(i, k)].clone() + * matrix[(k, k)].clone(); + } + + if matrix[(j, j)] == T::zero() { + matrix[(i, j)] = T::zero(); + } else { + matrix[(i, j)] = l_ij / matrix[(j, j)].clone(); + } + + // Zero out the upper triangular part. + matrix[(j, i)] = T::zero(); + } + } + + Some(Self(matrix)) } } diff --git a/tests/linalg/ldl.rs b/tests/linalg/ldl.rs index 271e91a1f..d51630db4 100644 --- a/tests/linalg/ldl.rs +++ b/tests/linalg/ldl.rs @@ -9,10 +9,17 @@ fn ldl_simple() { Complex::new(-1.0, -0.5), Complex::new(2.0, 0.0), Complex::new(-1.0, 0.0), Complex::zero(), Complex::new(-1.0, 0.0), Complex::new(2.0, 0.0)); - let ldl = m.lower_triangle().ldl().unwrap(); + let ldl = m.ldl().unwrap(); + + println!("{:}", &m); + println!("{:}", ldl.l_matrix()); + println!("{:}", ldl.d()); // Rebuild - let p = ldl.l * ldl.d_matrix() * ldl.l.adjoint(); + let p = ldl.l_matrix() * ldl.d_matrix() * ldl.l_matrix().adjoint(); + + + println!("{:}", &p); assert!(relative_eq!(m, p, epsilon = 3.0e-12)); } @@ -28,7 +35,7 @@ fn ldl_partial() { let ldl = m.lower_triangle().ldl().unwrap(); // Rebuild - let p = ldl.l * ldl.d_matrix() * ldl.l.adjoint(); + let p = ldl.l_matrix() * ldl.d_matrix() * ldl.l_matrix().adjoint(); assert!(relative_eq!(m, p, epsilon = 3.0e-12)); } @@ -59,7 +66,7 @@ fn ldl_non_sym_panic() { let ldl = m.ldl().unwrap(); // Rebuild - let p = ldl.l * ldl.d_matrix() * ldl.l.transpose(); + let p = ldl.l_matrix() * ldl.d_matrix() * ldl.l_matrix().transpose(); assert!(relative_eq!(m, p, epsilon = 3.0e-16)); } From 973110782e7529c56793366bc1132eb0cc085bc3 Mon Sep 17 00:00:00 2001 From: "Andreas J. Weiss" Date: Mon, 13 Apr 2026 14:37:02 -0400 Subject: [PATCH 6/8] ldl docs and unchecked --- src/linalg/decomposition.rs | 2 +- src/linalg/ldl.rs | 73 ++++++++++++++++++++++++++++++------- tests/linalg/ldl.rs | 19 +++------- 3 files changed, 66 insertions(+), 28 deletions(-) diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index 583588a77..82ebaa166 100644 --- a/src/linalg/decomposition.rs +++ b/src/linalg/decomposition.rs @@ -281,7 +281,7 @@ impl> Matrix { Hessenberg::new(self.into_owned()) } - /// Attempts to compute the LDL decomposition of this matrix. + /// Attempts to compute the strictly diagonal LDL / LDLT decomposition of this matrix. /// /// The input matrix `self` is assumed to be symmetric and this decomposition will only read /// the lower-triangular part of `self`. diff --git a/src/linalg/ldl.rs b/src/linalg/ldl.rs index 217171054..dd3770516 100644 --- a/src/linalg/ldl.rs +++ b/src/linalg/ldl.rs @@ -7,6 +7,9 @@ use crate::dimension::Dim; use simba::scalar::ComplexField; /// The LDL / LDL^T factorization of a Hermitian matrix A = LDL^T where L is a lower unit-triangular matrix and D is diagonal matrix. +/// +/// This implementation referes to the strictly diagonal LDL algorithm. +/// For the more general Bunch-Kaufman block-diagonal pivoting method we refer to the LBL factorization. #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize-no-std", @@ -32,21 +35,21 @@ impl LDL where DefaultAllocator: Allocator + Allocator, { - /// Returns the diagonal elements as a vector. + /// Returns the diagonal elements as a matrix. #[must_use] - pub fn d(&self) -> OVector { - self.0.diagonal() + pub fn d(&self) -> OMatrix { + OMatrix::from_diagonal(&self.0.diagonal()) } - /// Returns the diagonal elements as a matrix. + /// Returns the diagonal elements as a vector. #[must_use] - pub fn d_matrix(&self) -> OMatrix { - OMatrix::from_diagonal(&self.0.diagonal()) + pub fn diagonal(&self) -> OVector { + self.0.diagonal() } /// Returns the lower triangular matrix. #[must_use] - pub fn l_matrix(&self) -> OMatrix { + pub fn l(&self) -> OMatrix { let mut l = self.0.clone(); l.column_iter_mut() @@ -61,16 +64,16 @@ where /// Returns the matrix L * sqrt(D). /// This function returns `None` if the square root of any of the values in the diagonal matrix D is not finite. /// - /// This function can be used to generate a lower triangular matrix as if it were generated by the Cholesky decomposition, without the requirement of positive definiteness. + /// This function can be used to generate a lower triangular matrix as if it were generated by the Cholesky decomposition. pub fn lsqrtd(&self) -> Option> { let n_dim = self.0.shape_generic().1; let lsqrtd: crate::Matrix>::Buffer> = &self - .l_matrix() + .l() * OMatrix::from_diagonal(&OVector::from_iterator_generic( n_dim, Const::<1>, - self.d().iter().map(|value| value.clone().sqrt()), + self.0.diagonal().iter().map(|value| value.clone().sqrt()), )); // Check for any non-finite numbers in lsqrtd and return None if necessary. @@ -94,7 +97,11 @@ where } } - matrix[(j, j)] = d_j; + if d_j.is_zero() { + return None; + } + + matrix[(j, j)] = d_j.clone(); for i in (j + 1)..matrix.ncols() { let mut l_ij = matrix[(i, j)].clone(); @@ -105,10 +112,48 @@ where * matrix[(k, k)].clone(); } - if matrix[(j, j)] == T::zero() { + matrix[(i, j)] = l_ij / matrix[(j, j)].clone(); + + dbg!(&matrix[(i, j)]); + + // Zero out the upper triangular part. + matrix[(j, i)] = T::zero(); + } + } + + Some(Self(matrix)) + } + + /// Computes the LDL / LDL^T factorization, allowing for zero pivots. If a zero pivot is detected the corresponding column & row is zero'd out. + /// + /// This factorization is only valid for singular matrices with a specific structure, e.g. singular matrices of the form A = xxT. + pub fn new_unchecked(mut matrix: OMatrix) -> Self { + for j in 0..matrix.ncols() { + let mut d_j: T = matrix[(j, j)].clone(); + + if j > 0 { + for k in 0..j { + d_j -= matrix[(j, k)].clone() + * matrix[(j, k)].clone().conjugate() + * matrix[(k, k)].clone(); + } + } + + matrix[(j, j)] = d_j.clone(); + + for i in (j + 1)..matrix.ncols() { + let mut l_ij = matrix[(i, j)].clone(); + + for k in 0..j { + l_ij -= matrix[(j, k)].clone().conjugate() + * matrix[(i, k)].clone() + * matrix[(k, k)].clone(); + } + + if d_j.is_zero() { matrix[(i, j)] = T::zero(); } else { - matrix[(i, j)] = l_ij / matrix[(j, j)].clone(); + matrix[(i, j)] = l_ij / d_j.clone(); } // Zero out the upper triangular part. @@ -116,6 +161,6 @@ where } } - Some(Self(matrix)) + Self(matrix) } } diff --git a/tests/linalg/ldl.rs b/tests/linalg/ldl.rs index d51630db4..48da18cac 100644 --- a/tests/linalg/ldl.rs +++ b/tests/linalg/ldl.rs @@ -1,4 +1,4 @@ -use na::{Complex, Matrix3}; +use na::{linalg::LDL, Complex, Matrix3}; use num::Zero; #[test] @@ -10,16 +10,9 @@ fn ldl_simple() { Complex::zero(), Complex::new(-1.0, 0.0), Complex::new(2.0, 0.0)); let ldl = m.ldl().unwrap(); - - println!("{:}", &m); - println!("{:}", ldl.l_matrix()); - println!("{:}", ldl.d()); // Rebuild - let p = ldl.l_matrix() * ldl.d_matrix() * ldl.l_matrix().adjoint(); - - - println!("{:}", &p); + let p = ldl.l() * ldl.d() * ldl.l().adjoint(); assert!(relative_eq!(m, p, epsilon = 3.0e-12)); } @@ -32,10 +25,10 @@ fn ldl_partial() { Complex::zero(), Complex::zero(), Complex::zero(), Complex::zero(), Complex::zero(), Complex::new(2.0, 0.0)); - let ldl = m.lower_triangle().ldl().unwrap(); + let ldl = LDL::new_unchecked(m.lower_triangle()); // Rebuild - let p = ldl.l_matrix() * ldl.d_matrix() * ldl.l_matrix().adjoint(); + let p = ldl.l() * ldl.d() * ldl.l().adjoint(); assert!(relative_eq!(m, p, epsilon = 3.0e-12)); } @@ -50,7 +43,7 @@ fn ldl_lsqrtd() { let chol= m.cholesky().unwrap(); let ldl = m.ldl().unwrap(); - + assert!(relative_eq!(ldl.lsqrtd().unwrap(), chol.l(), epsilon = 3.0e-16)); } @@ -66,7 +59,7 @@ fn ldl_non_sym_panic() { let ldl = m.ldl().unwrap(); // Rebuild - let p = ldl.l_matrix() * ldl.d_matrix() * ldl.l_matrix().transpose(); + let p = ldl.l() * ldl.d() * ldl.l().transpose(); assert!(relative_eq!(m, p, epsilon = 3.0e-16)); } From 091d1a1a45d3073541ea30566d459f93e47a13c3 Mon Sep 17 00:00:00 2001 From: "Andreas J. Weiss" Date: Tue, 14 Apr 2026 11:12:32 -0400 Subject: [PATCH 7/8] ldl determinant, solve, solve_mut --- src/linalg/ldl.rs | 45 +++++++++++++++++++++++++++++++++++++++++---- tests/linalg/ldl.rs | 23 +++++++++++++++++++++-- 2 files changed, 62 insertions(+), 6 deletions(-) diff --git a/src/linalg/ldl.rs b/src/linalg/ldl.rs index dd3770516..8b228af56 100644 --- a/src/linalg/ldl.rs +++ b/src/linalg/ldl.rs @@ -4,9 +4,10 @@ use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; use crate::base::{Const, DefaultAllocator, OMatrix, OVector}; use crate::dimension::Dim; +use crate::{Matrix, Storage}; use simba::scalar::ComplexField; -/// The LDL / LDL^T factorization of a Hermitian matrix A = LDL^T where L is a lower unit-triangular matrix and D is diagonal matrix. +/// The LDL / LDL^T factorization of a Hermitian matrix A = LDL^T where L is a lower unit-triangular matrix and D is a diagonal matrix. /// /// This implementation referes to the strictly diagonal LDL algorithm. /// For the more general Bunch-Kaufman block-diagonal pivoting method we refer to the LBL factorization. @@ -41,6 +42,13 @@ where OMatrix::from_diagonal(&self.0.diagonal()) } + /// Computes the determinant of the decomposed matrix. + pub fn determinant(&self) -> T { + self.diagonal() + .iter() + .fold(T::one(), |acc, next| acc * next.clone()) + } + /// Returns the diagonal elements as a vector. #[must_use] pub fn diagonal(&self) -> OVector { @@ -114,8 +122,6 @@ where matrix[(i, j)] = l_ij / matrix[(j, j)].clone(); - dbg!(&matrix[(i, j)]); - // Zero out the upper triangular part. matrix[(j, i)] = T::zero(); } @@ -124,7 +130,7 @@ where Some(Self(matrix)) } - /// Computes the LDL / LDL^T factorization, allowing for zero pivots. If a zero pivot is detected the corresponding column & row is zero'd out. + /// Computes the LDL / LDL^T factorization, allowing for zero pivots. If a zero pivot is detected the corresponding column & row is filled with zeros. /// /// This factorization is only valid for singular matrices with a specific structure, e.g. singular matrices of the form A = xxT. pub fn new_unchecked(mut matrix: OMatrix) -> Self { @@ -163,4 +169,35 @@ where Self(matrix) } + + /// Solves the linear system A * x = b using this factorization. + pub fn solve(&self, b: &Matrix) -> Option> + where + S: Storage, + DefaultAllocator: Allocator, + { + let mut result = b.clone_owned(); + + if self.solve_mut(&mut result) { + Some(result) + } else { + None + } + } + + /// Solves the linear system A * x = b in place, overwriting `b` with the solution. + pub fn solve_mut(&self, b: &mut OMatrix) -> bool + where + DefaultAllocator: Allocator, + { + if self.determinant() == T::zero() { + return false; + } + + self.l().solve_lower_triangular_unchecked_mut(b); + self.d().solve_lower_triangular_unchecked_mut(b); + self.l().transpose().solve_upper_triangular_unchecked_mut(b); + + true + } } diff --git a/tests/linalg/ldl.rs b/tests/linalg/ldl.rs index 48da18cac..8399ff024 100644 --- a/tests/linalg/ldl.rs +++ b/tests/linalg/ldl.rs @@ -1,4 +1,4 @@ -use na::{linalg::LDL, Complex, Matrix3}; +use na::{linalg::LDL, Complex, Matrix2, Matrix3, Vector2, QR}; use num::Zero; #[test] @@ -15,11 +15,12 @@ fn ldl_simple() { let p = ldl.l() * ldl.d() * ldl.l().adjoint(); assert!(relative_eq!(m, p, epsilon = 3.0e-12)); + assert!(relative_eq!(m.determinant(), ldl.determinant(), epsilon = 3.0e-12)); } #[test] #[rustfmt::skip] -fn ldl_partial() { +fn ldl_unchecked() { let m = Matrix3::new( Complex::new(2.0, 0.0), Complex::zero(), Complex::zero(), Complex::zero(), Complex::zero(), Complex::zero(), @@ -31,6 +32,7 @@ fn ldl_partial() { let p = ldl.l() * ldl.d() * ldl.l().adjoint(); assert!(relative_eq!(m, p, epsilon = 3.0e-12)); + assert!(ldl.determinant() == Complex::zero()); } #[test] @@ -47,6 +49,23 @@ fn ldl_lsqrtd() { assert!(relative_eq!(ldl.lsqrtd().unwrap(), chol.l(), epsilon = 3.0e-16)); } +#[test] +#[rustfmt::skip] +fn ldl_solve() { + let m = Matrix2::new(8.0, 2.0, 2.0, 4.0); + let b = Vector2::new(5.0, 1.0); + + println!("{}", m); + + let ldl = m.ldl().unwrap(); + let chol = m.cholesky().unwrap(); + + let sol_ldl = ldl.solve(&b).unwrap(); + let sol_chol = chol.solve(&b); + + assert!(relative_eq!(sol_ldl, sol_chol, epsilon = 3.0e-16)); +} + #[test] #[should_panic] #[rustfmt::skip] From 9bdf4374568478622b7b0bb8f4bf1c40cc7121bd Mon Sep 17 00:00:00 2001 From: "Andreas J. Weiss" Date: Fri, 24 Apr 2026 15:29:51 -0400 Subject: [PATCH 8/8] fixed statement in the docs --- src/linalg/ldl.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linalg/ldl.rs b/src/linalg/ldl.rs index 8b228af56..26f3a71e2 100644 --- a/src/linalg/ldl.rs +++ b/src/linalg/ldl.rs @@ -132,7 +132,7 @@ where /// Computes the LDL / LDL^T factorization, allowing for zero pivots. If a zero pivot is detected the corresponding column & row is filled with zeros. /// - /// This factorization is only valid for singular matrices with a specific structure, e.g. singular matrices of the form A = xxT. + /// This factorization is only valid for singular matrices with a specific structure, e.g. covariance matrices. pub fn new_unchecked(mut matrix: OMatrix) -> Self { for j in 0..matrix.ncols() { let mut d_j: T = matrix[(j, j)].clone();