Variational Model Selection for Sparse GaussianProcess Regression
Michalis K. Titsias
School of Computer Science,
University of Manchester, UKmtitsias@cs.man.ac.uk

Abstract
Sparse Gaussian process methods that use inducing variables require the selection of the inducing
inputs and the kernel hyperparameters. We introduce a variational formulation for sparse approximations that jointly infers the inducing inputs and the kernel hyperparameters by maximizing a lower
bound of the true log marginal likelihood. The key property of this formulation is that the inducing
inputs are defined to be variational parameters which are selected by minimizing the Kullback-Leibler
divergence between the variational distribution and the exact posterior distribution over the latentfunction values. We apply this technique to regression and we compare it with other approaches in
the literature.

1

Introduction

Gaussian processes (GPs) are stochastic processes that, in the context of Bayesian statistics, can be used
as non-parametric priors over real-valued functions that can be combined with data to give posterior
processes over these functions (O’Hagan, 1978; Wahba, 1990). In machine learning GPs offer a Bayesiankernel-based framework for solving supervised learning tasks such as regression and classification; see e.g.
(Rasmussen and Williams, 2006).
However, the application of GP models is intractable for large datasets because the time complexity
scales as O(n3 ) and the storage as O(n2 ) where n is the number of training examples. To overcome
this limitation, many approximate or sparse methods have been proposed in the literature (Williams andSeeger, 2001; Smola and Bartlett, 2001; Csato and Opper, 2002; Lawrence et al., 2002; Seeger et al.,
2003; Schwaighofer and Tresp, 2003; Keerthi and Chu, 2006; Snelson and Ghahramani, 2006; QuiñoneroCandela and Rasmussen, 2005; Walder et al., 2008). These methods construct an approximation basedon a small set of m support or inducing variables that allow the reduction of the time complexity fromO(n3 ) to O(nm2 ). They mainly differ in the strategies they use to select the inducing inputs which are
typically selected from the training or test examples. Snelson and Ghahramani (2006) allow the inducing
variables to be considered as auxiliary pseudo-inputs that are inferred along with kernel hyperparameters
using continuous optimization.
The selection of inducing variables and kernel hyperparameters in a sparse GP method can be approached as a model selection problem. Therefore, the most relevant criterion for solving such problem
is an approximation to the exact (intractable) marginal likelihood that can be maximized over inducing
1

inputs and hyperparameters. Existing state-of-the-art methods (Snelson and Ghahramani, 2006; Seeger
et al., 2003) derive such approximations by modifying the likelihood function or the GP prior (QuiñoneroCandela and Rasmussen, 2005) and then computing the marginal likelihood of the modified model. Thisapproach turns the inducing inputs into additional kernel hyperparameters. While this can increase flexibility when we fit the data, it can also lead to overfitting when we optimize with respect to all unknown
hyperparameters. Furthermore, fitting a modified model to the data is not so rigorous approximation
procedure since there is no distance or divergence between the exact and the modified model that is
minimized.
In this paper we introduce a variational method for sparse GP models that jointly selects the inducing
inputs and the hyperparameters by maximizing a lower bound to the exact marginal likelihood. In thisformulation we do not modify the likelihood or the GP prior in the training examples. Instead we follow
the standard variational approach (Jordan et al., 1999) according to which a variational distribution isused to approximate the exact posterior over the latent function values. The important difference between
this formulation and previous methods is that here the inducing inputs are defined to be variationalparameters which are selected by minimizing the Kullback-Leibler (KL) divergence. This allows i) to
avoid overfitting and ii) to rigorously approximate the exact GP model by minimizing a divergence
between the sparse model and the exact one. The selection of the inducing inputs and hyperparametersis achieved either by assuming pseudo-inputs and applying continuous optimization over all unknown
quantities, similarly to (Snelson and Ghahramani, 2006), or by using a variational EM algorithm whereat the E step we greedily select the inducing inputs from the training data and at the M step we update
the hyperparameters. In contrast to previous greedy approaches, e.g. (Seeger et al., 2003), our scheme
monotonically increases the optimized objective function.
We apply the variational method to regression with additive Gaussian noise and we compare its
performance to training schemes based on the projected process marginal likelihood (Seeger et al., 2003;
Csato and Opper, 2002) and the sparse pseudo-inputs marginal likelihood (Snelson and Ghahramani,
2006).
Our method is most closely related to the variational sparse GP method described in (Csato, 2002;
Csato and Opper, 2002; Seeger, 2003) that is applied to GP classification (Seeger, 2003). The maindifference between our formulation and these techniques is that we maximize a variational lower boundin order to select the inducing inputs, while these methods use variational bounds for estimating only
the kernel hyperparameters. This paper is a revised and extended version of (Titsias, 2009).

2

Gaussian process regression and sparse methods

A GP is a collection of random variables {f (x)|x ∈ X }, where X is an index or input set, for which any
finite subset follows a Gaussian distribution. A GP is fully specified by the mean function m(x) and thecovariance or kernel function k(x, x′ ) defined bym(x) = E[f (x)],(1)

k(x, x′ ) = E[(f (x) − m(x))(f (x′ ) − m(x′ ))].

(2)

In the context of regression and statistical learning, a GP can be used as a non-parametric prior over areal-valued function which can be combined with data to give a posterior over the function. Suppose that2

we wish to estimate a real-valued function f (x), where for simplicity the input set X is taken to be theRD-dimensional real space D . We shall call f (x) the unobserved or latent function. We further assume
that the mean function of the GP prior is zero and the covariance function is specified through a set of
hyperparameters θ. Suppose we collect a training dataset, {(xi , yi )}ni=1 , consisting of n noisy realizationsof the latent function where each scalar yi is obtained by adding Gaussian noise to f (x) at input xi , i.e.
yi = fi + ǫi , ǫi ∼ N (0, σ 2 ), fi = f (xi ).
We denote by X the n × D matrix of all training inputs, y = [y1 . . . yn ]T the vector of all outputs andf = [f1 . . . fn ]T the corresponding vector of the training latent function values. The marginalization
property of GPs allows us to simplify the (initially infinite dimensional) prior so that after marginalizingout all function points not associated with the data, we obtain a n-dimensional Gaussian distribution,p(f ) = N (f |0, Knn ), where 0 denotes the n-dimensional zero vector and Knn is the n × n covariancematrix obtained by evaluating the kernel function on the observed inputs. The joint probability modelof observed output and latent variables data can be written asp(y, f ) = p(y|f )p(f ),where p(y|f ) is the likelihood function and p(f ) the GP prior. Notice that this is a conditional model
since we condition on the observed inputs X. Nevertheless, for the sake of clarity we omit reference to
X and the hyperparameters1 throughout the paper. The training dataset induce a posterior process overthe latent function f (x) which (because of the Gaussian likelihood) is also a GP specified by a posteriormean function and a posterior covariance function. We can easily work out the functional form of thesefunctions, e.g. by augmenting the marginal probability p(y) of the observed outputs by two latent function
values f (x) and f (x′ ), and find thatmy (x) = Kxn (σ 2 I + Knn )−1 y,(3)ky (x, x′ ) = k(x, x′ ) − Kxn (σ 2 I + Knn )−1 Knx′ .
Here, Kxn = [k(x, x1 ) . . . k(x, xn )] is an n-dimensional row vector of kernel function values between x andT. Any query related to the posterior GP can be answered by the above
the training inputs and Knx = Kxn
mean and covariance functions. For instance, the Gaussian posterior distribution p(f |y) on the training
latent variables f is computed by evaluating eq. (3) at the inputs X. Similarly the prediction of the outputy∗ = f∗ + ǫ∗ at some unseen input x∗ is described by p(y∗ |y) = N (y∗ |my (x∗ ), ky (x∗ , x∗ ) + σ 2 ). The
posterior GP depends on the values of the hyperparameters (θ, σ 2 ) which can be estimated by maximizing
the log marginal likelihood given bylog p(y) = log[N (y|0, σ 2 I + Knn )].

