AutoGP: Exploring the Capabilities and Limitations of Gaussian Process Models

INTRODUCTION

Recent advances in deep learning (dl; LeCun et al., 2015) have revolutionized the application of machine learning in areas such as computer vision (Krizhevsky et al., 2012), speech recognition (Hinton et al., 2012) and natural language processing (Collobert and Weston, 2008). Although certain kernel-based methods have also been successful in such domains (Cho and Saul, 2009; Mairal et al., 2014), it is still unclear whether these methods can indeed catch up with the recent dl breakthroughs.

Aside from the benefits obtained from using compositional representations, we believe that the main components contributing to the success of dl techniques are: (i) their scalability to large datasets and efficient computation via gpus; (ii) their large representational power; and (iii) the use of well-targeted objective functions for the problem at hand.

In the kernel world, Gaussian process (gp; Rasmussen and Williams, 2006) models are attractive because they are elegant Bayesian nonparametric approaches to learning from data. Nevertheless, besides the limitations intrinsic to local kernel machines (Bengio et al., 2005), it is clear that gp-based methods have not fully explored the desirable criteria highlighted above.

Firstly, with regards to (i) scalability, despite recent advances in inducing-variable approaches and variational inference in gp models (Titsias, 2009; Hensman et al., 2013, 2015a; Dezfouli and Bonilla, 2015), the study of truly large datasets in problems other than regression and the investigation of gpu-based acceleration in gp models are still under-explored areas. We note that these issues are also shared by non-probabilistic kernel methods such as support vector machines (svms; Scholkopf and Smola, 2001).

Furthermore, concerning (ii) their representational power, kernel methods have been plagued by the overuse of very limited kernels such as the squared exponential kernel, also known as the radial-basis-function (rbf) kernel. Even worse, in some gp-based approaches, and more commonly in svms, these have been limited to having a single length-scale (a single bandwidth in svm parlance) shared across all inputs instead of using different length-scales for each dimensions. The latter approach is commonly known in the gp literature as automatic relevance determination (ard; Mackay, 1994).

Finally, and perhaps most importantly, gp methods have largely ignored defining (iii) well-targeted objective functions, and instead focused heavily on optimizing the marginal likelihood (Rasmussen and Williams, 2006, Ch. 5). This includes variational approaches that optimize a lower bound of the marginal likelihood (Titsias, 2009). Although different approaches for hyperparameter learning in gp models are discussed by Rasmussen and Williams (2006, §5.4.2), and indeed carried out by Sundararajan and Keerthi (2001); Sundararajan et al. (2007); Sundararajan and Keerthi (2008) and more recently by Vehtari et al. (2016), their performance results are somewhat limited by their disregard for other directions of improvement (scalability and greater representational power) mentioned above.

In this work we push the capabilities of gp-models, while also investigating their limitations, by addressing the above issues jointly. In particular,

we develop scalable and statistically efficient inference methods for gp models using gpu computation and stochastic variational inference along with the reparameterization trick (Kingma and Welling, 2014);

we investigate the flexibility of models using ard kernels for high-dimensional problems, as well as “deep” arc-cosine kernels (Cho and Saul, 2009); and

we study the impact of employing a leave-one-out-based objective function on hyperparameter learning.

As we rely on automatic differentiation (Baydin et al., 2015) for the implementation of our model in TensorFlow (Abadi et al., 2015), we refer to our approach as autogp. While we thoroughly evaluate our claimed contributions on various regression and classification problems, our most significant results show that autogp:

has superior performance to all previously reported gp-based methods on the standard mnist dataset;

achieves comparable performance to previous kernel-based methods using the rectangles-image dataset; and

breaks the 1% error-rate barrier in gp models using the mnist8m dataset, showing along the way the scalability of our approach at unprecedented scale for gp models (8 million observations) in multiclass classification problems.

Overall, autogp represents a significant breakthrough in kernel methods in general and in gp models in particular, bridging the gap between deep learning approaches and kernel machines.

Related Work

