Gaussian Processes

Introduction

Throughout our study of statistics, we have worked with parametric models: models with a fixed number of parameters \(\boldsymbol{\theta}\) that does not grow with the training set size \(N\). For example, deep neural networks learn a function approximation \(f(\boldsymbol{x}; \boldsymbol{\theta})\) where \(\dim(\boldsymbol{\theta})\) is determined by the architecture. Such models can overfit when \(N\) is small or underfit when the model capacity is insufficient. We now consider a fundamentally different approach: nonparametric models, whose effective complexity adapts to the amount of available data.

The Gaussian Process (GP) is a fundamental example of a nonparametric Bayesian model. Rather than placing a prior distribution over a finite-dimensional parameter vector, we place a prior directly over functions \(p(f)\), and update it after observing data to obtain the posterior \(p(f \mid \mathcal{D})\). This function-space view builds on two key prerequisites: the multivariate normal distribution, whose conditioning identities we will use to derive the predictive posterior, and positive definite matrices, which ensure that our covariance structures are well-defined.

The key insight is to think of a function as an infinite-dimensional vector. Consider a Gaussian random vector \(\boldsymbol{f} = [f_1, f_2, \ldots, f_N]^\top\), characterized by its mean \(\boldsymbol{\mu} = \mathbb{E}[\boldsymbol{f}]\) and covariance \(\boldsymbol{\Sigma} = \operatorname{Cov}[\boldsymbol{f}]\).

Now, let \(f : \mathcal{X} \rightarrow \mathbb{R}\) be a function evaluated at a finite set of input points \[ \boldsymbol{X} = \{\boldsymbol{x}_n \in \mathcal{X}\}_{n=1}^N, \] and define \[ \boldsymbol{f}_X = [f(\boldsymbol{x}_1), f(\boldsymbol{x}_2), \ldots, f(\boldsymbol{x}_N)]^\top. \] If \(\boldsymbol{f}_X\) follows a joint multivariate Gaussian distribution for any finite set of input points, then \(f\) is said to follow a Gaussian process.

In practice, the mean function is often set to zero, \(m(\boldsymbol{x}) = 0\), since the data can be centered and the GP's expressiveness comes primarily from the covariance function.

For any finite collection of input points \(\boldsymbol{X} = \{\boldsymbol{x}_1, \boldsymbol{x}_2, \ldots, \boldsymbol{x}_N\}\), the corresponding function values follow a multivariate normal distribution: \[ p(\boldsymbol{f}_X \mid \boldsymbol{X}) = \mathcal{N}(\boldsymbol{f}_X \mid \boldsymbol{\mu}_X, \boldsymbol{K}_{X,X}), \] where \(\boldsymbol{\mu}_X = [m(\boldsymbol{x}_1), m(\boldsymbol{x}_2), \ldots, m(\boldsymbol{x}_N)]^\top\) and \((\boldsymbol{K}_{X,X})_{ij} = \mathcal{K}(\boldsymbol{x}_i, \boldsymbol{x}_j)\). When \(\boldsymbol{K}_{X,X}\) is strictly positive definite this is the multivariate normal distribution with the density above. When it is only positive semi-definite the distribution is degenerate, concentrated on a proper affine subspace, and has no density. The regression formulas below add observation noise, which makes the only matrix they invert invertible. The predictive distribution at the test inputs can still be degenerate in this sense, for example when two test inputs coincide.

Mercer Kernels

The kernel function defines the notion of similarity between input points, and hence the covariance between the function values there. In a Gaussian process, the choice of kernel is crucial since it determines the smoothness, periodicity, and other structural properties of the functions drawn from the GP prior.

While the kernel function intuitively measures similarity between inputs, not every symmetric function can serve as a covariance function in a Gaussian process. We now make precise the class of admissible kernels.

Note that this condition does not require each individual kernel value \(\mathcal{K}(\boldsymbol{x}_i, \boldsymbol{x}_j)\) to be nonnegative. Only the quadratic form built from the Gram matrix \(\boldsymbol{K}\) must be nonnegative.

