Gaussian Processes and Bayesian Regression

Gaussian process regression describes a distribution over possible functions, then updates it using observations. Predictions include a mean, a variance, and correlations between predictions at different inputs. Those uncertainties depend on the prior and observation model; they are not an automatic guarantee that intervals will be well calibrated on new data. This article works through Gaussian-noise regression, where the update has a closed form, and examines how the kernel, fitted hyperparameters, and computational budget affect it.

A prior over functions

A Gaussian process (GP) specifies that the function values at any finite set of inputs are jointly Gaussian. Its mean function \(m_0(x)\) describes the prior center, and its kernel gives \(\operatorname{Cov}(f(x),f(x’))=k(x,x’)\). For an RBF kernel, nearby function values are strongly correlated; a periodic kernel can link values one cycle apart. Bayesian linear regression is a finite-feature example: if \(f(x)=\phi(x)^\top\beta\) and \(\beta\sim\mathcal N(0,\Sigma_\beta)\), then \(k(x,x’)=\phi(x)^\top\Sigma_\beta\phi(x’)\). Here \(\phi(x)\) contains chosen features, such as an intercept and the original input, and \(\Sigma_\beta\) is the prior covariance of their coefficients.

Assume observations \(y_i=f(x_i)+\varepsilon_i\), with independent zero-mean Gaussian noise of variance \(\sigma_n^2\), independent of the function. For fixed kernel and noise parameters, define \(K=k(X,X)+\sigma_n^2I\), an \(n\times n\) noisy training covariance. For \(q\) query points, \(K_*=k(X,X_*)\) has shape \(n\times q\), and \(K_{**}=k(X_*,X_*)\) has shape \(q\times q\). The posterior mean is \(m_0(X_*)+K_*^\top K^{-1}(y-m_0(X))\), and the latent-function covariance is \(K_{**}-K_*^\top K^{-1}K_*\). Non-Gaussian likelihoods generally require another inference method.

For one observation with zero prior mean, suppose its noisy variance is 1.25, the query’s prior variance is 1, and their covariance is 0.8. Observing \(y=1\) gives query mean \(0.8/1.25=0.64\) and variance \(1-0.8^2/1.25=0.488\). Correlation transfers information from the observation and reduces uncertainty. The code below performs this update for 12 observations and seven queries with \(m_0=0\). Its kernel helper is specific to one-dimensional inputs. Run the two body code blocks in order.

import numpy as np
from scipy.linalg import solve_triangular
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel

rng = np.random.default_rng(0)
X_tr = rng.uniform(-3, 3, (12, 1))
y_tr = np.sin(X_tr.ravel()) + 0.1 * rng.normal(size=12)
X_te = np.linspace(-5, 5, 7).reshape(-1, 1)
ell, sf2, sn2 = 1.0, 1.0, 0.01

k = lambda A, B: sf2 * np.exp(-0.5 * ((A - B.T) ** 2) / ell ** 2)
K = k(X_tr, X_tr) + sn2 * np.eye(len(X_tr))
K_s, K_ss = k(X_tr, X_te), k(X_te, X_te)

L = np.linalg.cholesky(K)                      # factor the noisy covariance once
alpha = solve_triangular(L.T, solve_triangular(L, y_tr, lower=True), lower=False)
mu = K_s.T @ alpha
v = solve_triangular(L, K_s, lower=True)
sd = np.sqrt(np.diag(K_ss) - np.sum(v ** 2, axis=0))

gp = GaussianProcessRegressor(
    kernel=ConstantKernel(sf2, "fixed") * RBF(ell, "fixed"),
    alpha=sn2, optimizer=None).fit(X_tr, y_tr)
mu_sk, sd_sk = gp.predict(X_te, return_std=True)

print("     x    by hand mu   sklearn mu   by hand sd   sklearn sd")
for i in range(7):
    print(f"  {X_te[i, 0]:5.2f} {mu[i]:12.6f} {mu_sk[i]:12.6f}"
          f" {sd[i]:12.6f} {sd_sk[i]:12.6f}")
