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)

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.

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.