The finite-sample condition above is all that a Gaussian process requires of its kernel. A deeper treatment belongs to the functional-analytic setting, where the same notion is studied under the name positive definite kernel. That treatment covers the connection to reproducing kernel Hilbert spaces and Mercer's integral-operator decomposition \(\mathcal{K}(\boldsymbol{x}, \boldsymbol{x}') = \sum_n \lambda_n e_n(\boldsymbol{x}) e_n(\boldsymbol{x}')\).

We often focus on stationary kernels, which depend only on the difference between inputs, \(\boldsymbol{r} = \boldsymbol{x} - \boldsymbol{x}'\), and not on their absolute locations. In many cases, only the Euclidean distance between inputs matters: \[ r = \|\boldsymbol{r}\| = \|\boldsymbol{x} - \boldsymbol{x}'\|. \]

One of the most widely used stationary kernels in practice is the radial basis function (RBF) kernel (also used in classification, where its random Fourier-feature approximation enables training that scales linearly in the number of samples): \[ \begin{align*} \mathcal{K}_\text{RBF}(r ; \ell) &= \exp \left( -\frac{r^2}{2\ell^2}\right) \\\\ &= \exp \left( -\frac{\| \boldsymbol{x} - \boldsymbol{x}'\|^2}{2\ell^2}\right) \end{align*} \] where \(\ell \gt 0\) is the length-scale parameter controlling how quickly the correlation between points decays with distance. This kernel is also known as the Gaussian kernel.

Furthermore, by replacing Euclidean distance with Mahalanobis distance, we obtain the generalized RBF kernel: \[ \mathcal{K}(\boldsymbol{r}; \boldsymbol{\Sigma}, \sigma^2) = \sigma^2 \exp \left( - \frac{1}{2} \boldsymbol{r}^\top \boldsymbol{\Sigma}^{-1} \boldsymbol{r}\right). \]

When \(\boldsymbol{\Sigma}\) is diagonal, that is, \(\boldsymbol{\Sigma} = \operatorname{diag}(\ell_1^2, \ldots, \ell_D^2)\), we obtain the Automatic Relevance Determination (ARD) kernel: \[ \mathcal{K}_\text{ARD}(\boldsymbol{r}; \boldsymbol{\Sigma}, \sigma^2) = \sigma^2 \exp \left( - \frac{1}{2} \sum_{d=1}^D \frac{1}{\ell_d^2 }r_d^2 \right) \] where \(\sigma^2\) is the overall variance and \(\ell_d \gt 0\) is the characteristic length scale of dimension \(d\), controlling the sensitivity of the function to that input. If a dimension \(d\) is irrelevant, letting \(\ell_d \to \infty\) removes the dependence on that input entirely.

While RBF and ARD kernels are infinitely smooth and often used as default choices, there are situations where we may want more flexible control over the smoothness of the functions. For example, in Bayesian optimization or when modeling rougher processes, it is useful to have kernels whose sample paths are only once or twice differentiable.

This motivates the Matérn kernels, a family of stationary kernels with an explicit smoothness parameter \(\nu \gt 0\). The family contains the RBF kernel as the limiting case \(\nu \to \infty\). By adjusting \(\nu\), we can interpolate between very rough kernels (like the exponential kernel) and infinitely smooth kernels (like the RBF).

For parameters \(\nu \gt 0\), \(\ell \gt 0\), the Matérn kernel is \[ \mathcal{K}_\text{Matérn}(r; \nu, \ell) = \frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}\,r}{\ell}\right)^\nu K_\nu\!\left(\frac{\sqrt{2\nu}\,r}{\ell}\right), \] where \(K_\nu\) denotes the modified Bessel function of the second kind, and the value at \(r = 0\) is \(1\), the limit of the right-hand side as \(r \to 0\).

The parameter \(\nu\) controls the smoothness of the resulting functions. Specifically, a zero-mean Gaussian process with this kernel is \(k\)-times mean-square differentiable if and only if \(\nu \gt k\). The sample paths themselves are almost surely \(\lceil \nu \rceil - 1\) times continuously differentiable, and not \(\lceil \nu \rceil\) times. This pathwise statement is stronger than mean-square differentiability and is stated here without proof. Common choices (up to the variance factor \(\sigma^2\)) admit closed forms:

• Exponential Kernel (Non-Differentiable):
\[ \mathcal{K}_\text{Matérn}(r; \frac{1}{2}, \ell ) = \exp \left(- \frac{r}{\ell}\right) \]
• Once Differentiable:
\[ \mathcal{K}_\text{Matérn}(r; \frac{3}{2}, \ell ) = \left(1 + \frac{\sqrt{3}r}{\ell}\right)\exp \left(- \frac{\sqrt{3}r}{\ell}\right) \]
• Twice Differentiable:
\[ \mathcal{K}_\text{Matérn}(r; \frac{5}{2}, \ell ) = \left(1 + \frac{\sqrt{5}r}{\ell} + \frac{5r^2}{3\ell^2}\right)\exp \left(- \frac{\sqrt{5}r}{\ell}\right) \]

In Gaussian process regression, the kernel can be multiplied by a positive amplitude parameter \(\sigma^2 \gt 0\) to scale the variance of the function values: \[ \mathcal{K}(r;\nu,\ell) \longrightarrow \sigma^2 \mathcal{K}(r;\nu,\ell). \] This scaling changes the magnitude of function fluctuations but does not affect the smoothness or correlation structure of the GP.

Some functions exhibit repeating or periodic behavior that is not captured by standard RBF or Matérn kernels. To model such patterns, we can use a periodic kernel, which explicitly encodes periodicity into the covariance function.

A common choice for a one-dimensional periodic kernel is: \[ \mathcal{K}_\text{per}(r; \ell, p) = \exp \Bigg( - \frac{2 \, \sin^2\big(\pi r / p\big)}{\ell^2} \Bigg), \] where \(p \gt 0\) is the period of the function and \(\ell \gt 0\) is the length-scale controlling how quickly correlations decay away from exact multiples of the period.

Periodic kernels are particularly useful in Gaussian process regression for modeling phenomena that repeat regularly over time or space, such as seasonal time series, cyclic patterns, or physical systems with inherent periodicity.

GP Regression

Suppose we observe a training data set \(\mathcal{D} = \{(\boldsymbol{x}_n, y_n)\}_{n=1}^N\), where each observation follows \[ y_n = f(\boldsymbol{x}_n) + \epsilon_n, \quad \epsilon_n \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma_y^2), \] and the noise \(\{\epsilon_n\}\), with \(\sigma_y^2 \gt 0\), is independent of the latent function \(f\). Under this model, the covariance of the noisy observations is \[ \begin{align*} \operatorname{Cov}[y_i, y_j] &= \operatorname{Cov}[f(\boldsymbol{x}_i), f(\boldsymbol{x}_j)] + \operatorname{Cov}[\epsilon_i, \epsilon_j] \\\\ &= \mathcal{K}(\boldsymbol{x}_i, \boldsymbol{x}_j) + \sigma_y^2 \, \mathbb{1}\{i = j\}, \end{align*} \] so that \[ \operatorname{Cov}[\boldsymbol{y} \mid \boldsymbol{X}] = \boldsymbol{K}_{X,X} + \sigma_y^2 \boldsymbol{I}_N. \]

