Gaussian Processes for Regression:
A Quick Introduction

M. Ebden, August 2008
Comments to mebden@gmail.com

 

1 Motivation

Figure 1 illustrates a typical example of a prediction problem: given some noisy observations of a dependent variable at certain values of the independent variable $x$, what is our best estimate of the dependent variable at a new value, $x_{\ast }$?

If we expect the underlying function $f(x)$ to be linear, and can make some assumptions about the input data, we might use a least-squares method to fit a straight line (linear regression). Moreover, if we suspect $f(x)$ may also be quadratic, cubic, or even nonpolynomial, we can use the principles of model selection to choose among the various possibilities.

Gaussian process regression (GPR) is an even finer approach than this. Rather than claiming $f(x)$ relates to some specific models (e.g. $f(x)=mx+c$), a Gaussian process can represent $f(x)$ obliquely, but rigorously, by letting the data ‘speak’ more clearly for themselves. GPR is still a form of supervised learning, but the training data are harnessed in a subtler way.

As such, GPR is a less ‘parametric’ tool. However, it’s not completely free-form, and if we’re unwilling to make even basic assumptions about $f(x)$, then more general techniques should be considered, including those underpinned by the principle of maximum entropy; Chapter 6 of Sivia and Skilling (2006) offers an introduction.

Figure 1: Given six noisy data points (error bars are indicated with vertical lines), we are interested in estimating a seventh at $x_{\ast }=0.2$.

2 Definition of a Gaussian Process

Gaussian processes (GPs) extend multivariate Gaussian distributions to infinite dimensionality. Formally, a Gaussian process generates data located throughout some domain such that any finite subset of the range follows a multivariate Gaussian distribution. Now, the $n$ observations in an arbitrary data set, $𝐲=\{y_1,\mathrm{\dots },y_n\}$, can always be imagined as a single point sampled from some multivariate ($n$-variate) Gaussian distribution, after enough thought. Hence, working backwards, this data set can be partnered with a GP. Thus GPs are as universal as they are simple.

Very often, it’s assumed that the mean of this partner GP is zero everywhere. What relates one observation to another in such cases is just the covariance function, $k(x,x^{\prime })$. A popular choice is the ‘squared exponential’,

$$k(x,x^{\prime })={\sigma }_f^2\exp \left[\frac{-{(x-x^{\prime })}^2}{2l^2}\right],$$ (1)

where the maximum allowable covariance is defined as ${\sigma }_f^2$ — this should be high for functions which cover a broad range on the $y$ axis. If $x\approx x^{\prime }$, then $k(x,x^{\prime })$ approaches this maximum, meaning $f(x)$ is nearly perfectly correlated with $f(x^{\prime })$. This is good: for our function to look smooth, neighbours must be alike. Now if $x$ is distant from $x^{\prime }$, we have instead $k(x,x^{\prime })\approx 0$, i.e. the two points cannot ‘see’ each other. So, for example, during interpolation at new $x$ values, distant observations will have negligible effect. How much effect this separation has will depend on the length parameter, $l$, so there is much flexibility built into (1).

Not quite enough flexibility though: the data are often noisy as well, from measurement errors and so on. Each observation $y$ can be thought of as related to an underlying function $f(x)$ through a Gaussian noise model:

$$y=f(x)+𝒩(0,{\sigma }_n^2),$$ (2)

something which should look familiar to those who’ve done regression before. Regression is the search for $f(x)$. Purely for simplicity of exposition in the next page, we take the novel approach of folding the noise into $k(x,x^{\prime })$, by writing

$$k(x,x^{\prime })={\sigma }_f^2\exp \left[\frac{-{(x-x^{\prime })}^2}{2l^2}\right]+{\sigma }_n^2\delta (x,x^{\prime }),$$ (3)