The majority of works related to scalable kernel machines primarily target issues pertaining to the storage and computational requirements of algebraic operations involving kernel matrices. Low-rank approximations are ubiquitous among these approaches, with the Nyström approximation being one of the most popular methods in this category (Williams and Seeger, 2001). Nevertheless, most of these approximations can be understood within the unifying probabilistic framework of Quiñonero-Candela and Rasmussen (2005).

The optimization of inducing inputs using Nyström approximations for gps made significant progress with the work on sparse-gp variational inference in Titsias (2009). This exposed the possibility to develop stochastic gradient optimization for gp models (Hensman et al., 2013) and handle non-Gaussian likelihoods (Nguyen and Bonilla, 2014; Dezfouli and Bonilla, 2015; Hensman et al., 2015a; Sheth et al., 2015), which has sparked interest in the implementation of scalable inference frameworks such as those in Hensman et al. (2016) and Deisenroth and Ng (2015). These developments greatly improved the generality and applicability of gp models to a variety of applications, and our work follows from this line of ideas.

Other independent works have focused on kernel design (Wilson and Adams, 2013; Cho and Saul, 2009; Mairal et al., 2014; Wilson et al., 2016). This has sparked some debate as to whether kernel machines can actually compete with deep nets. Evidence suggests that this is possible; notable examples include the work by Huang et al. (2014); Lu et al. (2014). These works also provide insights into the aspects that make kernel methods less competitive to deep neural networks, namely lack of application specific kernels, and scalability at the price of poor approximations. These observations are corroborated by our experiments, which report how the combination of these factors impacts performance of gps.

Because gps are probabilistic kernel machines, it is natural to target the optimization of the marginal likelihood. However, alternative objective functions have also been considered, in particular the leave-one-out cross-validation (loo-cv) error. Originally this was done by Sundararajan and Keerthi (2001) for regression; later extended by Sundararajan et al. (2007) to deal with sparse gp formulations; and broadened by Sundararajan and Keerthi (2008); Vehtari et al. (2016) to handle non-Gaussian likelihoods. While the results in these papers seem inconclusive as to whether this can generally improve performance, this may allow to better accommodate for incorrect model specifications (Rasmussen and Williams, 2006, §5.4.2). Motivated by this observation, we explore this direction by extending the variational formulation to optimize a loo-cv error. We also note that none of these previous loo-cv approaches can actually handle large datasets.

Gaussian Process Models

We are interested in supervised learning problems where we are given a training dataset of input-output pairs $\mathcal{D}=\{\mathbf{x}_{n},\mathbf{y}_{n}\}_{n=1}^{n}$, with $\mathbf{x}_{n}$ being a $D$-dimensional input vector and $\mathbf{y}_{n}$ being a $P$-dimensional output. Our main goal is to learn a probabilistic mapping from inputs to outputs so that for a new input $\mathbf{x}_{\star}$, we can estimate the probability of its associated label $p(\mathbf{y}_{\star}|\mathbf{x}_{\star})$.

To this end, we follow a Gaussian process modeling approach (gp; Rasmussen and Williams, 2006), where latent functions are assumed to be distributed according to a gp, and observations are modeled via a suitable conditional likelihood given the latent functions. A function $f_{j}$ is said to be distributed according to a Gaussian process with mean function $\mu_{j}(\mathbf{x})$ and covariance function $\kappa_{j}(\mathbf{x},\mathbf{x}^{\prime};\boldsymbol{\theta}_{j})$, which we denote with $f_{j}\sim\mathcal{GP}(\mu_{j}(\mathbf{x}),\kappa_{j}(\mathbf{x},\mathbf{x}^{\prime};\boldsymbol{\theta}_{j}))$, if any subset of function values follow a Gaussian distribution.

Since we are dealing with multiple outputs $\{\mathbf{y}_{n}\}$, we follow the standard practice of considering $Q$ underlying latent functions $\{f_{j}\}_{j=1}^{Q}$ which are drawn independently from zero-mean gp priors along with i.i.d. conditional likelihood models. Such a modeling approach encompasses a large class of machine learning problems, including gp-regression, gp-classification and modeling of count data (Rasmussen and Williams, 2006).

