Gaussian Processes: a versatile data science method that packs infinite dimensions

Introduction¶

One of the most basic problems in data science goes something like this:

I want to fit a curve to data that is flexible enough to capture non-linear relationships (i.e. go through known data points without any error) but I also want confidence intervals for predictions at unknown points.

Most methods give you one or the other. You can go with linear/polynomial regresssions and have confidence intervals, but they won’t perfectly fit your data. Or you can go with splines, LOESS and other non-parametric methods, but then you won’t get confidence intervals out of the box (you still can estimate them with bootstrapping but it’s not always straightforward).

It was in this context that I learned about Gaussian Processes (GPs) first. GPs give you both a perfect interpolating function over known data points while also providing variance estimates at unknown data points. Here’s a very basic (and recognizable) illustration of a GP in action.

Once I learned about Gaussian processes, I started seeing references to them everywhere. Want to optimize hyperparameters? GPs can help. Reading an academic paper doing some cool stuff? Oh look, GPs are involved, too. That’s how I decided to spend some more time to understand them (a bit better), and they did not disappoint.

Cool stuff is meant to be shared, and that is how this post/notebook was born. In it, you will find:

1. Brief intro into how GPs work
2. GP regression as a way to interpolate in multi-dimensional space
3. Importance of kernels in GPs
4. GPs for time-series forecasting
5. GP regressions with multiple output variables 🤯
6. GP regressions as surrogate functions and Bayesian optimization

Be forewarned: it’s a loooong one!

Data used in this notebook¶

One of the earliest applications of Gaussian Process is known as kriging, named after Danie G. Krige, a mining enginner from South Africa. Staying close to the roots, this notebook uses weather forecasts in Lithuania to illustrate GPs (and it makes for pretty charts, too!). An extract of the dataset is displayed below. This notebook uses GPy library. If you are reading it in post format, some of the code cells are excluded for brevity. You can find the entire notebook and associated data on GitHub.

How do GPs work?¶

Parametric algorithms typically fit a fixed number of parameters to data (think linear regression and its coefficients). Non-parametric methods, on the other hand, typically use data points themselves as “parameters” (think K-nearest neighbours where classification depends solely on existing data points).

Gaussian Processes fall somewhere in between. Instead of fitting a fixed number of parameters, they fit a distribution of functions to data, where these functions are drawn from a multivariate normal distribution (which is what enables having confidence intervals out of the box), while the functions themselves are additionally parameterized by existing data points. That’s not an easy to digest sentence, I know… One of the best definitions I found is from Julia’s GP package:

Gaussian processes are a family of stochastic processes which provide a flexible nonparametric tool for modelling data. A Gaussian Process places a prior over functions, and can be described as an infinite dimensional generalisation of a multivariate Normal distribution. Moreover, the joint distribution of any finite collection of points is a multivariate Normal. This process can be fully characterised by its mean and covariance functions, where the mean of any point in the process is described by the mean function and the covariance between any two observations is specified by the kernel. Given a set of observed real-valued points over a space, the Gaussian Process is used to make inference on the values at the remaining points in the space.

In some ways, GPs thus incorporate both parametric and non-parametric worlds – they have fixed parameters (the mean and the variance of the prior), but they also use data as parameters (as the posterior of the functions is effectively drawn from n-dimensional Gaussian distribution, where n is the number of sample points).

If you are interested to understand GPs deeper (and I definitely suggest that), I highly recommend a tutorial on YouTube by Richard Turner and a visual exploration of Gaussian processes by three researchers from University of Konstanz.

Finally, I cannot leave out two wonderful memes by Richard McElreath…

If you are wondering how one can fit a “infinite-dimensional” normal distribution, the answer lies in the fact that a marginal distribution of a normal distribution is normal, just like is a joint and a conditional distribution of two normal distributions. Richard Turner does a really great job at illustrating these dynamics in his tutorial. If you don’t watch it now, make sure at least to put it in your bookmarks/open tabs list – it’s a really good one!

For all their coolness, GPs Achilles heel is computational complexity. They scale with amount of data as any other non-parametric methods, with most implementations having $O(n^3)$. Improving that is an active line of research, whereas for every day data scientists concerned with speed can choose between GPU-enabled GPyTorch and efficiency focused celerite libraries. Just keep in mind that GPs are not really an option if you have 1,000,000s of observations.

1. Gaussian Process Regression – interpolations¶

Let’s start with the classical application of Gaussian Process in geostatistics. Imagine you are building an application that requires knowing air temperature at any given location in Lithuania. What you have, though, is just the measurements at weather stations. Here’s how it looks visually.