where $\delta (x,x^{\prime })$ is the Kronecker delta function. (When most people use Gaussian processes, they keep ${\sigma }_n$ separate from $k(x,x^{\prime })$. However, our redefinition of $k(x,x^{\prime })$ is equally suitable for working with problems of the sort posed in Figure 1. So, given $n$ observations $𝐲$, our objective is to predict $y_{\ast }$, not the ‘actual’ $f_{\ast }$; their expected values are identical according to (2), but their variances differ owing to the observational noise process. e.g. in Figure 1, the expected value of $y_{\ast }$, and of $f_{\ast }$, is the dot at $x_{\ast }$.)

To prepare for GPR, we calculate the covariance function, (3), among all possible combinations of these points, summarizing our findings in three matrices:

(4)

(5)

Confirm for yourself that the diagonal elements of $K$ are ${\sigma }_f^2+{\sigma }_n^2$, and that its extreme off-diagonal elements tend to zero when $x$ spans a large enough domain.

3 How to Regress using Gaussian Processes

Since the key assumption in GP modelling is that our data can be represented as a sample from a multivariate Gaussian distribution, we have that

(6)

where $\mathrm{T}$ indicates matrix transposition. We are of course interested in the conditional probability $p(y_{\ast }|𝐲)$: “given the data, how likely is a certain prediction for $y_{\ast }$?”. As explained more slowly in the Appendix, the probability follows a Gaussian distribution:

$$y_{\ast }|𝐲\sim 𝒩(K_{\ast }{K}^{-1}𝐲,K_{\ast \ast }-K_{\ast }{K}^{-1}{K}_{\ast }^{\mathrm{T}}).$$ (7)

Our best estimate for $y_{\ast }$ is the mean of this distribution:

$${\overline{y}}_{\ast }=K_{\ast }{K}^{-1}𝐲,$$ (8)

and the uncertainty in our estimate is captured in its variance:

$$\mathrm{var}(y_{\ast })=K_{\ast \ast }-K_{\ast }{K}^{-1}{K}_{\ast }^{\mathrm{T}}.$$ (9)

We’re now ready to tackle the data in Figure 1.

1. There are $n=6$ observations $𝐲$, at

   We know ${\sigma }_n=0.3$ from the error bars. With judicious choices of ${\sigma }_f$ and $l$ (more on this later), we have enough to calculate a covariance matrix using (4):

   From (5) we also have $K_{\ast \ast }=1.70$ and

2. From (8) and (9), ${\overline{y}}_{\ast }=0.95$ and $\mathrm{var}(y_{\ast })=0.21$.

3. Figure 1 shows a data point with a question mark underneath, representing the estimation of the dependent variable at $x_{\ast }=0.2$.

We can repeat the above procedure for various other points spread over some portion of the $x$ axis, as shown in Figure 2. (In fact, equivalently, we could avoid the repetition by performing the above procedure once with suitably larger $K_{\ast }$ and $K_{\ast \ast }$ matrices. In this case, since there are 1,000 test points spread over the $x$ axis, $K_{\ast \ast }$ would be of size 1,000 $\times$ 1,000.) Rather than plotting simple error bars, we’ve decided to plot ${\overline{y}}_{\ast }\pm 1.96\sqrt{\mathrm{var}(y_{\ast })}$, giving a 95% confidence interval.

Figure 2: The solid line indicates an estimation of $y_{\ast }$ for 1,000 values of $x_{\ast }$. Pointwise 95% confidence intervals are shaded.

4 GPR in the Real World

The reliability of our regression is dependent on how well we select the covariance function. Clearly if its parameters — call them $𝜽=\{l,{\sigma }_f,{\sigma }_n\}$ — are not chosen sensibly, the result is nonsense. Our maximum a posteriori estimate of $𝜽$ occurs when $p(𝜽|𝐱,𝐲)$ is at its greatest. Bayes’ theorem tells us that, assuming we have little prior knowledge about what $𝜽$ should be, this corresponds to maximizing $\log p(𝐲|𝐱,𝜽)$, given by

$$\log p(𝐲|𝐱,𝜽)=-\frac{1}{2}{𝐲}^{\mathrm{T}}{K}^{-1}𝐲-\frac{1}{2}\log |K|-\frac{n}{2}\log 2\pi .$$ (10)

