Add LDL decomposition

src/linalg/decomposition.rs

@@ -1,7 +1,7 @@
 use crate::storage::Storage;
 use crate::{
     Allocator, Bidiagonal, Cholesky, ColPivQR, ComplexField, DefaultAllocator, Dim, DimDiff,
-    DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, LBLT, LU, Matrix, OMatrix, QR, RealField,
+    DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, LBLT, LDL, LU, Matrix, OMatrix, QR, RealField,
     SVD, Schur, SymmetricEigen, SymmetricTridiagonal, U1, UDU,
 };
 
@@ -294,6 +294,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 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`.
+    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,203 @@
+#[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 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 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.
+#[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 matrix.
+    #[must_use]
+    pub fn d(&self) -> OMatrix<T, D, D> {
+        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<T, D> {
+        self.0.diagonal()
+    }
+
+    /// Returns the lower triangular matrix.
+    #[must_use]
+    pub fn l(&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.
+    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()
+            * OMatrix::from_diagonal(&OVector::from_iterator_generic(
+                n_dim,
+                Const::<1>,
+                self.0.diagonal().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();
+                }
+            }
+
+            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();
+
+                for k in 0..j {
+                    l_ij -= matrix[(j, k)].clone().conjugate()
+                        * matrix[(i, k)].clone()
+                        * matrix[(k, k)].clone();
+                }
+
+                matrix[(i, j)] = l_ij / matrix[(j, j)].clone();
+
+                // 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 filled with zeros.
+    ///
+    /// This factorization is only valid for singular matrices with a specific structure, e.g. covariance matrices.
+    pub fn new_unchecked(mut matrix: OMatrix<T, D, D>) -> 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 / d_j.clone();
+                }
+
+                // Zero out the upper triangular part.
+                matrix[(j, i)] = T::zero();
+            }
+        }
+
+        Self(matrix)
+    }
+
+    /// Solves the linear system A * x = b using this factorization.
+    pub fn solve<M: Dim, S>(&self, b: &Matrix<T, D, M, S>) -> Option<OMatrix<T, D, M>>
+    where
+        S: Storage<T, D, M>,
+        DefaultAllocator: Allocator<D, M>,
+    {
+        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<M: Dim>(&self, b: &mut OMatrix<T, D, M>) -> bool
+    where
+        DefaultAllocator: Allocator<D, M>,
+    {
+        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
+    }
+}

src/linalg/mod.rs

@@ -17,7 +17,11 @@ pub mod givens;
 mod hessenberg;
 pub mod householder;
 mod inverse;
+<<<<<<< ldldecomp
+mod ldl;
+=======
 mod lblt;
+>>>>>>> main
 mod lu;
 mod permutation_sequence;
 mod pow;
@@ -40,7 +44,11 @@ pub use self::cholesky::*;
 pub use self::col_piv_qr::*;
 pub use self::full_piv_lu::*;
 pub use self::hessenberg::*;
+<<<<<<< ldldecomp
+pub use self::ldl::*;
+=======
 pub use self::lblt::*;
+>>>>>>> main
 pub use self::lu::*;
 pub use self::permutation_sequence::*;
 pub use self::qr::*;

tests/linalg/ldl.rs

@@ -0,0 +1,127 @@
+use na::{linalg::LDL, Complex, Matrix2, Matrix3, Vector2, QR};
+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();
+
+    // Rebuild
+    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_unchecked() {
+    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 = LDL::new_unchecked(m.lower_triangle());
+
+    // Rebuild
+    let p = ldl.l() * ldl.d() * ldl.l().adjoint();
+
+    assert!(relative_eq!(m, p, epsilon = 3.0e-12));
+    assert!(ldl.determinant() == Complex::zero());
+}
+
+#[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]
+#[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]
+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() * 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 lblt;
 mod lu;
 mod pow;