LKJ prior throwing positivity exceptions

Modeling
1

I am trying to use an LKJ prior. I have simplified my model down to the following. It’s so simple that I feel like I must be doing something wrong, yet it does not work for me:

parameters {
  cholesky_factor_corr[2] L;
}
model {
  L ~ lkj_corr_cholesky(1.0);
}

I’m calling this lkj_test.stan. I run this through CmdStanPy. Again, after simplifying things as much as possible, here is exactly what I run:

import cmdstanpy

model = cmdstanpy.CmdStanModel(stan_file="lkj_test.stan")
fit = model.sample({}, chains=1)
print(fit.summary())

The result is:

16:09:46 - cmdstanpy - INFO - CmdStan start processing
chain 1: 100%|█████████████████████████████████████████████████████████████████████| 2000/2000 [00:00<00:00, 134301.53it/s, (Sampling completed)]
                                                                                                                                                 
16:09:46 - cmdstanpy - INFO - CmdStan done processing.
16:09:46 - cmdstanpy - WARNING - Non-fatal error during sampling:
Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in 'lkj_test.stan', line 5, column 2 to column 29)
        Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in 'lkj_test.stan', line 5, column 2 to column 29)
        Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in 'lkj_test.stan', line 5, column 2 to column 29)
        Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in 'lkj_test.stan', line 5, column 2 to column 29)
Consider re-running with show_console=True if the above output is unclear!
            Mean      MCSE    StdDev       MAD        5%       50%       95%  ESS_bulk  ESS_tail  ESS_bulk/s    R_hat
lp__   -0.609607  0.047758  0.855861  0.365846 -2.327080 -0.260325 -0.001084   394.434   410.966     98608.5  1.00091
L[1,1]  1.000000       NaN  0.000000  0.000000  1.000000  1.000000  1.000000       NaN       NaN         NaN      NaN
L[1,2]  0.000000       NaN  0.000000  0.000000  0.000000  0.000000  0.000000       NaN       NaN         NaN      NaN
L[2,1] -0.022324  0.031549  0.572107  0.707185 -0.912745 -0.009924  0.890539   301.995   331.944     75498.8  1.00651
L[2,2]  0.788447  0.011914  0.225666  0.169100  0.312378  0.877953  0.999458   394.434   410.966     98608.5  1.00091

While I don’t see anything wrong with the final L, the exceptions are worrying. I looked at the samples, and they looked like valid Cholesky factorizations of correlation matrices; but anything invalid should have been rejected, so I don’t think that’s any evidence of anything. When I run my actual model, I get the same exception many dozens of times.

As an experiment, I also tried

parameters {
  corr_matrix[2] M;
}
model {
  M ~ lkj_corr(1.0);
}

Calling this lkj_test2.stan and running it precisely as before results in:

16:21:25 - cmdstanpy - INFO - CmdStan start processing
16:24:32 - cmdstanpy - INFO - CmdStan start processing
chain 1: 100%|██████████████████████████████████████████████████████████████████████| 2000/2000 [00:00<00:00, 78904.82it/s, (Sampling completed)]
                                                                                                                                                 
16:24:32 - cmdstanpy - INFO - CmdStan done processing.
16:24:32 - cmdstanpy - WARNING - Non-fatal error during sampling:
Exception: lkj_corr_lpdf: Correlation matrix is not positive definite. (in 'lkj_test2.stan', line 5, column 2 to column 20)
        Exception: lkj_corr_lpdf: Correlation matrix is not positive definite. (in 'lkj_test2.stan', line 5, column 2 to column 20)
        Exception: lkj_corr_lpdf: Correlation matrix is not positive definite. (in 'lkj_test2.stan', line 5, column 2 to column 20)
        Exception: lkj_corr_lpdf: Correlation matrix is not positive definite. (in 'lkj_test2.stan', line 5, column 2 to column 20)
Consider re-running with show_console=True if the above output is unclear!
            Mean      MCSE    StdDev       MAD        5%       50%       95%  ESS_bulk  ESS_tail  ESS_bulk/s    R_hat
