Thoroughly confused - I have a dataset w/ n=564 observations, p=12 features, and a single target.

Applying a Linear kernel as such works fine (p=2 for simplicity and debugging at the moment):

```
features = FEATURES[0:2]
with pm.Model() as m5:
eta_2 = pm.InverseGamma("eta2", 1, 1)
c = pm.Normal("c", 0, 1, shape=2)
cov_se = eta_2**2 * pm.gp.cov.Linear(len(features), c=c)
gp = pm.gp.Marginal(cov_func=cov_se)
# noise
sigma = pm.InverseGamma("sigma", 1, 1)
y_ = gp.marginal_likelihood("y", X=X[:,0:2], y=y, sigma=sigma)
with m5:
mp = pm.find_MAP()
mp
```

However, instantiating a Polynomial kernel gives the following error: `ValueError: Incompatible Elemwise input shapes [(564, 564), (1, 2)]`

```
features = FEATURES[0:2]
with pm.Model() as m6:
eta_2 = pm.InverseGamma("eta2", 1, 1)
c = pm.Normal("c", 0, 1, shape=2)
offset = pm.Normal("offset", 0, 1, shape=2)
d = pm.Normal("d", 0, 1, shape=2)
cov_se = eta_2**2 * pm.gp.cov.Polynomial(input_dim=len(features), c=c, offset=offset, d=d)
gp = pm.gp.Marginal(cov_func=cov_se)
# noise
sigma = pm.InverseGamma("sigma", 1, 1)
y_ = gp.marginal_likelihood("y", X=X[:,0:2], y=y, sigma=sigma)
```

PyMC version 5.6.0