Gaussian process regression (RBF)

Background

A Gaussian Process is a distribution over functions: rather than fitting parameters, it places a prior over functions through a kernel (covariance) and conditions on observed data to obtain a posterior. GP regression predicts, at any test input, a Gaussian with a mean and a variance — giving not just a value but calibrated uncertainty. It is a go-to model for small-data regression and Bayesian optimisation.

Problem statement

Implement gaussian_process_predict(X_train, y_train, X_test, length_scale, sigma, noise) that returns the GP posterior mean and standard deviation at the test points, using the RBF (squared-exponential) kernel:

$$k(x, x') = \sigma^2 \exp\!\Big(-\frac{\lVert x - x'\rVert^2}{2\,\ell^2}\Big)$$

With $K = k(X, X) + \text{noise}\cdot I$, $\;K_* = k(X, X_*)$, $\;K_{**} = k(X_*, X_*)$:

$$\mu_* = K_*^\top K^{-1} y, \qquad \Sigma_* = K_{**} - K_*^\top K^{-1} K_*$$

Return $\mu_*$ and the per-point standard deviations $\sqrt{\operatorname{diag}(\Sigma_*)}$.

Input

• X_train — np.ndarray (n, d): training inputs.
• y_train — np.ndarray (n,): training targets.
• X_test — np.ndarray (m, d): test inputs.
• length_scale — float: RBF length scale $\ell$.
• sigma — float: kernel output scale.
• noise — float: observation-noise variance added to the diagonal of $K$.

Output

Returns a tuple (mu, std): mu is np.ndarray (m,) of posterior means and std is np.ndarray (m,) of posterior standard deviations.

Examples

Example 1

Input:  X_train=[[0],[2],[4]], y_train=[0,1,0], X_test=[[0],[100]],
        length_scale=1.0, sigma=1.0, noise=1e-8
Output: mu ~= [0.0, 0.0], std ~= [~0.0, 1.0]

Explanation: at $x=0$ (a training point) the mean equals the observed $y=0$ with near-zero uncertainty; at $x=100$, far from all data, the RBF covariance vanishes, so the prediction reverts to the prior mean $0$ with the prior std $\sigma = 1$.

Constraints

• Use the RBF kernel; add noise to the diagonal of the train–train covariance before inverting.
• Posterior mean is $\mu_* = K_*^\top K^{-1} y$; the variance is the diagonal of $K_{**} - K_*^\top K^{-1} K_*$.
• Return per-point standard deviations (square root of the variance, clipped at $0$).
• Pure numpy (np.linalg.inv / np.linalg.solve); tests compare with atol=1e-3.

Notes

• As noise -> 0 the GP interpolates the training data exactly (the mean passes through every observation); larger noise smooths the fit.
• The length scale $\ell$ controls wiggliness: small $\ell$ gives fast-varying functions, large $\ell$ gives smooth, slowly varying ones.