Consequently, when realized at the observed data, our probabilistic model is given by:

One commonly used covariance function is the so-called squared exponential (or rbf kernel):

The main inference problems in gp models are (i) estimation of the posterior over the latent functions $p(\mathbf{f}|\mathbf{X},\mathbf{y},\boldsymbol{\theta})$ and subsequent estimation of the predictive distribution $p(\mathbf{y}_{\star}|\mathbf{x}_{\star},\mathbf{X},\mathbf{y},\boldsymbol{\theta})$; and (ii) inference over the hyperparameters $\boldsymbol{\theta}$.

For general likelihood models, both inference problems are analytically intractable, as they require the computation of non-trivial high-dimensional integrals. Furthermore, even for the simplest case when the conditional likelihood is Gaussian, evaluating the posterior and hyperparameter estimation is computationally prohibitive, as they entail time and memory requirements of $\mathcal{O}(N^{3})$ and $\mathcal{O}(N^{2})$ respectively.

In addition to dealing with general likelihood models and scalability issues, the performance of gp models (and kernel machines in general), is highly dependent on the covariance function (or kernel) used. Most previous work on gp methods and svms in the machine learning literature, use the squared exponential kernel in Equation (3). Furthermore, especially in high dimensions, this kernel is constrained to its isotropic version. As pointed out by Bengio et al. (2005), such an approach is severely limited by the curse of dimensionality when learning complex functions.

Finally, the estimation of hyper-parameters based solely on marginal likelihood optimization can be very sensitive to model misspecification, with predictive approaches such as those in Sundararajan and Keerthi (2001) being seemingly more appealing for this task, especially in problems such as classification, where we are ultimately interested in having lower error rates and better calibrated predictive probabilities.

In the following sections, we show how to deal with the above issues, namely non-Gaussian likelihood models and scalability; more flexible kernels; and better objective functions for hyperparameter learning, with the endmost goal of improving performance significantly over current gp approaches.

Automated Variational Inference

In this section we detail our method for scalable and statistically efficient inference in Gaussian process models with general likelihoods as specified by Equations (1) and (2). We adopt the variational framework of Dezfouli and Bonilla (2015), which we prefer over the mcmc method of Hensman et al. (2015b) as we avoid the sampling overhead which is especially significant in very large datasets.

In order to have a scalable inference framework, we augment our prior with $M$ inducing variables $\{\mathbf{u}_{\cdot j}\}$ per latent process and corresponding inducing inputs $\{\mathbf{Z}_{j}\}$ such that,

where $\mathbf{K}_{\mathbf{u}\mathbf{v}}^{j}$ is the covariance matrix obtained by evaluating the covariance function $\kappa_{j}(\cdot,\cdot;\boldsymbol{\theta})$ at all pairwise columns of matrices $\mathbf{U}$ and $\mathbf{V}$. We note that this prior is equivalent to that defined in Equation (1), which is obtained by integrating out $\mathbf{u}$ from the joint $p(\mathbf{f},\mathbf{u})$. However, as originally proposed by Titsias (2009), having an explicit representation of $\mathbf{u}$ will allow us to derive a scalable variational framework without additional assumptions on the train or test conditional distributions (Quiñonero-Candela and Rasmussen, 2005).

We now define our approximate joint posterior distribution as:

where $p(\mathbf{f}|\mathbf{u})$ is the conditional prior given in Equation (4) as variational parameters. $q(\mathbf{u}|\boldsymbol{\lambda})$ is our variational posterior with $\boldsymbol{\lambda}=\{\pi_{k},\mathbf{m}_{kj},\mathbf{S}_{kj}\}$. Note that we assume a variational posterior in the form of a mixture for added flexibility compared to using a single distribution.

2 Evidence Lower Bound

Such a definition of the variational posterior allows for simplification of the log-evidence lower bound ($\mathcal{L}_{\text{elbo}}$) such that no $\mathcal{O}(N^{3})$ operations are required,