#      x    by hand mu   sklearn mu   by hand sd   sklearn sd
#   -5.00     0.208036     0.208036     0.974594     0.974594
#   -3.33     0.204936     0.204936     0.258724     0.258724
#   -1.67    -1.130532    -1.130532     0.164573     0.164573
#    0.00     0.065411     0.065411     0.157043     0.157043
#    1.67     0.977530     0.977530     0.064665     0.064665
#    3.33     0.478568     0.478568     0.470942     0.470942
#    5.00     0.091014     0.091014     0.994337     0.994337

The manual calculation and scikit-learn agree to six decimals. The reported standard deviation is about 0.065 at \(x=1.67\) and 0.994 at \(x=5\). These are uncertainties about the latent function \(f(x)\), conditional on the fixed parameters. For a new noisy observation with independent variance 0.01, add 0.01 to the latent variance before taking the square root; the first value then becomes about 0.119. In this model, alpha=sn2 adds noise to the training covariance but predict(return_std=True) does not add that noise to the query variance. A WhiteKernel component contributes to the query variance as part of the kernel, so the two API choices require different interpretation.

The code factors \(K=LL^\top\) once and uses triangular solves to apply \(K^{-1}\) without forming the inverse. A valid kernel matrix is positive semidefinite; the strictly positive noise variance here makes \(K\) positive definite in exact arithmetic. Roundoff and poorly conditioned parameters can still cause numerical problems. Cholesky exploits symmetry and costs roughly half the leading factorization work of general LU, but that is not a universal twofold wall-clock speedup. The same factor gives \(\log\det K=2\sum_i\log L_{ii}\), used below for the marginal likelihood.

Where the uncertainty goes

import numpy as np

X_q = np.array([[0.0], [2.0], [4.0], [8.0], [20.0]])
K_q = k(X_tr, X_q)
v_q = solve_triangular(L, K_q, lower=True)
mu_q = K_q.T @ alpha
sd_q = np.sqrt(np.diag(k(X_q, X_q)) - np.sum(v_q ** 2, axis=0))

for i in range(len(X_q)):
    nearest = np.abs(X_tr.ravel() - X_q[i, 0]).min()
    print(f"x={X_q[i, 0]:5.1f}  nearest training point {nearest:5.2f} away"
          f"   mean {mu_q[i]:+8.4f}  sd {sd_q[i]:.4f}")
print(f"prior sd {np.sqrt(sf2):.4f}")
# x=  0.0  nearest training point  0.26 away   mean  +0.0654  sd 0.1570
# x=  2.0  nearest training point  0.10 away   mean  +0.8626  sd 0.0672
# x=  4.0  nearest training point  1.39 away   mean  +0.3601  sd 0.8419
# x=  8.0  nearest training point  5.39 away   mean  +0.0000  sd 1.0000
# x= 20.0  nearest training point 17.39 away   mean  +0.0000  sd 1.0000
# prior sd 1.0000

For this zero-mean RBF model, predictions approach the prior as the query becomes far from all training inputs: the cross-covariances approach zero, so the mean approaches zero and the latent standard deviation approaches 1. The displayed values at 8 and 20 are rounded. At five length scales the RBF correlation is \(e^{-25/2}\approx3.7\times10^{-6}\), small but nonzero; there is no exact finite-distance cutoff.

This reversion is a property of a decaying kernel with the stated mean, not of every GP. A periodic kernel can remain strongly correlated with distant observations, and a linear kernel can continue a trend. Bayesian linear regression also supplies predictive uncertainty, so uncertainty estimation is not unique to GPs. At fixed GP hyperparameters, the posterior covariance depends on input locations and noise levels, not on the observed target values. A surprising target does not by itself widen these bands; refitting hyperparameters can change them.

