issues with Cholesky decomposition

Hello guys. I have a question with Cholesky decompostion. I am running a code for MCMC sampling in high dimensional linear regression. You know that this means I have to perform Cholesky decomposition many times. In some cases I get error by non-positive definite covariance matrix. However, when I rerun the code again with same dataset, the problem can disappear. I really have no clue why this happens. Any suggestions? Thanks!

[Walter Roberson]

It is possible to get a small bit of round-off error in building your matrix. For any given matrix, M, the way around it is:

M = (M+M.')/2;

Answer 1

[Christine Tobler]

This can happen if your matrix is close to symmetric positive semi-definite (meaning the smallest eigenvalue is around machine epsilon compared to the largest eigenvalue). The Cholesky function does not support symmetric positive semi-definite input.

If this doesn't happen every time you run the code, it must mean that the matrices you're passing to Cholesky are not exactly equal between runs - one of the other operations used must give slightly different results between runs. The function chol guarantees that if you call it twice with the exact same input within the same MATLAB session, it will return the exact same output (unless you're changing some very specific MATLAB settings in between).

How are you generating your covariance matrix? If it's by forming C = M'*M, you could instead compute the QR decomposition of M: M = Q*R, M'*M = R'*Q'*Q*R = R'*R (using that Q'*Q is the identity matrix for the QR decomposition). This is numerically more robust; note that if C is close to positive semi-definite, the last row(s) of R will likely be close to machine epsilon relative to the other entries of R. You may need to be careful about using R depending on the next step of your algorithm.

[Yu Bai]

Hi Christine,

Indeed matrices passing to Cholesky are not the same in each run, since this is in a MCMC algorithm covariance matrix is obtained from drawing from some distributions. Let me check your way to see it works. Thanks!

I run my code 674 times for forecast evaluation. Initially, there are 10 cases I get an error of not positive definite. After the second run, problems for 9 cases are solved and the last one is OK if I rerun the code for the third time.

Best,

Yu