diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index 957d38ab1..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, + SymmetricEigen, SymmetricTridiagonal, LDL, LU, QR, SVD, U1, UDU, }; /// # 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: ComplexField, + 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..d029b71f4 --- /dev/null +++ b/src/linalg/ldl.rs @@ -0,0 +1,102 @@ +#[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::ComplexField; + +/// 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 hermitian c 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(); + } + } + + d[j] = d_j; + + for i in j..n { + let mut l_ij = p[(i, j)].clone(); + + for k in 0..j { + l_ij -= l[(j, k)].clone().conjugate() * l[(i, k)].clone() * d[k].clone(); + } + + if d[j] == T::zero() { + 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..6f3571882 --- /dev/null +++ b/tests/linalg/ldl.rs @@ -0,0 +1,108 @@ +use na::{Complex, Matrix3}; +use num::Zero; + +#[test] +#[rustfmt::skip] +fn ldl_simple() { + 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), + Complex::zero(), Complex::new(-1.0, 0.0), 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_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( + 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(); + + assert!(relative_eq!(ldl.cholesky_l(), chol.l(), 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;