Gaussian Process Regression Formulation

Gaussian Processes serve as priors over functions. When applying this framework to regression, the core idea is straightforward: if an underlying function is assumed to originate from a Gaussian Process prior, and observed data points are noisy versions of this function, then predicting the function's value at new points becomes a standard Bayesian inference problem that involves multivariate Gaussians.

The Setup: Observations and the Latent Function

Recall that a Gaussian Process defines a distribution over functions. We usually specify it with a mean function $m(\mathbf{x})$ and a covariance function (kernel) $k(\mathbf{x}, \mathbf{x}')$. For simplicity, we often start by assuming a zero mean function, $m(\mathbf{x}) = 0$.

In a regression setting, we have a training dataset consisting of $N$ input points $X = \{\mathbf{x}_1, \dots, \mathbf{x}_N\}$ and their corresponding observed target values $\mathbf{y} = \{y_1, \dots, y_N\}^T$. The fundamental assumption in GP regression is that these observed values $y_i$ are noisy realizations of an underlying latent function $f$ evaluated at the inputs $\mathbf{x}_i$. Specifically, we assume:

$$y_i = f(\mathbf{x}_i) + \epsilon_i$$

where $f \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}'))$ and the noise term $\epsilon_i$ is typically assumed to be independent and identically distributed Gaussian noise with zero mean and variance $\sigma_n^2$:

$$\epsilon_i \sim \mathcal{N}(0, \sigma_n^2)$$

Let $\mathbf{f} = [f(\mathbf{x}_1), \dots, f(\mathbf{x}_N)]^T$ be the vector of latent function values at the training inputs. According to the GP prior, this vector follows a multivariate Gaussian distribution:

$$\mathbf{f} \sim \mathcal{N}(\mathbf{m}, K)$$

where $\mathbf{m}$ is the vector of prior means $[m(\mathbf{x}_1), \dots, m(\mathbf{x}_N)]^T$ and $K$ is the $N \times N$ covariance matrix (or Gram matrix) evaluated at all pairs of training inputs: $K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)$.

Combining the GP prior on $\mathbf{f}$ and the Gaussian noise model, the marginal distribution of the observed targets $\mathbf{y}$ also becomes Gaussian:

$$\mathbf{y} | X \sim \mathcal{N}(\mathbf{m}, K + \sigma_n^2 I)$$

Here, $I$ is the $N \times N$ identity matrix, and the noise variance $\sigma_n^2$ is added to the diagonal elements of the covariance matrix $K$, reflecting the independent noise added to each observation. This distribution $p(\mathbf{y} | X)$ is often called the marginal likelihood or evidence.

Predicting New Points

Our goal is to predict the function value $f_*$ at a new test input point $\mathbf{x}_*$. We might also want to predict at multiple test points $X_* = \{\mathbf{x}_{*1}, \dots, \mathbf{x}_{*M}\}$. Let $\mathbf{f}_* = [f(\mathbf{x}_{*1}), \dots, f(\mathbf{x}_{*M})]^T$ be the vector of latent function values at these test points.

The core property of GPs is that any finite set of function values drawn from a GP follows a multivariate Gaussian distribution. Therefore, the joint distribution of the observed training targets $\mathbf{y}$ and the latent function values $\mathbf{f}_*$ at the test points, conditioned on the inputs $X$ and $X_*$, is also Gaussian:

$$\begin{pmatrix} \mathbf{y} \\ \mathbf{f}_* \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} \mathbf{m} \\ \mathbf{m}_* \end{pmatrix}, \begin{pmatrix} K(X, X) + \sigma_n^2 I & K(X, X_*) \\ K(X_*, X) & K(X_*, X_*) \end{pmatrix} \right)$$

Let's break down the covariance matrix components:

• $K(X, X)$: The $N \times N$ matrix of covariances between all pairs of training inputs (we previously called this $K$).
• $K(X, X_*)$: The $N \times M$ matrix of covariances between training inputs and test inputs. $K(X, X_*)_{ij} = k(\mathbf{x}_i, \mathbf{x}_{*j})$.
• $K(X_*, X)$: The $M \times N$ matrix of covariances between test inputs and training inputs. Note that $K(X_*, X) = K(X, X_*)^T$.
• $K(X_*, X_*)$: The $M \times M$ matrix of covariances between all pairs of test inputs. $K(X_*, X_*)_{ij} = k(\mathbf{x}_{*i}, \mathbf{x}_{*j})$.

And the mean vectors:

• $\mathbf{m}$: The $N \times 1$ vector of prior means at training inputs.
• $\mathbf{m}_*$: The $M \times 1$ vector of prior means at test inputs.

