Additive Gaussian Processes

Questions
1

Hi all,

I was playing with additive Gaussian processes similar to documentation here:
https://docs.pymc.io/Gaussian_Processes.html?highlight=gaussian%20processes

but I am having issues while sampling from the conditional distribution for f1_star and f2_star.

My notebook can be found here:

with pm.Model() as model:
    sigma1 = pm.HalfCauchy("η1", beta=5)
    sigma2 = pm.HalfCauchy("η2", beta=5)
    ls1 = pm.HalfCauchy('ls1', beta=10)
    ls2 = pm.HalfCauchy('ls2', beta=10)
    cov1 = sigma1 * pm.gp.cov.ExpQuad(2, ls=ls1, active_dims=[0])
    cov2 = sigma2 * pm.gp.cov.ExpQuad(2, ls=ls2, active_dims=[1]) 
    gp1 = pm.gp.Marginal(cov1)
    gp2 = pm.gp.Marginal(cov2)
    gp = gp1 + gp2
    sigma = pm.HalfCauchy("sigma", beta=3)
    f = gp.marginal_likelihood("f", X_2D_train, y, sigma)
    trace = pm.sample(1000, chains=4)

with model:
    f_star = gp.conditional("f_star", X_2D_star)

TypeError: For compute_test_value, one input test value does not have the requested type.

The error when converting the test value to that variable type:
Wrong number of dimensions: expected 1, got 2 with shape (1000, 1000).

X_2D_star has shape (1000,2). So, I think I am using additive GP wrong. Can someone please help me with this?

Thanks,
Anurag

2

Hi all,

I wonder if anyone has some tips on how to do additive GP with pymc3 and sample from the conditional distributions.

Thanks,
Anurag

3

What happens if you use gp.predict? Does the following work?

with model:
    pred = gp.predict(X_2D_star)
4

Thanks for your reply @ckrapu.

It works as long as the shape of X_2D_star is the same as the shape in X_2D_train. However, that may not be the case in the most general case. Also, the output is a list with 2 arrays of shape (1000, 1000). How do I interpret these 2 arrays?

Thanks,
Anurag

5

I’m a little stuck here, though I think this thread could be very helpful. @bwengals, any idea? Here’s a minimal reproducible example

With regard to the output of gp.predict, the first output should be the predictive mean vector and the output is the covariance matrix.

1 Like
6

Thanks @ckrapu.

Thats how I interpreted it but was confused because both elements of the list had shape (1000, 1000) in my case. Shouldn’t the mean be a vector of (1000, 1)?

Cheers,
Anurag

7

Apologies, I missed that detail. There is clearly something fishy going on with the broadcasting as the predictive mean should only have one nontrivial dimension.

8

Yes it should, something is fishy here. I can find some time to look into this.

2 Likes
9

Sorry about the delay. So I ran through the notebook and got it to run by specifying the cov_func argument.

gp1 = pm.gp.Marginal(cov_func=cov1)
gp2 = pm.gp.Marginal(cov_func=cov2)

So this is a bit of a gotcha, but if you check the docstring:

Init signature:
pm.gp.Marginal(
    mean_func=<pymc3.gp.mean.Zero object at 0x7fa7d8182130>,
    cov_func=<pymc3.gp.cov.Constant object at 0x7fa7d8182160>,
)

so the mean function was set with the cov function.

2 Likes