1d examples . . . . . . .

Introduction

Overview

Machine Learning as an inductive problem. Machine Learning (ml) is seen as an application of Artificial Intelligence. The use of learning algorithms equips systems with the ability to automatically acquire helpful information from experience without being explicitly programmed. Depending on the contextual scenarios, ml has been categorized into different approaches, e.g. supervised, unsupervised, semi-supervised or by reinforcement. In this dissertation, we mainly focus on supervised learning problems where pairs of input and outputs are collected to learn a mapping functions from an input space to an output space. From the available observations, we wish to derive a function that models the underlying mapping from the input data (covariates) to labels (or target values); from the function, we can then make predictions for all possible input values. It is obvious that the problem at hand is inductive.

The approaches for learning the mappting function in a given task can be grouped into two categories: parametric and non-parametric.

Parametric Modeling. Traditionally, we can use parametric machine learning algorithms to deal with supervised problems. This kind of modeling restricts the underlying mapping to a family of functional forms which is parameterized by a finite set of parameters. It also implies that no matter how much data is fed to a parametric model, it will not change its mind about how many parameters it needs [START_REF] Russell | Artificial Intelligence: A Modern Approach[END_REF]. Such parametric models often perform inefficiently if the functional form is inadequate to represent the actual unknown underlying correlation between inputs and its labels. One may be tempted to employ a flexible functional form, e.g. we can assume the parametric function is the one obtained by a neural network, but this runs into the danger of overfitting, so that we can obtain a good fit to training data, but perform badly in predictions.

Introduction

Non-parametric Modeling. In contrast to parametric model accompanying with a specified functional form, algorithms using free-form mapping functions are classified as non-parametric machine learning algorithms, such as k-Nearest Neighbor [START_REF] Cover | Nearest neighbor pattern classification[END_REF], Decision Trees [START_REF] Quilan | Decision Trees and Multi-Valued Attributes[END_REF] or Support Vector Machine [START_REF] Cortes | Support Vector Networks[END_REF] or Kernel methods [START_REF] Hofmann | Kernel methods in machine learning[END_REF]. Non-parametric feature extraction algorithms have more advantages than parametric ones and are well suited for non-normally distributed data along with being able to extract more features than the classic linear discriminant analysis [START_REF] Russell | Artificial Intelligence: A Modern Approach[END_REF][START_REF] Yang | A nonparametric feature extraction and its application to nearest neighbor classification for hyperspectral image data[END_REF].

In general, such non-parametric models possessing an infinite set of parameters are capable of fitting any complicated functional form. Nevertheless, it also implies that the number of labeled data required by non-parametric approaches to estimate the mapping function is greater than the parametric model with a finite set of parameters. Therefore, non-parametric models are easy prone to overfitting, especially when labeled data is scarce.

Scalability of Non-parametric approaches. In the era of big data, nonparametric models are promising solutions allowing to learn complicated patterns from data. Nevertheless, the computational complexity of non-parametric approaches depends on the training size. For example, the training phase of Kernel Support Vector Machine [START_REF] Cortes | Support Vector Networks[END_REF] involves solving a quadratic problem which generally suffers cubic time complexity with respect to data size. Consider K-Nearest Neighbor [START_REF] Cover | Nearest neighbor pattern classification[END_REF] as another example; it is a non-parametric lazy learning algorithm which does not require an explicit training phase. However, K-Nearest Neighbor makes prediction on unseen data as a vote by using all the training data. Generally, the testing phase of these methods requires linear time complexity to data size. Hence, the application of non-parametric models to large-scale problems is hindered by their poor scalability.

Needs of Predictive Uncertainty. The problem of enhancing the safety of decision-making system by acting on the model's prediction in an informed manner has obtained a significant attention from the machine learning community [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF]. Predictive uncertainty quantification has a crucial role to strengthen the safety of an AI system [START_REF] Amodei | Concrete problems in ai safety[END_REF] by acting on the model's prediction in an informed manner. This is essential to applications where the consequence of an error is serious, such as in autonomous vehicle control and medical, financial and legal fields. Hence, accurate fitting capabilities are no longer the most important aspects for evaluating the model's effectiveness. Evaluation of predictive uncertainty. Evaluating the quality of predictive uncertainties is challenging as the ground-truth uncertainty estimates are unknown. Being motivated by practical applications, there are two aspects that are able to examine the plausibility of predictive uncertainty. The first notion of quality of predictive uncertainty concerns calibration [START_REF] Dawid | The well-calibrated bayesian[END_REF][START_REF] Degroot | The comparison and evaluation of forecasters[END_REF], which measures the discrepancy between subjective forecast and (empirical) long-run frequencies. Traditionally, the quality of calibration can be numerically assessed by proper scoring rules [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF], such as the Brier score [START_REF] Brier | VERIFICATION OF FORECASTS EXPRESSED IN TERMS OF PROBABILITY[END_REF]. Secondly, the quality of predictive uncertainty is also obtainable using out-of-distribution examples [START_REF] Hendrycks | A baseline for detecting misclassified and outof-distribution examples in neural networks[END_REF]. For example, if a model is trained on one dataset, but is tested on a completely different dataset, the predictive uncertainty returned by the model should be high, as testing points would be distant from training points. Recently, the works of approximation of predictive uncertainty based upon ensemble learning are robust to calibration as well as the scenarios of data shift (Lakshminarayanan et al., 2017). Alternatively, a plethora of works revolves around the Bayesian formalism [START_REF] Bernardo | Bayesian Theory[END_REF] with the aim of adapting neural networks to encompass predictive uncertainty and give them a probabilistic flavor [START_REF] Mackay | Bayesian Methods for Adaptive Models[END_REF]Graves, 2011;[START_REF] Liu | Structured and efficient variational deep learning with matrix gaussian posteriors[END_REF][START_REF] Blundell | Weight Uncertainty in Neural Network[END_REF].

Gaussian Processes. As alluded earlier, an ideal modeling approach in the era of big data should possess not only a powerful fitting capability but also a firm mechanism to determine predictive uncertainty. Bayesian nonparametric approaches are ideal candidates due to their advantages over flexibility and calibrated predictive uncertainty. The philosophies and motivations of this area have been well discussed by a number of authors [START_REF] Hjort | Bayesian Nonparametrics. Cambridge Series in Statistical and Probabilistic Mathematics[END_REF][START_REF] Ghosh | Bayesian nonparametrics[END_REF][START_REF] Ghahramani | Bayesian non-parametrics and the probabilistic approach to modelling[END_REF]. Gaussian Processes (gps) [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] are an attractive way of doing non-parametric Bayesian modeling. A Gaussian Process is a collection of 1. Introduction random variables indexed by a variable in the input domain, such that every subset of those random variables has a multivariate normal distribution.

Thanks to the properties of the multivariate normal distribution, given observations, gps are able to make inferences as well as predictive uncertainties with a firm mathematical background. In addition to providing uncertainty in predictions, there are also compelling reasons to use gps, such as the gps can represent a rich family of functions; also, gps are protected from overfitting with an appropriate prior on hyperparameters. In practice, gps achieve stateof-the-art results in a wide spectrum of applications including robotics [START_REF] Ko | Gp-bayesfilters: Bayesian filtering using gaussian process prediction and observation models[END_REF][START_REF] Deisenroth | Pilco: A model-based and dataefficient approach to policy search[END_REF], geostatistics [START_REF] Diggle | Model-based Geostatistics[END_REF], numerics [START_REF] Briol | Probabilistic integration: A role in statistical computation?[END_REF], active sensing [START_REF] Guestrin | Near-optimal sensor placements in gaussian processes[END_REF] and optimization [START_REF] Snoek | Practical bayesian optimization of machine learning algorithms[END_REF].

Deep Gaussian Processes. A shallow gp is defined by a mean and covariance/kernel function. Kernel functions hold a crucial role as it not only encodes our assumptions as well as the desired flexibility into the functions we wish to learn. Thus, enhancing the expressiveness of kernel functions are able to boost the gps' power. A Deep Gaussian Process dgp [START_REF] Damianou | Deep Gaussian Processes[END_REF] which is a hierarchical composition of multiple gps, comes to a rescue of the limitation of the representational power of a single-layer gp.

dgps is more flexible than a standard gp, just as deep neural networks are more powerful than a Multilayer Perceptron with one hidden layer. In contrast to models constructed by with a highly parameterized functional form, dgps learn a hierarchical representation with very few hyperparameters to optimize.

Extensions and Open Problems

In this section, I introduce the extensions and open problems of (Deep) Gaussian Processes which will appear in the dissertation.

Combination of Neural Networks and Gaussian Processes. In 1996, [START_REF] Neal | Bayesian Learning for Neural Networks (Lecture Notes in Statistics)[END_REF] showed that Bayesian Neural Networks with infinitely many hidden units converged to Gaussian Processes (gps) with a particular kernel function. Speaking theoretically, Gaussian Processes were viewed as an interpretable alternative to neural networks. However, in practice, the power of gps are restricted by the limitations of the kernel function. By contrast, neural networks are able to automatically discover meaningful representations in high-dimensional data by learning multiple layers of highly adaptive basis functions MacKay (1998); [START_REF] Hinton | A fast learning algorithm for deep belief nets[END_REF]; [START_REF] Bengio | Learning deep architectures for ai[END_REF].

Extensions and Open Problems

Despite the impressive expressiveness, neural networks access predictive uncertainties via sampling using approaches [START_REF] Welling | Bayesian Learning via Stochastic Gradient Langevin Dynamics[END_REF]; Gal and Ghahramani (2016a); Lakshminarayanan et al. (2017). Unlike neural networks, gps directly capture predictive uncertainties with a firm mathematical background. Another advantage of gps over neural networks is that the prior knowledge about the properties of mapping function, e.g. smoothness, differentiability or periodicity, can be added by specifying an appropriate kernel function.

As neural networks and gps have particular strengths, the question of what the best paradigm, e.g. kernel methods in general (Gaussian Processes in particular) and neural networks) is become irrational. Instead, it is more sensible to think about the idea of combining the advantages of each approach. There are several works about the combinations of convolutional neural networks and gps on image recognition, e.g. substituting gps for the last fully connected layers [START_REF] Bradshaw | Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks[END_REF]; Wilson et al. (2016) or introducing convolutions in the calculation of the covariance between images van der [START_REF] Van Der Wilk | Convolutional Gaussian Processes[END_REF].

Evaluation of Predictive Uncertainty of Probabilistic Models. As alluded in the introduction section, predictive uncertainty can be evaluated by inspecting the calibration and out-of-distribution samples. The majority of works accessing predictive uncertainty on NN involve with Bayesian formalism [START_REF] Mackay | Bayesian Methods for Adaptive Models[END_REF]; Graves (2011); [START_REF] Liu | Structured and efficient variational deep learning with matrix gaussian posteriors[END_REF]; [START_REF] Blundell | Weight Uncertainty in Neural Network[END_REF]. Along a similar vein, combining cnn and gp is an intuitive way to add probabilistic flavor to cnn [START_REF] Bradshaw | Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks[END_REF]; Wilson et al. (2016); van der [START_REF] Van Der Wilk | Convolutional Gaussian Processes[END_REF]. Intuitively, the motivation to impose these Bayesian treatments into neural networks is to do a better quantification of uncertainty compared to plain neural networks. Nevertheless, analyzing Bayesian Neural Networks and the combination of neural networks and gps on predictive uncertainty has not been conducted carefully.

Inducing point-based approximation. GPs [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] are well-known because of the predictive uncertainties with a firm mathematical background. Despite being able to underpin a range of algorithms for supervised and unsupervised learning, the application of gps is hindered to the large-scale problems due to the burden of computational and storage cost. Assuming that the input dimensionality D is significantly less than the number of observations N , gps require the complexities of O (N 3 ) and O (N 2 ) for computation and storage. These costs are sourced from linear algebraic operation with the N × N kernel matrix. To improve the scalability 1. Introduction of gps, we must employ a technique accelerating the computation involving the kernel matrix. Almost works discussing the scalable gp have focused on the low-rank approximation of kernel matrix using inducing points [START_REF] Lawrence | Fast Sparse Gaussian Process Methods: The Informative Vector Machine[END_REF]Seeger et al., 2003;[START_REF] Snelson | Sparse Gaussian Processes using Pseudoinputs[END_REF]Naish-Guzman and Holden, 2007;[START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF][START_REF] Hensman | Gaussian processes for big data[END_REF]Wilson and Nickisch, 2015;Hensman et al., 2015a). Using M inducing points to obtain an approximation to the kernel matrix, the computational and storage costs are contracted to O (M 3 ) and O (M 2 ) respectively. It is obvious that inducing point-based approaches lead to a remarkable development on the scalability of gps if M is significantly less than N .

Recently, it has been shown that it is possible to obtain an arbitrarily good approximation for a certain class of GP models (i.e. conjugate likelihoods, concentrated distribution for the training data) with M growing more slowly than N . However, the general case remains elusive and it is still possible that the required value for M may exceed a certain computational budget. To employ a large number of inducing points without exploding the computational cost, these inducing inputs are arranged into a structure such that the resulting kernel matrix allows for the application of fast linear algebra, and the entries of the kernel matrix evaluated at the training inputs are approximated through interpolation via sparse matrices. A well-known example for this line of work was introduced by Wilson et al Wilson and Nickisch (2015), namely Kernel Interpolation for Scalable Structured gps (kiss-gp). The applicability of kiss-gpon higher-dimensional problems has been addressed in Wilson et al. (2015) by means of low-dimensional projections. A more recent extension allows for a constant-time variance prediction using Lanczos methods [START_REF] Pleiss | Constant-time predictive distributions for Gaussian processes[END_REF]. The limitation of these approaches is that inducing inputs must abide by the Kronecker structure due to computational acceleration. This leads to the partial restriction on the freedom of the optimization of inducing inputs.

Outline and Contributions of Thesis

The content of this thesis is organized as follows:

• Chapter 2 starts with a brief introduction to Gaussian Processes (gps).

We also investigate state-of-the-art techniques for dealing with the notorious limitation of gps on time and storage complexity as well as the flexibility of kernel function. In this text, these approaches is grouped into 1.3. Outline and Contributions of Thesis three main categories of approximations, namely inducing point-based approximations, structure exploiting approximations, random featurebased approximations are discussed. This chapter is intended to equip the reader with the background knowledge required for apprehending the underlying concepts presented in this thesis, and clarify how our contributions fit within the landscape of existing research on Gaussian process inference;

• Chapter 3 covers the first primary contribution of this thesis. The study expresses a thorough investigation of the calibration properties of Bayesian Convolutional Neural Networks (cnns) . Along a similar vein, independently of the works on Bayesian cnns, there are other attempts to impose a probabilistic formalism to cnns by integrating cnns with gps. Previous work on combining cnns with gps has been developed under the assumption that the predictive probabilities of these models are well-calibrated. We show that, in fact, current combinations of cnns and gps are miscalibrated. We propose a novel combination that considerably outperforms previous approaches to this aspect, while achieving state-of-the-art performance on image classification tasks.

• As alluded earlier, inducing point-based idea are a well-known approach to mitigate the computational bottleneck of gps in the large-scale problems. However, this solution still suffers cubic time complexity to the number of inducing points. Wilson et al Wilson and Nickisch (2015) propose to employ the Kronecker structure on inducing inputs to accelerate the approximation of covariance matrices. The trick also accompanies with significant restrictions on inducing inputs. Besides, the approach only performs well on low-dimensional datasets (Wilson and Nickisch, 2015). In Chapter 4, we address one limitation of sparse GPs, which is due to the challenge in dealing with a large number of inducing variables without imposing a special structure on the inducing inputs. In particular, we introduce a novel hierarchical prior, which imposes sparsity on the set of inducing variables. The study enables the possibility to use sparse GPs using a large number of inducing points without incurring a prohibitive computational cost.

• Finally, in Chapter 5, we summarize the contributions presented in this thesis. We conclude the thesis by a discussion to an outlook on possible extensions and future work.

Gaussian Processes for Big Data

Increasing the scalability and representational power of models without compromising performance is a core problem in machine learning. As emphasized in the introduction to this dissertation, the scalability of Gaussian Processes to training size and dimensionality is significantly limited by algebraic operations, which discourages their application to datasets having more than a few thousands of examples or high-dimensional covariates. Additionally, the flexibility of Gaussian Processes is possibly weakened by the need to choose a kernel functions, which might lead to difficulties in learning the intricate patterns concealed in the data. This chapter is a literature review on the developments of gps in both aspects, which involve the major contributions of the thesis.

Overview

Gaussian Processes (henceforth gps) which are powerful non-parametric Bayesian models can yield sensible predictions with a small number of available observations. However, it is notorious that gps suffer from high complexity in terms of both computation and storage with respect to training size N , i.e. O (N 3 ) and O (N 2 ) respectively, so they not the primary choice in datasets with a massive number of data points. To broaden the application of gps to larger datasets, there is plenty of ideas in the literature that have been proposed and analyzed. According to the groupings mentioned in (Liu et al., 2018b), these approaches are categorized into global and local approximations. While the former approximate the full kernel matrix by a global distillation, the latter abide to the divide-and-conquer concept and make predictions using a local subset of training data. We further split global approximations into subcategories: Inducing Point-Based Approximation and Random Feature-Based Approximation.

Gaussian Processes

As alluded earlier, a modeling approach in the era of big data should possess not only a powerful fitting capability but also a firm mechanism on predictive uncertainty. Bayesian nonparametrics is obviously an ideal candidate as it offers flexibility as well as calibrated predictive uncertainties. The philosophy and motivation of this area have been well discussed by a number of authors [START_REF] Ghosh | Bayesian nonparametrics[END_REF]; [START_REF] Hjort | Bayesian Nonparametrics. Cambridge Series in Statistical and Probabilistic Mathematics[END_REF]; [START_REF] Ghahramani | Bayesian non-parametrics and the probabilistic approach to modelling[END_REF].

Gaussian Process (gps) [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] are an attractive way of doing non-parametric Bayesian modeling in supervised learning problems.

Firstly, I succinctly introduce Gaussian Processes Regression (gpr) which is the simplest way to describe gps.

Gaussian Processes for Regression

Given a dataset D of N examples, D = {(x n , y n ) | n = 1, . . . , N }, where
x n denotes the n-th input vector (covariates) and y n denotes the n-th scalar output or target; the column vector inputs for all N cases are aggregated in the D × N design matrix X, and the outputs are collected in the vector y, so we can write D = (X, y). We would like to specify a function y representing the correlation between inputs and its targets, i.e. y n = y (x n ).

From a generative perspective, the observable labels y (x n ) are modeled via an appropriate conditional likelihood p (y (x n ) | f (x n )), where f is the latent function which can also be perceived as the intermediate representation of function y. In regression, the conditional likelihood is intuitively often assumed to be a Gaussian with mean of f and variance of σ 2 n , i.e.

p (y | f , σ 2 n ) = N (y | f , σ 2 n I).
In general, the function f can be viewed as

a column vector f , i.e. f = [f (x 1 ) , . . . , f (x N )]
p (f | X, ζ, θ) = N (m X , K X ) . (2.1)
where m X are column N -dimensional vector containing mean values at N

covariates, i.e. m X = [m (x 1 | ζ) , . . . , m (x N | ζ)] T ; and K X is a N × N sy-
metric and positive semi-definite matrix representing the correlation between latent random variables each other,

[K X ] ij = k (x i , x j | θ).
Hyper-parameter optimization. For convenience sake, we introduce ψ the set of all parameters involving mean function's parameters ζ, kernel parameters' θ and variance of likelihood σ 2 n , i.e. ψ = (σ 2 n , ζ, θ). Given dataset D,

Gaussian Processes Regressors are fitted to D by optimizing hyper-parameter ψ using the logarithm marginal likelihood function, log p (y | X, ψ). In general, the marginal likelihood can be found by marginalizing over latent random variables f .

p (y | X, ψ) = p y | f , σ 2 n p (f | X, ζ, θ) df . (2.2)
Thanks to the Gaussian likelihood p (y | f , σ 2 n ), we can derive an analytic form for the marginal likelihood as the Gaussian likelihood and Gaussian prior are conjugate to each other.

p (y | X, ψ) = N y | m X , K X + σ 2 n I . (2.3)
Setting K σ 2 n as K X + σ 2 n I, the logarithm marginal likelihood is written as:

log [p (y | X, ψ)] = - 1 2 log K σ 2 n - 1 2 (y -m X ) T K -1 σ 2 n (y -m X ) - N 2 log 2π.
(2.4) The quadratic form appearing in this expression corresponds to the model fit term of the GPR, advocating parameter settings that fit the data well. In contrast, the log determinant term penalizes overly complex models that are characterized by kernel matrices which are diagonally dominant, indicating little interaction between observations. It follows that the optimal parameters ψ OP T are identified by maximizing this objective function using iterative gradient ascent.

Prediction. Generally, gps governs the distribution of a finite-dimensional vector including latent values at a set of covariates using a multivariate normal distribution. Therefore, the joint distribution of training latent values, f , and the testing latent values, f * , according to the gp prior is: Remind that, in regression, the likelihood of observable targets given training latent values are intuitively assumed to be a Gaussian with the variance of σ 2 n ,

f f * ∼ N m X m X * , K X K X,X * K X * ,X K X * (2.5) If X
p (y n | f n ) = N (y n | f n , σ 2 n ).
It means that the functions for observable targets can be modeled as a noisy version of latent function f a Gaussian noise with