Simply run your favourite multivariate optimization algorithm (e.g. conjugate gradients, Nelder-Mead simplex, etc.) on this equation and you’ve found a pretty good choice for $𝜽$; in our example, $l=1$ and ${\sigma }_f=1.27$.

It’s only “pretty good” because, of course, Thomas Bayes is rolling in his grave. Why commend just one answer for $𝜽$, when you can integrate everything over the many different possible choices for $𝜽$? Chapter 5 of Rasmussen and Williams (2006) presents the equations necessary in this case.

Finally, if you feel you’ve grasped the toy problem in Figure 2, the next two examples handle more complicated cases. Figure 3(a), in addition to a long-term downward trend, has some fluctuations, so we might use a more sophisticated covariance function:

$$k(x,x^{\prime })=\sigma _f^{}{}_1{}^2\exp \left[\frac{-{(x-x^{\prime })}^2}{2l_1^2}\right]+\sigma _f^{}{}_2{}^2\exp \left[\frac{-{(x-x^{\prime })}^2}{2l_2^2}\right]+{\sigma }_n^2\delta (x,x^{\prime }).$$ (11)

The first term takes into account the small vicissitudes of the dependent variable, and the second term has a longer length parameter ($l_2\approx 6l_1$) to represent its long-term trend. Covariance functions can be grown in this way ad infinitum, to suit the complexity of your particular data.

Figure 3:  Estimation of $y_{\ast }$ (solid line) for a function with (a) short-term and long-term dynamics, and (b) long-term dynamics and a periodic element. Observations are shown as crosses.

The function looks as if it might contain a periodic element, but it’s difficult to be sure. Let’s consider another function, which we’re told has a periodic element. The solid line in Figure 3(b) was regressed with the following covariance function:

$$k(x,x^{\prime })={\sigma }_f^2\exp \left[\frac{-{(x-x^{\prime })}^2}{2l^2}\right]+\exp \{-2{\sin }^2[\nu \pi (x-x^{\prime })]\}+{\sigma }_n^2\delta (x,x^{\prime }).$$ (12)

The first term represents the hill-like trend over the long term, and the second term gives periodicity with frequency $\nu$. This is the first time we’ve encountered a case where $x$ and $x^{\prime }$ can be distant and yet still ‘see’ each other (that is, $k(x,x^{\prime })\not\approx 0$ for $x\gg x^{\prime }$).

What if the dependent variable has other dynamics which, a priori, you expect to appear? There’s no limit to how complicated $k(x,x^{\prime })$ can be, provided $K$ is positive definite. Chapter 4 of Rasmussen and Williams (2006) offers a good outline of the range of covariance functions you should keep in your toolkit.

“Hang on a minute,” you ask, “isn’t choosing a covariance function from a toolkit a lot like choosing a model type, such as linear versus cubic — which we discussed at the outset?” Well, there are indeed similarities. In fact, there is no way to perform regression without imposing at least a modicum of structure on the data set; such is the nature of generative modelling. However, it’s worth repeating that Gaussian processes do allow the data to speak very clearly. For example, there exists excellent theoretical justification for the use of (1) in many settings (Rasmussen and Williams (2006), Section 4.3). You will still want to investigate carefully which covariance functions are appropriate for your data set. Essentially, choosing among alternative functions is a way of reflecting various forms of prior knowledge about the physical process under investigation.

5 Discussion

We’ve presented a brief outline of the mathematics of GPR, but practical implementation of the above ideas requires the solution of a few algorithmic hurdles as opposed to those of data analysis. If you aren’t a good computer programmer, then the code for Figures 1 and 2 is at github.com/mebden/GPtutorial, and more general code can be found at www.gaussianprocess.org/gpml.

