Lab 1: Gaussian Process Regression¶

Gaussian Process Summer School 2019¶

This lab is designed to introduce Gaussian processes in a practical way, illustrating the concepts introduced in the first two lectures. The key aspects of Gaussian process regression are covered: the covariance function (aka kernels); sampling a Gaussian process; and the regression model. The notebook will introduce the Python library GPy$^\dagger$ which handles the kernels, regression and optimisation of hyperparameter, allowing us to easily access the results we want.

The level of this notebook is aimed at introductory, as the background of attendees is diverse, and so cover a wide range of basic GP concepts. There are seven exercises to complete, the difficulty of which varies, but you should aim to complete all during the lab session. The notebook will not be marked and we will provide answers to the exercises after the lab session.

In addition, there is a second optional notebook with extra work. The content of this is more advanced, so completion is at your discretion.

$^\dagger$GPy: A Gaussian process framework in Python (since 2012). Available from http://github.com/SheffieldML/GPy

1. Getting started¶

[Binder] If you are using Binder, all libraries required are installed on the Binder image.

[Local] If you are running this notebook on a local machine, make sure that GPy is already installed on your machine. You should be using Python 3.5 or 3.6. A set of instructions for setting up your environment are available from the GPSS site.

First, we need to setup our notebook with the libraries we are going to use. We will use numpy for maths functionality, pyplot for plotting, and of course GPy for Gaussian processes.

The documentation for GPy is available at gpy.readthedocs.io. We will be using GPy to define our kernels, and regression. Note that GPy also contains plotting utilities, but we will not use these in this lab.

Covariance functions, aka kernels¶

We will define a covariance function, from hereon referred to as a kernel, using GPy. The most commonly used kernel in machine learning is the Gaussian-form radial basis function (RBF) kernel. It is also commonly referred to as the exponentiated quadratic or squared exponential kernel – all are equivalent.

The definition of the (1-dimensional) RBF kernel has a Gaussian-form, defined as:

It has two parameters, described as the variance, $\sigma^2$ and the lengthscale $\mathscr{l}$.

In GPy, we define our kernels using the input dimension as the first argument, in the simplest case input_dim=1 for 1-dimensional regression. We can also explicitly define the parameters, but for now we will use the default values:

We can see from the above table that our kernel has two parameters, variance and lengthscale, both with value 1.0. There is also information on the constraints and priors on each parameter, but we will look at this later.

Visualising the kernel¶

We can visualise our kernel in a few different ways. We can plot the shape of the kernel by plotting $k(x,0)$ over some sample space $x$ which, looking at the equation above, clearly has a Gaussian shape. This describes the covariance between each sample location and $0$.