lp__   -0.663346  0.082399  1.006740  0.414771 -2.299060 -0.307206 -0.003113   425.331   212.841     53166.3  1.00199
M[1,1]  1.000000       NaN  0.000000  0.000000  1.000000  1.000000  1.000000       NaN       NaN         NaN      NaN
M[1,2] -0.005443  0.031301  0.587459  0.766037 -0.930505 -0.008561  0.880960   265.173   212.546     33146.6  1.00977
M[2,1] -0.005443  0.031301  0.587459  0.766037 -0.930505 -0.008561  0.880960   265.173   212.546     33146.6  1.00977
M[2,2]  1.000000       NaN  0.000000  0.000000  1.000000  1.000000  1.000000       NaN       NaN         NaN      NaN

Again, the final M and the individual samples look okay to me, but the exceptions are worrying. Also again, my actual model produces the same error, but many more times.

My question is: Am I doing something wrong? Or is this just the expected behavior when sampling from the LKJ distribution? If it is the expected behavior, is there some way I can suppress these exceptions (but not others)?

I am using Python 3.14.4 and CmdStanPy 1.3.0.

2

I don’t have a lot of experience with Cholesky matrices, but I have seen similar messages pop up when I invoke the Gamma distribution. There, Stan insists on y > 0 for \Gamma( y | \alpha, \beta ) , so I usually define my inputs with something like real<lower=0> y;. Placing hard constraints on a variable can cause havok with MCMC algorithms, due to tonnes of rejected samples; Stan tries to solve this by internally sampling an unconstrained value then transforming it into a constrained one, adjusting the gradient appropriately. For positive-only values, that transform is y = \exp( \dot y ). Problem solved, right?

transformed data {
   print( "exp( -747 ) = ", exp( -747 ), "." );
}
// STDOUT: exp( -747 ) = 0.

Alas, floating point values only have a certain amount of precision associated with them, so in reality the hard constraint has merely shifted location. That’s still a huge improvement, most of the time. Stan’s default initialization routine picks random unconstrained values between -2 and 2, the gradient steers it away from those problematic negative values, and you never see an exception.

So if you are seeing an exception, NUTS must be drifting into the Danger Zone. Unfortunately, there are quite a few potential reasons for this: you may have a broken model that favours values near zero for certain constrained variables, you may have a poorly-constrained model that doesn’t prevent the sampler from drifting, you may have not have enough data to steer the sampler away, you may be running across a corner case within Stan’s internal routines, or it may be some combination of all of the above. Even if any of the above apply, that doesn’t necessarily mean your samples are unrepresentative. Most MCMC algorithms reject far more than a dozen samples during their runtime, and yet we still consider the sample to be representative of the true distribution.

There’s really no simple solution here, you’ve got to roll up your sleeves and figure out the reason(s) yourself. One thing that may help: if I add the following to your toy example …

transformed parameters {
    if( L[1,1] == 0. )
        reject( "ERROR: invalid value in L[1,1]! Unconstrained value = ", lower_bound_unconstrain( L[1,1], 0. ), "." );
    if( L[2,2] == 0. )
        reject( "ERROR: invalid value in L[2,2]! Unconstrained value = ", lower_bound_unconstrain( L[2,2], 0. ), "." );
    if( L[2,1] == 0. )
        reject( "ERROR: invalid value in L[2,1]! Unconstrained value = ", lower_bound_unconstrain( L[2,1], 0. ), "." );
}

… I get different error messages.

Exception: ERROR: invalid value in L[2,2]! Unconstrained value = -inf.
Exception: ERROR: invalid value in L[2,2]! Unconstrained value = -inf.
Exception: ERROR: invalid value in L[2,2]! Unconstrained value = -inf.
Exception: ERROR: invalid value in L[2,2]! Unconstrained value = -inf.

The ERRORs forL[1,1] and L[2,1]never pop up. I figure this could be Stan constraining the value ofL[2,2]based on the contents of other elements in the matrix, and with no data or prior to steer it that occasionally leads to problematic values. The -inf is almost certainly misleading here, the underlying code is probably just \log\left( \exp\left( \dot y \right) \right) without compensation for floating point rounding. Still, it was worth showing off lower_bound_unconstrain().

3

