A Second Introduction to Gaussian Processes for Function Approximation

One does not simply get GPs in one sitting. This is particularly true if you are an engineer or computer scientist and haven’t taken a course on stochastic processes or anything of that nature.

My goal with this second 0: As a first introduction, this may make for slightly tedious reading: too thin on the motivation and thick in the details. For more easy-going introductions, see this talk from David MacKay, this one from Rich Turner or the introductory section of the GPML book. ^0[0] introduction to GPs is to provide a self-contained treatment of GPs for the reader who’s sat through a first course or lecture on the subject, yet still feels like they’re missing the full picture.

There are a few essential preliminaries to be familiar with, without which the going can be rather tough. 1: I made this mistake when first learning about GPs: I was able to apply the formulae in a rote manner and code up the model. Yet I didn’t feel I understood why the equations worked. ^1[1] These include several Gaussian identities, the definition of a stochastic process, among a few others. Happily, these preliminaries aren’t difficult to grasp, and we will cover them as they come up 2: If you’d like to cover them preemptively, then begin by taking a read through my previous article on the Gaussian distribution, followed by the this article’s Appendix covering the measure theory. ^2[2] . With a good hold on these key ideas, the rest of the GP framework follows almost effortlessly.

Gaussian Processes

The terminology of measure theory proves to be rather valuable in succinctly expressing the key ideas behind GPs. 3: Often, introductions to GPs will avoid this for the sake of accessibility; choosing analogies to parametric function approximation and other visualisations instead. While intuitive and accessible, they might lack the precision that can be afforded with just a little bit of basic measure theory. ^3[3]

If the words probability space or Borel sigma algebra don’t mean very much to you, then take a look through the Appendix before proceeding, in which I give a quick back-alley treatment of the subject, covering just the bits you need to know.

To establish notation, we begin with the definition of a stochastic process.

Stochastic Processes

Definition (Stochastic process). Consider a probability space $\mathcal{S}=(\mathrm{\Omega },\mathcal{F},\mathrm{P})$ and a measurable space $(E,\mathcal{E})$. A stochastic process $f$ is a collection of $E$-valued random variables, all defined on the same probability space $\mathcal{S}$, where this collection is indexed by some set $\mathcal{X}$.

Note: to avoid unnecessary abstraction at this point, you may wish to substitute real-valued random variables; $(E,\mathcal{E})=(\mathbb{R},\mathcal{B}(\mathbb{R}))$ and $d$-dimensional real-valued indices $\mathcal{X}={\mathbb{R}}^d$ every time these come up; since these will closely resemble the GP regression case later.

Let ${\text{f}}_x$ or $f(x)$ denote the random variable in the collection $f$ indexed by $x$. That is, by the definition of a random variable, $f(x)$ is itself a function $f(x):\mathrm{\Omega }\to E$ of some sample $\omega \in \mathrm{\Omega }$. Outcomes of a stochastic process are individually functions of both the index of the random variable $x\in \mathcal{X}$, as well as a sample $\omega \in \mathrm{\Omega }$. Hence, it is common to write a stochastic process as the following mapping:

$$f:\mathcal{X}\times \mathrm{\Omega }\to E.$$

In this sense, the fact that $f$ accepts an index parameter $x$, in addition to the sample parameter $\omega$ that is ordinarily the only parameter to a random variable, allows us to see $f$ as a random function or a distribution over functions. The result is that we can sample functions from our stochastic process just as we would sample values from a random variable.

To do this, we fix some $\omega \in \mathrm{\Omega }$ at which we will evaluate all the random variables in our collection $f$ that we decide to look at (i.e. that we decide to index). That is, the sampled function is

$$f(\cdot ,\omega ):\mathcal{X}\to E,$$

where by selecting some $\omega$, we have ‘fixed’ the random function 4: In practice, when using stochastic processes, we almost never fix some $\omega$ unless we are plotting sampled functions from the process. Instead, it is usually more interesting to consider the distribution (and moments) of the function points, as well as evaluating the likelihood of new function points under the random function. ^4[4] , which now behaves like any ordinary deterministic function. The indices $x\in \mathcal{X}$ can be interpreted more directly as ‘function inputs’ which are mapped to single points in the output space $E$.

A stochastic process is characterised by the joint probability distribution of $f$ on arbitrary finite subsets of $\mathcal{X}$. That is, if ${\mathbf{X}}_n\subseteq \mathcal{X}$ is a finite $n$-element vector of indices, then we obtain the following random vector

$${\mathbf{f}}_{{\mathbf{X}}_n}=({\text{f}}_{x_1},\dots ,{\text{f}}_{x_n})=(f(x_1),\dots ,f(x_n))=f({\mathbf{X}}_n).$$

More formally, we refer to these as finite-dimensional distributions.