We’ve merely scratched the surface of a powerful technique (MacKay, 1998). First, although the focus has been on one-dimensional inputs, it’s simple to accept those of higher dimension. Whereas $x$ would then change from a scalar to a vector, $k(x,x^{\prime })$ would remain a scalar and so the maths overall would be virtually unchanged. Second, the zero vector representing the mean of the multivariate Gaussian distribution in (6) can be replaced with functions of $x$. Third, in addition to their use in regression, GPs are applicable to integration, global optimization, mixture-of-experts models, unsupervised learning models, and more — see Chapter 9 of Rasmussen and Williams (2006). The next tutorial will focus on their use in classification.

Appendix

Imagine a data sample $𝐝$ taken from some multivariate Gaussian distribution with zero mean and a covariance given by matrix $D$. Now decompose $𝐝$ arbitrarily into two consecutive subvectors $𝐚$ and $𝐛$ — in other words, writing $𝐝\sim 𝒩(\mathrm{𝟎},D)$ would be the same as writing

(13)

where $A$, $B$, and $C$ are the corresponding bits and pieces that make up $D$.

Interestingly, the conditional distribution of $𝐛$ given $𝐚$ is itself Gaussian-distributed. If the covariance matrix $D$ were diagonal or even block diagonal, then knowing $𝐚$ wouldn’t tell us anything about $𝐛$: specifically, $𝐛|𝐚\sim 𝒩(\mathrm{𝟎},B)$. On the other hand, if $C$ were nonzero, then some matrix algebra leads us to

$$𝐛|𝐚\sim 𝒩(C{A}^{-1}𝐚,B-C{A}^{-1}{C}^{\mathrm{T}}).$$ (14)

The mean, $C{A}^{-1}𝐚$, is known as the ‘matrix of regression coefficients’, and the variance, $B-C{A}^{-1}{C}^{\mathrm{T}}$, is the ‘Schur complement of $A$ in $D$’.

In summary, if we know some of $𝐝$, we can use that to inform our estimate of what the rest of $𝐝$ might be, thanks to the revealing off-diagonal elements of $D$.

 

Gaussian Processes for Classification:
A Quick Introduction

M. Ebden, August 2008
Prerequisite reading: Gaussian Processes for Regression

 

1 Overview

As mentioned in the previous document, GPs can be applied to problems other than regression. For example, if the output of a GP is squashed onto the range $[0,1]$, it can represent the probability of a data point belonging to one of say two types, and voilà, we can ascertain classifications. This is the subject of the current document.

The big difference between GPR and GPC is how the output data, $𝐲$, are linked to the underlying function outputs, $𝐟$. They are no longer connected simply via a noise process as in (2) in the previous document, but are instead now discrete: say $y=1$ precisely for one class and $y=-1$ for the other. In principle, we could try fitting a GP that produces an output of approximately $1$ for some values of $x$ and approximately $-1$ for others, simulating this discretization. Instead, we interpose the GP between the data and a squashing function; then, classification of a new data point $x_{\ast }$ involves two steps instead of one:

1. Evaluate a ‘latent function’ $f$ which models qualitatively how the likelihood of one class versus the other changes over the $x$ axis. This is the GP.

2. Squash the output of this latent function onto $[0,1]$ using any sigmoidal function, $\pi (f)=\text{prob}(y=1|f)$.

Writing these two steps schematically,

data, $x_{\ast }$ $\stackrel{\text{GP}}{\to }$ latent function, $f_{\ast }|x_{\ast }$ $\stackrel{\text{sigmoid}}{\to }$ class probability, $\pi (f_{\ast })$.

The next section will walk you through more slowly how such a classifier operates. Section 3 explains how to train the classifier, so perhaps we’re presenting things in reverse order! Section 4 handles classification when there are more than two classes.

Before we get started, a quick note on $\pi (f)$. Although other forms will do, here we will prescribe it to be the cumulative Gaussian distribution, $\mathrm{\Phi }(f)$. This $S$-shaped function satisfies our needs, mapping high $f$ into $\pi (f)\approx 1$, and low $f$ into $\pi (f)\approx 0$.

