# Correlation between Gaussian processes

**URL:** <https://discourse.mc-stan.org/t/correlation-between-gaussian-processes/21251>\
**Category:** Modeling\
**Tags:** gaussian-process, covariance, techniques\
**Created:** [March 12, 2021, 10:09am UTC](https://discourse.mc-stan.org/t/correlation-between-gaussian-processes/21251 "2021-03-12T10:09:58Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![bernibra](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bernibra/32/10918_2.png) [@bernibra](https://discourse.mc-stan.org/u/bernibra)\
**Post date:** [March 12, 2021, 10:09am UTC](https://discourse.mc-stan.org/t/correlation-between-gaussian-processes/21251/1 "2021-03-12T10:09:58Z")

</div>

Hey there,

I recently posted for the first time in this forum, and I had such quick and great responses that I can’t help it but ask another question regarding Gaussian processes that has been bugging me for a while (and for which I can’t seem to find anything relevant online). I was hoping someone could shed some light on it or point me towards the right direction.

Imagine that I have two variables \alpha\_i and \beta\_i for which I have basic information on the differences between indices i. In particular, I want to assume that for parameters \alpha\_i those differences are characterised by a distance matrix D^{\alpha}, and for parameters \beta\_i the same are characterised by a distance matrix D^{\beta}. As a starting point, I can write those such that their priors are independent multivariate distributions such that

\alpha\_i \sim \text{MVNormal}\left(\bar{\alpha}, K^{\alpha}\right),  
\beta\_i \sim \text{MVNormal}\left(\bar{\beta}, K^{\beta}\right),  
\bar{\alpha}, \bar{\beta}\sim \text{Normal}\left(0,1\right),  
K^{\alpha}\_{ij} = \eta\_{\alpha}\,\exp(- \rho\_{\alpha} {D^{\alpha}\_{ij}}^2) +\delta\_{ij}\,\sigma\_{\alpha},  
K^{\beta}\_{ij} = \eta\_{\beta}\,\exp(- \rho\_{\beta} {D^{\beta}\_{ij}}^2) +\delta\_{ij}\,\sigma\_{\beta},  
\sigma\_{\alpha}, \sigma\_{\beta}, \eta\_{\alpha}, \eta\_{\beta} \sim \text{Exponential}\left(1\right)  
\rho\_{\alpha}, \rho\_{\beta}\sim \text{Exponential}\left(0.5\right)

Assume, however, that I am expecting some correlation between \alpha\_i and \beta\_i. I can certainly do a post hoc analysis of this correlation by comparing their posterior distributions; but, I wonder if there is a way to include this correlation in the model, directly. Discussing things with people that know much more than me, I have been suggested multiple ways to do this. One way is by modelling \alpha\_i and \beta\_i together such that:

\begin{bmatrix} \alpha\_i \\ \beta\_i \end{bmatrix} \sim \text{MVNormal}\left( \begin{bmatrix} \bar{\alpha} \\ \bar{\beta} \end{bmatrix}, K\right),  
K=\begin{bmatrix} K^{\alpha} & \nu\\ \nu & K^{\beta} \end{bmatrix},

where \nu is the correlation between \alpha\_i and \beta\_i. Unfortunately, I have found little success doing this with simulated data. This could well be a result of me messing up things along the process, but I have a feeling that my covariance structure K is missing something. More specifically, I believe that there should be a factor multiplying \nu that accounts for the covariance between alphas, and the covariance between betas (i.e. K^{\alpha} and K^{\beta}).

Alternatively, I have been suggested to account for this correlation in the prior distribution of the corresponding hyperparameters (i.e. \rho, \sigma and \eta). That is, writing their prior distributions as a multivariate distribution with several \nu measuring the pairwise correlation between hyperparameters. However, I don’t necessarily see how I can then understand the overall correlation between \alpha\_i and \beta\_i via the correlation of their hyperparameters. Other interesting suggestions that I was given included modelling \alpha\_i and \beta\_i as correlated via the LKJ and then impose spatial autocorrelation on the residuals following a Gaussian process that include the distance matrices D, but I haven’t had a chance to try this yet.

Not sure if any of this makes sense or if I am missing something trivial or already discussed in the forum/documentation, but it would be wonderful to hear some ideas on how to approach this.

Thanks in advance for your time!

---

<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:** [March 12, 2021, 10:18am UTC](https://discourse.mc-stan.org/t/correlation-between-gaussian-processes/21251/2 "2021-03-12T10:18:21Z")

</div>

I just woke up and I’m making coffee…can you write it out as a matrix variate normal? There’s an easy conversion back to multivariate normal and then investigate that correlation matrix. See [Matrix normal distribution - Wikipedia](https://en.m.wikipedia.org/wiki/Matrix_normal_distribution)

---

<div class="post-metadata">

**Author:** ![bernibra](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bernibra/32/10918_2.png) [@bernibra](https://discourse.mc-stan.org/u/bernibra)\
**Post date:** [March 12, 2021, 10:39am UTC](https://discourse.mc-stan.org/t/correlation-between-gaussian-processes/21251/3 "2021-03-12T10:39:14Z")

</div>

I must admit is the first time I hear of matrix variate normal distributions. I’ll have a look and see where it takes me. Thanks a lot!

---

<div class="post-metadata">

**Author:** ![mike-lawrence](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mike-lawrence/32/59_2.png) [@mike-lawrence](https://discourse.mc-stan.org/u/mike-lawrence)\
**Post date:** [March 15, 2021, 11:25am UTC](https://discourse.mc-stan.org/t/correlation-between-gaussian-processes/21251/4 "2021-03-15T11:25:03Z")

</div>

How about modelling things with an intercept & difference parameterization such that:

\alpha\_i = I\_i + \frac{D\_i}{2}  
\beta\_i = I\_i - \frac{D\_i}{2}  
I\_i \sim \text{MVNormal}\left(\bar{I}, K^{I}\right),  
D\_i \sim \text{MVNormal}\left(\bar{D}, K^{D}\right),  
etc…

---

<div class="post-metadata">

**Author:** ![bernibra](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bernibra/32/10918_2.png) [@bernibra](https://discourse.mc-stan.org/u/bernibra)\
**Post date:** [March 15, 2021, 12:47pm UTC](https://discourse.mc-stan.org/t/correlation-between-gaussian-processes/21251/5 "2021-03-15T12:47:44Z")

</div>

Thank you for the suggestion @mike-lawrence. How would the correlation be defined in your example? Would the K^I and K^D be defined the same as K^{\alpha} and K^{\beta}, respectively?

---

<div class="post-metadata">

**Author:** ![mike-lawrence](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mike-lawrence/32/59_2.png) [@mike-lawrence](https://discourse.mc-stan.org/u/mike-lawrence)\
**Post date:** [March 15, 2021, 1:00pm UTC](https://discourse.mc-stan.org/t/correlation-between-gaussian-processes/21251/6 "2021-03-15T13:00:25Z")

</div>

Yup! So, taking the time to copy your post’s formulas and doing a find-and-replace, the full thing would be:

\alpha\_i = I\_i + \frac{D\_i}{2}  
\beta\_i = I\_i - \frac{D\_i}{2}  
I\_i \sim \text{MVNormal}\left(\bar{I}, K^{I}\right),  
D\_i \sim \text{MVNormal}\left(\bar{D}, K^{D}\right),  
I\_i \sim \text{MVNormal}\left(\bar{I}, K^{I}\right)  
D\_i \sim \text{MVNormal}\left(\bar{D}, K^{D}\right)  
\bar{I}, \bar{D}\sim \text{Normal}\left(0,1\right)  
K^{I}\_{ij} = \eta\_{I}\,\exp(- \rho\_{I} {D^{I}\_{ij}}^2) +\delta\_{ij}\,\sigma\_{I}  
K^{D}\_{ij} = \eta\_{D}\,\exp(- \rho\_{D} {D^{D}\_{ij}}^2) +\delta\_{ij}\,\sigma\_{D}  
\sigma\_{I}, \sigma\_{D}, \eta\_{I}, \eta\_{D} \sim \text{Exponential}\left(1\right)  
\rho\_{I}, \rho\_{D}\sim \text{Exponential}\left(0.5\right)

---

<div class="post-metadata">

**Author:** ![mike-lawrence](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mike-lawrence/32/59_2.png) [@mike-lawrence](https://discourse.mc-stan.org/u/mike-lawrence)\
**Post date:** [March 15, 2021, 1:08pm UTC](https://discourse.mc-stan.org/t/correlation-between-gaussian-processes/21251/7 "2021-03-15T13:08:52Z")

</div>

Though be sure to look at the current recommendations for priors for the GP parameters; I think exponential at least for the \rho parameters might be expected to behave badly because it’s concentrating mass in regions of the parameter space where the data can’t distinguish between values very strongly. I always get confused about whether folks’ definition of \rho translates colloquially to “lengthscale” or “inverse-lengthscale”/“wiggliness” (latter being my term :) ), but for both you want to stay away from values very close to zero and very large values. I like a `lengthscale ~ weibull(2,x)` for this reason, where it’s peaked at about .8 and x is the max distance in the data (or 1 if you simply scale the data so the max distance is 1).

---

<div class="post-metadata">

**Author:** ![Matias\_Guzman\_Naranjo](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/matias_guzman_naranjo/32/4977_2.png) [@Matias\_Guzman\_Naranjo](https://discourse.mc-stan.org/u/Matias_Guzman_Naranjo)\
**Post date:** [May 23, 2025, 12:03pm UTC](https://discourse.mc-stan.org/t/correlation-between-gaussian-processes/21251/8 "2025-05-23T12:03:12Z")

</div>

I am currently working on exactly this same problem. Did you have any success? Thanks!

---

<div class="post-metadata">

**Author:** ![Matias\_Guzman\_Naranjo](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/matias_guzman_naranjo/32/4977_2.png) [@Matias\_Guzman\_Naranjo](https://discourse.mc-stan.org/u/Matias_Guzman_Naranjo)\
**Post date:** [May 23, 2025, 12:18pm UTC](https://discourse.mc-stan.org/t/correlation-between-gaussian-processes/21251/9 "2025-05-23T12:18:06Z")

</div>

This is, I think, what brms does for correlated phylogenetic effects. However, if I understand it correctly (and I probably don’t), this requires an identical covariance matrix for both GPs, no?