Definition (Finite-dimensional distribution). Consider a probability space $(\mathrm{\Omega },\mathcal{F},\mathrm{P})$, measurable space $(E,\mathcal{E})$ and a stochastic process $f:\mathcal{X}\times \mathrm{\Omega }\to E$. The finite-dimensional distributions (FDDs) of $f$ are the (joint, $n$-dimensional) probability distributions ${\mathrm{P}}_{x_1,\dots ,x_n}^f$ on the product space $E^n$ defined by

$${\mathrm{P}}_{x_1,\dots ,x_n}^f(S)\doteq \mathrm{P}(\{\omega \in \mathrm{\Omega }:(f(x_1,\omega ),\dots ,f(x_n,\omega ))\in S\}),$$

where $S\in E^n$ and $n\in \mathbb{N}$.

While the above is sufficiently abstract to entertain a variety of Gaussian processes (e.g. over indices $\mathcal{X}$ which are graphs or strings), for most straightforward applications of GP regression, the following will suffice:

• We consider real-valued random variables, that is $(E,\mathcal{E})=(\mathbb{R},\mathcal{B}(\mathbb{R}))$ in which case each random variable ${\text{f}}_x$ indexed by some $x\in \mathcal{X}$ maps some element of the sample space to the reals: ${\text{f}}_x:\mathrm{\Omega }\to \mathbb{R}$.
• The collection of random variables that makes up the stochastic process can be denoted $\{f(x):x\in \mathcal{X}\}$, where in GP regression we often have $\mathcal{X}={\mathbb{R}}^n$, for $n\ge 1$. In this case, since $\mathcal{X}$ is uncountable, then so is this collection of random variables, and we call the stochastic process a continuous stochastic process.

Modelling Functions with GPs

For the purposes of modelling functions with Gaussian processes 5: or any stochastic process for that matter… ^5[5] , we interpret the index set $\mathcal{X}$ as the domain or inputs to our function.

For instance, for time-series prediction the indices may be a set of time values we have in a training dataset, as well as time values we’d like to make predictions at. For Bayesian optimisation of an ML model’s hyperparameters, the indices may be the set of hyperparameters we’ve already tried, as well as the ones we’re interested in trying next. For plotting a 2d function on a graph, the indices are the real values along the x-axis that we’d like to plot the function at.

More generally, the indices are the vector-valued covariates for our prediction problem 6: It is easiest to think about regression for now, although we will see later how GPs can be used for classification by changing the likelihood. ^6[6] for which we have training observations 7: Note that we needn’t have observed any observations to use a GP; indeed Bayesian methods shine for how well they do a-priori and when we have observed very few data points ^7[7] as well as those for which we’d like to make predictions (obtain a predictive distribution).

Single Sampled Function

To see how the GP corresponds to a collection of random variables at the points of interest (i.e. train and test points) which are in some sense consistent with one another, we will build up to this one index at a time, working with a simple 1d function $f:\mathbb{R}\to \mathbb{R}$, having already fixed some value in the sample space, ${\omega }_1\in \mathrm{\Omega }$ 8: The effect of this ${\omega }_1$ is to ‘pin down’ the random function, which we may now think of as a deterministic function, whose points can be described by a Gaussian with a particular mean and covariance. ^8[8] .

Note that in this section we are only interested in seeing how GPs can be used to represent functions. We leave the form of $f$ unspecified, as well as how we actually find the parameters of the GP’s FDD; which will be the subject of the next section.

Picking an initial index value ${\mathbf{x}}_1\in \mathcal{X}$ gives a single element subset of the GP’s infinite collection of real-valued random variables, $\{{\text{f}}_{{\mathbf{x}}_1}\}\subset \{{\text{f}}_{\mathbf{x}}:\mathbf{x}\in \mathcal{X}\}$. The random variable in this singleton set is univariate Gaussian

$$\begin{array}{}\text{(1)}& {\text{f}}_{{\mathbf{x}}_1}\sim \mathcal{N}({\mu }_1,{\mathrm{\Sigma }}_1).\end{array}$$

Having already fixed an element of the sample space, ${\omega }_1$, we can evaluate this random variable ${\text{f}}_{{\mathbf{x}}_1}:\mathrm{\Omega }\to \mathbb{R}$ at ${\omega }_1$ to get a real value which we may plot as the first function point / ‘pixel’ of our sampled random function:

Adding a second index ${\mathbf{x}}_2$ to the index set, $\{{\mathbf{x}}_1,{\mathbf{x}}_2\}\in \mathcal{X}$, we now obtain a two-element set of random variables $\{{\text{f}}_{{\mathbf{x}}_1},{\text{f}}_{{\mathbf{x}}_2}\}\subset \{{\text{f}}_{\mathbf{x}}:\mathbf{x}\in \mathcal{X}\}$. By the definition of a GP as a stochastic process with Gaussian FDDs, these two random variables form a Gaussian random vector distributed as

$${\left[\begin{array}{cc}{\text{f}}_{{\mathbf{x}}_1}& {\text{f}}_{{\mathbf{x}}_2}\end{array}\right]}^{\mathrm{\top }}=\mathcal{N}(\mathit{\mu },\mathbf{\Sigma }),$$

