# Lkj\_corr\_lpdf: Correlation matrix is not positive definite

**URL:** <https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995>\
**Category:** General\
**Created:** [July 6, 2023, 7:59pm UTC](https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995 "2023-07-06T19:59:18Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![ramiro](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ramiro/32/7587_2.png) [@ramiro](https://discourse.mc-stan.org/u/ramiro)\
**Post date:** [July 6, 2023, 7:59pm UTC](https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995/1 "2023-07-06T19:59:18Z")

</div>

I am working with a multi normal distribution, my model of interest is a true multi-normal distribution but I am starting with something simple and build from it. The model is the following:

y \sim \mathcal{N}(0,\Sigma)

where the covariance matrix \Sigma is just a diagonal matrix

\Sigma \sim \begin{pmatrix} \sigma^2 & 0 \\ \ldots \\ 0 & \sigma^2 \end{pmatrix}

and

\sigma \sim \mathcal{G}(11,0.3)

I am trying to infer \sigma.

This can work without the multi-normal of course, but I am trying to get it to work with this to learn from it and slowly add more complexity. If I introduce the `lkj_corr` prior I start getting errors. Below is the stan code:

```plaintext
data {
  int<lower=0> N;
  vector[N] y;
}

parameters {
  real<lower=0> sigmaD;
  corr_matrix[N] Omega;
}

transformed parameters {
  vector<lower=0>[N] sigma;
  for (i in 1:N) {
    sigma[i] = sigmaD^2;
  }
}

model {
  sigmaD ~ gamma(11,0.3);
  Omega ~ lkj_corr(3);
  y ~ multi_normal(rep_vector(0,N), quad_form_diag(Omega, sigma));
}

```

But I get errors such as these:

“Rejecting initial value:  
Error evaluating the log probability at the initial value.  
Exception: lkj\_corr\_lpdf: Correlation matrix is not positive definite.  
Initialization between (-2, 2) failed after 100 attempts.  
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.  
Initialization failed.”

I don’t quite understand why this is happening. Can someone explain? Any suggestions?

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [July 6, 2023, 10:11pm UTC](https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995/2 "2023-07-06T22:11:06Z")

</div>

@spinkney any insights?

---

<div class="post-metadata">

**Author:** ![billWalker](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/billwalker/32/97_2.png) [@billWalker](https://discourse.mc-stan.org/u/billWalker)\
**Post date:** [July 7, 2023, 7:33am UTC](https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995/3 "2023-07-07T07:33:05Z")

</div>

I think the guidance is to always use the cholesky factorized version, `lkj_corr_cholesky`, for better stability.

See [https://mc-stan.org/docs/stan-users-guide/multivariate-hierarchical-priors.html#optimization-through-cholesky-factorization](https://mc-stan.org/docs/stan-users-guide/multivariate-hierarchical-priors.html#optimization-through-cholesky-factorization) in the guide.

---

<div class="post-metadata">

**Author:** ![andre.pfeuffer](https://avatars.discourse-cdn.com/v4/letter/a/e480ec/32.png) [@andre.pfeuffer](https://discourse.mc-stan.org/u/andre.pfeuffer)\
**Post date:** [July 7, 2023, 7:46am UTC](https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995/4 "2023-07-07T07:46:33Z")

</div>

@billWalker all correct! Agreed. one should better work on the cholesky form.

The problem arise with increasing likelihood if the dimension of the correlation matrix  
increases. In the code above the dimension D = N, assigning each datapoint it’s own  
variance which is cannot be estimated.

Even if D \< N and D is large (\>30) the sampler may run into problems. I figured two  
ways around it:

One may set in the sampler statement `init = 0`, which is a quick hack, but not  
recommended. Another is to do your own initialization. See Stan manuals for details.  
Set the initialization of the correlation matrix or better the cholesky form to  
to Identity matrix same dimension, eg.

```plaintext
initf2 <- function(chain_id = 1) {
# cat("chain_id =", chain_id, "\n")
list(
# ...
   Omega = diag(D)
)

```

---

<div class="post-metadata">

**Author:** ![ramiro](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ramiro/32/7587_2.png) [@ramiro](https://discourse.mc-stan.org/u/ramiro)\
**Post date:** [July 7, 2023, 2:33pm UTC](https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995/5 "2023-07-07T14:33:30Z")

</div>

Thank you very much Bill, however, I don’t think the cholesky will save this, it seems from @andre.pfeuffer’s response that it has more to do with the size of the correlation matrix.

---

<div class="post-metadata">

**Author:** ![ramiro](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ramiro/32/7587_2.png) [@ramiro](https://discourse.mc-stan.org/u/ramiro)\
**Post date:** [July 7, 2023, 2:50pm UTC](https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995/6 "2023-07-07T14:50:37Z")

</div>

Thank you very much!! I see that if the matrix is too big this multivariate normal or cholesky will have problems, I might be able to find a way around this for my problem though!.

---

<div class="post-metadata">

**Author:** ![spinkney](https://avatars.discourse-cdn.com/v4/letter/s/dec6dc/32.png) [@spinkney](https://discourse.mc-stan.org/u/spinkney)\
**Post date:** [July 7, 2023, 2:56pm UTC](https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995/7 "2023-07-07T14:56:52Z")

</div>

You can try using the onion method to construct the matrix. This works for sizes much larger than 30.

```stan
functions {
  array[] vector create_shapes(int K, real eta) {
    real alpha = eta + (K - 2) / 2.0;
    array[2] vector[K - 1] shape;
    
    shape[1, 1] = alpha;
    shape[2, 1] = alpha;
    for (k in 2 : (K - 1)) {
      alpha -= 0.5;
      shape[1, k] = k / 2.0;
      shape[2, k] = alpha;
    }
    
    return shape;
  }
  
  matrix lkj_onion(int K, row_vector l, vector R2, data array[] vector shape) {
    matrix[K, K] L = rep_matrix(0, K, K); // cholesky_factor corr matrix
    {
      int start = 1;
      int end = 2;
      
      L[1, 1] = 1.0;
      L[2, 1] = 2.0 * R2[1] - 1.0;
      L[2, 2] = sqrt(1.0 - square(L[2, 1]));
      for (k in 2 : (K - 1)) {
        int kp1 = k + 1;
        row_vector[k] l_row = segment(l, start, k);
        real scale = sqrt(R2[k] / dot_self(l_row));
        L[kp1, 1 : k] = l_row[1 : k] * scale;
        L[kp1, kp1] = sqrt(1.0 - R2[k]);
        start = end + 1;
        end = start + k - 1;
      }
    }
    return L;
  }
}
data {
  int<lower=0> K; // dim of corr mat
  real<lower=0> eta; // lkj param
}
transformed data {
  array[2] vector[K - 1] shapes = create_shapes(K, eta);
}
parameters {
  row_vector[choose(K, 2) - 1] l; // do NOT init with 0 for all elements
  vector<lower=0, upper=1>[K - 1] R2; // first element is not really a R^2 but is on (0,1)  
}
transformed parameters {
  cholesky_factor_corr[K] L = lkj_onion(K, l, R2, shapes);
}
model {
  target += std_normal_lpdf(l);
  target += beta_lpdf(R2 | shapes[1], shapes[2]);
}
generated quantities {
  matrix[K, K] Sigma = multiply_lower_tri_self_transpose(L);
}

```

---

<div class="post-metadata">

**Author:** ![jmh530](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jmh530/32/4579_2.png) [@jmh530](https://discourse.mc-stan.org/u/jmh530)\
**Post date:** [July 11, 2023, 3:18pm UTC](https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995/8 "2023-07-11T15:18:18Z")

</div>

That looks nice. What’s the intuition of the R2 parameter? Is that controlling the average level of correlation or partial correlations? I also don’t have a good sense of what the create\_shapes function is doing.

For the OP, I’ve run into problems with the LKJ prior over the years. My problems were associated with problems where the variables are multivariate normal but with high correlations (like in a factor model). So the correlation parameters tend to be correlated in weird non-linear ways that Stan has difficulty with.

---

<div class="post-metadata">

**Author:** ![defjaf](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/defjaf/32/413_2.png) [@defjaf](https://discourse.mc-stan.org/u/defjaf)\
**Post date:** [July 11, 2023, 10:41pm UTC](https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995/9 "2023-07-11T22:41:12Z")

</div>

I am possibly missing something crucial , but I think you are vastly overcomplicating the problem. If you really know your correlation matrix to be diagonal, with the same variance along that diagonal, it’s just the product of identical univariate Gaussians from which you’re just trying to estimate the variance. Isn’t this just a textbook problem?

You don’t need to work with the Cholesky factorisation: you can do it by hand, and it’s just the matrix with sigma along the diagonal. But indeed you can write the whole problem as a simple sampler for sigma.

More likely I am missing a subtlety.

---

<div class="post-metadata">

**Author:** ![ramiro](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ramiro/32/7587_2.png) [@ramiro](https://discourse.mc-stan.org/u/ramiro)\
**Post date:** [July 12, 2023, 1:40pm UTC](https://discourse.mc-stan.org/t/lkj-corr-lpdf-correlation-matrix-is-not-positive-definite/31995/10 "2023-07-12T13:40:17Z")

</div>

Thank you for your answer!! As mentioned:

> [@ramiro](#):
>
> This can work without the multi-normal of course, but I am trying to get it to work with this to learn from it and slowly add more complexity

My “real problem” is much more complicated, and the only reason that I am doing it like this is just so that I can get it to work in a simple case and then build from it. What I am aiming for is not a diagonal matrix, but rather a hidden markov model where there a multi-normal distribution with covariation between consecutive observations, but we have to start somewhere!!

Thank you again very much for your wanting to help!