Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 10 additions & 2 deletions src/linalg/cholesky.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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() };
Expand Down
38 changes: 38 additions & 0 deletions tests/linalg/cholesky.rs
Original file line number Diff line number Diff line change
Expand Up @@ -178,3 +178,41 @@ macro_rules! gen_tests(

gen_tests!(complex, RandComplex<f64>);
gen_tests!(f64, RandScalar<f64>);

#[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<f64> 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());
}
Loading