variance of σ 2 n , y (x n ) = f (x n ) + ε, where ε follows N (ε | 0, σ 2 n ).
Assuming additive independent identically distributed Gaussian noise with variance σ 2 n , the prior on the noisy observations becomes:

y ∼ N y | m X , K X + σ 2 n I . (2.6)
We can write the joint distribution of the observed target values and the function values at the test locations under prior as:

y f * ∼ N m X m X * , K X + σ 2 n I K X,X * K X * ,X K X * (2.7)
To get the posterior distribution over function, we need to restrict this joint prior distribution to contain only those functions which agree with the observed data points. By virtue of the nice properties of the multivariate normal distribution, the operation of eliminating those violating the available observations is extremely simple, corresponding to conditioning the joint Gaussian prior distribution on the observations to give:

f * | X * , X, y, ψ ∼ N µ f * , Σ f * , where (2.8) µ f * = m X * + K X * ,X K X + σ 2 n I -1 (f -m X )
, and, (2.9)

Σ f * = K X * -K X * ,X K X + σ 2 n I -1 K X,X * .
(2.10)

Once again, thanks to Gaussian likelihood with noise variance of σ 2 n , the predictive distribution p (y * | X, y, ψ) turns out:

p (y * | X * , X, y, ψ) = p (y * | f * ) p (f * | X * , X, y, ψ) df * = N y * | µ f * , Σ f * + σ 2 n I (2.11)

Covariance function

In gps or any kernel machine learning methods, the notion of similarity between data points is crucial as the predictions are made based upon these similarities. Under the Gaussian process view, a covariance matrix specified by a kernel function defines nearness or similarity between latent random variables by using inputs. Therefore, it is able to encode our assumptions about the function which we wish to learn through. It is uncertain whether an arbitrary matrix of input pair x i and x j will be a valid kernel function or not.

The first purpose of the section is to show the properties and construction of a valid covariance function. In addition, examples of some commonly-used covariance functions in this dissertation are also given.

Construction and properties. The covariance matrix of the is constructed from a kernel function k of an input pair. Consider a gps for the sequence of N latent values, the dimensionality of the covariance matrix of gps is N × N , and the element at i-th row and j-column of the covariance matrix is kernel function values of x i and x j , k (x i , x j ). In general, the kinds of kernel function for all examples x i and x j in an input space X ⊂ R D :

k (x i , x j ) = φ (x i ) , φ (x j ) .
(2.12)

, where φ is a non-linear (or linear) map from the input space X to the feature space F, and ., . is an inner product. Due to being computed by the inner product, a kernel function must be symmetric and also satisfy the Cauchy-Schwartz inequality:

k (x i , x j ) = k (x j , x i ) , and k 2 (x i , x j ) ≤ k (x i , x i ) k (x j , x j ) .
(2.13) Practically, the kernel function k is usually specified directly, thus implicitly defining the map φ and the feature space F. Therefore, a kernel function is stated to be valid if it guarantees the existence of the feature space. Mercer [START_REF] Mercer | Functions of positive and negative type and their connection with the theory of integral equations[END_REF] showed that a necessary and sufficient condition for a symmetric function k (., .) to be a kernel is that it be positive definite. This means that for any set of x 1 , . . . , x N and any set of real numbers λ 1 , . . . , λ N , the function k must satisfy:

∀x 1 , . . . , x N ∈ X , ∀λ 1 , . . . , λ N ∈ R, N i,j=1 λ i λ j k (x i , x j ) ≥ 0. (2.14)
In summary, a symmetric positive definite function constructs a valid covariance matrix in kernel methods. As the positive definiteness possesses pleasant algebraic properties, a new kernel can be created from existing valid kernels.

Introducing a 1 and a 2 are positive real numbers, and k 1 and k 2 are valid kernels, a new kernel can be manipulated using a weighted summation or multiplication:

k (x i , x j ) = a 1 k 1 (x i , x j ) + a 2 k 2 (x i , x j ) . (2.15) k (x i , x j ) = k 1 (x i , x j ) k 2 (x i , x j ) .
(2.16)

k RBF (r) = exp - r 2 2l 2 .
(2.17)

, with positive parameter l defines the characteristic length-scale which indicating the complexity of underlying latent functions.

• The Matérn Covariance Function.

The Matérn class of covariance functions is given by

k Matérn (r) = 2 1-ν Γ (ν) √ 2νr l ν K ν √ 2νr l .
(2.18)

, with positive parameters ν and l, and K ν is a modified Bessel function [START_REF] Abramowitz | Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables[END_REF]. The most interesting cases of Matérn class for machine learning are ν = 3/2 and ν = 5/2, for which with D input dimensionality. It turns out that the term r/l in the isotropic kernel is replaced using a quadratic form. For example, the rbf kernel can be rewritten as:

k Matérn 3/2 (r) = 1 + √ 3r l exp √ 3r l , (2.19) k Matérn 5/2 (r) = 1 + √ 5r l + 5r 2 3l 2 exp - √ 5r l , ( 2
k rbf (x i , x j ) = exp - 1 2 (x i -x j ) T Λ -1 (x i -x j ) .
(2.21)

, where Λ = Diag [l 2 1 , . . . , l 2 D ] T . This interpretation of the lengthscales allows for automatic relevance determination whereby relevant features in the data are weighted by their corresponding lengthscale parameter. This can also be seen as an implicit form of feature selection [START_REF] Mackay | Bayesian interpolation[END_REF].

Non-Gaussian Likelihoods

Recall that in gp regression the Gaussian likelihood p (y | f ) is conjugate to the Gaussian prior p (f ). Thus, it is possible to calculate the marginal likelihood and carrying out inference in gp regression analytically. In contrast, these calculations are analytically intractable in gp models with a non-Gaussian likelihood. There is a plethora of approaches to deal with the problem, including the Laplace approximation method [START_REF] Williams | Bayesian classification with Gaussian processes[END_REF], expectation propagation [START_REF] Minka | Expectation Propagation for approximate Bayesian inference[END_REF], sparse approximation employing online learning schemes [START_REF] Lawrence | Fast Sparse Gaussian Process Methods: The Informative Vector Machine[END_REF][START_REF] Csató | Sparse on-line gaussian processes[END_REF] and methods attempting to characterize the full posterior [START_REF] Murray | Elliptical slice sampling[END_REF][START_REF] Filippone | A comparative evaluation of stochastic-based inference methods for gaussian process models[END_REF]Hensman et al., 2015b). As the prerequisite backgrounds for proposed models which will be introduced in the next chapters do not significantly depend on the techniques of approximating posterior with non-Gaussian likelihood, the discussion about the non-Gaussian likelihood or gps classification will not be provided in this manuscript.

Limitations of Gaussian Processes

Scalability. Theoretically, gps is an ideal approach for the supervised scenario in the era of big data. However, the scalability of gps is limited on small datasets including a few thousands of data points due to linear algebraic operations requiring large computational complexity. Having considered the optimization of gps hyper-parameters ψ, the problem of gps scalability is revealed. As alluded in section 2.2.1, the process of fitting gps regressors given a dataset can be done by using a gradient-based method with the target function of marginal likelihood p (y | X, ψ, y). Take gps regression with zero mean prior as an example, the gradients of marginal likelihood with respect to parameter ψ i is computed as:

∂ log [p (y | X, θ)] ψ i = - 1 2 Tr K -1 σ 2 n ∂K σ 2 n ∂ψ i + 1 2 y T K -1 σ 2 n ∂K σ 2 n ∂ψ i K -1 σ 2 n y. (2.22)
The computation of gradients involves with solving the linear system, i.e.

K -1

σ 2 n y where K σ 2 n is N × N covariance matrix with additive noise and y is Ndimensional column vector of outputs, where N is the number of data points.

Practically, this linear system is solved by using Cholesky decomposition to factorize the symmetric positive definite kernel matrix K σ 2 n into LL T , where

L is a lower triangular matrix. Generally, factorization with Cholesky decomposition necessitates O (N 3 ) operations. The calculation of the trace terms appearing in gradient formula also need O (N 3 ) operations. In the progress of computing the gradients, the lower triangular matrix L must be cached.

Therefore, the storage cost of the training phase is O (N 2 ).

Besides the cubic complexity in the training phase, the computational cost of gps inference also depends on the training size. On inspection of the predictive distribution given from equation 2.8 to 2.10, we can observe that evaluating this expression also involves the inversion of N × N kernel matrix. Theoretically, the computational cost for gp inference is also O (N 3 ). However, in practice, the inversion of K σ 2 n and the vector which is the multiplication of

K -1 σ 2
n and ym X can be recorded after the training phase. Therefore, the computational costs of predictive mean and predictive variance at an unseen data point are O (N ) and O (N 2 ). As discussed, the likelihood mapping the latent values to observation is not obligated to be a Gaussian, as the case of classification. Under these conditions, the computation of the marginal likelihood as well as inference is no longer analytic, and further approximations are required, and computational budgets required in these case is identical to gp regression. Due to the dependence of computational complexity on training size, gps are hindered to large-scale problems.

To strengthen the scalability, while retaining the desired prediction quality, a large number of scalable gps have been proposed. According to (Liu et al., 2018b), these scalable approaches are sorted into two main categories: local and global approximation. Local approximations arising from the divide-andconquer concept focus on the local area of input spaces (Gramacy and Lee, 2007;Yuksel et al., 2012b;[START_REF] Masoudnia | Mixture of experts: A literature survey[END_REF][START_REF] Rasmussen | Infinite mixtures of gaussian process experts[END_REF][START_REF] Associates | Variational inference for infinite mixtures of gaussian processes with applications to traffic flow prediction[END_REF][START_REF] Hinton | Training products of experts by minimizing contrastive divergence[END_REF][START_REF] Deisenroth | Distributed gaussian processes[END_REF][START_REF] Rullière | Nested kriging predictions for datasets with large number of observations[END_REF]Liu et al., 2018a). Whereas global approximations replace kernel matrix K X by a compact representation reducing the burden of computation. The substitution is done through global distillation which can be accomplished by several ways, e.g. use a small subset of training data [START_REF] Chalupka | A framework for evaluating approximation methods for Gaussian process regression[END_REF], or remove uncorrelated entries in K X using sparse kernel [START_REF] Gneiting | Compactly supported correlation functions[END_REF], or employ low-rank representation (Nyström approximation) [START_REF] Hensman | Gaussian processes for big data[END_REF][START_REF] Pratola | A unifying view of sparse approximate Gaussian process regression[END_REF][START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF]Wilson and Nickisch, 2015).

Representational power. Kernel functions hold a crucial role as it not only encodes our assumptions as well as the desired flexibility into the functions we wish to learn. Concerning the representational capability, kernel-based methods possibly lose their power as very limited kernels such as rbf kernel sharing a single length-scale across input are overused, e.g. in some gp-based approaches and, especially in Support Vector Machine (svm) . Having been encouraged by the achievement of deep architectures, there have been several

attempts to build kernel-based method that mimic deep neural networks, for example, multilayer arc-cosine kernel (Cho and Saul, 2009) which is built by successive kernel compositions, and kernel function at each layer are defined via an integral representation, or convolutional multilayer kernels [START_REF] Mairal | Convolutional kernel networks[END_REF] which are built by concatenations of convolutional layers, and the compact representation of the kernel are learned in a data-dependent manner.

Another approach to enhance the flexibility of kernel methods is to use its deep architecture, e.g. Deep gps [START_REF] Damianou | Deep Gaussian Processes[END_REF]Salimbeni and Deisenroth, 2017;[START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF].

Inducing Point Approximations

Prior approximation

Main idea. As mentioned in section 2.2.4, the computational bottleneck of Gaussian Processes (gps) stems from the algebraic operation of the full kernel matrix that appears in the prior distribution. Intuitively, the idea of employing the approximations to these true priors accelerating the computations come to a rescue for the problem of scalability. In this approach, the joint prior p (f * , f ) is modified in ways that reduces the computational cost.

Here, f * and f are the latent values at training points X and testing points X * respectively. For clarity, it is useful to derive the exact expression for the joint prior before discussing about the particular approaches employing the idea.

Without loss of generality, the mean of all priors is set to zero. Introducing the auxiliary random variables u, which are latent values at inducing inputs Z, the joint prior p (f * , f ) is expressed by marginalizing out u from the joint prior p (f * , f , u).

p (f * , f ) = p (f * , f , u) du = p (f * , f | u) p (u) du (2.23)
Due to the consistency of gps, all probabilistic components appearing in equation 2.23, i.e. the joint prior p (f * , f , u) and the conditional prior p (f * , f | u)

and the prior p (u) are Gaussian densities. Introducing f as the general latent values for both training and testing points, we can rewrite the joint prior as:

Joint prior: p (f * , f ) = p f = p f , u du = p f | u p (u) du (2.24) Prior: p (u) = N (u | 0, K Z,Z ) (2.25) Conditional: p f | u = N f | K X,Z K -1 Z,Z u, K X, X -Q X, X (2.26)
Here, X generally indicates training inputs X and testing inputs X * . Assuming that A and B are the matrices constructed by concatenating covariates likewise X and X * , we define K A,B as a cross covariance matrix whose element in the i, j position is the covariance between the i-th covariate in A and j-th covariate in B. We also introduce the shorthand no-

tation Q A,B = K A,Z K -1 Z,Z K Z,B
which can be seen as an approximation to K A,B using inducing inputs Z. For simplicity, we use the Gaussian likelihood

p (y | f ) = N (y | f , σ 2 n I). The predictive latent distributions p (f * | y) can be
written in a closed-form using Gaussian density:

