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.
This commit is contained in:
parent
5aa6033a3b
commit
cea3bdc8e5
|
@ -140,13 +140,13 @@ where
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Computes the determinant of the decomposed matrix.
|
/// Computes the determinant of the decomposed matrix.
|
||||||
pub fn determinant(&self) -> N {
|
pub fn determinant(&self) -> N::SimdRealField {
|
||||||
let dim = self.chol.nrows();
|
let dim = self.chol.nrows();
|
||||||
let mut prod_diag = N::one();
|
let mut prod_diag = N::one();
|
||||||
for i in 0..dim {
|
for i in 0..dim {
|
||||||
prod_diag *= unsafe { *self.chol.get_unchecked((i, i)) };
|
prod_diag *= unsafe { *self.chol.get_unchecked((i, i)) };
|
||||||
}
|
}
|
||||||
prod_diag * prod_diag
|
prod_diag.simd_modulus_squared()
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -7,6 +7,7 @@ macro_rules! gen_tests(
|
||||||
use na::dimension::{U4, Dynamic};
|
use na::dimension::{U4, Dynamic};
|
||||||
use na::{DMatrix, DVector, Matrix4x3, Vector4};
|
use na::{DMatrix, DVector, Matrix4x3, Vector4};
|
||||||
use rand::random;
|
use rand::random;
|
||||||
|
use simba::scalar::ComplexField;
|
||||||
#[allow(unused_imports)]
|
#[allow(unused_imports)]
|
||||||
use crate::core::helper::{RandScalar, RandComplex};
|
use crate::core::helper::{RandScalar, RandComplex};
|
||||||
|
|
||||||
|
@ -83,18 +84,20 @@ macro_rules! gen_tests(
|
||||||
fn cholesky_determinant(n in PROPTEST_MATRIX_DIM) {
|
fn cholesky_determinant(n in PROPTEST_MATRIX_DIM) {
|
||||||
let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
|
let m = RandomSDP::new(Dynamic::new(n), || random::<$scalar>().0).unwrap();
|
||||||
let lu_det = m.clone().lu().determinant();
|
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();
|
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]
|
#[test]
|
||||||
fn cholesky_determinant_static(_n in PROPTEST_MATRIX_DIM) {
|
fn cholesky_determinant_static(_n in PROPTEST_MATRIX_DIM) {
|
||||||
let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
|
let m = RandomSDP::new(U4, || random::<$scalar>().0).unwrap();
|
||||||
let lu_det = m.clone().lu().determinant();
|
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();
|
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]
|
#[test]
|
||||||
|
|
Loading…
Reference in New Issue