Because a Gaussian process defines a joint Gaussian distribution over any finite collection of function values, we can express both the noisy training outputs and the noise-free test values together. For test inputs \(\boldsymbol{X}_* = \{\boldsymbol{x}_{*,m}\}_{m=1}^M\), let \(\boldsymbol{f}_* = [f(\boldsymbol{x}_{*,1}), \ldots, f(\boldsymbol{x}_{*,M})]^\top\) and \(\boldsymbol{\mu}_* = [m(\boldsymbol{x}_{*,1}), \ldots, m(\boldsymbol{x}_{*,M})]^\top\). Define the cross-covariance and test-test covariance blocks by \((\boldsymbol{K}_{X,*})_{ij} = \mathcal{K}(\boldsymbol{x}_i, \boldsymbol{x}_{*,j})\) and \((\boldsymbol{K}_{*,*})_{ij} = \mathcal{K}(\boldsymbol{x}_{*,i}, \boldsymbol{x}_{*,j})\). Then \[ \begin{bmatrix} \boldsymbol{y} \\ \boldsymbol{f}_* \end{bmatrix} \sim \mathcal{N} \Bigg( \begin{bmatrix} \boldsymbol{\mu}_X \\ \boldsymbol{\mu}_* \end{bmatrix}, \begin{bmatrix} \boldsymbol{K}_{X,X} + \sigma_y^2 \boldsymbol{I}_N & \boldsymbol{K}_{X,*} \\ \boldsymbol{K}_{X,*}^\top & \boldsymbol{K}_{*,*} \end{bmatrix} \Bigg). \]