for some $\mathit{\mu }\in {\mathbb{R}}^2$ and $\mathbf{\Sigma }\in {\mathbb{R}}^{2\times 2}$, whose exact form we will cover in the next section. For now, note that the fact that ${\text{f}}_{{\mathbf{x}}_1}$ and ${\text{f}}_{{\mathbf{x}}_2}$ are elements of the same joint Gaussian distribution is the sense in which the predictions at ${\mathbf{x}}_1$ and ${\mathbf{x}}_2$ are consistent. Marginalising 9: Marginalising a Gaussian simply reduces to picking the relevant sub-vector and sub-matrix of the mean and covariance. See this section in my previous article on the Gaussian for more on this. ^9[9] this 2-dimensional FDD further gives ${\text{f}}_{{\mathbf{x}}_1}\sim \mathcal{N}({\mu }_1,{\mathrm{\Sigma }}_1)$ and ${\text{f}}_{{\mathbf{x}}_2}\sim \mathcal{N}({\mu }_2,{\mathrm{\Sigma }}_2)$, where the former is exactly the univariate Gaussian from Equation $\text{1}$.

Once again, having already fixed ${\omega }_1\in \mathrm{\Omega }$, we can evaluate both ${\text{f}}_{{\mathbf{x}}_1}$ and ${\text{f}}_{{\mathbf{x}}_2}$ at ${\omega }_1$ which gives us two real-valued function points / ‘pixels’ which lie on the same sampled function line:

This does not offer us a very high ‘resolution’ for plotting our function. Indexing using 100 additional points, spaced from $[-10,10]$, evaluating the resulting 100-dimensional Gaussian random vector at ${\omega }_1$ and plotting these points lets us see the rest of the sampled function more clearly:

Hence indexing results in a Gaussian random vector, and sampling from this gives us consistent realisations of the random function evaluated at the indices / inputs.

Distributions over Functions

If we now refrain from picking a single element of the sample space $\mathrm{\Omega }$, we can consider the distribution over functions.

Working up to this once more from a single index, the univariate Gaussian FDD at ${\mathbf{x}}_1$ is ${\text{f}}_{{\mathbf{x}}_1}\sim \mathcal{N}({\mu }_1,{\mathrm{\Sigma }}_1)$. This is just a univariate Gaussian distribution, which we will suggestively plot vertically on the function plot, along with five sampled function points from this univariate Gaussian:

Continuing as we did before by adding a second index to the index set, ${\mathbf{x}}_2$, this yields a second marginal ${\text{f}}_{{\mathbf{x}}_2}$ of the now 2-dimensional joint, which we will also plot vertically at the corresponding index location ${\mathbf{x}}_2$. Under the density, we have also plotted 5 samples using the same $\{{\omega }_1,\dots ,{\omega }_5\}$ as were used to draw the samples at ${\mathbf{x}}_1$.

Continuing in this way by adding more and more indices to the index set, we can draw out the full function lines 10: Note that when using GPs, we never actually do this unless we’re trying to visualise the GP. Instead, we only need to add the training points and test points to the index set. ^10[10] for our five elements of the sample space $\{{\omega }_1,\dots ,{\omega }_5\}\subset \mathrm{\Omega }$, for some index / design matrix $\mathbf{X}\in {\mathbb{R}}^{n\times 1}$.

You’ll notice that we stopped visualising the distribution over functions at the various indices using vertical Gaussian densities at each input location, since this quickly becomes messy. Visualising the distribution by drawing a bunch of functions (as we have done above, using 5 samples) is equally messy. Rather, to represent the broadness of these vertical Gaussians at various inputs and communicate our uncertainty about the function values at those points, it is common to plot ${\mu }_{\mathbf{x}}\pm 2{\sigma }_{\mathbf{x}}$ for all $\mathbf{x}\in \mathbf{X}$.

For the GP prior that we have been sampling from in the plots above (i.e. the GP without conditioning any training data) this is usually fairly boring; with a constant mean and standard deviation / uncertainty at all points, which correspond to how we have defined the prior (see next section).

Once we have conditioned on training data however, we may observe varying degrees of uncertainty in our random function at different points away from the data.

In general, to define a Gaussian it suffices to define its first two moments: the mean and covariance. For any $n$ indices $\mathbf{X}=\{{\mathbf{x}}_i\in \mathcal{X}{\}}_{i=1}^n$, the resulting FDD is a $n$-dimensional Gaussian, parametrised by $\mu \in {\mathbb{R}}^n$ and $\mathbf{\Sigma }\in {\mathbb{R}}^{n\times n}$. In the next section, we will see how to obtain these parameters for a given GP prior, a training dataset, and some test points of interest.

Defining the GP Prior

When we place a Gaussian process prior on a function that we’re interested in modelling, we define a way of obtaining the first two moments of the FDD; in other words, we define a mean function $m:\mathcal{X}\to \mathbb{R}$ and a positive definite covariance (or kernel) function $k:\mathcal{X}\times \mathcal{X}\to \mathbb{R}$.

