Cholesky factorization for slightly non-Hermitian matrices

[caldwellshane]

I wrote a function to rotate the coordinates of a 2D MvNormal distribution, and in the course of testing it I discovered the Cholesky factorization throws an error depending on the input with bit-level sensitivity. Is this a bug? For my use case, the workaround of wrapping the input in Hermitian() doesn't work, because I can't pass a Hermitian to the MvNormal constructor, which is calling cholesky for me.

Example:

julia> VERSION
v"1.1.0"

julia> rot1 = [1 -1; 1 1] / sqrt(2)
2×2 Array{Float64,2}:
 0.707107  -0.707107
 0.707107   0.707107

julia> rot2 = [cos(π/4) -sin(π/4); sin(π/4) cos(π/4)]
2×2 Array{Float64,2}:
 0.707107  -0.707107
 0.707107   0.707107

julia> cholesky(rot1 * [25.0 0.0; 0.0 4.0] * rot1')
Cholesky{Float64,Array{Float64,2}}
U factor:
2×2 UpperTriangular{Float64,Array{Float64,2}}:
 3.80789  2.75744
  ⋅       2.62613

julia> cholesky(rot2 * [25.0 0.0; 0.0 4.0] * rot2')
ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.
Stacktrace:
 [1] checkpositivedefinite(::Int64) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/factorization.jl:11
 [2] #cholesky!#97(::Bool, ::Function, ::Array{Float64,2}, ::Val{false}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/cholesky.jl:182
 [3] #cholesky#101 at ./none:0 [inlined]
 [4] cholesky at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/cholesky.jl:275 [inlined] (repeats 2 times)
 [5] top-level scope at none:0

Related issues from Discourse:

• https://discourse.julialang.org/t/multivariate-normal-with-positive-semi-definite-covariance-matrix/3029
• https://discourse.julialang.org/t/is-this-a-bug-with-cholesky/16970
• https://discourse.julialang.org/t/numerical-stability-of-cholesky-factorization/20063

[stevengj]

This is not a problem of "numerical stability". The issue is, what should the cholesky function do if the input is only nearly Hermitian? Should it simply discard the non-Hermitian part, assuming it to be an error? If so, what should be the threshold?

This was discussed in JuliaLang/LinearAlgebra.jl#182, and the basic conclusion was that the user should explicitly indicate if the non-Hermitianness of a matrix should be ignored by doing cholesky(Hermitian(matrix)).

(The documentation of cholesky could probably be clearer about this.)

[stevengj]

Closing as a duplicate of JuliaLang/LinearAlgebra.jl#182.

[caldwellshane]

Thank you @stevengj.

For posterity: The fix to the problem that led me here is to cast to Symmetric in

MvNormal(μ, Symmetric(Σ))

rather than to Hermitian as I suggested in the OP. (Makes sense.) This solution is also found in the first Discourse issue, linked above.

[david-vicente]

(The documentation of cholesky could probably be clearer about this.)

Please do this. Just lost a bunch of hours trying to implement some Gaussian Processes code that just wouldn't work. I ended up re implementing in Python to find out that what I wanted was indeed possible and finally found this issue here. This could have been avoided if the documentation was a bit more clear :(

[stevengj]

@4aHxKzD, see #35320.