In order to compute the log-evidence lower bound in Equation (4.2) and its gradients, we use Jensen’s inequality to bound the KL term (as in Dezfouli and Bonilla, 2015) and estimate the expected likelihood term using Monte Carlo (mc).

3 The Reparameterization Trick

One of the key properties of the log-evidence lower bound in Equation (4.2) is its decomposition as a sum of expected likelihood terms on the individual observations. This allows for the applicability of stochastic optimization methods and, therefore, scalability to very large datasets.

In contrast, if one is willing to relax the constraints of a truly black-box approach by having access to the implementation of the conditional likelihood and its gradients, then a simple trick proposed by Kingma and Welling (2014) can significantly reduce the variance in the gradients and make the stochastic optimization of $\mathcal{L}_{\text{elbo}}$ much more practical. This has come to be known as the reparameterization trick.

We now present explicit expressions of our estimated $\mathcal{L}_{\text{elbo}}$, focusing on an individual expected likelihood term in Equation (4.2). Thus, the individual expectations can be estimated as:

where $S$ is the number of mc samples and $f^{(k,i)}_{nj}$ is an element in the vector of latent functions $\mathbf{f}_{n\cdot}^{(k,i)}$.

4 Mini-Batch Optimization

We can now obtain unbiased estimates of the gradients of $\mathcal{L}_{\text{elbo}}$ to use in mini-batch stochastic optimization. In particular, for a mini-batch $\Omega$ of size $B$ we have that the estimated gradient is given by:

where $\boldsymbol{\eta}\in\{\boldsymbol{\lambda},\boldsymbol{\theta}\}$; $\boldsymbol{\theta}$ now includes not only the covariance parameters, but also the inducing inputs across all latent processes $\mathbf{Z}_{j}$. The corresponding gradients are obtained through automatic differentiation using TensorFlow.

5 Full Approximate Posterior

An important characteristic of our approach is that, unlike most previous work that uses stochastic optimization in variational inference along with the reparameterization trick (see e.g. Kingma and Welling, 2014; Dai et al., 2016), we do not require a fully factorized approximate posterior to estimate our gradients using only simple scalar transformations of uncorrelated Gaussian samples, as given in Equation (14). Indeed, the posterior over the latent functions is fully correlated and is given in the Appendix. This result is a property inherited from the original framework of Nguyen and Bonilla (2014) and Dezfouli and Bonilla (2015), who showed that, even when having a full approximate posterior, only expectations over univariate Gaussians are required in order to estimate the expected log likelihood term in $\mathcal{L}_{\text{elbo}}$.

Flexible Kernels

We now turn our attention to increase the flexibility of gp models through the investigation of kernels beyond the isotropic rbf kernel. This is a complementary direction for exploring the capabilities of gp models, but it benefits from the results of the previous section in terms of large-scale inference and computation of gradients via (automatic) back-propagation.

For the rbf kernel, we limit our experiments to the automatic relevance determination (ard) setting in Equation (3), including problems of large input dimensionality such as mnist and cifar10. Furthermore, given the interesting results reported by Cho and Saul (2009), we also investigate their arc-cosine kernel, which was proposed as a kernel that mimics the computation of neural networks. We give the mathematical details of the arc-cosine kernel in the Appendix.

Leave-One-Out Learning

Here we focus on the average leave-one-out log predictive probability for hyper-parameter learning, as an alternative to optimization of the marginal likelihood. This can be particularly useful in problems such as classification where one is mainly interested in having lower error rates and better calibrated predictive probabilities.

This objective function is obtained by leaving one datapoint out of the training set and computing its log predictive probability when training on the other points; the resulting values are then averaged across all the datapoints. Interestingly, in our gp model this can be computed without the need for training $N$ different models, as the resulting expression is given by:

where the approximation stems from using the variational marginal posterior $q(\mathbf{f}_{n\cdot}|\mathcal{D},\boldsymbol{\theta})$ instead of the true marginal posterior $p(\mathbf{f}_{n\cdot}|\mathcal{D},\boldsymbol{\theta})$ for datapoint $n$. We also note that we have made explicit the dependency of the posterior on all the data. The derivation of this expression is given in the Appendix. Since $\mathcal{L}_{\text{oo}}(\boldsymbol{\theta})$ contains an expectation with respect to the marginal posterior we can also estimate it via mc sampling.