This is convenient since we often don’t know which indices $\mathbf{X}$ we’ll choose to look at when choosing the GP prior 11: i.e. where our training points will land, and where we will need to evaluate the function at test time ^11[11] . With these two functions, we can postpone the calculation of the FDD’s first two moments until we obtain the indices we need.

Supposing that our indices or covariates live in the space $\mathcal{X}={\mathbb{R}}^d$, and we select / observe $n$ of them $\mathbf{X}\in {\mathbb{R}}^{n\times d}$, the resulting prior FDD can be found by the element-wise evaluation of these functions at each of the indicies; giving us our mean vector $m(\mathbf{X})=\mathit{\mu }\in {\mathbb{R}}^n$ and covariance or Gram matrix $\mathbf{K}$, where ${\text{K}}_{i,j}=k({\mathbf{x}}_i,{\mathbf{x}}_k)$. The positive definite constraint on the kernel function says that the Gram matrices it produces satisfy ${\mathbf{v}}^{\mathrm{\top }}\mathbf{K}\mathbf{v}>0$, which is required to use these as covariance matrices. Drawing $n$-dimensional samples from $\mathcal{N}(\mathit{\mu },\mathbf{K})$ results in (prior) function samples, as plotted in the penultimate digram above.

The mean and covariance functions are the main levers we can pull to influence the expression of our prior beliefs about the function of interest. In most cases, using a zero-mean function $m(\mathbf{x})\doteq \mathbf{0}$ is a reasonable choice 12: although this is not always the case; a non-zero mean functions can be extremely useful, if not essential, in certain settings such as Bayesian optimisation, time-series predictions and deep Gaussian processes ^12[12] leaving the characterisation of the function down to the choice of kernel function. Different kernel functions can influence things such as the smoothness or periodicity of the resulting functions—for a nice taxonomy of some common kernels, see this page by David Duvenaud.

Some popular kernels include the following, where $r=|\mathbf{x}-{\mathbf{x}}^{\prime }|$ represents a pairwise distance between points 13: See this Fragment for more on computing pairwise distances between datapoints. ^13[13]

Note that for the ArcCos kernel, we have $\cos \theta =\frac{\mathbf{x}\cdot {\mathbf{x}}^{\prime }}{|\mathbf{x}||{\mathbf{x}}^{\prime }|}$. The $\sigma$ and $\ell$ terms are hyperparameters which influence the properties of each kernel; we discuss how to set these later.

In summary, we define a Gaussian process prior $f:{\mathbb{R}}^n\times \mathrm{\Omega }\to \mathbb{R}$, which we denote as

$$f\sim \mathcal{G}\mathcal{P}(m,k),$$

by defining the functions $m(\cdot )$ and $k(\cdot ,\cdot )$, which in turn determine the first two moments of the resulting FDDs:

$$\begin{array}{rl}m(\mathbf{x})& =\mathbb{E}[f(\mathbf{x})],\\ k(\mathbf{x},{\mathbf{x}}^{\prime })& =\mathbb{E}[(f(\mathbf{x})-\mu (\mathbf{x}){)}^{\mathrm{\top }}(f({\mathbf{x}}^{\prime })-\mu ({\mathbf{x}}^{\prime }))].\end{array}$$

Conditioning on Data

Defining a GP prior and drawing samples from it is not usually very useful in itself: rather we seek to make predictions about new inputs by updating or conditioning our beliefs about the latent function values on observed (training) data points. In this section, we will see how we can obtain the updated parameters of the FDD and parametrise the posterior predictive / posterior process; as was done to produce the last plot in the previous section.

For clarity of exposition, we will first assume a noiseless regression setup 14: These might come up when modelling the output of a computer simulation, or something for which we have extremely precise measurements. ^14[14] .

Consider a training set of $n$ input and function-value pairs, $\mathcal{D}=\{({\mathbf{x}}_i,f_i){\}}_{i=1}^n$, whose $d$-dimensional inputs we collectively denote as $\mathbf{X}\in {\mathbb{R}}^{n\times d}$ and where the function points are denoted as $\mathbf{f}=f(\mathbf{X})\in {\mathbb{R}}^n$. Further consider a set of $m$ new indices / test points, which we denote ${\mathbf{X}}_{\star }\in {\mathbb{R}}^{m\times d}$ at which we’d like to make predictions ${\mathbf{f}}_{\star }=f({\mathbf{X}}_{\star })$.

With these input points defined we can compute the GP’s FDD at $\{\mathbf{X},{\mathbf{X}}_{\star }\}$ by evaluating the mean and covariance functions at these input points. This gives us a joint distribution over the function values at the training and test locations $p(\mathbf{f},{\mathbf{f}}_{\star })$. Writing these parameters as a partitioned mean vector and covariance matrix, we obtain the distribution over the Gaussian random vector $\left[\begin{array}{c}\mathbf{f},{\mathbf{f}}_{\star }\end{array}\right]$ as