I’m not sure that explains the problem; but I tried digging through the source, and the behavior I’m seeing looks impossible, so I must be misunderstanding something.

I found ‎lkj_corr_cholesky_rng in math/stan/math/prim/prob/lkj_corr_cholesky_rng.hpp. When K = 2, this creates an array CPCs with only a single entry. It sets alpha = eta + 0.5 and then runs a pair of nested loops. Both loops execute only once. In the only loop iteration, alpha gets updated to eta and CPCs[0] is set to 2.0 * beta_rng(alpha, alpha, rng) - 1.0. For my choice of \eta = 1, this is a \operatorname{Beta}(1, 1)-distributed random variable shifted and scaled to [-1, 1], i.e., it’s uniform on [-1, 1].

As the last step of ‎lkj_corr_cholesky_rng, CPCs is sent to read_corr_L in stan/math/prim/fun/read_corr_L.hpp. This creates an array acc with a single entry and a 2\times 2 matrix L. It first sets L to zero. Then, if I’m reading this right, it sets (using C++ indexing) L[0,0] = 1.0, L[1,0] = CPCs[0], acc[0] = 1.0 - CPCs[0].square(), and L[1,1] = sqrt(acc[0]). That is, if CPCs[0] is c, it produces the matrix L = \begin{pmatrix} 1 & 0 \\ c & \sqrt{1 - c^2} \end{pmatrix}. This is a valid Cholesky factorization of a correlation matrix; indeed, it’s easy to calculate LL^T = \begin{pmatrix} 1 & c \\ c & 1 \end{pmatrix}.

The error messages your suggestion yields seem informative with regards to what variable is causing the exceptions, but not with regards to why. The -infs appearing in your error messages look erroneous to me. Ultimately, lower_bound_unconstrain calls lb_free in stan/math/prim/constraint/lb_free.hpp. This function assumes that it’s been passed exponentiated values, so it takes a log. The fact that it’s reporting -inf means that it’s actually complaining about L[2,2] being zero.

The exception I’m seeing is coming from lkj_corr_cholesky_lpdf in stan/math/prim/prob/lkj_corr_cholesky_lpdf.hpp. Early on, this calls check_cholesky_factor in stan/math/prim/err/check_cholesky_factor.hpp. The check that fails is the last one in the function, check_positive applied to the diagonal. And that function is in stan/math/prim/err/check_positive.hpp, and all it does is check that what it’s passed is (strictly) greater than zero.

L[2,2] is a square root of a non-negative number. Mathematically, it’s almost surely positive. Numerically, there’s some tiny but positive probability that c is 1.0 and that we get zero exactly, but that’s extraordinarily unlikely unless \eta \ll 1. Which it’s not. That leaves me still confused: I can’t see how this is happening, and yet it is.

4

This is almost surely an overflow issue where the value is overflowing to 1 in which causes the diagonal to be 0 and triggers the warning. As @kyle.hofmann noted the mathematics say this is impossible but with finite-precision computing this happens.

In the original model, the dimension of the correlation matrix is only 2. I almost never use the built-in Cholesky factor of correlation matrices for this. Since you can easily parameterize it with 1 parameter yourself as

real<lower=-1, upper=1> phi;

We, the Stan devs, need to update the parameterization to something like Updated cholesky corr parameterization testing which shouldn’t have this issue.

5

I found ‎lkj_corr_cholesky_rng in math/stan/math/prim/prob/lkj_corr_cholesky_rng.hpp.

I’m not sure why you’re looking there. The PRNG sampling routines might be relevant during initialization, but if check_positive() on line 15 was the trigger then adding a transformed parameters section wouldn’t have changed the exception message. Line 22 in lkj_corr_cholesky_lpdf.hpp is more likely.

The -infs appearing in your error messages look erroneous to me. Ultimately, lower_bound_unconstrain calls lb_free in stan/math/prim/constraint/lb_free.hpp. This function assumes that it’s been passed exponentiated values, so it takes a log. The fact that it’s reporting -inf means that it’s actually complaining about L[2,2] being zero.

Yep, I said as much here:

The -inf is almost certainly misleading here, the underlying code is probably just log⁡(exp⁡(˙𝑦)) without compensation for floating point rounding.

While that might have been a bit too subtle, the conditional was ( L[2,2] == 0. ) .

That leaves me still confused: I can’t see how this is happening, and yet it is.

Ah, my inexperience with Cholesky correlation matrices is showing. As you and @spinkney pointed out, L[2,1] must be between -1 and 1. Thus Stan uses this transformation to go from unconstrained to constrained values …

y = \left( 1 + \exp(-\dot y) \right)^{-1}

… which has precision issues on two sides. To make what @spinkney said a bit more concrete:

    print( "inv_logit( -747 ) = ", inv_logit( -747 ), "." );  // STDOUT: inv_logit( -747 ) = 0.
    print( "inv_logit( 747 ) = ", inv_logit( 747 ), "." );    // STDOUT: inv_logit( 747 ) = 1.

Again, whether that’s an issue or not depends on your model. Does it require more precision than an IEEE double-precision float gives you? If so, then you’ve either got to get creative or ditch Stan.

6

I assumed that it would be relevant, but I didn’t put the pieces together to see that it wasn’t related to where the exception is actually being thrown. Sorry about that; this is my first time using Stan.

The transformation between unconstrained and constrained seems to happen at stan/math/prim/constraint/cholesky_corr_constrain.hpp#L17. This calls corr_constrain from stan/math/prim/constraint/corr_constrain.hpp#L27, and corr_constrain uses \tanh as its transformation. So not quite the logit and its inverse (it can’t be those because they map between [0, 1] and \mathbf{R}, while this function needs to map between [-1, 1] and \mathbf{R}), but different only by some shifting and scaling. And in standard double precision arithmetic, \tanh(19) = 1.0, so we don’t have to wander very far before we’re in trouble.

I’m still rather concerned that I get these exceptions every time even in the trivial case of my original post. My actual data has non-trivial correlation (which is a priori expected; it’s one of those “when people like something they do it more” phenomena), so the sampler is going to gravitate towards larger values of the correlation coefficient, and I guess that makes numerical issues more likely. My plan, after I figured this out, was to ask about translating the LKJ distribution so that I could put its peak somewhere other than the identity. But I’m starting to think it might be safer to switch to a Wishart or maybe inverse Wishart distribution. I don’t have a strong preference among them, and for this application, I suspect they’ll all work fine in the end.

Okay, that’s good to know. I didn’t know whether parameterizing by hand like that would be liable to cause issues elsewhere. Thanks!

7

Welcome! This is a pretty ambitious place to start.

This shouldn’t be a problem unless the correlations get close to 1 or the variances get close to 0. More information will tend to keep the results away from the boundaries as the posteriors concentrate away from uniform.

No. Yes.

The rejections will cause duplicated states in the sampler, but that’s not the only cause.

We intentionally label these as “Non-fatal error” because it’s recoverable. What’s usually going on is that early in sampling you will get these exceptiosn or you will get them occassionally during sampling. They can indicate there might be bias if they’re all due to the same cause and show up as divergences in the same region of the posterior.

Negative values themselves are necessary in order to get values less than 1. The problem starts when the absolute value of unconstrained variable exceeds 350 or so and starts underflowing 0 (negative) or overflowing to +infinity (positive).

This isn’t necessarily a broken model so much as our sampler not being up to the job. @andrewgelman likes to recommend things like half normal priors on scale parameters in hierarchical models. If the data is consistent with complete pooling, this can run into numerical difficulties compared to a “zero avoiding” prior like lognormal. So while Bayes wants a really small value in these cases, the exp() transform can be a problem.

In another topic here, @aseyboldt posted a more stable altnerative he calls “bisoftplus”. Like we did for sum-to-zero and simplex, we may wind up changing our default built-in transform as a result.

If L[2,2] is bad, you wont see messages for L[2,1] because rejections stop where they’re at.

Right, you take the correlation phi, and the scales sigma[1] and sigma[2] and create