p (f * | y) = N f * |K X * ,X K X,X + σ 2 n I -1 y, K X * ,X * -K X * ,X (K X,X ) -1 K X,X * + σ 2 n I .
(

p (f * , f ) ≈ q (f * , f ) = q (f * | u) q (f | u) p (u) du.
(2.28)

Following the unifying view mentioned by Quiñonero Candela and Rasmussen (2005), a particular algorithm complying with the idea of prior approximation corresponds to different additional assumptions about the two approximate inducing conditionals q (f | u) and q (f * | u) appearing in the approximation defined in 2.28. The method pic [START_REF] Snelson | Local and global sparse Gaussian process approximations[END_REF] mentioned at the end of the section is also an extension of the idea by using another way to approximate the joint prior p (f * , f ).

Subset of Data. The most straightforward approach to reduce the computational burden of gps, which stems from the inverse of the kernel matrix K X , is to work on subsets of the data (henceforth sod ), D sod for the whole training points, D, i.e. simply speaking, we use K Xsod instead of K X . By restricting the number of data point M in X sod to be less than the total number of observations, N , the computational cost will decrease from O (N 3 ) to O (M 3 ). In case X sod is specified in an appropriate manner, the approaches of sod will produce reasonable predictive distributions. Otherwise and most often, sod yields overconfident predictions. On the inspection of the selection of D sod , one could, for example, randomly choose M data points, use clustering techniques to divide the training data to M subsets and then choose the centroids as representative for all the whole data sets, or employ online learning scheme with criteria based on information theoretic principles, i.e.

differential entropy [START_REF] Lawrence | Fast Sparse Gaussian Process Methods: The Informative Vector Machine[END_REF], to choose active data points sequentially.

Turning to the unifying view mentioned above, u and f are replaced by f sod which are the latent values of subset input D sod . sod also uses the true testing conditional distribution instead of its approximation, i.e.

q (f * | u) = p (f * | f sod ).
The joint prior turns out to be:

p (f * , f ) → p (f * , f sod ) = p (f * | f sod ) p (f sod ) (2.29)
Subset of Regressors. According to the study on Subset of Regressors (sor)

of [START_REF] Silverman | Some aspects of the spline smoothing approach to nonparametric regression curve fitting[END_REF] and [START_REF] Wahba | The biasvariance tradeoff and the randomized gacv[END_REF], [START_REF] Smola | Sparse greedy gaussian process regression[END_REF] have adjusted sor for a sparse approximation to Gaussian Processes Regression.

sor assumes that there is a deterministic relationship between latent values, i.e. f * and f , and inducing variables u. This correlation can be represented as a Gaussian distribution with zero covariance as follows:

q sor f | u = N f | K X,Z K -1 Z,Z u, 0 .
(2.30)

Substituting q sor f | u to the Equation 2.28, the approximated joint prior is:

q sor (f , f * ) = N f f * 0 0 , Q X,X Q X,X * Q X * ,X Q X * ,X * (2.31)
From the approximated joint prior and the Gaussian likelihood, we can obtain the approximated predictive latent distribution:

q sor (f * | y) = N f * |Q X * ,X Q X,X + σ 2 n I -1 y, Q X * ,X * -Q X * ,X Q X,X + σ 2 n I -1 Q X,X * (2.32)
Having observed the true predictive latent distributions p (f * | y) defined in Equation 2.27, the approximated predictive latent distributions q sor (f * | y) are identical with p (f * | y), except that the covariance K has been substituted by Q. Therefore, sor approximation operates as an exact Gaussian

Processes with the covariance matrix K sor defined by the kernel function

k sor (x i , x j ) = K x i ,Z K -1 Z,Z K Z,x j .
Deterministic Training Conditional. According to the analysis of [START_REF] Williams | Observations on the nyström method for gaussian process prediction[END_REF], sor can yield negative predictive variances due to the approximation of the full covariance matrix using the Nyström method. In order to avoid these nonsensical predictive variances, Seeger and Williams (2003)

p (y | f ) = N y | f , σ 2 n I ≈ q dtc (y | u) = N y | K X,Z K -1 Z,Z u, σ 2 n I (2.33)
This approach uses the point estimate to variational distribution over training latent value similarly to sor, it remains to use the exact test conditional defined in 2.26.

q dtc (f | u) = q sor (f | u) = N f | K X,Z K -1 Z,Z u, 0 .
(2.34)

q dtc (f * | u) = p (f * | u) = N f * | K X * ,Z K -1 Z,Z u, K X * ,X * -Q X * ,X * . (2.35)
Another difference between sor and dtc is indicated in the joint prior. While sor uses Q X,X * to govern the relation between testing points, dtc use the exact full covariance matrix K X * ,X * .

q dtc (f , f * ) = N f f * 0 0 , Q X,X Q X,X * Q X * ,X K X * ,X * .
(2.36)

The predictive distribution of dtc is similar to sor, but Q X * ,X * is replaced by K X * ,X * : Integrating sgpp into the unifying framework, we can observe clearly the differences in the formalism between sgpp and sor and dtc. While the likelihood variance of dtc is characterized by only the noise variance, the likelihood variance of sgpp also takes into account the residual difference between Diag (K X,X ) and Diag (Q X,X ). sgpp assumes that the auxiliary variables u induces the relation of training latent variables f . Due to this assumption, sgpp can be called Fully Independent Training Conditional (fitc) approximation.

q dtc (f * | y) = N f * |Q X * ,X Q X,X + σ 2 n I -1 y, K X * ,X * -Q X * ,X Q X,X + σ 2 n I -1 Q X,
The approximation to the likelihood as well as the variational distribution of training and testing latent values given u also relies on the projection as in dtc, but the predictive variances are more sophisticated than dtc.

p (y | f ) ≈ q fitc (y | u) = N y | K X,Z K -1 Z,Z u, Diag [K X,X -Q X,X ] + σ 2 n I . (2.38) q fitc (f | u) = N n=1 p (f n | u) = N n=1 N f n | K xn,Z K -1 Z,Z u, K xn,xn -Q xn,xn .
(2.39)

q fitc (f * | u) = p (f * | u) = N f * | K X * ,Z K -1 Z,Z u, K X,X -Q X,X .
(2.40)

By introducing A as a square matrix, the operator Diag (A) constructs a diagonal matrix whose elements are taken from the diagonal line of A. The approximation to joint prior q fitc (f , f * ) is similar to dtc, except for the covariance matrix governing the relation of training latent variables. While dtc uses Q X,X in q dtc (f , f * ) defined in the equation 2.36, fitc also uses Q X,X in the approximation to joint prior, but remain the true kernel value at diagonal elements.

q fitc (f , f * ) = N f f * 0 0 , Q X,X + Λ Q X,X * Q X * ,X K X * ,X * (2.41) , where Λ = Diag [K X,X -Q X,X ].
From the joint prior defined in 2.41, the predictive distribution of fitc or sgpp turns out:

q fitc (f * | y) = N f * |Q X * ,X Q X,X + Λ + σ 2 n I -1 y, K X * ,X * -Q X * ,X Q X,X + Λ + σ 2 n I -1 Q X,X * . (2.42)
Observe the approximation to joint prior in fitc defined in the equation 2.41, we realize that the training and testing covariance are constructed heterogeneously. Therefore, the approximation fitc does not comply with the strict definition of gps where the covariance for all points must be computed by identical manners. In contrast, if the assumption of conditional independence given active points is extended to the testing case, fitc turns into another approach which is logically called Fully Independent Conditional (fic) . fic is equivalent to Gaussian Processes with the covariance function

k fic (x i , x j ) = k sor (x i , x j ) + δ i,j (k (x i , x j ) -k sor (x i , x j ))
, where δ i,j is Kronecker delta function. The prior and predictive distribution implied by fic is:

q fic (f , f * ) = N f f * 0 0 , Q X,X + Λ Q X,X * Q X * ,X Q X * ,X * + Λ * (2.43) q fic (f * | y) = N f * |Q X * ,X Q X,X + Λ + σ 2 n I -1 y, Q X * ,X * + Λ * -Q X * ,X Q X,X + Λ + σ 2 n I -1 Q X,X * .
(

q pitc (f | u) = N f | K X,Z K -1 Z,Z u, Λ (2.45) q pitc (f * | u) = p (f * | u) = N K X * ,Z K -1 Z,Z u, K X * ,X * -Q X * ,X * . (2.46) , where Λ = bkdiag [K X,X -Q X,X
] is a block diagonal matrix that is not clearly specified in Quiñonero Candela and [START_REF] Pratola | A unifying view of sparse approximate Gaussian process regression[END_REF]. An intuitive blocking structure is to group training points using clustering techniques as mentioned in [START_REF] Snelson | Local and global sparse Gaussian process approximations[END_REF]. Similar to fitc, the approximation to joint prior of pitc is defined as:

q pitc (f , f * ) = N f f * 0 0 , Q X,X + Λ Q X,X * Q X * ,X K X * ,X * (2.47)
The approximation to the predictive distribution of pitc is identical to fitc defined in equation 2.42, except for the substitution of a block diagonal matrix.

q pitc (f * | y) = N f * |Q X * ,X Q X,X + Λ + σ 2 n I -1 y, K X * ,X * -Q X * ,X Q X,X + Λ + σ 2 n I -1 Q X,X * .
(2.48)

As argued by [START_REF] Snelson | Local and global sparse Gaussian process approximations[END_REF], predictions obtained by pitc are empirically identical to fitc and fic given a specified set of active positions and hyper-parameters. They have speculated that mean predictions of pitc are still just a weighted sum of basis functions centered on the same inducing inputs as in fitc or fic, and the blocking structures on training covariance only changes the weights slightly. In addition, the structure of covariance of pitc defined in equation 2.47 means that the pitc approximation is not equivalent to a Gaussian Processes with a particular kernel function. To solve these problems, [START_REF] Snelson | Local and global sparse Gaussian process approximations[END_REF] relax the assumption of conditional independence between training and testing latent variables given inducing variable, i.e. do not approximate p (f ,

f * | u) by q (f | u) q (f * | u).
They treat the training and testing inputs in the same manner, and put them into S blocks using clustering techniques. Consider a single testing input x * which are assigned to block B S , then the joint prior are approximate as:

p (f , f * ) = p (f , f * | u) p (u) du ≈ p (f B S , f * | u) S-1 s=1 p (f Bs | u) p (u) du.
(2.49) According to the approximation, the assumption of using partial independence applies to both training and testing points. Therefore, this idea is logically called Partially Independence Conditional (pic) . Thanks to the relaxation of conditional independence, pic corresponds to a Gaussian Process with covariance matrix K pic .

K pic (x i , x j ) = Q (x i , x j ) + ψ (x i , x j ) [K (x i , x j ) -Q (x i , x j )] .
(2.50)

, where ψ (x i , x j ) = 1 if x i and x j are in the same block 0 otherwise.

(2.51)

The predictive distribution implied by pic is identical to the exact predictive distribution, except for the alternation of K by K pic .

Posterior Approximations

Weakness of Prior Approximations.

As alluded to previously, the aforementioned algorithms complying with the idea of prior approximation operate as an exact Gaussian Processes with a particular kernel function or an approximation to covariance matrices, i.e. sor [START_REF] Smola | Sparse greedy gaussian process regression[END_REF] and fic [START_REF] Snelson | Variable noise and dimensionality reduction for sparse gaussian processes[END_REF] and pic [START_REF] Snelson | Local and global sparse Gaussian process approximations[END_REF]. Suppose we would like to employ M inducing variables which are latent values at some auxiliary inputs Z to approximate the gp prior.

The quality of these sparse approximations depends on the optimization of Z and hyper-parameters, i.e. kernel's parameters and variance noise (for Gaussian likelihood). An approximation to the true marginal likelihood defined in 2.1 allows us to select Z and other hyper-parameter using a gradient-based iterative method. For example, consider a zero-mean gp, Projected Process approximation (pp) (Seeger et al., 2003) and Sparse Gaussian Processes using Pseudo-points (sgpp) [START_REF] Snelson | Variable noise and dimensionality reduction for sparse gaussian processes[END_REF] following the idea of prior approximation replace the logarithm of gp marginal likelihood

F gp = log [p (y | X)
] by F pp and F sgpp :

gp: F gp = log N y | 0, K X,X + σ 2 n I .
(2.52)

pp : F pp = log N y | 0, Q X,X + σ 2 n I .
(2.53)

sgpp : F sgpp = log N y | 0, Q X,X + Λ + σ 2 n I . (2.54)
where we recall that

Q X,X = K X,Z K -1 Z,Z K Z,X is the Nyström approximation to K X,X using inducing inputs Z. Λ = Diag [K X,X -Q X,X ]
is the difference on diagonal elements between the true kernel values and approximated ones.

Observe F pp and F sgpp , the covariance of approximate marginal likelihood are parameterized by inducing inputs Z. While it is tempting to think that the introduction of Z in kernel function improves the representational power of approximate gps, the highly-parameterized form probably lead to an overfitting problem because the continuous optimization of F pp and F sgpp with respect to Z does not reliably approximate true gp.

Main idea of Posterior Approximations.

In order to deal with the lack of consideration of the convergence between true gp and approximate ones, we intuitively find Z by minimizing the distance of the approximated predictive distributions produced by the inducing points and the true ones.

Further, the idea also allows us to access the divergence between the true gp and sparse approximation to gp. Since both of predictive distribution p (f * | y)

and posterior p (f | y) are conditional prior given observations, the selection of Z based upon the idea is equivalent to minimize the Kullback-Leibler divergence between the approximate posterior q (f ) and the true posterior p (f | y).

Starting from the true gp conditional prior over arbitrary auxiliary variable v given observations y, we construct the approximate gp posterior using M inducing points. We can express the conditional prior p (v | y) by integrating out inducing variables u and training latent values f as follows:

p (v | y) = p (v | u, f ) p (f | u, y) p (u | y) df du.
(2.55)

By assuming u completely capture the relation between v and f , it holds that

p (v | u, y) = p (v | u).
Thanks to the assumption of conditional independence, the variable f only appears in p (f | u, y), and therefore, f is canceled out as p (f | u, y) df = 1. Subsequently, the above p (v | y) can be written as q (v):

q (v) = p (v | u) q (u) du = q (v, u) du.
(2.56)

where

q (v) = p (v | y) and q (u) = p (u | y).
Practically speaking, it is infeasible to find inducing variables u which are sufficient statistics for the parameters f . Thus, q (v) should be expected as an approximation to p (v | y).

Subsequently, the q (u) can be represented by a parameterized form.

Since the joint variable [z, y] T and [z, u] T follow a gp, the conditional priors of p (z | y) and p (z | u) are also Gaussian densities. Intuitively, q (z)

which is the approximation to p (z | y) should be also a Gaussian. Thanks to the equation 2.56 and the conjugacy properties, we see that assuming q (u) a variational Gaussian distribution defined by a mean vector m and covariance matrix S turns q (v) to be a Gaussian. Introducing X as the indices of v, we can express q (v) under a closed form:

q (v) = N Am, K X, X + A (S -K Z,Z ) A .
(2.57)

where A = K X,Z K -1 Z,Z . Since v is an arbitrary variable representing all latent function values, q (v) can be perceived as an approximation to gp posterior q (f ) or predictive distributions q (f * ).

Turning to the idea of posterior approximation, all parameters θ, e.g. inducing inputs or hyper-parameters, are selected to minimize the Kullback-Leibler divergence between the approximate posterior q (f ) and the true posterior p (f | y). Equivalently, we can minimize a distance between the aug-

mented variational posterior q (f , u) defined in equation 2.56, i.e. q (f , u) = p (f | u) φ (u) and the augmented true posterior p (f , u | y). θ * = arg min θ KL [q (f , u) || p (f , u | y)] (2.58)
Taking further analysis, we see that KL [q (f , u) || p (f , u | y)] can be represented as:

KL [q (f , u) || p (f , u | y)] = log [p (y)] -E q(f ,u) log p (f , u, y) q (f , u) . (2.59)
Since log [p (y)] is constant for q (f , u), learning all parameters can be inferred by maximizing F q defined as follows:

F q = E q(f ,u) log p (f , u, y) q (f , u) = p (f | u) q (u) log p (y | f ) p (u) q (u) df du.
(2.60) Sparse Variational Gaussian Processes Regression. The most wellknown representative following the idea of posterior approximation is proposed by [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF] using variational inference technique [START_REF] Blei | Variational inference: A review for statisticians[END_REF].

In the approach, the lower bound to marginal likelihood is developed from F q defined in 2.60. To derive a tighter bound, they firstly maximize the bound F q by analytically solving for the optimal choice of the variational distribution q * (u). By differentiating 2.60 with respect to q (u) and using

Gaussian likelihood p (y | f ) = N (y | f , σ 2 n I
), the optimal q * (u) is derived as follows:

q * (u) = N (u | m * , S * )
, where, (2.61)

m * = σ 2 n K Z,Z K Z,Z + σ 2 n K Z,X K X,Z -1 K Z,X y (2.62) S * = K Z,Z K Z,Z + σ 2 n K Z,X K X,Z -1 K Z,Z (2.63)
By replacing q * (u) into the bound F q defined in Equation 2.60, we obtain the lower bound of Sparse Gaussian Processes for Regression (sgpr) proposed by [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF]:

F sgpr = log N y | 0, σ 2 n I + Q X,X - 1 2σ 2 n Tr (Λ) . (2.64)
where we recall that Q

X,X = K X,Z K -1 Z,Z K Z,X is the Nyström approximation to K X,X using inducing inputs Z, and Λ = Diag [K X,X -Q X,X ]
is the difference on diagonal elements between the true kernel values and approximated ones. Observe the approximation to gp marginal likelihood of the approach of Projected Process Approximation (pp) or Deterministic Training Conditions (dtc) defined in 2.53, we can rewrite F sgpr in terms of F pp :

F sgpr = F pp - 1 2σ 2 n T r (Diag (K X,X -Q X,X )) . (2.65)
It is obvious that sgpr differs dtc only by trace term, which have a significant impact on the inference. Intuitively, the T r (Diag (K X,X -Q X,X )) represents the total variance of predicting the latent variables f given u. When maximizing the bound F sgpr , the positive trace term should be decreased, and, in particular, the fact of the trace is zero means that u recover the full gp.

Therefore, the trace term not only seeks to deliver a good inducing set but also prevents sgpr from overfitting.

Stochastic Variational Inference for Gaussian Processes. A downside of sgpr proposed by [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF] is that the computational and storage cost depends on the training size N linearly. On the inspection of the bound Processes. This approach is, therefore, abbreviated by svi . While Titsias' bound are derived by replacing q (u) by optimal distribution for inducing variable q * (u) defined in equation 2.61, svi [START_REF] Hensman | Gaussian processes for big data[END_REF] parameterize the variational distribution q (u) by a Gaussian density N (m, S). Substituting N (m, S) for q (u) in the general bound F q defined in the equation 2.60, the bound F svi are obtained as follows:

F SGP R defined
F svi = N i=1 log N y n | K x i ,Z K -1 Z,Z m, σ 2 n - 1 2σ 2 n K x i ,Z K -1 Z,Z SK -1 Z,Z K Z,x i - 1 2σ 2 n (K x i ,x i -Q x i ,x i ) -KL (q (u) || p (u)) .
(2.66)

Due to the Gaussian form of q (u), the KL term can be expressed analytically with the computational cost of O (M 3 ). The most important property of F svi is that it can be written as a sum of N terms, each of them corresponds to one pair of input and output (x i , y i ). This allows us to perform stochastic gradient ascent by using a mini-batch I as follows:

F svi ≈ N |I| (x i ,y i )∈I log N y i | K x i ,Z K -1 Z,Z m, σ 2 n - 1 2σ 2 n K x i ,Z K -1 Z,Z SK -1 Z,Z K Z,x i - 1 2σ 2 n (K x i ,x i -Q x i ,x i ) -KL (q (u) || p (u)) .
(2.67)

Apart from accelerating the computation cost by applying stochastic variational inference, the factorization over training examples allows the combina-tion of svi and non-Gaussian likelihood. As a consequence, a more general approach of svi has also proposed by Hensman et al. (2015a), which is called Scalable Variational Gaussian Processes (svgp) . The bound F svgp can be derived easily by rewriting F q in equation 2.60:

F svgp = F q = p (f | u) q (u) log p (y | f ) p (u) q (u) df du = E q(f ) log p (y | f ) -KL (q (u) || p (u)) = N i=1 E q(f i ) log p (y i | f i ) -KL (q (u) || p (u)) ≈ N |I| (x i ,y i )∈I E q(f i ) log p (y i | f i ) -KL (q (u) || p (u)) .
(2.68)

, where q (f i ) is calculated as q (v) defined in equation 2.57. In case the likelihood p (y i | f i ) is Gaussian, the variational expectation term can be expressed analytically. In general, the one-dimensional integrals of the log-likelihood can be computed by Gauss-Hermite quadrature as in Hensman et al. (2015a).

Further Improvement. The approaches of posterior approximation [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF][START_REF] Hensman | Gaussian processes for big data[END_REF]Hensman et al., , 2015a) can be further improved in various ways. The first direction is to apply a Bayesian treatment to all kernel hyperparameters rather than optimizing them, which is prone to overfitting (Titsias and Lazaro-Gredilla, 2013;Hensman et al., 2015b;[START_REF] Yu | Stochastic variational inference for bayesian sparse gaussian process regression[END_REF]. Another extension is to allow to work with a non-Gaussian posterior, e.g. mixture of

Gaussians (Nguyen and Bonilla, 2014a), or a free-form posterior (Hensman et al., 2015b).

Structure Exploiting Approximations

Main idea. Consider a gp with Gaussian likelihood p

(y i | f i ) = N (y i | f i , σ 2 n ),
the gradients of logarithm of the marginal likelihood p (y | X) with respect to an arbitrary trainable parameter ψ is as follows:

∂ log [p (y | X, θ)] ψ = - 1 2 Tr K -1 σ 2 n ∂K σ 2 n ∂ψ + 1 2 y T K -1 σ 2 n ∂K σ 2 n ∂ψ K -1 σ 2 n y. (2.69)
where

K σ 2 n = K X,X + σ 2 n I.
Traditionally, the Cholesky decomposition is ap- [START_REF] Golub | Matrix computations[END_REF].

plied to factorize K σ 2 n = LL T which cost O (N 3 ) (Golub
Therefore, the computational problems of these gradients start to arise when N exceed a few thousands. It is possible to enhance the scalability of the computations by imposing a special algebraic structure on the kernel matrix 2. Gaussian Processes for Big Data K X,X . A well-known approach following the idea of structural exploitation is to use Kronecker product with the assumption of grid-structure dataset and tensor product kernel [START_REF] Saatçi | Scalable Inference for Structured Gaussian Process Models[END_REF][START_REF] Gilboa | Scaling Multidimensional Inference for Structured Gaussian Processes[END_REF].

Grid-structured data. According to the exposition of Chapter 5 in Saatçi's dissertation, we assume all input points X are located on a Cartesian grid, i.e.

X = X 1 × • • • × X D (2.70)
, where X d represents all distinct input locations along dimension d, and operator × indicates the Cartesian product between vectors. The number of elements of the vector X d is generally arbitrary, i.e. we can say that

X d ∈ R G d
where G d is the size of vector X d . The definition of Cartesian product entails that X is restricted to contain exactly D d=1 G d points which are put on the D-dimensional Cartesian grid. Though a grid-structured data can enable the computational acceleration, the number of data points grows exponentially with dimensions, and, consequently, the limitation of the computational resource is quickly reached. Therefore, the speed-up of gp using the idea of grid-structured data is feasible with few dimensions. For example, following [START_REF] Saatçi | Scalable Inference for Structured Gaussian Process Models[END_REF], the applicability of gp on multidimensional grid data is restricted with the dimension which is less than 8. Despite severely suffering from the curse of dimensionality, this structure arises naturally in several spatial-temporal problems such as climate modeling, where the input points generally denote latitude and longitude coordinates that can be further augmented with some periodically spaced time dimension. Multimedia such as images and videos are also likely to inherently have such structure.

Tensor product kernel. In this section, the covariance functions are assumed to be tensor product kernels, which compute the covariance as a separable product over dimensions. Introducing two D-dimensional covariates x i and x j belonging to the grid-structure input space X mentioned above, the covariance between x i and x j can be written as:

k (x i , x j ) = D d=1 k d (x i,d , x j,d ) .
(2.71) where x i,d ∈ X d is the d-th element of input x i and k d (., .) is any symmetric positive definite function which is described in section 2.2.2.

Algebraic advantages of the Kronecker method. Introducing A as m × n matrix and B as p × q matrix, the Kronecker product of A and B, denoted by A ⊗ B, is an mp × nq matrix having the following form:

A ⊗ B =    a 1,1 B . . . a 1,n B . . . . . . . . . a m,1 B . . . a m,n B    (2.72)
, more explicitly:

A⊗B =                   
                  
(2.73) For the sake of clarity, we mention the basic properties of Kronecker product with square matrices, which is helpful for a forthcoming explanation.

Bilinearity: A ⊗ (B + C) = A ⊗ B + A ⊗ C (2.74) Associativity: (A ⊗ B) ⊗ C = A ⊗ (B ⊗ C) (2.75)
Mixed-product property:

(A ⊗ B) (C ⊗ D) = AC ⊗ BD (2.76) Inverse: (A ⊗ B) -1 = A -1 ⊗ B -1 (2.77) Transpose: (A ⊗ B) T = A T ⊗ B T (2.78) Trace: Tr (A ⊗ B) = Tr (A) Tr (B) (2.79) Determinant: det (A ⊗ B) = (det A) D A (det B) D B (2.80) Vec: V ec CXB T = (B ⊗ C) V ec (X) . (2.81)
, where D A and D B are dimensions of matrices A and B. Introducing X as m-by-n matrix, V ec (X) denotes flatten operator yielding mn-dimensional vector.

Thanks to the assumptions of grid-structured data and tensor product kernel, the full covariance matrix for points on the grid can be evaluated by Kronecker product of a sequence of kernels:

K X,X = K 1 (X 1 , X 1 ) ⊗ K 2 (X 2 , X 2 ) ⊗ • • • ⊗ K D (X D , X D ) (2.82)
, where K d is G d × G d covariance matrix of the vector of scalar input locations X d . In order to see how Kronecker product gain the benefit in gp computations, we remind the logarithm of the marginal likelihood of a zero-mean gp with the Gaussian likelihood p These reductions come from the fast computation of eigendecomposition K X,X = QVQ T , where V is the diagonal matrix constructed by corresponding eigenvalues v i , and Q is the square matrix whose i-th column is the eigenvector q i of K X,X , and therefore, Q is guaranteed to be an orthogonal matrix, and consequently, Q -1 = Q T . Since K X,X can be expressed by Kronecker product, the computation for matrices Q and V is accelerated by separately carrying out the eigendecomposition of

(y i | f i ) = N (y i | f i , σ 2 n ): log p (y | X) = - 1 2 log K X,X + σ 2 n I - 1 2 y T K X,X + σ 2 n I -1 y - N 2 log 2π.
K 1 (X 1 , X 1 ) , . . . , K D (X D , X D ). Denoting Q d as matrix containing eigenvectors of K d (X d , X d ) and V d as a diagonal ma- trix of eigenvalues of K d (X d , X d ), i.e. K d (X d , X d ) = Q d V d Q T d , matrix Q
and V can be expressed as Kronecker products by using the Mixed-product property defined at 2.76. Actually, V d and V are diagonal matrices, and V are constructed by concatenating diagonal elements of V d .

Q = Q 1 ⊗ • • • ⊗ Q D , and V = Diag Diag (V 1 ) T , . . . , Diag (V D ) T (2.84)
Due to the V ec property mentioned at 2.81, the fast matrix vector multiplication are enabled by using the Algorithm kron-mvn mentioned in Saatçi's dissertation [START_REF] Saatçi | Scalable Inference for Structured Gaussian Process Models[END_REF].

In order to analyze the complexity of kron_mvm conveniently, I assume all

= G D Output: α = ⊗ D d=1 A d b. 1: α ← b. 2: for d ← D to 1 do 3: X ← reshape (α, G, N/G); 4: Z ← A d X; 5: α ← vec Z T 6: end for is repeated D times, therefore, the ultimate budget is O DN 1+ 1 D .
Turning to the matrix-vector multiplication (K X,X + σ 2 n I) -1 y appearing in the logarithm of marginal likelihood, it can be rewritten in terms of Q and V. Thanks to the property of T ranspose defined at 2.78, (K X,X + σ 2 n I) -1 y can be further represented by Kronecker product:

K X,X + σ 2 n I -1 y = Q V + σ 2 n I -1 Q T y = (Q 1 ⊗ • • • ⊗ Q D ) V + σ 2 n I -1 Q T 1 ⊗ • • • ⊗ Q T D y (2.85)
With the above expression, the matrix vector multiplication can be solved efficiently using the following steps:

α ← kron_mvm Q T 1 , . . . , Q T D , y α ← V + σ 2 n I -1 α α ← kron_mvm ([Q 1 , . . . , Q D ] , α) (2.86)
where kron_mvm is a procedure detailed in Algorithm 1. Remind that V is diagonal matrix containing the eigenvalues of block covariances

{K d (X d , X d )} D d=1 ,
the matrix V + σ 2 n I is also diagonal, and, therefore, its inversion can be computed with linear complexity.

The fast eigendecomposition of K X,X also speeds up the computation of log-

arithm of determinant of K X,X + σ 2 n I. Denoting v 1 , . . . , v N the diagonal ele- ments of V, we know that {v i } N i=1 are eigenvalues of K X,X as K X,X = QVQ T .
Due to the definition of eigenvalue, there is a relation between matrix K X,X , its eigenvalue v i and its corresponding eigenvector q i : K X,X q i = v i q i , then

(K X,X + σ 2 n I) q i = v i q i + σ 2 n Iq i = (v i + σ 2 n ) q i . Therefore, it can be de- rived that if v i is an eigenvalue of K X,X then v i + σ 2 n is also an eigenvalue of K X,X + σ 2 n I. In consequence, the logarithm of determinant of K X,X + σ 2 n I is reduced from O (N 3 ) to O DN 3 D which are the cost for eigendecomposition of D matrices {K d (X d , X d )} D d=1 . log K X,X + σ 2 n I = i log v i + σ 2 n .
(2.87)

Structural Kernel Interpolation. Despite their impressive computational acceleration of Kronecker-based methodology presented above, the main limitation of this approach is the restriction of grid-structured data. However, most datasets will not satisfy this requirement, making the application of such techniques narrow. In order to relax the condition of having observations at all possible input locations in the grid, there are several attempts such that missing observations and incomplete grid are also permitted [START_REF] Flaxman | Fast kronecker inference in gaussian processes with non-gaussian likelihoods[END_REF][START_REF] Wilson | Fast kernel learning for multidimensional pattern extrapolation[END_REF]. Ultimately, Wilson and Nickisch (2015) have extended the concept of the Kronecker method to a general scenario with the proposal of Kernel Interpolation for Scalable Structured Gaussian Processes (kiss-gp) . This method constrains that the set of inducing positions Z constructs a complete multidimensional grid. Consider D-dimensional problems and introduce Z d as a vector containing distinct inducing locations along with dimension d, we again define Z as Cartesian product of Z 1 , . . . , Z D , i.e.

Z = Z 1 ⊗ • • • ⊗ Z D .
Similarly, the size of Z is M = D d=1 G d where G d is the number of elements in the vector Z d . By utilizing the tensor product kernel defined above, the Kronecker idea enables the fast algebraic operations of the covariance matrix of inducing points. Nevertheless, setting a massive M could be problematic due to the time-consuming operations associated with the cross-covariance K X,Z between design matrix X and inducing locations Z. For example, the full covariance matrix K X,X which can be expressed by Nyström approximation K X,Z K -1 Z,Z K Z,X dominate the computations with quadratic complexity to M , i.e. O (N M 2 ). Instead of computing directly K X,Z , it is estimated by interpolating on the M × M covariance matrix K Z,Z . For example, if we would like to estimate k (x, z j ), for point x and inducing input z j , we can start by finding the two inducing points z a and z b which are the two closest to x.

Then, we can estimate k (x, z j ) by k (x,

z j ) = wk (z j , z a ) + (1 -w)k (z j , z b ),
where w and 1 -w are represented the relative distance from x to z a and z b .

Generally, the cross covariance K X,Z between design matrix X and inducing points Z can be interpolated by:

K X,Z ≈ KX,Z = WK Z,Z (2.88)
While M is expected very large in the scenario, W is constraint to be extremely sparse. There are several options to construct matrix W based upon various strategy including (i) local linear interpolation where each row of W contains only 2 non-zero entries or (ii) local cubic interpolation for greater accuracy with 4 non-zero elements per row.

As a consequence, from the Nyström approximation to full covariance K X,X , a further estimation can be obtained by substituting KX,Z for K X,Z . This general approach of approximation is so called Structured Kernel Interpolation (ski) .

K X,X ≈ K X,Z K -1 Z,Z K Z,X ≈ KX,Z K -1 Z,Z KZ,X = WK Z,Z W K ski (2.89)
By exploiting the fast Kronecker matrix-vector multiplications mentioned above, the overall complexity of learning

gp is O DM 1+ 1 D computations and O N + DM 2 D
storage. Nonetheless, this approach also introduces additional design choices, such as determining the optimal density of the interpolation point grid, which require further fine-tuning than the relatively more straightforward inducing point methods. In general, the grid density is expected to be heavily dependent on the choice of the kernel since more expressive kernels are likely to require a greater number of interpolation points and less sparse W. In summary, the combination of ski and Kronecker algebraic structure results in the method kiss-gp.

Random Feature Approximations

As highlighted earlier, the inducing point-based approximation is a well-known approach for improving gps' scalability. In these methods, a small number of pairs of inducing inputs and outputs are learned to define a new gps, which is expected to be close as possible to the gps, and the computational and storage cost now depend on the number of inducing points. These approaches are appropriate for locally complex functions. Intuitively, most inducing inputs would be located in regions where the function is complex, while the rest would be placed in regions where the function is simpler. Highly complex functions cannot be modeled well with these inducing point-based approaches.

In order to capture complex behaviors at a global level and improve the scalability of gps, random feature-based approximations were proposed by [START_REF] Lázaro-Gredilla | Sparse Spectrum Gaussian Process Regression[END_REF] and [START_REF] Gal | Improving the Gaussian Process Sparse Spectrum Approximation by Representing Uncertainty in Frequency Inputs[END_REF], which relies on spectral representations of kernel functions. For this kind of approximation, we only consider stationary gps whose covariance functions are written as a function of the distance between observations, i.e. k (x, x ) = k (xx ) = k (r).

The spectral density for non-stationary kernels can be found in Remes et al. (2017). The concept of spectral expressiveness and random feature expansions are discussed here because these concepts are essential in the next chapter 2

where we propose a novel combination of cnns and gps approximated with random features.

Spectral Representations

Kernel trick and its problem. Kernel methods are a class of algorithms enabling the operations in an infinite-dimensional feature space, which leads to an enhancement of representational power. This is materialized by observing that inference for these methods is expressed through inner products between test points and input points, e.g. svm [START_REF] Cortes | Support Vector Networks[END_REF].

Thanks to this observations and Mercer's theorem, we can implicitly define the transformation from the original space to the infinite-dimensional space by specifying the kernel function between points. This is the so-called kernel trick. However, the weakness of these methods is that algorithms needs to evaluate the kernel function between all pairs of datapoints. Consequently, large training sets incur large computational and storage costs.

Dual representation of a stationary kernel. Rahimi and Recht (2008) proposed an idea to define a transformation of the input space enabling a numerical approximation to kernel values without suffering a prohibitive cost.

Due to the significant impact on research communities working on kernelbased models such as support vector machines, kernel ridge regression, and ultimately gps, this seminal work is considered to be one of the most influential papers published in the previous decade. Their work is inspired by Bochner's theorem [START_REF] Rudin | Fourier analysis on groups[END_REF] which states that any continuous shift-

invariant normalized covariance function k (x i , x j ) = k (x i -x j
) is said to be positive definite if and only if it can be rewritten as the Fourier transform of some non-negative measure p (ω). The spectral density s (ω) can be constructed from k (r) and vice versa through Wiener-Khintchin theorem: Kernel Name Kernel function k (r) Spectral density p (ω)

k (r) = F -1 {p (ω)} = +∞ -∞ p (ω) exp iω T r dω (2.90) p (ω) = F {k (r)} = 1 2π +∞ -∞ k (r) exp -iω T r
Gaussian exp - ||r|| 2 2 2 (2π) -D 2 exp -||ω|| 2 2 Matérn 1/2 σ 2 exp - ||r|| 1 l 2 σ 2 l 1 l 2 + ||ω|| 2 2 -1 Laplacian exp (-||r|| 1 ) D d 1 π(1+ω 2 d )
Approximation of rbf Kernel using Random Fourier Features. Generally, we consider the rbf kernel parameterized by θ = (σ 2 , l 1 , ..., l D ) and its corresponding spectral density can be found using Equation 2.91 as follows:

k rbf (x i , x j | θ) = k rbf (x i -x j | θ) = k rbf (r | θ) = σ 2 exp D d=1 r 2 d l d .
(2.92)

p rbf (ω) = N ω | 0, Λ -1
, where Λ = Diag (l 1 , . . . , l D ) .

(2.93)

From equation 2.90, the kernel function can be rewritten as the expectation under the density p rbf (ω).

k rbf (x i , x j | θ) = k (r | θ) = σ 2 E p(ω) exp iω T r = σ 2 E p(ω) cos ω T r + i sin ω T r .
(2.94)

As sin (.) is an odd function, i.e. sin (-x) = -sin (x), the imaginary term can be canceled out from the expectation in Equation 2.94. Further, the kernel function can be approximated using N RF spectral samples ω from density function p (ω).

k rbf (x i , x j | θ) = k rbf (r | θ) = σ 2 E prbf(ω) cos r T ω ≈ σ 2 N RF N RF r=1 cos r T ω(r)
(2.95) Replacing r by x i -x j into equation 2.95, we can express the approximation of kernel function by an inner product representation:

k rbf (x i , x j | θ) ≈ σ 2 N RF N RF r=1 cos x T i ω(r) -x T j ω(r) = φ rbf (x i ) T φ rbf (x j ) , (2.96)
where φ rbf (x) is known as random features of x for rbf kernel, which is defined as follows:

φ rbf (x) = σ 2 N RF cos x T ω(1) , . . . , cos x T ω(N RF ) , sin x T ω(1) , . . . , sin x T ω(N RF ) T (2.97)
Approximation of order-one arc-cosine Kernel using Random Features. In addition to working with rbf, we also consider order-one arccosine covariance which is a prevalent kernel function.

k arc (x i , x j | θ) = σ 2 π Λ -1 2 x i Λ -1 2 x j [sin (α) + (π -α) cos (α)] , (2.98)
where θ = (σ, Λ = Diag (l 2 1 , . . . , l 2 D )) and α is the angle between Λ -1 2 x i and Λ -1 2 x j . Let H (.) be the Heaviside function. Following Cho and Saul (2009), this covariance can be written under an integral form:

k arc (x i , x j | θ) = 2σ 2 H ω T x i ω T x i H ω T x j ω T x j ×N (ω | 0, I) dω.
(2.99) The convenient integral representation allows for a Monte Carlo approximation obtaining a low-rank approximation to the covariance matrix involving Rectified Linear Unit (relu) activation (Cho and Saul, 2009).

φ arc (x) = 2σ 2 N RF max 0, x T ω(1) , . . . , max 0, x T ω(N RF ) T (2.100)
2.4.2 Random featured-based Gaussian Processes.

In this section, I firstly present a well-known study of approximation of gps using random features, which is proposed by [START_REF] Lázaro-Gredilla | Sparse Spectrum Gaussian Process Regression[END_REF].

The key novel idea is to sparsify the spectral representation of gps.

Sparse Spectrum Gaussian Process Regression. As alluded earlier, Gaussian Processes Regression (gpr) is introduced in function-space view.

Here, we remind that, by considering the dataset D = {X, y} and Gaussian

likelihood p (y | f ) = N (y | f , σ 2 n I), the predictive distributions p (y * | x * , D)
and the logarithm of the marginal likelihood log (y | θ) given parameters θ are expressed as follows:

p (y * | x * , D) = N y * | µ * , σ 2 * , where µ * = K x * ,X K X,X + σ 2 n I -1 y σ 2 * = σ 2 n + K x * -K x * ,X K X,X + σ 2 n I -1 K X,x * (2.101) log p (y | θ) = - N 2 log (2π) - 1 2 K X,X + σ 2 n I - 1 2 y T K X,X + σ 2 n I -1 y
(2.102) Computing the gradients of logarithm of the marginal likelihood with respect to related parameters requires the cubic cost to training size, i.e. O (N 3 ), which is unacceptable for large-scale data sets. In order to avoid the prohibitive cost, Lázaro-Gredilla et al. ( 2010) have employed the approximation of the covariance matrix using spectral representation. Consider for example ard kernel (a stationary anisotropic squared exponential covariance function):

k ard (x i , x j | θ) = k ard (r = x i -x j | θ) = σ 2 0 exp - 1 2 r T Λ -1 r , (2.103)
where Λ = Diag ([l 2 1 , . . . , l 2 D ]). Based on the dual representation of the stationary kernel mentioned above, we can approximate the k ard (x i , x j ) using N RF spectral samples, and express the approximation as an inner product:

k ard (x i , x j ) ≈ σ 2 0 N RF N RF r=1 cos r T ω(r) = σ 2 0 N RF φ (x i ) T φ (x j ) , (2.104)
where ω(r From the transformation φ (.), we construct 2N RF by N matrix of random features Φ X = [φ (x 1 ) , . . . , φ (x N )]. Now, the full kernel matrix K X can be approximated as follows:

) ∼ p (ω) = N ω | 0, Λ -1 ,
K X ≈ σ 2 0 N RF Φ T X Φ X (2.106)
Replacing the kernel matrix by this approximation in equation 2.101 and 2.102, we obtain the spectral approximation of predictive distribution with mean µ * and variance σ 2 * :

µ * ≈ φ (x * ) A -1 Φ X y, and σ 2 * ≈ σ 2 n + σ 2 n φ (x * ) T A -1 φ (x * ) , (2.107) where A = Φ X Φ T X + N RF σ 2 n σ 2 0 I.
Similarly, we also obtain the approximate logarithm of the marginal likelihood:

log p (y | θ) ≈ -y T y -y T Φ T X A -1 Φ X y / 2σ 2 n - 1 2 log |A| + N RF log N RF σ 2 n σ 2 0 - N 2 log 2πσ 2 n (2.108)
Since this method approximates kernel matrices using the spectral density, it is called the Sparse Spectrum Gaussian Process (ssgp). Model selection can be done by optimizing jointly the logarithm of the marginal likelihood defined in 2.108 with respect to spectral points ω(r) and hyperparameters θ. The computational cost for each training step of ssgp algorithm is O (N N 2 RF ). In terms of making prediction for each test point, the cost is O (N RF ) for the predictive mean and O (N 2 RF ) for the predictive variance.

Extensions of ssgp. [START_REF] Gal | Improving the Gaussian Process Sparse Spectrum Approximation by Representing Uncertainty in Frequency Inputs[END_REF] show that the original ssgp model's have a tendency of overfitting. They have presented a Variational Sparse Spectrum approximation to the Gaussian Processes (vssgp) that allows one to integrate out the set of spectral samples Ω = ω(1) , . . . , ω(N RF ) .

The model vssgp is shown to yield better calibrated uncertainty estimates accompanying predictions, and a procedure for deriving the optimal weights analytically is given for the Gaussian likelihood case. Other approaches of applying variational inference on ssgp are featured in [START_REF] Tan | Variational inference for sparse spectrum gaussian process regression[END_REF] and [START_REF] Hoang | A generalized stochastic variational bayesian hyperparameter learning framework for sparse spectrum gaussian process regression[END_REF]. Besides, efficient random feature maps have also been proposed to accelerate the computation and reduce the storage cost, such as the Fastfood approximation [START_REF] Le | Fastfood -Approximating Kernel Expansions in Loglinear Time[END_REF] and Orthogonal Random Features (Yu et al., 2016). as the softmax [START_REF] Jacobs | Adaptive mixture of local expert[END_REF] and probit function [START_REF] Geweke | Smoothly mixing regressions[END_REF]. More general, it can be extended to a tree-structured hierarchical architecture [START_REF] Jordan | Hierarchical mixtures of experts and the em algorithm[END_REF].

Local Approximation

The application of gps mixture experts for big data scenarios must deal with various problems. For example, the question of determining the number of local experts can be dealt with by Akaike information criterion [START_REF] Huang | Estimating mixture of gaussian processes by kernel smoothing[END_REF], or the synchronously balancing criterion (Zhao et al., 2015a). Another problem is on the reduction of computational cost, which includes several re- which can be reduced to O (nm 2 ) using hard-cut EM (Nguyen and Bonilla, 2014b;[START_REF] Nguyen | Variational inference for infinite mixtures of sparse gaussian processes through kl-correction[END_REF].

Calibrating Deep Convolutional

Gaussian Processes

The wide adoption of Convolutional Neural Networks (cnns) in applications where decision-making under uncertainty is fundamental, has brought a great deal of attention to the ability of these models to accurately quantify the uncertainty in their predictions. Previous work on combining cnns with Gaussian processes (gps) has been developed under the assumption that the predictive probabilities of these models are well-calibrated. In this paper we show that, in fact, current combinations of cnns and gps are miscalibrated.

We propose a novel combination that considerably outperforms previous approaches on this aspect, while achieving state-of-the-art performance on image classification tasks.

Introduction

The wide adoption of Convolutional Neural Networks (cnns) in increasingly popular pieces of technology such as self driving cars and medical imaging, where decision-making under uncertainty is fundamental, has brought attention to the ability of these learning architectures to accurately quantify the uncertainty in their predictions (Kendall and Gal, 2017;[START_REF] Gal | Bayesian Convolutional Neural Networks with Bernoulli Approximate Variational Inference[END_REF]. In short, the reliability of predictive probabilities of learning algorithms can be evaluated through the analysis of their calibration [START_REF] Flach | Classifier Calibration[END_REF]. In particular, a classifier is well calibrated when its output offers an accurate account of the probability of a given class, i.e. when it predicts a given class label with probability p that matches the true proportion p of test points belonging to that class.

The calibration properties of standard classifiers and neural networks have been studied in the literature [START_REF] Kull | Beta calibration: a well-founded and easily implemented improvement on logistic calibration for binary classifiers[END_REF][START_REF] Niculescu-Mizil | Predicting Good Probabilities with Supervised Learning[END_REF], which has shown that classifiers that use the standard cross-entropy loss are generally well calibrated. Perhaps surprisingly, modern cnns, which are a particular case of deep neural networks (dnns), have been found to be miscalibrated, and the depth of convolutional filters is the main factor affecting calibration [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF]. The work in [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF] shows that regularization, implemented through weight decay, improves calibration and that, ultimately, simple methods such as post-calibration [START_REF] Platt | Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods[END_REF] can be an effective remedy for most calibration issues of cnns.

Alternatively, Bayesian cnns [START_REF] Gal | Bayesian Convolutional Neural Networks with Bernoulli Approximate Variational Inference[END_REF] where convolutional filters are inferred using Bayesian inference techniques, seem like perfect candidates to model uncertainty in these architectures in a principled way.

However, while Bayesian cnns have been shown to be effective in obtaining state-of-the-art performance in image classification tasks, we are not aware of studies that show their calibration properties. Hence, our first contribution is to investigate the calibration properties of Bayesian cnns.

Along a similar vein, independently of the works on Bayesian cnns, there have been other attempts to give a probabilistic flavor to cnns by combining them with Gaussian processes (gps, [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]). Most of these approaches can be seen as a way to parameterize a cnn-based covariance for gps, and the aim is to learn end-to-end both the filters and the gps (see, e.g., [START_REF] Bradshaw | Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks[END_REF]; Wilson et al. (2016)). A crucial aspect that the literature has overlooked, however, is that methods that combine cnns and gps suffer from the same issues of miscalibration that characterize modern cnns. Therefore, the second contribution of this paper is to show that current combinations of cnns and gps are miscalibrated.

Consequently, as our third contribution, we propose a novel combination of cnns and gps that is indeed well-calibrated, while being simple to implement. In particular, we propose to replace the fully connected layers of cnns with gps that we approximate with random features [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF][START_REF] Lázaro-Gredilla | Sparse Spectrum Gaussian Process Regression[END_REF]. Due to this approximation, the resulting model becomes a Bayesian cnn with a nonlinear transformation applied to the convolutional features. Building on the connection between variational inference and dropout, we apply Monte Carlo dropout (mcd, (Gal and Ghahramani, 2016a)) to carry out joint inference over the filters and the approximate gps, thus obtaining an end-to-end learning method for the proposed model, which we call cnn+gp(rf). The resulting approach is characterized by a number of attractive features: (i) it is well calibrated, given that it uses the multinomial likelihood and the filters are regularized using Bayesian inference techniques;

(ii) it is as scalable as state-of-the-art cnns, in so much as it can be trained using mini-batch updates and can exploit GPU and distributed computing;

(iii) unlike other works that combine cnns and gps, it is as easy to implement as standard cnns, as it leverages the equivalence of gps approximated with random features and Bayesian dnns [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF][START_REF] Gal | Improving the Gaussian Process Sparse Spectrum Approximation by Representing Uncertainty in Frequency Inputs[END_REF][START_REF] Neal | Bayesian Learning for Neural Networks (Lecture Notes in Statistics)[END_REF], and the connections between dropout and variational inference (Gal and Ghahramani, 2016a). We extensively validate these properties in a variety of image classification tasks.

Our final contribution extends the above framework by replacing the last layer of cnns with Deep gps [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF] and by proposing the use of structured random features to obtain faster and more compact gp approximations [START_REF] Le | Fastfood -Approximating Kernel Expansions in Loglinear Time[END_REF]Yu et al., 2016). In all, our proposal considerably improves on classification accuracy compared to previous combinations of cnns and gps (e.g., ∼88% on cifar10 and ∼67% on cifar100, all without data augmentation), while being competitive with state-of-the-art cnns; we are not aware of other gp works that approach these results. Crucially, we achieve these performance without compromising on calibration, again considerably improving on previous approaches that combine cnns and gps.

Related Work

Calibration of Convolutional Networks: The issue of calibration of classifiers in machine learning was popularized in the 90's with the use of support vector machines for probabilistic classification [START_REF] Platt | Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods[END_REF]. Calibration techniques aim to learn a transformation of the output using a validation set in order for the transformed output to give a reliable account of the actual probability of class labels [START_REF] Flach | Classifier Calibration[END_REF]; interestingly, calibration can be applied regardless of the probabilistic nature of the untransformed output of the classifier. Popular calibration techniques include Platt scaling [START_REF] Platt | Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods[END_REF] and isotonic regression [START_REF] Zadrozny | Transforming Classifier Scores into Accurate Multiclass Probability Estimates[END_REF].

Classifiers based on Deep Neural Networks (dnns) have been shown to be well-calibrated [START_REF] Niculescu-Mizil | Predicting Good Probabilities with Supervised Learning[END_REF]. The reason is that the optimization of the cross-entropy loss promotes calibrated output. The same loss is used in Platt scaling and it corresponds to the correct multinomial likelihood for class labels. Recent sudies on the calibration of cnns, which are a particular case of dnns, however, show that depth has a negative impact on calibration, despite the use of a cross-entropy loss, and that regularization improves the calibration properties of classifiers [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF]. et al., 2017), the left and central panels show the same curve for cgp.

Combinations of Conv Nets and Gaussian Processes: Thinking of

Bayesian priors as a form of regularization, it is natural to assume that Bayesian cnns can "cure" the miscalibration of modern cnns. Despite the abundant literature on Bayesian dnns [START_REF] Neal | Bayesian Learning for Neural Networks (Lecture Notes in Statistics)[END_REF][START_REF] Mackay | Bayesian methods for backpropagation networks[END_REF], far less attention has been devoted to Bayesian cnns (Gal and Ghahramani, 2016a), and the calibration properties of these approaches have not been investigated.

Several approaches have proposed the combination of cnns and gps as a means to give a probabilistic character to cnns. Most of these works are based on ideas developed in the context of manifold gps [START_REF] Calandra | Manifold Gaussian Processes for regression[END_REF], where inputs are transformed using some parametric transformation.

In these works, the parametric transformation is based on convolutional layers, and scalability to large data is achieved through the use of ideas drawn from the literature on scalable gps, for example the Stochastic Variational Deep Kernel Learning (svdkl) approach in Wilson et al. (2016). In contrast, the work on hybrid gps and dnns (gpdnn, [START_REF] Bradshaw | Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks[END_REF]) combines cnns and gps using an inducing point approximation. Other recent approaches that aim to introduce convolutions in the calculation of the covariance between images include the work in van der Wilk et al. ( 2017), which proposes a way to construct covariances between domains/patches, mimicking the computations in cnns.

In this work, we propose an alternative way to combine cnns and gps, where gps are approximated using random features expansions (Rahimi and Recht, 2008;[START_REF] Lázaro-Gredilla | Sparse Spectrum Gaussian Process Regression[END_REF]. The random feature expansion approximation amounts to replacing the orginal kernel matrix with a low-rank approximation, turning gps into Bayesian linear models. Combining this with cnns leads to a particular form of Bayesian cnns, much like gps and dgps are particular forms of Bayesian dnns [START_REF] Duvenaud | Avoiding pathologies in very deep networks[END_REF]Gal and Ghahramani, 2016a;[START_REF] Neal | Bayesian Learning for Neural Networks (Lecture Notes in Statistics)[END_REF]. Inference in Bayesian cnns is intractable and requires some form of approximation. In this work, we draw on the interpretation of dropout as variational inference, employing the so-called Monte Carlo Dropout (mcd, (Gal and Ghahramani, 2016a)) to obtain a practical way of combining cnns and gps.

On calibration of Convolutional GPs

Consider a Q-class image classification task where X denotes a set of N images

x i ∈ R px×py (1 ≤ i ≤ n),
and Y is the matrix consisting of the corresponding one-hot encoded labels y i stacked by row. We can use various metrics to determine the quality of a classifier, and here we focus in particular on calibration.

Let g(x) be the output of a classifier for an input image x. To compute the calibration properties of a classifier, consider a partitioning of the test set X * into disjoint sets {X 1 , . . . , X M }, such that each subset X m contains the inputs yielding predictions in the range ( m-1 M , m M ]. Hence, the confidence associated with each subset X m is characterized by the midpoint of its corresponding range, i.e. conf(X m ) = m-0.5 M . Then, the accuracy acc(X m ) for each subset can be evaluated as follows:

1 |X m | x * ∈Xm δ (arg max(y * ) -arg max(g(x * ))) , (3.1)
where δ(x) is equal to one if x = 0, and zero otherwise.

In what follows, we use reliability diagrams to assess calibration, where we plot accuracy as a function of confidence for the subsets {X 1 , . . . , X M }. For a perfectly calibrated classifier, we expect acc(X m ) = conf(X m ) for all m, with deviations implying that the class probabilities are either underestimated or overestimated. A useful summary statistics that can be extracted from reliability diagrams is the Expected Calibration Error (ece), which is the average of the absolute difference between accuracy and confidence weighted according to its size:

ece = M m=1 |X m | |X * | |acc(X m ) -conf(X m )| . (3.2)
Another metric that measures the accuracy in predicting class probabilities is the brier score which takes into account the factors of calibration, resolution and uncertainty [START_REF] Murphy | A new vector partition of the probability score[END_REF]. It is defined as the squared distance between labels and outputs averaged across classes and test points:

brier = 1 N test x * ∈X * 1 Q Q k=1 ((y * ) k -(g(x * )) k ) 2 . (3.3)
In figure 3.1, we report the reliability diagrams of three state-of-the-art combinations of cnns and gps, i.e gpdnn approach in [START_REF] Bradshaw | Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks[END_REF], cgp in van der Wilk et al. ( 2017) and svdkl in Wilson et al. (2016). These approaches are applied to the cifar10 and cifar100 data sets with various convolutional structures. Note that the lines for cgp in the sub-figure of cifar10-LeNet and cifar10-resnet are identical because there is no equivalent cnn architecture in cgp. All of reliability diagrams for these methods and ours can be found in the supplemental material.

The results indicate that current approaches that combine cnns and gps are miscalibrated, with a tendence of being overconfident in predictions. This is an important and perhaps surprising finding, because one of the motivations to combine cnns with gps is to do better quantification of uncertainty compared to plain cnns. In the experiments section we report more extensively on the calibration of these classifiers, as well as illustrating other performance metrics. These considerations call for the study of better ways to combine cnns and gps to recover calibration while attempting to improve on standard metrics such as error rate and test log-likelihood. The next section illustrates our proposal that achieves this goal.

Proposed Method

In the proposed model, the labels Y i• are assumed to be conditionally independent given a set of corresponding latent variables F i• , i.e. we consider the likelihood p(Y|F)

= N i=1 p(Y i• |F i• )
, where the latent variables F are realizations of a set of Q functions f j (x) at the input images x 1 , . . . , x n , i.e.,

(F) ij = f j (x i ) for j = 1, . . . , Q. Each individual p(Y i• |F i• ) is multinomial
with probabilities obtained using a softmax transformation of the latent variables. In this work we focus on functions f j (x) that are modeled using gps;

note that extension to dgps is actually easy to consider in our framework, as we show in the experiments.

Due to the gp modeling assumption, the latent function values F •j comprising (f j (x 1 ), . . . , f j (x n )) are jointly Gaussian with p(F

•j |X, θ) ∼ N (0, K),
where K is the covariance matrix. The entries of the covariance matrix pooling. This is one of the key successes of cnn models that allows for the learning of their filters, which we exploit for the end-to-end learning of our model.

K = {k (x i , x j |θ)} i,j ,
Inference in this model requires being able to characterize the posterior over all or a selected group of model parameters, but this posterior is analytically intractable and thus computationally prohibitive [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. In the remainder of this paper, we build on previous work on scalable inference for gps and dgps with random features [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF] to obtain an approximation to the proposed model that can be learned end-to-end.

Random Feature Expansions

Naïve inference in gp models requires algebraic operations with K that would cost O(n 3 ) in time. Popular approaches to recover tractability use low-rank approximations of the kernel matrix. Among this family of low-rank approximations, we choose to work with random feature approximations [START_REF] Lázaro-Gredilla | Sparse Spectrum Gaussian Process Regression[END_REF][START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF]. The reason is that they offer a number of possible extensions to speedup computations (e.g., using structured approximations [START_REF] Le | Fastfood -Approximating Kernel Expansions in Loglinear Time[END_REF]Yu et al., 2016)) and increase the complexity of the model (e.g., considering Deep gps [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF])); we elaborate on this in the experiments section. In random feature expansions, the kernel matrix is replaced by a low-rank approximation K ≈ ΦΦ , with Φ ∈ R n×m and m n. This approximation suggests the construction of a Bayesian linear model to approximate the gp latent variables as F = ΦW.

Using p(W ij ) = N (W ij |0, 1) it is straightforward to show that the covariance of each of the latent functions F •j is indeed an approximation to K, as

Φ arc = 2σ 2 N RF max (0, C(X|Ψ) Ω) . (3.5)
In this expression, we have defined C(X|Ψ) as the matrix resulting from the application of convolutional layers to the image training set X and Ω is obtained by stacking N RF samples from p(ω) = N ω|0, Λ -1 by column.

Note that in the case of a popular Radial Basis Function (rbf) covariance, it is possible to obtain a similar random feature approximation, where the relu activation is replaced by trigonometric functions; see Rahimi and Recht (2008) and the supplement for details.

End-to-end learning

Inference in the proposed model is intractable due to the likelihood that is not conjugate to the gp prior. Further complications stem from the need to infer kernel parameters, which include convolutional parameters, and the need to be able to scale to large data. Our aim is to carry out inference within a consistent framework that is characterized by simplicity, as described next.

We start by introducing an approximate posterior over W, Ω and Ψ, that we denote as q(W, Ω, Ψ). Following standard variational inference arguments, we can define an operative way to obtain these approximate posteriors. The log-marginal likelihood L = log [p(Y|X, θ] can be bounded by the sum of an expected log-likelihood term and a negative Kullback-Leibler (KL) divergence term as follows:

L ≥E q(W,Ω,Ψ) (log [p (Y|X, W, Ω, Ψ, θ)]) -KL [q (W, Ω, Ψ) p (W, Ω, Ψ)] . (3.6)
Variational inference amounts to optimizing the lower bound above with respect to q(W, Ω, Ψ) and any other parameters of interest.

We have now a number of options on the form for the approximate posteriors q(W, Ω, Ψ). In previous works on variational inference for dnns, it has been proposed to define the approximating distributions to be Gaussian and factorized across parameters [START_REF] Kingma | Auto-Encoding Variational Bayes[END_REF]Graves, 2011). The drawback of this is that it doubles the number of parameters. Alternatively, we can rely on the connections between dropout and variational inference (Gal and Ghahramani, 2016a,b) which is drawn by assuming the posterior of W, Ω and Ψ as a mixture of two Gaussian distributions (see supplement).

From this connection, we are able to obtain an easier approximate inference scheme, which is also known as Monte Carlo Dropout (mcd). Focusing on the weights for now, the connection with dropout is apparent if we rewrite

W = M w Diag[z w ] (3.7)
with (z w ) i ∼ Bernoulli(π w ). The reparameterization introduces variational parameters M w (one for each weight in W) and a vector of binary variables that can switch on or off the columns of the weight matrix with probability π w . A similar reprameterization can be done for the convolutional parameters Ψ and matrices of random feature Ω, introducing M ψ , M Ω and π ψ , π Ω . The optimization of the lower bound wrt all variational parameters requires being able to evaluate the expectation and the KL term in (3.16).

In mcd, the KL term in (3.16) can be approximated following Gal and Ghahramani (2016a), obtaining a regularization term involving the squared-norm of the parameters

KL [q (W, Ω, Ψ) p (W, Ω, Ψ)] ≈ π w 2 M w 2 + π Ω 2 M Ω 2 + π ψ 2 M ψ 2 (3.8)
The expectation in (3.16), instead, can be unbiasedly estimated using Monte

Carlo and also considering a mini-batch of size m:

N m 1 N MC N MC i=1 k∈Im log p y k |x k , W (i) , Ω (i) , Ψ (i) , θ (3.9)
with W (i) , Ω (i) , Ψ (i) ∼ q(W, Ω, Ψ), and I m is a set of m indices to select a mini-batch of training points (Graves, 2011). This doubly-stochastic approximation is differentiable wrt variational parameters when the Bernoulli variables are fixed.

The approximate objective can now be optimized in the same vein as in standard back-propagation with dropout, noting that dropout is applied to W, Ω and to convolutional parameters Ψ. What changes, however, is the interpretation of the procedure as stochastic variational inference, whereby the Bernoulli variables are resampled at each iteration. A practical implication is in the way we compute the predictive distribution, which has a probabilistic flavor as follows:

p(y * |x * , X, θ) ≈ p(y * |W, Ω, Ψ, x * , X, θ)q(W, Ω, Ψ)dWdΩdΨ, (3.10)
and can be approximated using Monte Carlo by resampling the Bernoulli variables. While mcd has been proposed for cnns in [START_REF] Gal | Bayesian Convolutional Neural Networks with Bernoulli Approximate Variational Inference[END_REF], in this work we extend it to the case of joint inference over convolutional parameters and the gp approximation in the cnn+gp(rf) model, thus obtaining a practical inference and prediction scheme, which combines cnns and gps.

Depth

Extensions

Structured random feature approximations: One of the advantages of the proposed model, compared to other gp approximations, is that it can exploit structured random feature expansions to accelerate computations and reduce the size of the approximate gp [START_REF] Le | Fastfood -Approximating Kernel Expansions in Loglinear Time[END_REF]Yu et al., 2016). In the random features approximation, random features are constructed by multiply-

ing Ω with the convolutional features. Without loss of generality, assuming

that Ω ∈ R m×d and c(x|Ψ) ∈ R d×1 , the cost of computing products Ωc(x|Ψ)

is O (md), while storing Ω requires O (md) storage. with G ij ∼ N (0, 1). One way to make computations cheaper is to replace the Gaussian matrix G with a pseudo-random alternative. The Structured Orthogonal Random Feature (sorf) approximation (Yu et al., 2016) approximates G through a series of Hadamard transformations of diagonal matrices

D i with elements randomly sampled from {-1, +1} or Rademacher distri- bution, that is G ≈ √ dHD 1 HD 2 HD 3
, where H is the normalized Walsh-Hadamard matrix. We refer to this variation of the model as cnn+gp(sorf).

Similarly to the other parameters, we infer the diagonal matrices D i using mcd. We denote by d i the diagonal of D i , i = 1, 2, 3. The mcd scheme (Gal and Ghahramani, 2016a,b) A dgp model represents a deep probabilistic nonparametric approach where the output of one gp at each layer is used as the input to the gp in the next layer [START_REF] Damianou | Deep Gaussian Processes[END_REF]. Extending the random feature approximation to dgps and the inference scheme presented here is straightforward;

see [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF] for details. The random feature approximation turns the dgp into a Bayesian dnn for which we can apply stochastic variational inference to infer model parameters. In the experiments section, we explore the possibility to stack a dgp on top of convolutional layers, and we show the impact of depth on performance.

Experiments

We carry out the experimental evaluation using popular benchmark datasets, such as mnist, cifar10 and cifar100 and with a number of popular cnn architectures based on LeNet and resnet (see table 3.1).

We report three state-of-the-art competitors combining cnns and gps, namely gpdnn [START_REF] Bradshaw | Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks[END_REF], svdkl (Wilson et al., 2016), and cgp (van der [START_REF] Van Der Wilk | Convolutional Gaussian Processes[END_REF]. We also report Bayesian cnns, as suggested in [START_REF] Gal | Bayesian Convolutional Neural Networks with Bernoulli Approximate Variational Inference[END_REF] and cnns with post-calibration as proposed in [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF], which we refer to as cnn+mcd and cnn+cal, respectively. For all the competing methods we used available implementations, adding the same cnn architecture to ensure a fair comparison. In all experiments, we use a batch-size m = 100 and the Adam optimizer with default learning rate [START_REF] Kingma | A method for stochastic optimization[END_REF]. In the methods that use mcd, we use a dropout rate of 0.5 for all parameters.

The results are reported in figure 3.2, where we have used different training sizes N , keeping the classes balanced. In the figure, we report the calibration measures that we have introduced earlier, namely ece and brier scores, and we also report the classification error rate (err) and the mean negative test log-likelihood (mnll). Compared to other combinations of cnns and gps, cnn+gp(rf) improves considerably on all metrics. It is interesting to see that our proposal is competitive with Bayesian cnns employing mcd, with only a marginal improvement on err and mnll in some configurations.

In temp it is necessary to leave out part of the data to perform post-calibration, which can be problematic in applications where obtaining labeled data is difficult or expensive. As a result, our proposal is considerably better, although temp is competitive in ece; this is expected given that this is the metric that is optimized after training.

The two variants of our approach, namely cnn+gp(rf) where we learn the frequencies Ω and cnn+gp(sorf) where we sample Ω from its prior, are comparable. This suggests that the extra level of complexity of learning the spectral frequencies does not lead to substantial gains in performance and that the structured random feature approximation yields satisfactory performance.

We also note that these results have been obtained by fixing the covariance parameters θ of the gp, as we found it to be unstable when learning these jointly with Ω. This might be the reason why these parameters were learned through cross-validation in Gal et al. (2017). In the supplement, we report the results obtained when learning θ and fixing Ω, which we found yielding similar performance as fixing θ. All these observations corroborate the hypothesis that most of the performance of cnn-based classification models is due to the convolutional layers.

In summary, figure 3.2 shows that our cnn+gp(rf) is the best strategy for calibrating these models compared to other approaches using gps. Furthermore, we found perhaps surprisingly that mcd has comparable performance.

In the supplementary material, we report results on gpdnn where we infer convolutional parameters using mcd, so as to gain insights as to whether most of the improvements in performance are due to this form of regularization. The results support the intuition that inferring these parameters yields improvements in calibration, but also that our cnn+gp(rf) still offers better performance.

Reliability diagrams

In figure 3.3, we report the reliability diagrams of all the methods studied in figure 3.1. The figure shows that temp, mcd and cnn+gp(rf) produce wellcalibrated predictions when using a shallow convolutional structure (LeNet).

For a deeper architecture (resnet), cnn+gp(rf) is slightly under-confident.

Compared to previous combinations of cnns and gps, our approach yields better reliability curves.

Extension with Deep gps

In figure 3.4, we report results varying the depth of a dgp on top of the convolutional layers; again, we learn the convolutional filters and the dgp end-to-end as discussed in the previous sections. We show results when applying our model to the whole cifar10 data set in the case of the shallow convolutional structure (table 3.1). We feed-forward the convolutional features to all layers of the dgp, in line with what suggested in the literature of dgps to avoid pathologies in the functions that can be modeled [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF][START_REF] Duvenaud | Avoiding pathologies in very deep networks[END_REF][START_REF] Neal | Bayesian Learning for Neural Networks (Lecture Notes in Statistics)[END_REF]. The results indicate that increasing the complexity of the model improves on all performance metrics, and worsen calibration, which however is still around 3% ece. This is in line with the intuition that increasing model complexity negatively impacts calibration.

Knowing when the model does not know

We report experiments showing the ability of our model to know when it does not know, following a similar experimental setup as in Lakshminarayanan et al. (2017). In this experiment we train our cnn+gp(rf) model on mnist and test on the not-mnist dataset, which contains images of letters from "A" to "J" in various typefaces. For this experiment, while we do not know the exact value that we should obtain for predictive probabilities, we expect to observe low entropy in the predictions when tesing on mnist and high entropy when predicting on not-mnist, indicating high uncertainty. The results are reported in figure 3.5, where we show the cumulative distribution of the entropy of predictive probabilities for two depths of the convolutional structure. In the figure, we compare our cnn+gp(rf) against one of the methods combining cnns and gps, that is gpdnn. In the figure, we also include results on cnns with post-calibration and Bayesian cnns inferred with mcd. Our approach is competitive with Bayesian cnns and it is considerably superior to post-calibration. This is especially true in the case of the resnet convolutional structure, where post-calibration still yields a large number of predictions with low uncertainty. Interestingly, gpdnn assigns large uncertainty to predictions on not-mnist, although with the deeper convolutional architecture it yields a large fraction of predictions with low entropy. We speculate that this due to the inducing point approximation of the gp, which nicely captures uncertainty away from training data except for test points which are closer to the training data.

Extension with the sorf

In table 3.2, we report further results comparing mcd with cnn+gp(sorf).

In this experiment, we use the alexnet structure [START_REF] Krizhevsky | ImageNet classification with deep convolutional neural networks[END_REF] on cifar10 and cifar100 datasets. The results in table 3.2 show improvements in using our model compared cnns with mcd. We attribute this to the fact that the gp approximated through sorf in place of the fully connected layer of alexnet reduces model parameters from 30 million to 2.3 million.

Mathematical details and other experiments

Random Feature Expansion of the rbf Covariance

We report here the expansion of the popular Radial Basis Function (rbf) covariance. Following the convolutional representation of images in our cnn+gp(rf) model, the rbf covariance is defined as:

k rbf (x i , x j |Ψ, θ) = σ 2 exp -(c(x i |Ψ) -c(x j |Ψ)) Λ -1 (c(x i |Ψ) -c(x j |Ψ)) , (3.12)
with θ = (σ, Λ = Diag( 2 1 , . . . , 2 d )). It is possible to express this covariance function as the Fourier transform of a non-negative measure p(ω) Rahimi and Recht (2008), where ω are the so-called spectral frequencies. It is straightforward to verify that p (ω) = N ω|0, Λ -1 . Stacking N RF Monte Carlo samples from p(ω) into Ω by column, we obtain

Φ rbf = σ 2 N RF [cos (C(X|Ψ) Ω) , sin (C(X|Ψ) Ω)] , (3.13)
where C(X|Ψ) denotes the matrix resulting from the application of convolutional layers to the image training set X, and the sin and cos functions are applied elementwise to their argument.

3.6.2 Variational Inference for the Proposed Model

3.6.2.1 cnn+gp(rf)
In cnn+gp(rf), the variational parameters we would like to optimize are M w , M ψ and M Ω . Our model parameters W, Ψ and Ω share an identical form for the approximate posterior and prior. Focusing on W, its elements have a standard normal prior, and we assume that the posterior q (W) is a mixture of two Gaussian distribution, which can be factorized over rows, governed by variational parameters M w :

q (W) = R r=1 q (W r ) , with q (W r ) = π w N M wr , σ 2 I D +(1-π w )N 0, σ 2 I D ,
(3.14) where π w ∈ [0, 1], σ 2 ≈ 0 and M wr ∈ R D . This form of posterior leads to the sampling procedure which characterizes dropout Gal and Ghahramani (2016a,b). Given the choice of σ 2 ≈ 0, W can be sampled by introducing Bernoulli variables

W = M w Diag[z w ] with (z w ) i ∼ Bernoulli(π w ), (3.15)
and similarly for Ψ and Ω.

All variational parameters are optimized to maximize the lower bound of marginal likelihood which is defined as follows

log [p(Y|X, θ] ≥E q(W,Ψ,Ω) (log [p (Y|X, W, Ψ, Ω, θ)]) -KL [q (W, Ψ, Ω) p (W, Ψ, Ω|θ)] (3.16)
The expectation in 3.16 can be unbiasedly estimated using Monte Carlo and also considering a mini-batch of size m

E q(W,Ψ,Ω) (log [p (Y|X, W, Ψ, Ω, θ)]) ≈ N m 1 N MC N MC i=1 k∈Im log p y k |x k , W (i) , Ψ (i) , Ω (i) , θ , (3.17)
where W (i) , Ψ (i) , Ω (i) is a sample from q(W, Ψ, Ω), and can be obtained via 3.15. I m is a set of m indices to select a mini-batch of training points. In classification, each individual p y k |x k , W (i) , Ψ (i) , Ω (i) , θ can be computed using a softmax transformation. The KL term can be approximated following Gal and Ghahramani (2016a), noting that the fact that we are treating Ω

variationally, gives rise to extra terms that involve the gp length-scale :

KL [q (W, Ψ, Ω) p (W, Ψ, Ω|θ)] ≈ π w 2 M w 2 + π ψ 2 M ψ 2 + 2 π Ω 2 M Ω 2 + N RF d log -2 (3.18) 3.6.2.2 cnn+gp(sorf)
In cnn+gp(sorf), our proposed variational inference scheme is similar to the one in cnn+gp(rf), except that Ω is replaced by l -1 √ N RF HD 1 HD 2 HD 3 , with length-scale l and D i = Diag (d i ) and H is the normalized Walsh-Hadamard matrix. Because d i is Rademacher distributed, the form of prior and posterior in mcd proposed by Gal and Ghahramani (2016a,b) is inadequate. Therefore, we use the prior

p ε (d i ) = N (d i |d * i , ε 2 I N RF ) with d * i sampled
from the Rademacher distribution and a small positive ε. The posterior q (d)

is also composed by two Gaussian distribution as in cnn+gp(rf) Gal and Ghahramani (2016a), we can approximate the KL term between q (d i ) and p

q (d i ) = N RF j=1 q [d i ] j , where q [d i ] j = π d N M [d i ] j , σ 2 + (1 -π d ) N [d * i ] j , σ 2 (3.19) with π d ∈ [0, 1] , σ 2 ≈ 0 and M d i ∈ R N RF . Following
(d i ) KL (q (d i ) p ε (d i )) ≈ π d 2ε 2 M d i -d * i 2 (3.20)
In terms of implementation, we do not apply mcd to d i -d * i but on d i directly.

According to this choice, each element in d i is sampled based on the variational parameters M d i -d * i as in 3.21. Thanks to this trick, the implementation of mcd scheme does not change for optimizing d i

d i = M d i -d * i + d * i , with probability π d d * i , otherwise (3.21)
In figure 3.6, we report some experimental results to illustrate the impact of optimizing d i . For cifar10-LeNet and cifar100-resnet, the optimization of sorf parameters outperforms the case where spectral frequencies are fixed in terms of err, mnll and brier. In the case of cifar10-resnet, the gains are marginal.

Optimization for covariance parameters

When using 3.18 to optimize all variational parameters pertaining to q (W, Ψ, Ω)

jointly with covariance θ we encountered some instabilities, and therefore we decided to report results when fixing the covariance parameters θ in our paper. For the case where Ω is not learned variationally we can simply draw Ω from the prior N (Ω •j |0, Λ -1 ) and consider the reparameterization:

Ω •j = Λ -1 2 ε, (3.22)
where [START_REF] Lázaro-Gredilla | Sparse Spectrum Gaussian Process Regression[END_REF]. This reparameterization allows for the update of covariance parameters θ fixing the randomness in the sampling from p(Ω|θ). The results comparing cnn+gp(sorf) when updating or fixing θ throughout optimization are reported in table 3.3. It is interesting to notice how fixing covariance parameters θ leads to comparable performance to the case where they are learned.

ε i ∼ N (ε i |0, 1) (Lázaro

Variational inference of filters in gpdnn

In this section we report results when applying variational inference on the weights in gpdnn [START_REF] Bradshaw | Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks[END_REF]. In order to do this, we implemented mcd for the convolutional parameters, similarly to what presented in the main paper for our cnn+gp(rf) model. The results in table 3.4 indicate that this improves the calibration and accuracy of gpdnn compared to optimizing the filters. In the case of a shallow convolutional architecture, the performance of cnn+gp(rf) and gpdnn are comparable, although in the deeper case cnn+gp(rf) achieves better performance. This supports the intuition that inferring convolutional parameters, ranther than optimizing them, leads to considerable improvements in calibration.

Reliability diagrams

In this section, we report the reliability diagram and histogram of predictive output for all methods with various datasets, i.e cifar10 and cifar100 and convolutional architectures, i.e LeNet and resnet. We use the best config- Having observed these figures, we see that regularizing convolutional filters has a huge impact on calibration. From figures 3.7, 3.8, 3.9 and 3.10 we see that cnns and the previous combinations of gps and cnns are miscalibrated.

From figure 3.12 and 3.13, instead, we see that Bayesian cnns improve the reliability of the prediction, which is comparable with post-calibration.

It seems that there is a correlation between the histogram of predictive output and the reliability line. When the histogram is skewed to the right, the corresponding classifier is poorly calibrated.

Local and Global Approximation of Gaussian Processes

Approximations to Gaussian processes (GPs) based on inducing variables, combined with variational inference techniques, enable state-of-the-art sparse approaches to infer GPs at scale through mini-batch-based learning. In this work, we address one limitation of sparse GPs, which is due to the challenge in dealing with a large number of inducing variables without imposing a special structure on the inducing inputs. In particular, we introduce a novel hierarchical prior, which imposes sparsity on the set of inducing variables.

We treat our model variationally, and we experimentally show considerable computational gains compared to standard sparse GPs when sparsity on the inducing variables is realized considering the nearest inducing inputs of a random mini-batch of the data. We perform an extensive experimental validation that demonstrates the effectiveness of our approach compared to the state-ofthe-art. Our approach enables the possibility to use sparse GPs using a large number of inducing points without incurring a prohibitive computational cost.

Introduction

Gaussian Processes (gps) [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]) offer a powerful framework to perform inference over functions; being Bayesian, gps provide rigorous uncertainty quantification and prevent overfitting. However, the applicability of gps on big datasets is hindered by their computational complexity of O (N 3 ), where N is the training size. This issue has fuelled a considerable amount of research towards scalable gp methodologies that operate on a set of inducing variables [START_REF] Pratola | A unifying view of sparse approximate Gaussian process regression[END_REF]. In the literature, there is a plethora of approaches that offer different treatments of the inducing variables [START_REF] Lawrence | Fast Sparse Gaussian Process Methods: The Informative Vector Machine[END_REF]Seeger et al., 2003;[START_REF] Snelson | Sparse Gaussian Processes using Pseudoinputs[END_REF]Naish-Guzman and Holden, 2007;[START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF][START_REF] Hensman | Gaussian processes for big data[END_REF]Wilson and Nickisch, 2015;Hensman et al., 2015a). Some of the more recent approaches, such as Scalable Variational Gaussian Processes (svgps) (Hensman et al., 2015a), allow for the application of gps to problems with millions of data points. In most applications of scalable gps, these are approximated using M inducing points (IPs) , which results in a complexity of O (M 3 ). It has been shown recently by [START_REF] Burt | Rates of convergence for sparse variational Gaussian process regression[END_REF] that it is possible to obtain an arbitrarily good approximation for a certain class of gp models (i.e.

conjugate likelihoods, concentrated distribution for the training data) with M growing more slowly than N . However, the general case remains elusive and it is still possible that the required value for M may exceed a certain computational budget. Our result contributes to strengthen our belief that sparsity does not only enjoy desirable theoretical properties, but it also constitutes an extremely computationally efficient method in practice.

In this work, we push the limits of scalability and effectiveness of sparse gps enabling a further reduction in complexity, which can be translated to higher Hierarchical priors are often applied in Bayesian modeling to achieve compression and to improve flexibility [START_REF] Molchanov | Variational dropout sparsifies deep neural networks[END_REF]Louizos et al., 2017). To the best of our knowledge, this work is the first to explore these ideas for the purposes of sparsifying the inducing set in sparse gps.

Related work and background

Sparse gps that operate on inducing inputs have been extensively studied in the last 20 years [START_REF] Csató | Sparse on-line gaussian processes[END_REF][START_REF] Lawrence | Fast Sparse Gaussian Process Methods: The Informative Vector Machine[END_REF][START_REF] Snelson | Sparse Gaussian Processes using Pseudoinputs[END_REF][START_REF] Pratola | A unifying view of sparse approximate Gaussian process regression[END_REF]Naish-Guzman and Holden, 2007). Many attempts on sparse gps specified inducing inputs by satisfying certain criteria that produce an informative set of inducing variables [START_REF] Csató | Sparse on-line gaussian processes[END_REF][START_REF] Lawrence | Fast Sparse Gaussian Process Methods: The Informative Vector Machine[END_REF]Seeger et al., 2003).

