Behaviour of Gibbs kernel seems off

I am trying to draw samples from a GP prior with Gibbs covariance - following the example in the API on the website if I just were to replace the tanh_func with a function of my own choosing (in this case I am sampling the warping function from another GP) – it is very clear that the resulting covariance K is not positive semi-definite.

(I think it might be because if you zoom into the plot for the warping function there are these feeble discontinuities – I am not sure why they arise as a monotonic transformation of a perfectly smooth function shouldn’t result in any discontinuities. )

import numpy as np
import pymc3 as pm 
import matplotlib.pyplot as plt
import theano.tensor as tt

X = np.linspace(-10, 10, 1000)[:,None]
y = (X**3 - X**2 + X)

lengthscale = 2
eta = 2.0
cov = eta**2 *, lengthscale)
K_base = cov(X).eval()

###### Samples from Gibbs covariance

def warp_gp_func(x, K):
    sample = pm.MvNormal.dist(mu=np.zeros(K.shape[0]), cov=K + np.eye(1000)*1e-5, shape=(len(X),)).random(size=1).T*tt.ones(shape=(1000,1))
    return sample.exp()

def tanh_func(x, ls1, ls2, w, x0):
    ls1: left saturation value
    ls2: right saturation value
    w:   transition width
    x0:  transition location.
    return (ls1 + ls2) / 2.0 - (ls1 - ls2) / 2.0 * tt.tanh((x - x0) / w)

ls1 = 0.05
ls2 = 0.6
w = 0.3
x0 = 1.0
#cov =, tanh_func, args=(ls1, ls2, w, x0))
cov =, warp_gp_func, args=(K_base))

#th = tanh_func(X, ls1, ls2, w, x0).eval()
wf = warp_gp_func(X, K_base).eval()
plt.figure(figsize=(14, 4))
plt.plot(X, wf)
plt.title("Lengthscale as a function of X")

K = cov(X).eval()

### Below block throws an error

plt.figure(figsize=(14, 4))
plt.plot(X, pm.MvNormal.dist(mu=np.zeros(K.shape[0]), cov=K + np.eye(1000)*1e-6, shape=(1000,)).random(size=3).T)
plt.title("Samples from the GP prior")


I’m not exactly sure what the problem is, but it seems like it might be coming from the warp function and not from the implementation of the Gibbs kernel?

One thing I tried was seeing if instead of exponentiating your samples if taking the inverse logistic would work, but this gave me the same problem. I have a sneaking suspicion it might be caused by exponentials amplifying errors somewhere.

While PyMC3’s pm.MvNormal will fail quite loudly if you try to pass a non positive semi-definite matrix as a covariance, I’ve found that NumPy’s np.random.multivariate_normal will still produce samples but will just warn you. I’m not sure if the mathematics breaks in a meaningful way This definitely breaks the mathematics but if I replace your last block of code with the following I’m able to draw samples:

plt.figure(figsize=(14, 4))
plt.plot(X, np.random.multivariate_normal(mean=np.zeros(K.shape[0]), cov=K + np.eye(1000)*1e-9, size=3).T)
plt.title("Samples from the GP prior")

Hope this helps!

1 Like

I wasn’t able to figure out why your example as written didn’t work. It looks like this is using exp(GP) for the lengthscale of another GP. Does this give something sensible for samples from the GP prior?:

import numpy as np
import pymc3 as pm 
import matplotlib.pyplot as plt
import theano.tensor as tt

X = np.linspace(-10, 10, 300)[:,None]  # i changed the 1000 to 300 to speed it up here

with pm.Model() as model:
    eta = 2.0
    lengthscale = 2.0
    base_cov = eta**2 *, lengthscale)
    base_gp =
    log_wf = base_gp.prior("log_wf", X=X)
    wf = pm.Deterministic("wf", pm.math.exp(log_wf))
    def warp_gp_func(x, wf):
        return wf
    cov =, warp_gp_func, args=(wf,))
    K = pm.Deterministic("K", cov(X))
    gp =
    f = gp.prior("f", X=X)

with model:
    tr = pm.sample(tune=1000, draws=100, cores=1)

The Gibbs covariance matrices here at least look symmetric. Is this kind of what you’re after?

1 Like

thank you, the behaviour of np.random.multivariate_normal is a bit worrying because even if I can sample …the fact that the matrix isn’t psd might break something downstream. :frowning:

thank you, this seems reasonable!