(4)

Once we have obtained point estimates for the hyperparameters (σ 2 , θ) by maximizing the above log
marginal likelihood, we can use these estimates in eq. (3) in order to make prediction in unseen input
points.
Although the above GP framework is elegant, it requires O(n3 ) computations as clearly we need to
1 A precise notation is to write p(y, f |X, σ 2 , θ) = p(y|f , σ 2 )p(f |X, θ).

3

invert a matrix of size n × n; once when we evaluate the prediction in eq. (3) and multiple times when
we maximize the marginal likelihood in eq. (4). Therefore, we need to consider approximate or sparse
methods in order to deal with large datasets. Advanced sparse methods use a small set of m function
points as support or inducing variables. This yields a time complexity that scales as O(nm2 ). Some
important issues in these methods involve the selection of the inducing variables and the hyperparameters.
For reviews of current approaches see chapter 8 in (Rasmussen and Williams, 2006) and (QuiñoneroCandela and Rasmussen, 2005). In section 3, we propose a variational framework to deal with theselection of the inducing variables and hyperparameters. Therefore, in the remaining of this section weanalyze the most relevant previous methods which are based on maximizing an approximate marginal
likelihood obtained by modifying the GP prior (Quiñonero-Candela and Rasmussen, 2005).
Suppose we wish to use m inducing variables to construct our sparse GP method. The inducing variables are latent function values evaluated at some inputs Xm . Xm can be a subset of the training inputs or
auxiliary pseudo-points (Snelson and Ghahramani, 2006). Learning Xm and the hyperparameters (θ, σ 2 )is the crucial problem we need to solve in order to obtain a sparse GP method. An approximation to thetrue log marginal likelihood in eq. (4) can allow us to infer these quantities. The current state-of-the-art
approximate marginal likelihood is given in the sparse pseudo-inputs GP method (SPGP) proposed in(Snelson and Ghahramani, 2006). A related objective function used in (Seeger et al., 2003) corresponds
to the projected process approximation (PP). These approximate log marginal likelihoods have the form
F = log[N (y|0, σ 2 I + Qnn )],(5)−1
where Qnn is an approximation to the true covariance Knn . In PP, Qnn = Knm Kmm
Kmn , i.e. the exact
covariance is replaced by the Nyström approximation. Here, Kmm is the m × m covariance matrix on