A second quick note, revisiting (6) and (7) in the first document: confirm for yourself that, if there were no noise (${\sigma }_n=0$), the two equations could be rewritten as

(1)

and

$$f_{\ast }|𝐟\sim 𝒩(K_{\ast }{K}^{-1}𝐟,K_{\ast \ast }-K_{\ast }{K}^{-1}{K}_{\ast }^{\mathrm{T}}).$$ (2)

2 Using the Classifier

Suppose we’ve trained a classifier from $n$ input data, $𝐱$, and their corresponding expert-labelled output data, $𝐲$. And suppose that in the process we formed some GP outputs $𝐟$ corresponding to these data, which have some uncertainty but mean values given by $\widehat{𝐟}$. We’re now ready to input a new data point, $x_{\ast }$, in the left side of our schematic, in order to determine at the other end the probability ${\pi }_{\ast }$ of its class membership.

In the first step, finding the probability $p(f_{\ast }|𝐟)$ is similar to GPR, i.e. we adapt (2):

$$p(f_{\ast }|𝐟)=𝒩(K_{\ast }{K}^{-1}\widehat{𝐟},K_{\ast \ast }-K_{\ast }{(K^{\prime })}^{-1}{K}_{\ast }^{\mathrm{T}}).$$ (3)

($K^{\prime }$ will be explained soon, but for now consider it to be very similar to $K$.) In the second step, we squash $f_{\ast }$ to find the probability of class membership, ${\pi }_{\ast }=\pi (f_{\ast })=\mathrm{\Phi }(f_{\ast })$. The expected value is

$${\overline{\pi }}_{\ast }=\int \pi (f_{\ast })p(f_{\ast }|𝐟)𝑑f_{\ast }.$$ (4)

This is the integral of a cumulative Gaussian times a Gaussian, which can be solved analytically. By Section 3.9 of Rasmussen and Williams (2006), the solution is:

$${\overline{\pi }}_{\ast }=\mathrm{\Phi }\left(\frac{{\overline{f}}_{\ast }}{\sqrt{1+\mathrm{var}(f_{\ast })}}\right).$$ (5)

An example is depicted in Figure 1.

Figure 1: (a) Toy classification dataset, where circles and crosses indicate class membership of the input (training) data $𝐲$, at locations $𝐱$. $𝐲$ is $1$ for one class and $-1$ for another, but for illustrative purposes we pretend $y=0$ instead of $-1$ in this figure. The solid line is the (mean) probability ${\overline{\pi }}_{\ast }=\text{prob}(y_{\ast }=1|x_{\ast })$, i.e. the ‘answer’ to our problem after successfully performing GPC. (b) The corresponding distribution of the latent function $f$, not constrained to lie between 0 and 1.

3 Training the GP in the Classifier

Our objective now is to find $\widehat{𝐟}$ and $K^{\prime }$, so that we know everything about the GP producing (3), the first step of the classifier. The second step of the classifier does not require training as it’s a fixed sigmoidal function.

Among the many GPs which could be partnered with our data set, naturally we’d like to compare their usefulness quantitatively. Considering the outputs $𝐟$ of a certain GP, how likely they are to be appropriate for the training data can be decomposed using Bayes’ theorem:

$$p(𝐟|𝐱,𝐲)=\frac{p(𝐲|𝐟)p(𝐟|𝐱)}{p(𝐲|𝐱)}.$$ (6)

Let’s focus on the two factors in the numerator. Assuming the data set is i.i.d.,

$$p(𝐲|𝐟)=\prod _{i=1}^{n}p(y_i|f_i).$$ (7)

Dropping the subscripts in the product, $p(y|f)$ is informed by our sigmoid function, $\pi (f)$. Specifically, $p(y=1|f)$ is $\pi (f)$ by definition, and to complete the picture, $p(y=-1|f)=1-\pi (f)$. A terse way of combining these two cases is to write $p(y|f)=\mathrm{\Phi }(yf)$.

