R coda "The leading minor of order 3 is not positive definite"

Question

Can someone explain to me what this error message means?

I have an MCMC sampling method called hitandrun (in unfinished package https://github.com/davidkane9/kmatching) that gives me a list of matrices which have columns as multivariate output of individual samples, here is an example of it in action:

> A = matrix(1, ncol = 3)
> b = 1
> ## gives me solutions of Ax = b (a.k.a x + y + z = 1)
> h = hitandrun(A,b, n=10, chains = 2)
> h
[[1]]
          [,1]      [,2]      [,3]      [,4]      [,5]       [,6]      [,7]      [,8]      [,9]
[1,] 0.1804431 0.3340590 0.4195820 0.2061222 0.3591085 0.09984353 0.6707110 0.3926639 0.1283919
[2,] 0.6135745 0.4256909 0.3619727 0.2918238 0.5057426 0.81919629 0.2368842 0.1178713 0.2666737
[3,] 0.2059824 0.2402501 0.2184453 0.5020541 0.1351489 0.08096018 0.0924048 0.4894647 0.6049344
         [,10]
[1,] 0.1322112
[2,] 0.4736057
[3,] 0.3941831

[[2]]
           [,1]      [,2]      [,3]      [,4]       [,5]       [,6]      [,7]       [,8]      [,9]
[1,] 0.32883534 0.1284182 0.1735151 0.2005726 0.94511422 0.61653717 0.5130324 0.33228224 0.2088865
[2,] 0.65868549 0.3066952 0.5182009 0.3065610 0.01214334 0.07007548 0.1191157 0.01137002 0.3311197
[3,] 0.01247917 0.5648866 0.3082840 0.4928664 0.04274244 0.31338735 0.3678519 0.65634774 0.4599938
          [,10]
[1,] 0.61412223
[2,] 0.32289039
[3,] 0.06298738

I wanted to see the Gelman-Rubin diagnostic of this data, to see how much I need to thin it, but when I put it in function, I get an obscure error, and I don't know what it means. Here is that:

> mclist = lapply(h, function(x) mcmc(t(x), thin = 5))
> gelman.diag(mclist)
Error in chol.default(W) : 
  the leading minor of order 1 is not positive definite

(I guess now it's order 1, but before it was order 3) Are there any coda experts on SO? I've tried to debug it, but it leads me to an internal function La_chol and I don't know what to do from there.

Answer 1

There seems to be something wrong with getting the multivariate estimate of the Gelman-Rubin diagnostic. Setting multivariate = FALSE fixes the problem and outputs a single variable estimate for each variable. However, most of my variables are correlated because of the nature of the problem I am trying to solve, so I wonder (and hope) that this is giving me an overestimate of the diagnostic.