Gaussian process regression

Lecture 1

Suzanne Aigrain, University of Oxford

LSST DSFP Session 4, Seattle, Sept 2017

• Lecture 1: Introduction and basics

  • Tutorial 1: Write your own GP code

• Lecture 2: Examples and practical considerations

  • Tutorial 3: Useful GP modules

• Lecture 3: Advanced applications

Why GPs?

• flexible, robust probabilistic regression and classification tools.

• applied across a wide range of fields, from finance to zoology.

• useful for data containing non-trivial stochastic signals or noise.

• time-series data: causation implies correlation, so noise always correlated.

• increasingly popular in astronomy [mainly time-domain, but not just].

Spitzer exoplanet transits and eclipses (Evans et al. 2015)

GPz photometric redshifts (Almosallam, Jarvis & Roberts 2016)

What is a GP?

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution.

Consider a scalar variable $y$, drawn from a Gaussian distribution with mean $\mu$ and variance $\sigma^2$:

$$ p(y) = \frac{1}{\sqrt{2 \pi} \sigma} \exp \left[ - \frac{(y-\mu)^2}{2 \sigma^2} \right]. $$

As a short hand, we write: $y \sim \mathcal{N}(\mu,\sigma^2)$.

Now let us consider a pair of variables $y_1$ and $y_2$, drawn from a bivariate Gaussian distribution. The joint probability density for $y_1$ and $y_2$ is:

$$ \left[ \begin{array}{l} y_1 \\ y_2 \end{array} \right] \sim \mathcal{N} \left( \left[ \begin{array}{l} \mu_1 \\ \mu_2 \end{array} \right] , \left[ \begin{array}{ll} \sigma_1^2 & C \\ C & \sigma_2^2 \end{array} \right] \right), $$

where $C = {\rm cov}(y_1,y_2)$ is the covariance between $y_1$ and $y_2$. The second term on the right hand side is the covariance matrix, $K$.

We now use two powerful identities of Gaussian distributions to elucidate the relationship between $y_1$ and $y_2$.

The marginal distribution of $y_1$ describes what we know about $y_1$ in the absence of any other information about $y_2$, and is simply:

$$ p(y_1)= \mathcal{N} (\mu_1,\sigma_1^2). $$

If we know the value of $y_2$, the probability density for $y_1$ collapses to the the conditional distribution of $y_1$ given $y_2$:

$$ p(y_1 \mid y_2) = \mathcal{N} \left( \mu_1 + C (y_2-\mu_2)/\sigma_2^2, \sigma_1^2-C^2\sigma_2^2 \right). $$

If $K$ is diagonal, i.e. if $C=0$, $p(y_1 \mid y_2) = p(y_1)$. Measuring $y_2$ doesn't teach us anything about $y_1$. The two variables are uncorrelated.

If the variables are correlated ($C \neq 0$), measuring $y_2$ does alter our knowledge of $y_1$: it modifies the mean and reduces the variance.

To make the relation to time-series data a bit more obvious, let's plot the two variables side by side, then see what happens to one variable when we observe (fix) the other.

Now consider $N$ variables drawn from a multivariate Gaussian distribution:

$$ \boldsymbol{y} \sim \mathcal{N} (\boldsymbol{m},K) $$

where $y = (y_1,y_2,\ldots,y_N)^T$, $\boldsymbol{m} = (m_1,m_2,\ldots,m_N)^T$ is the mean vector, and $K$ is an $N \times N$ positive semi-definite covariance matrix, with elements $K_{ij}={\rm cov}(y_i,y_j)$.

A Gaussian process is an extension of this concept to infinite $N$, giving rise to a probability distribution over functions.

This last generalisation may not be obvious conceptually, but in practice only ever deal with finite samples.

Textbooks

A good, detailed reference is Gaussian Processes for Machine Learning by C. E. Rasmussen & C. Williams, MIT Press, 2006.

The examples in the book are generated using the Matlab package GPML.

A more formal definition

