# Constraining the trace of a covariance matrix
**URL:** https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815
**Category:** Modeling
**Created:** [July 10, 2018, 8:26pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815 "2018-07-10T20:26:45Z")
**Posts on this page:** 13
**Page:** 1
### Author:  [@ihincks](https://discourse.mc-stan.org/u/ihincks)
#### Post date: [July 10, 2018, 8:26pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/1 "2018-07-10T20:26:45Z")
I am interested in an inference problem where the posterior distribution is over positive semi-definite matrices whose trace is equal to 1. I have not begun working on this, but before I do, I would be interested to hear any opinions from those experienced users about how best to implement this, as in, what to put in the `parameters` section of the stan program. I understand there are constrained native types for cholesky factors, correlations, and covariance matrices.
I once wrote pair of inverse functions that transform between \mathbb{R}^{d^2-1} and the set of d\times d unit-trace positive semi-definite operators with the constraint based on stick-breaking ([see here](https://gist.github.com/ihincks/0b636988959f1226c2107e7a78516aff)).
Would implementing this in stan be a useful approach, or is there a simpler way? Here is a blurb copied from the docstring of `vec_to_state()`:
> The unit trace condition is met by using a stick-breaking like trick, where
> a given value of `x` doesn’t directly store a cholesky factor matrix element,
> but rather, it stores (after being shrunk) the fraction of the remaining
> budget of the trace until `1`. Note that the trace of a positive operator
> is the square Hilbert-Schmidt norm of its cholesky factor, so this is a
> 2-norm version of stick breaking.
Thanks for any help.
---
### Author:  [@James\_Savage](https://discourse.mc-stan.org/u/James_Savage)
#### Post date: [July 10, 2018, 9:27pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/2 "2018-07-10T21:27:12Z")
Could you just define the square of the scale vector to be a simplex, give it a Dirichlet prior, and proceed as normal?
Something like
```
data {
int P;
}
transformed data {
vector[P] alpha = rep_vector(1.0/P, P);
}
parameters {
simplex[P] sqr_scale;
corr_matrix[P] Omega;
}
transformed parameters {
matrix[P, P] Sigma = quad_form_diag(Omega, sqrt(sqr_scale));
}
model {
target += dirichlet_lpdf(sqr_scale | alpha);
target += lkj_corr_lpdf(Omega | 3);
}
```
In R then
```
compiled_model <- stan_model("trace_sums_to_1.stan")
fit <- sampling(compiled_model, data = list(P = 5), chains = 1)
sum(diag(matrix(unlist(as.data.frame(fit, pars = "Sigma")[1000,]), 5, 5)))
> 1
```
This puts no weight on scales of 1 or 0, which might not be desirable.
Cheers,
Jim
---
### Author:  [@ihincks](https://discourse.mc-stan.org/u/ihincks)
#### Post date: [July 10, 2018, 11:35pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/3 "2018-07-10T23:35:33Z")
Thanks, this is exactly what I am looking for! (though I will need to muck around a bit to simulate complex matrices).
---
### Author:  [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [July 11, 2018, 5:35am UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/4 "2018-07-11T05:35:07Z")
> [@ihincks](#):
>
> so this is a 2-norm version of stick breaking.
That’s exactly what we do with correlation matrices under the hood, by the way. The full definition with Jacobian is in the manual in the chapter on constrained parameter transforms.
> [@James\_Savage](#):
>
> Could you just define the square of the scale vector to be a simplex, give it a Dirichlet prior, and proceed as normal?
Neat!
---
### Author:  [@ihincks](https://discourse.mc-stan.org/u/ihincks)
#### Post date: [July 11, 2018, 3:49pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/5 "2018-07-11T15:49:39Z")
> That’s exactly what we do with correlation matrices under the hood, by the way.
Cool. It seems the main difference is that cholesky stick-breaking operates per-row, while trace constraining stick-breaking operates on all entries together.
While reading through, I found a minor typo in the subsection “Correlation Matrix Inverse Transform”: It reads that tanh maps into (0,1), rather than (-1,1).
---
### Author:  [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)
#### Post date: [July 11, 2018, 7:19pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/6 "2018-07-11T19:19:24Z")
> [@ihincks](#):
>
> While reading through, I found a minor typo in the subsection
Nice, Bob keeps an issue open for changes like this: [next manual, 2.18(++) · Issue #15 · stan-dev/docs · GitHub](https://github.com/stan-dev/stan/issues/2552) , so go ahead and post it there!
---
### Author:  [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [July 12, 2018, 11:06pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/7 "2018-07-12T23:06:40Z")
> [@ihincks](#):
>
> a minor typo in the subsection “Correlation Matrix Inverse Transform”: It reads that tanh maps into (0,1), rather than (-1,1).
You’re right. I’ll fix for next release.
> [@bbbales2](#):
>
> an issue open for changes like this: [https://github.com/stan-dev/stan/issues/2552](https://github.com/stan-dev/stan/issues/2552)
Thanks for following up with an issue, @ihincks—that way, it won’t get lost.
---
### Author:  [@ihincks](https://discourse.mc-stan.org/u/ihincks)
#### Post date: [July 17, 2018, 3:12pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/8 "2018-07-17T15:12:42Z")
I have another fundamental question about the same problem. My prior distribution is based on a Ginibre ensemble, parameterized by `K` and `D` as follows:
```
data {
int D; // Dimension
int K; // Number of Ginibre subsystems
}
parameters {
matrix[D,K] X;
}
transformed parameters {
matrix[D,D] rho;
rho = X * X';
rho /= trace(rho);
}
model {
// Ginibre DxK prior
for (idx_row in 1:D) {
X[idx_row] ~ normal(0,1);
}
}
```
This is very convenient to implement, just iid normals on `X`. However, say `D=K`, then there are more parameters than necessary: `X` requires `D*D` real numbers, but `rho` is symmetric with trace 1 and thus requires `D*(D+1)/2-1` parameters. My guess is that this is a bad thing to do because, depending on the likelihood added to the model, it might become multimodal. Is this correct? Is it better to use a 1-1 and deal with the jacobian, or will I be okay with this more convenient-to-implement parameterization?
Thanks again.
(Note, I am still interested in the complex valued problem where X has real and imaginary components sampled iid from N(0,1), but I have stated it as above to be cleaner. James’ trick does not work as cleanly in the complex case and I think I would need to write the transform and its jacobian from scratch, and figure out the prior on rho, etc.)
---
### Author:  [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [July 17, 2018, 3:22pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/9 "2018-07-17T15:22:10Z")
> [@ihincks](#):
>
> Is it better to use a 1-1 and deal with the jacobian,
I think this should be OK. It’s what we do for generating points on a sphere. There are only N-1 degrees of freedom for a sphere embedded in N dimensions, but we use N parameters, give each a normal distribution, then project to the sphere. That provides a uniform distribution on the points on the sphere at the cost of some unit normals. This algorithm comes from [Marsaglia 1972](https://projecteuclid.org/download/pdf_1/euclid.aoms/1177692644).
By the way, `to_vector(X) ~ normal(0, 1)` should be more efficient than the loop since we can share a lot of computations across the rows. The copies have to happen anyway to pull the rows out as Eigen stores column major, so for big matrices, it’s even more of a win.
---
### Author:  [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [July 17, 2018, 3:27pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/10 "2018-07-17T15:27:58Z")
> [@ihincks](#):
>
> My guess is that this is a bad thing to do because, depending on the likelihood added to the model, it might become multimodal.
Sorry, missed that the first go `round.
Yes, this can be a problem with things like matrix multiplication because (-x)^2 = x^2.
In the sampling on a sphere case, there’s only a single line from the origin that results in the same point on a sphere, so nothing gets multimodal.
Without a way to control the multimodality of the parameters `X`, it can be problematic in the posterior monitoring `X`. On the other hand, it might not matter in that `rho` won’t wind up being multimodal.
Also, I’m not sure what distribution you get from `X * X'` given a normal distribution on the `X` components. What was your intended distribution on `rho`? The best thing to do would be to check whether that’s what you get with this coding.
---
### Author:  [@BFiles](https://discourse.mc-stan.org/u/BFiles)
#### Post date: [July 17, 2018, 5:00pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/11 "2018-07-17T17:00:17Z")
> [@Bob\_Carpenter](#):
>
> Yes, this can be a problem with things like matrix multiplication because (-x)^2 = x^2.
>
> In the sampling on a sphere case, there’s only a single line from the origin that results in the same point on a sphere, so nothing gets multimodal.
>
> Without a way to control the multimodality of the parameters `X`, it can be problematic in the posterior monitoring `X`. On the other hand, it might not matter in that `rho` won’t wind up being multimodal.
I do not want to hijack this thread, but I have been looking for a clarification of a point here. Continuing the above example, if we only care about rho, and `X` is for some reason of no interest, do we need to worry about the posterior monitoring on `X`? In other words, if we get `Rhat=1` for `rho`, can we safely ignore `Rhat>1` for `X`?
Note: if more details are needed about my specific case, I can start a new thread to avoid derailing this one further.
---
### Author:  [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [July 17, 2018, 5:17pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/12 "2018-07-17T17:17:09Z")
> [@BFiles](#):
>
> if we only care about rho, and `X` is for some reason of no interest, do we need to worry about the posterior monitoring on `X`
The answer is complicated, as it depends on the properties of particular transforms.
In some cases, no, we don’t need to worry about it. For example, you can get away with this to generate a uniform prior in (0, 1) on `alpha`:
```no-highlight
parameters {
real alpha_raw;
}
transformed parameters {
real alpha = fabs(alpha_raw);
}
model {
alpha ~ normal(0.2, 0.05)
}
```
If the posterior for `alpha` concentrates away from 0, as in this example, the posterior on `alpha_raw` will be multimodal. This should be OK for `alpha` in this simple case. Here’s what the output looks like for 4 chains and default settings:
```no-highlight
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
theta 0.00 0.14 0.21 -0.28 -0.20 0.01 0.20 0.29 2 4.42
theta_pos 0.20 0.00 0.05 0.10 0.17 0.20 0.23 0.30 1779 1.00
lp__ -1.23 0.01 0.69 -3.18 -1.41 -0.96 -0.78 -0.73 2188 1.00
```
By inspection, the 95% intervals look good—they should be about +/- 1.96 sd and they’re (0.1, 0.3).
You can see from Rhat and n\_eff that `theta` didn’t mix—two or more chains got stuck in different modes.
The way our adaptation works within a single chain actually mitigates the most problematic aspect of transforming a multimodal prior, namely that it’s likely that if `alpha` is heavily constrained, then `alpha_raw` won’t try to switch modes in the posterior during warmup, so adaptation will be OK. If a chain would explore both modes, then we run into a problem with adaptation trying to account for them both, which won’t be nearly as efficient.
---
### Author:  [@ihincks](https://discourse.mc-stan.org/u/ihincks)
#### Post date: [July 17, 2018, 8:17pm UTC](https://discourse.mc-stan.org/t/constraining-the-trace-of-a-covariance-matrix/4815/13 "2018-07-17T20:17:52Z")
Thanks for the very helpful info and the example, I will give it some more thought. The posterior on rho should be unimodal, and all of my tests are looking very much like your example in terms of Rhat and n\_eff.
> [@BFiles](#):
>
> if we only care about rho, and `X` is for some reason of no interest, do we need to worry about the posterior monitoring on `X`
Glad this is of broader interest to some, this is exactly what I am wondering too.
> [@Bob\_Carpenter](#):
>
> `to_vector(X) ~ normal(0, 1)` should be more efficient
I suspected there must be a better way…thanks!