From cea3bdc8e5fe1569eac3544888f613dd470396bf Mon Sep 17 00:00:00 2001 From: Jim Turner Date: Fri, 9 Apr 2021 13:55:40 -0400 Subject: [PATCH] Fix Cholesky::determinant for Complex elements The previous implementation was correct only for real elements. The Cholesky decomposition is `L L^H`, so the determinant is `det(L) * det(L^H)`. Since `L` is a triangular matrix, `det(L)` is the product of the diagonal elements of `L`. Since `L^H` is triangular and its diagonal elements are the conjugates of the diagonal elements of `L`, `det(L^H)` is the conjugate of `det(L)`. So, the overall determinant is the product of the diagonal elements of `L` times its conjugate. --- src/linalg/cholesky.rs | 4 ++-- tests/linalg/cholesky.rs | 7 +++++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 6bfbeb81..82f3ad86 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -140,13 +140,13 @@ where } /// Computes the determinant of the decomposed matrix. - pub fn determinant(&self) -> N { + pub fn determinant(&self) -> N::SimdRealField { let dim = self.chol.nrows(); let mut prod_diag = N::one(); for i in 0..dim { prod_diag *= unsafe { *self.chol.get_unchecked((i, i)) }; } - prod_diag * prod_diag + prod_diag.simd_modulus_squared() } } diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index 99893bd4..df69b8b1 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -7,6 +7,7 @@ macro_rules! gen_tests( use na::dimension::{U4, Dynamic}; use na::{DMatrix, DVector, Matrix4x3, Vector4}; use rand::random; + use simba::scalar::ComplexField; #[allow(unused_imports)] use crate::core::helper::{RandScalar, RandComplex}; @@ -83,18 +84,20 @@ macro_rules! gen_tests( fn cholesky_determinant(n in PROPTEST_MATRIX_DIM) { let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap(); let lu_det = m.clone().lu().determinant(); + assert_relative_eq!(lu_det.imaginary(), 0., epsilon = 1.0e-7); let chol_det = m.cholesky().unwrap().determinant(); - prop_assert!(relative_eq!(lu_det, chol_det, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(lu_det.real(), chol_det, epsilon = 1.0e-7)); } #[test] fn cholesky_determinant_static(_n in PROPTEST_MATRIX_DIM) { let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap(); let lu_det = m.clone().lu().determinant(); + assert_relative_eq!(lu_det.imaginary(), 0., epsilon = 1.0e-7); let chol_det = m.cholesky().unwrap().determinant(); - prop_assert!(relative_eq!(lu_det, chol_det, epsilon = 1.0e-7)); + prop_assert!(relative_eq!(lu_det.real(), chol_det, epsilon = 1.0e-7)); } #[test]