A Gaussian process is completely specified by its mean function and covariance function.

We define the mean function $m(x)$ and the covariance function $k(x,x)$ of a real process $y(x)$ as $$ \begin{array}{rcl} m(x) & = & \mathbb{E}[y(x)], \\ k(x,x') & = & \mathrm{cov}(y(x),y(x'))=\mathbb{E}[(y(x) − m(x))(y(x') − m(x'))]. \end{array} $$

A very common covariance function is the squared exponential, or radial basis function (RBF) kernel $$ K_{ij}=k(x_i,x_j)=A \exp\left[ - \Gamma (x_i-x_j)^2 \right], $$ which has 2 parameters: $A$ and $\Gamma$.

We then write the Gaussian process as $$ y(x) \sim \mathcal{GP}(m(x), k(x,x')) $$

Here we are implicitly assuming the inputs $x$ are one-dimensional, e.g. $x$ might represent time. However, the input space can have more than one dimension. We will see an example of a GP with multi-dimensional inputs later.

The prior

Now consider a finite set of inputs $\boldsymbol{x}$, with corresponding outputs $\boldsymbol{y}$.

The joint distribution of $\boldsymbol{y}$ given $\boldsymbol{x}$, $m$ and $k$ is $$ \mathrm{p}(\boldsymbol{y} \mid \boldsymbol{x},m,k) = \mathcal{N}( \boldsymbol{m},K), $$

where $\boldsymbol{m}=m(\boldsymbol{x})$ is the mean vector,

and $K$ is the covariance matrix, with elements $K_{ij} = k(x_i,x_j)$.

Test and training sets

Suppose we have an (observed) training set $(\boldsymbol{x},\boldsymbol{y})$.

We are interested in some other test set of inputs $\boldsymbol{x}_*$.

The joint distribution over the training and test sets is $$ \mathrm{p} \left( \left[ \begin{array}{l} \boldsymbol{y} \\ \boldsymbol{y}_* \end{array} \right] \right) = \mathcal{N} \left( \left[ \begin{array}{l} \boldsymbol{m} \\ \boldsymbol{m}_* \end{array} \right], \left[ \begin{array}{ll} K & K_* \\ K_*^T & K_{**} \end{array} \right] \right), $$

where $\boldsymbol{m}_* = m(\boldsymbol{x}_*)$, $K_{**,ij} = k(x_{*,i},x_{*,j})$ and $K_{*,ij} = k(x_i,x_{*,j})$.

For simplicity, assume the mean function is zero everywhere: $\boldsymbol{m}=\boldsymbol{0}$. We will consider to non-trivial mean functions later.

The conditional distribution

The conditional distribution for the test set given the training set is: $$ \mathrm{p} ( \boldsymbol{y}_* \mid \boldsymbol{y},k) = \mathcal{N} ( K_*^T K^{-1} \boldsymbol{y}, K_{**} - K_*^T K^{-1} K_* ). $$

This is also known as the predictive distribution, because it can be use to predict future (or past) observations.

More generally, it can be used for interpolating the observations to any desired set of inputs.

This is one of the most widespread applications of GPs in some fields (e.g. kriging in geology, economic forecasting, ...)

Adding white noise

Real observations always contain a component of white noise, which we need to account for, but don't necessarily want to include in the predictions.

If the white noise variance $\sigma^2$ is constant, we can write $$ \mathrm{cov}(y_i,y_j)=k(x_i,x_j)+\delta_{ij} \sigma^2, $$

and the conditional distribution becomes $$ \mathrm{p} ( \boldsymbol{y}_* \mid \boldsymbol{y},k) = \mathcal{N} ( K_*^T (K + \sigma^2 \mathbb{I})^{-1} \boldsymbol{y}, K_{**} - K_*^T (K + \sigma^2 \mathbb{I})^{-1} K_* ). $$

In real life, we may need to learn $\sigma$ from the data, alongside the other contribution to the covariance matrix.

We assumed constant white noise, but it's trivial to allow for different $\sigma$ for each data point.

Single-point prediction

Let us look more closely at the predictive distribution for a single test point $x_*$.

It is a Gaussian with mean $$ \overline{y}_* = \boldsymbol{k}_*^T (K + \sigma^2 \mathbb{I})^{-1} \boldsymbol{y} $$ and variance $$ \mathbb{V}[y_*] = k(x_*,x_*) - \boldsymbol{k}_*^T (K + \sigma^2 \mathbb{I})^{-1} \boldsymbol{k}_*, $$ where $\boldsymbol{k}_*$ is the vector of covariances between the test point and the training points.

Notice the mean is a linear combination of the observations: the GP is a linear predictor.

It is also a linear combination of covariance functions, each centred on a training point: $$ \overline{y}_* = \sum_{i=1}^N \alpha_i k(x_i,x_*), $$ where $\alpha_i = (K + \sigma^2 \mathbb{I})^{-1} y_i$.

The likelihood

The likelihood of the data under the GP model is simply: $$ \mathrm{p}(\boldsymbol{y} \,|\, \boldsymbol{x}) = \mathcal{N}(\boldsymbol{y} \, | \, \boldsymbol{0},K + \sigma^2 \mathbb{I}). $$

This is a measure of how well the model explains, or predicts, the training set.

In some textbooks this is referred to as the marginal likelihood.

This arises if one considers the observed $\boldsymbol{y}$ as noisy realisations of a latent (unobserved) Gaussian process $\boldsymbol{f}$.

The term marginal refers to marginalisation over the function values $\boldsymbol{f}$: $$ \mathrm{p}(\boldsymbol{y} \,|\, \boldsymbol{x}) = \int \mathrm{p}(\boldsymbol{y} \,|\, \boldsymbol{f},\boldsymbol{x}) \, \mathrm{p}(\boldsymbol{f} \,|\, \boldsymbol{x}) \, \mathrm{d}\boldsymbol{f}, $$ where $$ \mathrm{p}(\boldsymbol{f} \,|\, \boldsymbol{x}) = \mathcal{N}(\boldsymbol{f} \, | \, \boldsymbol{0},K) $$ is the prior, and $$ \mathrm{p}(\boldsymbol{y} \,|\, \boldsymbol{f},\boldsymbol{x}) = \mathcal{N}(\boldsymbol{y} \, | \, \boldsymbol{0},\sigma^2 \mathbb{I}) $$ is the likelihood.

Parameters and hyper-parameters

The parameters of the covariance and mean function as known as the hyper-parameters of the GP.

This is because the actual parameters of the model are the function values, $\boldsymbol{f}$, but we never explicitly deal with them: they are always marginalised over.

Conditioning the GP...

...means evaluating the conditional (or predictive) distribution for a given covariance matrix (i.e. covariance function and hyper-parameters), and training set.

Training the GP...

...means maximising the likelihood of the model with respect to the hyper-parameters.

The kernel trick

Consider a linear basis model with arbitrarily many basis functions, or features, $\Phi(x)$, and a (Gaussian) prior $\Sigma_{\mathrm{p}}$ over the basis function weights.

One ends up with exactly the same expressions for the predictive distribution and the likelihood so long as: $$ k(\boldsymbol{x},\boldsymbol{x'}) = \Phi(\boldsymbol{x})^{\mathrm{T}} \Sigma_{\mathrm{p}} \Phi(\boldsymbol{x'}), $$ or, writing $\Psi(\boldsymbol{x}) = \Sigma_{\mathrm{p}}^{1/2} \Phi(\boldsymbol{x})$, $$ k(\boldsymbol{x},\boldsymbol{x'}) = \Psi(\boldsymbol{x}) \cdot \Psi(\boldsymbol{x'}), $$

Thus the covariance function $k$ enables us to go from a (finite) input space to a (potentially infinite) feature space. This is known as the kernel trick and the covariance function is often referred to as the kernel.

Non-zero mean functions

In general (and in astronomy applications in particular) we often want to use non-trivial mean functions.

To do this simply replace $\boldsymbol{y}$ by $\boldsymbol{r}=\boldsymbol{y}-\boldsymbol{m}$ in the expressions for predictive distribution and likelihood.

The mean function represents the deterministic component of the model

• e.g.: a linear trend, a Keplerian orbit, a planetary transit, ...

The covariance function encodes the stochastic component.

• e.g.: instrumental noise, stellar variability

Covariance functions

The only requirement for the covariance function is that it should return a positive semi-definite covariance matrix.

The simplest covariance functions have two parameters: one input and one output variance (or scale). The form of the covariance function controls the degree of smoothness.

The squared exponential

The simplest, most widely used kernel is the squared exponential: $$ k_{\rm SE}(x,x') = A \exp \left[ - \Gamma (x-x')^2 \right]. $$ This gives rise to smooth functions with variance $A$ and inverse scale (characteristic length scale) $A$ and output scale (amplitude) $l$.

The Matern family

The Matern 3/2 kernel $$ k_{3/2}(x,x')= A \left( 1 + \frac{\sqrt{3}r}{l} \right) \exp \left( - \frac{\sqrt{3}r}{l} \right), $$ where $r =|x-x'|$.

It produces somewhat rougher behaviour, because it is only differentiable once w.r.t. $r$ (whereas the SE kernel is infinitely differentiable). There is a whole family of Matern kernels with varying degrees of roughness.

The rational quadratic kernel

is equivalent to a squared exponential with a powerlaw distribution of input scales $$ k_{\rm RQ}(x,x') = A^2 \left(1 + \frac{r^2}{2 \alpha l} \right)^{-\alpha}, $$ where $\alpha$ is the index of the power law.

This is useful to model data containing variations on a range of timescales with just one extra parameter.

Periodic kernels...

...can be constructed by replacing $r$ in any of the above by a periodic function of $r$. For example, the cosine kernel: $$ k_{\cos}(x,x') = A \cos\left(\frac{2\pi r}{P}\right), $$ [which follows the dynamics of a simple harmonic oscillator], or...

...the "exponential sine squared" kernel, obtained by mapping the 1-D variable $x$ to the 2-D variable $\mathbf{u}(x)=(\cos(x),\sin(x))$, and then applying a squared exponential in $\boldsymbol{u}$-space: $$ k_{\sin^2 {\rm SE}}(x,x') = A \exp \left[ -\Gamma \sin^2\left(\frac{\pi r}{P}\right) \right], $$ which allows for non-harmonic functions.

Combining kernels

Any affine tranform, sum or product of valid kernels is a valid kernel.

For example, a quasi-periodic kernel can be constructed by multiplying a periodic kernel with a non-periodic one. The following is frequently used to model stellar light curves: $$ k_{\mathrm{QP}}(x,x') = A \exp \left[ -\Gamma_1 \sin^2\left(\frac{\pi r}{P}\right) -\Gamma_2 r^2 \right]. $$

Example: Mauna Kea CO$_2$ dataset

(From Rasmussen & Williams textbook)

2 or more dimensions

So far we assumed the inputs were 1-D but that doesn't have to be the case. For example, the SE kernel can be extended to D dimensions...

using a single length scale, giving the Radial Basis Function (RBF) kernel: $$ k_{\rm RBF}(\mathbf{x},\mathbf{x'}) = A \exp \left[ - \Gamma \sum_{j=1}^{D}(x_j-x'_j)^2 \right], $$ where $\mathbf{x}=(x_1,x_2,\ldots, x_j,\ldots,x_D)^{\mathrm{T}}$ represents a single, multi-dimensional input.

or using separate length scales for each dimension, giving the Automatic Relevance Determination (ARD) kernel: $$ k_{\rm ARD}(\mathbf{x},\mathbf{x'}) = A \exp \left[ - \sum_{j=1}^{D} \Gamma_j (x_j-x'_j)^2 \right]. $$