$$\left[\begin{array}{c}\mathbf{f}\\ {\mathbf{f}}_{\star }\end{array}\right]\sim \mathcal{N}(\left[\begin{array}{c}\mathit{\mu }\\ {\mathit{\mu }}_{\star }\end{array}\right],\left[\begin{array}{cc}\mathbf{K}& {\mathbf{K}}_{\star }\\ {\mathbf{K}}_{\star }^{\mathrm{\top }}& {\mathbf{K}}_{\star \star }\end{array}\right]),$$

where $\mathbf{K}=k(\mathbf{X},\mathbf{X})\in {\mathbb{R}}^{n\times n}$, ${\mathbf{K}}_{\star }=k(\mathbf{X},{\mathbf{X}}_{\star })\in {\mathbb{R}}^{n\times m}$, ${\mathit{\mu }}_{\star }=m({\mathbf{X}}_{\star })\in {\mathbb{R}}^m$ and so forth.

Given this joint model, we now note that while we know the function values at the training inputs $\mathbf{f}$ (these are in our training set), we don’t know the function values at the test points ${\mathbf{f}}_{\star }$—this is what we would like to infer.

To make predictions, we therefore seek the predictive distribution of function values at the test points, conditioned on the observed training dataset, $p({\mathbf{f}}_{\star }|\mathbf{f})$ 15: or, using more explicit notation, $p({\mathbf{f}}_{\star }|\mathcal{D})$ or even $p({\mathbf{f}}_{\star }|{\mathbf{X}}_{\star },\mathbf{f},\mathbf{X})$. However, it is common to leave out the indices for notational brevity, with the understanding that $\mathbf{f}\doteq f(\mathbf{X})$. ^15[15] .

Owing to the convenient properties of the Gaussian, this Gaussian conditional distribution has a closed-form solution, which I give a full derivation of in this section of my previous article on the Gaussian distribution. If you would find it satisfying to know how GPs work end-to-end, I encourage you to read through the proof. Otherwise, I have provided the relevant results below which you may use ad-hoc:

Hence the Gaussian random vector ${\mathbf{f}}_{\star }$ conditioned on our training points is distributed as $p({\mathbf{f}}_{\star }|{\mathbf{X}}_{\star },\mathbf{f},\mathbf{X})=\mathcal{N}({\mathbf{f}}_{\star };{\mathit{\mu }}_{\star |\mathbf{X}},{\mathbf{\Sigma }}_{\star |\mathbf{X}})$ with parameters

$$\begin{array}{}\text{(4)}& {\mathit{\mu }}_{\star |\mathbf{X}}& ={\mathit{\mu }}_{\star }+{\mathbf{K}}_{\star }^{\mathrm{\top }}{\mathbf{K}}^{-1}(\mathbf{f}-\mathit{\mu })\text{(5)}& {\mathbf{\Sigma }}_{\star |\mathbf{X}}& ={\mathbf{K}}_{\star \star }-{\mathbf{K}}_{\star }^{\mathrm{\top }}{\mathbf{K}}^{-1}{\mathbf{K}}_{\star },\end{array}$$

where we have used Equations $\text{2}$ and $\text{3}$ for the mean and covariance, respectively.

These are the equations that can be used 16: I refrain from saying ‘that were used’, since in practice these equations prove to be numerically unstable. I discuss some methods for practical inversion of kernel matrices in this section of another article on Bayesian linear regression. ^16[16] to plot the predictive GP conditioned on some training data:

Notice how the mean of the posterior perfectly interpolates the training data, and the predictive uncertainty ‘pinches’ around the training points. This is an artefact of the noiseless modelling assumption. In most cases when modelling observed phenomena this is unrealistic, and this may further have undesired effects on the predicted function. For instance if two similar inputs have differing observations (e.g. due to corrupting observation noise on a sensor), then the predicted function (distribution) will oscillate wildly around that point to match the data:

Noise-Corrupted Observations

To rectify the above, we now assume that instead of observing noiseless function values $f({\mathbf{x}}_n)$ in our training set, we rather observe noise corrupted versions thereof; ${\text{y}}_n=f({\mathbf{x}}_n)+{ϵ}_n$, for ${ϵ}_n\sim \mathcal{N}(0,{\sigma }_y^2)$.

Note that this unimodal, Gaussian noise modelling assumption is required in the exact GP setting to remain tractable, however in a subsequent article we will see how we can relax this constraint and use different likelihoods.

Due to our assumption that the noise has zero mean (i.e. is centred at the latent function values $\mathbf{f}$) and i.i.d. (i.e. it has a diagonal covariance), incorporating this modelling assumption is fairly straightforward; only requiring us to update the diagonal entries of the covariance matrix at the training locations 17: we don’t update the test locations since we are concerned with inferring the latent function values ${\mathbf{f}}_{\star }$, not noisy versions thereof! ^17[17] :