A different treatment has been proposed by [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF], which involves formulating the selection of inducing inputs as optimization of a variational lower bound to the marginal likelihood. The variational framework was later expanded so that stochastic optimization can be admitted, thus improving scalability for regression [START_REF] Hensman | Gaussian processes for big data[END_REF] and classification (Hensman et al., 2015a). In a more recent work [START_REF] Panos | Fully Scalable Gaussian Processes using Subspace Inducing Inputs[END_REF] scalability is addressed in terms of the dimensionality of the input. All the aforementioned methodologies share a computational complexity of O (M 3 ). Although there have been some attempts in the literature to infer the appropriate number of inducing points as well as the inducing inputs (Pourhabib et al., 2014a;[START_REF] Burt | Rates of convergence for sparse variational Gaussian process regression[END_REF], a large number of inducing variables is desirable in improving the approximation to the posterior. In this work we present a methodology that builds on the svgp framework (Hensman et al., 2015a) and reduces its complexity, thus increasing the potential of sparse gp application on even larger datasets and with a larger set of inducing variables.

A different approach to scalable gps was introduced by Wilson and Nickisch (2015), namely Kernel Interpolation for Scalable Structured gps (kiss-gp).

This line of work involves arranging a large number of inducing inputs into a grid structure; this allows one to scale to very large datasets by means of fast linear algebra. The applicability of kiss-gp on higher-dimensional problems has been addressed by Wilson et al. (2015) by means of low-dimensional projections. A more recent extension allows for a constant-time variance prediction using Lanczos methods [START_REF] Pleiss | Constant-time predictive distributions for Gaussian processes[END_REF]. Our work takes a different approach by keeping the gp prior intact, and by imposing sparsity on the set of inducing variables.

Local approximation of gps inspired by the the concept of divide-and-conquer is also a practical solution to implement scalable gps [START_REF] Kim | Analyzing nonstationary spatial data using piecewise gaussian processes[END_REF]Urta-sun and Darrell, 2008;Datta et al., 2016;?;?) which allows gps to work on large-scale datasets. In our work, we use neighbour information in a different way, by incorporating it in a certain hierarchical structure of the auxiliary variables through a variational scheme.

