Dimension of GPy.core.gp.posterior_samples_f()

[lawrennd]

Just using this function in some deep GP sampling.

First, I think there's a couple of bugs. On line 567, do we really want self.num_data here? Shouldn't it be X.shape[1].

But further, the defined behaviour of this function seems pretty nasty. It returns a dx n x size ndarray. Unless d==1 when it returns an n x size array.

For example, see the code below, where d is output_dim

import numpy as np
import GPy
out_dim = 3
x = np.random.normal(size=(10, 2))
y = np.random.normal(size=(10, out_dim))
model = GPy.models.GPRegression(x, y)
xt = np.random.normal(size=(12, 2))
y_samp = model.posterior_samples(xt)
print("Output Dim: ", out_dim)
print("model.posterior_sample dimension: ", y_samp.shape)

out_dim=1
x = np.random.normal(size=(10, 2))
y = np.random.normal(size=(10, out_dim))
model = GPy.models.GPRegression(x, y)
xt = np.random.normal(size=(12, 2))
y_samp = model.posterior_samples(xt)
print("Output Dim: ", out_dim)
print("model.posterior_sample dimension: ", y_samp.shape)

I'm pretty uncomfortable with this behaviour. It's very difficult to handle downstream.

As standard in GPy we set the first dimension of any matrix to the number of data. I think this function should return an n x d x size array.

And the singleton dimension of d should never be squashed. Perhaps for size=1 then that should be squashed. But actually, that's still ugly.

In my local fork, I'll make this change, but in the meantime it would be good to hear about potential downstream affects. I'm sure there are some in plotting functions, but what about elsewhere?

[lawrennd]

As a follow up, does it make any sense to sample this without the full covariance? I don't think so, but there may be a use case I'm missing.