The Predictive Distribution

We are interested in the posterior predictive distribution $p(\mathbf{f}_* | X, \mathbf{y}, X_*)$, which tells us the distribution of function values at the test points given the observed training data. Since the joint distribution is Gaussian, this conditional distribution is also Gaussian.

Using standard results for conditional Gaussian distributions, we can derive the mean and covariance of this predictive distribution. Assuming a zero prior mean function ($m(\mathbf{x}) = 0$, thus $\mathbf{m}=\mathbf{0}$ and $\mathbf{m}_*=\mathbf{0}$) for simplicity first:

The posterior predictive mean is:

$$\mathbb{E}[\mathbf{f}_* | X, \mathbf{y}, X_*] = K(X_*, X) [K(X, X) + \sigma_n^2 I]^{-1} \mathbf{y}$$

The posterior predictive covariance is:

$$\text{Cov}[\mathbf{f}_* | X, \mathbf{y}, X_*] = K(X_*, X_*) - K(X_*, X) [K(X, X) + \sigma_n^2 I]^{-1} K(X, X_*)$$

If we use a non-zero prior mean function $m(\mathbf{x})$, the equations become slightly more general:

Posterior Predictive Mean:

$$\mathbb{E}[\mathbf{f}_* | X, \mathbf{y}, X_*] = \mathbf{m}_* + K(X_*, X) [K(X, X) + \sigma_n^2 I]^{-1} (\mathbf{y} - \mathbf{m})$$

Posterior Predictive Covariance:

$$\text{Cov}[\mathbf{f}_* | X, \mathbf{y}, X_*] = K(X_*, X_*) - K(X_*, X) [K(X, X) + \sigma_n^2 I]^{-1} K(X, X_*)$$

Notice that the posterior covariance does not depend on the observed targets $\mathbf{y}$, only on the inputs $X$, $X_*$, the kernel $k$, and the noise variance $\sigma_n^2$.

Interpretation of the Predictive Equations

• Predictive Mean: The predictive mean $\mathbb{E}[\mathbf{f}_* | \dots]$ starts with the prior mean $\mathbf{m}_*$ and adds a correction term. This correction term is a linear combination of the centered observed targets $(\mathbf{y} - \mathbf{m})$. The weights of this combination depend on the covariances between the test points $X_*$ and the training points $X$, modulated by the inverse of the noisy covariance matrix of the training data. Intuitively, training points that are "closer" (as defined by the kernel $k$) to a test point $\mathbf{x}_*$ will have a stronger influence on its predicted value.
• Predictive Covariance: The predictive covariance $\text{Cov}[\mathbf{f}_* | \dots]$ starts with the prior covariance $K(X_*, X_*)$ and subtracts a positive semi-definite matrix. This subtraction represents the reduction in uncertainty about $\mathbf{f}_*$ gained from observing the training data $\mathbf{y}$. The reduction is larger when the training points provide more information about the test points (i.e., when $K(X_*, X)$ is large). The diagonal elements of this covariance matrix, $\text{Var}[f_{*i} | \dots]$, give us the predictive variance for each individual test point $\mathbf{x}_{*i}$. This variance is typically smaller near the training data points and larger further away, reflecting increased uncertainty where we lack data.

The ability to obtain not just point predictions (the mean) but also a full measure of uncertainty (the variance/covariance) is a significant advantage of GP regression.

Example of 1D GP regression. The blue points are training data. The red line is the posterior predictive mean. The shaded red area represents +/- two standard deviations from the mean, visualizing the model's uncertainty. Uncertainty is lower near the training points and increases further away.

Computational Challenges

The primary computational bottleneck in standard GP regression is the calculation of the inverse term $[K(X, X) + \sigma_n^2 I]^{-1}$. This involves inverting an $N \times N$ matrix, which typically requires $O(N^3)$ computational time and $O(N^2)$ memory to store the matrix $K$.

• Training/Inference: Computing the inverse takes $O(N^3)$.
• Prediction: Once the inverse (or its Cholesky decomposition) is computed and stored, predicting the mean for $M$ test points takes $O(NM)$ time, and predicting the full covariance takes $O(NM^2)$. Predicting the variance only for $M$ test points takes $O(NM)$ time after the initial $O(N^3)$ investment.

This $O(N^3)$ complexity limits the applicability of exact GP regression to datasets with typically fewer than a few thousand data points. We will discuss methods to scale GPs to larger datasets later in this chapter.

Now that we have the core formulation for regression, the next steps involve understanding how to choose the kernel function $k(\mathbf{x}, \mathbf{x}')$ and how to set its parameters (hyperparameters) along with the noise variance $\sigma_n^2$.