Gaussian Process Regression

Gaussian Process Regression

GP regression estimates an unknown function $f:\Omega \to \mathbb{R}$ from noisy observations $y_i=f(x_i)+{ϵ}_i$ with ${ϵ}_i\stackrel{\text{iid}}{\sim }\mathcal{N}(0,{\sigma }^2)$, by placing a Gaussian process prior $f\sim \text{GP}(0,K)$ on the function itself. Because the joint distribution of observations and test values is multivariate Gaussian, the posterior is a GP with closed-form mean and kernel.

Given test point $x$ and Gram matrix $\mathbf{K}$ with ${\mathbf{K}}_{ij}=K(x_i,x_j)$ and ${\mathbf{k}}_i=K(x,x_i)$:

$$f(x)\mid \mathbf{y}\sim \mathcal{N}\left(\underset{\text{posterior mean}}{\underbrace{{\mathbf{k}}^{\top }(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{y}}},\underset{\text{posterior variance}}{\underbrace{K(x,x)-{\mathbf{k}}^{\top }(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{k}}}\right).$$

Setting ${\sigma }^2=0$ gives the interpolation case; with noise, observations are allowed to disagree with the estimate and repeated measurements at the same $x_i$ become informative. The marginal likelihood
$p(\mathbf{y}\mid \sigma ,\mathbf{\theta })=\mathcal{N}(\mathbf{y};0,{\mathbf{K}}_{\mathbf{\theta }}+{\sigma }^2\mathbf{I})$
is used to learn kernel hyperparameters by empirical Bayes or full posterior sampling.

Computational cost: the naive posterior computation is $O(n^3)$ from inverting the $n\times n$ matrix $\mathbf{K}+{\sigma }^2\mathbf{I}$. Scalability beyond $n\approx 10^4$ requires inducing-point or low-rank approximations (not in scope here). When the kernel has a linear component with improper $\mathcal{N}(0,C)$ priors, the Sherman-Morrison-Woodbury identity makes the $C\to \infty$ limit well-defined and computable.

Prerequisites

• 3.10 Gaussian Processes — parent index.
• 3.10.1 Kernels — kernel design and algebra.
• 6.2 Multivariate Gaussian Distribution — conditioning formula for jointly Gaussian vectors.
• 3.7 Marginal Likelihood — used for hyperparameter estimation.

Intuition

GP regression is the nonparametric analog of 3.8.1 Bayesian Linear Regression. Instead of putting a prior on a finite-dimensional coefficient vector $\mathbf{\beta }$, put a prior directly on the function $f$. The rest follows: the posterior of $f$ given noisy observations is a new GP, the posterior mean is a kernel-weighted smoother, and the posterior variance grows smoothly between observation points — small near data, large far from data.

Two kinds of uncertainty get combined. The noise variance ${\sigma }^2$ captures how much the $y_i$ deviate from $f(x_i)$. The kernel variance $K(x,x)$ captures how much prior uncertainty there is about $f(x)$ itself. The regularized-kernel matrix $\mathbf{K}+{\sigma }^2\mathbf{I}$ cleanly combines both: the ${\sigma }^2\mathbf{I}$ term is the exact same diagonal loading that makes the linear-regression MAP into ridge regression, and here it plays the analogous role of keeping the matrix invertible and damping overfitting.

Unlike interpolation, observations are averaged rather than exactly matched. At an observed input $x_i$, the posterior mean of $f(x_i)$ need not equal $y_i$ — it typically pulls $y_i$ toward the GP prior mean (zero) and toward neighboring observations. Repeated noisy measurements at the same $x_i$ are informative: each reduces posterior variance at $x_i$ as more information accumulates.

Properties

(Posterior is a GP) Closed-form posterior mean and kernel.

The posterior process $f\mid \mathbf{y}$ is itself a Gaussian process with mean function and covariance kernel

$${\mu }_{\text{post}}(x)=\mathbf{k}(x{)}^{\top }(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{y},$$ $$K_{\text{post}}(x,x^{\prime })=K(x,x^{\prime })-\mathbf{k}(x{)}^{\top }(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{k}(x^{\prime }),$$

where $\mathbf{k}(x{)}_i=K(x,x_i)$. Because the posterior is itself a GP, any functional of $f$ (derivative, integral, linear predictor) is also Gaussian with computable mean and variance. This closure under conditioning is what makes GPs so clean to work with analytically.

(Posterior Mean is a Kernel Smoother) Linear combination of observed $y$ values.

The posterior mean ${\mu }_{\text{post}}(x)=\mathbf{k}(x{)}^{\top }(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{y}$ is a linear combination of the observations with weights depending only on inputs and kernel, not on $\mathbf{y}$. Writing $\mathbf{\alpha }=(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{y}$, the predictor is ${\mu }_{\text{post}}(x)={\sum }_i{\alpha }_{i}K(x,x_i)$ — a weighted sum of kernel functions centered at the training inputs. This is formally identical to the predictor in kernel ridge regression, which reveals that GP regression and kernel ridge regression are the same algorithm with different probabilistic interpretations.

(Posterior Variance Ignores Observations) Depends only on inputs and kernel.

The posterior variance $K(x,x)-{\mathbf{k}}^{\top }(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{k}$ depends on the inputs $x_1,\dots ,x_n$ and the test point $x$, but not on the observed $y$ values. This is because the conditional covariance of a Gaussian depends only on the covariance matrix, not on the observed values. A counterintuitive consequence: the posterior variance can be computed before any data is collected. This is the basis for Bayesian optimization and active learning strategies, which choose new query points $x_{\text{new}}$ to minimize posterior variance before looking at $y_{\text{new}}$.

(Hyperparameter Estimation via Marginal Likelihood) Optimize evidence, not held-out loss.

The kernel depends on hyperparameters $\mathbf{\theta }$ (length-scale, signal variance, etc.) and the noise level ${\sigma }^2$. These are estimated by maximizing (or placing priors on) the log-marginal likelihood:

$$\log p(\mathbf{y}\mid \mathbf{\theta },\sigma )=-\frac{1}{2}{\mathbf{y}}^{\top }({\mathbf{K}}_{\mathbf{\theta }}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{y}-\frac{1}{2}\log |{\mathbf{K}}_{\mathbf{\theta }}+{\sigma }^2\mathbf{I}|-\frac{n}{2}\log 2\pi .$$

The first term rewards fit, the second penalizes complexity — the log-determinant is large when the kernel allows high-variance functions, providing automatic regularization. Unlike cross-validation, marginal likelihood does not require held-out data and uses all $n$ observations. Full details and computational tricks in 3.7 Marginal Likelihood.

(Sherman-Morrison-Woodbury for Infinite-Variance Priors) Linear components with $C\to \infty$ priors are handled analytically.

When the kernel has a linear component $K(s,t)=C(1+st)+{\tau }^2{K}_I(s,t)$ with $C\to \infty$ (flat priors on intercept and slope), direct computation is unstable. Using the Woodbury identity

$$(C\mathbf{X}{\mathbf{X}}^{\top }+\mathbf{\Sigma }{)}^{-1}={\mathbf{\Sigma }}^{-1}-{\mathbf{\Sigma }}^{-1}\mathbf{X}(C^{-1}{\mathbf{I}}_d+{\mathbf{X}}^{\top }{\mathbf{\Sigma }}^{-1}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }{\mathbf{\Sigma }}^{-1}$$

and sending $C\to \infty$, the posterior mean for the combined kernel reduces to

$$\hat{f}(x)=(1,x)({\mathbf{X}}^{\top }{\mathbf{A}}_{\gamma }^{-1}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }{\mathbf{A}}_{\gamma }^{-1}\mathbf{y}+{\gamma }^2{\mathbf{k}}_I^{\top }({\mathbf{A}}_{\gamma }^{-1}-{\mathbf{A}}_{\gamma }^{-1}\mathbf{X}({\mathbf{X}}^{\top }{\mathbf{A}}_{\gamma }^{-1}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }{\mathbf{A}}_{\gamma }^{-1})\mathbf{y},$$

where ${\mathbf{A}}_{\gamma }=\mathbf{I}+{\gamma }^2{\mathbf{K}}_I$ and $\gamma =\tau /\sigma$. This expression depends only on the identifiable ratio $\gamma$, not on $\tau$ and $\sigma$ separately. The first term is the ordinary-least-squares fit projected onto the linear component; the second is the IBM-smoothed residual.

(Computational Cost is $O(n^3)$) Dominated by the matrix solve.

Computing $(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{y}$ requires a Cholesky factorization of the $n\times n$ matrix, which costs $O(n^3)$ time and $O(n^2)$ memory. This makes naive GP regression impractical beyond $n\approx 10^4$. Standard scaling strategies (not in scope here): inducing-point approximations ($O(n{m}^2)$ with $m$ inducing points), structured kernel interpolation, random-feature approximations, or Kronecker methods when inputs lie on a grid. For the problem sizes in STAT 238 (a few hundred observations), the $O(n^3)$ cost is inconsequential.

(Separable-Noise Predictive Distribution) Predictions account for both signal and noise variance.

To predict a new noisy observation $y_{\ast }=f(x_{\ast })+{ϵ}_{\ast }$ rather than the underlying $f(x_{\ast })$, add ${\sigma }^2$ to the posterior variance:

$$y_{\ast }\mid \mathbf{y}\sim \mathcal{N}\left({\mathbf{k}}^{\top }(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{y},K(x_{\ast },x_{\ast })-{\mathbf{k}}^{\top }(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{k}+{\sigma }^2\right).$$

The noise variance cleanly separates: predictions of the underlying function are narrower than predictions of new observations by exactly ${\sigma }^2$. This is the GP analog of the prediction-interval / confidence-interval distinction in linear regression.

Derivation

Start from the joint distribution of the noisy observations and the latent value $f(x)$ at a test point. Under the prior, $f$ is a GP with mean zero and kernel $K$, so

$$\left(\begin{array}{c}\mathbf{y}\\ f(x)\end{array}\right)\sim \mathcal{N}\left(0,\left(\begin{array}{cc}\mathbf{K}+{\sigma }^2\mathbf{I}& \mathbf{k}\\ {\mathbf{k}}^{\top }& K(x,x)\end{array}\right)\right),$$

using $\text{Cov}(y_i,y_j)=K(x_i,x_j)+{\sigma }^2{\delta }_{ij}$ and $\text{Cov}(y_i,f(x))=K(x_i,x)$. Applying the standard Gaussian conditioning formula

$$\left(\begin{array}{c}\mathbf{a}\\ b\end{array}\right)\sim \mathcal{N}\left(0,\left(\begin{array}{cc}\mathbf{A}& \mathbf{c}\\ {\mathbf{c}}^{\top }& d\end{array}\right)\right)⟹b\mid \mathbf{a}\sim \mathcal{N}({\mathbf{c}}^{\top }{\mathbf{A}}^{-1}\mathbf{a},d-{\mathbf{c}}^{\top }{\mathbf{A}}^{-1}\mathbf{c})$$

with $\mathbf{a}=\mathbf{y}$, $b=f(x)$, $\mathbf{A}=\mathbf{K}+{\sigma }^2\mathbf{I}$, $\mathbf{c}=\mathbf{k}$, $d=K(x,x)$ gives the headline result. The derivation extends immediately to joint prediction at multiple test points by replacing the scalar $b$ with a test vector and adjusting matrix dimensions.

Examples

(Easy) Predict a GP with Brownian motion kernel between two observations.

Observe $y_1=1.2$ at $x_1=0.2$ and $y_2=2.0$ at $x_2=0.8$, both with noise $\sigma =0$. Use the scaled BM plus constant kernel $K(s,t)=C+{\tau }^2\min (s,t)$ with $C\to \infty$ and $\tau =1$. Predict $f(0.5)$.

Solution. With $C\to \infty$, the posterior mean becomes the linear interpolant (Example 2.1 in L22):

$$\hat{f}(0.5)=\frac{0.8-0.5}{0.8-0.2}(1.2)+\frac{0.5-0.2}{0.8-0.2}(2.0)=\frac{1}{2}(1.2)+\frac{1}{2}(2.0)=\mathrm{1.6.}$$

Outside $[x_1,x_2]$, the posterior mean is constant-extrapolated: $\hat{f}(0.1)=1.2$, $\hat{f}(0.9)=2.0$. This confirms that the Brownian motion kernel induces linear interpolation inside the observed range and constant extrapolation outside.

(Medium) Hyperparameter fitting via marginal-likelihood grid.

Given $n$ observations and the scaled-IBM-plus-linear-term kernel with ratio $\gamma$, sketch how to learn $\gamma$ and $\sigma$ from the data.

Solution. Place a prior $p(\gamma )\propto 1/\gamma$ on $[{\gamma }_{\text{low}},{\gamma }_{\text{high}}]$ (log-uniform). After integrating out $\sigma$ analytically (inverse-Gamma conjugacy) and letting $C\to \infty$ via Sherman-Morrison-Woodbury, the marginal posterior of $\gamma$ is

$$p(\gamma \mid \mathbf{y})\propto \frac{p(\gamma )|{\mathbf{A}}_{\gamma }{|}^{-1/2}|{\mathbf{X}}^{\top }{\mathbf{A}}_{\gamma }^{-1}\mathbf{X}{|}^{-1/2}}{{\left({\mathbf{y}}^{\top }[{\mathbf{A}}_{\gamma }^{-1}-{\mathbf{A}}_{\gamma }^{-1}\mathbf{X}({\mathbf{X}}^{\top }{\mathbf{A}}_{\gamma }^{-1}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }{\mathbf{A}}_{\gamma }^{-1}]\mathbf{y}\right)}^{(n/2)-1}}$$

with ${\mathbf{A}}_{\gamma }=\mathbf{I}+{\gamma }^2{\mathbf{K}}_I$. Workflow:

1. Choose a grid ${\gamma }_1,\dots ,{\gamma }_G$ (log-spaced between underfit and overfit extremes).
2. Evaluate $\log p({\gamma }_k\mid \mathbf{y})$ on the grid.
3. Sample $\gamma$ from the discrete grid posterior.
4. For each $\gamma$-sample, draw $\sigma$ from the inverse-Gamma conditional.
5. For each $(\gamma ,\sigma )$ pair, draw $\mathbf{\beta }$ and the posterior function values via the closed-form Gaussian formulas.
   This produces posterior samples of the function — full Bayesian inference without MCMC.

(Hard) Two noisy measurements at the same input reduce posterior variance.

Observe $y_1=0.9$ and $y_2=1.1$ both at $x=0.5$, with noise $\sigma =0.3$. Use an RBF kernel with ${\sigma }_f^2=1$ and any length-scale $\ell >0$. Compute the posterior mean and variance of $f(0.5)$.

Solution. Both observations are at the same input, so the Gram matrix is

$$\mathbf{K}+{\sigma }^2\mathbf{I}=\left(\begin{array}{cc}{\sigma }_f^2+{\sigma }^2& {\sigma }_f^2\\ {\sigma }_f^2& {\sigma }_f^2+{\sigma }^2\end{array}\right)=\left(\begin{array}{cc}1.09& 1\\ 1& 1.09\end{array}\right).$$

Inverting gives $(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}=\frac{1}{1.09^2-1}\left(\begin{array}{cc}1.09& -1\\ -1& 1.09\end{array}\right)\approx \frac{1}{0.1881}\left(\begin{array}{cc}1.09& -1\\ -1& 1.09\end{array}\right)$. The vector $\mathbf{k}=(1,1{)}^{\top }$ (both at $x=0.5$). Posterior mean:

$${\mu }_{\text{post}}={\mathbf{k}}^{\top }(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{y}=\frac{0.09}{0.1881}(y_1+y_2)\approx 0.4785\times 2.0=0.957,$$

close to the mean of the two observations (1.0) but shrunk slightly toward the prior mean (0). Posterior variance:

$$K(0.5,0.5)-{\mathbf{k}}^{\top }(\mathbf{K}+{\sigma }^2\mathbf{I}{)}^{-1}\mathbf{k}=1-\frac{2\times 0.09}{0.1881}\approx 1-0.957=0.043,$$

smaller than the posterior variance from a single observation ${\sigma }_f^2{\sigma }^2/({\sigma }_f^2+{\sigma }^2)\approx 0.0826$. Replication at the same input did tighten the posterior — as expected.

See Also

• 3.10 Gaussian Processes — parent index.
• 3.10.1 Kernels — kernel design and algebra.
• 3.10.3 Integrated Brownian Motion — the IBM kernel used in the STAT 238 running example.
• 3.8.1 Bayesian Linear Regression — parametric cousin; GP regression is the nonparametric limit.
• 3.7 Marginal Likelihood — engine for hyperparameter estimation.
• 6.2 Multivariate Gaussian Distribution — conditioning formula.