Upper and Lower Bounds on the Learning Curve for Gaussian Processes.
In this paper we introduce and illustrate non-trivial upper and lower bounds on the learning curves for one-dimensional Guassian Processes. The analysis is carried out emphasising the effects induced on the bounds by the smoothness of the random process described by the Modified Bessel and the Squared Exponential covariance functions. We present an explanation of the early, linearly-decreasing behavior of the learning curves and the bounds as well as a study of the asymptotic behavior of the curves. The effects of the noise level and the lengthscale on the tightness of the bounds are also discussed.
Introduction
A fundamental problem for systems learning from examples is to estimate
the amount of training samples needed to guarantee satisfactory generalisation
capabilities on new data. This is of theoretical interest but also of vital
practical importance; for example, algorithms which learn from data should
not be used in safety-critical systems until a reasonable understanding of
their generalisation capabilities has been obtained. In recent years several
authors have carried out analysis on this issue and the results presented
depend on the theoretical formalisation of the learning problem.
Approaches to the analysis of generalisation include those based on
asymptotic expansions around optimal parameter values (e.g. AIC (Akaike,
1974), NIC (Murata et al., 1994)); the Probably Approximately Correct
convergence approaches (e.g.
Vapnik, 1995); and Bayesian methods.
The PAC and uniform convergence methods are concerned with frequentist-
style confidence intervals derived from randomness introduced with respect
to the distribution of inputs and noise on the target function. A central
concern in these results is to identify the flexibility of the hypothesis class F
to which approximating functions belong, for example, through the Vapnik-Chervonenkis
dimension of F . Note that these bounds are independent
of the input and noise densities, assuming only that the training and test
samples are drawn from the same distribution.
The problem of understanding the generalisation capability of systems
can also be addressed in a Bayesian framework, where the fundamental
assumption concerns the kinds of function our system is required to model.
In other words, from a Bayesian perspective we need to put priors over
target functions. In this context learning curves and their bounds can be
analysed by an average over the probability distribution of the functions. In
this paper we use Gaussian priors over functions which have the advantage
of being more general than simple linear regression priors, but they are
more analytically tractable than priors over functions obtained from neural
networks.
Neal (1996) has shown that for fixed hyperparameters, a large class of
neural network models will converge to Gaussian process priors over functions
in the limit of an infinite number of hidden units. The hyperparameters
of the Bayesian neural network define the parameters of the corresponding
Gaussian Process (GP). Williams (1997) calculated the covariance functions
of GPs corresponding to neural networks with certain weight priors
and transfer functions.
The investigation of GP predictors is motivated by the results of Rasmussen
(1996), who compared the performances obtained by GPs to those obtained
by Bayesian neural networks on a range of tasks. He concluded that GPs
were at least as good as neural networks. Although the present study deals
with regression problems, GPs have also been applied to classification problems
(e.g. Barber and Williams, 1997).
In this paper we are mainly concerned with the analysis of upper and
lower bounds on the learning curve of GPs. A plot of the expected generalisation
error against the number of training samples n is known as a learning
curve. There are many results available concerning leaning curves under different
theoretical scenarios. However, many of these are concerned with the
asymptotic behaviour of these curves, which is not usually of great practical
importance as it is unlikely that we will have enough data to reach the
asymptotic regime. Our main goal is to explain some of the early behaviour
of learning curves for Gaussian processes.
The structure of the paper is as follows. GPs for regression problems
are introduced in Section 2. As will be shown, the whole theory of GPs is
based on the choice of the prior covariance function C p (x; x 0 in Section
3 we present the covariance functions we have been using in this study.
In Section 4 the learning curve of a GP is introduced. We present some
properties of the learning curve of GPs as well as some problems may arise
in evaluating it. Upper and lower bounds on the learning curve of a GP in a
non-asymptotic regime are presented in Section 5. These bounds have been
derived from two different approaches: one makes use of main properties of
the generalisation error, whereas the other is derived from an eigenfunction
decomposition of the covariance function. The asymptotic behaviour of the
upper bounds is also discussed.
A set of experiments have been run in order to assess the upper and lower
bounds of the learning curve. In Section 6 we present the results obtained
and investigate the link between tightness of the bounds and the smoothness
of the stochastic process modelled by a GP. A summary of the results and
some open questions are presented in the last Section.
Gaussian Processes
A collection of random variables fY (x) jx 2 Xg indexed by a set X defines a
stochastic process. In general the domain X might be R d for some dimension
d although it could be even more general. A joint distribution characterising
the statistics of the random variables gives a complete description of the
stochastic process.
A GP is a stochastic process whose joint distribution is Gaussian; it is
fully defined by giving a Gaussian prior distribution for every finite subset
of variables.
In the following we concentrate to the regression problem assuming that
the value of the target function t (x) is generated from an underlying function
y (x) corrupted by Gaussian noise with mean 0 and variance oe 2
. Given a
collection of n training data D
(where each t i is the
observed output value at the input point x i ), we would like to determine
the posterior probability distribution p (yjx; D n ).
In order to set up a statistical model of the stochastic process, the set
of n random variables
modelling the function values
at respectively, is introduced. Similarly t is the collection of
target values
denote the set of training inputs
We also denote with ~ y the vector whose components are
y and the test value y at the point x. The distribution p (~ yjx; D n ) can
be inferred using Bayes' theorem. In order to do so, we need to specify a
prior over functions as well as evaluate the likelihood of the model and the
evidence for the data.
A choice for a prior distribution of the stochastic vector ~ y is a Gaussian
prior distribution:
\Gamma2
~
y
This is a prior as it describes the distribution of the true underlying values
without any reference to the target values t. The covariance matrix \Sigma can
be partitioned as
The element (K
is the covariance between the i-th and the j-th training
points, i.e. (K
y
y
. The components
of the vector k (x) are the covariances of the test point with all the
training data
is the covariance of the test
point with itself.
A GP is fully specified by its mean E [y covariance function
Below we set -
this is a valid assumption provided that any known offset or trend in the
data has been removed. We can also deal with - (x) 6= 0, but this introduces
some extra notational complexity. A discussion about the possible choices
of the covariance function C p (x; x 0 ) is given in Section 3. For the moment
we note that the covariance function is assumed to depend upon the input
variables (x; x 0 ). Thus the correlation between function values depends
upon the spatial position of the input vectors; usually this will be chosen so
that the closer the input vectors, the higher the correlation of the function
values.
The likelihood relates the underlying values of the function to the target
data. Assuming a Gaussian noise corrupting the data, we can write the
likelihood as
I. The likelihood refers to the stochastic variables representing
the data; so t; y 2 R n
and\Omega is an n \Theta n matrix.
Given the prior distribution over the values of the function p
Bayes' rule specifies the distribution p (~ yjx; D n ) in terms of the likelihood
of the model p (tjy) and the evidence of the data p (D n ) as
Given such assumptions, it is a standard result (e.g. Whittle, 1963) to derive
the analytic form of the predictive distribution marginalising over y. The
predictive distribution turns out to be y (x) - N
where the
mean and the variance of the Gaussian function are
The most probable value -
y (x) is regarded as the prediction of the GP on the
test point x; K is the covariance matrix of the targets t:
I. The
estimate of the variance oe 2
(x) of the posterior distribution is considered
as the error bar of -
y (x). In the following, we always omit the subscript - y in
, taking it as understood. Since the estimate 1 is a linear combination
of the training targets, GPs are regarded as linear smoother (Hastie and
Tibshirani, 1990).
3 Covariance functions
The choice of the covariance function is a crucial one. The properties of
two GPs, which differ only in the choice of the covariance function, can be
remarkably diverse. This is due to the r-ole of the covariance function which
has to incorporate in the statistical model the prior belief about the underlying
function. In other words the covariance function is the analytical
expression of the prior knowledge about the function being modelled. A misspecified
covariance function affects the model inference as it has influence
on the evaluation of Equations 1 and 2.
Formally every function which produces a symmetric, positive semi-definite
covariance matrix K for any set of the input space X can be chosen
as covariance function. From an applicative point of view we are interested
only in functions which contain information about the structure of the underlying
process being modelled.
The choice of the covariance function is linked to the a priori knowledge
about the smoothness of the function y (x) through the connection between
the differentiability of the covariance function and the mean-square differentiability
of the process. The relation between smoothness of a process
and its covariance function is given by the following theorem (see e.g. Adler,
exists and is finite at (x; x), then the stochastic
process y (x) is mean square differentiable in the i-th Cartesian direction at
x. This theorem is relevant as it links the differentiability properties of the
covariance function with the smoothness of the random process and justifies
the choice of a covariance function depending upon the prior belief about
the degree of smoothness of y (x).
In this work we are mainly concerned with stationary covariance func-
tions. A stationary covariance function is translation invariant (i.e. C p (x; x
depends only upon the distance between two data points.
In the following, the covariance functions we have been using are presented.
In order to simplify the notation, we consider the case
The stationary covariance function squared exponential (SE) is defined
as
where - is the lengthscale of the process. The parameter - defines the
characteristic length of the process, estimating the distance in the input
space in which the function y (x) is expected to vary significantly. A large
value of - indicates that the function is almost constant over the input space,
whereas a small value of the lengthscale designates a function which varies
rapidly. The graph of this covariance function is shown by the continuous
line in Figure 1. As the SE function has infinitely many derivatives it gives
rise to smooth random processes (y (x) posses mean-square differentiability
up to order 1).
It is possible to tune the differentiability of a process, introducing the
modified Bessel covariance function of order k (MB k ). It is defined as
a i
exp
where K - (\Delta) is the modified Bessel function of order - (see e.g. Equation
8:468 in Gradshteyn and Ryzhik, 1993), with
Below we set the constant - such that C p 1. The factors a k are
constants depending on the order - of the Bessel function. Mat'ern (1980)
shows that the functions MB k define a proper covariance. Stein (1989) also
noted that the process with covariance function MB k is
differentiable.
In this study we deal with modified Bessel covariance function of orders
We note that MB 1 corresponds to the Ornstein-Uhlenbeck covariance function
which describes a process which is not mean square differentiable.
If k !1, the MB k behaves like the SE covariance function; this can be
easily shown by considering the power spectra of MB k and SE which are
and S se (!) / - exp
Since
lim
the MB k behaves like SE for large k, provided that - is rescaled accordingly.
Modified Bessel covariance functions are also interesting because they
describe Markov processes of order k. Ihara (1991) defines Y (x) to be a
strict sense Markov process of order k if it is differentiable
at every x 2 R and if P (Y
states that a Gaussian process is a
1 Note that the definition of a Markov process in discrete and continuous time is rather
different. In discrete time, a Markov process of order k depends only on the previous
k times, but in continuous time the dependence is on the derivatives at the last time.
However, function values at previous times clearly allow approximate computation of
Markov process of order k in the strict sense if and only if it is an autoregressive
model of order k (AR(k)) with a power spectrum (in the Fourier
domain) of the form
Y
As the power spectrum of MB k has the same form of the power spectrum of
an AR(k) model, the stochastic process whose covariance function is MB k
is a strict sense k-ple Markov process. This characteristic of the MB k covariance
functions is important as it ultimately affects the evaluation of the
generalisation error (as we shall see in Section 6).

Figure

2 shows the graphs of four (discretised) random functions generated
using the MB k covariance functions (with and the SE func-
tion. We note how the smoothness of the random function specified is dependent
of the choice of the covariance function. In particular, the roughest
function is generated by the Ornstein-Uhlenbeck covariance function (Figure
whereas the smoothest one is produced by the SE (Figure 2(d)). An intermediate
level of regularity characterises the functions of Figures 2(b) and
2(c), corresponding to MB 2 and MB 3 respectively. Note that the number
of zero-level upcrossings in [0; 1] (denoted N u ) is only weakly dependent on
the order of the process. For MB 2 and MB 3 E[N u
derivatives (e.g. via finite differences) and thus one would expect that in the continuous-time
situation the previous k process values will contain most of the information needed
for prediction at the next time. Note that for the Ornstein-Uhlenbeck process Y
depends only on the previous observation Y (t).
and (
respectively (see Papoulis (1991) eqn 16-7 for details). For
the SE process E[N u As the Ornstein-Uhlenbeck process is
non-differentiable, the formula given for E[N u ] cannot be applied in this
case.
Learning curve for Gaussian processes
A learning curve of a model is a function which relates the generalisation
error to the amount of training data; it is independent of the test points as
well as the locations of the training data and depends only upon the amount
of data in the training set. The learning curve for a GP is evaluated from
the estimation of the generalisation error averaged over the distribution of
the training and test data.
For regression problems, a measure of the generalisation capabilities of
a GP is the squared difference E g
between the target value on a test
point x and the prediction made by using Equation 1:
The Bayesian generalisation error at a point x is defined as the expectation
of
Dn (x; t) over the actual distribution of the stochastic process t:
. Under the assumption that the data set is actually
generated from a GP, it is possible to read Equation 2 as the Bayesian
generalisation error at x given training data D n . To see this, let us consider
the (n 1)-dimensional distribution of the target values at x 1
x. This is a zero-mean multivariate Gaussian. The prediction at the test
point x is -
I. Hence the expected
generalisation error at x is given by
\Theta

\Theta
\Theta tt T

where we have used
\Theta tt T
K. Equation 5 is identical
to oe 2
(x) as given in Equation 2 with the addition of the noise variance oe 2
(since we are dealing with noisy data). The variance of
can also be calculated (Vivarelli, 1998).
The covariance matrix pertinent for these calculations is the true prior;
if a GP predictor with a different (incorrect) covariance function is used, the
expression for the generalisation error becomes
c
where the indices c and i denote the correct and incorrect covariance functions
respectively. It can be shown (Vivarelli, 1998) that this is always larger
than Equation 5.
Another property of the generalisation error can be derived from the
following observation: adding more data points never increases the size of
the error bars on prediction (oe 2
n (x)). This can be proved using
standard results on the conditioning of a multivariate Gaussian (see
Vivarelli, 1998). It can also be understood by the information theoretic
argument that conditioning on additional variables never increases the entropy
of a random variable. Considering t (x) to be the random variable,
we observe that its distribution is Gaussian, with variance independent of
t (although the mean does depend on t). The entropy of a Gaussian is2 log
. As log is monotonic, the assertion is proved. This argument
is an extension of that in (Qazaz et al., 1997), where the inequality
was derived for generalized linear regression.
Dn (x), a similar inequality applies also to the Bayesian
generalisation errors and hence
This remark will be applied in Section 5 for evaluating upper bounds on the
learning curve.
Equation 5 calculates the generalisation error at a point x. Averaging
Dn (x) over the density distribution of the test points p (x), the expected
generalisation error E
Dn is
For particular choices of p (x) and C p (x) the computation of this expression
can be reduced to a n \Theta n matrix computation as E x
\Theta
\Theta k (x) k T (x)
. We also note that Equation 7 is independent
of the test point x but still depends upon the choice of the training data
D n . In order to obtain a proper learning curve for GP, E g
Dn needs to be
averaged 2 over the possible choices of the training data D n . However, it is
very difficult to obtain the analytical form of E g for a GP as a function of
n. Because of the presence of the k T Equation 5, the
matrix K and vector k (x) depend on the location of the training points:
the calculations of the averages with respect to the data points seems very
hard. This motivates looking for upper and lower bounds on the learning
curve for GP.
5 Bounds on the learning curve
For the noiseless case, a lower bound on the generalisation error after n
observations is due to Michelli and Wahba (1981). Let be the
ordered eigenvalues of the covariance function on some domain of the input
space X . They showed that E g (n) -
a bound on the learning curve for the noisy case; since the bound uses
observations consisting of projections of the random function onto the first
eigenfunctions, it is not expected that it will be tight for observations
which consist of function evaluations.
Other results that we are aware of pertain to asymptotic properties of
(n). Ritter (1996) has shown that for an optimal sampling of the input
space, the asymptotics of the generalisation error is O
Hansen (1993) showed that for linear regression models it is possible to average over
the distribution of the training sets.
a random process which obeys to the Sacks-Ylvisaker 3 conditions of order s
(see Ritter et al., 1995 for more details on Sacks-Ylvisaker conditions). In
general, the Sacks-Ylvisaker order of the MB k covariance function is
1. For example an MB 1 process has hence the generalisation error
shows a n \Gamma1=2 asymptotic decay. In the case that X ae R, the asymptotically
optimal design of the input space is the uniform grid.
Silverman (1985) proved a similar result for random designs. Haussler
and Opper (1997) have developed general (asymptotic) bounds for the expected
log-likelihood of a test point after seeing n training points.
In the following we introduce upper and lower bounds on the learning
curve of a GP in a non-asymptotic regime. An upper bound is particularly
useful in practice as it provides an (over)estimate of the number of
examples needed to give a certain level of performance. A lower bound is
similarly important because it contributes to fix the limit which can not be
outperformed by the model.
The bounds presented are derived from two different approaches. The
first approach makes use of the particular form assumed by the generalisation
error at x (E g
(x)). As the error bar generated by one data point
is greater than that generated by n data points, the former can be considered
as an upper bound of the latter. Since this observation holds for the variance
due to each one the data points, the envelope of the surfaces generated by
Loosely speaking, a stochastic process possessing s mean-square derivatives but not
is said to satisfy the Sacks-Ylvisaker conditions of order s.
the variances due to each data point is also an upper bound of oe 2
n (x). In
particular as oe 2
Dn (x) (cf. Equation 5), the envelope is an upper
bound of the generalisation error of the GP. Following this argument, we
can assert that an upper bound on E g
Dn (x) is the one generated by every
GP trained with a subset of D n . The larger the subset of D n the tighter the
bound.
The two upper bounds we present differ in the number of training points
considered in the evaluation of the covariance: the derivation of the one-point
upper bound E u
1 (n) and the two-point upper bound E u
2 (n) are presented
in Section 5.1 and Section 5.2 respectively. Section 5.3 reports the
asymptotic expansion of E u
1 (n) in terms of - and oe 2
- .
The second approach is based on the expansion of the stochastic process
in terms of the eigenfunctions of the covariance function. Within this
framework, Opper proposed bounds on the training and generalisation error
(Opper and Vivarelli, 1999) in terms of the eigenvalues of C p (x; x 0 ); the
lower bound E l (n) obtained is presented in Section 5.4.
In order to have tractable analytical expressions, all the bounds have
been derived by introducing three assumptions:
i The input space X is restricted to the interval [0; 1];
ii The probability density distribution of the input points is uniform:
iii The prior covariance function C p (x; x 0 ) is stationary.
5.1 The one-point upper bound E u
For the derivation of the one-point upper bound, let us consider the error
bar generated by one data point x i . Since C
Equation 2 becomes
For x far away from the training point x i , oe 2
the confidence on
the prediction for a test point lying far apart from the data point x i is quite
low as the error bar is large. The closer x to x i , the smaller the error bar on
y (x). When
Irrespective
of the value of C p (0), r varies from 0 to 1. As normally C p (0) AE oe 2
and thus oe 2
- . So far we have not used any hypothesis concerning
the dimension of the variable x, thus this observation holds regardless the
dimension of the input space.
The effect of just one data point helps in introducing the first upper
bound. The interval [0; 1] is split up in n subintervals
\Theta
a

(where a
=2 and b
centred around the i-th
data point x i , with a
Let us consider the i-th training point and the error bar oe 2
by x i . When x 2
\Theta
a
1 this relation is illustrated in Figure
3, where the envelope of the surfaces of the errors due to each datapoint
(denoted by E g
(x)) is an upper bound of the overall generalisation error.
Since we are dealing with positive functions, an upper bound of the expected
generalisation error on the interval
\Theta
a
can be written as
a i
a i
where p (x) is the distribution of the test points. Summing up the contributions
coming from each training datapoint in both sides of Equation 8 and
setting
a i
a i
The interval where the contribution of the variance due to x i contributes to
Equation 8 is also shown in Figure 3.
Under the assumption of the stationarity of the covariance function,
integrals such as those in the right hand side of Equation 9 depend only
upon differences of adjacent training points (i.e. x
The right hand side of Equation 9 can be rewritten as
a i
a i
dx
I
I
where
I
Equation 11 can be derived changing the variables in the two integrals of
Equation Equation 11 is
an upper bound on E g
and still depends upon the choice of the training
data D n through the interval of integration. We note that the arguments
of the integrals I (\Delta) in Equation 11 are the differences between adjacent
training points. Denoting those differences with , we can model
their probability density distribution by using the theory of order statistics
(David, 1970). Given an uniform distribution of n training data over the
interval [0; 1], the density distribution of the differences between adjacent
points is p . Since this is true for all the differences ! i
we can omit the superscript i and thus the expectation of the integrals in
Equation 11 over p (!) is
I
I
I (! n )
. Both the integrals
can be calculated following a similar procedure. Let us consider
where the second line has been obtained integrating by parts. The last line
follows from the fact that [I (!)
We are now able to write an upper bound on the learning curve as
The calculations of the integrals in the above expression are straightforward
though they involve the evaluation of hyper-geometric functions (because of
the As the evaluation of such functions is computationally
intensive, we found preferable to evaluate Equation 14 numerically.
5.2 The two-points upper bound E u(n)
The second bound we introduce is the natural extension of the previous idea,
using two data points rather than one. By construction, we expect that it
will be tighter than the one introduced in Section 5.1.
Let us consider two adjacent data points x i and x i+1 of the interval [0; 1],
with By the same argument presented in the previous section,
the following inequality holds:
2 (x) is the variance on the prediction -
y (x) generated by the data
points x i and x i+1 . Similarly to Equation 9, summing up the contributions
of both sides of Equation 15 we get an upper bound on the generalisation
error:
where we have defined
After some calculations (see Appendix A) we obtain
where
I 1
(!). The calculation of the integrals with respect to !
in E u
2 (n) is complicated by the determinant \Delta (!) in the denominator and by
the distribution n so we preferred to evaluate them numerically
as we did for E u
5.3 Asymptotics of the upper bounds
From Equation 14, an expansion of E u
1 (n) in terms of - and oe 2
- in the limit
of a large amount of training data can be obtained. The expansion depends
upon the covariance function we are dealing with. Expanding the covariance
function around 0, the asymptotic form of E u
1 (n) for MB 1 is
n-
whereas for the functions MB 2 , MB 3 and SE it is
The asymptotic value of E u
depends neither on the lengthscale of
the process nor on the order of the covariance function MB k for k - 1 but
is a function of the ratio r:
lim
As we pointed out in Section 5.1, this is the minimum generalisation error
achievable by a GP when it is trained with just one datapoint. The n !1
scenario corresponds to the situation in which every test point is close to
a datapoints. As mentioned at the beginning of this Section, the asymptotics
of the learning curve for the MB k and SE covariance functions are
O
\Delta and O
respectively. Although the expansions of
decay asymptotically faster than the learning curves, they reach an
asymptotic plateau oe 2
- . We also note that the asymptotic values
get closer to the true noise level when r - 1, i.e. for the unrealistic
case oe 2
The smoothness of the process enters into the asymptotics through a
factor O
This factor affects the rate of approach to the asymptotic value oe 2
of E u
1 (n). We notice that larger lengthscales and noise levels increase the
rate of decay of E u
1 (n) to the asymptotic plateau.
The asymptotic form of E u
2 (n) for the MB 1 , MB 2 , MB 3 and SE covariance
functions is (Vivarelli, 1998)
a
where the value of a depends upon the choice of the covariance function and
(0). Similarly to the expansion of E u
1 (n), the decay rate of
2 (n) is faster than the asymptotic decay of the actual learning curves but
it reaches an asymptotic plateau of
lim
It is straightforward to verify that the asymptotic plateau of E u
2 (n) is lower
than the one of E u
1 (n) and that it corresponds to the error bar estimated
by a GP with two observations located at the test point.
5.4 The lower bound E l (n)
Opper (Opper and Vivarelli, 1999) proposed a bound on the learning curve
and on the training error based on the decomposition of the stochastic process
y (x) in terms of the eigenfunctions of the covariance C p (x; x 0 ).
Denoting with ' k set of functions satisfying
the integral equation
Z
the Bayesian generalisation error E
(where
y (x) is the true underlying stochastic function and - y (x) is the GP predic-
tion) can be written in terms of the eigenvalues of C p (x; x 0 ). In particular,
after an average over the distribution of the input data, E g (D n ) can be
written as E g (D n

, where   is the infinite dimension
diagonal matrix of the eigenvalues and V is a matrix depending on
the training data, i.e. V
By using Jensen's inequality, it is possible to show that a lower bound of
the learning curve and an upper bound of the training error is (Opper and
In this paper we mean to compare this lower bound to the actual learning
curve of a GP. As our bounds are on t rather than y, we must add oe 2
- to the
expression obtained in Equation 23 giving an actual lower bound of
6 Results
As we pointed out in Section 4, the analytic calculation of the learning curve
of a GP is infeasible. Since the generalisation error
is a complicated function of the training data (which are inside the elements
of k (x) and K \Gamma1 ), it is problematic to perform an integration over the
distribution of the training points. For comparing the learning curve of
the GP with the bounds we found, we need to evaluate the expectation
of the integral in Equation 25 over the distribution of the data: E
EDn
\Theta
Dn
. An estimate of E g (n) can be obtained using a Monte Carlo
approximation of the expectation. We used 50 generations of training data,
sampling uniformly the input space [0; 1]. For each generation, the expected
generalisation error for a GP has been evaluated using up to 1000 datapoints.
Using the 50 generations of training data, we can obtain an estimate of the
learning curve E g (n) and its 95% confidence interval.
Since this study is focused on the behaviour of bounds on learning
curve on GP, we assume the true values of the parameters of the GP are
known. So we chose the value of the constant - for the covariance functions
Equation 4) such that C p
allowed the lengthscale - and the noise level oe 2
- to assume several values
To begin with, we study how the smoothness of a process affects the
behaviour of the learning curve. The empirical learning curves of Figure
4 have been obtained for processes whose covariance functions are MB 1 ,
0:1. We can notice that all the
learning curves exhibit an initial linear decrease. This can be explained
considering that without any training data, the generalisation error is the
maximum allowable by the model (C
- ). The introduction
of a training point x 1 creates a hole on the error surface: the volume of
the hole is proportional to the value of the lengthscale and depends on the
covariance function. The addition of a new data point x 2 will have the effect
of generating a new hole in the surface. With such a few data points it is
likely that the two data lie down far apart one from the other, giving rise
to two distinct holes. Thus the effect that a small dataset exerts to pull
down the error surface is proportional to the amount of training points and
explains the initial linear trend.
Concerning the asymptotic behaviour of the learning curves, we have
verified that they agree with the theoretical analysis carried out by Ritter
(1996). In particular, a log-log plot of the learning curves with a MB k
covariance function shows an asymptotic behaviour as O
. A
similar remark applies to the SE covariance function, with an asymptotic
decay rate of O
(Opper, 1997). We have also noted that the
smoother the process described by the covariance function the smaller the
the amount of training data needed to reach the asymptotic regime.
The behaviour of the learning curves is affected also by the value of the
lengthscale of the process and by the noise level and this is illustrated in

Figure

7. The learning curves shown in Figure 5(a) have been obtained for
the MB 1 covariance function setting the noise level oe 2
0:1 and varying the
values of the parameters Intuitively, Figure 5(a) suggests
that decreasing the lengthscale stretches the early behaviour of the learning
curve and the approach to the asymptotic plateau lasts longer; this is due
to the effect induced by different values of the lengthscale which stretch or
compress the input space. We have verified that rescaling the amount of
data n by the ratio of the two lengthscales, the two curves of Figure 5(a)
lay on top of each other.
The variation of the noise level shifts the learning curves from the prior
value C p (0) by an offset equal to the noise level itself (cf. Equation 5);
in order to see any significant effect of the noise on the learning curve,

Figure

5(b) shows a log-log graph of E
obtained for a stochastic
process with MB 3 covariance function, setting
. We can notice two main effects. The noise variance affects
the actual values of the generalisation error since the learning curve obtained
with high noise level is always above the one obtained with a low noise
level. A second effect concerns the amount of data necessary to reach the
asymptotic regime. The learning curve characterised by an high noise level
needs fewer datapoints to attain to the asymptotic regime.
Stochastic processes with different covariance functions and different values
of lengthscales and noise variance behave in a similar way.
In the following we discuss the results in two main subsections: results
about the bounds E u
2 (n) are presented in Section 6.1, whereas
the lower bound of Section 5.4 is shown in Section 6.2. As the results we
obtained for these experiments show common characteristics, we show the
bounds of the learning curve obtained by setting
6.1 The upper bounds E u(n) and E u(n)
Each graph in Figure 6 shows the empirical learning curve with its confidence
interval and the two upper bounds E u
(n). The curves are shown
for the MB 1 , MB 2 , MB 3 and the SE covariance functions.
For a limited amount of training data it is possible to notice that the upper
error bar associated to EDn [E g (n)] lies above the actual upper bounds.
This effect is due to the variability of the generalisation error for small data
sets and suggests that the bounds are quite tight for small n. The effect
disappears for large n, when the estimate of the generalisation error is less
sensitive to the composition of the training set.
As expected, the two-point upper bound E u
2 (n) is tighter than the one-point
upper bound E u
We note that the tightness of the upper bound depends upon the covariance
function, being tighter for rougher processes (such as MB 1 ) and getting
worse for smoother processes. This can be explained by recalling that covariance
functions such as the MB k correspond to Markov processes of order k
(cf. Section 3). Although the Markov process is actually hidden by the presence
of the noise, E g (n) is still more dependent on training data lying close
to the test point x than on more distant points. Since the bounds E u
calculated by using only local information (namely the
closest datapoint to the test point, or the closest datapoints to the left and
right, respectively), it is natural that the more the variance at x depends on
local data points, the tighter the bounds become.
For instance, let us consider MB 1 , the covariance function of a first order
Markov process. For the noise-free process, knowledge of data-points lying
beyond the the left and right neighbours of x does not reduce the generalisation
error at x 4 . Although in the noisy case more distant data-points
4 This is because the process values at the training points and test point form a Markov
chain, and knowledge of the process values to the left and right of the test point "blocks"
reduce the generalisation error (because of the term oe 2
- in the covariance
matrix K), it is likely that local information is still the most important.
The bounds on the learning curves computed for MB 2 and MB 3 confirm
this remark, as they are looser than for MB 1 . For the SE covariance function,
this effect still holds and is actually enlarged.
In Section 5.3 we have shown that the asymptotic behaviour of the bound
depends on the covariance function, being O
O
plots of the upper bounds confirm the
analysis carried out in Section 5.3, where we showed that E u
approach asymptotic plateaux. In particular, E u
tends to oe 2
O
tends to
The quality of the bounds for processes characterised by different length-
scales and different noise levels are comparable to the ones described so far:
the tightness of E u
still depend on the smoothness of the
process. As explained at the beginning of this section, a variation of the
lengthscale has the same effect of a rescaling in the number of training data.
This can be observed explicitly in the asymptotic analysis of Equations
and 19, where the decay rate depends on the factor n-.
For a fixed covariance function, we note that the bounds are tighter for
lower noise variance; this is due to the fact that the lower the noise level the
better the hidden Markov process manifests itself. For smaller noise levels
the influence of more remote observations.
the learning curve becomes closer to the bounds because the generalisation
error relies on the local behaviour of the processes around the test data; on
the contrary, a larger noise level hides the underlying Markov Process thus
loosening the bounds.
6.2 The bound E l (n)
We have also run experiments computing the lower bound we obtained from
Equation 24 for processes generated by the covariance priors MB 1 , MB 2 ,
MB 3 and SE .
Equation 24 shows that the evaluation of E l (n) involves the computation
of an infinite sum of terms; we truncated the series considering only those
terms which add a significant contribution to the sums, i.e. j k =oe 2
" is the machine precision. Since each contribution in the series is positive,
the quantity computed is still a lower bound of the learning curve.

Figure

7 shows the results of the experiment in which we set
0:1. The graphs of the lower bound lies below the empirical learning
curve, being tighter for large amount of data; in particular for the smoothest
processes with large amount of data, the 95% confidence intervals lay below
the actual lower bound.
For the lower bound tends to the noise level oe 2
- . As with the
empirical learning curve, log-log plots of E l
y (n) show an asymptotic decay to
zero as O(n \Gamma(2k\Gamma1)=2k ) and O
\Delta for the MB k and the SE covariance
functions, respectively.
The graphs of Figure 7 show also that the tightness of the bound depends
on the smoothness of the stochastic process; in particular smooth processes
are characterised by a tight lower bound on the learning curve E g (n). This
can be explained by observing that E l (n) is a lower bound on the learning
curve and an upper bound of the training error. The values of smooth
functions do not have large variation between training points and thus the
model can infer better on test data; this reduces the generalisation error
pulling it closer to the training error. Since the two errors sandwich the
bound of Equation 24, E l (n) becomes tight for smooth processes.
We can also notice that the tightness of the lower bound depends on the
noise level, becoming tight for high the noise level and loose for small noise
level. This is consistent with a general characteristic of E l (n) which is monotonically
decreasing function of the noise variance (Opper and Vivarelli,
1999).
In this paper we have presented non-asymptotic upper and lower bounds for
the learning curve of GPs. The theoretical analysis has been carried out for
one-dimensional GPs characterised by several covariance functions and has
been supported by numerical simulations.
Starting from the observation that increasing the amount of training
data never worsens the Bayesian generalisation error, an upper bound on
the learning curve can be estimated as the generalisation error of a GP
trained with a reduced dataset. This means that for a given training set the
envelope of the generalisation errors generated by one and two datapoints
is an upper bound of the actual learning curve of the GP. Since the expectation
of the generalisation error over the distribution of the training data is
not analytically tractable, we introduced the two upper bounds E u
1 (n) and
2 (n) which are amenable to average over the distribution of the test and
training points. In this study we have evaluated the expected value of the
future directions of research should also deal with the evaluation of
the variances.
In order to highlight the behaviour of the bounds with respect to the
smoothness of the stochastic process, we investigated the bounds for the
modified Bessel covariance function of order k (describing stochastic processes
differentiable) and the squared exponential
function (describing processes mean square-differentiable up to the order
1).
The experimental results have shown that the learning curves and their
bounds are characterised by an early, linearly decreasing behaviour; this is
due to the effect exerted by each datapoint in pulling down the surface of
the prior generalisation error. We also noticed that the tightness of the
bounds depends on the smoothness of the stochastic processes. This is due
to the facts that the bounds rely on subsets of the training data (i.e. one
or two datapoints) and the modified Bessel covariance functions describe
Markov processes of order k; although in our simulations the Markovian
processes were hidden by noise, the learning curves depend mainly on local
information and our bounds become tighter for rougher processes.
We also investigated the behaviour of the curves with respect to the
variation of the correlation lengthscale of the process and the variance of
the noise corrupting the stochastic process. We noticed that the lengthscale
stretches the behaviour of the curves effectively rescaling the number of
training data. As the noise level has the effect of hiding the underlying
Markov process, the upper bounds become tighter for smaller noise variance.
The expansion of the bounds in the limit of large amount of data highlights
an asymptotic behaviour depending upon the covariance function;
approaches the asymptotic plateau as O
(for the MB 1 covariance
and as O
for smoother processes; the rate of decay
to the plateau of E u
2 (n) is O
. Numerical simulations supported our
analysis.
One limitation of our analysis is the dimension of the input space; the
bounds have been made analytically tractable by using order statistics results
after splitting up the one dimensional input space of the GP. In higher
dimensional spaces the partition of the input space can be replaced by a
Voronoi tessellation that depends on the data D n but averaging over this
distribution appears to be difficult. One can suggest an approximate evaluation
of the upper bounds by an integration over a ball whose radius depends
upon the number of examples and the volume of the input space in which
the bound holds. In any case we expect that the effect due to larger input
dimension is to loosen the upper bounds. We note that recent work by (Sol-
lich, 1999) has derived some good approximations to the learning curve, and
that his methods apply in more than one dimension 5 .
We also ran some experiments by using the lower bound proposed by
Opper, based on the knowledge of the eigenvalues of the covariance function
of the process. Since the bound E l (n) is also an upper bound on the training
error, we observed that the bound is tighter for smooth processes, when the
learning curve becomes closer to the training error. Also the noise can vary
the tightness of E l (n); a low noise level loosens the lower bound. Unlike the
upper bounds, the lower bound can be applied also in multivariate problems,
as it is easily extended to high dimension input space; however it has been
verified (Opper and Vivarelli, 1999) that the bound becomes less tight in
input space of higher dimension.


Appendix

A: The two-points upper bound E u
In this Appendix we derive Equation 17 starting from Equation 16.
We start by calculating oe 2
(x). As the covariance matrix generated by
two data points is a 2 \Theta 2 matrix, it is straightforward to evaluate oe 2
Considering the two training data x i and x i+1 , the covariance matrix of the
5 The reference to Sollich (1999) was added when the manuscript was revised in April
1999.
GP is
From the evaluation of the determinant of K as
As the covariance vector for the test point x is k
the variance assumes the form
Changing variables in the covariances C p
(as
turns out that the upper bound
generated by oe 2
2 (x) in the interval
\Theta
(when i 6= 0; n), is
I 1
where
I 1
(-) d- and I 2
It is noticeable that, similarly to Equation 11, also the integrals I 1 (\Delta), I 2 (\Delta)
and the determinant \Delta
depend upon the length of the interval
of integration We evaluate the contributions to the upper
bound over the intervals
\Theta 0; x 1
and [x n ; 1] by integrating the variance oe 2
generated by x 1 and x n over
\Theta 0; x 1
and [x n ; 1] respectively. Hence the right
hand side of Equation 16 can be rewritten as
I 1
I
where I (\Delta) is defined in Equation 12.
Equation 26 is still dependent on the distribution of the training data
because it is a function of the distances between adjacent training points
. Similarly to Equation 11, we obtain an upper bound independent of
the training data by integrating Equation 13 over the distribution of the
differences
\GammaC

Acknowledgments

This research forms part of the "Validation and Verification of Neural Net-work
Systems" project funded jointly by EPSRC (GR/K 51792) and British
Aerospace. We thank Dr. Manfred Opper, and Dr. Andy Wright of BAe for
helpful discussions. We also thank the anonymous referees for their comments
which have helped improve this paper. F. V. was supported by a
studentship from British Aerospace.


--R

The Geometry of Random Fields.
A new look at statistical model identification.
Gaussian processes for Bayesian classification via hybrid Monte Carlo.
Order Statistics.
Table of Integrals
Stochastic linear learning: Exact test and training error averages.
Generalized Additive Models.
Mutual information
Information Theory.
Design problems for optimal surface interpolation.
Network information criterion-determining the number of hidden units for artificial neural network models
Bayesian Learning for Neural Networks.
Lecture Notes in Statistics 118.
Regression with gaussian processes: average case per- formance
General bounds on Bayes errors for regression with Gaussian Processes.


An Upper Bound on the Bayesian Error Bars for Generalized Linear Regression.
Evaluation of Gaussian Processes and Other Methods for Non-linear Regression
Almost optimal differentiation using noisy data.
Multivariate integration and approximation for random fields satisfying Sacks- Ylvisaker conditions
Some aspects of the spline smoothing approach to non-parametric regression curve filtering
Learning Curves for Gaussian Processes.
A theory of the learnable.
The Nature of Statistical Learning Theory.
Studies on generalisation in Gaussian processes and Bayesian neural networks.
Prediction and regulation by linear least square methods.
Computing with infinite networks.




Figures 6(a)
Figure 7: Figures 7(a)
--TR

--CTR
Peter Sollich , Anason Halees, Learning curves for Gaussian process regression: approximations and bounds, Neural Computation, v.14 n.6, p.1393-1428, June 2002
gaussian processes;bounds;generalisation error;bayesian inference