Alternatively, we can construct a full covariance matrix, $\mathbf{K}_{xx} \triangleq k(x,x')$ with samples $x = x'$. The resulting GP prior is a multivariate normal distribution over the space of samples $x$: $\mathcal{N}(\mathbf{0}, \mathbf{K}_{xx})$. It should be evident then that the elements of the matrix represents the covariance between respective points in $x$ and $x'$, and that it is exactly $\sigma^2[=1]$ in the diagonal.

We can show this using pyplot to plot the vector $k(x,0)$ and the matrix $k(x,x')$ using k.K($\cdot$ , $\cdot$):

Setting the kernel parameters¶

Looking at the above definition of the RBF kernel, we can see that the parameters, i.e. variance and lengthscale, control the shape of the covariance function and therefore the value of the covariance between points $x$ and $x'$.

We can access the value of the kernel parameters in GPy and manually set them by calling k.param_name, e.g. k.lengthscale or k.variance for the RBF kernel. The following example demonstrates how the value of the lengthscale affects the RBF kernel:

Exercise 1¶

(a) What is the effect of the lengthscale parameter on the covariance function?

(b) Change the code used above to plot the covariance function showing the effects of the variance on the covariance function. Comment on the effect.

2. Types of covariance function¶

There are many different covariance functions already implemented in GPy. Aside from the RBF kernel, there are others such as the following:

• Exponential
• Matern32
• Matern52
• Brownian
• Bias
• Linear
• StdPeriodic
• Cosine
• PeriodicMatern32,

Note: when defining these, all are preceded by GPy.kern. The following are some examples of the Matérn 5/2 and Cosine kernels, compared with the RBF kernel:

Not every kernel has the same set of parameters. Some kernels are not parameterised by a lengthscale, for example, like the Linear kernel which only has a list of variances corresponding to each linear component

Likewise, not every kernel is stationary. In the case of the Gaussian RBF, or Matérn kernels, the kernel can be written $\kappa(x,x') = f(x-x')$, however this is not true for, e.g., the Brownian motion covariance function, which is defined as $k(x,x') = \min(x,x')$

3. Combining covariance functions¶

Exercise 2¶

(a) A matrix, $\mathbf{K}$, is positive semi-definite if the matrix inner product is greater than or equal to zero, $\mathbf{x}^\text{T}\mathbf{K}\mathbf{x} \geq 0$, regardless of the values in $\mathbf{x}$. Given this, it should be easy to see that the sum of two positive semi-definite matrices is also positive semi-definite. In the context of Gaussian processes, this is the sum of two covariance functions. What does this mean from a modelling perspective?

(b) What about the element-wise product of two covariance functions? If we define $k(\mathbf{x}, \mathbf{x}') = k_1(\mathbf{x},\mathbf{x}')k_2(\mathbf{x},\mathbf{x}')$, then is $k(\mathbf{x},\mathbf{x}')$ a valid covariance function?

Combining kernels in GPy¶

We can easily combine kernels using GPy using the + and * operators, respectively denoting addition and product of kernels.

Summing kernels¶

An example of adding kernels is shown here. We create a new kernel that is the sum of an RBF and a Matern 5/2 kernel.

We can visualise our kernel sum to see the resulting effect. It should be fairly clear that the result is simply the sum of evaluations of the respective kernels for each sample point.

Multiplying two kernels¶

We also demonstrate here the effect of multiplying two kernels. Here, we multiply an RBF and Periodic kernel, effectively encapsulating the periodicity into an RBF window:

4. Sampling from a Gaussian Process¶

A Gaussian process provides a prior over some infinite-dimensional function, defined by a mean function and covariance function

When we sample from the covariance function, $k$, to create a matrix over some sample space, we are creating a matrix of values that describe the covariance between sample points. Since it is not possible to sample every single point in an infinite dimensional function, we have to sample a finite subset of the input domain. Let $\mathbf{X}$ denote some sample inputs, and $\mathbf{K}$ the covariance matrix, with elements $K_{ij} = k(\mathbf{X}_i,\mathbf{X}_j)$, then we can describe the prior over $f(\mathbf{X})$ as a (finite-dimensional) normal distribution with covariance $\mathbf{K}$. As such, we can easily create samples of $f$ which, for a good choice of $\mathbf{X}$, are representative of the true function.

We can also sample from the kernel prior by creating a covariance matrix over a sample space and sampling from a zero-mean multivariate normal distribution with covariance $\mathbf{K}$. Below are examples of different kernels with different parameters, including composite kernels.

These samples are from the Gaussian process prior made up of the covariance function and a zero mean. After GP regression, the fitted posterior can also be sampled in this manner, to get samples of the fitted function.

Exercise 3¶

Can you identify the covariance function used to generate the following samples?

5. Gaussian Process Regression¶

We will now use our Gaussian process prior with some observed data to form a GP regression model.

Suppose we have a data model for which we only have noisy observations, $y = f(x) + \epsilon$ at some small number of sample locations, $\mathbf{X}$. Here, we set up an example function

We will first fit a Gaussian process using the exact equations.

A Gaussian process regression model using a Gaussian RBF covariance function can be defined first by setting up the kernel:

And then combining it with the data to form a Gaussian process regression model, with $\mathbf{X}^*$ representing any new inputs (imagine $\mathbf{f}^*$ approximates $f(\mathbf{X}^*)$):

where $ \mathbf{m} = \mathbf{K}_{*x}(\mathbf{K}_{xx} + \sigma^2\mathbf{I})^{-1}\mathbf{y}$ and $\mathbf{C} = \mathbf{K}_{**} - \mathbf{K}_{*x}(\mathbf{K}_{xx} + \sigma^2\mathbf{I})^{-1}\mathbf{K}_{*x}^\text{T} $ and covariance matrices are defined by evaluations of the kernel functions: $\mathbf{K}_{xx} = k(\mathbf{X}, \mathbf{X})$; $\mathbf{K}_{*x} = k(\mathbf{X}^*, \mathbf{X})$; and $\mathbf{K}_{**} = k(\mathbf{X}^*,\mathbf{X}^*)$.

Here we define a quick plotting utility function for our GP fits. There are a number of plotting options available in GPy, but we will use the below method, which plots the mean and 95% confidence fit of a GP for a given input $\mathbf{X}^*$. Optionally, we will allow it to plot the initial training points.

We can also save effort and time by to do Gaussian process regression using GPy, by creating a GP regression model with sample points $(\mathbf{X}, \mathbf{Y})$ and the Gaussian RBF kernel:

We can use GPy's regression and prediction tools, which should give the same result as our basic implementation:

It can be clearly seen that this is the same fit as our above model. However, using GPy gives flexibility and ease of use for extending the capabilities of the fitting, including use of different kernels, optimising parameters and solving more complicated problems, including classification. We also don't need to write explicit equations and manually creating covariances matrices.

Exercise 4¶

(a) What do you think of this initial fit? Does the prior given by the GP seem to be adapted?

(b) The parameters of the model can be editted much like those of the kernel. For example,

m.Gaussian_noise = 0.1

or

m.rbf.variance = 2.0

Change the values of the parameters to try and obtain a better fit of the GP. You can recalculate the updated mean and covariance after changing the values by calling m.predict_noiseless as above.

Note: changing the original kernel k will also affect the model parameters due to how Python connects objects, but this is not a reliable way of setting the parameters, so you should adjust the kernel parameters via the model m as described

(c) Given that we can obtain the mean and covariance of the GP fit, we can also sample the GP posterior as a multivariate Gaussian process. This can be done as in Section 4, where we sampled the priors as defined by the kernels, i.e. with np.random.multivariate_normal. Obtain 10 samples from the GP posterior and plot them alongside the data. Try to simulate noisy measurements using m.predict (rather than m.predict_noiseless).

Remember to get the full covariance matrix, using full_cov=True, and note that to make the mean vector 1-D (for sampling a multivariate normal), you need np.random.multivariate_normal(mean[:,0], Cov))

6. Covariance Function Parameter Estimation¶

As discussed in the lectures, the values of kernel parameters can be estimated by maximising the likelihood of the observations. This is useful to optimise our estimate of the underlying function, without eye-balling parameters to get a good fit.

In GPy, the model objects such as GPRegression, have parameter optimisation functionality. We can call this as following:

We can see the selected parameters in the model table above. The regression fit with the optimised parameters can be plotted:

Parameter constraints¶

We can see in the above model that the regression model is fit to the data, as the optimiser has minimised the noise effect in the model, Gaussian_noise.variance$ = 4.804\times10^{-8}$. If we know, or can reasonably approximate, the variance of the observation noise $\epsilon$, we can fix this parameter for the optimiser, using fix, which in the case of the above is $0.01$. We can also limit the values that the parameters take by adding constraints. For example, the variance and lengthscale can only be positive, so calling constrain_positive, we can enforce this (note that this is the default constraint for GP regression anyway).

We can see our constraints in the corresponding column in the above table, where "+ve" means we are constrained to positive values, and fixed means the optimiser will not try and optimise this parameter. We can see here that the variance of the noise in the model is unchanged by the optimiser. Looking at the resulting plot, we can see that we have a much more reasonable confidence in our estimate, and that the true function is hard to distinguish from samples drawn from our fit, indicating that we have reasonable approximation of the true function given noisy observations.

Using the prior knowledge of the noise in the data has given us a reasonably good approximation of the true function. The samples from the GP demonstrate provide fits to the observed data and roughly follow the shape of the true function $f$, particularly in the range for which we have samples.

Exercise 5¶

The function we used is a sum of sinusoids, and therefore has inherent periodicity. However, we only have samples for a single period, so this is not directly seen in the domain of $[0, 1]$.

(a) Try predicting the function on the range of $[0, 2]$ instead (Xnew = np.linspace(0., 2., 100)[:,None]), and compare the fitted GP posterior with the true function $f$.

Note, you shouldn't need to retrain the model or adjust the hyperparameters.

(b) Comment on the fit of the GP, and the uncertainty in regions where we have no observations. Is the GP still a good fit? How might we produce a better fit, for example, if we knew $f(x)$ had a periodic nature?

7. Real World Example¶

We'll consider now a classic real world example using atmospheric CO$_2$ observations from the Mauna Loa Observatory, Hawaii. The dataset is a good demonstration of the predictive power of Gaussian processes, and we will use it to show how we can encode our prior beliefs to combine kernels.

First, we need to download the dataset, which can be done through a utility provided by GPy.

You may want to subsample the data to save time during the labs

Run the following to reduce the datasize by only using every other training point:

Naive GP regression¶

First, we will try to fit a basic RBF to our data, as we have used in previous examples

It is possible to make the model fit the data near perfectly by minimising the variance of the Gaussian noise in the likelihood and fixing the kernel variance.

But this has no predictive power, and we have really just overfitted to the training data

GP Regression with combined covariance functions¶

Taking a look at the training data, we can see a number of features that occur in the data. There is a clear periodic trend that is yearly, and an approximately linear trend. We can use this prior information in our choice of kernel to give some meaning to the GP fit.

First, we will look at the linear trend. It should be obvious that the overall trend (ignoring the periodicity) can be described approximately by $f(x) \approx a + bx$. To embed this as a covariance function, we can use the Linear covariance functions, which adds a linear trend to the covariance.

Exercise 6¶

(a) Create a Linear kernel with reasonable estimates of the parameters that represent the trend?

(b) How might we encode this trend in the mean estimate of the GP ?

(c) Create a GP regression model using the kernels to fit the data. Comment on good is the fit and the predictive power of the model?

Periodicity¶

It's there is a seasonal trend over the year, and that a simple linear fit cannot capture this information. However, we can add this using our choice of kernel by adding a StdPeriodic kernel to our regression model. It's evident to the data that the period is yearly, so a period of $\omega=1$ is a sensible choice for our initial parameter:

From the plot we can see that while we have maintained the periodicity in our prediction, there is some deviation in the amplitude of each period. Likewise, the data is not strictly linear. We can embed model these non-linearities and deviations by adding an Exponential kernel. First, in the summative kernel, but also by multiplying a Gaussian RBF by the periodic kernel. For some intiution on how this allows for deviation of the periodicity, see the corresponding sample plot in Section 4.

Note that the Exponential kernel describes an exponential decay in covariance with distance between points, and does not model an exponential trend in the overall data. The RBF kernel describes a smoother (Gaussian) decay, which may be a better model for changes in the periodic signal.

We can see that with even with our initial (mostly default) parameter choices, the GP fit to the training data is representative of the trends we are observing. However, the prediction of the test data is relatively poor (though we still have an upward trend). The obvious step now is to optimise our kernel parameters.

Optimisation of the GP in itself is imperfect, often because the likelihood we are maximising can be multimodal, or flat, and so it can get stuck in local-maxima. It is always important to sanity-check the GP fit when optimising the parameters, to mitigate problems that could occur as a result (such as the minimised noise example we saw earlier). One of the ways to avoid the problem of local maxima is to reinitialise the optimiser with different starting locations, and take the maximum of these outputs. This is possible in GPy using optimize_restarts(n), which will optimise the parameters $n$ times and take the best estimate.

We can see from the output of the optimiser that there are a number of different minima in the negative log-likelihood (corresponding to maxima in the likelihood), implying that there are several modes. However, it's possible to get a good fit of the function with a GP with great predictive power. You might try allowing the optimiser to run for more iterations when you have more time.

We can also look at the effects of each kernel on the GP model, which will show the features that it represents in the fit. For example, we expect the Bias to capture the offset, and the Exponential kernel to extract most of the non-linearity (that's not explained by the periodicity). We can view the effects of each of the components and the combination of them to see how they deconstruct the training data. This shows what features are being learned by which kernel.

It's clear that while our final fit is very good, some of the kernel effects are not as clear cut as we expected. For example, much of the upward trend is captured by the periodic combination, likely due to its ability to vary non-linearly. However, we can see that, for example, the Exponential kernel has a clear purpose in adjusting the GP fit to the data.

Exercise 7¶

What improvements might be made to the GP model, for example in terms combinations of kernels and parameters, to get a better fit, especially for the prediction of our testing data?

8. Multiple Inputs¶

Typically, we will have data that is more than one-dimensional. Gaussian processes regression is scalable in dimension, and everything used above can be used with multiple inputs. The following problem is reasonably straight-forward, a simple extension of the regression problem we've used to 2-dimensional input.

Consider the toy example: $$ f = (x_1, x_2) \mapsto \sin(x_1)\sin(x_2) $$ $$ y = f(\mathbf{X}) + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2) $$

Now, our input space is made up of two coordinates, $(x_1, x_2)$. Suppose we have some random noisy observations of the function, which are 1-dimensional in nature:

Our observations then are simply our sample locations propagated through $f$ with some Gaussian noise

We can fit a Gaussian process over the observations by selecting a 2-D kernel. For example, the RBF kernel has a 2-D (really N-D) form, so we can use this as our kernel choice. Here we need to set both the input_dim to 2, and specify the dimensions of the data we want the kernel to act over. In this case, it is the first 2 (of 2) dimensions, which in Python are indexed by 0 and 1.

We can see here that the fit with the GP is generally good (this is quite a simple example, after all), and that functionally, the setup is almost identical to the 1-D case using GPy. We can also observe from the variance surface that we have the most uncertainty where there are no samples, which matches our expectations and the observation of the power of GPs in 1-D.

Notice that in the model preview, we only have a single lengthscale that is based on both dimensions. This is because by default, the RBF kernel is isotropic, which in our case is fine. However, in reality, it is often desirable to have different parameters corresponding to different dimensions. There are two approaches to doing this in GPy. The first is to enable different anisotropy within the kernel itself, by using the ARD keyword, e.g.

k = GPy.kern.RBF(2, ARD=True)

Alternatively, we can use kernel composition to create two kernels acting seperately on different dimensions. An alternative construction of the anisotropic RBF is to construct a kernel from two 1-D RBF kernels active on $x_1$ and $x_2$ respectively:

k = GPy.kern.RBF(1, active_dims=[0]) * GPy.kern.RBF(1, active_dims=[1])

We can also use this construction to allow for different covariance assumptions in different dimensions. In the following example, we use an RBF to model $x_1$ and a Matérn 3/2 kernel for $x_2$:

k = GPy.kern.RBF(1, active_dims=[0]) * GPy.kern.Matern32(1, active_dims=[1])

Footnote¶

Using Gaussian processes for machine learning gives you a powerful tool for learning latent functions in a Bayesian fashion, especially with a low number of observations. However, full GP regression can be quite time consuming, given it scales $\mathscr{O}(n^3)$ with the number of training points, $n$. There are methods for dealing with this: for example, sparse approximation or representing the regression problem as a state-space model. Each of these comes with their own assumptions and drawbacks. In tomorrow's session, there will be a lecture on computational efficiency in GPs and the next lab will have examples of some of these approaches.

When using GPs in a computational setting, there are a wide range of tools that can be used depending on your preference of programming language and libraries. GPy, the library we use in these labs, is a well established Python library that was developed here at the University of Sheffield, and is maintained by members here and at Amazon Research Cambridge. Some other libraries that are commonly used for Gaussian processes include:

Credit¶

This notebook was written by Wil Ward. It adapted from notebooks by Rich Wilkinson and Neil Lawrence. Additional writing and feedback was provided by Maria Skoularidou and Chunchao Ma.