Scalable Variational Gaussian Processes

Consider a supervised learning problem with inputs X = (x 1 , . . . , x N ) associated with labels y = (y 1 , . . . , y N ) . Given a set of latent variables f = (f 1 , . . . , f N ) , gp models assume that labels are stochastic realizations based on f and a likelihood function p(y | f ). In svgps, the set of inducing points is characterized by inducing inputs Z = (z 1 , . . . , z M ) and inducing variables u = (u 1 , . . . , u M ) . Regarding f and u, we have the following joint prior:

p(f , u) = N 0, K X K X,Z K Z,X K Z , (4.1)
where K X , K Z and K X,Z are covariance matrices evaluated at the inputs indicated by the subscripts. The posterior over inducing variables is approximated by a variational distribution q (u) = N (u | m, S), while keeping the above and by applying Jensen's inequality:

exact conditional p(f | u) intact, that is q(f , u) = p(f | u)q(u).
E q(f ) log p (y | f ) -KL (q (u) p (u)) . (4.2)
The approximate posterior q (f ) can be computed by integrating out u: q (f ) = q (u) p (f |u) du. Thanks to the Gaussian form of q (u), q (f ) can be computed analytically:

q (f ) = N (f | Am, K X + A (S -K Z ) A) , (4.3)
where A = K X,Z K -1 Z . When the likelihood factorizes over training points, the lower bound can be re-written as:

N i=1 E q(f i ) [log p (y i | f i )] -KL (q (u) p (u)) . (4.4)
Each term of the one-dimensional expectation of the log-likelihood can be computed by Gauss-Hermite quadrature for any likelihoods (and analytically for the Gaussian likelihood). The KL (q (u) p (u)) term can be computed analytically given that q (u) and p (u) are both Gaussian. To maintain positivedefiniteness of S and perform unconstrained optimization, S is parametrized as S = LL T , with L lower triangular.

Sparse within Sparse Gaussian Processes

We present a novel formulation of sparse gps, which permits the use of a random subset of the inducing points with little loss in performance. We introduce a set of binary random variables w ∈ {0, 1} M to govern the inclusion of inducing inputs Z and the corresponding variables u. We then employ these random variables to define a hierarchical structure on the prior as follows:

p (u | w) = N (0, D w K Z D w ) , (4.5)
where D w = Diag (w) , and w ∼ p (w). Although the marginalized prior p(u)

is not Gaussian, it is possible to use the joint p(u, w) = p(u | w) p(w) within a variational scheme. We thus consider a random subset of the inducing points during the evaluation of the prior in the variational scheme that follows; no inducing points are permanently removed. Regarding p(w), we consider an implicit distribution: its analytical form is unknown, but we can draw samples from it. Later, we will consider p(w) based on the nearest inducing inputs to random mini-batches of data.

p (f , u | w) = N     0 0 0   ,   K X K X,Z I 0 K Z I ,X K Z I 0 0 0 0     (4.6)
The rows and columns of u J can simply be ignored. Regardless of the value of w, the conditional f , u I | w is always a Gaussian marginal, as it is a subset of Gaussian variables. The marginalized p(f , u) = p(f , u | w) p(w)dw is mixture of Gaussian densities, where the marginal over f is the same for every component of the mixture.

The effect on f is demonstrated in Of course, although the prior remains unchanged, that is not the case for the posterior approximation. It is well known that the choice of inducing inputs has an effect on the variational posterior [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF][START_REF] Burt | Rates of convergence for sparse variational Gaussian process regression[END_REF]. Our choice to impose a hierarchical structure to the inducing variables through w effectively changes the model compared to svgp, and we adapt the variational scheme accordingly.

Lower bound on marginal likelihood

By introducing u, w and using Jensen's inequality, the lower bound on log p (y)

can be obtained as follows E q(u,w) log p (y |u, w ) -KL (q (u, w) p (u, w)) , (4.7)

where we choose the variational distribution q to reflect the hierarchical structure of the prior, i.e. q (u, w) = q(u | w) p(w). This choice enforces sparsity over the approximate posterior q; the variational parameters are shared among the conditionals q(u | w), for which we assume:

q(u | w) = N (u |D w m, D w SD w ) (4.8)
By maximizing the variational bounds that follow, we impose a q that performs well under a sparsified inducing set. We continue by applying Jensen's inequality on p (y |u, w ), obtaining:

log p (y |u, w ) ≥ E p(f |u,w) log p (y |f ) (4.9)
We can now substitute (4.9) into (4.7), obtaining a bound where we expand q(u, w) as q (u |w ) p (w). By making this assumption, we obtain the following evidence lower bound L elbo :

N n=1 E p(w) E q(u|w ) E p(fn|u,w ) log p (y n |f n ) - 1 N KL q (u |w ) p (u |w ) (4.10)
Recall that p (w) is implicit: although we do not make any particular assumptions about its analytical form, we can draw samples from it. Using MC sampling from p (w), we can obtain the approximation Lelbo : (4.11) where w(n) is sampled from p (w).

N n=1 E q(u| w(n) ) E p(fn|u, w(n) ) log p (y n |f n ) - 1 N KL q u w(n) p u w(n) ,
4.3.1.0.1 Sampling from the set of inducing points.

Recall that any sample w from p (w) is a binary vector, i.e. w ∈ {0, 1} M . In case all elements of w are set to one, our approach recovers the original svgp with computational cost of O (M 3 ) coming from computing p (f n |u, w = 1)

and KL (q (u |w ) p (u |w )) in the elbo. When a wi is set to zero, the entries of the i-th row and i-th column of the covariance matrix in p (u |w ) and q (u |w ) are zero. This means that the i-th variable becomes unnecessary, so we get rid of i-th row and column in these matrices, and also eliminate the i-th element in mean vectors of q (u |w ) and p (u |w ). This is equivalent to selecting a set of active inducing points in each training iteration.

H-nearest inducing inputs

Despite the fact that p(w) is an implicit distribution, we have been able to define and calculate a variational bound, assuming we can sample from p(w).

We shall now describe our sampling strategy, which relies on neighbor information of random mini-batches.

In order to explain the idea conveniently, we introduce Z H x as the set of Hnearest inducing inputs. Intuitively, the prediction for an unseen data x using