$$\mathrm{Cov}[\mathbf{y}|\mathbf{X}]=\mathbf{K}+{\sigma }_y^2{\mathbf{I}}_n.$$

Hence, our model for the joint density of our observations (which we now denote $\mathbf{y}$) and the latent function values at the test inputs ${\mathbf{f}}_{\star }$ becomes:

$$\begin{array}{}\text{(6)}& \left[\begin{array}{c}\mathbf{y}\\ {\mathbf{f}}_{\star }\end{array}\right]\sim \mathcal{N}(\left[\begin{array}{c}\mathit{\mu }\\ {\mathit{\mu }}_{\star }\end{array}\right],\left[\begin{array}{cc}\mathbf{K}+{\sigma }_y^2\mathbf{I}& {\mathbf{K}}_{\star }\\ {\mathbf{K}}_{\star }^{\mathrm{\top }}& {\mathbf{K}}_{\star \star }\end{array}\right]).\end{array}$$

Conditioning on the observations $\mathbf{y}$ proceeds exactly as in the noiseless case (Equations $\text{2}$ and $\text{3}$) with $\mathbf{K}+{\sigma }_y^2\mathbf{I}$ substituted for $\mathbf{K}$, giving the posterior predictive distribution in the noisy case as:

$$\begin{array}{}\text{(7)}& p({\mathbf{f}}_{\star }|{\mathbf{X}}_{\star },\mathbf{y},\mathbf{X})& =\mathcal{N}(\mathbf{f};{\mathit{\mu }}_{\star |\mathbf{X}},{\mathbf{\Sigma }}_{\star |\mathbf{X}})\text{(8)}& {\mathit{\mu }}_{\star |\mathbf{X}}& ={\mathit{\mu }}_{\star }+{\mathbf{K}}_{\star }^{\mathrm{\top }}(\mathbf{K}+{\sigma }_y^2\mathbf{I}{)}^{-1}(\mathbf{y}-\mathit{\mu })\text{(9)}& {\mathbf{\Sigma }}_{\star |\mathbf{X}}& ={\mathbf{K}}_{\star \star }-{\mathbf{K}}_{\star }^{\mathrm{\top }}(\mathbf{K}+{\sigma }_y^2\mathbf{I}{)}^{-1}{\mathbf{K}}_{\star }.\end{array}$$

Setting the noise level to ${\sigma }_y^2=0.2$, we can now see that the mean function values now merely pass close to the training datapoints, without the hard requirement to interpolate through them:

Due to this behaviour, these simple, single-layer GPs are often referred to as smoothing devices.

Appendix A: The Measure Theory You Need To Know (and no more)

This Appendix collects the bits of basic measure theory which are useful for understanding GPs. It reads like a list of definitions, which you should aim to commit to memory.

Definition (Sample Space). A sample space $\mathrm{\Omega }$ is the set of all outcomes of a probabilistic experiment. More generally, a measurable space, commonly denoted $E$, is the set of things which we are interested in measuring. A sample space is a measurable space 18: The former arising when talking about probability, with the latter arising more generally in measure theory. ^18[18] .

Example. For a single toss of a coin, the sample space is $\mathrm{\Omega }=\{{\omega }_1,{\omega }_2\}=\{H,T\}$.

Sigma Algebras

Definition (Sigma-algebra, event space). These are the same thing, with the latter being a probabilistic interpretation of the former. A $\sigma$-algebra $\mathcal{E}$ is a collection of ‘nice’ (i.e. not infinite, tractable) subsets of a measurable space. If the measurable space is finite and discrete, then its power set $\mathcal{P}(E)$ can be used as a valid sigma algebra. In the probabilistic interpretation, for a finite sample space $\mathrm{\Omega }$, the event space (which we usually denote $\mathcal{F}$) can be taken as a subset of $\mathcal{P}(\mathrm{\Omega })$; elements of which are called events. Elements of a sigma-algebra must follow the following regularity conditions:

• $\mathbb{E}\in \mathcal{E}$; the sample space itself is an event (the sure event)
• $A\in \mathcal{E}⟹A^c\in \mathcal{E}$; it is closed under complementation (the complement of an event is an event)
• $A_1,A_2,\cdots \in \mathcal{E}⟹\bigcup _{i=1}^{\mathrm{\infty }}{A}_i\in \mathcal{E}$; it is closed under countable unions (the union of countably many events is an event)

If a measurable space $E$ is countable, then this is nice because we can always find the power set $\mathcal{E}=\mathcal{P}(E)$ and the measure can be defined on all possible subsets. However, for ‘bigger’ spaces, we have to be more careful, since the set of all subsets if often too large to work with. For such inconvenient, infinite event spaces, we usually define $\sigma$-algebras in terms of a smaller set, known as the generating sets.

Definition (Generator of a $\sigma$-algebra). Let $E$ be a set, and $\mathcal{A}\subseteq \mathcal{P}(E)$ be a collection of subsets of $E$. We define

