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,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<T, D>: Serialize, OMatrix<T, D, D>: Serialize"))
+)]
+#[cfg_attr(
+    feature = "serde-serialize-no-std",
+    serde(bound(
+        deserialize = "OVector<T, D>: Deserialize<'de>, OMatrix<T, D, D>: Deserialize<'de>"
+    ))
+)]
+#[derive(Clone, Debug)]
+pub struct LDL<T: ComplexField, D: Dim>
+where
+    DefaultAllocator: Allocator<D> + Allocator<D, D>,
+{
+    /// The lower triangular matrix resulting from the factorization
+    pub l: OMatrix<T, D, D>,
+    /// The diagonal matrix resulting from the factorization
+    pub d: OVector<T, D>,
+}
+
+impl<T: ComplexField, D: Dim> Copy for LDL<T, D>
+where
+    DefaultAllocator: Allocator<D> + Allocator<D, D>,
+    OVector<T, D>: Copy,
+    OMatrix<T, D, D>: Copy,
+{
+}
+
+impl<T: ComplexField, D: Dim> LDL<T, D>
+where
+    DefaultAllocator: Allocator<D> + Allocator<D, D>,
+{
+    /// 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<T, D, D>) -> Option<Self> {
+        let n = p.ncols();
+
+        let n_dim = p.shape_generic().1;
+
+        let mut d = OVector::<T, D>::zeros_generic(n_dim, Const::<1>);
+        let mut l = OMatrix::<T, D, D>::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<T, D, D> {
+        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<T, D, D> {
+        OMatrix::from_diagonal(&self.d)
+    }
+}

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

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;