# Multiple Gaussian processes without for loop?

**URL:** https://discourse.pymc.io/t/multiple-gaussian-processes-without-for-loop/8606
**Category:** Questions
**Tags:** gaussian\_process
**Created:** [January 17, 2022, 11:58am UTC](https://discourse.pymc.io/t/multiple-gaussian-processes-without-for-loop/8606 "2022-01-17T11:58:58Z")
**Posts on this page:** 7
**Page:** 1

<div class="post-metadata">

### Author: ![ysfoo](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ysfoo/32/7303_2.png) [@ysfoo](https://discourse.pymc.io/u/ysfoo)
#### Post date: [January 17, 2022, 11:58am UTC](https://discourse.pymc.io/t/multiple-gaussian-processes-without-for-loop/8606/1 "2022-01-17T11:58:58Z")

</div>

I have a model involving a multinomial likelihood with multiple GPs as my prior (the GPs should have different covariance hyperparameters). I’ve heard that for loops are usually frowned upon when defining a model. Is there a way to define multiple GPs while avoiding a for loop? How ‘bad’ is it if it can’t be avoided?

Edit:

```auto
import numpy as np
import pymc3 as pm
import theano.tensor as tt

N = 20
H = 8
# assume X is N x 1 array and y is N x H array of multinomial counts
X = np.arange(N)
y = np.ones((N,H))

with pm.Model() as model:
    ls = pm.Lognormal('ls', mu=1, sigma=0.5, shape=H)
    gps = [pm.gp.Latent(cov_func=pm.gp.cov.ExpQuad(1, ls[i])) for i in range(H)]
    p = tt.nnet.softmax(tt.stack([gp.prior('f'+str(i), X=X[:,None])
                                  for i, gp in enumerate(gps)], axis=1))
    counts = pm.Multinomial('y', p=p, n=np.sum(y, axis=1), shape=(N,H))

```

---

<div class="post-metadata">

### Author: ![ericmjl](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ericmjl/32/196_2.png) [@ericmjl](https://discourse.pymc.io/u/ericmjl)
#### Post date: [January 17, 2022, 3:52pm UTC](https://discourse.pymc.io/t/multiple-gaussian-processes-without-for-loop/8606/2 "2022-01-17T15:52:00Z")

</div>

Hi @ysfoo! Thanks for chiming in, and welcome!

Given the nature of us nerds here, I think it helps the community to answer your question if you _show_ us the kind of code you’re trying to avoid. Would you be open to editing your post with a minimum workable example that can be copy/pasted and executed in a Jupyter notebook? That will really help us grok what you’re trying to do.

---

<div class="post-metadata">

### Author: ![ysfoo](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ysfoo/32/7303_2.png) [@ysfoo](https://discourse.pymc.io/u/ysfoo)
#### Post date: [January 18, 2022, 9:11am UTC](https://discourse.pymc.io/t/multiple-gaussian-processes-without-for-loop/8606/3 "2022-01-18T09:11:26Z")

</div>

Thanks for the welcome! I’ve added a MWE now.

---

<div class="post-metadata">

### Author: ![twiecki](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/twiecki/32/6930_2.png) [@twiecki](https://discourse.pymc.io/u/twiecki)
#### Post date: [January 19, 2022, 6:35am UTC](https://discourse.pymc.io/t/multiple-gaussian-processes-without-for-loop/8606/4 "2022-01-19T06:35:09Z")

</div>

You can do multi-dimensional GPs quite easily: [Mean and Covariance Functions — PyMC3 3.11.4 documentation](https://docs.pymc.io/en/v3/pymc-examples/examples/gaussian_processes/GP-MeansAndCovs.html#Two-(and-higher)-Dimensional-Inputs)

---

<div class="post-metadata">

### Author: ![ysfoo](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ysfoo/32/7303_2.png) [@ysfoo](https://discourse.pymc.io/u/ysfoo)
#### Post date: [January 19, 2022, 6:50am UTC](https://discourse.pymc.io/t/multiple-gaussian-processes-without-for-loop/8606/5 "2022-01-19T06:50:49Z")

</div>

Hmm, I don’t think that’s quite what I’m looking for. That shows multi-dimensional inputs, but my inputs are in one dimension, with multiple outputs.

I also know that PyMC implements Kronecker-structure covariance, but I don’t think I can have multiple length scales in that case (like the MWE).

---

<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 19, 2022, 7:44am UTC](https://discourse.pymc.io/t/multiple-gaussian-processes-without-for-loop/8606/6 "2022-01-19T07:44:44Z")

</div>

There’s nothing too bad with that loop, you are just automating the creation of 8 separate variables. Loops are generally discouraged in the case where you can use a single variable tensor instead of 100s individual variables, but your case might fall more in the “premature optimization” bucket.

---

<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: [January 19, 2022, 9:14pm UTC](https://discourse.pymc.io/t/multiple-gaussian-processes-without-for-loop/8606/7 "2022-01-19T21:14:39Z")

</div>

I think the way you’re doing this with the for loop is exactly right here. Since your covariances don’t share hyperparameters, there aren’t any efficiency tricks possible, so automating the creation of multiple independent GPs with a for loop is 💯