The second factor in the numerator is $p(𝐟|𝐱)$. This is related to the output of the first step of our schematic drawing, but first we’re interested in the value of $p(𝐟|𝐱)$ which maximizes the posterior probability $p(𝐟|𝐱,𝐲)$. This occurs when the derivative of (6) with respect to $𝐟$ is zero, or equivalently and more simply, when the derivative of its logarithm is zero. Doing this, and using the same logic that produced (10) in the previous document, we find that

$$\widehat{𝐟}=K\nabla \log p(𝐲|\widehat{𝐟}),$$ (8)

where $\widehat{𝐟}$ is the best $𝐟$ for our problem. Unfortunately, $\widehat{𝐟}$ appears on both sides of the equation, so we make an initial guess (zero is fine) and go through a few iterations. The answer to (8) can be used directly in (3), so we’ve found one of the two quantities we seek therein.

The variance of $𝐟$ is given by the negative second derivative of the logarithm of (6), which turns out to be ${(K^{-1}+W)}^{-1}$, with $W=-\nabla \nabla \log p(𝐲|𝐟)$. Making a Laplace approximation, we pretend $p(𝐟|𝐱,𝐲)$ is Gaussian distributed, i.e.

$$p(𝐟|𝐱,𝐲)\sim q(𝐟|𝐱,𝐲)=𝒩(\widehat{𝐟},{(K^{-1}+W)}^{-1}).$$ (9)

(This assumption is occasionally inaccurate, so if it yields poor classifications, better ways of characterizing the uncertainty in $𝐟$ should be considered, for example via expectation propagation.)

Now for a subtle point. The fact that $𝐟$ can vary means that using (2) directly is inappropriate: in particular, its mean is correct but its variance no longer tells the whole story. This is why we use the adapted version, (3), with $K^{\prime }$ instead of $K$. Since the varying quantity in (2), $𝐟$, is being multiplied by $K_{\ast }{K}^{-1}$, we add $K_{\ast }{K}^{-1}\text{cov}(𝐟)K^{-1}{K}_{\ast }^{\mathrm{T}}$ to the variance in (2). Simplification leads to (3), in which $K^{\prime }=K+W^{-1}$.

With the GP now completely specified, we’re ready to use the classifier as described in the previous section.

GPC in the Real World
As with GPR, the reliability of our classification is dependent on how well we select the covariance function in the GP in our first step. The parameters are $𝜽=\{l,{\sigma }_f\}$, one fewer now because ${\sigma }_n=0$. However, as usual, $𝜽$ is optimized by maximizing $p(𝐲|𝐱,𝜽)$, or (omitting $𝜽$ on the righthand side of the equation),

$$p(𝐲|𝐱,𝜽)=\int p(𝐲|𝐟)p(𝐟|𝐱)𝑑𝐟.$$ (10)

This can be simplified, using a Laplace approximation, to yield

$$p(𝐲|𝐱,𝜽)=-\frac{1}{2}{\widehat{𝐟}}^{\mathrm{T}}{K}^{-1}\widehat{𝐟}+\log p(𝐲|\widehat{𝐟})-\frac{1}{2}\log (|K|\cdot |K^{-1}+W|).$$ (11)

This is the equation to run your favourite optimizer on, as performed in GPR.

4 Multi-Class GPC

We’ve described binary classification, where the number of possible classes, $C$, is just two. In the case of $C>2$ classes, one approach is to fit an $f$ for each class. In the first of the two steps of classification, our GP values are concatenated as

$$𝐟={(f_1^1,\mathrm{\dots },f_n^1,f_1^2,\mathrm{\dots },f_n^2,\mathrm{\dots },f_1^C,\mathrm{\dots },f_n^C)}^{\mathrm{T}}.$$ (12)

Let $𝐲$ be a vector of the same length as $𝐟$ which, for each $i=1,\mathrm{\dots },n$, is $1$ for the class which is the label and $0$ for the other $C-1$ entries. Let $K$ grow to being block diagonal in the matrices $K^1,\mathrm{\dots },K^C$. So the first change we see for $C>2$ is a lengthening of the GP. Section 3.5 of Rasmussen and Williams (2006) offers hints on keeping the computations manageable.

