Error in Chol.Default(Cxx):The Leading Minor of Order Is Not Positive Definite

Error in chol.default(Cxx) : the leading minor of order is not positive definite

The error you are seeing occurs when some of the eigenvectors of the matrix you are trying to operate on are not positive (typically they'll be zero, or below some very small threshold); this means, essentially, that your data are too noisy/small to estimate a full covariance matrix.

Regularizing means (approximately) adding a penalty term to push your estimates away from zero (in this case, pushing your matrices away from having non-positive eigenvectors). If your regularization parameters (lambda1, lambda2) are too small, then you'll get the error. Since your grid1 and grid2 sequences start from zero or very small values, rCCA will choke for these too-small values.

Try setting your grid1 and grid2 sequences to start at a larger value, e.g.

grid1 <- grid2 <- seq(0.05, 0.2, length=5)

Getting Error in chol.default(K) : the leading minor of order 5 is not positive definite with betareg

Fitting a beta distribution to the data in category 1 will be very challenging with three observations being essentially zero. With rounding to five digits: 0.00000, 0.00000, 0.00000, 0.00320, 0.00610, 0.01500. It's not clear to me whether this category should be modeled in the same way as the others.

In category 4 there is another observation that is numerically zero although the other observations are somewhat larger: 0.00000, 0.01000, 0.02500, 0.03100, 0.03400, 0.04100.

Omitting category 1 at least allows to estimate the model without numerical problems. Whether or not the asymptotic inference is a good approximation for two parameters from six observations per group is another question. The precision does not seem to be the same across groups though.

betareg(value ~ category | 1, data = df, subset = category != "c1")
## Call:
## betareg(formula = value ~ category | 1, data = df, subset = category != 
##     "c1")
## 
## Coefficients (mean model with logit link):
## (Intercept)   categoryc3   categoryc4   categoryc5  
##      0.2634      -2.2758      -4.4627      -1.0206  
## 
## Phi coefficients (precision model with log link):
## (Intercept)  
##       2.312  
betareg(value ~ category | category, data = df, subset = category != "c1")
## Call:
## betareg(formula = value ~ category | category, data = df, subset = category != 
##     "c1")
## 
## Coefficients (mean model with logit link):
## (Intercept)   categoryc3   categoryc4   categoryc5  
##      0.2566      -2.6676      -4.0601      -0.9784  
## 
## Phi coefficients (precision model with log link):
## (Intercept)   categoryc3   categoryc4   categoryc5  
##      2.0849       3.5619      -0.2308      -0.1376  

coda gelman.diag(): Error in chol.default(W) : the leading minor of order nn is not positive definite

In gelman.diag(), the MPSRF is calculated by:

> coda::gelman.diag <-
function (x, confidence = 0.95, transform = FALSE, autoburnin = FALSE,
        multivariate = TRUE)
{
#A lot of code ...

Niter <- niter(x)
Nchain <- nchain(x)
Nvar <- nvar(x)
#W is the variance matrix of the chain
S2 <- array(sapply(x, var, simplify = TRUE), dim = c(Nvar,
                                                   Nvar, Nchain))
W <- apply(S2, c(1, 2), mean)
#A lot of code ...

if (Nvar > 1 && multivariate) {
CW <- chol(W)
#This is max eigenvalue from W^-1*B.
emax <- eigen(
  backsolve(CW, t(backsolve(CW, B, transpose = TRUE)), transpose = TRUE),
  symmetric = TRUE, only.values = TRUE)$values[1]
}

I edited the gelman.diag() by removing the Cholesky decomposition that gave the error, and by adding W and B to the list to be returned. This allows me to see why W cannot undergo Cholesky decomposition.

my.gelman.diag <- function(x,
                           confidence = 0.95,
                           transform = FALSE, 
                           autoburnin = FALSE,
                           multivariate = TRUE
){
  x <- as.mcmc.list(x)
  if (nchain(x) < 2)
    stop("You need at least two chains")
  if (autoburnin && start(x) < end(x)/2)
    x <- window(x, start = end(x)/2 + 1)
  Niter <- niter(x)
  Nchain <- nchain(x)
  Nvar <- nvar(x)
  xnames <- varnames(x)
  if (transform)
    x <- gelman.transform(x)
  x <- lapply(x, as.matrix)
  S2 <- array(sapply(x, var, simplify = TRUE),
              dim = c(Nvar, Nvar, Nchain)
  )
  W <- apply(S2, c(1, 2), mean)
  xbar <- matrix(sapply(x, apply, 2, mean, simplify = TRUE),
                 nrow = Nvar, ncol = Nchain)
  B <- Niter * var(t(xbar))
  if (Nvar > 1 && multivariate) {  #ph-edits 
    # CW <- chol(W)
    #    #This is W^-1*B.
    # emax <- eigen(
    #  backsolve(CW, t(backsolve(CW, B, transpose = TRUE)), transpose = TRUE),
    # symmetric = TRUE, only.values = TRUE)$values[1]
    emax <- 1
    mpsrf <- sqrt((1 - 1/Niter) + (1 + 1/Nvar) * emax/Niter)
  }  else {
    mpsrf <- NULL
  }

  w <- diag(W)
  b <- diag(B)
  s2 <- matrix(apply(S2, 3, diag), nrow = Nvar, ncol = Nchain)
  muhat <- apply(xbar, 1, mean)
  var.w <- apply(s2, 1, var)/Nchain
  var.b <- (2 * b^2)/(Nchain - 1)
  cov.wb <- (Niter/Nchain) * diag(var(t(s2), t(xbar^2)) - 2 *
                                    muhat * var(t(s2), t(xbar)))
  V <- (Niter - 1) * w/Niter + (1 + 1/Nchain) * b/Niter
  var.V <- ((Niter - 1)^2 * var.w + (1 + 1/Nchain)^2 * var.b +
              2 * (Niter - 1) * (1 + 1/Nchain) * cov.wb)/Niter^2
  df.V <- (2 * V^2)/var.V
  df.adj <- (df.V + 3)/(df.V + 1)
  B.df <- Nchain - 1
  W.df <- (2 * w^2)/var.w
  R2.fixed <- (Niter - 1)/Niter
  R2.random <- (1 + 1/Nchain) * (1/Niter) * (b/w)
  R2.estimate <- R2.fixed + R2.random
  R2.upper <- R2.fixed + qf((1 + confidence)/2, B.df, W.df) *
    R2.random
  psrf <- cbind(sqrt(df.adj * R2.estimate), sqrt(df.adj * R2.upper))
  dimnames(psrf) <- list(xnames, c("Point est.", "Upper C.I."))
  out <- list(psrf = psrf, mpsrf = mpsrf, B = B, W = W) #added ph
  class(out) <- "gelman.diag"
  return( out )
}

Applying my.gelman.diag() to short.mcmc.list:

> l <- my.gelman.diag(short.mcmc.list, autoburnin = FALSE, multivariate = TRUE)
> W <- l$W  #Within-sequence variance
> B <- l$B  #Between-sequence variance

An investigation into W shows that W is indeed positive definite, but its eigen-values are near 0 and hence it fails.

> evals.W <- eigen(W, only.values = TRUE)$values
> min(evals.W)
[1] 1.980596e-16

Indeed, increasing the tolerance shows that W is indeed positive definite.

> matrixNormal::is.positive.definite(W, tol = 1e-18)
[1] TRUE

So in reality, W is near linear dependency...

> eval <- eigen(solve(W)%*%B, only.values = TRUE)$values[1]

Error in solve.default(W) :
system is computationally singular: reciprocal condition number = 7.10718e-21

So in fact, removing the last two columns of W makes it now linear independent/positive definite. This indicates that there are correlated parameters within the chain, and the number of parameters can be reduced.

> evals.W[198]
[1] 9.579362e-05
> matrixNormal::is.positive.definite(W[1:198, 1:198])
[1] TRUE
> chol.W <- chol(W)

W is the within-variance of Markov chain, a measure of how each value in the state is different from the mean. If W is near singular, the variance is small, and thus the chain does not vary much. It is a slow-moving chain. The user should investigate this using trace plots, and possibly reducing the number of parameters. The chain may also be too short; if the chain is longer, the values within the chain may be different enough so that W is linearly independent.

To avoid the function from crashing, I suggest to use purrr::possibly() to substitute a missing value instead of throwing an archaic error.

>  pgd <- purrr::possibly(coda::gelman.diag, list(mpsrf = NA), quiet = FALSE)
> pgd(short.mcmc.list, autoburnin = FALSE, multivariate = TRUE)
Error: the leading minor of order 199 is not positive definite
[1] NA

See Brooks and Gelman, 1992 for more information.

Reference:
Gelman, A and Rubin, DB (1992) Inference from iterative simulation using multiple sequences, Statistical Science, 7, 457-511.

3D PCA plotting error: the leading minor of order 3 is not positive definite

When I run the current version of the code, I get the error you saw here:

> ellipse_base <- ellipse3d(cov(BASE_for3d),level = 0.95, centre=center_base)
Error in chol.default(cov) : 
  the leading minor of order 3 is not positive definite

If I look at the eigenvalues of that matrix, I see that the message is correct:

> eigen(cov(BASE_for3d))
eigen() decomposition
$values
[1] 1.454977e+01 1.433775e+00 8.521084e-16

$vectors
           [,1]       [,2]       [,3]
[1,] 0.81230516  0.5213444  0.2614581
[2,] 0.04362183  0.3927276 -0.9186197
[3,] 0.58159906 -0.7576048 -0.2962727

Notice that the smallest eigenvalue is approximately zero.

The BASE_for3d object has 3 points in it. They all lie in a plane (as any 3 points must). Thus the covariance matrix for them is rank 2, not rank 3. You need to increase BASE to include at least 4 rows not all in a plane, not just 3. For example,

BASE <- c("SR1.4", "SR2.2", "CI1.4", "SR1.5")

This results in the following:

BASE_for3d <- for_plot3d(BASE)
center_base <- center4plot(BASE_for3d)
ellipse_base <- ellipse3d(cov(BASE_for3d),level = 0.95, centre=center_base)
shade3d(ellipse_base, alpha = 0.2)

---

Related Topics