Z H
x is a good approximation of the prediction using all M inducing points, that is Z M x . This can be verified by looking at the predictive mean, which is expressed as a linear combination of kernel functions evaluated between training points and a test point, as in Eq. ( 4.3). The majority of the contribution is given by the inducing points with the largest kernel values, so we can use this as a criterion to establish whether an inducing input is "close" to an input vector (the effect of different kernels on the definition of nearest neighbors is explored in the supplement ELL ← 0 and KL ← 0.

4:

Sample mini-batch I of size n from D.

5:

for (x i , y i ) ∈ I do 6:

Find Z H x i , i.e. the H-nearest Z to x i .

7:

Compute w (x i ) using Z H x i as in (4.12)

8:

Extract m w(x i ) and S w(x i ) from m and S.

9:

Compute q (f i |w (x i )) as in (4.13).

10:

ELL ← ELL + E q(f i |w(x i ) ) log p (y i |f i ).
11:

KL ← KL + KL q u w(x i ) p u w(x i )

12:

end for 13:

Lelbo ← N n ELL -1 n KL.

14:

Update θ using the derivative of Lelbo . 15: end while test time, however, the inputs of interest are not random; we need to describe the predictive distribution in terms of the deterministic function w (x). In fact, if we would like to approximate the predictive distribution at x n using H-nearest inducing inputs to x, i.e. Z H xn , then w (x) = w

(1)

x ...w We extract the relevant elements using w (x); for the mean, we have m w(x i ) = D w(x i ) m, and for the covariance we select the appropriate rows and columns using S w(x i ) = D w(x i ) SD w(x i ) . The approximate posterior over f i given w (x i ), i.e. q (f i |w (x i )) is:

N f i |A x i m w(x i ) , K x i + A x i S w(x i ) -K Z H x i A x i , (4.13) where A x i = K x i ,Z H x i K -1 Z H x i
. One-dimensional regression example. We visualize the posterior distribution for a synthetic dataset generated on a one-dimensional input space.

We execute svgp and swsgp, and depict the posterior distributions of these two methods by showing the predictive means (orange lines) and the 95% credible intervals (shaded areas) in Figure 4.2. We consider identical settings for the two methods (i.e. 128 inducing points, kernel parameters, likelihood variance) and a neighbor area of 16 for swsgp; a full account of the setup can be found in the supplement. We see that although the models are different, the predictive distributions appear remarkably similar. A more extensive evaluation follows in Section 4.4.

Complexity

The computational cost of swsgp is dominated by lines 6, 8 and 9 in Algorithm 2. For each data point (x i , y i ) in mini-batch I, we need to find the H nearest inducing neighbors Z H

x i for n points in line 6, where n = |I|; this contributes to the worst-case complexity by O (nM H).

In line 8, we extract relevant parameters from m and S. We focus on the cost of extracting S w(x i ) from S. Similar to svgp (Section 4.2.1), we consider S = LL T , where L is lower triangular. We extract L w(x i ) = D w(x i ) L which contains the rows of L that correspond to the Cholesky decomposition of S w(x i ) = L w(x i ) L T w(x i ) . The computational complexity of selecting the variational parameters is O (nM H 2 ).

Finally, the computation of approximating the predictive distribution in line 9 requires O (nH 3 ). The overall complexity for swsgp in the general case is O (nM H + nM H 2 + nH 3 ), which is a significant improvement over the O (M 3 ) complexity of standard svgp, assuming that n, H M . If we choose S to be diagonal, the total complexity reduces to O (nM H + nH 3 ); if we additionally consider Z to be fixed, the computational cost is O(nH 3 ). In the experiments of Section 4.4 we also explore these settings.

Experiments

In this section, we conduct experiments to evaluate swsgp on a variety of experimental conditions. We denote our approach by swsgp-M-H, where M inducing points are used and H determines how many neighbors are selected.

We introduce svgp-M, svgp-H and svgp-M-H as competitors; svgp-M and svgp-H are using M and H inducing points, respectively. svgp-M-H, instead, refers to svgp using M inducing points at training time and H-nearest inducing inputs at test time.

The comparison is carried out on some UCI data sets for regression and classification, i.e., powerplant, kin8nm, naval, eeg, credit, and spam. We also consider larger scale data sets, such as mnist and the airline data. We use the Matérn-5 /2 kernel in all cases except for the airline dataset, where the sum of a Matérn-3 /2 and a linear kernel is used, similar to Hensman et al. (2015a). All models are trained using the Adam optimizer [START_REF] Kingma | A method for stochastic optimization[END_REF] with a learning rate of 0.001 and a mini-batch size of 64. The likelihood for regression and binary classification are set to Gaussian and probit function, respectively. All models are trained over 100, 000 iterations except for the airline data set where models are trained for one million iterations. In regression tasks, we report the test root mean squared error (rmse) and the test mean negative log-likelihood (mnll), whereas we report the test error rate (err) and mnll in classification tasks. The results are averaged over three folds.

Increasing the number of neighbors

We begin our experimental evaluation by investigating the behavior of swsgp with respect to H. In Figure 4.

Increasing the number of inducing points

In this set of experiments, we show that the performance swsgp improves when increasing the total number of inducing points, while keeping the number of active inducing points H fixed. We first illustrate this on the banana We also investigate the impact of increasing H and M simultaneously. In each regression and classification data set, we test swsgp with H = 4, 8 and M = 8, 16, 32, 64. The results shown in Fig. 4.6 indicate that using a small H is not detrimental to performance when M is large. In addition, swsgp with a small H is comparable or better than svgp in almost all cases. 4.4.4 Large-scale problems with a huge number of IPs

We showcase a large-scale classfication problem, where we illustrate that swsgp enables the possibility to use sparse gps with a massive number of inducing points without incurring a prohibitive computational cost. We employ the airline data set, featuring 5 million training points. We test swsgp with M = 100, 000 inducing points. We attempted to run svgp with such a large M without success (out of memory in a system with 32GB of RAM).

Therefore, as a baseline we report the results of svgp with the configuration in Hensman et al. (2015a).

In swsgp, we impose a diagonal matrix S in the variational distribution By setting H and the mini-batch size n to 100 and 16 respectively, in about 24 hours of training we could run swsgp-100,000-100 for one million iterations.

q (u | w),
The err and mnll of swsgp-100,000-100 evaluated on the test set are 21% and 0.48, respectively, while the err and mnll of svgp-200 published in Hensman et al. (2015a) are about 34% and 0.61, respectively. To the best of our knowledge, swsgp is the first to enable sparse gps with such a large set of inducing points without imposing a grid structure on the inducing inputs.

We conclude by reporting comparisons with other gp-based models. In particular, we compare against the Stochastic Variational Deep Kernel Learning (svdkl) (Wilson et al., 2016) and the Deep gp approximated with random features [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF] As we discuss in the paper, the selection of H-nearest inducing points Z H x is made by using the kernel as a proxy to the concept of distance. Intuitively, a kernel defines the similarity between two points in the input space, which is more formally expressed as correlation. The kernel implicitly defines a kind of distance that we use to determine the active neiborhood. Thus, the selected neiborhood is dominated by the inducing points with largest kernel values.

In the main paper, we have used different versions of the Matérn kernel. We shall now explore the effect of our neiborhood-selection strategy on a number of different kernels, both stationary and non-stationary. We apply swsgp on the banana data-set using different heuristics for the H-nearest inducing points selection. Let K-swsgp denote what is essentially the vanilla version of our method, where the kernel-based heuristic is used as a proxy to distance. In the case of the rbf kernel, K-swsgp essentially corresponds to the Euclidean distance. We also examine a random-based heuristic (R-swsgp) in which H-nearest inducing points are randomly chosen. In all cases, we set M and H as 32 and 8 respectively. We also compare against svgp with M of 32.

In Fig. Regarding the kernels rbf, arc-cosine-0 and arc-cosine-1, our method (K-swsgp) seems to be virtually identical to svgp. The advantages of Kswsgp over svgp are shown when using polynomial-3. It is highly possible that the flexibility of variational distribution over inducing variables, i.e. q (u), in swsgp is the main reason for this difference.

Further visualizations on 1d examples

We demonstrate swsgp on one-dimensional regression problem. We have generated a synthetic data-set by sampling inputs x i from the interval [-2, 2]; the targets have been computed as y i = sin(12x i ) + 0.66 cos(25x i ) + ε, where ε is additive Gaussian noise with variance 0.1. We also show that the performance of swsgp improves when increasing the total number of inducing points while keeping the number of active inducing points H fixed. We intuitively expect that a larger the total number of inducing points should translate to a more accurate model. In these experiments, the size of neighbor area is fixed to 4, i.e. H = 4, and the total number of inducing points are varies from 4 to 64. We see that the sequence of the predictive means in Fig. 4.9 are more and more accurate from left to right.

Although we are using a small neighbor area, our model is improved when increasing the total number of inducing points.

Conclusions

Sparse approaches that rely on inducing points have met with success in reducing the complexity of gp regression and classification. However, these methods are limited by the number of inducing inputs that is required to obtain an accurate approximation of the true gp model. A large number of inducing inputs is often necessary in cases of very large datasets, which marks the limits of practical applications for most gp-based approaches.

In this work, we further improve the computational gains of sparse gps.

We proposed swsgp, a novel methodology that imposes a hierarchical and sparsity-inducing effect on the prior over the inducing variables. This has been realized as a conditional gp given a random subset of the inducing points, which is defined as the nearest neighbors of random mini-batches of data. We have developed an appropriate variational bound which can be estimated in an unbiased way by means of mini-batches. We have performed an extensive experimental campaign that demonstrated the superior scalability properties of swsgp compared to the state-of-the-art.

Conclusion

The models and techniques presented in this thesis are unified by the overarching goal of improving the calibration and scalability of Gaussian Processes.

We conclude this thesis by summarizing the principal themes and contributions presented in the preceding chapters, with particular emphasis on their significance in the context of complementary work in this direction of research.

This is followed by a brief outlook on possible avenues for future work where we indicate how one might go about achieving these objectives.

Themes and Contributions

In this thesis, we primarily investigated the following themes in relation to Gaussian processes:

• Well-calibrated deep convolutional probabilistic model. Developing models which are able to provide accurate predictions and reliable uncertainties has been a long-standing research topic attracting significant attention from machine learning community. Deep cnns that have accomplished stateof-the-art results in a range of tasks have been illustrated to be miscalibrated, the depth of architecture are the main factor affecting calibration [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF]. Thinking of Bayesian priors as a form of regularization, it is natural to assume that Bayesian cnns are an appropriate treatment for the problem of miscalibration of modern cnns. Independently of the works on Bayesian cnns implemented by Monte Carlo Dropout [START_REF] Gal | Bayesian Convolutional Neural Networks with Bernoulli Approximate Variational Inference[END_REF], there have been other attempts to give a probabilistic flavor to cnns by combining them with Gaussian processes (Wilson et al., 2016;[START_REF] Bradshaw | Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks[END_REF][START_REF] Van Der Wilk | Convolutional Gaussian Processes[END_REF]. To the best of our knowledge, prior to our work there were no studies showing calibration properties of these Bayesian cnns approaches. Hence, in Chapter 3, we investigated the calibration properties of Bayesian treatment on cnns. Perhaps surprisingly, the results indicated that current combinations of cnns and gps are miscalibrated, with a tendency of being overconfident in predictions. Consequently, by extending the random feature expansion approximation for dgps [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF], we proposed a novel combination of cnns and gps which is well-calibrated, and we validated it through several experimental results on image classification tasks.

Furthermore, our model was extended by replacing the last fully-connected layers of cnns with Deep gps [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF] and by employing structured random features to obtain faster and more compact gp approximations [START_REF] Le | Fastfood -Approximating Kernel Expansions in Loglinear Time[END_REF]Yu et al., 2016).

• Combination of global and local approximation. Gaussian Processes [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] offer a powerful statistical framework for inference on functions. However, the applicability of gps on big datasets is hindered by the prohibitive complexity depending on training size N . Due to the rigorous uncertainty quantification of gps, the inducing point-based sparse approximation of gps have been extensively studied [START_REF] Snelson | Sparse Gaussian Processes using Pseudoinputs[END_REF][START_REF] Pratola | A unifying view of sparse approximate Gaussian process regression[END_REF][START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF]Hensman et al., 2015a). The state-of-the-art approaches, e.g Scalable Variational Gaussian Processes (Hensman et al., 2015a), allows for the application of gps to large-scale problems with a small number of inducing points M . As shown recently by [START_REF] Burt | Rates of convergence for sparse variational Gaussian process regression[END_REF], it is possible to obtain an arbitrarily good approximation for a certain class of gp models with M growing more slowly than N . However, in general, it is still possible that the required value for M may exceed a certain computational budget. In Chapter 4, by imposing a sparsity-inducing structure on the prior over the inducing variables and by carrying out a variational formulation of this model, we pushed the limits of scalability and effectiveness of sparse gps enabling a further reduction of computational complexity. Our experimental results showed that the use of unprecedented number of inducing points led to higher accuracy on airline which is a dataset with millions data points. In addition, we showed that our proposed model is able to know what it does not know by yielding sensible predictive uncertainties.

Future work

Beyond the discussion featured in this thesis, the themes explored in this body of work not only motivate immediate extensions for improvements, but also set the foundations for broader long-term objectives. In this section, we expand upon the directions for future work which we believe to be particularly pertinent to ongoing developments in both the theoretical and practical aspects of machine learning using gps. We partition this discussion into the overarching themes of (i) studying calibration properties of gps regression; and (ii) proposing more elegant mixtures of cnns and gps; and (iii) adapting the state-of-the-art scalable gps to online machine learning.

Calibrated GP regression

In addition to improving the scalability of gps, producing reliable predictive uncertainties is also a primary goal for the application of gps in the era of big data, especially when gps are components of larger decision-making pipelines.

This aspect can be evaluated by analyzing calibration properties mentioned in Chapter 3. While the reliability of the predictive uncertainties of Bayesian [START_REF] Platt | Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods[END_REF][START_REF] Guo | On Calibration of Modern Neural Networks[END_REF] or training with adversarial samples (Lakshminarayanan et al., 2017).

Elegant mixtures of CNNs and GPs

While studying on Bayesian cnn, we have realized that combining cnns and gps does not generally improve the performance of standard gps. We speculate that the kernel's parameterization with a high number of parameter increases the risk of overfitting, and leads to overconfident tendency in predictions. As shown in Chapter 3, the Bayesian treatment on convolutional parameters enhances not only model's generalization but also model's calibration. However, the improvements of our approach carries a great computational cost due to repeated feed-forward procedure. This limit can serve as a motivation for investigating new approximation methods for scalable inference in gp models and combinations with cnns.

Generally, the Bayesian flavor in the mixtures of cnns and gps can be strengthened by applying a full Bayesian treatment. For example, following our works in Chapter 3, the proposed models can be further improved by applying a Bayesian treatment on priors of parameters, which would result in the optimization of dropout rates of convolutional hyperparameters [START_REF] Kingma | A method for stochastic optimization[END_REF][START_REF] Molchanov | Variational dropout sparsifies deep neural networks[END_REF]Louizos et al., 2017). Along a similar vein, independent works replacing the fully-connected layers of cnns by gps (Wilson et al., 2016;[START_REF] Bradshaw | Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks[END_REF][START_REF] Tran | Calibrating deep convolutional gaussian processes[END_REF], while Deep Convolutional Gaussian Processes (DCGPs) proposed by [START_REF] Blomqvist | Deep convolutional gaussian processes[END_REF] substitutes gps for convolutional filters. Another interesting approach could be applying the Bayesian formulation mentioned in Titsias and Lazaro-Gredilla (2013) together with sparsity inducing priors (Louizos et al., 2017;[START_REF] Molchanov | Variational dropout sparsifies deep neural networks[END_REF] on DCGPs, a procedure of learning architecture is proposed, which not only accelerates computations but also allows one to approximately integrate out kernel hyperparameters, such as length-scales.

Adaptability to online machine learning

According to the extensive literature review in Liu et al. (2018b), local approximations are common approaches to implement scalable statistical inference systems. The uses of local approximations require to define the localization of experts, which directly affects to the assignments of data points to local experts. Likewise, in chapter 4, our proposal named Sparse-within-sparse Gaussian Processes (swsgp) perceived as a combination of global and local approximations also relies on the way to select active inducing points for each inputs. On offline tasks, swsgp was shown to be effective in terms of accuracy and complexity. With the application of online machine learning wherein training sets are constantly evolving, the selection of active inducing points based upon spatial or temporal distance, which is implemented in swsgp, may ignore the information related to periodic patterns. In such scenario, a kernelbased distance seems to be more appropriate because the kernel intuitively determines the correlation between two points in the input space. More general, by perceiving the selection of active inducing points as a gating function, the input-dependent Dirichlet Process [START_REF] Rasmussen | Infinite mixtures of gaussian process experts[END_REF] and Polya urn distribution (Meeds and Osindero, 2006)

  Source of Predictive Uncertainty. Predictive uncertainty is a conflation of several separate factors: model uncertainty, data uncertainty and distributional uncertainty. Model uncertainty or epistemic uncertainty represents the uncertainty in the estimate of model's parameter given the training data. This uncertainty can be explained away given enough data. Data uncertainty or aleatoric uncertainty comes from the complication in the observations, such as class overlap, label noise, input-dependent noise. As this kind of uncertainty accompanying the nature of data, it is irreducible even if more data are collected. Distributional uncertainty appears due to the mismatch between the training and testing distribution.

  and X * include N training points and N * testing points, respectively, then m X and m X * contain N and N * values of the mean function at X and X * ; and K (X, X * ) denotes the N × N * matrix of the covariances evaluated at all pairs of training and testing points, and similarly for the other covariance matrices.

  .20) Automatic Relevance Determination Kernel. The kernel functions mentioned above are called isotropic where the flexibility of kernel function is indicated by a lengthscale parameter, l. To enhance the flexibility of kernel function, we augment D length-scale parameters, l 1 , . . . , l D accompanying

  .44) Partially Independent (Training) Conditional. Having compared the predictive distribution of dtc and fitc defined in equation 2.37 and 2.42, it is obvious that fitc is an improvement of dtc by remaining the exact diagonal elements of the covariance matrix. Relying on the unifying framework, Quiñonero Candela and Rasmussen (2005) have proposed a further improved approximation compared to fitc by extending the training conditional to have a block of diagonal covariance and remaining the exact testing covariance as defined in equation 2.26. Due to the usage of diagonal block covariance on training conditional, the approximation is called Partially Independent Training Conditionals (pitc) .

  in equation 2.64, each training iteration of sgpr requires the budget of O (N M 2 ) for computation and O (N M ). These costs come from the linear algebraic operation appearing the computation of Q X,X , i.e. the matrix inversion of K -1 Z,Z and the matrix multiplication K X,Z K -1 Z,Z K Z,X . Though the reduction of sgpr on computation and memory requirement are impressive, these demands are quickly prohibitive for big data, where the training size N reaches to many millions or billions. In order to overcome the dependency of complexities on training size, Hensman et al. (2013) have employed Stochastic Variational Inference on Gaussian

aaaa

  11 b 11 a 11 b 12 . . . a 11 b 1q . . . . . . a 1n b 11 a 1n b 12 . . . a 1n b 1q a 11 b 21 a 11 b 22 . . . a 11 b 2q . . . . . . a 1n b 21 a 1n b 22 . . . a 1n b 2q . 11 b p1 a 11 b p2 . . . a 11 b pq . . . . . . a 1n b p1 a 1n b p2 . . . a 1n b pq . m1 b 11 a m1 b 12 . . . a m1 b 1q . . . . . . a mn b 11 a mn b 12 . . . a mn b 1q a m1 b 21 a m1 b 22 . . . a m1 b 2q . . . . . . a mn b 21 a mn b 22 . . . a mn b 2q . m1 b p1 a m1 b p2 . . . a m1 b pq . . . . . . a mn b p1 a mn b p2 . . . a mn b pq

  (2.83) It is infeasible to access the logarithm marginal likelihood of gp regression on X containing N = D d=1 G d points due to the computational bottlenecks from the algebraic operations, i.e. the inversion and matrix-vector multiplication (K X,X + σ 2 n I) -1 y and logarithm of determinant log |K X,X + σ 2 n I|. The original computational and storage cost are O (N 3 ) and O (N 2 ) respectively. Due to the nice properties of Kronecker product, the complexity of learning and inference turns out O DN 1+ 1 D and O DN 2 D for storage. In the next section, I will analyze and explain why Kronecker product can lead to the improvements.