Let’s use Gaussian Process regression to interpolate. We will fit a GP using coordinates as features and air temperature as the dependent variable.

And, here it is – now we have predictions over the entire country. Looks pretty, doesn’t it?

But how good is it?¶

Pretty, sure, but are these predictions any good? Let’s test it out. I will split the dataset into a train and a test set and measure the prediction accuracy by mean absolute error (MAE).

Looks pretty good to me! It’s off by 0.1 degree (Celsius) on average in a dataset where temperature ranges from $-7.9$ to $-0.3$ degrees.

What is somewhat confusing is that its training error is not zero. GPs are interpolating functions, so this should only happen when there are multiple observations at the same input points (in which case the default GP kernel would end up predicting the mean of the observations), and there are no such cases in the dataset we are looking at. I suspect that this happens either because of numerical artifacts in the fitting process or there may be very sharp differences in temperature in nearby locations that the default GP kernel is not able to perfectly capture. Conceptually, that number should be zero, though.

So far, we looked only at mean predictions. Remember that one of the unique aspects of a Gaussian Process is that we also get confidence (variance) over the predictions. Let’s see how accurate the confidence intervals are.

It seems that in our case, the confidence intervals produced by Gaussian Process are actually even conservative!

2. Kernels in Gaussian Processes¶

You may have noticed that I have used all defaults of the GPy regresssion module, and may be curious if an even better performance may be achieved by tuning some of the hyperparameters. By far the most important hyperparameter in Gaussian Processes is the kernel function that defines the correlation between two data points (the default is the radial basis function, $R(h) = e^{-\theta h^2}$). GPy has a number of kernels implemented, and even allows defining custom kernels. They have an entire notebook dedicated to exploring kernels which is definitely worthwhile skimming through.

Let’s see what we can achieve in our interpolation task by using other kernels. Specifically, let’s try out:

• Matern52 kernel (Matern kernels are a generalization of RBF kernels with an additional smoothing parameter)
• A combination of Matern52 and a linear kernel. I choose this particular combination by reasoning that there may be a linear temperature trend across latitude/longitude that may be worth trying to capture.

And, indeed – the performance is better; using kernels allowed reducing test MAE by 20%. So, if there is one thing to remember – make sure to test different kernels (and see if you can reason what kernel combinations may be most appropriate given your data).

3. Forecasting with GPs¶

The flexibility of kernel functions is what makes GPs useful for time series forecasting, too. However, you won’t get far with a default RBF kernel as its predictions converge to zero (or mean, if a bias kernel is added) very quickly once you step out of a domain observed during training. Here’s an illustration using the 1D example from earlier.

I will illustrate a better way to do it using one of the city temperature forecasts. We will take the first 48 hours of predictions in one of the Lithuanian cities and see if we can generate a reasonable forecast for it.

Given the cyclical nature of temperature during a day, let’s use a periodic exponential kernel combined with a linear kernel that should capture the day-over-day trend.

Here’s the initial result.

The predictions appear quite reasonable – GP definitely picked up the cyclical trends and there seem to be a linear downwards trend, too. But there’s something fishy. We see that GPs (again!) do not perfectly interpolate existing data. Furthermore, what is with the variances? They seem to be constant over time which is weird, too. What is going on?

In this case, the answer is that the GP optimization process is stochastic, and a single attempt did not yield good results. If we rerun the optimization process multiple (say, 20) times, we will see different results (see below). This is another take away about GPs – in most cases, you may need to re-initialize the fitting process multiple times to get what you are looking for.

At the same time, looking below you may reasonably ask: is this not overfitting? I would tend to agree, actually! How do you go about reducing overfitting if you have that problem? Recall that GPs are defined over a prior of functions – it is possible to provide stronger priors as well as introduce constraints over some of the fitting parameters. I haven’t tried it out myself, but if you are interested, you can find more detail in one of the GPy tutorials.

GPs and multiple dependent variables 🤯¶

Another cool aspect of Gaussian Processes is that they can model relationships between multiple output variables. I have to admit that the maths behind how it works exactly is beyond me, though I suspect it must have to do with the fact that a joint distribution of multiple normal distributions is a normal distribution, so one can model the joint distribution as a single output and then marginalize it as needed.