The posterior predictive distribution at the test points \(\boldsymbol{X}_*\) is obtained by conditioning on the observed data: \[ p(\boldsymbol{f}_* \mid \mathcal{D}, \boldsymbol{X}_*) = \mathcal{N}\big( \boldsymbol{f}_* \mid \boldsymbol{\mu}_{* \mid X}, \boldsymbol{\Sigma}_{* \mid X} \big), \] where \[ \boldsymbol{\mu}_{* \mid X} = \boldsymbol{\mu}_* + \boldsymbol{K}_{X,*}^\top \big(\boldsymbol{K}_{X,X} + \sigma_y^2 \boldsymbol{I}_N \big)^{-1} (\boldsymbol{y} - \boldsymbol{\mu}_X), \] \[ \boldsymbol{\Sigma}_{* \mid X} = \boldsymbol{K}_{*,*} - \boldsymbol{K}_{X,*}^\top \big(\boldsymbol{K}_{X,X} + \sigma_y^2 \boldsymbol{I}_N \big)^{-1} \boldsymbol{K}_{X,*}. \]

These formulas are the conditional mean and covariance of the block \(\boldsymbol{f}_*\) of the multivariate normal distribution above, given the block \(\boldsymbol{y}\). We state this conditioning rule without proof. Here \(\boldsymbol{K}_{X,X} + \sigma_y^2 \boldsymbol{I}_N\) is invertible because \(\boldsymbol{K}_{X,X}\) is positive semi-definite and \(\sigma_y^2 \gt 0\). The posterior mean \(\boldsymbol{\mu}_{* \mid X}\) is \(\boldsymbol{\mu}_*\) plus a linear combination of the centered observations \(\boldsymbol{y} - \boldsymbol{\mu}_X\), with weights \(\boldsymbol{K}_{X,*}^\top (\boldsymbol{K}_{X,X} + \sigma_y^2 \boldsymbol{I}_N)^{-1}\) determined by the kernel. The posterior covariance \(\boldsymbol{\Sigma}_{* \mid X}\) starts from the prior covariance \(\boldsymbol{K}_{*,*}\) and subtracts the information gained from the observations. Since the subtracted term is positive semi-definite, conditioning never increases the predictive variance, and the more data we observe near a test point, the smaller the posterior uncertainty there tends to be.

The computational bottleneck is the inversion of the \(N \times N\) matrix \((\boldsymbol{K}_{X,X} + \sigma_y^2 \boldsymbol{I}_N)\), which costs \(O(N^3)\) in time and \(O(N^2)\) in storage. This cubic scaling limits standard GP regression to moderate dataset sizes (typically \(N \lesssim 10^4\)). For larger datasets, sparse approximations (inducing point methods), random Fourier features, and structured kernel interpolation provide scalable alternatives.