# Using LKJCorr together with MvNormal

**URL:** <https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606>\
**Category:** version agnostic\
**Created:** [January 9, 2024, 3:09pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606 "2024-01-09T15:09:40Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![Velochy](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/velochy/32/7602_2.png) [@Velochy](https://discourse.pymc.io/u/Velochy)\
**Post date:** [January 11, 2024, 11:55am UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/21 "2024-01-11T11:55:12Z")

</div>

We are talking of a linear transformation here between a LKJ corr matrix and LKJ Cov matrix (multiplying both rows and columns with the sd vector). Therefore, if we took ` _LKJCholeksyCovRV_logp` and just removed the logp\_sd from the final sum, shouldnt’t it give us the logp for the cholesky factorization of the pure correlation matrix?

In which case I would structure things so that LKJCorrRV returns the cholesky factorization and LKJCorr allows you to parametrize if you want matrix or its cholesky factorization.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 11, 2024, 12:03pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/22 "2024-01-11T12:03:19Z")

</div>

> [@Velochy](#):
>
> We are talking of a linear transformation here between a LKJ corr matrix and LKJ Cov matrix

Is LKJ Cov matrix the LKJCholeskyCov one? If that is the case, you should be able to invert the linear transformation deterministically after you get draws from the LKJCholeskyCov. My guess is that you cannot.

If you were referring to another object ignore the message (and let me know which one it is)

---

<div class="post-metadata">

**Author:** ![Velochy](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/velochy/32/7602_2.png) [@Velochy](https://discourse.pymc.io/u/Velochy)\
**Post date:** [January 11, 2024, 12:07pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/23 "2024-01-11T12:07:17Z")

</div>

This is in effect the workaround I’m currently using, actually.

But I was thinking about [ENH: LKJCorr should return matrix rather than vector · Issue #7074 · pymc-devs/pymc · GitHub](https://github.com/pymc-devs/pymc/issues/7074) and how to get LKJCorr to actually output either the cholesky factorization or corr as matrix, without needing to sample sd\_dist in the first place.

Anyways, I unfortunately think I’m underqualified to actually do this as I can just very, very barely keep up with what you are saying unfortunately. So I think I’m tapping out.

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [January 11, 2024, 12:09pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/24 "2024-01-11T12:09:02Z")

</div>

Alternative suggestion then: a thin wrapper around `LKJCorr` that behaves like the wrapper around `LKJCholeskyCov`, that basically applies @Velochy 's original helper function to the outputs of `LKJCorr` if `chol=True`. I think handling the packing/unpacking and doing an “unnecessary” cholesky factorization if you want to do non-centered parameterization is already quite nice.

This would avoid having to worry about any logp stuff.

---

<div class="post-metadata">

**Author:** ![Velochy](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/velochy/32/7602_2.png) [@Velochy](https://discourse.pymc.io/u/Velochy)\
**Post date:** [January 11, 2024, 12:25pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/25 "2024-01-11T12:25:52Z")

</div>

Problem is, my original helper does not provide the cholesky factorization, just the full matrix.

My personal workaround is to sample from LKJCholeskyCov and then scale it back by it’s own sd vector i.e.

```auto
chol, corr, sd = LKJCholeskyCov(...., compute_corr=True)
corr_chol = (chol.T/sd).T

```

But this still requires specifying a distribution for sd\_dist and sampling from it spuriously

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [January 11, 2024, 12:28pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/26 "2024-01-11T12:28:02Z")

</div>

Something like this doesn’t work?

```auto
corr_flat = pm.LKJCorr(...)
corr = unpacking_helper(corr_lat)
chol_corr = pt.linalg.cholesky(corr)

```

---

<div class="post-metadata">

**Author:** ![Velochy](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/velochy/32/7602_2.png) [@Velochy](https://discourse.pymc.io/u/Velochy)\
**Post date:** [January 11, 2024, 12:33pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/27 "2024-01-11T12:33:09Z")

</div>

No it does 😃 I was actually not aware of pt.linalg as google did not provide any matches with “pytensor cholesky”. Thank you @jessegrabowski

In this case, you are right, the thin wrapper seems like a good plan, and one I can probably manage to implement and do a PR for. I’ll at least give it a try 🙂

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [January 11, 2024, 12:36pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/28 "2024-01-11T12:36:06Z")

</div>

Yeah google doesn’t index our docs very well, sadly. But things in pytensor are organized to follow numpy/scipy as closely as possible, so you can always try to look in the corresponding place. The [API docs](https://pytensor.readthedocs.io/en/latest/library/index.html) are here but they need some love ☹

Tag me on the PR and I’ll be happy to give it a review

---

<div class="post-metadata">

**Author:** ![Velochy](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/velochy/32/7602_2.png) [@Velochy](https://discourse.pymc.io/u/Velochy)\
**Post date:** [January 13, 2024, 2:15pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/29 "2024-01-13T14:15:22Z")

</div>

Turns out there might still be some issues though.

I tired running

```auto
    c_arr = pm.LKJCorr(name,eta=eta,n=n)
    tri = pt.zeros( (n,n) )
    tri = pt.subtensor.set_subtensor(tri[np.triu_indices(n,1)],c_arr)
    cmat = tri + tri.T + pt.diag(np.ones(n))
    corr_chol = pt.linalg.cholesky(cmat)

```

in the context of a slightly more complex model and got the following error:

> 

```auto
LinAlgError: 7-th leading minor of the array is not positive definite
Apply node that caused the error: Cholesky{lower=True, destructive=False, on_error='raise'}(Add.0)
Toposort index: 226
Inputs types: [TensorType(float64, shape=(None, None))]
Inputs shapes: [(12, 12)]
Inputs strides: [(96, 8)]
Inputs values: ['not shown']
Outputs clients: [[Dot22(Cholesky{lower=True, destructive=False, on_error='raise'}.0, Composite{...}.0), Dot22Scalar(Cholesky{lower=True, destructive=False, on_error='raise'}.0, Composite{...}.0, 3.0), Dot22Scalar(Cholesky{lower=True, destructive=False, on_error='raise'}.0, Composite{...}.0, 3.0)]]

```

You guys have any thoughts on how this might have come about? Because looking at the matrix with .eval, it looks like a perfectly fine correllation matrix.

Also, in the same situation

```auto
chol, corr, sd = LKJCholeskyCov(...., compute_corr=True)
corr_chol = (chol.T/sd).T

```

works just fine. Reverted to that for the time being.

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [January 13, 2024, 4:23pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/30 "2024-01-13T16:23:05Z")

</div>

Can you post a MVE that reproduces the bug? It samples fine if I am just trying to estimate the correlations, like this:

```python
import pymc as pm
import pytensor.tensor as pt
import numpy as np
import arviz as az

# Generate Data
n = 5
eta = 1

corr_values = pm.LKJCorr.dist(n=n, eta=eta)
corr = pt.zeros((n, n))
corr = pt.set_subtensor(corr[np.triu_indices(n, k=1)], corr_values)
corr = corr + corr.T + pt.identity_like(corr)

chol_corr = pt.linalg.cholesky(corr)

chol = pm.draw(chol_corr, 1)
y = pm.draw(pm.MvNormal.dist(mu=0, chol=chol), 100)

# PyMC Model
with pm.Model() as m:
    
    corr_values = pm.LKJCorr('corr_values', n=n, eta=eta)
    corr = pt.zeros((n, n))
    corr = pt.set_subtensor(corr[np.triu_indices(n, k=1)], corr_values)
    corr = corr + corr.T + pt.identity_like(corr)

    chol_corr = pt.linalg.cholesky(corr)
    y_hat = pm.MvNormal('y_hat', mu=0, chol=chol_corr, observed=y)
    idata = pm.sample()

```

Without seeing a model, my hypothesis is that the sampler is proposing illegal values that should just be discarded because they have `-np.inf` logp, but are causing a computation problem for you before that can happen (see [here](https://discourse.pymc.io/t/variable-lag-and-slicing/12003/7) for a similar situation in a totally unrelated context). If I’m right, the fix is just to pass `on_error='nan'` to `pt.linalg.cholesky`. That will prevent it from erroring out if a proposal is not PSD, and the proposal will be rejected anyway so it’s no harm no foul.

---

<div class="post-metadata">

**Author:** ![Velochy](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/velochy/32/7602_2.png) [@Velochy](https://discourse.pymc.io/u/Velochy)\
**Post date:** [January 13, 2024, 11:02pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/31 "2024-01-13T23:02:08Z")

</div>

It seems the minimal reproducible example is the same code you sent me, but with a higher n (say 10 or 20). I actually got the error first time I ran it as it was (with n=5), but with smaller n it works very often but the error rate very clearly increases as n does.

I also tried setting the on\_error=‘nan’ and while the sampling runs then, it is clear from looking at the trace it actually fails miserably, so it is not really a fix but just sweeps the error under the rug.

Do you have any more ideas as to what it might be. I’m pretty stumped.

---

<div class="post-metadata">

**Author:** ![Velochy](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/velochy/32/7602_2.png) [@Velochy](https://discourse.pymc.io/u/Velochy)\
**Post date:** [January 14, 2024, 12:04am UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/32 "2024-01-14T00:04:14Z")

</div>

Having thought about this for a bit:

Isn’t the problem that in backwards sampling, the logp for LKJCorr in no way selects for positive semidefiniteness. The interval transform only transforms all the values to within [-1,1] range and then LKJCorr logp applies the determinant-based distribution which already assumes a semi-definite matrix, but neither of these really nudges the numbers to form a semi-definite matrix.

This is consistent with failures increasing as n gets larger, as PSD matrices become increasingly less likely. As for why the probability does not drop off a cliff after 2 or 3, my guess is that the interval transform steers away from values close to -1 and 1 in favor of more middling values which lets the ones on the diagonal dominate for a while.

This is also consistent with the sampling failing if instead of chol=chol\_corr you use cov=corr in the MvNormal parametrization. Same reason - it ends up proposing a non-PSD matrix. In this case, the error message contains the initial value though, so it is easy to check:

```auto
vec = [-0.39718353, 0.94680841, -0.41832789, 0.67318716, 0.74000446,
       -0.91103117, 0.17496013, 0.05499186, -0.60507651, -0.37989516,
        0.00281214, -0.83842367, 0.04520434, 0.51635372, -0.26586628,
        0.96277075, 0.66443148, -0.46217671, 0.77692162, -0.20452818,
        0.99873618, -0.23998811, 0.65304258, 0.98352337, 0.64847332,
       -0.08343583, 0.52405396, 0.15585091, 0.31938468, 0.90564176,
       -0.35090953, -0.43284892, -0.21282992, -0.48646797, -0.50581992,
       -0.06065079, -0.57660361, 0.90999947, 0.07498183, -0.01768012,
       -0.54075488, 0.8193952 , -0.13612187, 0.29677601, 0.07794663]
n=10
corr = np.zeros((n, n))
corr[np.triu_indices(n, k=1)] = vec
corr = corr + corr.T + np.eye(n)
print(np.linalg.eig(corr)[1].min())

```

indeed returns -0.7785569784890056 i.e. a negative eigenvalue =\> not PSD

So it looks like LKJCorr as it currently stands is fundamentally broken, i.e. it does not sample from the LKJ distribution as it does not respect the PSD requirement.

I think the way to fix this would be to move to sampling cholesky decomposition (as in LKJCholeskyCorr) as the existance of it already implies the matrix is PSD hence selecting for it.

Although - at that point - maybe it makes sense to just rewrite LKJCorr as a thin wrapper around LKJCholeskyCov that just returns the ‘corr’ part?

Thoughts @ricardoV94 @jessegrabowski ?

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 14, 2024, 7:32am UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/33 "2024-01-14T07:32:33Z")

</div>

> [@Velochy](#):
>
> Although - at that point - maybe it makes sense to just rewrite LKJCorr as a thin wrapper around LKJCholeskyCov that just returns the ‘corr’ part?

I’m not sure about this part. In the end you will need the logp, and for sampling efficiency, a transform that unconstrains LKJ Corr matrices. It’s unclear to me how much logic would be shared between these methods and the ones LKJCholeskyCov uses.

Also LKJCholeskyCov has more parameters than the Corr counterpart

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 14, 2024, 7:39am UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/34 "2024-01-14T07:39:37Z")

</div>

We may need something like this: [tfp.bijectors.CorrelationCholesky &nbsp;|&nbsp; TensorFlow Probability](https://www.tensorflow.org/probability/api_docs/python/tfp/bijectors/CorrelationCholesky)

I was quite surprised that we were using just an IntervalTransform, so much I missed the regression described in [this issue](https://github.com/pymc-devs/pymc/issues/7002) when we forbid univariate transforms on multivariate rvs.

Apparently it was indeed not sufficient?

Edit: the Stan pages are pretty informative as usual:

> **[10.12 Cholesky factors of correlation matrices | Stan Reference Manual](https://mc-stan.org/docs/reference-manual/cholesky-factors-of-correlation-matrices-1.html)**
>
> Stan reference manual specifying the syntax and semantics of the Stan programming language.

> **[10.11 Cholesky factors of covariance matrices | Stan Reference Manual](https://mc-stan.org/docs/reference-manual/cholesky-factors-of-covariance-matrices-1.html)**
>
> Stan reference manual specifying the syntax and semantics of the Stan programming language.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 14, 2024, 8:23am UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/35 "2024-01-14T08:23:57Z")

</div>

> [@Velochy](#):
>
> Same reason - it ends up proposing a non-PSD matrix. In this case, the error message contains the initial value though, so it is easy to check:

We do have a check in the logp that should reject non PSD values. The problem is the logp of MvNormal will raise before we get the chance of rejecting that proposal, unless you don’t set that `on_error` flag.

> <https://github.com/pymc-devs/pymc/blob/6f8f9eef9e91999680877e3ed48dd903011de79c/pymc/distributions/multivariate.py#L1635>

Actually our check is to just try to do a cholesky and see if it raises :D. Maybe we should use the eigenvalues as a check? I guess that’s slower, or maybe not sufficient?

> <https://github.com/pymc-devs/pymc/blob/6f8f9eef9e91999680877e3ed48dd903011de79c/pymc/distributions/multivariate.py#L834>

Even if we reject non PSD values, we should try to have a transform that guarantees invalid values are never proposed, or sampling will suck like you noticed.

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [January 14, 2024, 10:16am UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/37 "2024-01-14T10:16:58Z")

</div>

> [@Velochy](#):
>
> `print(np.linalg.eig(corr)[1].min())`

The 0 return is the eigenvalues, you are checking the eigenvectors here. Better to just use `np.linalg.eigvals` to avoid this problem. But your point is still correct, and there are negative eigenvalues.

I believe the problem is in the transformation, because when you just forward sample from the model (e.g. with `pm.sample_prior_predictive`), all the correlation matrices will be PSD.

> [@ricardoV94](#):
>
> Actually our check is to just try to do a cholesky and see if it raises :D. Maybe we should use the eigenvalues as a check? I guess that’s slower, or maybe not sufficient?

For now I think it’s fine how we do it. In the future if we have special rewrites for eigenvalues on different types of matrices it might be faster to check eigenvalues. Eigenvalues can also be imaginary, and right now the imaginary component of returns get set to 0 during casting to `float64`, so we’d miss negative eigenvalues of the form `0-xj`

---

<div class="post-metadata">

**Author:** ![Velochy](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/velochy/32/7602_2.png) [@Velochy](https://discourse.pymc.io/u/Velochy)\
**Post date:** [January 14, 2024, 10:39am UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/38 "2024-01-14T10:39:47Z")

</div>

I’m taking an engineering / product management approach here.

LKJCorr seems to have been broken for a while without anyone noticing. Presumably because people using it have only used small matrices (n=2 or 3 is all I could find googling on forums).

Now. We could spend time implementing the logp. Probably at least a day of my time, probably less for either of you, but I’m still guessing hours.

Or - we could just rewrite LKJCorr as a thin wrapper around LKJCholeskyCov by providing it with sd\_dist=pm.LogNormal.dist() (or any other well-behaving distribution, really) and just marginalize it out by returning corr\_chol = (chol.T/sd).T (or returned corr directly).  
We would get correct behavior. We would not significantly effect existing users (as n=2 and n=3 cases are so small adding extra variables to sample over has a very marginal overhead ). And most notably - this would take a fraction of the time to accomplish.

Anyways - you guys are core team so you are better positioned to make such decisions. From my perspective though, it seems like an easy win as opposed to a hard but ideal solution to a problem no-one except me is having.

---

<div class="post-metadata">

**Author:** ![Velochy](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/velochy/32/7602_2.png) [@Velochy](https://discourse.pymc.io/u/Velochy)\
**Post date:** [January 14, 2024, 11:16am UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/39 "2024-01-14T11:16:58Z")

</div>

To maybe better explain what I mean, I took the 15 minutes required to write that thin wrapper:

> <https://github.com/velochy/pymc/blob/70c17eb87ddcffdc269b5446cf5451c7beb3abf4/pymc/distributions/multivariate.py#L1545>

I’m not saying it is the ideal solution. But it might be the pragmatic one.

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 14, 2024, 12:24pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/40 "2024-01-14T12:24:59Z")

</div>

You are discarding the variances right? That’s still inefficient because you are sampling more parameters than you actually need. I don’t know about this specific case, but that can sometimes introduce bad indeterminacies in the posterior (unless the extra parameters are completely orthogonal)

It may be less inefficient than the not completely constrained LKJCorr we have right now, but not great either.

To be pedantic, the LKJCorr is not broken, everything it does seems correct AFAICT, we just don’t have a transformation that allows NUTS to sample in a completely unconstrained space, which means that for some models you will be getting divergences.

That’s something we definitely should address.

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [January 14, 2024, 12:58pm UTC](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606/41 "2024-01-14T12:58:12Z")

</div>

I believe the TFP implementation has everything we need, I’m looking at it now. @Velochy could you open an issue on github explaining that LKJCorr breaks during sampling?

[Previous page](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606.md?page=1)

[Next page](https://discourse.pymc.io/t/using-lkjcorr-together-with-mvnormal/13606.md?page=3)