posterior_samples hangs when using multioutput GPs with input dimension > 1

Hello,

I am trying to draw posterior samples from a multi output GP which has a two dimensional input and a two dimensional output. I can call predict() on the trained model just fine, but it appears that posterior_samples() hangs (it never returns), even if I'm requesting one sample only. If the input has dimension 1, the model works fine.

I'm using GPy version 1.10.0 on Python 3.8.10.

Here is a code that demos the issue (the last line is the problem), I hope I'm not doing something wrong.

import numpy as np
import numpy.random as rng
import GPy
import itertools

N = 50

# generate synthetic data with 2d inputs and 2d outputs
X = np.random.rand(N,2)
k = GPy.kern.RBF(1, lengthscale=1)
f1 = rng.multivariate_normal(np.zeros(N), k.K(X)).reshape(N,1)
f2 = rng.multivariate_normal(np.zeros(N), k.K(X)).reshape(N,1)
f = np.hstack((f1, f2))
Y = f + rng.randn(*f.shape) * 0.1

# use icm multioutput gp
n_inputs = 2
n_outputs = 2

icm_kernel = GPy.util.multioutput.ICM(input_dim=n_inputs,
                                      num_outputs=n_outputs, 
                                      kernel=GPy.kern.RBF(n_inputs))

m = GPy.models.GPCoregionalizedRegression([X for i in range(n_outputs)], 
                                          [Y[:,[i]] for i in range(n_outputs)],
                                           kernel=icm_kernel)

m.optimize()

# put the test points in the right format for the multioutput gp
Xstar = np.array(list(itertools.product(np.linspace(-3, 3, 100), np.linspace(-3, 3, 100))))
ones = np.ones(Xstar.shape[0])[:,np.newaxis]
stacked_Xstar = np.vstack([np.hstack([Xstar,ones*i]) for i in range(n_outputs)])
noise_dict = {'output_index': stacked_Xstar[:,[-1]].astype(int)}

# predict and reshape outputs at test points -- works just fine
mu, var = m.predict(stacked_Xstar, Y_metadata=noise_dict)
mu = mu.reshape(n_outputs, Xstar.shape[0]).T

# draw a single posterior sample ... doesn't work when input dimension > 1 :(
samples = m.posterior_samples(stacked_Xstar, Y_metadata=noise_dict, size=1)

SheffieldML / GPy

posterior_samples hangs when using multioutput GPs with input dimension > 1 #973