A zero prior mean can be appropriate after centering or modeling residuals around a baseline. With a specified nonzero mean, the update uses residuals \(y-m_0(X)\) and adds the mean back at the query points. Fitting that baseline introduces uncertainty of its own. The mean, kernel, noise assumptions, and treatment of their parameters jointly determine what extrapolation says. Pointwise intervals such as mean ±1.96 standard deviations have a Gaussian posterior interpretation under this model; they are not simultaneous bands for the whole curve or automatic frequentist coverage guarantees.

Choosing a covariance structure

Kernel	Covariance assumption	Possible use
RBF	Mean-square derivatives of all orders	Very smooth variation
Matérn 3/2 or 5/2	One or two mean-square derivatives	Rougher variation than RBF
Periodic	Exact repetition at its specified period	A stable repeating cycle
Linear	Gaussian uncertainty over linear coefficients	Bayesian linear regression
Sum	Additive independent components	Trend plus a cycle
Periodic × RBF	Correlation between cycles decays with separation	A cycle whose shape evolves

The smoothness in the table is mean-square differentiability: derivatives exist as limits in mean squared error under the GP prior. RBF permits derivatives of every order; Matérn 3/2 and 5/2 permit one and two respectively. These are different structural assumptions, not a ranking of kernels. RBF can work well for smooth signals, while a Matérn kernel may better represent rougher variation. Compare plausible kernels using predictions and uncertainty checks that match the application.

Sums of kernels describe sums of independent GP components, such as a trend plus a periodic signal. Products give another valid covariance: multiplying a periodic kernel by an RBF kernel weakens correlation between matching phases as cycles get farther apart. That permits the cycle’s shape to evolve. With constant component variances, the product still has constant marginal variance, so it does not specifically model a drifting amplitude envelope. Also, multiplying two GP sample functions does not generally produce a Gaussian process. The product here is a covariance construction.

Fitting the hyperparameters

Kernel hyperparameters can be selected by cross-validation or by marginal likelihood. In Gaussian regression, integrating out the training function values gives \(\log p(y\mid X,\theta)=-\tfrac12r^\top K^{-1}r-\tfrac12\log\det K-\tfrac n2\log(2\pi)\), where \(r=y-m_0(X)\) and \(\theta\) collects the kernel and noise parameters. The quadratic term measures how plausible the residual vector is under the covariance, while the determinant normalizes the Gaussian density. Neither term is simply training prediction error. Maximizing this expression gives a point estimate of \(\theta\); it does not integrate uncertainty over \(\theta\).

Each fixed parameter choice can be scored using one factorization on all training rows, without a validation split. Optimizing that score usually evaluates many parameter choices and repeats the factorization; restarts add work. The objective can have local optima, and multiple starts can help but do not guarantee the global maximum. Predictive validation remains useful, especially when comparing assumptions or checking intervals. The returned posterior at an optimized parameter value is conditional on that value.

Dense inference and its costs

For standard dense exact regression, training factors an \(n\times n\) matrix, costing \(O(n^3)\) time and \(O(n^2)\) memory. Hyperparameter optimization repeats this work as parameters change. Prediction reuses the stored factor and solved coefficients: a mean prediction needs \(n\) kernel evaluations and a weighted sum; obtaining its variance also needs a triangular solve costing \(O(n^2)\). The factorization is not repeated for every query. A full covariance for many queries adds storage and computation beyond separate marginal standard deviations.

One dense float64 matrix occupies 800 MB at 10,000 rows and 80 GB at 100,000 rows, before factors and work arrays. Wall-clock time and feasibility depend on hardware, implementation, and the number of parameter evaluations. Many inducing-point methods use \(m\ll n\) auxiliary locations and reduce a batch training calculation to roughly \(O(nm^2+m^3)\); costs differ across algorithms. Approximation quality depends on the number and placement of those locations and the covariance structure.

GPs are useful candidates when observations are costly, structural assumptions are informative, and prediction uncertainty affects decisions. In Bayesian optimization, for example, many acquisition rules use the predictive mean and uncertainty to choose the next evaluation. The value of that uncertainty still depends on checking the model. A numerically exact Gaussian update can be too confident under a poorly specified prior or noise model.

References

• Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.