How to check if a matrix is positive definite or not

[Jiang Du]

Hi,

I'm in a process of building a large spatial covariance matrix with lots of blocks. There is a correlation parameter (\gamma) for the off diagonal blocks which need to be constrained such that the resulting covariance matrix is positive definite. 

Unfortunately, there is no closed form to get the support of \gamma to ensure the positive definiteness. The support is between -1 and 1 with holes inside. What I did in R is during the sampling, if non-positive definiteness is found (catch the error message from a Cholesky decomposition), I will set the joint log density to the smallest number in log scale (i.e. zero density). Then that proposed \gamma will be declined and keep the previous \gamma. 

I'm interested in implementing my model in Rstan to check my own coding (it's so easy to make some small bugs when programming my sampling procedure with such a complicated model, and the convergence is not good so far). But I got stuck when dealing with theses issues.

1. I notice that I can calculate the eigenvalues of a matrix and then check if all positive( or simply the first one, since it's been ordered). But I'm wondering if there be some machine error going on if I'm using this method.

2. Also, if I find a good way to check positive definiteness, is there a way for me to set the joint density to 0 (let stan stay with the previous sample and move on)? Or do I have other ways to address my unknown support issue?

Thanks,

Jiang

[Krzysztof Sakrejda]

You _could_ check and then reject ("reject()"?) when it is not positive definite but that would kill sampling efficiency.  It's better to construct it so that it is positive definite but I don't think you give us enough detail to make any useful suggestions. Krzysztof

 

Thanks,

Jiang

[Ben Goodrich]

On Thursday, May 12, 2016 at 4:25:29 PM UTC-4, Jiang Du wrote:

Unfortunately, there is no closed form to get the support of \gamma to ensure the positive definiteness. The support is between -1 and 1 with holes inside.

In cases like these, the prior often does not integrate to a constant over the support. If it does integrate to a constant, then you could just declare the thing as a cov_matrix in transformed parameters and it will reject any indefinite ones (at considerable cost to the efficiency). If it doesn't integrate to a constant, then you are not drawing from a well-defined posterior distribution.