A

  d have the same dimensions. Similarly, the algorithm kron_mvm also works with matrices {A} D d=1 with various sizes, i.e. A d ∈ R G d ×G d . Consider the iterative steps appear in the loop, the computational cost mainly relies on the matrix multiplication A d X which requires O (N G) or O N 1+ 1 D . The loop Algorithm 1 Fast Matrix Vector Multiplication with Kronecker Product -kron_mvm. Input: G-by-G matrices A 1 , . . . , A D , N -dimensional vector b where N

  and we define φ(x) as a column vector of length 2N RF containing the evaluation of the m pairs of trigonometric functions at x. φ (x) = cos x T ω(1) , . . . , cos x T ω(N RF ) , sin x T ω(1) , . . . , sin x T ω(N RF ) T(2.105)

  Inducing point-based and random feature-based approximations of gps are implemented based on a global distillation, and they are commonly used to approximate gps. However, these approaches require the computational and storage costs which are determined by auxiliary variables, i.e. number of inducing points or spectral samples. An alternative class of methods for improving the scalability of gps is to follow the divide-and-conquer idea, which focuses on the local subsets of training data. According to the literature survey conducted by(Liu et al., 2018b), in this text, we opt to split the approach of local approximation into two groups: Separate-Local-Experts and Ensemble-Local-Experts.Separate-Local-Experts. Intuitively speaking, there is almost no dependence between two points which are distant from each other. Thus, the prediction at an unseen input can be made sensibly by using localized experts with an acceptable computational cost. For example,[START_REF] Kim | Analyzing nonstationary spatial data using piecewise gaussian processes[END_REF] andDatta et al. (2016) assume that a local expert model completely governs prediction at inputs inside its corresponding area. Simply, these approaches firstly parti-tion the input space, then all local experts are trained based on these disjoint subsets, and then the inference at x * can be made by an appropriate local expert. By introducing M i as a local expert which is responsible for the subregion Ω i and D i as the subset of data located inside Ω i , we mathematically state that the predictive distribution at x * can be approximated by using a subset of data D i , i.e. p (y * | D, x * ) ≈ p (y i | M i , D i , x * ). The partition on input space can be made by some clustering algorithms, e.g. Voronoi tessellations[START_REF] Kim | Analyzing nonstationary spatial data using piecewise gaussian processes[END_REF], and tree techniques[START_REF] Vasudevan | Gaussian process modeling of large scale terrain[END_REF][START_REF] Pratola | A unifying view of sparse approximate Gaussian process regression[END_REF]. By restricting the number of data points of a local model to M , there are N/M local gps where N is training size. Learning all independent gps experts requires a cost of O (N M 2 ).Instead of grouping data points into disjoint subsets statically before training local gps experts, an alternative approach is to select a neighborhood subsets D * around x * , and train a particular expert M * to make the prediction at x * . For example,[START_REF] Urtasun | Sparse probabilistic regression for activityindependent human pose inference[END_REF] employ a dynamic partition to choose m 0 neighbor points around x * , resulting in O (n t m 3 0 ) complexity that relies on the test size n t . The primary problem of the approach is the concept of the neighborhood set D * around x * . The most straightforward way is to use geometric closeness criteria for selection, i.e. the selected points should be close to x * . However, the approach is not optimal due to these closest points convey redundant information. Thus, there are several gp-based methods which have been employed to sequentially update the neighborhood set (Gramacy, 2016;[START_REF] Gramacy | Speeding up neighborhood search in local gaussian process prediction[END_REF][START_REF] Gramacy | Adaptive design and analysis of supercomputer experiments[END_REF][START_REF] Gramacy | Local gaussian process approximation for large computer experiments[END_REF].While improving significantly the scalability and enjoying the capability of capturing non-stationary features due to the localized structure, Separate-Local-Experts yields discontinuous predictions on the boundaries of subregions, which is illustrated inLiu et al. (2018b). To alleviate the discontinuity problem, the patched gps[START_REF] Park | Efficient computation of gaussian process regression for large spatial data sets by patching local gaussian processes[END_REF][START_REF] Park | Patchwork kriging for large-scale gaussian process regression[END_REF] restricts that two adjacent local GPs are patched to share the nearly identical predictions on the boundary. However, it possibly yields non-sensible predictive variances, and are only available in low dimensional space(Pourhabib et al., 2014b;[START_REF] Park | Patchwork kriging for large-scale gaussian process regression[END_REF]. Another problem of Separable LocalExperts is to suffer from poor generalization since it misses the long-term spatial correlations. To address the generalization issue, we can restrict that all local expert use the same hyperparameters[START_REF] Deisenroth | Distributed gaussian processes[END_REF], or combine local and global approximation of gps as mentioned in[START_REF] Snelson | Local and global sparse Gaussian process approximations[END_REF].Ensemble-Local-Experts. An alternative solution to mitigate the problems raised by Separable Local Experts is to use the model averaging strategy, which is accomplished by an ensemble of local experts. The approach combines various local gps possessing individual hyperparameters for enhancing accuracy and reliability(Yuksel et al., 2012a;[START_REF] Masoudnia | Mixture of experts: A literature survey[END_REF]. Mathematically, Ensemble-Local-Experts can be expressed as a mixture of M Gaussian model, where the weight for each component can be seen as a gating function of covariates, which often takes a parametric form such

  search directions. The first one is to the localization of experts. This can be accomplished by Expectation Maximization (ME) algorithm, wherein the data points are assigned to local experts through Maximum a Posterior in E-step(Nguyen and Bonilla, 2014b; Zhao et al., 2015b;[START_REF] Chen | A precise hard-cut em algorithm for mixtures of gaussian processes[END_REF], and subsequently, the optimization in M-step only operates on small subsets of data. The second one is to combine global approximation with local experts. When using m inducing points for each local gps that is responsible for n samples, the complexity for training M experts is intuitively O (nm 2 M ),

Figure 3 . 1 -

 31 Figure3.1 -Reliability diagrams for three state-of-the-art combinations of cnns and gps, i.e gpdnn[START_REF] Bradshaw | Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks[END_REF]), cgp (van der Wilk et al., 2017), svdkl(Wilson et al., 2016) applied to cifar10 and cifar100 data sets with LeNet and resnet architectures. See table3.1 for details on the convolutional architectures that we apply to cifar10 and cifar100. Because it is not possible to specify the convolutional structure in cgp (van der[START_REF] Van Der Wilk | Convolutional Gaussian Processes[END_REF], the left and central panels show the same curve for cgp.

  are specified by a covariance (kernel) function k (with hyperparameters θ) and this form is shared across output dimensions, although this can be relaxed and allow for a different k for the Q outputs.Instead of applying the gp modeling directly to the images, we propose to employ a transformation c(x|Ψ) using convolutional layers, where Ψ denotes the parameters of such layers. The vector-valued function c(x|Ψ) is differentiable as it implements a series of differentiable operations, such as convolutions and

Figure 3

 3 Figure3.2 -Comparison of our cnn+gp(rf) and cnn+gp(sorf) with existing combinations of cnns with gps, and with Bayesian cnns and postcalibrated cnns. All performace metrics are defined so that the lower the better.

Figure 3

 3 Figure 3.3 -Reliability diagrams of our cnn+gp(rf) in comparison with existing combinations of cnns with gps, and with Bayesian cnns and postcalibrated cnns.

Figure 3

 3 Figure 3.4 -Performance of the proposed model when varying the depth of the dgp on top of a resnet convolutional structure on cifar10 dataset. Note that the scale of y-axes indicates that the metrics change only slightly when increasing the depth of the dgp.

  Figure3.5 -Cumulative distribution function plot of predictive entropies when the models trained on mnist are tested on mnist and not-mnist. We report results for two different depths of the convolutional structure. not-mnist dataset available at http://yaroslavvb.blogspot.fr/2011/ 09/notmnist-dataset.html

Figure 3

 3 Figure 3.6 -Impact of optimization of sorf parameters

  uration for cgp according to the implementation released by the Authors. In each figure, rows correspond with the dataset and convolutional architecture, while the column refer to the training size. After the training phase, all mod-

  els are evaluated on the entire testing set. The number of bins used to draw the reliability diagram is 20. In each subfigure, the dashed line indicates perfect calibration. The horizontal axis is the softmax output ranging from 0 to 1. The vertical axis indicates accuracy rate for the red line or frequency for the green bars. The red dot is the real average accuracy at each bin, while the line segments at the red dots refer to the standard deviation of the accuracies. The green bar is the average frequency histogram at each bin of softmax values. The experiments of gpdnn, cgp, mcd-cifar10-LeNet and cnn+gp(rf) are repeated three times.

  accuracy by considering a larger set of inducing variables. The idea is to operate on a subset of H inducing points during training and prediction, with H M , while maintaining a sparse approximation with M inducing variables. We formalize our strategy by imposing a sparsity-inducing structure on the prior over the inducing variables and by carrying out a variational formulation of this model. This extends the original svgp framework and enables mini-batch-based optimization for the variational objective. We then consider ways to select the set of H inducing points based on neighbor information; at training time, for a given mini-batch, we activate H out of Minducing variables considering the nearest inducing inputs to the samples in the mini-batch, whereas at test time we select inducing variables corresponding to the inducing inputs which are nearest to the test data-points. We name our proposal Sparse within a Sparse gp (swsgp). swsgp is characterized by a number of attractive features: (i) it improves significantly the prediction quality using a small number of neighboring inducing inputs, and (ii) it accelerates the training phase, especially when the total number of inducing points becomes large. We extensively validate these properties on a variety of regression and classification tasks. We also showcase swsgp on a large scale classification problem where we set M = 100, 000; we are not aware of other approaches that can handle such a large set of inducing inputs without imposing some special structure on them (e.g., grid) or without considering one-dimensional inputs.

  The variational parameters m and S, as well as the inputs Z, are optimized by maximizing a lower bound on the marginal likelihood p(y | X) = p(y|f )p(f | X)df . The lower bound on log p (y | X) can be obtained by considering the form of q(f , u)

  on the prior over f Our strategy simply assumes a certain structure on the auxiliary variables, but it has no effect on the prior over f ; the latter remains unchanged. Let I and J bet the sets of indices such that w I = 1 and w J = 0. Given an appropriate ordering, the conditional u | w is effectively the element-wise product [u I , u J ] = u • w. This reduces the variances and covariances of some elements of u to zero yielding a distribution of this form:

Figure 4

 4 Figure 4.1 -The choice of inducing points does not affect the prior samples drawn from p(f ). Left: visualizations of f | u, w for different samples of w. Right: comparison of the marginalised (w.r.t. u, w) prior over f , against the true p(f ).K X,Z I K -1 Z I K Z I ,X. These conditionals can be seen for different samples of u, w in the left side of Figure4.1, while in the right side we compare the marginalized prior over f against the true gp prior.

  z m ∈ Z H x 0 else, with m = 1, ..., M (4.12)

Figure 4 . 2 -

 42 Figure 4.2 -Visualization of posterior distribution of svgp and swsgp. In both cases, we consider 128 inducing points; in terms of our scheme (swsgp) we use 16 neighbors.

Figure 4

 4 Figure 4.4 -Evaluation of swsgp on high-dimensional data sets with increasing H. The black up-triangles are for svgp with M inducing points, the cyan down-triangles are for svgp with H inducing points, the red circles are for svgp training with M inducing points and the prediction at an unseen data x are made by Z Hx , and the green squares are for swsgp. In these experiments, M is set to 64 and H varies from 4 to 32. Horizontal axis shows various configurations of H. The standard deviation of the error metrics over the different folds is represented by vertical bars; they are very small for most configurations.

Figure 4

 4 Figure 4.6 -Evaluation of swsgp on high-dimensional data sets with increasing M . The black up-triangles are for svgp with M inducing points. The green stars and plus are for swsgp with H of 4 and 8 respectively. In these experiments, M varies from 4 to 64, as shown on horizontal axes. The standard deviation of the error metrics over the different folds is represented by vertical bars; they are very small for most configurations.

  and we fix the position of the inducing inputs during training. By fixing the inducing inputs, we can operate with pre-computed information about which inducing inputs are neighbors of training inputs. Thanks to these settings, swsgp's training phase requires O (nH 3 ) operations only, where n is the mini-batch size. Due to the appropriate choice of H and n, and the computational cost being independent of M , unlike svgp, we can successfully run swsgp with M = 100, 000.

4. 5 Figure 4

 54 Figure4.7swsgp on various kernels and strategies for selecting the Hnearest inducing points. As we discuss in the paper, the selection of H-nearest inducing points Z Hx is

  (4.7), we visualize the contours of classifiers of svgp and swsgp with various configurations. Clearly, R-swsgp does not work, i.e. the contours are discontinuous and the locations of contours does not makes sense.

Figure 4

 4 Figure 4.8swsgp is applied on a one-dimensional data set, where M is fixed to 64 and H is increased gradually from 4 to 64. The red dots are inducing positions; the black crosses are testing points; the green line refers to predictive means.

M = 4 ,

 4 Figure 4.9swsgp is applied on 1d. The red dots are inducing positions. The black crosses are testing samples. The green lines are predictive means. The title of each sub-figures shows M and corresponding rmse.

  cnns on classification tasks has been analyzed[START_REF] Guo | On Calibration of Modern Neural Networks[END_REF];Lakshminarayanan et al. (2017);[START_REF] Tran | Calibrating deep convolutional gaussian processes[END_REF], the calibration of gp-based regression methods has not been considered carefully. As mentioned in[START_REF] Kuleshov | Accurate uncertainties for deep learning using calibrated regression[END_REF], the calibration property of regressors is evaluated by their predictive interval. A regressor is stated to be calibrated if p-percent credible intervals contain the true outcomes p-percent of the time. Starting with the novel vision about reliable regressors, investigating calibration properties of gps on regression promises to be interesting. Some potential candidates reinforcing the model's calibration may be inspired by the preceding works, for example post-calibration by Platt scaling

  Stationary covariance function.A stationary covariance function of x i and x j only depends on Euclidean distance of x i and x j , i.e. k (x i , x j ) = k S (r),

	where r = || x i -x j || 2 2 . Thus, it is invariant to translations in the input
	space. This kind of kernel are commonly-used because, intuitively, it is a
	basic similarity assumption that points with inputs x which are close are
	likely to have similar target values y, and thus training points that are near to
	a test point should be informative about the prediction at that point. Next,
	we mention two commonly-used isotropic kernel functions. The covariance
	functions are given in a normalized form where k(0) = 1; we can multiply k
	by a (positive) constant σ 2 f to get any desired process variance.
	• Squared Exponential Covariance Function.

The square exponential function or Radial Basis Function (rbf) kernel has the form:

  art combinations of cnns and gps, i.e gpdnn[START_REF] Bradshaw | Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks[END_REF]), cgp (van der Wilk et al., 2017), svdkl(Wilson et al., 2016) applied to cifar10 and cifar100 data sets with LeNet and resnet architectures. See table3.1 for details on the convolutional architectures that we apply to cifar10 and cifar100. Because it is not possible to specify the convolutional structure in cgp (van der Wilk

Table 3

 3

		Data set	cnn architecture	cnn name
	Shallow mnist	2 Conv Layers + 2 Fully connected	LeNet
	Shallow cifar10	2 Conv Layers + 3 Fully connected	LeNet
	Deep	cifar10	30 Conv Layers + 1 Fully connected resnet
	Deep	cifar100 150 Conv Layers + 1 Fully connected resnet

.1cnn architectures considered in this work. The same architectures are used in gpdnn and svdkl by replacing the fully connected layers with gps, while cgp does not explicitly use a convolutional structure.

  assumes an L 2 regularization which implies a zero-mean Gaussian prior, which is inappropriate for d i as it is Rademacher distributed. We propose to bypass this limitation by applying mcd to a diagonal element is sampled based on the variational parameters M d i -d *

			i
	d i =	M d i -d * i + d * i , with probability π d d * i , otherwise	(3.11)
	Convolutional Networks with Random-Feature-Expanded Deep gps:

reparameterization of d i . In particular, denoting by d * i ∈ {-1, +1} d the initialized values of d i , we apply mcd to d i -d * i . According to this choice, each

Table 3

 3

	blogspot.fr/2011/

.2 -Comparison between cnn+gp(sorf) and mcd with alexnet architecture on cifar10 and cifar100.

Table 3 .

 3 3 -Results on the proposed cnn+gp(sorf) when fixing or learning covariance parameters θ. All results were obtained on mnist, cifar10, and cifar100 without subsampling the data. Please refer to table 1 in the main paper for details on the convolutional structure corresponding to SHALLOW and DEEP.

			SHALLOW
			mnist		cifar10
	Metrics	Fixed	Learned	Fixed	Learned
	err	0.006	0.005	0.203	0.192
	mnll	0.018	0.018	0.610	0.584
	ece	0.002	0.003	0.015	0.010
	brier	0.009	0.008	0.288	0.271
			DEEP
			cifar10	cifar100
	Metrics	Fixed	Learned	Fixed	Learned
	err	0.113	0.115	0.352	0.359
	mnll	0.348	0.355	1.264	1.287
	ece	0.051	0.054	0.050	0.054
	brier	0.170	0.173	0.466	0.478

Table 3 .

 3 4 -Results on the proposed cnn+gp(sorf) vs gpdnn when inferring convolutional parameters using mcd. All results were obtained on mnist, cifar10, and cifar100 without subsampling the data. Please refer to table 1 in the main paper for details on the convolutional structure corresponding to SHALLOW and DEEP.

			SHALLOW
		mnist	cifar10
	Metrics cnn+gp(rf)	gpdnn	cnn+gp(rf)	gpdnn
	err	0.005	0.005	0.172	0.172
	mnll 0.014	0.019	0.535	0.531
	ece	0.004	0.005	0.012	0.012
	brier 0.0071	0.008	0.245	0.244
			DEEP
		cifar10	cifar100
	Metrics cnn+gp(rf)	gpdnn	cnn+gp(rf)	gpdnn
	err	0.111	0.190	0.351	0.820
	mnll 0.344	0.675	1.255	8.606
	ece	0.051	0.036	0.050	0.527
	brier 0.168	0.278	0.466	1.268

  ). With this intuition, p (w) becomes a deterministic function w (x) indicating which inducing inputs are activated. For mini-batch-based training, the value of w remains random, as it depends on the elements x that are selected in the random mini-batch; this materializes the sampling from the implicit distribution p(w). The maximization of the elbo in the setting described is summarized in Algorithm 2 (swsgp). At Algorithm 2 Sparse within sparse gp (swsgp).

Input: D, H, M . Result: The optimum of trainable parameters θ. 1: Initialize θ, i.e. kernel's parameters, Z, m and S. 2: while stopping criteria is False do 3:

  3, we examine swsgp on a two-dimensional Figure4.3 -Visualization of swsgp on banana data sets with increasing H. The total number of inducing points M is fixed to 64, while the size of neighbor area H varies from 4 to 64. The red dots represent the inducing inputs. The orange and blue dots are training points from two different classes. The black lines are the contours of a classifier where the predictive mean is 0.5. This suggests that including neighbor information at prediction time, combined with the use of a larger set of inducing points alone is not enough to obtain competitive performance, and that only thanks to the sparsity-inducing prior over latent variables, this yields improvements.Crucially, the performance obtained by swsgp are comparable with those obtained by svgp-M, while at each iteration only a subset of H out of M inducing points are updated, carrying a significant complextity reduction.

	H = 4	H = 16	H = 64

classification data set (banana), where M is fixed to 64 and H is increased from 4 to 64. In general, these boundaries remain sensible across the whole range of values of H, suggesting that swsgp is able to work and converge well even though H is significantly less than M . We also observe that the contours of the classifier become smoother as H is increasing.

We then test swsgp on other data sets with larger dimensional inputs. In these experiments, H is gradually increased to M . For powerplant, kin8nm, naval, eeg, credit and spam, M is set to 64, and for mnist and airline, M is set to 512. In Fig.

(4

.4), we see that swsgp-M-H consistently outperforms svgp-M-H and svgp-H.

  . In the former, kiss-gp is trained on top of a deep neural network which is optimized during training, and in the latter the layers of a deep gp are approximated as parametric models using random feature expansions. Both competitors feature mini-batch-based learning, so this represents a challenging test for swsgp. The results in Tab. 4.2 show that swsgp is comparable with these competitors. We believe that this is a remarkable result obtained by our shallow swsgp, supporting the conclusions of previous works showing that advances in kernel methods can result in performance which are competitive with deep learning approaches (see, e.g., We set the number of local experts to 64, and we use the same number of inducing points for swsgp (with H either 4 or 8). As the size of powerplant and kin8nm are approximately 7000, we set the number of training points governed by a local expert to 100. For the local gp approaches, we choose 64 locations in the input space using the K-means algorithm, and for each location we choose 100 neighboring points; we then train the corresponding local gp expert. Regarding the testing phase, inductive gps simply rely on the nearest local experts to an unseen point x * . Whereas for transductive gps, we use 100 neighbors of x * and the nearest local expert to make predictions. In table 4.3, we summarize rmse and mnll for all methods; swsgp clearly outperforms the local gp approaches in terms of mnll. Table 4.3 -Comparison with Local gp approximations.

	Method	powerplant	kin8nm
		rmse | mnll	rmse | mnll
	Rudi et al. (2017)). swsgp-64-4	4.27 | 2.41		0.11 | -1.27
	swsgp-64-8	4.24 | 2.40		0.10 | -1.38
	Inductive GPs	9.93 | 38.38		0.13 | -0.40
	Method Transductive GPs	6.17 | 18.78	Data set	rmse mnll 0.09 | -0.65
	swsgp-64-4	powerplant	4.29 2.42
	kiss-gp	powerplant	11.26 5.78
	swsgp-100k-100	airline	0.21	0.48
	svdkl		airline	0.22	0.46
	Deep gp random features	airline	0.21	0.46
	Table 4.2 -Comparison of swsgp, kiss-gp (Wilson and Nickisch, 2015),
	svdkl (Wilson et al., 2016) and Deep gps random features (Cutajar et al.,
	2017)
	4.4.5 Comparison to Local GPs

We finally demonstrate that swsgp behaves differently from other approaches that use local approximations of gps. We consider two well-established approaches of local gps proposed by

[START_REF] Kim | Analyzing nonstationary spatial data using piecewise gaussian processes[END_REF]

and

[START_REF] Urtasun | Sparse probabilistic regression for activityindependent human pose inference[END_REF]

. Following

Liu et al. (2018b)

, we shall refer to these methods as Inductive gps and Transductive gps, respectively. We run all methods on two regression data sets: powerplant and kin8nm.

  can automatically infer which inducing points are necessary from data. Another problem in the scenario of online machine learning is to define a scheme for removing unnecessary inducing points. This can be done simply by eliminating the oldest ones. More elegantly, the frameworks proposed by McIntire et al. (2016); Bijl et al. (2016) could be employed.

			Index
	dgp, 4	2961	svgp, 27
	dtc, 19	2962	svi, 26
	fic, 21	2963	svm, 16
	gps, 3	2964	vssgp, 38
	IPs, 67 kiss-gp, 32	2965	CNNs, 7
	ml, 1	2966	fitc, 20
	pic, 22
	pitc, 21	2967	GPR, 9
	rbf, 13	2968	GPs, 9
	ski, 33 sod, 17		ME, 40
	sor, 18

2969

sgpp, 19

cov(F

In this work, we focus in particular on the order-one arc-cosine kernel (Cho and Saul, 2009) k (1) arc (x i , x j |Ψ, θ) =

[sin(α) + (πα) cos(α)] ,

(3.4)

where θ = (σ, Λ = Diag( 2 1 , . . . , 2 d )) and α is the angle between Λ -1 2 c(x i |Ψ)

and Λ -1 2 c(x j |Ψ).

The arc-cosine covariance has a convenient integral representation that allows for a Monte Carlo approximation, obtaining a low-rank approximation to the covariance matrix involving Rectified Linear Unit (relu) activations (Cho and Saul, 2009) Post calibration, mcd and cnn+gp(rf) (our method) are able to yeild calibrated classification.

Conclusions

Despite the considerable interest in combining cnns with gps, little attention has been devoted to understand the implications in terms of the ability of these models to accurately quantify the level of uncertainty in predictions. This is the first work that highlights the issues of calibration of these models, showing that gps cannot cure the issues of miscalibration in cnns. We have proposed a novel combination of cnns and gps where the resulting model becomes a particular form of a Bayesian cnn for which inference using variational inference is straightforward. However, our results also indicate that combining cnns and gps does not generally improve the performance of standard cnns.

This can serve as a motivation for investigating new approximation methods for scalable inference in gp models and combinations with cnns.

Running time