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..217171054 --- /dev/null +++ b/src/linalg/ldl.rs @@ -0,0 +1,121 @@ +#[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; + +/// 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 = "OMatrix: Serialize")) +)] +#[cfg_attr( + feature = "serde-serialize-no-std", + serde(bound(deserialize = "OMatrix: Deserialize<'de>")) +)] +#[derive(Clone, Debug)] +pub struct LDL(OMatrix) +where + DefaultAllocator: Allocator + Allocator; + +impl Copy for LDL +where + DefaultAllocator: Allocator + Allocator, + OMatrix: Copy, +{ +} + +impl LDL +where + DefaultAllocator: Allocator + Allocator, +{ + /// Returns the diagonal elements as a vector. + #[must_use] + pub fn d(&self) -> OVector { + self.0.diagonal() + } + + /// Returns the diagonal elements as a matrix. + #[must_use] + pub fn d_matrix(&self) -> OMatrix { + OMatrix::from_diagonal(&self.0.diagonal()) + } + + /// Returns the lower triangular matrix. + #[must_use] + pub fn l_matrix(&self) -> OMatrix { + let mut l = self.0.clone(); + + l.column_iter_mut() + .enumerate() + .for_each(|(idx, mut column)| { + column[idx] = T::one(); + }); + + l + } + + /// 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. + pub fn lsqrtd(&self) -> Option> { + let n_dim = self.0.shape_generic().1; + + 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()), + )); + + // 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) + } + } + + /// 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/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..d51630db4 --- /dev/null +++ b/tests/linalg/ldl.rs @@ -0,0 +1,115 @@ +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.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); + + 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_matrix() * ldl.d_matrix() * ldl.l_matrix().adjoint(); + + assert!(relative_eq!(m, p, epsilon = 3.0e-12)); +} + +#[test] +#[rustfmt::skip] +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), + 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.lsqrtd().unwrap(), 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_matrix() * ldl.d_matrix() * ldl.l_matrix().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;