# Extending Gaussian Process functionality: Coregion and beyond

**URL:** <https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479>\
**Category:** Development\
**Created:** [October 21, 2017, 11:55pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479 "2017-10-21T23:55:28Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![jordan-melendez](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jordan-melendez/32/399_2.png) [@jordan-melendez](https://discourse.pymc.io/u/jordan-melendez)\
**Post date:** [October 21, 2017, 11:55pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/1 "2017-10-21T23:55:28Z")

</div>

Let’s discuss what steps would be required to extend the current GP api in some meaningful ways, given the [upcoming](https://github.com/pymc-devs/pymc3/pull/2635) `MatrixNormal` distribution. Some possibilities include

- `Coregion` kernels (as in [GPflow](https://github.com/GPflow/GPflow/blob/69364c2b79e1e6dd0fe7419377a044c0bdedb4f5/gpflow/kernels.py#L664) or [GPy](http://gpy.readthedocs.io/en/deploy/_modules/GPy/kern/src/coregionalize.html#Coregionalize), which uses a [mixed noise](http://gpy.readthedocs.io/en/deploy/GPy.likelihoods.html#module-GPy.likelihoods.mixed_noise) likelihood. Does that have any benefit over `MatrixNormal`?)
- `Marginal2D`
- Multiple curve observations of same GP ([works](https://discourse.pymc.io/t/multiple-uncertain-function-observations-of-the-same-gaussian-process/400) without `MatrixNormal`, but needs testing)
- I personally enjoy `plot_gp_dist` and would like to see if it could be modified to handle the above points
- Other ideas?

These features would help in my research, so this is partly a selfish plea for your help, but I really do think that PyMC3 would benefit from these additions. This post is just meant to spur discussion on these points.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [October 22, 2017, 6:44am UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/2 "2017-10-22T06:44:09Z")

</div>

These are all great suggestions. One question: what do you mean by `Marginal2D`? We can already do `Marginal2D` by unpacking the `y` (latent) into 1D.

---

<div class="post-metadata">

**Author:** ![jordan-melendez](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jordan-melendez/32/399_2.png) [@jordan-melendez](https://discourse.pymc.io/u/jordan-melendez)\
**Post date:** [October 22, 2017, 4:35pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/3 "2017-10-22T16:35:52Z")

</div>

That is a good question. That was an idea, along with `Latent2D`, raised by @bwengals in [this](https://discourse.pymc.io/t/multiple-uncertain-function-observations-of-the-same-gaussian-process/400/15?u=jordan-melendez) reply. I wasn’t entirely clear on it either, so maybe he can elaborate.

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [October 24, 2017, 5:17pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/4 "2017-10-24T17:17:06Z")

</div>

Sorry for my slowness, been a busy week… should clear up soon.

I completely agree that PyMC3 would benefit from these additions. There have already been a few people on discourse asking about stuff that could be helped by this. The whole point is to benefit peoples research and work!

For the reasons @junpenglao said, probably not the best name. How about `MarginalKron` and `LatentKron` to emphasize that the implementation relies on a Kronecker structured covariance matrix? I think the implementation of `LatentKron` should follow `Latent` nearly exactly, except with appropriate input args and `MatrixNormal` instead of `MvNormal`. `MarginalKron` will be the same when the covariance structure is $K\_1 \otimes K\_2$. If we want to also have the structure $K\_1 \otimes K\_2 + \sigma^2 I$, where the noise is diagonal across the whole covariance matrix, then we need to use the eigendecomposition instead of Cholesky ([ie Ch 5 here](http://mlg.eng.cam.ac.uk/pub/pdf/Saa11.pdf)).

I need to think a bit more about a `Coregion` covariance to say more, but from what I understand I don’t think there should be too much difficulty.

---

<div class="post-metadata">

**Author:** ![jordan-melendez](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jordan-melendez/32/399_2.png) [@jordan-melendez](https://discourse.pymc.io/u/jordan-melendez)\
**Post date:** [October 26, 2017, 3:30pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/5 "2017-10-26T15:30:33Z")

</div>

@bwengals that thesis is a good resource. Does theano provide anything like the Kronecker matrix-vector and Kronecker matrix-matrix products, defined as `kron_mvprod` and `kron_mmprod` respectively in the thesis?

I assume that it would be nice to allow for an arbitrary number of Kronecker products, right? If so, will we need something more than just a `MatrixNormal`? Possibly a `KroneckerNormal` or something that takes advantage of `kron_mvprod` or `kron_mmprod` described above when computing the exponential and some determinant tricks for the prefactor.

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [October 26, 2017, 5:14pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/6 "2017-10-26T17:14:49Z")

</div>

Unfortunately theano does not. Would have to implement it.

I think it would be useful, though i bet that a vast majority of cases are covered by MatrixNormal. If youre making a spatial model, you just need x and y. Maybe you could have time too? But past that its harder to anticipate usecases. The amount of work to implement the general version might be more than the utility it brings. Im for sticking with MatrixNormal for now, and we can see if theres demand for the general case?

---

<div class="post-metadata">

**Author:** ![jordan-melendez](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jordan-melendez/32/399_2.png) [@jordan-melendez](https://discourse.pymc.io/u/jordan-melendez)\
**Post date:** [October 27, 2017, 9:32pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/7 "2017-10-27T21:32:40Z")

</div>

You are probably right, and I am not opposed to sticking with `MatrixNormal`. I may be in the 0.1% that could benefit from the generic case, so I may dig deeper and see what I can do, but who knows if anything useful for PyMC3 will come of it. For reference, [this page](https://github.com/jameshensman/VFF/blob/master/VFF/kronecker_ops.py) implements some nice Kronecker tricks for GPflow and also defines a Kronecker GP class, so some subset of this may be of relevance here.

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [October 27, 2017, 9:43pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/8 "2017-10-27T21:43:39Z")

</div>

Oh! nice find. It’s certainly relevant. Ill look more at it. I thought for your case you just needed MatrixNormal? Do you have more stuff coming that would benefit from the generic case?

---

<div class="post-metadata">

**Author:** ![jordan-melendez](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jordan-melendez/32/399_2.png) [@jordan-melendez](https://discourse.pymc.io/u/jordan-melendez)\
**Post date:** [October 27, 2017, 10:01pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/9 "2017-10-27T22:01:59Z")

</div>

So in my case I have multiple realizations of GPs, some of which have 2D inputs. The total matrix is then the Kronecker product of the covariance between curves with the Kronecker product of the two input covariances. Using `MatrixNormal` would work for this if I unpack the `y` into 1D, but the speed increases from the general case would allow me to cover more of the input space. Who knows, it could be premature optimization on my end!

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [October 27, 2017, 10:06pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/10 "2017-10-27T22:06:50Z")

</div>

Are your 2D inputs on a grid? If not then you couldn’t use Kronecker structure anyway. You can just as well do something like $k(x1, x2 ; x1’, x2’) \otimes k(x3, x3’)$, where the x3 input would be indexing your different GP realizations, and the x1, x2 inputs would be the 2D inputs.

---

<div class="post-metadata">

**Author:** ![jordan-melendez](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jordan-melendez/32/399_2.png) [@jordan-melendez](https://discourse.pymc.io/u/jordan-melendez)\
**Post date:** [October 27, 2017, 10:17pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/11 "2017-10-27T22:17:17Z")

</div>

I can choose my inputs to be on a grid and can make them as fine or as course as needed. I essentially have access to continuous functions that are very cheap to sample from. Your suggestion would work, it would just involve inverting a larger matrix than is necessary.

---

<div class="post-metadata">

**Author:** ![jordan-melendez](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jordan-melendez/32/399_2.png) [@jordan-melendez](https://discourse.pymc.io/u/jordan-melendez)\
**Post date:** [November 21, 2017, 10:37pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/12 "2017-11-21T22:37:29Z")

</div>

I am resurrecting this discussion because I have made some progress on these fronts:

- Kronecker functionality is [here](https://github.com/jordan-melendez/pymc3/blob/KroneckerGP/pymc3/math.py#L26)! If Theano wizards can improve on `kron_dot`, etc., I’m all ears.
- Thus, a `KroneckerNormal` is now possible and is also [here](https://github.com/jordan-melendez/pymc3/blob/KroneckerGP/pymc3/distributions/multivariate.py#L1298).
- `MarginalKron` [too](https://github.com/jordan-melendez/pymc3/blob/KroneckerGP/pymc3/gp/gp.py#L797), though it requires testing and is not quite finished (currently getting a weird error when `jobs != 1` which is `AttributeError: Can't pickle local object 'EighGrad. __init__.<locals>.<lambda>'`).
- And `Coregion` kernel is [here](https://github.com/jordan-melendez/pymc3/blob/KroneckerGP/pymc3/gp/cov.py#L521), (shamelessly stolen from GPy) but hasn’t been tested.

Tests and further improvements are on their way.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [November 22, 2017, 6:14am UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/13 "2017-11-22T06:14:25Z")

</div>

> [@jordan-melendez](#):
>
> MarginalKron too, though it requires testing and is not quite finished (currently getting a weird error when jobs != 1 which is `AttributeError: Can't pickle local object 'EighGrad. __init__.&lt;locals&gt;.&lt;lambda&gt;')`.

That’s usually because you have a lambda function in the code and it can not be pickle for njobs, eg see [https://github.com/pymc-devs/pymc3/blob/master/pymc3/gp/gp.py#L640-L642](https://github.com/pymc-devs/pymc3/blob/master/pymc3/gp/gp.py#L640-L642)

cc @bwengals

---

<div class="post-metadata">

**Author:** ![jordan-melendez](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jordan-melendez/32/399_2.png) [@jordan-melendez](https://discourse.pymc.io/u/jordan-melendez)\
**Post date:** [November 22, 2017, 12:16pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/14 "2017-11-22T12:16:35Z")

</div>

I don’t see any lambdas, but it looks like the problem only occurs when using eigenvalue decomposition, which happens inside `KroneckerNormal`. It could be somewhere in [this](https://github.com/jordan-melendez/pymc3/blob/KroneckerGP/pymc3/distributions/multivariate.py#L1359) block, but I haven’t managed to fix it yet.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [November 22, 2017, 12:27pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/15 "2017-11-22T12:27:32Z")

</div>

Yeah that is likely. Submit a PR when you are ready, maybe @aseyboldt or @ferrine know a better way to do that instead of `*map`

---

<div class="post-metadata">

**Author:** ![ferrine](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ferrine/32/10_2.png) [@ferrine](https://discourse.pymc.io/u/ferrine)\
**Post date:** [November 22, 2017, 12:44pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/16 "2017-11-22T12:44:40Z")

</div>

Don’t think problem is here. It is [here](https://github.com/Theano/theano/blob/d395439aec5a6ddde8ef5c266fd976412a5c5695/theano/tensor/nlinalg.py#L378)

---

<div class="post-metadata">

**Author:** ![jordan-melendez](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jordan-melendez/32/399_2.png) [@jordan-melendez](https://discourse.pymc.io/u/jordan-melendez)\
**Post date:** [November 22, 2017, 2:37pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/17 "2017-11-22T14:37:04Z")

</div>

Good catch @ferrine. Eigendecomposition is required to add diagonal noise to the covariance matrix. Could someone submit a PR to Theano to switch to `partial` or would we have to find a workaround on our end?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [November 22, 2017, 4:23pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/18 "2017-11-22T16:23:14Z")

</div>

What we did before is work out a fix locally and submit a PR to theano, eg see here [https://github.com/pymc-devs/pymc3/blob/master/pymc3/distributions/dist\_math.py#L240-L241](https://github.com/pymc-devs/pymc3/blob/master/pymc3/distributions/dist_math.py#L240-L241)

---

<div class="post-metadata">

**Author:** ![jordan-melendez](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jordan-melendez/32/399_2.png) [@jordan-melendez](https://discourse.pymc.io/u/jordan-melendez)\
**Post date:** [November 25, 2017, 2:22pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/19 "2017-11-25T14:22:47Z")

</div>

PR for Kronecker stuff started [here](https://github.com/pymc-devs/pymc3/pull/2731).

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [November 25, 2017, 9:58pm UTC](https://discourse.pymc.io/t/extending-gaussian-process-functionality-coregion-and-beyond/479/20 "2017-11-25T21:58:43Z")

</div>

I’m really looking forward to seeing this stuff in action. There’ve been several questions that have come up on discourse where this sort of stuff would have helped people out.

@junpenglao that pickling fix replacing the lambda with functools.partial was made by @aseyboldt, I believe.