Cholesky_factor_corr initializing with invalid values

Developers
1

Thanks to those who have answered my previous questions, it is great that there is an active and helpful community for Stan.

Today I have potentially what might be a bug with cholesky_factor_corr. Basically when the size of the correlation matrix gets large, e.g. 40 x 40, I run into issues with the initial values for Cholesky factor correlation parameters:

Rejecting initial value:
  Error evaluating the log probability at the initial value.
Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite.  last conditional variance is 6.66134e-015.  (in 'model3d5066d27a9a_0b961540b66c542856a20ec83237a5e1' at line 18)

I have a simple example in R below:

stan.code <- "
data {
    int K;
    int R;
    vector[K] draws[R];
    vector[K] sigma;
    vector[K] mu_0;
}
parameters {
    cholesky_factor_corr[K] L_omega;
}
transformed parameters {
    matrix[K, K] Sigma;
    Sigma = quad_form_diag(tcrossprod(L_omega), sigma);
}
model {
    draws ~ multi_normal(mu_0, Sigma);
}
"

n.draws <- 200
n.variables <- 40

draws <- matrix(rnorm(n.variables * n.draws),
                nrow = n.draws,
                ncol = n.variables)

stan.data <- list(
    K = n.variables,
    R = n.draws,
    draws = draws,
    sigma = rep(1, n.variables),
    mu_0 = rep(0, n.variables)
)

fit <- stan(model_code = stan.code,
            data = stan.data,
            iter = 100,
            chains = 1)

If n.variables is reduced to 20 or I set init = 0 then everything works fine, so it appears to be an issue with initialization of the Cholesky factor correlation parameter when the matrix is large. Right now my solution is to pass in an identity matrix as the initial value. Has anyone reported this before?

2

Yes, but it is not really a bug, just a numerical problem. You could also do init_r = some number less than the default of 2.

3

Thanks, but that would affect all parameters and the number to use needs to be found with trial and error and may not work with different data. Are there plans to fix this numerical problem? Also, what are downsides of setting init = 0 all the time, apart from not being able to randomize the starting point?

4

The unconstrained parameters are all initialized on [-2,2], which if anything is too big.

No, as it is not really a problem with the code but a fundamental limitation of how computers do math.

For this model, that is about it. If you do anything with unit_vectors, the sampling won’t start if init = 0.

1 Like
5

Sometimes these problems can be fixed.

What’s failing is the Eigen library impelemntation of Cholesky factorization for badly conditioned matrices. The bad conditioning comes up from our initializations. We don’t know how to initialize robustly in general as it depends on the scale of the problem, which we can’t really know in advance.

Our code’s not well set up to change initialization for a single variable type just for initialization. We could tinker with trying to make the transforms more robust, but we don’t really know how to do that—it’s an open problem.

2 Likes
6

Or the diagonal of the Cholesky factor underflows to zero and Eigen is doing the right thing by saying it is not positive definite.

7

Right. If the true value can’t be represented in floating point without underflowing, Eigen has no choice.

If it happens on an intermediate calculation, there’s then the issue of trying to fix it at the possible cost of speed for non-boundary cases and increased doc and maintenance burden.

8

Hi there, I’m facing the same issue when initializing a Kronecker Gaussian Process model, wonder if there has been any update to this open problem yet?

9

One approach is to try initialising using an identity matrix, that’s worked well for me in the past

2 Likes