The second change is that the (merely one-dimensional) cumulative Gaussian distribution is no longer sufficient to describe the squashing function in our classifier; instead we use the softmax function. For the $i$th data point,

$$p(y_i^c|{𝐟}_i)={\pi }_i^c=\frac{\exp (f_i^c)}{{\sum }_{c^{\prime }}\exp (f_i^{c^{\prime }})}$$ (13)

where ${𝐟}_i$ is a nonconsecutive subset of $𝐟$, viz. ${𝐟}_i=\{f_i^1,f_i^2,\mathrm{\dots },f_i^C\}$. We can summarize our results with $𝝅=\{{\pi }_1^1,\mathrm{\dots },{\pi }_n^1,{\pi }_1^2,\mathrm{\dots },{\pi }_n^2,\mathrm{\dots },{\pi }_1^C,\mathrm{\dots },{\pi }_n^C\}$.

Now that we’ve presented the two big changes needed to go from binary- to multi-class GPC, we continue as before. Setting to zero the derivative of the logarithm of the components in (6), we replace (8) with

$$\widehat{𝐟}=K(𝐲-\widehat{𝝅}).$$ (14)

The corresponding variance is ${(K^{-1}+W)}^{-1}$ as before, but now $W=\text{diag}(𝝅)-\mathrm{\Pi }{\mathrm{\Pi }}^{\mathrm{T}}$, where $\mathrm{\Pi }$ is a $Cn\times n$ matrix obtained by stacking vertically the diagonal matrices $\text{diag}({𝝅}^c)$, if ${𝝅}^c$ is the subvector of $𝝅$ pertaining to class $c$.

With these quantities estimated, we have enough to generalize (3) to

$$p(f_{\ast }^c|𝐟)=𝒩(K_{\ast }^c{(K^c)}^{-1}{\widehat{𝐟}}^c,\text{diag}(K_{\ast \ast })-{(K_{\ast }^c)}^{\mathrm{T}}{(K^c+{(W^c)}^{-1})}^{-1}{(K_{\ast }^c)}^{\mathrm{T}}),$$ (15)

where $f_{\ast }^c$, $K_{\ast }^c$, and $W^c$ represent the class-relevant information only. Finally, (11) is replaced with

$$p(𝐲|𝐱,𝜽)=-\frac{1}{2}{\widehat{𝐟}}^{\mathrm{T}}{K}^{-1}\widehat{𝐟}+{𝐲}^{\mathrm{T}}\widehat{𝐟}-\sum _{i=1}^n\log \left[\sum _{c=1}^C\exp {\widehat{f}}_i^c\right]-\frac{1}{2}\log \left(|K|\cdot |K^{-1}+W|\right).$$ (16)

We won’t present an example of multi-class GPC, but hopefully you get the idea.

5 Discussion

As with GPR, classification can be extended to accept $x$ values with multiple dimensions, while keeping most of the mathematics unchanged. Other possible extensions include using the expectation propagation method in lieu of the Laplace approximation as mentioned previously, putting confidence intervals on the classification probabilities, calculating the derivatives of (16) to aid the optimizer, or using the variational Gaussian process classifiers described in MacKay (1998), to name but four extensions.

Second, we repeat the Bayesian call from the previous document to integrate over a range of possible covariance function parameters. This should be done regardless of how much prior knowledge is available — see for example Chapter 5 of Sivia and Skilling (2006) on how to choose priors in the most opaque situations.

Third, we’ve again spared you a few practical algorithmic details; computer code is available at www.gaussianprocess.org/gpml, with examples.

 

Gaussian Processes for Dimensionality Reduction:
A Quick Introduction

M. Ebden, August 2015

Prerequisite reading: Gaussian processes for regression

 