Equation (16) immediately suggests an alternating optimization scheme where we estimate the approximate posterior $q(\mathbf{f}_{n\cdot}|\mathcal{D},\boldsymbol{\theta})$ through optimization of $\mathcal{L}_{\text{elbo}}$, as described in §4, and then learn the hyperparameters via optimization of $\mathcal{L}_{\text{oo}}$. Algorithm 1 illustrates such an alternating scheme for hyperparameter learning in our model using the leave-one-out objective and vanilla stochastic gradient descent (sgd), where $\widehat{\mathcal{L}}_{\text{elbo}}$ and $\widehat{\mathcal{L}}_{\text{oo}}$ denote estimates of the corresponding objectives when using mini-batches.

Experiments

We evaluate autogp across various datasets with number of observations ranging from $12,000$ to $8.1$M and input dimensionality between $31$ and $3072$. The aim of this section is (i) to demonstrate the effectiveness of the re-parametrization trick in reducing the number of samples needed to estimate the expected log likelihood; (ii) to evaluate non-isotropic kernels across a large number of dimensions; (iii) to assess the performance obtained by loo-cv hyperparameter learning; and (iv) to analyze the performance of autogp as a function of time. Details of the experimental set-up can be found in the supplementary material.

For our first experiment we look at the behavior of our model as we vary the number of samples used to estimate the expected log likelihood. We trained our model on the sarcos dataset (Vijayakumar and Schaal, 2000), an inverse dynamics problem for a seven degrees-of-freedom anthropomorphic robot arm.

We used the Gaussian process regression network (gprn) model of Wilson et al. (2012) as our likelihood function. Bonilla et al. (2016) found that gprns require a large number of samples ($10,000$) to yield stable optimization of the $\mathcal{L}_{\text{elbo}}$ in their method. This is most likely due to the increased variance induced by multiplying latent samples together. The high variance of this likelihood model makes it an ideal candidate for us to evaluate how the performance of autogp varies with the number of samples.

As shown in Figure 1, increasing the number of samples decreased the number of epochs it took for our model to converge; however, our model always converged to the same optimal value. This is in stark contrast to savigp (Bonilla et al., 2016) which was able to converge to the optimal value of $0.0195$ with $10,000$ samples, but converged to sub-optimal values ($>0.05$) as the number of samples was lowered. This shows that our mc estimator is significantly more stable than savigp’s black-box estimator, without requiring any variance reduction methods. In practice, reducing the number of samples leads to an improved runtime. A gradient update using $10,000$ samples takes around $0.17$ seconds, whereas an update using $10$ samples only takes around $0.03$ seconds, which makes up for the extra epochs needed to converge.

2 Kernel Performance

In this section we compare the performance of autogp across two different kernels on a high-dimensional dataset. We trained our model on the rectangles-image dataset (Larochelle et al., 2007) which is a binary classification task that was created to compare shallow models (e.g. svms), and deep learning architectures. The task involves recognizing whether a rectangle has a larger width or height. Randomized images were overlayed behind and over the rectangles, which makes this task particularly challenging.

We compare the multilayer arc-cosine kernel (arc-cosine) described by Cho and Saul (2009) with an rbf kernel with automatic relevance determination (rbf-ard). arc-cosine was devised to mimic the architecture of deep neural networks through successive kernel compositions, and has shown good results when combined with non-Bayesian kernel machines (Cho and Saul, 2009). Unlike arc-cosine, rbf-ard is commonly used with Gaussian processes (Hensman et al., 2013). However, rbf-ard has not been applied to large-scale datasets with a large input dimensionality due to limitations in scalability. Given that our implementation uses automatic differentiation, and is thus able to efficiently compute gradients, we can use rbf-ard with no noticeable performance overhead.