$$\sigma (\mathcal{A})=\{A\subseteq E:A\in \mathcal{A}\text{ for all }\sigma \text{-algebras }\mathcal{E}\text{ that contain }\mathcal{A}\}.$$

Put otherwise, $\sigma (\mathcal{A})$ is the smallest sigma algebra that contains $\mathcal{A}$. This is known as the $\sigma$-algebra generated by $\mathcal{A}$.

Definition (Borel $\sigma$-algebra). Consider a real measurable space $E=\mathbb{R}$ and the generator $\mathcal{A}=\{U\subseteq \mathbb{R}:U\text{ is open}\}$. Then $\sigma (\mathcal{A})$ is known as the Borel $\sigma$-algebra, which is not the set of all subsets of $\mathbb{R}$ (which would be an uncountably infinite object, and hence hard to measure).

We can perhaps more intuitively, and entirely equivalently, define this as $\overline{\mathcal{A}}=\{(a,b):a<b,a,b\in \mathrm{Q}\}$—the set of tuples of ordered rational numbers. Then $\sigma (\overline{\mathcal{A}})$ is also the Borel $\sigma$-algebra. We denote the Borel $\sigma$-algebra as $\mathcal{B}(\mathbb{R})$.

Measures

Definition (Measurable space). A pair $(E,\mathcal{E})$ of a set and its $\sigma$-algebra (for instance a sample space $\mathrm{\Omega }$ and an event space $\mathcal{F}$) is called a measurable space. This is because we can define a function that ‘measures’ the size of the elements of $\mathcal{F}$. In fact, any element of a sigma algebra $B\in \mathcal{F}$ is called a measurable set.

Definition (Measure). A measure on a measurable space $(E,\mathcal{E})$ is a function $\mu :\mathcal{E}\to [0,+\mathrm{\infty })$, such that

1. $\mu (\mathrm{\varnothing })=0$
2. It follows countable additivity: for any disjoint sequence $(A_n)$ in $\mathcal{E}$, then

$$\mu \left(\bigcup _n{A}_n\right)=\sum _{n=1}^{\mathrm{\infty }}\mu (A_n).$$

Example (Mass function). Let $E$ be any countable set, and $\mathcal{E}=\mathcal{P}(E)$ its power set (i.e. the set of all subsets of $E$). A mass function is any function $m:E\to [0,+\mathrm{\infty })$. We can define a measure by setting

$$\mu (A)=\sum _{x\in A}m(x).$$

In particular, if we put $m(x)=1$ for all $x\in E$, then we obtain the counting measure.

Definition (Probability measure). Let $(\mathrm{\Omega },\mathcal{F})$ be a measurable space with the property that $\mu (\mathrm{\Omega })=1$. Then we call $\mu$ a probability measure and usually denote it as $\mathrm{P}$ instead of $\mu$. Note that $\mathrm{P}$ must also satisfy the additional property that it maps to $[0,1]$ rather than $[0,+\mathrm{\infty })$.

Definition (Lebesgue measure). The Lebesgue measure is the unique Borel measure $\mu$ on $\mathbb{R}$ with $\mu ([a,b])=b-a$. Note that the Lebesgue measure is not a finite measure, since $\mu (\mathbb{R})=\mathrm{\infty }$, however it is a $\sigma$-finite measure (defined below).

We often denote the Lebesgue measure as $\lambda$ (with its $n$-dimensional generalisation on $\mathcal{B}({\mathbb{R}}^n)$ denoted ${\lambda }^n$).

Definition ($\sigma$-finite measure). Let $(E,\mathcal{E})$ be a measurable space, and $\mu$ a measure. We say that $\mu$ is $\sigma$-finite if there exists a sequence $(E_n)$ in $\mathcal{E}$ such that $\bigcup _n{E}_n=E$, and $\mu (E_n)<\mathrm{\infty }$ for all $n$.

Notice how the measure takes elements of $\mathcal{F}$ and not elements of the sample space: a common example of this is the probability space $(\mathbb{R},\mathcal{B}(\mathbb{R}),\mathrm{P})$, where $\mathcal{B}(\mathbb{R})$ is a Borel sigma-algebra and $\mathrm{P}$ is the Lebesgue measure. We always evaluate a probability density function on an interval (i.e. an element of $\mathcal{B}(\mathbb{R})$) and never directly on a real number!

Probability

Definition (Probability space). A probability space is a triple $(\mathrm{\Omega },\mathcal{F},\mathrm{P})$, composed of

• the sample space $\mathrm{\Omega }$,
• an event space $\mathcal{F}$ which is a $\sigma$-algebra on $\mathrm{\Omega }$, and
• a probability measure $\mathrm{P}$ on $(\mathrm{\Omega },\mathcal{F})$.