parameters {
  vector<lower=0>[2] sigma;
  real<lower=-1, upper=1> phi;
}
transformed parameters {
  matrix[2, 2] Sigma;
  {
    vector[2] var = sigma^2;
    real cov = sigma[1] * sigma[2] * phi;
    Sigma = [[var[1], cov], [cov, var[2]];
  }
}
8

It’s the ambitious stuff that gives me enough of a kick to learn something new. But the learning curve has been more gentle than I expected!

Okay, great. Thanks, that helps.

Oh, that’s quite interesting. I like the fact that it’s basically linear for positive x and doesn’t go to zero too quickly for negative x. I looked for a better parameterization, but I didn’t find one. (My attempt was: The minimal polynomial of y = x + \sqrt{x^2 + 1} is y^2 - 2xy - 1. Over \mathbf{C} it defines a Riemann surface which double covers the sphere and is branched over \pm i. It’s a rational algebraic curve, so it has a rational parameterization. For example, fix a point on the curve like (0, 1) and consider the line of slope t through (0, 1), i.e., y = tx + 1; we can solve this pair of equations to express any one of x, y, or t in terms of the other two. But when we express y in terms of t, we find that y(t) is a parabola, and that’s not a helpful way to parameterize the positive reals for this application. And when we express x in terms of t, half the slopes correspond to the wrong branch, i.e,. to y = x - \sqrt{x^2 + 1}, so we end up with half the slopes being disallowed. That’s kinda like where we started with half the real numbers being disallowed, so it isn’t an improvement either. Consequently this idea seems to make things no better.)

In any case, I am busily sampling as we speak. A test run on a subset of my data produced good-looking results, so I’m feeling pretty good right now. Thanks again for your help!

9

I’m afraid I don’t know what minimal polynomials are or where t came from.

From a change-of-variables perspective, this should be fine. Let me be super pedantic about what we’re proposing.

In case anyone wants to implement, here are the details in the style of the Reference Manual chapter on transforms. Suppose we have a variable Y \in (0, \infty) and define

\qquad X = f(Y) = (Y^2 - 1) / (2 \cdot Y).

The inverse is

\qquad Y = f^{-1}(X) + \sqrt{X^2 + 1}.

Given the density of the constrained Y, p_Y(y), the unconstrained density of X is

\qquad p_X(x) = p_Y(f^{-1}(x)) \cdot \left| \nabla f^{-1}(x) \right|,

which reduces to

\qquad p_X(x) = p_Y\!\left(x + \sqrt{x^2 + 1}\right) \cdot \left(1 + \dfrac{x}{\sqrt{x^2 + 1}}\right).

10

Ah, well, you can safely ignore what I said! I was hoping to find something that behaved similarly but had a different formula. There is nothing wrong with the given one, but I thought I had an idea, so I worked it out. The idea failed—the result didn’t solve the original problem of mapping between (0, \infty) and \mathbf{R}. So while I found it to be a fun exercise, it’s not something anyone else ought to spend serious time with!

(But in case you’re curious: The minimal polynomial of x + \sqrt{x^2 + 1} is the smallest degree polynomial with leading coefficient 1 such that plugging in x + \sqrt{x^2 + 1} gives zero. To write down that polynomial, we introduce a new variable y. The polynomial y^2 - 2xy - 1 is quadratic in y, has leading coefficient 1, and if we substitute y = x + \sqrt{x^2 + 1}, we get zero. There’s no degree 1 polynomial with this property because of the square root, so y^2 - 2xy - 1 is the minimal polynomial of x + \sqrt{x^2 + 1}. If you’re familiar with field theory and Galois theory, this is the minimal polynomial of x + \sqrt{x^2 + 1} over the field \mathbf{R}(x) of rational functions with real coefficients. This polynomial turns up in what you wrote, because solving y^2 - 2xy - 1 = 0 for x gives x = (y^2 - 1)/(2y) just like you had; minimality says that this expression for x in terms of y is, in a certain limited sense, as simple as possible. A degree 2 polynomial is small enough that this is quite special, and the theory of Riemann surfaces made me hope that I could do something fun and clever. That was where t came in; but alas, it was not meant to be!)

Oh, and a last follow-up: So far, my results look great! Thanks for developing Stan; without it, what I wanted to do would have been insurmountable, and instead it has been easy.