We trained our model using $10$, $200$, and $1000$ inducing points across both kernels. We experimented with various different depths and degrees for arc-cosine and found that $3$ layers and a degree of $1$ worked best with autogp. As such, we ran all arc-cosine experiments with these settings.

For reference, Cho and Saul (2009) report an error rate of $22.36\%$ using svms with arc-cosine. Larochelle et al. (2007) report error rates of $24.04\%$ on svms with a simple rbf kernel and $22.50\%$ with a deep belief network with $3$ hidden layers (DBN).

Figure 2 shows that autogp is competitive with previous approaches, except when using $10$ inducing points. Our results also show that rbf-ard performs comparatively (and sometimes better) than arc-cosine. This is likely because arc-cosine is better suited to deep architectures, such as the multilayer kernel machine of Cho and Saul (2009).

3 LOO-CV Hyperparameter Learning

In this section we analyze our approach when learning hyperparameters via optimization of the leave-one-out objective ($\widehat{\mathcal{L}}_{\text{oo}}$), as described in §6, and compare it with the standard variational method that learns hyperparameters via sole optimization of the variational objective ($\widehat{\mathcal{L}}_{\text{elbo}}$). To this end, we use the mnist dataset using $10$, $200$, and $1000$ inducing points and the rbf-ard kernel across all settings.

Our results indicate that optimizing the leave-one-out objective leads to a significant performance gain. While Figure 3 shows the corresponding error rates, we refer the reader to the supplementary material for similar results on negative log probabilities (nlps). In summary, at $1000$ inducing points our loo-cv approach attains an error rate of $1.55\%$ and a mean nlp of $0.061$; At $200$ inducing points it obtains an error rate of $1.71\%$ and a mean nlp of $0.063$. This represents a $40\%$–$45\%$ decrease in mean nlp and error rate from the equivalent model that only optimized the variational objective. Furthermore, our results outperform the work of Hensman et al. (2015b), Gal et al. (2014), and Bonilla et al. (2016) who reported error rates of $1.96\%$, $5.95\%$ and $2.74\%$ respectively, that are the best results reported on this data set in the literature of GPs.

We have shown that our approach reduces the error rate on the test set down to $1.55\%$ without artificially extending the mnist dataset through transformations, or doing any kind of pre-processing of the data. These results, while not on par with state-of-the-art vision results, show that Gaussian process models can perform competitively with non-Bayesian methods, while solving the harder problem of full posterior estimation. For example, standard convolutional neural networks achieve an error rate of $1.19\%$ without any pre-processing (Benenson, 2013). We will show in the following section how we can further improve our results by artificially extending the training data on mnist through distortions and affine transformations.

4 AutoGP All Together

We finish by evaluating the performance of autogp using all its main components simultaneously, namely statistical efficiency, flexible kernels and loo-cv hyperparameter learning, using three datasets: rectangles-image, cifar10, and mnist8m. Our goal is to gain further insights as to how our model compares against other works in the machine learning community. As before, we optimized our model with $200$ inducing points and used the rbf-ard kernel across all experiments.

We begin by looking at the rectangles-image dataset which we have already investigated in Section 7.2. Optimizing $\widehat{\mathcal{L}}_{\text{oo}}$ slightly improves performance in this case, with a reduction in error rate from $24.06\%$ to $24.04\%$ and a reduction in mean nlp from $0.495$ to $0.485$, when using 200 inducing points. The error rate can be further reduced to $23.6\%$ when using 1,000 inducing points (with corresponding nlp of 0.478). As a reference, performance of recent deep architectures is: deep trans-layer autoencoder networks (${\textsc{er}}=13.01\%$, Zhu et al., 2015), PCANet (${\textsc{er}}=13.94\%$, Chan et al., 2014), and two-layer stacked contractive auto-encoders (${\textsc{er}}=21.54\%$, Rifai et al., 2011). Bruna and Mallat (2013) have reported results using invariant scattering convolutional networks with an error rate of $8.01\%$.

