Sequential sampling from Gaussian conditional not working

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        
    else:
        # 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.