Cholesky erroneously produces (wrong) results for non-positive-definite complex matrices #1536
Description
See the following example:
https://play.rust-lang.org/?version=stable&mode=debug&edition=2024&gist=03236050db62b06a5a5a35b5438b0bad
The matrix here is
┌ ┐
│ -4 0.000000000000000122 0.000000000000000122 0 0 0 │
│ 0.000000000000000122 -4 2 0 0 0 │
│ 0.000000000000000122 2 -4 0 0 0 │
│ 0 0 0 -4 0.000000000000000122 0.000000000000000122 │
│ 0 0 0 0.000000000000000122 -4 2 │
│ 0 0 0 0.000000000000000122 2 -4 │
└ ┘
which is not positive definite. When the datatype of the matrix is f64, Cholesky::new correctly returns None. However, if I provide the same matrix but where the entries have the datatype Complex<f64>, Cholesky::new gives a result, which it shouldn't - if a matrix is not positive definite over the reals it certainly isn't positive definite over the complex plane. Furthermore, this result is not correct, as if I calculate
┌ ┐
│ 4+0i -0.000000000000000122+0.00000000000000000000000000000001494069094959771i -0.000000000000000122+0.00000000000000000000000000000001494069094959771i 0+0i 0+0i 0+0i │
│ -0.000000000000000122-0.00000000000000000000000000000001494069094959771i 4+0i -2+0.00000000000000024492935982947064i 0+0i 0+0i 0+0i │
│ -0.000000000000000122-0.00000000000000000000000000000001494069094959771i -2-0.00000000000000024492935982947064i 6.000000000000001+0i 0+0i 0+0i 0+0i │
│ 0+0i 0+0i 0+0i 4+0i -0.000000000000000122+0.00000000000000000000000000000001494069094959771i -0.000000000000000122+0.00000000000000000000000000000001494069094959771i │
│ 0+0i 0+0i 0+0i -0.000000000000000122-0.00000000000000000000000000000001494069094959771i 4+0i -2+0.00000000000000024492935982947064i │
│ 0+0i 0+0i 0+0i -0.000000000000000122-0.00000000000000000000000000000001494069094959771i -2-0.00000000000000024492935982947064i 6.000000000000001+0i │
└ ┘
which is not my original matrix (note that ignoring the floating point issues in the imaginary parts, it isn't even the negative of my original matrix as the last entry on the diagonal is 6)
I believe the problem here is that the only way that we are checking for positive definiteness in the Cholesky implementation is by seeing whether the non-zero entries on the diagonal produce a square root - but for complex numbers this is always true!
This source says that the square roots need to particularly be real and strictly positive for a complex matrix, so adding that to the check may be the solution. That said, I can't seem to find a more rigorous source for this. Edit: it is exercise 1.4.65 of David S. Watkins' Fundamentals of Matrix Computations.