First two things to try: would definitely recommend using the default NUTS instead of Metropolis, and set cores=1 in pm.sample whenever using the GP submodule.
Also, the model @kejzlarv uses looks like:
gp1 = pm.gp.Marginal(mean_func = mean_fcn, cov_func = cov_fcn)
gp2 = pm.gp.Marginal(mean_func = mean_fcn, cov_func = cov_fcn)
y_1 = gp1.marginal_likelihood("y1", X=input_1, y=y1, noise=σ)
y_2 = gp2.marginal_likelihood("y2", X=input_2, y=y2, noise=σ)
which is a little different than what you have. Maybe these things help?