From fdcd522956eb1a511b83f8a50624b0b0be9a285e Mon Sep 17 00:00:00 2001 From: Diogo Ribeiro Date: Fri, 19 Jun 2026 20:18:47 +0100 Subject: [PATCH 1/2] fix: Cholesky::new returns None for non-positive-definite complex matrices MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit For complex types, try_sqrt always succeeds (every complex number has a square root), so non-PD matrices were silently producing garbage results. Fix by calling try_sqrt on the real part of each diagonal pivot instead — matching the behaviour already correct for f64. Fixes #1536 --- src/linalg/cholesky.rs | 8 ++++++-- tests/linalg/cholesky.rs | 25 +++++++++++++++++++++++++ 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 4d667fbc1..5343b4010 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -235,10 +235,14 @@ where } let sqrt_denom = |v: T| { - if v.is_zero() { + let re = v.clone().real(); + let im = v.imaginary(); + + if re <= T::RealField::zero() || !im.is_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..8737e4a53 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -178,3 +178,28 @@ 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 with a non-real diagonal pivot 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()); +} From 79d9ae56589315b3a166c6f1ab3c23f1dfc89b96 Mon Sep 17 00:00:00 2001 From: Diogo Ribeiro Date: Fri, 19 Jun 2026 22:22:07 +0100 Subject: [PATCH 2/2] fix: discard roundoff imaginary part of complex Cholesky pivots The previous fix rejected any diagonal pivot with a non-zero imaginary part, but floating-point roundoff gives legitimate Hermitian positive-definite matrices tiny non-zero imaginary parts on their diagonal pivots, causing valid decompositions to fail. Take the real part of the pivot (the imaginary part is a roundoff artifact) and require it to be strictly positive instead. Also run cargo fmt. --- src/linalg/cholesky.rs | 12 ++++++++---- tests/linalg/cholesky.rs | 31 ++++++++++++++++++++++--------- 2 files changed, 30 insertions(+), 13 deletions(-) diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs index 5343b4010..9431f2aa5 100644 --- a/src/linalg/cholesky.rs +++ b/src/linalg/cholesky.rs @@ -235,10 +235,14 @@ where } let sqrt_denom = |v: T| { - let re = v.clone().real(); - let im = v.imaginary(); - - if re <= T::RealField::zero() || !im.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; } diff --git a/tests/linalg/cholesky.rs b/tests/linalg/cholesky.rs index 8737e4a53..ac7100bb1 100644 --- a/tests/linalg/cholesky.rs +++ b/tests/linalg/cholesky.rs @@ -190,16 +190,29 @@ fn cholesky_non_pd_complex_returns_none() { 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), - ]); + 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 with a non-real diagonal pivot 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), - ]); + // 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()); }