Dealing with non-positive definite matrices in Cholesky decomposition: R vs. Julia

General Usage
1

Hi all,

While porting some R code into Julia, I have encountered many cases where cholesky in Julia returns an error because the matrix is not positive definite. However, in R, I do not run into these cases. I was wondering why this happens and if there is a workaround. For example, is there a tolerance I can adjust or some other workaround if this is due to float point error? Here is a MWE:

R

x = matrix(c(59.30480616243846, -174.3444811471312, -174.34448114713118, 512.7624424006956), 2, 2)
chol(x)

R Output

         [,1]       [,2]
 [1,] 7.700961 -22.639314
 [2,] 0.000000   0.473195

Julia

using LinearAlgebra
x = [59.30480616243846 -174.3444811471312; -174.34448114713118 512.7624424006956]
cholesky(x)

** Julia Output **

ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.

2

I forgot to note that I saw a tolerance value in the docs, but I am unsure what is a reasonable value to use if that is one way to workaround the problem.

3

If you don’t care about losing the differences between (i,j) and (j,i) entries you can wrap in Symmetric:

julia> x2 = Symmetric(x)
2×2 Symmetric{Float64, Matrix{Float64}}:
   59.3048  -174.344
 -174.344    512.762

julia> cholesky(x2)
Cholesky{Float64, Matrix{Float64}}
U factor:
2×2 UpperTriangular{Float64, Matrix{Float64}}:
 7.70096  -22.6393
  ⋅         0.473195
1 Like
4

Thank you. This might work for my case. I can also print out a warning or error if the degree of asymmetry is too high. In this case, it is very low:

julia> sum(abs.(x .- x'))
5.684341886080802e-14

I’m open to alternative ideas if this could be problematic for some reason.

5

Let me add that you can consider modifying the function that generates those data so that it automatically returns exactly symmetric matrices. Chances are that you are doing 2x the work somewhere, because you are computing those off-diagonal elements twice. I’m not sure if it is easy to do, but it’s something worth considering.

2 Likes