Multiple Gaussian processes without for loop?

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?


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 = [, 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))

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.


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

1 Like

You can do multi-dimensional GPs quite easily: Mean and Covariance Functions — PyMC3 3.11.4 documentation

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).

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.


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 :100: