Error in Gaussian Process Implementation

I have been trying to define a Gaussian Process from hand by simply writing down the required formulas and apply them to an toy example. Here is what my (noise-free), observed data looks like:

x_observed = np.arange(-5, 5, 0.2).reshape(-1, 1)
y_observed = np.sin(x_train)
fig, ax = plt.subplots(figsize=(12, 8))
ax.scatter(x_train, y)


Moreover, I used the exponential quadratic kernel function:

def kernel_func(X1, X2, l=1.0, sigma=1.0):

sqdist = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)
return sigma**2 * np.exp(-0.5 / l**2 * sqdist)


and defined the gaussian process:

def gp(x_observed, y_observed, x_new):

k = kernel_func(x_observed, x_observed)
K_star_star = kernel_func(x_new, x_new)
K_star = kernel_func(x_observed, x_new)
K_inv = np.linalg.inv(k)
mu = K_star.T @ K_inv @ y_observed
sigma = K_star_star - K_star.T @ K_inv @ K_star

return mu, sigma


Creating new data, which I want a prediction for:

X_new = np.array([-3.5, -3, -2, -1, 1]).reshape(-1, 1)


I can get the mean value for each $$x_i$$ in X_new, by calling the gaussian process function:

mu_s, cov_s = gp(x_observed, y_observed, X_new)


Now, I did expect that the point prediction, that is the mean, for some of the $$x_i$$ in X_new are equal to some values in y_observed, since we already know the target value, i.e. the values were already observed. However, looking at the observed y values and mu_s :

for i, x in  enumerate(X_new):
print(f'target function {np.sin(x)}, prediction {mu_s[i]}')

target function [0.35078323], prediction [-6.27205575]
target function [-0.14112001], prediction [-1.03358472]
target function [-0.90929743], prediction [2.2573257]
target function [-0.84147098], prediction [3.20067827]
target function [0.84147098], prediction [-7.19285408]


I either have a missunderstanding of Gaussian Processes or an error in my code, but I could not find any bugs or similar

• aren't you missing a transpose in the equation for mu ? (K_star.T instead of K_star) May 6, 2023 at 10:20
• I think this depends on the transpose used in the kernel function. I could also add the transpose, like you suggest, but then I have to adjust how I pass the input of the kernel function. The result should be the same.I will adjust it in the code though to avoid confusion. May 6, 2023 at 10:30
• Right. Anyway your understanding is correct so this must be a small bug somewhere. I would suggest testing the code with fewer observations and plotting the results, this might give you a hint about where things go wrong May 6, 2023 at 10:40

The problem is the ill conditioning of your kernel matrix. These are the singular values of your kernel matrix k:

As you can see, many of them are numerically zero. This leads to nonsense when you compute np.linalg.inv.

You have two options. The most common is to simply add a scaled identity matrix to your kernel matrix with some small scaling value, like $$10^{-6}$$:

K_inv = np.linalg.inv(k + 1e-6 * np.eye(k.shape[0]))


Another option is to notice that each time you explicitly compute a matrix inverse, Householder rolls around in his grave. We can allow him to rest easier by avoiding explicit inverse computation, and instead computing linear system solves. Like so:

def gp(x_observed, y_observed, x_new):

k = kernel_func(x_observed, x_observed)
K_star_star = kernel_func(x_new, x_new)
K_star = kernel_func(x_observed, x_new)
mu = K_star.T @ np.linalg.solve(k, y_observed)
sigma = K_star_star - K_star.T @ np.linalg.solve(k, K_star)

return mu, sigma


This also leads to good predictions:

• Great answer! I'm surprised numpy does not issue a warning when trying to inverse close to singular matrices May 6, 2023 at 16:13
• you saved my weekend, ty! May 6, 2023 at 16:30
• @J.Delaney - the issue is that "close" is in the eye of the beholder. Any value you picked would simultaneously get some people annoyed at nuisance warnings and others thinking they had a (nonexistent) safety net because they didn't get warnings.
– TLW
May 6, 2023 at 18:20

@John Madden gave a good answer pointing to the root cause of the problem, but adding to it, you should not invert that matrix directly in the first place. Matrix inversion is generally inefficient and not recommended for all kinds of applications. This also applies to Gaussian processes. An efficient algorithm is given by Rasmussen (2006) as Algorithm 2.1 on p 19:

\begin{align} L &= \operatorname{cholesky}(K + \sigma^2 I) \\ \alpha &= L^\top \backslash (L \backslash y) \\ \mu &= K_{*}^\top \alpha \\ v &= L \backslash K_{*} \\ \Sigma^2 &= K_{**} - v^\top v \\ \end{align}

where $$K = k(X, X)$$, $$K_{*} = k(X, X_*)$$ and $$K_{**} = k(X_*, X_*)$$.

• may be worth emphasizing that these are not any linear system solves that youre invoking, but lower diagonal ones which have complexity only quadratic in matrix dimension. BTW: do you think this approach is superior to the direct linear system solves I recommend? I've never actually compared the two, I could see it going either way... but anyways yes it's true this is the most popular approach, especially because determinant calculation is linear time once you've done a cholesky decomposition (in case you want to do marginal MLE). May 8, 2023 at 13:49