the inducing inputs, Knm is the n × m cross-covariance matrix between training and inducing points−1−1TKmn , i.e. the Nyström
Kmn ] + Knm Kmm. In SPGP, Qnn = diag[Knn − Knm Kmmand Kmn = Knmapproximation is corrected to be exact in the diagonal. By contrasting eq. (4) with (5), it is clear that F
is obtained by modifying the GP prior. This implies that the inducing inputs Xm play the role of extra
kernel hyperparameters (similar to θ) that parametrize the covariance matrix Qnn . However, because theprior has changed, continuous optimization of F with respect to Xm does not reliably approximate theexact GP model because there is no any distance between the modified and the exact GP model that is
minimized. Further, since F is heavily parametrized with the extra hyperparameters Xm , overfitting canoccur especially when we jointly optimize over (Xm , θ, σ 2 ). Despite all that, the flexibility introduced by
the extra hyperparameters Xm can often be advantageous. For instance, unlike the exact GP model, theSPGP model can fit heteroscedastic noise in the output data (Snelson and Ghahramani, 2006).
In the next section, we propose a formulation for sparse GP regression that follows a different philosophy to what eq. (5) implies. Rather than modifying the exact GP model and maximizing the marginal
likelihood of the modified model, we minimize the KL divergence between the exact posterior GP anda variational approximation. The inducing inputs Xm become now variational parameters which are
selected so as the KL divergence is minimized.

3

Variational learning of inducing inputs

We wish to define a sparse method that directly approximates the mean and covariance functions in
eq. (3) that characterize the posterior GP. This posterior GP can be also described by the predictive
4

Gaussian distributionp(z|y) =Zp(z|f )p(f |y)df ,where z is any finite set of function points, p(z|f ) denotes the conditional GP prior and p(f |y) is theposterior distribution over the training latent function values. Suppose that we wish to approximate theabove Bayesian integral by using a small set of m auxiliary inducing variables fm evaluated at the pseudoinputs Xm , which are independent from the training inputs. We further assume that fm are function
values drawn from the same GP prior2 as the training function values f . By using the augmented joint
model p(y|f )p(z, f , fm ), where p(z, fm , f ) is the GP prior jointly expressed over the function values z, f
and fm , we can equivalently write p(z|y) asp(z|y) =Zp(z|fm , f )p(f |fm , y)p(fm |y)df dfm .

(6)This expanded way of writing the predictive distribution is rather instructive since it indicates what one
should expect from a good set of inducing variables and in fact it reveals the properties of an optimal
set. More precisely, suppose that fm is a sufficient statistic for the parameter f in the sense that z and f
are independent given fm , i.e. it holds p(z|fm , f ) = p(z|fm ) for any z. The above can be written as
q(z)

=

Z

=

Z

p(z|fm )p(f |fm )φ(fm )df dfm
Zp(z|fm )φ(fm )dfm = q(z, fm )dfm ,(7)

where q(z) = p(z|y) and φ(fm ) = p(fm |y). Here, p(f |fm ) = p(f |fm , y) is true since y is a noisy version off and because of the assumption we made that any z is conditionally independent from f given fm 3 . The
predictive distribution in eq. (7) requires reference to only m function points, the inducing variables fm ,and it can be computed in O(nm2 ) time; see eq. (8) below. However, in practise the assumption of fm
being a sufficient statistic is unlikely to hold and we should expect q(z) to be only an approximation to
the exact predictive distribution p(z|y). In such case, and in order to think about how to optimize thequality of the approximation, we can let φ(fm ) be a “free” variational Gaussian distribution, where in
general φ(fm ) 6= p(fm |y), that depends on a mean vector µ and a covariance matrix A. Notice also thatthe quality of the approximation will crucially depend on the locations Xm of the inducing variables. By
using eq. (7), we can read off the mean and covariance functions of the approximate posterior GP andobtain:−1mqy (x) = Kxm Kmm
µ,

(8)

−1
−1
−1
kyq (x, x′ ) = k(x, x′ ) − Kxm Kmm
Kmx′ + Kxm Kmm
AKmm
Kmx′ .

The above defines the general form of the sparse posterior GP which is tractably computed in O(nm2 )time. The question that now arises is how can we select the φ distribution, i.e. (µ, A), and the inducing
inputs Xm . Notice that different sparse GP approaches such as the PP and SPGP methods, also called
DTC and FITC in the unified view of Quiñonero-Candela and Rasmussen (2005), use different forms forthe variational distribution φ and follow different strategies for the selection of the inducing inputs andhyperparameters.
2 More general inducing variables are defined in section 6.
R
3 From p(z|f , y) = R p(y|f )p(z,fm ,f )df and by using the fact p(z|f , f ) = p(z|f ), the result follows.
mmmp(y|f )p(z,fm ,f )df dz

5

Next we apply a variational inference method that allows to jointly specify the quantities (Xm , φ) andcrucially treat the inducing inputs as variational parameters which are rigorously selected by minimizing
the KL divergence. As will be shown our method regarding the form of the φ distribution leads to
the PP prediction (Csato and Opper, 2002; Seeger, 2003), while it gives a novel way of specifying theinducing inputs Xm and the hyperparameters (σ 2 , θ) by maximizing a lower bound to the exact marginal
likelihood.

3.1

The variational lower boundTo select the inducing inputs Xm , we intend to apply variational inference in an augmented probability
space that involves both the training latent function values f and the pseudo-input inducing variables
fm . Particularly, the initial joint model p(y, f ) is augmented with the variables fm to form the modelp(y, f , fm ) = p(y|f )p(f |fm )p(fm ),

(9)

−1
−1
Kmn ). As a confm , Knn −Knm Kmm
where the conditional GP prior is given by p(f |fm ) = N (f |Knm Kmm
R
sequence of the marginalization property of the Gaussian process according to which p(f ) = p(f , fm )dfm ,the two models are equivalent in terms of doing exact inference, i.e. computing the posterior p(f |y) andthe marginal likelihood p(y). For instance, the log marginal likelihood in eq. (4) can be equivalently
written aslog p(y) = logZp(y|f )p(f |fm )p(fm )df dfm .

(10)

However, the augmented probability model is more flexible in terms of doing approximate inference since
it contains a set of parameters, that is the inducing inputs Xm , which are somehow arbitrary. Moreprecisely, these parameters do not affect the exact GP model, p(y, f ), because p(f ) is not changing by
varying the values of Xm despite the fact that p(f |fm ) and p(fm ) do change. Therefore, Xm are not model
parameters (as (σ 2 , θ) are) and by applying approximate inference in the augmented probability model
we can turn them into variational parameters. This is what we do next.
We want to approximate the true posterior distribution p(f , fm |y) = p(f |fm , y)p(fm |y) by introducing
a variational distribution q(f , fm ) and minimizing the KL divergence:
KL(q(f , fm )||p(f , fm |y)) =Z

q(f , fm ) log

q(f , fm )df dfm .
p(f , fm |y)(11)This is based on the standard variational approach widely used in machine learning (Jordan et al., 1999).
To specify the form of the variational distribution q(f , fm ), we follow the arguments exposed earlier insection 3. Particularly, for an optimal setting of the inducing variables, the exact posterior p(f , fm |y)factorizes as p(f , fm |y) = p(f |fm )p(fm |y). This tells us that the variational distribution must satisfy thesame factorization as well, in order for the minimization of the KL divergence to search for these optimal
inducing variables. Thus,q(f , fm ) = p(f |fm )φ(fm ),(12)where φ(fm ) is an unconstrained variational distribution over fm and p(f |fm ) is the conditional GP prior.
Notice also that the form of this distribution directly follows from eq. (7). To determine the variational
quantities (Xm , φ), we minimize the KL divergence in eq. (11), which is equivalently expressed as the6maximization of the following variational lower bound on the true log marginal likelihood:log p(y) ≥ FV (Xm , φ)p(y|f )p(f |fm )p(fm )df dfm ,p(f |fm )φ(fm ) log
p(f |fm )φ(fm )ZZp(fm )=
φ(fm )
p(f |fm ) log p(y|f )df + log
dfm ,φ(fm )=

Z(13)

where the term p(f |fm ), in the first line inside the log, cancels out4 . We can firstly maximize the boundby analytically solving for the optimal choice of the variational distribution φ. To do this, we firstly
compute the integral
log G(fm , y)

=

Z

=



1
log N (y|α, σ 2 I) − 2 T r(Knn − Qnn ).
2σp(f |fm ) log p(y|f )df
(14)

−1
−1
Kmn . Notice that α and Knn − Qnn are thefm and Qnn = Knm Kmm
Here, α = E[f |fm ] = Knm Kmmmean vector and covariance matrix, respectively, of the conditional GP prior p(f |fm ). Eq. (13) is writtenasZ
G(fm , y)p(fm )FV (Xm , φ) = φ(fm ) log
(15)
dfm .
φ(fm )
We can now maximize the bound with respect to the distribution φ, without computing the optimal
distribution, called φ∗ , itself. This is done by reversing the Jensen’s inequality, i.e. moving the log
outside of the integral, which gives1
FV (Xm ) = log N (y|0, σ 2 I + Qnn ) − 2 T r(Knn − Qnn ),
2σ

(16)

−1
where Qnn = Knm Kmm
Kmn . More details of the derivation of this bound are given in the appendix A.

The novelty of the above objective function is that it contains a regularization trace term: − 2σ1 2 T r(Knn −Qnn ). This clearly differentiates FV from all marginal likelihoods, described by eq. (5), that were previously applied to sparse GP regression. We will analyze the trace term shortly.
The quantity in eq. (16) is computed in O(nm2 ) time and is a lower bound of the true log marginal
likelihood for any value of the inducing inputs Xm . Further maximization of the bound can be achievedby optimizing over Xm and optionally over the number of these variables. Notice that the inducing
inputs determine the flexibility of the variational distribution q(f , fm ) = p(f |fm )φ∗ (fm ) since by tuning
Xm we adapt both p(f |fm ) and the underlying optimal distribution φ∗ . To compute this optimal φ∗ , wedifferentiate eq. (15) with respect to φ(fm ) without imposing any constraints on the functional form ofφ(fm ) apart from being a distribution. This gives (see appendix A):φ∗ (fm ) = N (fm |µ, A),(17)

where µ = σ −2 Kmm ΣKmn y, A = Kmm ΣKmm and Σ = (Kmm + σ −2 Kmn Knm )−1 . This now fully
specifies our variational GP and we can substitute (µ, A) in eq. (8) in order to make predictions in
unseen input points. Clearly, the predictive distribution is exactly the one used by the projected process
4 ThisRcancellationis actually not needed as we can arrive faster at eq. (13) by writing log p(y)
R
p(f

)

p(y|f )p(f |f

)df

≥

m
m
dfm and then apply once again Jensen’s inequality to lower bound the termφ(fm ) log
φ(fm )
Rlog p(y|f )p(f |fm )df with respect to p(f |fm ).

7

approximation (PP) that has been previously proposed in (Csato and Opper, 2002; Seeger et al., 2003).
Thus, as far as the predictive distribution is concerned the above method is equivalent to PP.
However, the variational method is very different to PP and SPGP as far as the selection of theinducing inputs and the kernel hyperparameters is concerned. This is because of the extra regularization
term that appears in the bound in eq. (16) and does not appear in the approximate log marginal likelihoods
used in PP (Seeger et al., 2003) and SPGP (Snelson and Ghahramani, 2006). As discussed in section 2,
for the latter objective functions, the role of Xm is to form a set of extra kernel hyperparameters. In
contrast, for the lower bound, the inputs Xm become variational parameters due to the KL divergencethat is minimized.
To look into the functional form of the bound, note that FV is the sum of the PP log likelihood and theregularization trace term − 12 σ −2 T r(Knn − Qnn ). Thus, FV attempts to maximize the PP log likelihoodand simultaneously minimize the trace T r(Knn − Qnn ). This trace represents the total variance of theconditional prior p(f |fm ) which also corresponds to the squared error of predicting the training latentR−1
values f from the inducing variables fm : p(f , fm )||Knm Kmm
fm − f ||2 df dfm . When the trace term is−1zero, the Nyström approximation is exact, i.e. Knn = Knm Kmm
Kmn , which means that the variational
distribution q(f , fm ) matches exactly the true posterior distribution. Note that the trace T r(Knn − Qnn )itself has been used as a criterion for selecting the inducing points from the training data in (Smola andSchölkopf, 2000) and is similar to the criterion used in (Lawrence et al., 2002).
When we maximize the variational lower bound, the hyperparameters (σ 2 , θ) are regularized. It iseasy to see how this is achieved for the noise variance σ 2 . At a local maxima, σ 2 satisfies:
1
σ =
n
2

Z

φ∗ (fm )||y − α||2 dfm +1
T r(Knn − Qnn ),
n

(18)

−1
where ||z|| denotes the Euclidean norm and α = E[f |fm ] = Knm Kmm
fm . This decomposition reveals

that the obtained σ 2 will be equal to the estimated “actual” noise plus a “correction” term that is theaverage squared error associated with the prediction of the training latent values f from the inducing
variables fm . Thus, the variational lower bound naturally prefers to set σ 2 larger than the “actual” noise
in a way that is proportional to the inaccuracy of the approximation.
So far we assumed that the inducing variables correspond to pseudo-inputs which are are selected byapplying gradient-based optimization. However, this can be difficult in high dimensional input spaces
as the number of variables that need to be optimized becomes very large. Further, the kernel function
might not be differentiable with respect to the inputs. see e.g. kernels defined in strings. In such cases
we can still apply the variational method by selecting the inducing inputs from the training inputs. An
important property of this discrete optimization scheme is that the variational lower bound monotonicallyincreases when we greedily select inducing inputs and adapt the hyperparameters. Next we discuss thisgreedy selection method.

3.2

Greedy selection from the training inputs

Let m ⊂ {1, . . . , n} be the indices of a subset of data that are used as the inducing variables. The training
points that are not part of the inducing set are indexed by n − m and are called the remaining points, e.g.
fn−m denotes the remaining latent function values. The variational method is applied similarly to thepseudo-inputs case. Assuming the variational distribution q(f ) = p(fn−m |fm )φ(fm ), we can express a variational bound that has the same form as the bound in eq. (16) with the only difference that the covariance8−1
Km(n−m) .
matrix in the trace term is now given by Cov(fn−m |fm ) = K(n−m)(n−m) − K(n−m)m Kmm

The selection of inducing variables among the training data requires a prohibitive combinatorial
search. A suboptimal solution can be based on a greedy selection scheme where we start with an empty
inducing set m = ∅ and a remaining set n − m = {1, . . . , n}. At each iteration, we add a training pointj ∈ J ⊂ n − m, where J is a randomly chosen working set, into the inducing set that maximizes theselection criterion ∆j .
It is important to interleave the greedy selection process with the adaption of the hyperparameters(σ , θ). This can be viewed as an EM-like algorithm; at the E step we add one point into the inducing
2

set and at the M step we update the hyperparameters. To obtain a reliable convergence, the approximate
marginal likelihood must monotonically increase at each E or M step. The PP and SPGP log likelihoods
do not satisfy such a requirement because they can also decrease as we add points into the inducing set.
In contrast, the bound FV is guaranteed to monotonically increase since now the EM-like algorithm is avariational EM. To clarify this, we state the following proposition.
Proposition 1. Let (Xm , fm ) be the current set of inducing points and m the corresponding set ofindices. Any point i ∈ n − m added into the inducing set can never decrease the lower bound.
Proof: Before the new point (fi , xi ) is added, the variational distribution is p(fn−m |fm )φ∗ (fm ) =p(fn−(m∪i) |fi , fm )p(fi |fm )φ∗ (fm ). When we add the new point, the term p(fi |fm )φ∗ (fm ) is replaced by
the optimal φ∗ (fi , fm ) distribution. This can either increase the lower bound or leave it invariant. A
more detailed proof is given in the appendix B.
A consequence of the above proposition is that the greedy selection process monotonically increases
the lower bound and this holds for any possible criterion ∆. An obvious choice is to use FV as thecriterion, which can be evaluated in O(nm) time for any candidate point in the working set J. Such a
selection process maximizes the decrease in the divergence KL(q(f )||p(f |y)).
Finally, we should point out that the selection of the inducing variables from the training data should
not necessarily be performed based on the greedy selection strategy discussed above. One could consider
different strategies such as adding groups of training points simultaneously into the active set or swapping
points between the active and the remaining sets. Any strategy for reselecting the active set simply
updates the variational distribution q(f ), and as long as this is done so that the bound increases, thewhole procedure leads to a valid variational EM algorithm.

4Comparison of the objective functions

In this section we compare the lower bound FV , the PP and the SPGP log likelihood in some toy problems.
All these functions are continuous with respect to (Xm , σ 2 , θ) and can be maximized using gradient-based
optimization.
Our working example will be the one-dimensional dataset5 considered in Snelson and Ghahramani(2006) that consists of 200 training points; see Figure 1. We train a sparse GP model using the squared
exponential kernel defined by σf2 exp(− 2ℓ12 ||xi − xj ||2 ). Since the dataset is small and the full GP model is
tractable, we compare the sparse approximations with the exact GP prediction. The plots in the first row
of Figure 1 show the predictive distributions for the three methods assuming 15 inducing inputs. The left
plot displays the mean prediction with two-standard error bars (shown as blue solid lines) obtained by the
maximization of FV . The prediction of the full GP model is displayed using dashed red lines. The middle
5 obtained from www.gatsby.ucl.ac.uk/∼snelson/.

9

Figure 1: The first row corresponds to 200 training points and the second row to 20 training points. The firstcolumn shows the prediction (blue solid lines) obtained by maximizing FV over the 15 pseudo-inputs and thehyperparameters. The full GP prediction is shown with red dashed lines. Initial locations of the pseudo-inputs
are shown on the top as crosses, while final positions are given on the bottom as crosses. The second column
shows the predictive distributions found by PP and similarly the third column for SPGP.

plot shows the corresponding solution found by PP and the right plot the solution found by SPGP.
The prediction obtained by the variational method almost exactly reproduces the full GP prediction.
The final value of the variational lower bound was −55.5708, while the value of the maximized truelog marginal likelihood was −55.5647. Further, the estimated hyperparameters found by FV match thehyperparameters found by maximizing the true log marginal likelihood. In contrast, training the sparsemodel using the PP log likelihood gives a poor approximation. The SPGP method gave a much moresatisfactory answer than PP although not as good as the variational method.
To consider a more challenging problem, we decrease the number of the original 200 training examples
by maintaining only 20 of them6 . We repeat the experiment above using exactly the same setup. The
second row of Figure 1, displays the predictive distributions of the three methods. The prediction of thevariational method is identical to the full GP prediction and the hyperparameters match those obtainedby full GP training. On the other hand, the PP log likelihood leads to a significant overfitting of the
training data since the mean curve interpolates the training points and the error bars are very noisy.
SPGP provides a solution that significantly disagrees with the full GP prediction both in terms of themean prediction and the errors bars. Notice that the width of the error bars found by SPGP varies a lot
in different input regions. This nonstationarity is achieved by setting σ 2 very close to zero and modelling−1
the actual noise by the heteroscedastic diagonal matrix diag[Knn − Knm Kmm
Kmn ]. The fact that thisdiagonal matrix (the sum of its elements is the trace T r(Knn − Qnn )) is large indicates that the full GP

model is not well approximated.
The reason PP and SPGP do not recover the full GP model when we optimize over (Xm , σ 2 , θ) is not
the local maxima. To clarify this point, we repeated the experiments by initializing the PP and SPGP
log likelihoods to optimal inducing inputs and hyperparameters values where the later are obtained byfull GP training. The predictions found are similar to those shown in Figure 1. A way to ensure that
6 The points were chosen from the original set according to the MATLAB command: X = X(1:10:end).

10the full GP model will be recovered as we increase the number of inducing inputs is to select them fromthe training inputs. This, however, turns the continuous optimization problem into a discrete one andmoreover PP and SPGP face the non-smooth convergence problem.
Regarding FV , it is clear from section 3 that by maximizing over Xm we approach the full GP
model in the sense of KL(q(f , fm )|p(f , fm |y)). Something less clear is that FV efficiently regularizes thehyperparameters (σ 2 , θ) so as overfitting is avoided. This is achieved by the regularization trace term:− 12 σ −2 T r(Knn − Qnn ). When this trace term is large because there are not sufficiently many inducing
variables, this term favours kernel parameters that give a smoother function. Also, when the trace termis large the decomposition in eq. (18) implies that σ 2 must increase as well. These properties are usefulfor avoiding overfitting and also imply that the prediction obtained by FV will tend to be smoother thanthe prediction of the full GP model. In contrast, the PP and SPGP log likelihoods can find more flexible
solutions than the full GP prediction which indicates that they are prone to overfitting.

5

Experiments

In this section we compare the variational lower bound (VAR), the projected process approximate log
likelihood (PP) and the sparse pseudo-inputs GP (SPGP) log likelihood in four real datasets. As abaseline technique, we use the subset of data (SD) method. For all sparse GP methods we jointly
maximize the alternative objective functions w.r.t. hyperparameters (θ, σ 2 ) and the inducing inputs Xm
using the conjugate gradients algorithm. Xm is initialized to a randomly chosen subset of training
inputs. In each run all methods are initialized to the same inducing inputs and hyperparameters. The2∗ −f∗ ||
performance criteria we use are the standardized mean squared error (SMSE), given by T1 ||yvar(y∗ ) , andthe standardized negative log probability density (SNLP) as defined in (Rasmussen and Williams, 2006).
Smaller values for both error measures imply better performance. In all the experiments we use thesquared-exponential kernel with varied length-scale.
Firstly, we consider the Boston-housing dataset, which consists of 455 training examples and 51
test examples. Since the dataset is small, full GP training is tractable. In the first experiment, we
fix the parameters (θ, σ 2 ) to values obtained by training the full GP model. Thus we can investigate
the difference of the methods solely on how the inducing inputs are selected. We rigorously compare
the methods by calculating the moments-matching divergence KL(p(f∗ |y)||q(f∗ )) between the true test
posterior p(f∗ |y) and each of the approximate test posteriors. For the SPGP method the approximate test
posterior distribution is computed by using the exact test conditional p(f∗ |fm ). Figure 2(a) show the KL
divergence as the number of inducing points increases. Means and one-standard error bars were obtainedby repeating the experiment 10 times. Note that only the VAR method is able to match the full GP
model; for around 200 points we closely match the full GP prediction. Interestingly, when the inducing
inputs are initialized to all training inputs, i.e. Xm = X, PP and SPGP still give a different solution
from the full GP model despite the fact that the hyperparameters are kept fixed to the values of thefull GP model. The reason this is happening is that they are not lower bounds to the true log marginal
likelihood and as shown in Figure 2(c) they become upper bounds. To show that the effective selection of
the inducing inputs achieved by VAR is not a coincidence, we compare it with the case where the inputs
are kept fixed to their initial randomly selected training inputs. Figure 2(b) displays the evolution of the
KL divergence for the VAR, the random selection plus PP (RSPP) and the SD method. Note that the
only difference between VAR and RSPP is that VAR optimizes the lower bound over the initial values of

11

60

60
50

KL(p||q)

50
40
30

40
30

20

20

10

10

0

0

−500

−1000
FullGP
VAR
PP
SPGP

0
100
200
300
400
Number of inducing variables

−1500

100
200
300
400
Number of inducing variables

(a)

(b)

0.2
FullGP
VAR
PP
SPGP

1

SNLP

SMSE

0.15

0.1

FullGP
VAR
PP
SPGP

1.5

100
200
300
400
Number of inducing variables

(c)

2

0.5
0
−0.5

400

Log marginal likelihood

KL(p||q)

VAR
RSPP
SD

70

Log marginal likelihood

VAR
PP
SPGP

70

200
0
−200
FullGP
VAR
PP
SPGP

−400
−1
0.05

100
200
300
400
Number of inducing variables

(d)

−1.5

100
200
300
400
Number of inducing variables

(e)

−600

100
200
300
400
Number of inducing variables(f)Figure 2: (a) show the KL divergence as the number of inducing variables increases for the VAR the PP and SPGP
methods. Similarly (b) show the divergence for the VAR, RSPP and SD methods. (c) displays the approximate
log marginal likelihoods; the true log marginal likelihood value is displayed by using the dotted horizontal line.
(d) and (e) show the SMSE and SNLP errors (obtained by joint learning hyperparameters and inducing-inputs)
against the number of inducing variables. (f) shows the corresponding log marginal likelihoods.

the inducing inputs, while RSPP just keep them fixed. Clearly RSPP significantly improves over the SD
prediction, and VAR significantly improves over RSPP.
In a second experiment, we jointly learn inducing variables and hyperparameters and compare themethods in terms of the SMSE and SNLP errors. The results are displayed in the second row of Figure
2. Note that the PP and SPGP methods achieve a much higher log likelihood value (Figure 2(f)) than
the true log marginal likelihood. However, the error measures clearly indicate that the PP log likelihood
significantly overfits the data. SPGP gives better SMSE error than the full GP model but it overfits
w.r.t. the SNLP error. The variational method matches the full GP model.
We now consider three large datasets: the kin40k dataset, the sarcos and the abalone datasets7that have been widely used before. Note that the abalone dataset is small enough so as we will be able
to train the full GP model. The inputs were normalized to have zero mean and unit variance on the
training set and the outputs were centered so as to have zero mean on the training set. For the kin40k
and the sarcos datasets, the SD method was obtained in a subset of 2000 training points. We vary thesize of the inducing variables in powers of two from 16 to 1024. For the sarcos dataset, the experimentfor 1024 was not performed since is was unrealistically expensive. All the objective functions were jointly
maximized over inducing inputs and hyperparameters. The experiment was repeated 5 times. Figure 3
shows the results.
From the plots in Figure 3, we can conclude the following. The PP log likelihood is significantly
prone to overfitting as the SNLP errors clearly indicate. However, note that in the kin40k and sarcos
7 kin40k: 10000 training, 30000 test, 8 attributes, ida.first.fraunhofer.de/ anton/data.html.
sarcos: 44, 484 training, 4, 449 test, 21 attributes, www.gaussianprocess.org/gpml/data/.
abalone: 3, 133 training, 1, 044 test, 8 attributes, www.liaad.up.pt/ ltorgo/Regression/DataSets.html.

12

0.05
SD
VAR
PP
SPGP

0.5
SD
VAR
PP
SPGP

0.045
0.04

0.1

SMSE

SMSE

SMSE

0.035
0.03
0.025

0.05

FullGP
VAR
PP
SPGP

0.48

0.02

0.46

0.44

0.42

0.015
0

200
400
600
800
Number of inducing variables

0.01

1000

1

−0.5
−1

0.4

500

−1
−1.5

−1.5

0

200
400
600
800
Number of inducing variables

1000

0.1

SD
VAR
PP
SPGP

−0.5

SNLP

SNLP

0

100
200
300
400
Number of inducing variables

0

SD
VAR
PP
SPGP

0.5

0

FullGP
VAR
PP
SPGP

0
−0.1

SNLP

0

−0.2
−0.3
−0.4

−2

−2
−2.5
0

200
400
600
800
Number of inducing variables

1000

−2.5

−0.5

0

100
200
300
400
Number of inducing variables

500

0

200
400
600
800
Number of inducing variables

1000

Figure 3: The first column displays the SMSE (top) and SNLP (bottom) errors for the kin40k dataset withrespect to the number of inducing points. The second column shows the corresponding plots for the sarcos
dataset and similarly the third column shows the results for the abalone dataset.

datasets, PP gave the best performance w.r.t. to SMSE error. This is probably because of the ability of
PP to interpolate the training examples that can lead to good SMSE error when the actual observationnoise is low. SPGP often has the worst performance in terms of the SMSE error and almost always thebest performance in terms of the SNLP error. In the abalone dataset, SPGP had significantly better
SNLP error than the full GP model. Since the SNLP error depends on the predictive variances, we
believe that the good performance of SPGP is due to its heteroscedastic ability. For example, in the
kin40k dataset, SPGP makes σ 2 almost zero and thus the actual noise in the likelihood is modelled−1by the heteroscedastic covariance diag[Knn − Knm Kmm
Kmn ]. The fact that the latter term is large
may indicate that the full GP model is not well approximated. Finally the variational method has good
performance. VAR never had the worst performance and it didn’t exhibit overfitting. The examples in
section 4, the Boston-housing and the abalone dataset indicate that the VAR method remains muchcloser to the full GP model than the other methods.

6

Adding “jitter” and more general inducing variables

In this section, and in order to address future research regarding sparse GP methods using inducing
variables, we discuss the use of more general inducing variables than those being just function points.
We do this by giving an example of how to slightly generalize the variational inference method so thatto numerically stabilize the lower bound and the sparse GP prediction. This involves the usual trick ofadding a “jitter” factor to the covariance matrix Kmm of the inducing points. We will show that this isactually rigorous and has a precise mathematical interpretation. Our discussion will have a more general
outlook.
To start with, observe that when we maximize the exact log marginal likelihood in eq. (4) (assuming
for now that we deal with small datasets) we need to invert the matrix Knn + σ 2 I. This is usually13carried out by computing the Cholesky decomposition of this matrix. While the kernel matrix Knn canbe (in computer precision) singular, i.e. it will have some zero eigenvalues, the matrix Knn + σ 2 willoften be strictly positive definite due to the addition of σ 2 in the diagonal of Knn . Thus, σ 2 naturally
acts as a “jitter” term that numerically stabilizes the computation of the Cholesky decomposition andsubsequently the whole optimization of the log marginal likelihood. This, however, is not true for thesparse approximations to the exact marginal likelihood, neither for the PP and SPGP log marginal
likelihoods in eq. (5) nor for the variational lower bound in eq. (16). The reason is that now we haveto explicitly invert the kernel matrix Kmm evaluated at the inducing inputs without adding a positive
number to the diagonal. To deal with that in a software implementation the common approach is to add
a small “jitter” in the diagonal of Kmm . But what is the consequence of this regarding the mathematicalproperties of the objective function that we optimize? Next we show that for the variational method
adding “jitter” leads to a generalized lower bound on the exact log marginal likelihood. In fact “jitter”
will turn out to be another variational parameter and we can optimize over it.
The variational method in section 3, applies inference in an expanded probability space. Particularly,
we started with the joint model p(y, f ), we augmented it to form p(y, f , fm ) and then fitted an augmented
variational distribution q(f , fm ) to the exact posterior p(f , fm |y). The key property that made this acorrect8 inference procedure is the consistency condition between the original and the augmented model:p(y, f ) =Z

p(y, f , fm )dfm ,Rwhich is true because fm are function values drawn from the GP prior, i.e. it holds p(f ) = p(f , fm )dfm .
Its rather interesting to observe now that there is nothing in the above augmentation procedure that saysthat the newly introduced random variables must be function points of exactly the same type as f . So
let us augment with a different type of inducing variables, denoted by λ, that are not precisely function
values drawn from the GP prior, but instead they are noisy versions of functions values:
λ = fm + ǫ, ǫ ∼ N (0, vIm ),(19)

where fm , as before, is a vector of m auxiliary GP function values evaluated at the pseudo-inputs Xmand v is a noise variance parameter. Note that when v = σ 2 , then λ can be considered as pseudo-output
data, however, this is rather restrictive and it is better to let v be a free parameter. The covarianceof λ is Kλλ = Kmm + vI and the cross covariance matrix between f and λ is Knλ = Knm . Thus
augmenting the Gaussian prior p(f ) with inducing variables λ we obtain p(f , λ) = p(f |λ)p(λ) which, ofR
course, satisfies the consistency condition p(f ) = p(f , λ)dλ. Therefore, exact inference in the modelp(y, f ) (e.g. computing the marginal likelihood p(y)) can be equivalently performed using the augmentedmodel
p(y, f , λ) = p(y|f )p(f |λ)p(λ).

(20)

However, doing approximate inference in the augmented model can be more flexible because the “free”
parameters (Xm , v) can be treated as variational parameters. Thus, following the method described in
section 3, and assuming the variational distribution q(f , λ) = p(f |λ)φ(λ), we compute a lower bound to
8 In the sense that when the KL divergence becomes zero, then the variational distribution allows to compute exactly
the true posterior p(f |y).

14

the exact log marginal likelihood:
FV (Xm , v)

=
−



log N (y|0, σ 2 I + Knm (Kmm + vIm )−1 Kmn )1
T r(Knn − Knm (Kmm + vIm )−1 Kmn ).
2σ 2

(21)

This only differs from the initial bound of eq. (16) in that the “jitter” term, that is v, has been added
to the diagonal of Kmm which makes the above bound slightly more general. To summarize, we justshowed that adding “jitter” in the diagonal of Kmm gives always a lower bound on log p(y). The sparseR
GP prediction in eq. (8) is also affected since now it is expressed through q(z|y) = p(z|λ)φ(λ)dλ(analogously to eq. (7)) where the role of approximate sufficient statistics is played by λ. Notice thatthe “jitter” v is actually a variational parameters as Xm is, so we can optimize over both of them. A
sensible optimization strategy is to initialize v to a large value. This will help to numerically stabilize

the optimization at the first crucial iterations and allow also to escape from local maxima. After few
iterations v will typically become close to zero and we can constrain it to be larger than a certain minimum“jitter” value. The fact that v will approach zero as we optimize is because the noisy-free latent function
values fm are more informative about f compared to their noisy counterparts fm + ǫ. More precisely,
the accuracy of the variational distribution q(f , λ) = p(f |λ)φ(λ) depends on how much informative (i.e.
deterministic) is the conditional GP prior p(f |λ) and among all λs is eq. (19) the most deterministic
conditional prior is obtained when v = 0.
We address now the issue of what random variables could be used as inducing variables. A generalanswer is that any linear functional involving the GP function f (x) can be an inducing variable. An useful
set of inducing variables, λ, should allow to decrease significantly the uncertainty in the conditional priorp(f |λ), which is part of the variational distribution. So good λs are those that can predict f with highcertainty. From practical point of view inducing variables must be such that the covariance matrix Kλλ
and the cross covariance matrix Knλ are computed both in closed-form and in a reasonable time so thatthe complexity of maximizing the variational lower bound will remain of order O(nm2 ).
We shall now draw a more general picture regarding the variational sparse GP method and themethods based on modifying the GP model. Assume we define a set of inducing variables λ for which
the Knλ and Kλλ are computed analytically. Then, we can augment the joint model similarly to eq. (20).
This involves the introduction of some parameters θλ , for example, when we use the “jitter” inducing
variables λ = fm + ǫ, these parameters are θλ = (Xm , v). The variational method has the property that
automatically turns the augmentation parameters θλ into variational parameters that are selected so that
the KL divergence is minimized. In contrast, for the SPGP or PP method, or any other method basedon the philosophy of modifying the GP prior or likelihood, the θλ parameters will be part of the model
hyperparameters that are tuned in order to fit the data and not to approximate the exact GP model.

7

Discussion

We proposed a variational framework for sparse GP regression that can reliably learn inducing inputs andhyperparameters by minimizing the KL divergence between the true posterior GP and an approximate
one. This approach starts with the full GP model and applies standard variational inference without
introducing any approximation to the likelihood or the GP prior. The variational distribution follows thefactorization that holds for optimally selected inducing variables. This method can be more generally

15

applicable. An interesting topic for the future is to apply this method to GP models that assume multiple
latent functions and consider also models with non-Gaussian likelihood functions, such as classification.
Furthermore, for the standard GP regression model discussed in this paper, the variational inference
method gave rise to the PP/DTC prediction (Csato and Opper, 2002; Seeger et al., 2003). However, theSPGP/FITC prediction is considered to be more accurate that the PP/DTC prediction (Snelson, 2007;
Quiñonero-Candela and Rasmussen, 2005). This arises the question: can we reformulate FITC so thatthe selection of the inducing inputs and hyperparameters can be achieved by maximizing a variational
lower bound. An attempt to define such a bound is given in the appendix C. This bound is derived bystarting with a full GP model that does not correspond to the standard GP regression, but instead it isa more flexible model that explains heteroscedastic noise in the output data.
Acknowledgements
I am grateful to Neil Lawrence for his help. Also, I wish to thank Magnus Rattray, Mauricio Alvarez,
David Luengo Garcia and Chris Williams for useful comments. This work is funded by EPSRC Grant No
EP/F005687/1 ”Gaussian Processes for Systems Identification with Applications in Systems Biology”.

References
Csato, L. (2002). Gaussian Processes – Iterative Sparse Approximations. PhD thesis, Aston University,UK.
Csato, L. and Opper, M. (2002). Sparse online Gaussian processes. Neural Computation, 14:641–668.
Jordan, M., Ghahramani, Z., Jaakkola, T., and Saul, L. (1999). An introduction to variational methods forgraphical models. Machine Learning, 37(2):183–233.
Keerthi, S. and Chu, W. (2006). A Matching Pursuit approach to sparse Gaussian process regression. In Weiss,Y., Schölkopf, B., and Platt, J., editors, Neural Information Processing Systems 18. MIT Press.
Lawrence, N. D., Seeger, M., and Herbrich, R. (2002). Fast sparse Gaussian process methods: the informative
vector machine. In Neural Information Processing Systems, 13. MIT Press.
O’Hagan, A. (1978). Curve fitting and optimal design for prediction. Journal of the Royal Statistical Society,
Series B, 40(1):1–42.
Quiñonero-Candela, J. and Rasmussen, C. E. (2005). A unifying view of sparse approximate Gaussian process
regression. Journal of Machine Learning Research, 6:1939–1959.
Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.
Schwaighofer, A. and Tresp, V. (2003). Transductive and inductive methods for approximate Gaussian process
regression. In Neural Information Processing Systems 15. MIT Press.
Seeger, M. (2003). Bayesian Gaussian Process Models: PAC-Bayesian Generalisation Error Bounds and Sparse
Approximations. PhD thesis, University of Edinburgh.
Seeger, M., Williams, C. K. I., and Lawrence, N. D. (2003). Fast forward selection to speed up sparse Gaussian
process regression. In Ninth International Workshop on Artificial Intelligence. MIT Press.
Smola, A. J. and Bartlett, P. (2001). Sparse greedy Gaussian process regression. In Neural Information Processing
Systems, 13. MIT Press.

16

Smola, A. J. and Schölkopf, B. (2000). Sparse greedy matrix approximations for machine learning. In International
Conference on Machine Learning.
Snelson, E. (2007). Flexible and efficient Gaussian process models for machine learning. PhD thesis, Gatsby
Computational Neuroscience Unit, University College London.
Snelson, E. and Ghahramani, Z. (2006). Sparse Gaussian process using pseudo-inputs. In Neural Information
Processing Systems, 13. MIT Press.
Titsias, M. K. (2009). Variational learning of inducing variables in sparse gaussian processes. In Twelfth International Conference on Artificial Intelligence and Statistics, JMLR: W and CP, volume 5, pages 567–574.
Wahba, G. (1990). Spline Models for Observational Data. Society for Industrial and Applied Mathematics, 59.
Walder, C., Kim, I. K., and Schölkopf, B. (2008). Sparse multiscale gaussian process regression. In ICML, pages
1112–1119.
Williams, C. K. I. and Seeger, M. (2001). Using the Nyström method to speed up kernel machines. In Neural
Information Processing Systems 13. MIT Press.

A

Variational lower bound and optimal distribution φ∗

The true log marginal likelihood is written as follows:Z
log p(y) = log p(y|f )p(f |fm )p(fm )df dfm .

(22)

Introducing the variational distribution q(f , fm ) = p(f |fm )φ(fm ) and applying Jensen’s inequality we obtain thelower bound:
FV (Xm , φ)==Zp(y|f )p(f |fm )p(fm )df dfmp(f |fm )φ(fm )Zp(y|f )p(fm )df dfm .
p(f |fm )φ(fm ) log
φ(fm )
p(f |fm )φ(fm ) log(23)

It can be written as
FV (Xm , φ) =Z

φ(fm )

Zp(f |fm ) log p(y|f )df + log

p(fm )
φ(fm )

ﬀdfm .

The integral involving f is computed as follows
Z
log G(fm , y) =p(f |fm ) log p(y|f )df

Z
h
iﬀ
n
1
T
T
T
2
=
p(f |fm ) − log(2πσ ) − 2 T r yy − 2yf + ff
df
2
2σ
i
h
1
n
= − log(2πσ 2 ) − 2 T r yyT − 2yαT + ααT + Knn − Qnn ,
2
2σ

(24)

(25)

−1
−1
where α = E[f |fm ] = Knm Kmm
fm and Qnn = Knm Kmm
Kmn . From the above we can recognize the expressionˆ˜
1log G(fm , y) = log N (y|α, σ 2 I) − 2 T r(Knn − Qnn ).
2σ

(26)

FV (Xm , φ) is now written as
FV (Xm , φ)=Z

φ(fm ) log

N (y|α, σ 2 I)p(fm )1dfm − 2 T r(Knn − Qnn ).
φ(fm )
2σ

17

(27)

We can now maximize this bound w.r.t. the distribution φ. The usual way of doing this is to take the derivative2(y|α,σ I)p(fm )w.r.t. φ(fm ), set to zero and obtain the optimal φ∗ (fm ) = R NN(y|α,σ. Then by substituting φ∗ (fm ) back2 I)p(f )dfmm

to eq. (27) we can compute the optimal bound w.r.t. to the distribution φ. However, since φ was not constrained
to belong to any restricted family of distributions, a faster and by far simpler way to compute the optimal boundis by reversing the Jensen’s inequality, i.e. moving the log outside of the integral in eq. (27). This gives
FV (Xm )Z

N (y|α, σ 2 I)p(fm )dfm −1
T r(Knn − Qnn )
2σ 2

=

log

=

ˆ
˜
1
log N (y|0, σ 2 I + Qnn ) − 2 T r(Knn − Qnn ).
2σ

(28)

The optimal distribution φ that gives rise to this bound is given byφ∗ (fm )∝=N (y|α, σ 2 I)p(fm )

ﬀ
1 T
1 T
1 −1
−1
−1
−1
c exp − fm (Kmm + 2 Kmm Kmn Knm Kmm )fm + 2 y Knm Kmm fm ,
2
σ
σ

(29)

where c is a constant. Completing the quadratic form we recognize the Gaussianφ∗ (fm ) = N (fm |σ −2 Kmm Σ−1 Kmn y, Kmm Σ−1 Kmm ),(30)

where Σ = Kmm + σ −2 Kmn Knm .

B

Detailed proof of Proposition 1

Let the inducing variables be a subset of the training examples indexed by m ⊂ {1, . . . , n}. We also use m todenote the number of these variables. The training points that are not part of the inducing set are indexed byn − m and are called the remaining points, e.g. fn−m denotes the remaining latent function values.
The variational distribution we use to bound the true log marginal likelihood is
qm (f ) = p(fn−m |fm )φm (fm ),(31)

where we introduced m as an index of q and φ to emphasize that this variational distribution assumes m inducing
variables. The variational bound takes the formZp(y|f )p(fn−m |fm )p(fm )log p(y) ≥p(fn−m |fm )φm (fm ) log
dfp(fn−m |fm )φm (fm )Zfp(y|f )p(fm )df=p(fn−m |fm )φ(fm ) log
φ(fm )f= FV (Xm , φ).

(32)

The maximum value of this bound that corresponds to the optimal φ∗ (fm ) is
ˆ
˜
1
−1
−1
FV (Xm ) = log N (y|0, σ 2 I + Knm Kmm
Kmn ) − 2 T r(K(n−m)(n−m) − K(n−m)m Kmm
Km(n−m) ).
2σ

(33)

Note that FV (Xm ) ≥ FV (Xm , φm ) for any distribution φm (fm ).
Proposition 1. Let (Xm , fm ) be the current set of inducing variables and m the corresponding set of indices.
Any point i ∈ n − m added into the inducing set can never decrease the lower bound in eq. (33).
Proof: For the added training input xi , we have to show that FV (Xm , xi ) ≥ FV (Xm ). Let fm+1 = (fm , fi )and Xm+1 = (Xm , xi ). The variational distribution qm (f ) that corresponds to the bound FV (Xm ) can be written18

in the formqm (f )=p(fn−(m+1) , fi |fm )φ∗m (fm ),=p(fn−(m+1) |fi , fm )p(fi |fm )φ∗m (fm ),(34)

where φ∗m (fm ) is the optimal choice for φm . This quantity has the same form with the variational distribution
em+1 (fm+1 ) = p(fi |fm )φ∗ (fm ). Notice that this φem+1 may not be optimal. Thus, by construction
qm+1 (f ) where φem+1 (fm+1 )) = FV (Xm ) with FV (Xm+1 , φem+1 ) computed by eq. (32). Computing the bound for theFV (Xm+1 , φ

optimal φ∗m+1 (fm+1 ) in order to maximize w.r.t that distribution we obtain

This completes the proof.

em+1 (fm+1 )) = FV (Xm ).
FV (Xm+1 ) ≥ FV (Xm+1 , φNote that strict inequality can hold when the optimal φ∗m+1 (fm+1 ) differs thanp(fi |fm )φ∗m (fm ).

C

Variational reformulation of the SPGP model

The input-constant noise GP model with likelihood p(y|f ) = N (y|f , σ 2 ) can be restrictive when the actual noise of
the data varies across the input space. A sparse GP method that can model input-dependent noise was proposed bySnelson and Ghahramani (2006). Assume a set of pseudo-inputs Xm and inducing latent function values fm , where−1
fm is drawn from the GP prior. The covariance matrix of the conditional prior p(f |fm ), i.e. Knn − Knm Kmm
Kmn ,can have non-stationary properties although the kernel function of the GP prior is stationary. The diagonal ofthis conditional matrix can be used to parametrize an input-dependent variance. The likelihood takes the formp(y|f ) = N (y|f , Λ),

(35)

−1
where Λ = σ 2 + diag[Knn − Knm Kmm
Kmn ]. The likelihood parameters are (σ 2 , Xm , θ) where θ are the kernelparameters of the GP prior. Note that the parameters θ are common for the prior and the likelihood.
The marginal likelihood can be written as
p(y)

=
=

Z

Z

p(y|f )p(f )df
p(y|f )p(f |fm )p(fm )df dfm(36)

where we augmented the GP prior on the training latent function f with the latent function values fm evaluatedat the pseudo-inputs Xm . We assume the variational distribution
q(f , fm ) = p(f |fm )φ(fm ).

(37)

Using this distribution we can bound the true log marginal likelihood and obtain the lower bound:
FV (Xm , φ)==Zp(y|f )p(fm )p(f |fm )φ(fm ) log
df dfmφ(fm )ﬀZZp(fm )dfm .
p(f |fm ) log p(y|f )df + log
φ(fm )
φ(fm )

19(38)

The integral involving f is computed as follows
Z
log G(fm , y) =p(f |fm ) log p(y|f )df

Z
h
iﬀ
n
1
1
=
p(f |fm ) − log(2π) − log |Λ| − T r Λ−1 (yyT − 2yf T + ff T ) df
2
2
2
i
h
1
1
n
−1
= − log(2π) − log |Λ| − T r Λ (yyT − 2yαT + ααT + Knn − Qnn ,
2
2
2

(39)

−1
−1
where α = E[f |fm ] = Knm Kmm
fm and Qnn = Knm Kmm
Kmn . The above quantity can also be written aslog G(fm , y) = log [N (y|α, Λ)] −1T r[Λ−1 (Knn − Qnn )].
2(40)

FV (Xm , φ) is now written as
FV (Xm , φ) =Z

φ(fm ) log

N (y|α, Λ)p(fm )1dfm − T r[Λ−1 (Knn − Qnn )].
φ(fm )
2

We can maximize this bound by reversing the Jensen’s inequality:
Z
1
FV (Xm ) = log N (y|α, Λ)p(fm )dfm − T r[Λ−1 (Knn − Qnn )]2
1= log [N (y|0, Λ + Qnn )] − T r[Λ−1 (Knn − Qnn )].
2

(41)

(42)

To compute the optimal distribution φ∗ we differentiate eq. (41) with respect to φ(fm ) and set to zero. This givesφ∗ (fm )∝=N (y|α, Λ)p(fm )

ﬀ
1 T
−1
−1
−1
−1
c exp − fm
(Kmm
+ Kmm
Kmn Λ−1 Knm Kmm)fm + yT Λ−1 Knm Kmm
fm .
2

(43)

From the above expression we can recognize the Gaussianφ∗ (fm ) = N (fm |Kmm Σ−1 Kmn Λ−1 y, Kmm Σ−1 Kmm )(44)

where Σ = Kmm + Kmn Λ−1 Knm . This distribution is exactly the one used by the SPGP model of Snelson andGhahramani (2006). Particularly, when we do predictions assuming the exact test conditional p(f∗ |fm ), the above
framework gives rises to the FITC prediction. The difference of the variational formulation with the SPGP/FITC
method is that the former starts with an input-dependent noise full GP model and derives a variational boundfor this model. The first term of the lower bound is the log marginal likelihood used by the SPGP model and thesecond term involves the trace of Λ−1 (Knn − Qnn ).

20