Bayesian Classification With Gaussian Processes.
AbstractWe consider the problem of assigning an input vector to one of m classes by predicting P(c|${\schmi x}$) for m. For a two-class problem, the probability of class one given ${\schmi x}$ is estimated by (y(${\schmi x}$)), where Gaussian process prior is placed on y(${\schmi x}$), and is combined with the training data to obtain predictions for new ${\schmi x}$ points. We provide a Bayesian treatment, integrating over uncertainty in y and in the parameters that control the Gaussian process prior; the necessary integration over y is carried out using Laplace's approximation. The method is generalized to multiclass problems (m > 2) using the softmax function. We demonstrate the effectiveness of the method on a number of datasets.
Introduction
We consider the problem of assigning an input vector x to one out of m classes by predicting P (cjx) for
classic example of this method is logistic regression. For a two-class problem, the probability
of class 1 given x is estimated by oe(w T x+ b), where \Gammay ). However, this method is not at all
"flexible", i.e. the discriminant surface is simply a hyperplane in x-space. This problem can be overcome, to
some extent, by expanding the input x into a set of basis functions fOE(x)g, for example quadratic functions of
the components of x. For a high-dimensional input space there will be a large number of basis functions, each
one with an associated parameter, and one risks "overfitting" the training data. This motivates a Bayesian
treatment of the problem, where the priors on the parameters encourage smoothness in the model.
Putting priors on the parameters of the basis functions indirectly induces priors over the functions that
can be produced by the model. However, it is possible (and we would argue, perhaps more natural) to put
priors directly over the functions themselves. One advantage of function-space priors is that they can impose
a general smoothness constraint without being tied to a limited number of basis functions. In the regression
case where the task is to predict a real-valued output, it is possible to carry out non-parametric regression
using Gaussian Processes (GPs); see, e.g. [25], [28]. The solution for the regression problem under a GP
prior (and Gaussian noise model) is to place a kernel function on each training data point, with coefficients
determined by solving a linear system. If the parameters ' that describe the Gaussian process are unknown,
Bayesian inference can be carried out for them, as described in [28].
The Gaussian Process method can be extended to classification problems by defining a GP over y, the
input to the sigmoid function. This idea has been used by a number of authors, although previous treatments
typically do not take a fully Bayesian approach, ignoring uncertainty in both the posterior distribution of y
given the data, and uncertainty in the parameters '. This paper attempts a fully Bayesian treatment of the
problem, and also introduces a particular form of covariance function for the Gaussian process prior which,
we believe, is useful from a modelling point of view.
The structure of the remainder of the paper is as follows: Section 2 discusses the use of Gaussian processes
for regression problems, as this is essential background for the classification case. In Section 3 we describe
the application of Gaussian processes to two-class classification problems, and extend this to multiple-class
problems in section 4. Experimental results are presented in section 5, followed by a discussion in section 6.
This paper is a revised and expanded version of [1].
Gaussian Processes for regression
It will be useful to first consider the regression problem, i.e. the prediction of a real valued output y
for a new input value x   , given a set of training data ng. This is of relevance because
our strategy will be to transform the classification problem into a regression problem by dealing with the
input values to the logistic transfer function.
A stochastic process prior over functions allows us to specify, given a set of inputs, x the distribution
over their corresponding outputs
We denote this prior over functions as P (y), and similarly, P (y   ; y) for the
joint distribution including y   . If we also specify P (tjy), the probability of observing the particular values
actual values y (i.e. a noise model) then we have that
Z
Z
Z
Hence the predictive distribution for y   is found from the marginalization of the product of the prior and
the noise model. Note that in order to make predictions it is not necessary to deal directly with priors over
function space, only n- or n + 1-dimensional joint densities. However, it is still not easy to carry out these
calculations unless the densities involved have a special form.
If P (tjy) and P (y   ; y) are Gaussian then P (y   jt) is a Gaussian whose mean and variance can be calculated
using matrix computations involving matrices of size n \Theta n. Specifying P (y   ; y) to be a multidimensional
Gaussian (for all values of n and placements of the points x   means that the prior over functions
is a Gaussian process. More formally, a stochastic process is a collection of random variables fY (x)jx 2 Xg
indexed by a set X . In our case X will be the input space with dimension d, the number of inputs.
A GP is a stochastic process which can be fully specified by its mean function its
covariance function C(x; x any finite set of Y -variables will have a
joint multivariate Gaussian distribution. Below we consider GPs which have -(x) j 0.
If we further assume that the noise model P (tjy) is Gaussian with mean zero and variance oe 2 I , then the
predicted mean and variance at x   are given by
y
e.g. [25]).
2.1 Parameterizing the covariance function
There are many reasonable choices for the covariance function. Formally, we are required to specify functions
which will generate a non-negative definite covariance matrix for any set of points From a
modelling point of view we wish to specify covariances so that points with nearby inputs will give rise to
similar predictions. We find that the following covariance function works well:
l
where x l is the lth component of x and is the vector of parameters
that are needed to define the covariance function. Note that ' is analogous to the hyperparameters in a
neural network. We define the parameters to be the log of the variables in equation (4) since these are
positive scale-parameters. This covariance function can be obtained from a network of Gaussian radial basis
functions in the limit of an infinite number of hidden units [27].
The w l parameters in equation 4 allow a different length scale on each input dimension. For irrelevant
inputs, the corresponding w l will become small, and the model will ignore that input. This is closely related
to the Automatic Relevance Determination (ARD) idea of MacKay [10] and Neal [15]. The v 0 variable
specifies the overall scale of the prior. v 1 specifies the variance of a zero-mean offset which has a Gaussian
distribution.
The Gaussian process framework allows quite a wide variety of priors over functions. For example, the
Ornstein-Uhlenbeck process (with covariance function C(x; x very rough sample paths
which are not mean-square differentiable. On the other hand the squared exponential covariance function of
equation 4 gives rise to an infinitely m.s. differentiable process. In general we believe that the GP method
is a quite general-purpose route for imposing prior beliefs about the desired amount of smoothness. For
reasonably high-dimensional problems, this needs to be combined with other modelling assumptions such as
ARD. Another modelling assumption that may be used is to build up the covariance function as a sum of
covariance functions, each one of which may depend on only some of the input variables (see section 3.3 for
further details).
2.2 Dealing with parameters
Given a covariance function it is straightforward to make predictions for new test points. However, in
practical situations we are unlikely to know which covariance function to use. One option is to choose
a parametric family of covariance functions (with a parameter vector ') and then either to estimate the
parameters (for example, using the method of maximum likelihood) or to use a Bayesian approach where a
posterior distribution over the parameters is obtained.
These calculations are facilitated by the fact that the log likelihood l = log P (Dj') can be calculated
analytically as
log 2-; (5)
where ~
Kj denotes the determinant of ~
K. It is also possible to express analytically the
partial derivatives of the log likelihood with respect to the parameters
@l
~
~
(see, e.g. [11]).
Given l and its derivatives with respect to ' it is straightforward to feed this information to an optimization
package in order to obtain a local maximum of the likelihood.
In general one may be concerned about making point estimates when the number of parameters is large
relative to the number of data points, or if some of the parameters may be poorly determined, or if there
may be local maxima in the likelihood surface. For these reasons the Bayesian approach of defining a prior
s

Figure

1: -(x) is obtained from y(x) by "squashing" it through the sigmoid function oe.
distribution over the parameters and then obtaining a posterior distribution once the data D has been seen is
attractive. To make a prediction for a new test point x   one simply averages over the posterior distribution
Z
For GPs it is not possible to do this integration analytically in general, but numerical methods may be used.
If ' is of sufficiently low dimension, then techniques involving grids in '-space can be used.
If ' is high-dimensional it is very difficult to locate the regions of parameter-space which have high
posterior density by gridding techniques or importance sampling. In this case Markov chain Monte Carlo
methods may be used. These work by constructing a Markov chain whose equilibrium distribution
is the desired distribution P ('jD); the integral in equation 7 is then approximated using samples from the
Markov chain.
Two standard methods for constructing MCMC methods are the Gibbs sampler and Metropolis-Hastings
algorithms (see, e.g., [5]). However, the conditional parameter distributions are not amenable to Gibbs
sampling if the covariance function has the form given by equation 4, and the Metropolis-Hastings algorithm
does not utilize the derivative information that is available, which means that it tends to have an inefficient
random-walk behaviour in parameter-space. Following the work of Neal [15] on Bayesian treatment of neural
networks, Williams and Rasmussen [28] and Rasmussen [17] have used the Hybrid Monte Carlo (HMC)
method of Duane et al [4] to obtain samples from P ('jD). The HMC algorithm is described in more detail
in

Appendix

D.
3 Gaussian Processes for two-class classification
For simplicity of exposition we will first present our method as applied to two-class problems; the extension
to multiple classes is covered in section 4.
By using the logistic transfer function to produce an output which can be interpreted as -(x), the
probability of the input x belonging to class 1, the job of specifying a prior over functions - can be transformed
into that of specifying a prior over the input to the transfer function, which we shall call the activation, and
denote by y (see Figure 1). For the two-class problem we can use the logistic function
\Gammay ). We will denote the probability and activation corresponding to input x i by - i and y i
respectively. Fundamentally, the GP approaches to classification and regression problems are similar, except
that the error model which is t - N(y; oe 2 ) in the regression case is replaced by t - Bern(oe(y)). The choice
of v 0 in equation 4 will affect how "hard" the classification is; i.e. if -(x) hovers around 0:5 or takes on the
extreme values of 0 and 1.
Previous and related work to this approach is discussed in section 3.3.
As in the regression case there are now two problems to address (a) making predictions with fixed
parameters and (b) dealing with parameters. We shall discuss these issues in turn.
3.1 Making predictions with fixed parameters
To make predictions when using fixed parameters we would like to compute -
R
requires us to find P (-   for a new input x   . This can be done by finding the distribution
is the activation of -   ) and then using the appropriate Jacobian to transform the distribution.
Formally the equations for obtaining P (y   jt) are identical to equations 1, 2, and 3. However, even if we use a
GP prior so that P (y   ; y) is Gaussian, the usual expression for
classification
data (where the t's take on values of 0 or 1), means that the marginalization to obtain P (y   jt) is no longer
analytically tractable.
Faced with this problem there are two routes that we can follow: (i) to use an analytic approximation
to the integral in equations 1-3 or (ii) to use Monte Carlo methods, specifically MCMC methods, to approximate
it. Below we consider an analytic approximation based on Laplace's approximation; some other
approximations are discussed in section 3.3.
In Laplace's approximation, the integrand P (y   ; yjt; ') is approximated by a Gaussian distribution
centered at a maximum of this function with respect to y   ; y with an inverse covariance matrix given by
\Gammarr log P (y   ; yjt; '). Finding a maximum can be carried out using the Newton-Raphson iterative method
on y, which then allows the approximate distribution of y   to be calculated. Details of the maximization
procedure can be found in Appendix A.
3.2 Integration over the parameters
To make predictions we integrate the predicted probabilities over the posterior P ('jt) / P (tj')P ('), as we
saw in 2.2. For the regression problem P (tj') can be calculated exactly using P
R P (tjy)P (yj')dy,
but this integral is not analytically tractable for the classification problem. Let
Using log
log 2-: (8)
By using Laplace's approximation about the maximum ~
y we find that
log
log 2-:
We denote the right-hand side of this equation by log P a (tj') (where a stands for approximate).
The integration over '-space also cannot be done analytically, and we employ a Markov Chain Monte
Carlo method. Following Neal [15] and Williams and Rasmussen [28] we have used the Hybrid Monte Carlo
(HMC) method of Duane et al [4] as described in Appendix D. We use log P a (tj') as an approximation for
log P (tj'), and use broad Gaussian priors on the parameters.
3.3 Previous and related work
Our work on Gaussian processes for regression and classification developed from the observation in [15]
that a large class of neural network models converge to GPs in the limit of an infinite number of hidden
units. The computational Bayesian treatment of GPs can be easier than for neural networks. In the
regression case an infinite number of weights are effectively integrated out, and one ends up dealing only
with the (hyper)parameters. Results from [17] show that Gaussian processes for regression are comparable
in performance to other state-of-the-art methods.
Non-parametric methods for classification problems can be seen to arise from the combination of two
different strands of work. Starting from linear regression, McCullagh and Nelder [12] developed generalized
linear models (GLMs). In the two-class classification context, this gives rise to logistic regression. The other
strand of work was the the development of non-parametric smoothing for the regression problem. Viewed
as a Gaussian process prior over functions this can be traced back at least as far as the work of Kolmogorov
and Wiener in the 1940s. Gaussian process prediction is well known in the geostatistics field (see, e.g. [3])
where it is known as "kriging". Alternatively, by considering "roughness penalties" on functions, one can
obtain spline methods; for recent overviews, see [25] and [8]. There is a close connection between the GP
and roughness penalty views, as explored in [9]. By combining GLMs with non-parametric regression one
obtains what we shall call a non-parametric GLM method for classification. Early references to this method
include [21] and [16], and discussions can also be found in texts such as [8] and [25].
There are two differences between the non-parametric GLM method as it is usually described and a
Bayesian treatment. Firstly, for fixed parameters the non-parametric GLM method ignores the uncertainty
in y   and hence the need to integrate over this (as described in section 3.1).
The second difference relates to the treatment of the parameters '. As discussed in section 2.2, given
parameters ', one can either attempt to obtain a point estimate for the parameters or to carry out an
integration over the posterior. Point estimates may be obtained by maximum likelihood estimation of ',
or by cross-validation or generalized cross-validation (GCV) methods, see e.g. [25, 8]. One problem with
CV-type methods is that if the dimension of ' is large, then it can be computationally intensive to search
over a region/grid in parameter-space looking for the parameters that maximize the criterion. In a sense
the HMC method described above are doing a similar search, but using gradient information 1 , and carrying
out averaging over the posterior distribution of parameters. In defence of (G)CV methods, we note Wahba's
comments (e.g. in [26], referring back to [24]) that these methods may be more robust against an unrealistic
prior.
One other difference between the kinds of non-parametric GLM models usually considered and our method
is the exact nature of the prior that is used. Often the roughness penalties used are expressed in terms of a
penalty on the kth derivative of y(x), which gives rise to a power law power spectrum for the prior on y(x).
There can also be differences over parameterization of the covariance function; for example it is unusual to
find parameters like those for ARD introduced in equation 4 in non-parametric GLM models. On the other
hand, Wahba et al [26] have considered a smoothing spline analysis of variance (SS-ANOVA) decomposition.
In Gaussian process terms, this builds up a prior on y as a sum of priors on each of the functions in the
decomposition
ff
y ff
The important point is that functions involving all orders of interaction (from univariate functions, which
on their own give rise to an additive model) are included in this sum, up to the full interaction term which
is the only one that we are using. From a Bayesian point of view questions as to the kinds of priors that are
appropriate is an interesting modelling issue.
There has also been some recent work which is related to the method presented in this paper. In section
3.1 we mentioned that it is necessary to approximate the integral in equations 1-3 and described the use of
Laplace's approximation.
Following the preliminary version of this paper presented in [1], Gibbs and MacKay [7] developed an alternative
analytic approximation, by using variational methods to find approximating Gaussian distributions
that bound the marginal likelihood P (tj') above and below. These approximate distributions are then used
to predict P (y   jt; ') and thus -
-(x   ). For the parameters, Gibbs and MacKay estimated ' by maximizing
their lower bound on P (tj').
It is also possible to use a fully MCMC treatment of the classification problem, as discussed in the recent
paper of Neal [14]. His method carries out the integrations over the posterior distributions of y and '
simultaneously. It works by generating samples from P (y; 'jD) in a two stage process. Firstly, for fixed ',
each of the n individual y i 's are updated sequentially using Gibbs sampling. This ``sweep'' takes time O(n 2 )
once the matrix K \Gamma1 has been computed (in time O(n 3 )), so it actually makes sense to perform quite a few
Gibbs sampling scans between each update of the parameters, as this probably makes the Markov chain mix
faster. Secondly, the parameters are updated using the Hybrid Monte Carlo method. To make predictions,
one averages over the predictions made by each
It would be possible to obtain derivatives of the CV-score with respect to ', but this has not, to our knowledge, been used
in practice.
4 GPs for multiple-class classification
The extension of the preceding framework to multiple classes is essentially straightforward, although notationally
more complex.
Throughout we employ a one-of-m class coding scheme 2 , and use the multi-class analogue of the logistic
function-the softmax function-to describe the class probabilities. The probability that an instance labelled
by i is in class c is denoted by - i
c , so that an upper index to denotes the example number, and a lower index
the class label. Similarly, the activations associated with the probabilities are denoted by y i
c . Formally, the
link function relates the activations and probabilities through
c
which automatically enforces the constraint
1. The targets are similarly represented by t i
c , and are
specified using a one-of-m coding.
The log likelihood takes the form
c , which for the softmax link function gives
c
As for the two class case, we shall assume that the GP prior operates in activation space; that is we specify
the correlations between the activations y i
c .
One important assumption we make is that our prior knowledge is restricted to correlations between the
activations of a particular class. Whilst there is no difficulty in extending the framework to include inter-class
correlations, we have not yet encountered a situation where we felt able to specify such correlations.
Formally, the activation correlations take the form,
hy i
c (12)
where K i;i 0
c is the element of the covariance matrix for the cth class. Each individual correlation matrix
K c has the form given by equation 4 for the two-class case. We shall use a separate set of parameters for
each class. The use of m independent processes to perform the classification is redundant, but forcing the
activations of one process to be (say) zero would introduce an arbitrary asymmetry into the prior.
For simplicity, we introduce the augmented vector notation,
where, as in the two-class case, y
c denotes the activation corresponding to input x   for class c; this notation
is also used to define t + and -+ . In a similar manner, we define y, t and - by excluding the values
corresponding to the test point x   . Let y
With this definition of the augmented vectors, the GP prior takes the form,
ae
oe
where, from equation 12, the covariance matrix K + is block diagonal in the matrices, K
m . Each
individual matrix K
c expresses the correlations of activations within class c.
As in the two-class case, to use Laplace's approximation we need to find the mode of P (y jt). The
procedure is described in Appendix C. As for the two-class case, we make predictions for -(x   ) by averaging
the softmax function over the Gaussian approximation to the posterior distribution of y   . At present, we
simply estimate this integral using 1000 draws from a Gaussian random vector generator.
That is, the class is represented by a vector of length m with zero entries everywhere except for the correct component
which contains 1.
5 Experimental results
When using the Newton-Raphson algorithm, - was initialized each time with entries 1=m, and iterated until
the mean relative difference of the elements of W between consecutive iterations was less than 10 \Gamma4 .
For the HMC algorithm, the same step size " is used for all parameters, and should be as large as possible
while keeping the rejection rate low. We have used a trajectory made up of leapfrog steps, which gave
a low correlation between successive states. The priors over parameters were set to be Gaussian with a mean
of \Gamma3 and a standard deviation of 3. In all our simulations a step size produced a low rejection rate
(! 5%). The parameters corresponding to the w l 's were initialized to \Gamma2 and that for v 0 to 0. The sampling
procedure was run for 200 iterations, and the first third of the run was discarded; this "burn-in" is intended
to give the parameters time to come close to their equilibrium distribution. Tests carried out using the
R-CODA package 3 on the examples in section 5.1 suggested that this was indeed effective in removing the
transients, although we note that it is widely recognized (see, e.g. [2]) that determining when the equilibrium
distribution has been reached is a difficult problem. Although the number of iterations used is much less
than typically used for MCMC methods it should be remembered that (i) each iteration involves
leapfrog steps and (ii) that by using HMC we aim to reduce the "random walk" behaviour seen in methods
such as the Metropolis algorithm. Autocorrelation analysis for each parameter indicated, in general, that
low correlation was obtained after a lag of a few iterations.
The MATLAB code which we used to run our experiments is available from
ftp://cs.aston.ac.uk/neural/willicki/gpclass/.
5.1 Two classes
We have tried out our method on two well known two class classification problems, the Leptograpsus crabs
and Pima Indian diabetes datasets 4 . We first rescale the inputs so that they have mean of zero and unit
variance on the training set. Our Matlab implementations for the HMC simulations for both tasks each
take several hours on a SGI Challenge machine (200MHz R10000), although good results can be obtained in
much less time. We also tried a standard Metropolis MCMC algorithm for the Crabs problem, and found
similar results, although the sampling by this method is somewhat slower than that for HMC.
The results for the Crabs and Pima tasks, together with comparisons with other methods (from [20] and
[18]) are given in Tables 1 and 2 respectively. The tables also include results obtained for Gaussian processes
using (a) estimation of the parameters by maximizing the penalised likelihood (found using 20 iterations of
a scaled conjugate gradient optimiser) and (b) Neal's MCMC method. Details of the set-up used for Neal's
method are given in Appendix E.
In the Leptograpsus crabs problem we attempt to classify the sex of crabs on the basis of five anatomical
attributes, with an optional additional colour attribute. There are 50 examples available for crabs of each sex
and colour, making a total of 200 labelled examples. These are split into a training set of 20 crabs of each sex
and colour, making 80 training examples, with the other 120 examples used as the test set. The performance
of our GP method is equal to the best of the other methods reported in [20], namely a 2 hidden unit
neural network with direct input to output connections, a logistic output unit and trained with maximum
likelihood (Network(1) in Table 1). Neal's method gave a very similar level of performance. We also found
that estimating the parameters using maximum penalised likelihood (MPL) gave similar performance with
less than a minute of computing time.
For the Pima Indians diabetes problem we have used the data as made available by Prof. Ripley, with
his training/test split of 200 and 332 examples respectively [18]. The baseline error obtained by simply
classifying each record as coming from a diabetic gives rise to an error of 33%. Again, ours and Neal's
GP methods are comparable with the best alternative performance, with an error of around 20%. It is
encouraging that the results obtained using Laplace's approximation and Neal's method are similar 5 . We
also estimated the parameters using maximum penalised likelihood, rather than Monte Carlo integration.
The performance in this case was a little worse, with 21.7% error, but for only 2 minutes computing time.
3 Available from the Comprehensive R Archive Network at http://www.ci.tuwien.ac.at.
4 Available from http://markov.stats.ox.ac.uk/pub/PRNN.
5 The performance obtained by Gibbs and MacKay in [7] was similar. Their method made 4 errors in the crab task (with
colour given), and 70 errors on the Pima dataset.
Method Colour given Colour not given
Neural Network(1) 3 3
Neural Network(2) 5 3
Linear Discriminant 8 8
Logistic regression 4 4
PP regression (4 ridge
Gaussian Process (Laplace 3 3
Approximation, HMC)
Gaussian Process (Laplace 4 3
Approximation, MPL)
Gaussian Process (Neal's method) 4 3

Table

1: Number of test errors for the Leptograpsus crabs task. Comparisons are taken from from Ripley (1996)
and Ripley (1994) respectively. Network(2) used two hidden units and the predictive approach (Ripley, 1993) which
uses Laplace's approximation to weight each network local minimum.
Method Pima Indian diabetes
Neural Network 75+
Linear Discriminant 67
Logistic Regression 66
PP regression (4 ridge functions) 75
Gaussian Mixture 64
Gaussian Process (Laplace 68
Approximation, HMC)
Gaussian Process (Laplace 69
Approximation, MPL)
Gaussian Process (Neal's method) 68

Table

2: Number of test errors on the Pima Indian diabetes task. Comparisons are taken from from Ripley
(1996) and Ripley (1994) respectively. The neural network had one hidden unit and was trained with
maximum likelihood; the results were worse for nets with two or more hidden units (Ripley, 1996).
Analysis of the posterior distribution of the w parameters in the covariance function (equation
be informative. Figure 5.1 plots the posterior marginal mean and 1 standard deviation error bars for each
of the seven input dimensions. Recalling that the variables are scaled to have zero mean and unit variance,
it would appear that variables 1 and 3 have the shortest lengthscales (and therefore the most variability)
associated with them.
5.2 Multiple classes
Due to the rather long time taken to run our code, we were only able to test it on relatively small problems,
by which we mean only a few hundred data points and several classes. Furthermore, we found that a full
Bayesian integration over possible parameter settings was beyond our computational means, and we therefore
had to be satisfied with a maximum penalised likelihood approach. Rather than using the potential and its
gradient in a HMC routine, we now simply used them as inputs to a scaled conjugate gradient optimiser
(based on [13]) instead, attempting to find a mode of the class posterior, rather than to average over the
posterior distribution.
We tested the multiple class method on the Forensic Glass dataset described in [18]. This is a dataset
of 214 examples with 9 inputs and 6 output classes. Because the dataset is so small, the performance is
Figure

2: Plot of the log w parameters for the Pima dataset. The circle indicates the posterior marginal mean
obtained from the HMC run (after burn-in), with one standard deviation error bars. The square symbol
shows the log w-parameter values found by maximizing the penalized likelihood. The variables are 1. the
number of pregnancies, 2. plasma glucose concentration, 3. diastolic blood pressure, 4. triceps skin fold
thickness, 5. body mass index, 6. diabetes pedigree function, 7. age. For comparison, Wahba et al (1995)
using generalized linear regression, found that variables 1, 2 5 and 6 were the most important.
estimated from using 10-fold cross validation. Computing the penalised maximum likelihood estimate of our
multiple GP method took approximately 24 hours on our SGI Challenge and gave a classification error rate
of 23.3%. As we see from Table 3, this is comparable to the best of the other methods. The performance of
Neal's method is surprisingly poor; this may be due to the fact that we allow separate parameters for each
of the y processes, while these are constrained to be equal in Neal's code. There are also small but perhaps
significant differences in the specification of the prior (see Appendix E for details).
6 Discussion
In this paper we have extended the work of Williams and Rasmussen [28] to classification problems, and have
demonstrated that it performs well on the datasets we have tried. We believe that the kinds of Gaussian
Method Forensic Glass
Neural Network (4HU) 23.8%
Linear Discriminant 36%
PP regression (5 ridge functions) 35%
Gaussian Mixture 30.8%
Decision Tree 32.2%
Gaussian Process (LA, MPL) 23.3%
Gaussian Process (Neal's method) 31.8%

Table

3: Percentage of test error for the Forensic Glass problem. See Ripley (1996) for details of the methods.
process prior we have used are more easily interpretable than models (such as neural networks) in which
the priors are on the parameterization of the function space. For example, the posterior distribution of the
ARD parameters (as illustrated in Figure 5.1 for the Pima Indians diabetes problem) indicates the relative
importance of various inputs. This interpretability should also facilitate the incorporation of prior knowledge
into new problems.
There are quite strong similarities between GP classifiers and support-vector machines (SVMs) [23]. The
SVM uses a covariance kernel, but differs from the GP approach by using a different data fit term (the
maximum margin), so that the optimal y is found using quadratic programming. The comparison of these
two algorithms is an interesting direction for future research.
A problem with methods based on GPs is that they require computations (trace, determinants and linear
solutions) involving n \Theta n matrices, where n is the number of training examples, and hence run into problems
on large datasets. We have looked into methods using Bayesian numerical techniques to calculate the trace
and determinant [22, 6], although we found that these techniques did not work well for the (relatively) small
size problems on which we tested our methods. Computational methods used to speed up the quadratic
programming problem for SVMs may also be useful for the GP classifier problem. We are also investigating
the use of different covariance functions and improvements on the approximations employed.

Acknowledgements

We thank Prof. B. Ripley for making available the Leptograpsus crabs, Pima Indian diabetes and Forensic
Glass datasets. This work was partially supported by EPSRC grant GR/J75425, Novel Developments in
Learning Theory for Neural Networks, and much of the work was carried out at Aston University. The
authors gratefully acknowledge the hospitality provided by the Isaac Newton Institute for Mathematical
Sciences (Cambridge, UK) where this paper was written up. We thank Mark Gibbs, David MacKay and
Radford Neal for helpful discussions, and the anonymous referees for their comments which helped improve
the paper.


Appendix

Maximizing case
We describe how to find iteratively the vector y + so that P (y + jt) is maximized. This material is also
covered in [8] x5.3.3 and [25] x9.2. We provide it here for completeness and so that the terms in equation 9
are well-defined.
the complete set of activations. By Bayes' theorem log
log As P (t) does not depend on y + (it is just
a normalizing factor), the maximum of P (y + jt) is found by maximizing \Psi + with respect to y + . Using
log
log 2- (14)
where K+ is the covariance matrix of the GP evaluated at x . \Psi is defined similarly in equation
8. K+ can be partitioned in terms of an n \Theta n matrix K, a n \Theta 1 vector k and a scalar k   , viz.
As y   only enters into equation 14 in the quadratic prior term and has no data point associated with it,
maximizing with respect to y + can be achieved by first maximizing \Psi with respect to y and then doing
the further quadratic optimization to determine y   . To find a maximum of \Psi we use the Newton-Raphson
iteration y new Differentiating equation 8 with respect to y we find
where the 'noise' matrix is given by This results in the iterative equation,
To avoid unnecessary inversions, it is usually more convenient to rewrite this in the form
Note that \Gammarr\Psi is always positive definite, so that the optimization problem is convex.
Given a converged solution ~ y for y, y   can easily be found using y
-), as
is the W with a
zero appended in the (n diagonal position. Given the mean and variance of y   it is then easy to find
R
the mean of the distribution of P (-   jt). In order to calculate the Gaussian integral over
the logistic sigmoid function, we employ an approximation based on the expansion of the sigmoid function
in terms of the error function. As the Gaussian integral of an error function is another error function, this
approximation is fast to compute. Specifically, we use a basis set of five scaled error functions to interpolate
the logistic sigmoid at chosen points 6 . This gives an accurate approximation (to to the desired integral
with a small computational cost.
The justification of Laplace's approximation in our case is somewhat different from the argument usually
put forward, e.g. for asymptotic normality of the maximum likelihood estimator for a model with a finite
number of parameters. This is because the dimension of the problem grows with the number of data points.
However, if we consider the "infill asymptotics" (see, e.g. [3]), where the number of data points in a bounded
region increases, then a local average of the training data at any point x will provide a tightly localized
estimate for -(x) and hence y(x) (this reasoning parallels more formal arguments found in [29]). Thus we
would expect the distribution P (y) to become more Gaussian with increasing data.


Appendix

B: Derivatives of log P a (tj') wrt '.
For both the HMC and MPL methods we require the derivative of l a log P a (tj') with respect to components
of ', for example ' k . This derivative will involve two terms, one due to explicit dependencies of l a =
log 2- on ' k , and also because a change in ' will cause a change in ~ y. However,
as ~ y is chosen so that r\Psi(y)j y= ~
@l a
\Gamma2
@ log jK
The dependence of jK \Gamma1 +W j on ~ y arises through the dependence of W on ~
-, and hence ~ y. By differentiating
~
-), one obtains
@~ y
and hence the required derivative can be calculated.


Appendix

Maximizing
Multiple-class case
The GP prior and likelihood, defined by equations 13 and 11, define the posterior distribution of activations,
jt). As in Appendix A we are interested in a Laplace approximation to this posterior, and therefore
need to find the mode with respect to y Dropping unnecessary constants, the multi-class analogue of
equation 14 is
c
exp y i
6 In detail, we used the basis functions erf(-x)) for These were used to interpolate oe(x) at
By the same principle as in Appendix A, we define \Psi by analogy with equation 8, and first optimize \Psi with
respect to y, afterwards performing the quadratic optimization of \Psi + with respect to y   .
In order to optimize \Psi with respect to y, we make use of the Hessian given by
where K is the mn \Theta mn block-diagonal matrix with blocks K c , m. Although this is in the same
form as for the two class case, equation 17, there is a slight change in the definition of the 'noise' matrix,
W . A convenient way to define W is by introducing the matrix \Pi which is a mn \Theta n matrix of the form
Using this notation, we can write the noise matrix in the form of a
diagonal matrix and an outer product,
As in the two-class case, we note that \Gammarr\Psi is again positive definite, so that the optimization problem is
convex.
The update equation for iterative optimization of \Psi with respect to the activations y then follows the
same form as that given by equation 18. The advantage of the representation of the noise matrix in equation
23 is that we can then invert matrices and find their determinants using the identities,
and
As A is block-diagonal, it can be inverted blockwise. Thus, rather than
requiring determinants and inverses of a mn \Theta mn matrix, we only need to carry out expensive matrix
computations on n \Theta n matrices. The resulting update equations for y are then of the same form as given
in equation 18, where the noise matrix and covariance matrices are now in their multiple class form.
Essentially, these are all the results needed to generalize the method to the multiple-class problem.
Although, as we mentioned above, the time complexity of the problem does not scale with the m 3 , but
rather m (due to the identities in equations 24, 25), calculating the function and its gradient is still rather
expensive. We have experimented with several methods of mode finding for the Laplace approximation. The
advantage of the Newton iteration method is its fast quadratic convergence. An integral part of each Newton
step is the calculation of the inverse of a matrix M acting upon a vector, i.e. M \Gamma1 b . In order to speed up
this particular step, we used a conjugate gradient (CG) method to solve iteratively the corresponding linear
system b. As we repeatedly need to solve the system (because W changes as y is updated), it saves
time not to run the CG method to convergence each time it is called. In our experiments the CG algorithm
was terminated when 1=n
The calculation of the derivative of log P a (tj') wrt ' in the multiple-class case is analogous to the two-class
case described in Appendix B.


Appendix

D: Hybrid Monte Carlo
HMC works by creating a fictitious dynamical system in which the parameters are regarded as position
variables, and augmenting these with momentum variables p. The purpose of the dynamical system is to
give the parameters "inertia" so that random-walk behaviour in '-space can be avoided. The total energy,
H , of the system is the sum of the kinetic energy, and the potential energy, E. The potential
energy is defined such that p('jD) / exp(\GammaE), i.e. We sample from the joint
distribution for ' and p given by P ('; p) / exp(\GammaE \Gamma K); the marginal of this distribution for ' is the
required posterior. A sample of parameters from the posterior can therefore be obtained by simply ignoring
the momenta.
Sampling from the joint distribution is achieved by two steps: (i) finding new points in phase space
with near-identical energies H by simulating the dynamical system using a discretised approximation to
Hamiltonian dynamics, and (ii) changing the energy H by Gibbs sampling the momentum variables.
Hamilton's first order differential equations for H are approximated using the leapfrog method which
requires the derivatives of E with respect to '. Given a Gaussian prior on ', log P (') is straightforward to
differentiate. The derivative of log P a (tj') is also straightforward, although implicit dependencies of ~
y (and
hence ~
-) on ' must be taken into account as described in Appendix B. The calculation of the energy can be
quite expensive as for each new ', we need to perform the maximization required for Laplace's approximation,
equation 9. This proposed state is then accepted or rejected using the Metropolis rule depending on the
final energy H   (which is not necessarily equal to the initial energy H because of the discretization).


Appendix

E: Simulation set-up for Neal's code
We used the fbm software available from
http://www.cs.utoronto.ca/~radford/fbm.software.html. For example, the commands used to run the
Pima example are
model-spec pima1.log binary
gp-gen pima1.log fix 0.5 1
mc-spec pima1.log repeat 4 scan-values 200 heatbath hybrid 6 0.5
gp-mc pima1.log 500
which follow closely the example given in Neal's documentation.
The gp-spec command specifies the form of the Gaussian process, and in particular the priors on the
parameters v 0 and the w's (see equation 4). The expression 0.05:0.5 specifies a Gamma-distribution prior
on v 0 , and x0.2:0.5:1 specifies a hierarchical Gamma prior on the w's. Note that a ``jitter'' of 0:1 is also
specified on the prior covariance function; this improves conditioning of the covariance matrix.
The mc-spec command gives details of the MCMC updating procedure. It specifies 4 repetitions of 200
scans of the y values followed by 6 HMC updates of the parameters (using a step-size adjustment factor of
0.5). gp-mc specifies that this is sequence is carried out 500 times.
We aimed for a rejection rate of around 5%. If this was exceeded, the stepsize reduction factor was
reduced and the simulation run again.



--R



Statistics for Spatial Data.
Hybrid Monte Carlo.
Bayesian Data Analysis.
Efficient Implementation of Gaussian Processes.
Variational Gaussian Process Classifiers.
Nonparametric regression and generalized linear models.
A correspondence between Bayesian estimation of stochastic processes and smoothing by splines.
Bayesian Methods for Backpropagation Networks.
Maximum likelihood estimation for models of residual covariance in spatial regression.
Generalized Linear Models.
A scaled conjugate gradient algorithm for fast supervised learning.
Monte Carlo Implementation of Gaussian Process Models for Bayesian Regression and Classification.
Bayesian Learning for Neural Networks.
Automatic Smoothing of Regression Functions in Generalized Linear Models.
Evaluation of Gaussian Processes and Other Methods for Non-linear Regres- sion
Pattern Recognition and Neural Networks.
Statistical aspects of neural networks.
Flexible Non-linear Approaches to Classification
Density Ratios
Bayesian numerical analysis.
The Nature of Statistical Learning Theory.
A Comparison of GCV and GML for Choosing the Smoothing Parameter in the Generalized Spline Smoothing Problem.
Spline Models for Observational Data.
Classification
Computing with infinite networks.
Gaussian processes for regression.
A Comparison of Kriging with Nonparametric Regression Methods.
--TR

--CTR
Christopher K. I. Williams, On a Connection between Kernel PCA and Metric Multidimensional Scaling, Machine Learning, v.46 n.1-3, p.11-19, 2002
S. S. Keerthi , K. B. Duan , S. K. Shevade , A. N. Poo, A Fast Dual Algorithm for Kernel Logistic Regression, Machine Learning, v.61 n.1-3, p.151-165, November  2005
Mrio A. T. Figueiredo, Adaptive Sparseness for Supervised Learning, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.25 n.9, p.1150-1159, September
Koji Tsuda, Propagating distributions on a hypergraph by dual information regularization, Proceedings of the 22nd international conference on Machine learning, p.920-927, August 07-11, 2005, Bonn, Germany
W. D. Addison , R. H. Glendinning, Robust image classification, Signal Processing, v.86 n.7, p.1488-1501, July 2006
Hyun-Chul Kim , Daijin Kim , Zoubin Ghahramani , Sung Yang Bang, Appearance-based gender classification with Gaussian processes, Pattern Recognition Letters, v.27 n.6, p.618-626, 15 April 2006
Yasemin Altun , Alex J. Smola , Thomas Hofmann, Exponential families for conditional random fields, Proceedings of the 20th conference on Uncertainty in artificial intelligence, p.2-9, July 07-11, 2004, Banff, Canada
Wei Chu , Zoubin Ghahramani, Preference learning with Gaussian processes, Proceedings of the 22nd international conference on Machine learning, p.137-144, August 07-11, 2005, Bonn, Germany
Balaji Krishnapuram , Alexander J. Hartemink , Lawrence Carin , Mario A. T. Figueiredo, A Bayesian Approach to Joint Feature Selection and Classifier Design, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.26 n.9, p.1105-1111, September 2004
Wei Chu , S. Sathiya Keerthi , Chong Jin Ong, Bayesian trigonometric support vector classifier, Neural Computation, v.15 n.9, p.2227-2254, September
Yasemin Altun , Thomas Hofmann , Alexander J. Smola, Gaussian process classification for segmenting and annotating sequences, Proceedings of the twenty-first international conference on Machine learning, p.4, July 04-08, 2004, Banff, Alberta, Canada
Hyun-Chul Kim , Jaewook Lee, Clustering Based on Gaussian Processes, Neural Computation, v.19 n.11, p.3088-3107, November 2007
Bart Bakker , Tom Heskes, Task clustering and gating for bayesian multitask learning, The Journal of Machine Learning Research, 4, p.83-99, 12/1/2003
Mark Girolami , Simon Rogers, Variational Bayesian multinomial probit regression with Gaussian process priors, Neural Computation, v.18 n.8, p.1790-1817, August 2006
Liefeng Bo , Ling Wang , Licheng Jiao, Feature Scaling for Kernel Fisher Discriminant Analysis Using Leave-One-Out Cross Validation, Neural Computation, v.18 n.4, p.961-978, April 2006
Volker Tresp, The generalized Bayesian committee machine, Proceedings of the sixth ACM SIGKDD international conference on Knowledge discovery and data mining, p.130-139, August 20-23, 2000, Boston, Massachusetts, United States
Malte Kuss , Carl Edward Rasmussen, Assessing Approximate Inference for Binary Gaussian Process Classification, The Journal of Machine Learning Research, 6, p.1679-1704, 12/1/2005
Lehel Csat , Manfred Opper, Sparse on-line Gaussian processes, Neural Computation, v.14 n.3, p.641-668, March 2002
Manfred Opper , Ole Winther, Gaussian Processes for Classification: Mean-Field Algorithms, Neural Computation, v.12 n.11, p.2655-2684, November 2000
Michael Lindenbaum , Shaul Markovitch , Dmitry Rusakov, Selective Sampling for Nearest Neighbor Classifiers, Machine Learning, v.54 n.2, p.125-152, February 2004
Volker Tresp, Scaling Kernel-Based Systems to Large Data Sets, Data Mining and Knowledge Discovery, v.5 n.3, p.197-211, July 2001
Charles A. Micchelli , Massimiliano A. Pontil, On Learning Vector-Valued Functions, Neural Computation, v.17 n.1, p.177-204, January 2005
Michael E. Tipping, Sparse bayesian learning and the relevance vector machine, The Journal of Machine Learning Research, 1, p.211-244, 9/1/2001
Balaji Krishnapuram , Lawrence Carin , Mario A. T. Figueiredo , Alexander J. Hartemink, Sparse Multinomial Logistic Regression: Fast Algorithms and Generalization Bounds, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.27 n.6, p.957-968, June 2005
T. Van Gestel , J. A. K. Suykens , G. Lanckriet , A. Lambrechts , B. De Moor , J. Vandewalle, Bayesian framework for least-squares support vector machine classifiers, Gaussian processes, and kernel fisher discriminant analysis, Neural Computation, v.14 n.5, p.1115-1147, May 2002
Zhihua Zhang , James T. Kwok , Dit-Yan Yeung, Model-based transductive learning of the kernel matrix, Machine Learning, v.63 n.1, p.69-101, April     2006
Gavin C. Cawley , Nicola L. C. Talbot, Preventing Over-Fitting during Model Selection via Bayesian Regularisation of the Hyper-Parameters, The Journal of Machine Learning Research, 8, p.841-861, 5/1/2007
Matthias Seeger, Pac-bayesian generalisation error bounds for gaussian process classification, The Journal of Machine Learning Research, 3, p.233-269, 3/1/2003
Arnulf B. A. Graf , Felix A. Wichmann , Heinrich H. Blthoff , Bernhard H. Schlkopf, Classification of Faces in Man and Machine, Neural Computation, v.18 n.1, p.143-165, January 2006
Ralf Herbrich , Thore Graepel , Colin Campbell, Bayes point machines, The Journal of Machine Learning Research, 1, p.245-279, 9/1/2001
Liam Paninski , Jonathan W. Pillow , Eero P. Simoncelli, Maximum Likelihood Estimation of a Stochastic Integrate-and-Fire Neural Encoding Model, Neural Computation, v.16 n.12, p.2533-2561, December 2004
Alexander J. Smola , Bernhard Schlkopf, Bayesian kernel methods, Advanced lectures on machine learning, Springer-Verlag New York, Inc., New York, NY,
gaussian processes;Markov chain Monte Carlo;classification problems;hybrid Monte Carlo;bayesian classification;parameter uncertainty