Continuing with the weather data examples, here’s an application to consider: Suppose that, due to historical reasons, not all stations can measure temperature and wind speed simultaneously. The ones on the east side of the country tend to measure only temperature, and the ones to the west – only wind speed. If we wanted to build models to predict temperature and wind speed at other locations, we could simply build two models independently, using the relevant station data for each each. But we may have reason to believe that temperature and wind speed are related, and capturing that relationship in the model may be beneficial. Let’s see if that is the case.

Below is a visual illustration – I’ll be using a subset of the weather stations in this example.

To establish the baseline, let’s first train a GP regressions for the two datasets (temperature-measuring stations and wind-measuring stations) separately and see how well the models predict the other measurement at the locations where that measurement is not available to the model (i.e. how well a model trained on wind-measuring stations fares in predicting wind speed at temperature-measuring locations and vice-versa).

We can see that the erorrs are not too bad, actually. Can a multi-output GP regression do better?

The GPy model type to be used for such scenarios is GPCoregionalizedRegression, and it has a somewhat different API than most of the other GPy models (you need to pass inputs as list for each of the output variable; predictions require an extra index column in the input array). I based the code below on their dedicated tutorial for coregionalized regressions.

The results for this dataset are… interesting. In case of temperature, the model performs clearly better, but its performance is actually slightly worse for wind speeds. The thing about wind speeds, however, is that the forecasts are available only on an integer basis. So these results may simply be a factor of underlying data and shouldn’t necessarily be interpreted as a “worse fit”.

And, yes, I know – one really shouldn’t fit a continuous function if the output variables are count-like.. But hey, you know what – it is possible to model Poisson-like variables with GPs too! Maybe another time.

GPs as surrogate functions and Bayesian optimization¶

So far, I explored use cases where GP is the main model we use for fitting. The flexibility of GPs, however, make them very popular candidates to be used as surrogate functions, too.

Imagine you have a model that takes a long time to run – perhaps it is something as simple as an XGBoost classifier on a relatively large dataset, or maybe it is a complex physics simulation of the universe that takes a day to run on a supercomputer cluster. In both cases, you can only make a few runs of them with different (hyper-) parameters, and you would like to find the optimal setting of those hyperparameters (e.g. tune a learning rate for XGBoost) or use the limited outputs to interpolate in the rest of input space.

In such scenarios, you can do the following:

• Run the expensive model a few times with carefully spaced out parameter values and obtain initial observations of the dependent variable (or a loss function, in case of an ML model)
• Then, fit a Gaussian Process over the obtained data points and treat it as an approximate (surrogate) function of the underlying model
• Use GP predictions to infer expected values of the true underlying model at unknown locations
• If optimization is of interest (e.g. finding XGBoost hyperparameters that lead to the best model performance), use the uncertainty of predictions (variance) to determine the next combination of parameters that are most promising to evaluate the true underlying model on.
• Rinse and repeat the steps above iteratively.

The iterative process of looking for an optimal value is what is also known as “exploration vs exploitation”. You choose to either explore (using variance information to determine locations with highest uncertainty) or exploit (by evaluating the model at locations close to the best currently known point). There are multiple ways how to determine the next point for evaluation (i.e. how to make the exploration vs. exploitation trade-off). One of them is known as the Expected Improvement maximization approach. I won’t delve into detail on it here – I highly recommend a post by Martin Krasser on Bayesian Optimization if you are curious about the details.

In this notebook, I won’t fit XGBoost models. Instead, I will illustrate the concept with a simple task. Let’s pretend our task is to find the coldest place in Lithuania as measured by chill factor (which happens to be a village called Kirmėliai, with wind chill factor of -11.58). Here’s how the distribution of chill factor looks geographically; you can see it is not a convex environment, and so using gradient-based approaches is unlikely to get us to the right place.

Instead, let’s use Bayesian optimization with Gaussian processes. There are a few python packages that support Bayesian optimization, most notably hyperopt and pyGPGO. Hyperopt uses tree of Parzen estimators (TPE) as the surrogate function, however. For this illustration, I chose to go with the aptly named BayesianOptimization package that uses Gaussian processes and supports Expected Improvement algorithm for iterative choices.

Let’s start by randomly choosing 4 points in Lithuania and evaluating the chill factor there (we’ll use the nearest station measure for any given point). Then, let’s explore additional 15 locations iteratively, where each subsequent location is chosen by Expected Improvement criteria using a Gaussian Process model fit on all observations collected so far.

What are the results? Turns out we were able to find the coldest place!

To get the intuition behind what happened, here’s a graphical illustration of the algorithm’s path. You can see how trade-offs between exploring unknown areas vs. exploiting the information collected so far about the maximum observed value changes as more and more data points are collected. Much better than any grid search!