Add LDL decomposition

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<T: ComplexField, D: Dim, S: Storage<T, D, D>> Matrix<T, D, D, S> {
         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<LDL<T, D>>
+    where
+        T: ComplexField,
+        DefaultAllocator: Allocator<D> + Allocator<D, D>,
+    {
+        LDL::new(self.into_owned())
+    }
+
     /// Computes the Schur decomposition of a square matrix.
     pub fn schur(self) -> Schur<T, D>
     where

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<T, D, D>: Serialize"))
+)]
+#[cfg_attr(
+    feature = "serde-serialize-no-std",
+    serde(bound(deserialize = "OMatrix<T, D, D>: Deserialize<'de>"))
+)]
+#[derive(Clone, Debug)]
+pub struct LDL<T: ComplexField, D: Dim>(OMatrix<T, D, D>)
+where
+    DefaultAllocator: Allocator<D> + Allocator<D, D>;
+
+impl<T: ComplexField, D: Dim> Copy for LDL<T, D>
+where
+    DefaultAllocator: Allocator<D> + Allocator<D, D>,
+    OMatrix<T, D, D>: Copy,
+{
+}
+
+impl<T: ComplexField, D: Dim> LDL<T, D>
+where
+    DefaultAllocator: Allocator<D> + Allocator<D, D>,
+{
+    /// Returns the diagonal elements as a vector.
+    #[must_use]
+    pub fn d(&self) -> OVector<T, D> {
+        self.0.diagonal()
+    }
+
+    /// Returns the diagonal elements as a matrix.
+    #[must_use]
+    pub fn d_matrix(&self) -> OMatrix<T, D, D> {
+        OMatrix::from_diagonal(&self.0.diagonal())
+    }
+
+    /// Returns the lower triangular matrix.
+    #[must_use]
+    pub fn l_matrix(&self) -> OMatrix<T, D, D> {
+        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<OMatrix<T, D, D>> {
+        let n_dim = self.0.shape_generic().1;
+
+        let lsqrtd: crate::Matrix<T, D, D, <DefaultAllocator as Allocator<D, D>>::Buffer<T>> = &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<T, D, D>) -> Option<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;
+
+            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))
+    }
+}

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::*;

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);
+}

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;