Suppose you wanted to learn a low-dimensional representation of a dataset for ease of interpretation, sort of like how your shadow is a simplified, two-dimensional representation of a three-dimensional object. Let the original data be matrix $𝐘$ with $n$ rows (observations) of $d$ dimensions each, and let the new representation be $𝐗$ with $n$ rows of $q$ dimensions each ().

To learn this low-dimensional representation, we’ll assume that for the $i$th dimension of $𝐘$ the $n$ elements in ${𝐲}_{:,i}$ are a sample from a Gaussian process based on the low-dimensional space. Specifically, we’ll use a zero-mean Gaussian process, $𝒢𝒫({\mathrm{𝟎}}_{n\times 1},k(𝐱,{𝐱}^{\prime }))$, keeping the same model for each dimension of $𝐘$. We’ll choose for $k(\cdot )$ the squared exponential kernel (a.ka. RBF kernel) in order to ensure that points close in $𝐗$ will be close in $𝐘$. Renaming the noise variance from ${\sigma }_n^2$ (in the first report) to $1/\beta$, the kernel is:

$$k(𝐱,{𝐱}^{\prime })={\sigma }^2\exp \left[\frac{-{|𝐱-{𝐱}^{\prime }|}^2}{2l^2}\right]+\frac{\delta (x,x^{\prime })}{\beta }.$$

The covariance matrix $𝐊$ for this GP is constructed as per (4) in the first report, and we have no need for a ${𝐊}_{\ast }$ or ${𝐊}_{\ast \ast }$ from (5) because there are no points in $𝐘$ to predict. The tuple of kernel parameters this time is $𝜽=\{\sigma ,l,\beta \}$.

We’ll further assume that the GPs behind each of the $d$ dimensions have been sampled independently. Therefore, the likelihood of the observed values in $𝐘$ is the product of $d$ independent GPs:

$$\text{p}(𝐘|𝐗,𝜽)=\prod _{i=1}^d\frac{\exp \left[-\frac{1}{2}{𝐲}_{:,i}^{\text{T}}{𝐊}^{-1}{𝐲}_{:,i}\right]}{{(2\pi )}^{n/2}{|𝐊|}^{1/2}}.$$

Then we can adjust $𝐗$ and $𝜽$ to maximize this likelihood. If we were only adjusting $𝐗$, it would be akin to determining the shadow which is most likely to correspond to your body’s shape. Because we’re adjusting $𝜽$ as well, imagine manipulating the light source and the surface for the shadow to appear on as well.

Figure 1 gives an example outcome of the most likely $𝐗$ for a certain $𝐘$, with $n=17$, $d=2$, and $q=1$. $𝐘$ was preprocessed to have zero mean and the same variance in each dimension. Using an optimizer based on scaled conjugate gradients, the optimal $𝜽$ and $𝐗$ were then found. The former is $\{\sigma ,l,\beta \}=\{1.05,3.3\times {10}^{-4},93\}$ and the latter is given in Figure 1(b). In this example, $𝐗$ is more than just a ‘shadow’ (simple projection) of $𝐘$: the above method of dimensionality reduction is powerful enough to have learned that the points in the left curve in Figure 1(a) should be grouped together while the points in the right curve form a separate group; no simple lamplight projection can achieve this. The price we pay for this flexibility is in the high number of parameters being fit: the total here is $3+17\times 1=20$, i.e. $𝜽$ plus the one-dimensional values in $𝐗$. Usually $𝐗$ are referred to as latent variables rather than parameters, and our overall setup is called the ‘Gaussian process latent variable model’ (GP-LVM). The technique was first presented by Neil Lawrence in 2004; he examined a set of oil-flow data in which $n=1000$ and $d=12$, which we reduce to $q=2$ in Figure 2.

This tutorial on the GP-LVM is provided by Mind Foundry, a technology spin-out from the University of Oxford.

Figure 1: (a) Seventeen data points, $𝐘$, each of two dimensions.
d (b) A one-dimensional projection, $𝐗$, of those points.

Figure 2: Although the GP-LVM algorithm uses unsupervised learning, the three classes of data in this originally 12-D dataset generally occupy different regions.