We now turn our attention to the mnist8m (Loosli et al., 2007) dataset, which artificially extends the mnist dataset to $8.1$ million training points by pseudo-randomly transforming existing mnist images. To the best of our knowledge, no other gp models have been shown to scale to such a large dataset. We were able to achieve an error rate of $0.89\%$ and a mean nlp of $0.033$, breaking the $1\%$ error-rate barrier in gp models. This shows that our model is able to perform on par with deep architectures with some pre-processing of the mnist dataset. We also note that we were able to rapidly converge given the dataset size. As shown in Figure 4 autogp achieves near optimal results within the first hour, and converges to the optimal value within 10 hours.

Finally, we look at cifar10, a multiclass image classification task and we find that results are less conclusive. We achieve an error rate of $44.95\%$ and a mean nlp of $1.33$, which is well below the state of the art, with some models achieving less than $10\%$ error rate. However, we note that we still compare favorably to svms which have been shown to achieve error rates of $57.7\%$ (Le et al., 2013). We also perform on par with older deep architectures (Rifai et al., 2011). We believe these results provide a strong motivation for investigating application-specific kernels (such as those studied by Mairal et al., 2014) in gp models.

CONCLUSIONS & DISCUSSION

We have developed autogp, an inference framework for gp models that pushes the performance of current gp methods by exploring three complementary directions: (i) scalable and statistically efficient variational inference; (ii) flexible kernels; and (iii) objective functions for hyperparameter learning alternative to the marginal likelihood.

We have shown that our framework outperforms all previous gp approaches on mnist; achieves comparable performance to previous kernel approaches on the rectangles-image dataset; and can break the $1\%$ error barrier using mnist8m. Overall, this represents a significant breakthrough in Gaussian process methods and kernel machines. While our results on cifar10 are well below the state-of-the-art achieved with deep learning, we believe we can further reduce the gap between kernel machines and deep architectures by using application-specific kernels such as those recently proposed by Mairal et al. (2014).

Appendix A Details of Full Posterior Distribution

The full posterior distribution over the latent functions is given by:

Appendix B Details of The Arc-Cosine Kernel

Cho and Saul (2009) define an arc-cosine kernel of degree $d$ and depth $0pt$ using the following recursion:

Appendix C Derivation of Leave-One-Out Objective

In this section we derive an expression for the leave-one-out objective and show that this does not require training of $N$ models. A similar derivation can be found in Vehtari et al. (2016). Let $\mathcal{D}_{\neg n}=\{\mathbf{X}_{\neg n},\mathbf{y}_{\neg n}\}$ be the dataset resulting from removing observation $n$. Then our leave-one-out objective is given by:

We now that the marginal posterior can be computed as:

and substituting this expression in Equation (25) we have

We see that the objective only requires estimation of the marginal posterior $p(\mathbf{f}_{n\cdot}|\mathcal{D},\boldsymbol{\theta})$, which we can approximate using variational inference, hence:

where $q(\mathbf{f}_{n\cdot}|\mathcal{D},\boldsymbol{\theta})$ is our approximate variational posterior.

Appendix D Additional Details of Experiments

The datasets used are described in Table 2.

We trained our model stochastically using the rmsprop optimizer provided by TensorFlow (Abadi et al., 2015) with a learning rate of $0.003$ and mini-batches of size $1000$. We initialized inducing point locations by using the k-means clustering algorithm, and initialized the posterior mean to a zero vector, and the posterior covariances to identity matrices. When jointly optimizing $\widehat{\mathcal{L}}_{\text{oo}}$ and $\widehat{\mathcal{L}}_{\text{elbo}}$, we alternated between optimizing each objective for $100$ epochs. Unless otherwise specified we used $100$ Monte-Carlo samples to estimate the expected log likelihood term.

All timed experiments were performed on a machine with an Intel(R) Core(TM) i5-4460 CPU, 24GB of DDR3 RAM, and a GeForce GTX1070 GPU with TensorFlow 0.10rc.

D.2 Additional Results

Figure 5 shows the nlp for our evaluation of the loo-cv-based hyperparameter learning. As with the error rates described in the main text, the nlp obtained with loo-cv are significantly better than those obtained with a purely variational approach.

References