Definition (Random variable). Consider two measurable spaces: the first is the pair of a sample space and event space $(\mathrm{\Omega },\mathcal{F})$, and the other is an arbitrary pair of a set (a measurable space) and its sigma algebra $(E,\mathcal{E})$. A random variable (also called a measurable function 19: Note that when we talk about a measurable function, we will more often use the notation $f$ rather than $X$. ^19[19] ) is a function $X:\mathrm{\Omega }\to E$ which maps elements in the sample space to elements in the arbitrary set $E$, such that for any $\mathcal{E}$-measurable set $B\in \mathcal{E}$, the preimage $X^{-1}(B)$ is an $\mathcal{F}$-measurable set:

$$X^{-1}(B)=\{\omega \in \mathrm{\Omega }:X(\omega )\in B\}\in \mathcal{F}\mathrm{\forall }B\in \mathcal{E}.$$

$◼$

Intuitively, this is useful because the actual sample space and event space of our experiments, for instance $E$ and $\mathcal{E}$ might be tricky to work with, and we’d rather use alternative ones that we’re more comfortable using such as $\mathrm{\Omega }$ and $\mathcal{F}$ (which, concretely, might be $\mathbb{R}$ and $\mathcal{B}(\mathbb{R})$). The property above guarantees that there is some sort of correspondence between the events in our original event space $\mathcal{F}$ and our transformed event space $\mathcal{E}$.

On terminology, we often say that an $E$-valued random variable is a measurable function $X:\mathrm{\Omega }\to E$.

Example. If $E=\mathbb{R}$ and $\mathcal{E}=\mathcal{B}(\mathbb{R})$ a real-valued random variable (sometimes $\mathbb{R}$-valued random variable) assigns a real number to each $\omega \in \mathrm{\Omega }$. More generally, we might take $E={\mathbb{R}}^n$ and $\mathcal{E}=\mathcal{B}({\mathbb{R}}^n)$, where the random variable now maps elements of the sample space to an $n$-dimensional Euclidean vector.

Definition (Probability distribution / law). Suppose we have a probability space $(\mathrm{\Omega },\mathcal{F},\mathrm{P})$ allowing us to assign probabilities to events in $\mathcal{F}$. Suppose that we also have another measurable space $(E,\mathcal{E})$ for which we don’t yet have a probability measure to measure events in $\mathcal{E}$. A probability distribution (also a law or even the pushforward of a measure, denoted ${\mu }_x=\mathrm{P}\circ X^{-1}$ or $X_{\mathrm{\#}}\mathrm{P}$, respectively) is an image measure: the composition of the inverse of the measurable function defining the random variable $X$ with the probability measure $\mathrm{P}$. It is effectively a mapping from sets in $\mathcal{E}$ into $[0,1]$:

$$X_{\mathrm{\#}}\mathrm{P}=(\mathrm{P}\circ X^{-1}):\mathcal{E}\to [0,1].$$

The key idea is that, given some set (read event) in $\mathcal{E}$, we can use a random variable $X$ to find the pre-image of this set in the event space $\mathcal{F}$, and then measure this set via the probability measure $\mathrm{P}$ defined on this set. The existence of the random variable $X:\mathrm{\Omega }\to \mathcal{E}$ allows us to measure events $A\in \mathcal{E}$ not in the original probability space.

Using more probabilistic notation, we would write

$$\mathrm{P}(X\in A)={\mu }_x(A)=\mathrm{P}(X^{-1}(A))=\mathrm{P}(\{\omega \in \mathrm{\Omega }:X(\omega )\in A\}).$$

The above works well for countable $E$ (that is, discrete random variables). For uncountable sample spaces, such as when $E=\mathbb{R}$, then ${\mu }_X$ is determined by its values on the intervals $(-\mathrm{\infty },y]$. We set

$$F_X(x)={\mu }_X((-\mathrm{\infty },x])=\mathrm{P}(X\le x),$$

where $F_X(\cdot )$ is the distribution function of the random variable $X$.

One common example of this is the measurable space $(E,\mathcal{E})=({\mathbb{R}}^n,\mathcal{B}({\mathbb{R}}^n))$, where we have

$$X_{\mathrm{\#}}\mathrm{P}:\mathcal{B}({\mathbb{R}}^n)\to [0,1].$$

Definition (Probability density function). Consider the measurable space $({\mathbb{R}}^n,\mathcal{B}({\mathbb{R}}^n))$, and suppose that we have the usual Lebesgue measure ${\lambda }^n$ on this measurable space. This allows us to assign the notion of ‘size’ to each set $A\in \mathcal{B}({\mathbb{R}}^n)$ based on the likelihood of it happening. To do this, we assign a density to every point in ${\mathbb{R}}^n$, and to compute the new measure of $A$ according to $X_{\mathrm{\#}}\mathrm{P}$, we simply take the average of its density at all points in $A$, according to the Lebesgue measure:

$$X_{\mathrm{\#}}\mathrm{P}(A)={\int }_{A}d{X}_{\mathrm{\#}}\mathrm{P}={\int }_{A}fd{\lambda }^n\mathrm{\forall }A\in \mathcal{B}({\mathbb{R}}^n).$$