diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 4d667fbc1..9431f2aa5 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -235,10 +235,18 @@ where } let sqrt_denom = |v: T| { - if v.is_zero() { + // The diagonal pivots of a Hermitian positive-definite matrix + // are real and strictly positive. For complex scalars `try_sqrt` + // always succeeds, so we explicitly take the real part (the + // imaginary part is a roundoff artifact) and require it to be + // positive — otherwise the matrix is not positive-definite. + let re = v.real(); + + if re <= T::RealField::zero() { return None; } - v.try_sqrt() + + re.try_sqrt().map(T::from_real) }; let diag = unsafe { matrix.get_unchecked((j, j)).clone() }; diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index bd135fb6a..ac7100bb1 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -178,3 +178,41 @@ macro_rules! gen_tests( gen_tests!(complex, RandComplex); gen_tests!(f64, RandScalar); + +#[test] +fn cholesky_non_pd_complex_returns_none() { + // Regression test for https://github.com/dimforge/nalgebra/issues/1536. + // A non-positive-definite matrix with Complex entries must return None, + // just as it does for f64 entries. The diagonal pivot elements during + // Cholesky factorization must be real and strictly positive; for complex + // types, try_sqrt always succeeds, so positivity must be checked explicitly. + use na::DMatrix; + use num_complex::Complex; + + // This 2x2 negative-definite matrix is not positive definite. + let m = DMatrix::from_row_slice( + 2, + 2, + &[ + Complex::new(-4.0_f64, 0.0), + Complex::new(0.0, 0.0), + Complex::new(0.0, 0.0), + Complex::new(-4.0_f64, 0.0), + ], + ); + assert!(na::Cholesky::new(m).is_none()); + + // A matrix whose diagonal pivot has a non-positive real part (here a purely + // imaginary entry, real part 0) is also not positive definite. + let m2 = DMatrix::from_row_slice( + 2, + 2, + &[ + Complex::new(0.0_f64, 1.0), + Complex::new(0.0, 0.0), + Complex::new(0.0, 0.0), + Complex::new(1.0_f64, 0.0), + ], + ); + assert!(na::Cholesky::new(m2).is_none()); +}