I'm trying to sequentially sample from a Gaussian Process prior.

The problem is that the samples eventually converge to zero or diverge to infinity.

I'm using the basic conditionals described e.g. here

Note: the kernel(X,X) function returns the squared exponential kernel with isometric noise.

Here is my code:

n = 32

x_grid = np.linspace(-5,5,n)

x_all = []
y_all = []
for x in x_grid:
    x_all = [x] + x_all
    X = np.array(x_all).reshape(-1, 1)
    # Mean and covariance of the prior
    mu = np.zeros((X.shape), np.float)
    cov = kernel(X, X)
    if len(mu)==1: # first sample is not conditional
        y = np.random.randn()*cov + mu        
        # condition on all previous samples
        u1 = mu[0]
        u2 = mu[1:]
        y2 = np.atleast_2d(np.array(y_all)).T
        C11 = cov[:1,:1] # dependent sample
        C12 = np.atleast_2d(cov[0,1:])
        C21 = np.atleast_2d(cov[1:,0]).T
        C22 = np.atleast_2d(cov[1:, 1:])
        C22_ = la.inv(C22)
        u = u1 + np.dot(C12, np.dot(C22_, (y2 - u2)))
        C22_xC21 = np.dot(C22_, C21)
        C_minus = np.dot(C12, C22_xC21) # this weirdly becomes larger than C!
        C = C11 - C_minus        
        y = u + np.random.randn()*C
    y_all = [y.flatten()[0]] + y_all

Here's an example with 32 samples, where it collapses:

enter image description here

Here's an example with 34 samples, where it explodes:

enter image description here

(for this particular kernel, 34 is the number of samples at which (or more) the samples start to diverge.

I've gone through this code so many times that I'm going blind - something must be fundamentally wrong with it but I just can't see it.

  • 1
    $\begingroup$ Can you post all the code so the result can be replicated? $\endgroup$ Jan 16 at 14:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.