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.994337The 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.0000For 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.
Exercises
1. Extrapolation with fixed and fitted periods. Fit RBF, Matérn, and periodic models to two periods of a noisy sine on [0,6]. Evaluate 24 points beyond 6 through 12, covering two further periods. Compare a period fixed at 3 with one optimized from an initial value of 3.
Report error, mean latent standard deviation, and the period actually used. What uncertainty is absent when a fitted period is treated as known?
Solution
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import (RBF, Matern, ExpSineSquared,
ConstantKernel as C)
rng = np.random.default_rng(0)
X = np.linspace(0, 6, 40).reshape(-1, 1)
y = np.sin(2 * np.pi * X.ravel() / 3) + 0.05 * rng.normal(size=40)
X_q = np.linspace(6.25, 12, 24).reshape(-1, 1) # 24 query points beyond the training interval
truth = np.sin(2 * np.pi * X_q.ravel() / 3)
for name, kern in (("RBF", C(1.0) * RBF(1.0)),
("Matern 1.5", C(1.0) * Matern(1.0, nu=1.5)),
("period fixed", C(1.0) * ExpSineSquared(1.0, 3.0, periodicity_bounds="fixed")),
("period fit", C(1.0) * ExpSineSquared(1.0, 3.0, periodicity_bounds=(2.5, 3.5)))):
gp = GaussianProcessRegressor(kernel=kern, alpha=0.05**2).fit(X, y)
mu, sd = gp.predict(X_q, return_std=True)
print(f"{name:12s} extrapolation MAE {np.abs(mu - truth).mean():.4f}"
f" mean sd {sd.mean():.4f}")
if name.startswith("period"):
print(f" period used {gp.kernel_.k2.periodicity:.6f}")
# RBF extrapolation MAE 0.8809 mean sd 1.3840
# Matern 1.5 extrapolation MAE 0.6501 mean sd 1.0653
# period fixed extrapolation MAE 0.0176 mean sd 0.0158
# period used 3.000000
# period fit extrapolation MAE 0.0319 mean sd 0.0158
# period used 2.994646The fixed-period model has MAE 0.0176 and mean latent standard deviation 0.0158 on this query grid; RBF gives 0.8809 and 1.3840. Noise variance is set to the generating value \(0.05^2=0.0025\) for every model. These compare error against the noiseless sine and uncertainty in the latent function, not in a future noisy measurement. Low average error and a reported standard deviation do not establish interval coverage from one run.
periodicity_bounds="fixed" keeps the first periodic model at period 3. Its prior identifies \(f(x+3)\) with \(f(x)\), so distant matching phases remain connected to observations. The fitted-period model instead searches within [2.5,3.5] and reaches about 2.994646. Its MAE is 0.0319, while its mean latent standard deviation still rounds to 0.0158. The two standard deviations omit uncertainty about the other optimized kernel parameters as well. In particular, the second one treats its estimated period as if it were known.
For this fixed-period prior, changing the true cycle or allowing it to drift would violate an assumption and could cause extrapolation errors with narrow bands. Their size would depend on the mismatch and forecast horizon. The RBF and Matérn fits instead lose correlation at long distances and approach their specified prior means and fitted prior variances. Extrapolation depends on the mean and noise model as well as the kernel.
To include uncertainty about the period, use a distribution over plausible periods and integrate or approximate the resulting predictive mixture. Merely optimizing the period does not do this. For a scalar query, the total variance includes both the average conditional variance and the variance of conditional means across parameter values. That distinction matters when future phase depends strongly on a small uncertainty in the period.
2. Reading the marginal likelihood. Evaluate five fixed length scales with signal variance 1 and noise variance 0.04. Then optimize only the length scale, followed by an optimization of all three parameters. Compare the likelihood terms and the search results.
The fixed sweep and the joint search explore different sets of models. Which comparison isolates the effect of length scale?
Solution
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel, ConstantKernel as C
rng = np.random.default_rng(1)
X = rng.uniform(-4, 4, (60, 1))
y = np.sin(X.ravel() * 1.5) + 0.2 * rng.normal(size=60)
for ell in (0.05, 0.2, 0.7, 2.0, 10.0):
kern = C(1.0, "fixed") * RBF(ell, "fixed") + WhiteKernel(0.04, "fixed")
gp = GaussianProcessRegressor(kernel=kern, optimizer=None).fit(X, y)
print(f"length scale {ell:6.2f}"
f" lml {gp.log_marginal_likelihood_value_:9.3f}"
f" quadratic {-0.5 * y @ gp.alpha_:9.3f}"
f" logdet {-np.log(np.diag(gp.L_)).sum():9.3f}")
only_ell = GaussianProcessRegressor(
kernel=C(1.0, "fixed") * RBF(1.0) + WhiteKernel(0.04, "fixed"),
n_restarts_optimizer=3, random_state=0).fit(X, y)
print(f"length only: {only_ell.kernel_}")
print(f" lml {only_ell.log_marginal_likelihood_value_:.3f}")
gp = GaussianProcessRegressor(kernel=C(1.0) * RBF(1.0) + WhiteKernel(0.1),
n_restarts_optimizer=3, random_state=0).fit(X, y)
print(f"optimized: {gp.kernel_}")
print(f" lml {gp.log_marginal_likelihood_value_:.3f}")
# length scale 0.05 lml -55.159 quadratic -11.857 logdet 11.834
# length scale 0.20 lml -24.978 quadratic -11.796 logdet 41.954
# length scale 0.70 lml 0.696 quadratic -15.348 logdet 71.181
# length scale 2.00 lml -26.562 quadratic -54.925 logdet 83.499
# length scale 10.00 lml -246.261 quadratic -281.475 logdet 90.350
# length only: 1**2 * RBF(length_scale=1.17) + WhiteKernel(noise_level=0.04)
# lml 4.733
# optimized: 0.993**2 * RBF(length_scale=1.15) + WhiteKernel(noise_level=0.0218)
# lml 8.523Of the five fixed candidates, length scale 0.7 has the highest log marginal likelihood, 0.696. Relative to it, 0.05 is about 55.9 nats lower and 10 is about 247.0 nats lower; a nat uses the natural logarithm. Optimizing only length scale finds about 1.17 with score 4.733. Allowing signal and noise variances to change gives length scale 1.15 and score 8.523. That latter score is not a point on the fixed-amplitude, fixed-noise curve, and five candidates alone do not locate its continuous maximum.
The printed quadratic and logdet columns are the first two terms of the formula; adding the constant \(-n\log(2\pi)/2\) recovers the score. At length scale 0.05, they are −11.857 and +11.834. At 0.7, they are −15.348 and +71.181: the quadratic term becomes less favorable, but the determinant contribution improves more. At length scale 10, the quadratic term falls to −281.475, outweighing a determinant contribution of +90.350. This explains these scores without equating the quadratic term to training fit or assuming the determinant term must be negative.
The joint fit estimates noise variance 0.0218, below the generating value 0.04. This one sample does not establish a general downward bias or identify whether fitted smoothness absorbed noise. Compare fitted signal scale, length scale, and noise together, and check predictive performance on suitable held-out data. A small fitted noise level does not by itself prove interpolation: additional diagonal noise, repeated conflicting observations, and numerical regularization can also matter. Here the GPR default alpha=1e-10 is a small extra diagonal term alongside WhiteKernel.
The marginal likelihood uses all sixty rows without a held-out split because it is the model’s density for the observed vector after integrating out latent function values. It is a selection criterion, not an independent estimate of prediction error. Each optimizer call may make many evaluations; n_restarts_optimizer=3 means the initial run plus three additional starts. The five fixed scores do not demonstrate a separate local optimum at very short length scales.
3. Measure dense fitting costs. Time repeated GP fits at 200, 400, 800, and 1600 points with fixed hyperparameters. Separately calculate matrix storage, including sizes too large to allocate in this experiment.
Distinguish measured timings from calculated byte counts. What additional memory does a fitted model retain?
Solution
import numpy as np, time
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
from threadpoolctl import threadpool_limits
rng = np.random.default_rng(2)
with threadpool_limits(limits=1):
for n in (200, 400, 800, 1600):
X = rng.uniform(-3, 3, (n, 1))
y = np.sin(X.ravel()) + 0.1 * rng.normal(size=n)
model = GaussianProcessRegressor(
kernel=C(1.0, "fixed") * RBF(1.0, "fixed"),
alpha=1e-2, optimizer=None)
model.fit(X, y)
t = time.perf_counter()
for _ in range(3):
model.fit(X, y)
ms = (time.perf_counter() - t) * 1000 / 3
print(f"n={n:5d} mean fit {ms:8.1f} ms")
for n in (200, 400, 800, 1600):
print(f"n={n:5d} kernel matrix {n * n * 8 / 2 ** 20:8.2f} MiB")
for n in (10_000, 100_000, 1_000_000):
print(f"n={n:9,d} kernel matrix {n * n * 8 / 2 ** 30:12,.1f} GiB")
# n= 200 mean fit 1.7 ms # varies by machine
# n= 400 mean fit 3.9 ms # varies by machine
# n= 800 mean fit 19.4 ms # varies by machine
# n= 1600 mean fit 127.3 ms # varies by machine
# n= 200 kernel matrix 0.31 MiB
# n= 400 kernel matrix 1.22 MiB
# n= 800 kernel matrix 4.88 MiB
# n= 1600 kernel matrix 19.53 MiB
# n= 10,000 kernel matrix 0.7 GiB
# n= 100,000 kernel matrix 74.5 GiB
# n=1,000,000 kernel matrix 7,450.6 GiBThe benchmark limits native thread pools to one thread, warms up each dataset size, and averages three subsequent fits. It times the full fit, including matrix construction, factorization, and solves, rather than an isolated Cholesky operation. In the recorded run, 400 to 1600 rows takes 3.9 to 127.3 ms, roughly 33 times longer. Pure cubic scaling for that fourfold increase would give 64 times; these measurements do not establish a cubic timing law by themselves. Constants, other operations, caching, and numerical-library behavior all contribute.
The storage lines are calculated from \(8n^2\) bytes for one dense float64 matrix and rounded for display; the large matrices are not allocated. At 100,000 rows, one matrix is about 74.5 GiB. Dense GPR also retains the Cholesky factor and uses intermediate arrays, so peak memory exceeds that single-matrix count. Whether memory or runtime becomes limiting first depends on the machine and workload.
Inducing-point approximations summarize covariance relationships through auxiliary locations, which can be optimized and need not be selected training rows. For methods with a batch cost near \(O(nm^2+m^3)\), the leading dependence on \(n\) is linear only while \(m\) is held fixed. Increasing the dataset or the complexity of the function may require more inducing locations to maintain accuracy. Minibatch variational methods and structured kernels offer other computational tradeoffs.
Use measured fitting and prediction costs to decide whether dense inference is practical. If an approximation is needed, compare both predictions and uncertainty with a denser or exact reference where feasible.
References
- Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.
Discover more from Insightful Data Lab
Subscribe to get the latest posts sent to your email.