13

role of the inducing inputs in gp models and their treatment as variational parameters or even hyper-parameters. Given their potential high dimensionality and that the typical number of inducing variables goes beyond hundreds/thousands, we argue that they should be treated simply as model variables and, therefore, having priors and carrying out efficient posterior inference over them is an important-although challenging-problem.

• Chapter 7 will finally present some concluding remarks and analyze some possible future trends.

( )

Power profiling during inference on the test set of CIFAR10 with a l e x n e t The task is repeated 16 consecutive times and profiling is carried out using the nvidia-smi tool. . . . Figure 4.14 Analysis of the parameter efficiency of w h v i for scalable

Gaussian process regression with variational inference. . . . Ablation study on the effect of performing posterior inference on different sets of variables. From s v g p, where the posterior is constrained to be Gaussian and the remaining parameters are point-estimated, to our proposal b s g p,

where we infer a free-form posterior for all Ψ = {u, θ, Z}.

We refer the reader to Table 6.1 for details on the methods (colors are matched). . . . . . . . . . . . . . . . . . . . . . . . . Figure 6.6

Analysis of different choices of objectives when used for optimization and sampling. We refer the reader to Table 6.1 for a description of the methods. . . . . . . . . . . . . . . . . The Deep Learning revolution has arguably reshaped the way we do machine learning, both from a methodological and from an application point of view (LeCun, [START_REF] Lecun | Deep learning[END_REF]. Throughout the last decade, the practical advancements and the theoretical understanding of deep learning models and practices has arguably reached a level of maturity such that it is the preferred choice for any practitioner

seeking simple yet powerful solutions to solve machine learning-related problems.

In few years the community made tremendous advances to develop better and more efficient learning algorithms and architectures. Combined with massive availability of data and implementations that can easily scale from small edge devises to huge data-center of cloud-serviced clusters of graphic processing units (g p us), deep learning quickly has gained a place in our society, from personal assistants [START_REF] Su | Plan-then-Generate: Controlled Data-to-Text Generation via Planning[END_REF][START_REF] Muralidharan | Leveraging User Engagement Signals For Entity Labeling in a Virtual Assistant[END_REF], autonomous driving cars [START_REF] Kendall | Learning to Drive in a Day[END_REF][START_REF] Kiran | Deep Reinforcement Learning for Autonomous Driving: A Survey[END_REF][START_REF] Yurtsever | A Survey of Autonomous Driving: Common Practices and Emerging Technologies[END_REF][START_REF] Grigorescu | A survey of deep learning techniques for autonomous driving[END_REF], automated medical diagnostics [START_REF] Dang | Vessel-CAPTCHA: An efficient learning framework for vessel annotation and segmentation[END_REF][START_REF] Chan | Representing and Denoising Wearable ECG Recordings[END_REF][START_REF] Wu | How to leverage the multimodal EHR data for better medical prediction?[END_REF] and drug synthesis [START_REF] Jumper | Highly accurate protein structure prediction with AlphaFold[END_REF][START_REF] Jiménez-Luna | Drug discovery with explainable artificial intelligence[END_REF]. At the same time, this level of dissemination raises questions on how much we blindly rely on these model's predictions, or how fair they are [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF][START_REF] Du | Fairness in Deep Learning: A Computational Perspective[END_REF]. Taking informed decisions is the key to build reliable and trustworthy systems and it is only possible with a sound modeling of the uncertainty, both in the data and in the model. This is even more critical for applications for which accuracy is not the only important performance metric and for which having highly calibrated confidence in the predictions is a strict system requirement.

1 .1

In an attempt at providing an historical review of neural networks, [START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF] traces back their origins to early XIX century [START_REF] Legendre | Nouvelles méthodes pour la détermination des orbites des cometes[END_REF][START_REF] Gauss | Theoria motus corporum coelestium in sectionibus conicis solem ambientium[END_REF][START_REF] Gauss | Theoria combinationis observationum erroribus minimis obnoxiae (Theory of the combination of observations least subject to error[END_REF], arguing that deep neural networks are in practice variants of linear regression methods. For what we consider nowadays neural networks we need to wait the 1950's, when early success stories involve an attempt at solving the learning problem [START_REF] Hebb | The Organization of Behavior[END_REF] and the formalization of the perceptron [START_REF] Rosenblatt | The perceptron: a probabilistic model for information storage and organization in the brain[END_REF][START_REF] Rosenblatt | Principles of Neurodynamics[END_REF]. Arguably the first important milestone for deep learning is backpropagation (also known as reverse mode differentation), an efficient procedure for error minimization through the combination of (stochastic) gradient descent [START_REF] Hadamard | Mémoire sur le problème d'analyse relatif à l'équilibre des plaques élastiques encastrées[END_REF][START_REF] Robbins | A Stochastic Approximation Method[END_REF] and the chain rule of differentiation [START_REF] Leibniz | Memoir using the chain rule[END_REF]. This was preliminarily discussed by [START_REF] Linnainmaa | The representation of the cumulative rounding error of an algorithm as a Taylor expansion of the local rounding errors[END_REF] but applied to learning of neural networks later by [START_REF] Rumelhart | Learning Internal Representations by Error Propagation[END_REF]. This simple method was extended in the following years to speed-up convergence [START_REF] Vogl | Accelerating the Convergence of the Back-Propagation Method[END_REF][START_REF] Battiti | Accelerated Backpropagation Learning: Two Optimization Methods[END_REF], by e. g. using second order derivative information through the diagonalization of the Hessian [START_REF] Becker | Improving the Convergence of Back-Propagation Learning with Second-Order Methods[END_REF]; P. [START_REF] Simard | Local Computation of the Second Derivative Information in a Multi-Layer Network[END_REF], by adapting the global step size (LeCun, P. Y. [START_REF] Simard | Local Computation of the Second Derivative Information in a Multi-Layer Network[END_REF][START_REF] Yu | Dynamic Learning Rate Optimization of the Backpropagation Algorithm[END_REF] or by computing individual learning rates for each parameter [START_REF] Jacobs | Increased rates of convergence through learning rate adaptation[END_REF][START_REF] Silva | Speeding Up Back-Propagation[END_REF]. In more recent years, these methods have become the backbone of optimizers like a dag r a d [START_REF] Duchi | Adaptive Subgradient Methods for Online Learning and Stochastic Optimization[END_REF], a da d e lta [START_REF] Zeiler | ADADELTA: An Adaptive Learning Rate Method[END_REF]), r m s p ro p [START_REF] Schaul | No more pesky learning rates[END_REF] and a da m [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF]. Having now the tool for learning (i. e.optimization), research moved to study new architectures and exploit new computational devices, like g p us. Convolutional neural networks (c n ns), for example, were introduced for improving the performance of deep networks on highly structured data (like images for classification tasks) [START_REF] Lecun | Backpropagation Applied to Handwritten Zip Code Recognition[END_REF][START_REF] Lecun | Handwritten digit recognition with a back-propagation network[END_REF][START_REF] Lecun | Gradient-Based Learning Applied to Document Recognition[END_REF] while long short-term memory (l s t m) networks [START_REF] Hochreiter | Long Short-Term Memory[END_REF], a more powerful class of recurrent neural networks (r n ns), were found to learn long term dependences in the data, which was previously impossible to do [START_REF] Pérez-Ortiz | Kalman filters improve LSTM network performance in problems unsolvable by traditional recurrent nets[END_REF]. By early 2000's after initial implementations of neural networks on g p us were proved successful with speed-up factors up to 20 times [START_REF] Chellapilla | High Performance Convolutional Neural Networks for Document Processing[END_REF][START_REF] Oh | GPU implementation of neural networks[END_REF], it was clear that the next milestone would come by solving the computational issues of these models and by providing open source implementations for heterogeneous computing devices. After effectively solving the classification problem of handwritten digits on mnist [START_REF] Ranzato | Efficient Learning of Sparse Representations with an Energy-Based Model[END_REF][START_REF] Ciresan | Multi-column deep neural networks for image classification[END_REF], [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] proposed to use a multi-gpu implementation of a deep cnn (later known as alexnet) and with this architecture they achieved the best result on imagenet, a popular computer vision benchmark. This was the first time that a deep learning approach was able to outperform manual feature extraction and classic computer vision algorithms.

Since then, accuracy on this benchmark is continuously increasing, as well as the size of the models implemented (see Figure 1.1). The availability of libraries with automatic differentation [START_REF] Baydin | Automatic Differentiation in Machine Learning: A Survey[END_REF] and highly efficient implementations, like Tensorflow (Abadi et al., 2015), PyTorch [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF] and JAX [START_REF] Bradbury | JAX: Composable Transformations of Python+NumPy Programs[END_REF], made deep learning the preferred choice for any practitioner seeking simple yet powerful solutions to solve machine learning related problems.

Examples of such wide adoption of deep learning come from cosmology [START_REF] Dieleman | Rotation-invariant convolutional neural networks for galaxy morphology prediction[END_REF][START_REF] Lanusse | Hybrid physical-deep learning model for astronomical inverse problems[END_REF][START_REF] Fussell | Forging new worlds: High-resolution synthetic galaxies with chained generative adversarial networks[END_REF], experimental physics [START_REF] Acciarri | Convolutional neural networks applied to neutrino events in a liquid argon time projection chamber[END_REF][START_REF] Aurisano | A convolutional neural network neutrino event classifier[END_REF] and neuroscience (J. X. [START_REF] Wang | Prefrontal cortex as a meta-reinforcement learning system[END_REF][START_REF] Yamins | Using goal-driven deep learning models to understand sensory cortex[END_REF][START_REF] Bellec | Long short-term memory and learning-to-learn in networks of spiking neurons[END_REF][START_REF] Zador | A critique of pure learning and what artificial neural networks can learn from animal brains[END_REF], but it has cross-fertilized other computer science fields, such as digital hardware design [START_REF] Kwon | Transfer learning for design-space exploration with high-level synthesis[END_REF][START_REF] Zhang | From homogeneous to heterogeneous: Leveraging deep learning based power analysis across devices[END_REF][START_REF] Lu | TP-GNN: A graph neural network framework for tier partitioning in monolithic 3D ICs[END_REF], data management systems [START_REF] Kraska | The Case for Learned Index Structures[END_REF][START_REF] Li | QTune: A Query-Aware Database Tuning System with Deep Reinforcement Learning[END_REF][START_REF] Liu | Cardinality Estimation Using Neural Networks[END_REF][START_REF] Marcus | Plan-Structured Deep Neural Network Models for Query Performance Prediction[END_REF] and materials science [START_REF] Zhang | Physics-informed deep learning for digital materials[END_REF][START_REF] Kim | Designing an Adhesive Pillar Shape with Deep Learning-Based Optimization[END_REF][START_REF] Zhang | Finite-Element-Based Deep-Learning Model for Deformation Behavior of Digital Materials[END_REF][START_REF] Wilt | Accelerating Auxetic Metamaterial Design with Deep Learning[END_REF].

.2

Deep learning models, like any other machine learning tools, are generally only a part of larger systems with agents interacting in complex environments. Take as a simple example the problem of autonomous driving (Figure 1.2): the system collects data from multiple sensors (video, Lidar, etc.), it controls the car by e.g. accelerating, braking or turning the steering wheel and it interacts with a chaotic environment made of other cars, obstacles and pedestrians. The output of a model can be used for taking actions, like performing an emergency breaking if another car suddenly changes lane. Generally, decision making systems combine beliefs and data to optimize an utility function. For example, in Figure 1.2 the agent is learning to turn left while minimizing jointly the risk of collision, the time required to complete the manoeuvre and the lateral acceleration and jerk. Unfortunately, real systems are hindered by uncertainty, whether due to partial observability, non-determinism, or a combination of the two. This means that for example we may never know

for certain what state the system would be in after taking a sequence of actions, as predicted by our model. As a consequence, utility theory needs to take into account the unidentifiability of the problem by means of probability theory, which can express the degree of belief an agent has of the model and the environment. Going back to the example above, it can happen that the other car is actually not changing lane but a different shadow is casted on the road in such a way that it looks like.

Quantification of the uncertainty becomes then essential for taking informed decisions. Deep learning models, while being very flexible and powerful models, they are usually not equipped with the ability to accurately estimate the uncertainty of their predictions. Learning algorithms for these models were typically designed through the lens of statistical learning theory by means of minimizing the expected risk defined as

R(f) = e(f(x), y)p(x, y)dxdy (1.1)
where f(•) defines a neural network function in the space of hypothesis F and e(•, •) is a properly defined error function. Due to the intractability of this formulation, in practice we rely on the empirical risk as proxy to the expected risk and it formulates the learning problem as

min f∈F R(f) = min f∈F 1 n n i=1 e(f(x i ), y i ) (1.2)
where x i and y i are the training data. Stability analysis of the learning algorithm with respect to changes in the training data suggest that in order to build models with good generalization properties, we need to impose regularization which controls the tradeoff between data fitting and model variance. For example, a generic regularization algorithm (e.g. Tikhonov regularizer) minimizes the empirical risk plus a stabilizer on the model parameters,

f n = arg min f∈F R(f) + λ f F (1.3)
By building a learning algorithm is this way we obtain a solution f n , which under some conditions, exhibit a bounded generalization error:

P R(f n ) -R(f n ) ε P sup f R(f) -R(f) ε (1.4)
Note that in practice this way of designing learning algorithms implicitly aims at finding the best solution, which equivalently minimizes some expected generalization error. Regression models return a single vector of predictions, while for classification the output represents a probability vector for the associated set of classes.

One way to reason about the confidence in the model's predictions is by treating the problem probabilistically. In practice, instead of having predictions which are point-estimated, we can rely on proper probability distributions around the model outputs. In this way, we can reason about the confidence in the predictions of the model and quantify the risk associated for a down-stream decision making process.

Bayesian machine learning offers a powerful and intuitive way to propagate belief using data. Figure 1.3 visually represents two machine learning tasks (regression and binary classification) which are learned with regularized loss (on the left) and probabilistically with Bayesian inference (on the right). Note how the additional information about prediction confidence at a certain test point can be used in downstream tasks to build more reliable decision making systems.

. 2 .Not all the uncertainties are the same

Before diving into more details of Bayesian inference and how it can be implemented in deep learning models, we need to briefly discuss the role of uncertainty in practice. Broadly speaking, the sources of uncertainty we can distinguish are twofold: we are either uncertain of the data or we are uncertain of the model and its parameters. The former, also known as aleatoric uncertainty, can be due to several factors, including for example noisy data acquisition systems. This uncertainty is irreducible, in the sense that in general we have no controls over it and it's independent on the amount of data available. The second kind of uncertainty is related to the model parameters and it is usually referred to as epistemic uncertainty. This uncertainty is due to the fact that there are up to an infinite amount of possible models that are able to explain a given dataset. As a consequence, we could be unsure about which model parameters we need to choose to predict with. Differently from the aleatoric uncertainty, the epistemic is in fact reducible: once we collect enough data to estimate accurately the model's parameters such uncertainty decreases (even down to zero, in the infinite data limit). On the other hand, by tautology for regions away from the data the model has no information and it will exhibit higher uncertainty. Note that this could be useful to detect out of distribution test data, i.e. situations for which the model is inferring on input points that are far away from training data.

In Figure 1.4, we have a simple example where the aleatoric and the epistemic uncertainties are decomposed and visualized. Note, for example, how the epistemic uncertainty behaves far away from the training data, where due to missing information the variance is larger than in the region with data.

.3

In this section, we briefly review some relevant literature for Bayesian neural networks, from their inception to more recent developments. First examples of Bayesian neural networks can be traces back to late '80s. Researchers at the AT&T Bell Laboratories involved in the general problem of learning neural networks from examples discussed the possibility of using Bayesian inference for achieving this goal. J. [START_REF] Denker | Large automatic learning, rule extraction, and generalization[END_REF] proposes to place an uniform prior distribution on the space of the weights and it also shows how this translates into the function space, effectively by building a network for each configuration of the weights.

Few years later [START_REF] Tishby | Consistent inference of probabilities in layered networks: predictions and generalizations[END_REF] derives a theoretical framework for applying Bayes rule to layered models by assigning a likelihood, a nonsingular prior on the "configuration space" (i.e. the weight space) and by deriving the posterior as a Gibbs canonical distributions on the ensemble of networks. As Yarin Gal pointed out in his thesis [START_REF] Gal | Uncertainty in Deep Learning[END_REF], this is possibly the first example of modern Bayesian neural network (bnn) formalized by means of Bayes rule. The framework of [START_REF] Tishby | Consistent inference of probabilities in layered networks: predictions and generalizations[END_REF] was put in practice by J. S. [START_REF] Denker | Transforming Neural-Net Output Levels to Probability Distributions[END_REF], who propose to use the recently developed back-propagation algorithm for the Laplace approximation of the posterior. D. J. C. [START_REF] Mackay | Bayesian Model Comparison and Backprop Nets[END_REF] discusses the problem of Bayesian model selection and the effect of the Occam factor for the choosing the parameters of a neural network. Using the Laplace approximation he also shows how the model evidence correlates with the test error and generalization. D. J. C. MacKay (1992) extends some of these results, showing how the model evidence can be used to select the number of hidden units. Approximations for Bayesian neural networks were not limited to the Laplace method only. In the same period, Neal (1992a) and [START_REF] Neal | Bayesian training of backpropagation networks by the hybrid Monte Carlo method[END_REF] (1987). In a follow-up work, [START_REF] Neal | Priors for infinite networks[END_REF] analyzes the infinite limit behaviour of Bayesian neural networks priors and he draws important connections with another Bayesian model, the Gaussian process (g p). All these results are collected in his thesis [START_REF] Neal | Bayesian Learning for Neural Networks[END_REF]Neal, 1994a). Variational inference (v i) has been also proposed as a scalable approximation for these models. [START_REF] Barber | Ensemble learning in Bayesian neural networks[END_REF] derives an objective which forms a lower bound of the marginal likelihood. They proposes to use full covariance matrices (differently from what previously done by G. E. [START_REF] Hinton | Keeping the Neural Networks Simple by Minimizing the Description Length of the Weights[END_REF] and a Gamma hyper-prior on the model parameters and they derive a procedure to perform tractable v i with free-form variational posteriors.

. 3 .Modern inference for Bayesian neural networks

More recently, [START_REF] Graves | Practical Variational Inference for Neural Networks[END_REF] proposes a practical way to perform variational inference in a scalable way on neural networks. He discusses several important aspects of the vi framework, like the role of the posterior approximation and how it relates to regularized loss training, how the prior parameters can be adapted during training and finally how Monte-Carlo sampling can be used to approximated the intractable expectation of the log-likelihood. This work was later extended by [START_REF] Blundell | Weight Uncertainty in Neural Network[END_REF] with the so-called Bayes by backprop, which shows how it is possible to use the reparameterization trick to generate unbiased estimates of the variational parameters' gradients. Alternatives can be found in the probabilistic backpropagation (Hernandez-Lobato and R. [START_REF] Hernandez-Lobato | Probabilistic Backpropagation for Scalable Learning of Bayesian Neural Networks[END_REF], in Stein variational inference (Q. [START_REF] Liu | Stein Variational Gradient Descent: A General Purpose Bayesian Inference Algorithm[END_REF] and in the Monte Carlo dropout (Gal and Ghahramani, 2016b). The latter found its fortune [START_REF] Kendall | What Uncertainties Do We Need in Bayesian Deep Learning for Computer Vision?[END_REF] in the embarrassingly simple implementation required to compute the output multiple times with different dropout masks, with nonetheless drawing important theoretical connections with variational inference (Y. [START_REF] Li | Dropout Inference in Bayesian Neural Networks with Alphadivergences[END_REF].

Research in v i for b n n has been also driven by the necessity of moving away from the simple fully factorized Gaussian assumption for the posterior. In this line of works, we can find attempts at parametrizing efficiently matrix variate posteriors [START_REF] Louizos | Structured and Efficient Variational Deep Learning with Matrix Gaussian Posteriors[END_REF][START_REF] Sun | Learning Structured Weight Uncertainty in Bayesian Neural Networks[END_REF][START_REF] Zhang | Noisy Natural Gradient as Variational Inference[END_REF][START_REF] Rossi | Walsh-Hadamard Variational Inference for Bayesian Deep Learning[END_REF]. Note that only around this time we start to see more complex neural network architectures (e.g. c n ns) being employed in the variational setting. [START_REF] Osawa | Practical Deep Learning with Bayesian Principles[END_REF]-building on top of the work of M. E. [START_REF] Khan | Fast and Scalable Bayesian Deep Learning by Weight-Perturbation in Adam[END_REF]-shows how variational inference with natural gradients can scale to ResNet-18 trained on ImageNet. Changing the perspective of v i, Sun, G. [START_REF] Sun | Functional Variational Bayesian Neural Networks[END_REF] proposes to perform variational inference in function space, by defining and optimizing the lower bound directly on the stochastic processes (functional output). Due to the specific kind of priors used, Sun, G. [START_REF] Sun | Functional Variational Bayesian Neural Networks[END_REF] only limit their analysis to regression tasks, leaving classification and other problems as a future work. In a recent work by [START_REF] Farquhar | Liberty or Depth: Deep Bayesian Neural Nets Do Not Need Complex Weight Posterior Approximations[END_REF], the Authors challenge the assumption that structured posteriors are better than the mean-field approximation. In essence, among other things, they speculate that for a given b n n with full covariance there exists a deeper b n n which can induce the same functional posterior in function-space.

As an alternative to approximating the posterior with variational inference, it is also possible to use other kind of approximations. With the Laplace approximation (la), for example, one would choose to fit a Gaussian distribution around the (local) maximum of the joint likelihood (D. J. C. [START_REF] Mackay | Choice of Basis for Laplace Approximation[END_REF][START_REF] Ritter | A Scalable Laplace Approximation for Neural Networks[END_REF]Immer, Bauer, et al., 2021). Recently this approximation has found new life also thanks to efficient implementations for computing the Hessian of neural networks [START_REF] Dangel | BackPACK: Packing more into Backprop[END_REF]E. A. Daxberger et al., 2021). [START_REF] Immer | Improving predictions of Bayesian neural nets via local linearization[END_REF], for example, shows how we can build a locally linearized model starting from the Laplace approximation, while [START_REF] Kristiadi | Being Bayesian, Even Just a Bit, Fixes Overconfidence in ReLU Networks[END_REF] discusses how the la at the last layer is enough to fix the issue of overconfidence in neural networks.

Instead of approximating the posterior, we could also generate samples from this intractable distribution using some gradient-based m c m c. The main limitation of such approaches reside in the expensive computations required to evaluate the gradients of the logarithm of the unnormalized posterior density in big data regimes [START_REF] Neal | MCMC Using Hamiltonian Dynamics[END_REF]. Unfortunately, as [START_REF] Betancourt | The Fundamental Incompatibility of Scalable Hamiltonian Monte Carlo and Naive Data Subsampling[END_REF] discusses and quantifies, mini-batching the data fundamentally breaks and compromises the scalability and accuracy of näive hmc. To solve this problem, T. [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF] proposes to add an additional friction term to counter-balance the effect of the noise introduced by subsampling the data. Notably, T. [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF] draws significant connections between the proposed Stochastic Gradient Hamiltonian Monte Carlo (s g h m c)

and the classic stochastic gradient descent (s g d), which has been leveraged by [START_REF] Springenberg | Bayesian Optimization with Robust Bayesian Neural Networks[END_REF] to estimate some of the hyper-parameters hmc requires.

To improve the performance and the exploration efficiency of s g h m c, R. in theory to converge to the exact posterior (Franzese et al., 2021b) and in practice to deliver high quality samples compared with classic h m c [START_REF] Izmailov | What Are Bayesian Neural Network Posteriors Really Like?[END_REF].

Finally, ensembles of models could be employed as Bayesian approximation [START_REF] Newton | Approximate Bayesian inference with the weighted likelihood bootstrap[END_REF][START_REF] Garipov | Loss Surfaces, Mode Connectivity, and Fast Ensembling of DNNs[END_REF]. [START_REF] Lakshminarayanan | Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles[END_REF] proposes deep ensembles as a simple way to train a set of independent models and combine their predictions to generate informative uncertainty estimation. In a more recent work, B. [START_REF] He | Bayesian Deep Ensembles via the Neural Tangent Kernel[END_REF] generalize this setup to infinitely wide neural networks using the neural tangent kernel (ntk). As an attempt to theoretically characterize the learning behavior of these methods, [START_REF] Milios | A Variational View on Bootstrap Ensembles as Bayesian Inference[END_REF] analyzes the problem, finding that training these models is equivalent to a Kullback-Leibler (kl) minimization akin to variational inference.

.4 '

Despite the fact that we have witnessed huge advancements in the last few years, arguably Bayesian deep learning is still in its infancy, especially from a practical point of view. Broadly speaking, we can summarize the current challenges in two categories, which are tightly entangled: inference and scalability. Bayesian inference, while offering an elegant formulation, cannot be solved exactly. The approximations required could decrease the quality of the predictions while certainly increasing the complexity of the entire machine learning pipeline. Variational inference, for example, increases the model parameters to be optimized by at least a factor of two, while sampling methods typically require orders of magnitude more gradients evaluations. Also, they require to store the samples of the parameters' posterior, which increases the memory footprint of these algorithms hundredfold. On top of this, we are implicitly augmenting the challenges by adding the problem of choosing the prior distribution. This fundamental building block of the Bayesian theory is proved empirically to be quite critical and yet difficult to assess its quality a-priori. On the other hand, we can fortunately find comfort from analyzing other Bayesian models, like g ps and deep Gaussian processs (d g ps), which are closely related to Bayesian neural networks.

.4 .Structure of this thesis

In this thesis, we will attempt at analyzing some of the aforementioned limitations and proposing some solutions. The rest of the thesis is structured as follows:

• Chapter 2 is dedicated to a quick review of Bayesian machine learning and to an introduction to approximations to intractable inference for Bayesian neural networks.

• Chapter 3 will discuss the challenges of initializing variational inference for deep neural networks. We found that this problem is more severe for complex architectures, such as deep convolutional neural networks, for which vi systematically converges to trivial solutions (posterior equal to the prior).

We also propose a novel initialization strategy for v i, whose performance will be assessed in various experiments.

• Chapter 4 proposes a new structured parameterization for variational inference, inspired by the wide literature on scalable kernel machines. The key operation within our proposal is the Walsh-Hadamard transform, and this is why we name our proposal Walsh-Hadamard variational inference (whvi). Unlike mean field vi, whvi induces a matrix-variate distribution to approximate the posterior over the weights, thus increasing flexibility at a log-linear cost in D. We derive expressions for the reparameterization and the local reparameterization tricks, showing that, the computational complexity

is reduced from O(D 2 ) to O(D log D).
• Chapter 5 shows how important is the role of the prior for improving the performance of Bayesian deep learning models and it will address the challenges to accurately selecting them. In more details, we seek to tune the prior distributions over b n ns parameters so that the induced functional priors exhibit interpretable properties, similar to shallow g ps. While b n n priors induce a regularization effect that penalizes large values for the network weights, a g p-adjusted prior induces regularization directly on the space of functions. We will show how we can use the Wasserstein distance between the distribution of bnn functions induced by a prior over their parameters, and a target g p prior.

• Chapter 6 will revisit model approximations and inference schemes for scalable Gaussian processes and deep Gaussian processes. We will discuss the Given the number of times in which an unknown event has happened and failed. [It is required to calculate] the chance that the probability of its happening in a single trial lies somewhere between any two degrees of probability that can be named.

P R O B A B I L I S T I C M E T H O D S F O R M A C H I N E L E A R N I N G
-Thomas [START_REF] Bayes | An essay towards solving a problem in the doctrine of chances[END_REF] In his notes, Bayes derives several theorems that will constitute the backbone of modern probability theory, with particular emphasis on conditional probabilities. He then proposes an argument for using a uniform distribution for the binomial parameter p(0 π 1) to represent absence of knowledge about it [START_REF] Edwards | Commentary on the Arguments of Thomas Bayes[END_REF]. In practice, this is the formulation of what we refer to as Bayes' rule or Bayes' theorem, as a way to solve the inverse conditional probability problem. This new change of prospective to modeling reality is arguably one of the most important milestone in history of Science and Mathematics. Among others, Sir Harold Jeffreys prizes the accomplishments of Bayes by recognizing Bayes' theorem to be "to the theory of probability what the Pythagorean theorem is to geometry" [START_REF] Jeffreys | Theory of Probability[END_REF].

This chapter builds the foundations upon which the rest of the thesis will lie.

Specifically, we start in § 2.1 with a review of basis concept of probability theory and the use of Bayes' rule in machine learning. To this goal, we present some important results of Bayesian inference approximations, notably variational inference (vi) and Markov chain Monte Carlo (mcmc) sampling.

.1

Bayes offers us a principled machinery to do probabilistic modeling: the way of reasoning about uncertainty, a structured procedure to include prior knowledge of the problem (or the lack of thereof) and a formal methodology to model selection are only some of the characteristics that makes Bayes' theorem and Bayesian inference very well suited for machine learning [START_REF] Ghahramani | Bayesian non-parametrics and the probabilistic approach to modelling[END_REF][START_REF] Ghahramani | Probabilistic Machine Learning and Artificial Intelligence[END_REF].

Let's dive a bit more on the theory of Bayesian inference. Assume a generic

In case of supervised learning, we assume both X and y to be known. The Bayesian framework can also be applied to unsupervised tasks and generative modeling.

parametric model f parameterized by some unknown parameters θ (i.e. f(•, θ))

and a collection of data y ∈ R N corresponding to some input points (2.1)

X = {x i | x i ∈ R D in } i=1
In order to keep the notation uncluttered, in the following chapters we might avoid to explicitly write the dependency on X; nonetheless, f is always computed at some input points x. Note also that everything remains valid even if the output dimensionality is greater then one. Depending on the context, y could be a vector or a matrix

.1 .1 Bayesian model selection

The denominator of Bayes' theorem is of particular interest. We refer to this as the model evidence or marginal likelihood, and its role can be better understood by explicitly conditioning the entire inference to a specific model hypothesis M, p(y i | y 1:i-1 ) = p(y 1 )p(y 2 | y 1 )p(y 3 | y 1:2 ) . . . p(y N | y 1:N-1 ) , (2.3)

p(θ | y, M) = p(y | f(X, θ), M)p(θ | M) p(y | M) , ( 2
This highlights how the prediction for the next point y i is conditioned to the previous ones y 1:i-1 , similarly to the leave-one-out cross-validation estimate of the likelihood [START_REF] Murphy | Machine learning -A Probabilistic Perspective[END_REF]. Intuitively, if a model it too complex it will overfit with few examples and will perform poorly with the next ones.

We can further apply Bayes' rule to condition models on the data,

p(M | y) = p(y | M)p(M) p(y, M) dM (2.4)
where the right hand-side should be used in case of discrete model space. Solving such inference problem with complex systems (e.g. climate models, biological cell models, epidemiology, etc.) is challenging and it's a topic of an active area of research. The common and simple approach to this problem of model selection is the maximization of the marginal likelihood (also known as Type-II maximum likelihood).

This can be further seen a maximum-a-posteriori (m a p) estimate of the posterior in Equation (2.4), when we place an uniform prior on the models p(M) ∝ 1. In practice, this is what happens when we optimize the kernel parameters for Gaussian processs (g ps) by maximization of the marginal likelihood.

The marginal likelihood has also an interpretation as generalization metric. Looking at the second axiom of probability theory, for a particular model hypothesis the marginal will sum to one over all possible datasets, e. g. y p(y | M). gives an illustration of this property. For sake of clarity, suppose the dataset y to be a continuous random variable. Here, we plot the marginal likelihood as a function of the dataset for three different models with increasing complexity and we also indicate the observed data y . M 1 is a simple model that performs well only on a very specific data that might not match the one that we observe. The low marginal corresponds to the low probability that indeed M 1 has generated y . On the contrary, M 3 is a very complex model that can predict many datasets; this means M 3 assigns its probability thinly and widely across all possible datamaking also this one unlikely to have generated y . Finally, M 2 is just right as it performs well on y with a reasonable confidence on its neighbor [START_REF] Murphy | Machine learning -A Probabilistic Perspective[END_REF][START_REF] Bishop | Pattern recognition and machine learning[END_REF].

.2

In the previous discussion, we kept the form of the model f(•) generic. In this section, we shall formalize the choice of neural networks, highlighting some important characteristics of this class of models and remarking the challenges for a Bayesian treatment. L) . . . The basic idea behind deep neural networks is to use composition of functions to learn hierarchical features representation of data. The ability of learning and extracting useful features can be exploited for down-stream tasks, like regression or classification.

W (1) W (2) W (
f (2) f (L) y . . .
Deep feedforward neural networks, also known as multilayer perceptrons (mlps), are one of the most simple and popular architectures in deep learning. In feedforward neural networks (see an example in § 2.2), the data is processed by a chain of transformations, without the presence of any feedback loops. Convolutional neural networks are another example of feedforward models, where, differently from mlps, weights are usually shared among different spatial patches of the input images. These architectures have proven to be very efficient in learning useful feature representations for highly structured data, like images, videos and text. Alternatively, architectures that implement loops are also known as deep recurrent neural networks and, despite their popularity for tasks like natural language processing (nlp), will not be discussed in this thesis.

20 -2 0 2 -2 0 2 f 1 f 2 1 hidden units -2 0 2 f 1 10 hidden units -2 0 2 f 1 100 hidden units -2 0 2 f 1
In nite limit GP Figure 2.3: Analysis of the distribution of samples from the prior for a shallow neural network with tanh activation. Priors on the parameters are assumed to be Gaussian with zero mean and standard deviations of 5 for the weights and biases of the rst layer, D -0.5 and 0.1 for the weights and the bias on the second layer. The two random variables f 1 and f 2 are corresponding to inputs x 1 = -0.2 and x 2 = 0.4, respectively. The plots show 500 samples, while the density in the background is estimated with KDE and 10 000 samples.

. .1 Bayesian Neural Networks: Parameterization, Prior and Inference

Let's start by considering a deep neural network (d n n) consisting of L layers, where the output of the l-th layer f l (x) is a function of the previous layer outputs f l-1 (x), as follows:

f l (x) = 1 D l-1 W l a(f l-1 (x)) + b l , l ∈ {1, ..., L}, (2.5)
where a(•) is a nonlinearity, b l ∈ R D l is a vector containing the bias parameters for layer l, and W l ∈ R D l ×D l-1 is the corresponding matrix of weights. In order to keep the notation uncluttered, we refer to the set of weights and bias parameters of a layer l as θ l = {W l , b l }, while the entirety of trainable network parameters is θ = {θ l } L l=1 . In order to simplify the presentation, we focus on fully-connected dnns; note that the weight and bias parameters of convolutional neural networks (c n ns) can be treated in a similar way. To ensure the asymptotic variance of the output neither explodes nor vanishes, the output can be normalized with D l-1 . For fully-connected layers, D l-1 is the dimension of the input, while for convolutional layers D l-1 is replaced with the filter size multiplied by the number of input channels.

The Bayesian treatment of neural networks dictates that a prior distribution p(θ) is placed over the parameters. Chapter 5 will be completely dedicated to the discussion on the role of prior distributions on Bayesian neural network (b n n). Nonetheless, it is important to comment on the asymptotic behavior of neural networks when we consider an infinite width limit. [START_REF] Neal | Bayesian Learning for Neural Networks[END_REF] observed that Recently this was also proved without the use of the central limit theory, by analyzing the effect of width in the prior marginal covariance [START_REF] Pleiss | The Limitations of Large Width in Neural Networks: A Deep Gaussian Process Perspective[END_REF].

"priors over network parameters can be defined in such a way that the corresponding priors over functions computed by the network reach reasonable limits as the number of hidden units goes to infinity". He then proves with the central limit theory that for shallow neural networks (i. e.one hidden layer) "a Gaussian prior for hidden-to-output weights results in a Gaussian process prior for functions", meaning that for any set of input points X, the random variables describing the functions computed at X have a joint Gaussian distribution [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]. In Figure 2.3 we show a pictorial representation of this phenomenon on a single hidden layer with tanh activation function and increasing width. This behavior has been observed in other deeper models (G. [START_REF] Matthews | Gaussian Process Behaviour in Wide Deep Neural Networks[END_REF][START_REF] Lee | Deep Neural Networks as Gaussian Processes[END_REF]; M. E. E. [START_REF] Khan | Approximate Inference Turns Deep Networks into Gaussian Processes[END_REF], for convolutional neural networks [START_REF] Garriga-Alonso | Deep Convolutional Networks as shallow Gaussian Processes[END_REF] and recurrent neural networks (G. [START_REF] Yang | Wide Feedforward or Recurrent Neural Networks of Any Architecture are Gaussian Processes[END_REF].

The learning problem is formulated as a transformation of a prior belief into a posterior distribution by means of Bayes' theorem. Given a dataset with N inputtarget pairs D = {X, y} def = {(x i , y i )} N i=1 , the posterior over θ is:

p(θ | y) = p(y | θ)p(θ) p(y) , (2.6)
where the likelihood depends on the output of the neural network at input X and can generally be defined as p(y | f(X; θ)). With the posterior inferred, we can make For keeping the notation clean, we will avoid to explicit condition on the input X.

prediction on new test points x , which can be done my marginalizing out the parameters from the posterior,

p(y | x , y) = p(y | θ)p(θ | y)dθ (2.7)
While looking remarkably easy, Equations (2.6) and (2.7) hide several challenges: (i) the non-conjugacy of the likelihood and the prior makes such that we don't know the form of the posterior; (ii) the normalization constant at the dominator requires to compute the integral on the space of the parameters of the joint likelihood. As a consequence, the posterior on the parameters of a b n n is always analytically intractable and making predictions on new test points isnot possible, unless we resort to approximations. In Figure 2.4 we visualize the challenges of Bayesian ReLu(w 0 • x -0.5) + ReLu(w 1 • x -0.3), where some parameters are xed for visualization purposes. On the left, the posterior on the two parameters w 0 and w 1 ;

on the right, the mean and the con dence intervals (±2σ) of the predictive posterior p(y | x , y). In both cases we approximate all required integrals with numerical integration.

inference in neural networks with a simple case where we can run numerical integration to solve the various integrals involved. The visualization of the density highlights the non-trivial multimodal shape of the posterior; this gives some intuition on the huge challenges that approximate inference techniques have to deal with.

.3 :

Variational inference (v i) is a classic tool to tackle intractable Bayesian inference [START_REF] Jordan | An Introduction to Variational Methods for Graphical Models[END_REF][START_REF] Blei | Variational Inference: A Review for Statisticians[END_REF]. vi casts the inference problem into an optimization-based procedure to compute a tractable approximation of the true posterior. In our setting, we have a probabilistic model p(y | f(X; θ)) with parameters θ, a prior distributions on them p(θ) and a set of observations {X, y}. In a nutshell, the general recipe of vi consists of (i) introducing a set Q of distributions;

(ii) defining a tractable objective that "measure" the distance between any arbitrary distribution q(θ) ∈ Q and the true posterior p(θ | y); and finally (iii) providing a programmatic way to find the distribution q(θ) that minimizes such distance (Figure 2.5 provides a simple diagram of this procedure). In practice, q(θ) has The KL divergence between two distributions q and p is KL [q(z) p(z)] = q(z) log q(z) p(z) dz.

some free parameters ν (also known as variational parameters), which are optimized such that the approximating distribution q(θ; ν) is as closer as possible to the true posterior p(θ | y). The classical formulation of vi chooses the Kullback-Leibler (kl) divergence as a metric of similarity between these two distribution, but extensions to alternative divergences and distances are possible as well (Y. [START_REF] Li | Rényi Divergence Variational Inference[END_REF][START_REF] Li | Dropout Inference in Bayesian Neural Networks with Alphadivergences[END_REF][START_REF] Dieng | Variational Inference via \chi Upper Bound Minimization[END_REF][START_REF] Wan | f-Divergence Variational Inference[END_REF]. We can derive the variational objective starting from the definition of the kl,

The color red highlights the quantities that we cannot compute.

KL [q(θ; ν) p(θ | y)] = E q(θ;ν) [log q(θ; ν) -log p(θ | y)] = (2.8)

= E q(θ;ν) [log q(θ; ν) -log p(y | θ) -log p(θ)] + log p(y)

Rearranging we have that log p(y) -KL [q(θ; ν) p(θ | y)] = E q(θ;ν) [log q(θ; ν) -log p(y | θ) -log p(θ)]

(2.9)

The r.h.s. of the equation defines our variational objective, also known as evidence lower bound (e l b o), that can be arranged as follows,

L elbo (ν) = E q(θ;ν) log p(y | θ) Model fitting term -KL [q(θ; ν) p(θ)]
Regularization term .

(2.10)

This formulation highlights the property of this objective, which is made of two components: the first one is the expected log-likelihood under the approximate posterior q and measures how the model fits the data. The second term, on the other hand, has the regularization effect of penalizing posteriors that are far from the prior as measured by the k l. We finally observe that minimizing the k l in Equation (2.8) is equivalent to find the variational parameters ν such that the elbo is maximized,

ν = arg max ν E q(θ;ν) log p(y | θ) -KL [q(θ; ν) p(θ)]
(2.11)

Before diving into the challenges of Equation (2.11), we shall spend a brief moment discussing the form of the approximating distribution q. One of the simplest and easier choice is the mean field approximation (G. E. [START_REF] Hinton | Keeping the Neural Networks Simple by Minimizing the Description Length of the Weights[END_REF],

where each variable θ i is taken to be independent with respect to the remaining θ -i . Effectively, this imposes a factorization of the posterior,

q(θ; ν) = K i=1 q(θ i ; ν i ) (2.12)
where ν i is the set of variational parameters for the parameter θ i . On top of this approximation, q(θ i ) is often chosen to be Gaussian,

q(θ i ) = N(µ i , σ 2 i ) (2.13)
Now, the collection of all means and variances {µ i , σ 2 i } K i=1 defines the set of variational parameters to optimize. In the next chapters we will discuss the limitations of simple approximations for Bayesian neural networks and how things can be improved.

For b n ns the analytic evaluation of the e l b o (and its gradients) is always untractable due the non-linear nature of the expectation of the log-likelihood under the variational distribution. Nonetheless, this can be easily estimated via Monte Carlo integration [START_REF] Metropolis | The Monte Carlo Method[END_REF], by sampling N MC times from

q ν , E q(θ;ν) log p(y | θ) ≈ 1 N MC N MC j=1 log p(y | θ j ) , with θ j ∼ q(θ; ν) (2.14)
In practice, this is as simple as re-sampling the weights and the biases for all the layers N MC times and computing the output for each new sample. Despite its simplicity, the estimation in Equation (2.14) is unbiased and its variance asymptotically shrinks with 1 N MC , independently of the dimensionality of θ. Equation (2.9) allows us to understand the relationship between the marginal likelihood and the (negative) evidence lower bound. In the ideal case of exact solution for the posterior, equivalent to having KL [q(θ; ν) p(θ | y)] = 0, the minimization of the e l b o is equivalent to the maximization of the marginal likelihood, with all the properties analyzed in § 2.1.1. In the general case, the e l b o is just a lower bound of the marginal and the gap between these two quantities depends on the level of approximation of q. To narrow this gap and therefore obtain a tighter lower bound, we need to increase the expressiveness of the variational posterior.

. 3 .1 Optimization of the

We now have a tractable objective that needs to be optimized with respect to the variational parameters ν. Very often the kl term is known, making its differentation trivial. On the other hand the expectation of the likelihood is not available, making the computation of its gradients more challenging. This differentiation falls in the category of taking gradients of a quadrature, i. e.∇ ν E q(θ;ν) f(θ). An unbiased estimate of such gradient can be derived as follows,

∇ ν E q f(θ) = E q [f(θ)∇ ν log q(θ)]
(2.15)

which can be estimated again by sampling from q. The full gradient of the e l b o can be computed starting from Equation (6.2) as follows,

∇ ν L e l b o = E q [p θ (y | θ)∇ ν log q ν (θ)] -∇ ν KL [q ν (θ) p(θ)] . (2.16)
This estimator has a very general formulation that make it applicable to a wide range of situations. In the literature, this is also refer as score function estimator or black-box variational inference (D. M. [START_REF] Blei | Variational Bayesian Inference with Stochastic Search[END_REF], due to the fact that it doesn't require to take the derivative of the model likelihood p(y | θ). As such, the score function estimator can be used also for non-differentiable likelihoods. On the other hand, this estimator shows huge variance and, unless remedies like control variates [START_REF] Ross | Simulation, Fourth Edition[END_REF] are taken, this estimator is impractical. Given that most of b n ns architectures (if not all) are easily differentiable, there are not practical uses for this formulation in Bayesian deep learning.

Alternatively, this problem can be solved using the so-called reparameterization trick [START_REF] Salimans | Fixed-Form Variational Posterior Approximation through Stochastic Linear Regression[END_REF][START_REF] Kingma | Auto-Encoding Variational Bayes[END_REF]. The reparameterization trick aims at constructing θ as an invertible function T of the variational parameters ν and of another random variable ε, so that θ = T(ε; ν). ε is chosen such [START_REF] Maddison | The Concrete Distribution: A Continuous Relaxation of Discrete Random Variables[END_REF].

that its marginal p(ε) does not depend on the variational parameters. With this parameterization, T separates the deterministic components of q from the stochastic ones, making the computation of its gradient straightforward. For a Gaussian distribution with mean µ and variance σ 2 , T corresponds to as simple scale-location transformation of an isotropic Gaussian noise,

θ ∼ N(µ, σ 2 ) ⇐⇒ θ = µ + σε with ε ∼ N(0, 1) .
(2.17)

This simple transformation ensures that p(ε) = N(0, 1) does not depends on the variational parameters ν = {µ, σ 2 }. The gradients of the e l b o can be therefore computed as

∇ ν L e l b o = E p(ε) ∇ θ log p(y | θ) | θ=T(ε;ν) ∇ ν T(ε; ν) -∇ ν KL [q(θ; ν) p(θ)] .
(2.18)

The gradient ∇ θ log p(y | θ) depends on the model and it can be derived with automatic differentation tools (Abadi et al., 2015;[START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF], while ∇ ν T(ε; ν) doesn't have any stochastic components and therefore can be known deterministically. Note that the reparameterization trick can be also used when the k l is not analitically available. In that case, we would end up with, [START_REF] Roeder | Sticking the Landing: Simple, Lower-Variance Gradient Estimators for Variational Inference[END_REF] argue that when we believe that q(θ; ν) ≈ p(y | θ), Equation (2.19) should be prefered over Equation (2.18) even if computing analitically the kl is possible. Note that this case is very unlikely for bnn posteriors, and that the additional randomness introduced by the Monte Carlo estimation of the k l could be harmful.

∇ ν L e l b o = E p(ε) [∇ θ log p(y | θ) + log q(θ; ν) -log p(θ)] θ=T(ε;ν) ∇ ν T(ε; ν) (2.19)
In case of large datasets and complex models, the formulation summarized in Equation (2.18) can be computationally challenging, due to the evaluation of the likelihood and its gradients N MC times. Assuming factorization of the likelihood, this quantity can be approximated using mini-batching [START_REF] Graves | Practical Variational Inference for Neural Networks[END_REF][START_REF] Hoffman | Stochastic Variational Inference[END_REF]. Recalling y as the set of labels of our dataset with N examples, by taking B ⊂ y as a random subset of y, the likelihood term can be estimated in an unbiased way as

p(y | θ) = p(y | f(X; θ)) = N i=1 p(y i | f(x i ; θ)) (2.20)
log p θ (y | θ) ≈ N M y i ∼B log p(y i | θ) . (2.21)
where M is the number of points in the minibatch. At the cost of increase "randomness", we can use Equation (2.18) to compute the gradients of the elbo with the minibatch formulation in Equation (2.21). Stochastic optimization, e.g. any version of stochastic gradient descent (s g d), will converge to a local optimum provided with a decreasing learning rate and sufficient gradient updates [START_REF] Robbins | A Stochastic Approximation Method[END_REF].

This procedure can be further improved by analyzing the magnitude of the noise for the stochastic gradients. Focusing on the likelihood in Equation (2.21), let's define i the likelihood contribution log p(y i | θ) from the datapoint y i . From here, we can derive the variance of this term by considering the minibatch to be drawn randomly with replacement from y [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF],

Var E q(θ;ν) log p(y | θ) = N 2 M 2   M i=0 Var [ i ] + 2 M i=0 M j=i+1 Cov i , j   ≈ N 2 1 M Var [ i ] + M -1 M Cov i , j (2.22)
where the Var [•] and the Cov [•] are taken with respect to the empirical distribution defined by the dataset and the distribution of θ. As we can see, the variance term decreases linearly with M but the total contribution does not scale with the number of points in the minibatch. This is due to the correlation term which is asymptotically independent on the minibatch size. In practice, this term can dominate the whole variance of the likelihood even for big minibatch. This can be alleviated by imposing independence between cross-terms in the likelihood, i.e. Cov i , j = 0. This is achieved by imposing independence between input points and the samples from the distribution q ν (θ); in practice, one would need to resample θ exactly N MC times for each point y i . This is generally computationally challenging but for b n ns it becomes as easy as moving the randomness from the weights to the pre-activations. Take as an example the setup of Figure 2.7 where for conveniance we defined A as the pre-activation map for the first hidden layer of a mini-batch of data X. Given our choice of factorized Gaussian on the weights

q(W ij ) = N(µ ij , σ 2 ij )
, we can write down the expression for the distribution over A, which has the following form,

q(A; ν) = M m=1 D 1 j=1 N D 0 k=1 X nk µ kj , D 0 k=1 X 2 nk σ 2 kj (2.23)
This new estimator is also known as the local reparameterization trick [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF], if otherwise stated it will be the default choice for variational inference on b n ns.

.4 29

x y

f (1) f (2) W X F = a( A XW)

.4

Hamiltonian Monte Carlo (hmc) methods provide an elegant and efficient way to sample from intractable distributions by defining proposals with probability one of being accepted in the classic Metropolis-Hastings framework. These proposals are generated by simulating an Hamiltonian system where the potential energy U is our intractable (log-)posterior densitylog p(θ | y) and the kinetic energy is parameterized with some auxilary momentum variables ρ. To propose samples, h m c simulates the following dynamics, derived from the theory of Hamiltonian systems (Hamilton-Jacobi equations):

dθ = M -1 ρdt (2.24) dρ = ∇U(θ)dt (2.25)
While the simulation of this continuous system is often intractable, the classic recipe calls for discretization using e.g. the leapfrog scheme and correction step using the classic MH acceptance/rejection framework. In this case, the MH works on the Hamiltonian function, which is defined as

H(θ, ρ) = U(θ) + 1 2 ρ M -1 ρ (2.26)
Algorithm 1 sketches the h m c sampler, including the steps required for the discretization of the dynamics.

Algorithm 1: Hamiltonian Monte Carlo

Input: Starting position θ (1) and step size ε for t = 1, 2 • • • do Resample momentum:

ρ (t) ∼ N(0, M) (θ 0 , ρ 0 ) = (θ (t) , ρ (t) )
Leapfrog discretization of Hamiltonian dynamics in Equation (2.24):

ρ 0 ← ρ 0 -ε 2 ∇U(θ 0 ) for i = 1 to m do θ i ← θ i-1 + εM -1 ρ i-1 ρ i ← ρ i-1 -ε∇U(θ i ) ρ m ← ρ m -ε 2 ∇U(θ m ) ( θ, ρ) = (θ m , ρ m ) Metropolis-Hastings correction: u ∼ Uniform(0, 1) ρ = e H( θ, ρ)-H(θ (t) ,ρ (t) ) if u < min(1, ρ), then θ (t+1) = θ
This way of proposing samples is very efficient but computationally challenging, as it requires to compute the gradients of the unnormalized log-posterior. While the prior component is usually not a problem, the evaluation of the gradient of the likelihood is, especially for large datasets. If we suppose the likelihood to factorize on the observations, we can approximate this quantity using mini-batches:

∇ U(θ) = - N M M i=1 ∇ log p(y i | θ) -∇ log p(θ) (2.27)
Using mini-batches, we can write the gradients as unbiased estimates of the original

U(θ) as follows, ∇ U(θ) = ∇U(θ) + N(0, V(θ))
(2.28)

where V(θ) is the covariance of the stochastic gradient noise, which depends on many factors including the parameters themselves.

A first solution to implement Stochastic Gradient Hamiltonian Monte Carlo (sghmc) is to replace the computation of the gradients in Algorithm 1 with Equation (2.27). Unfortunately, this is not enough: as T. [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF] proves, the simulation of the resulting dynamics would make the stationary distribution no longer invariant (i.e., sampling is performed on a distribution which is not the original p(θ | y)).

To solve this problem, T. [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF] proposes to modify the dynamics by introducing an additional friction term C. With this addition, the discretized system of equation becomes,

∆θ = εM -1 ρ (2.29) ∆ρ = ε∇ U(θ) -εCM -1 ρ + N(0, 2(C -B)ε) (2.30)
where B is an estimation of the stochastic gradient noise covariance.

From a practical point of view, this system of dynamics is still challenging to work

with. Namely, one should choose the friction term C, perform an estimation of the noise gradient B and choose the mass matrix M and the step-size ε. While the friction term and the step-size are high model and dataset dependent, the other two quantities can be estimated during the burn-in phase [START_REF] Springenberg | Bayesian Optimization with Robust Bayesian Neural Networks[END_REF]. For the mass matrix, we can leverage the connection between sghmc and sgd. Results in stochastic optimization literature [START_REF] Duchi | Adaptive Subgradient Methods for Online Learning and Stochastic Optimization[END_REF][START_REF] Tieleman | RMSProp: Divide the gradient by a running average of its recent magnitude[END_REF] discuss how the robustness of s g d can be improved by normalizing the gradient by its magnitude. For s g h m c this is equivalent to choosing,

M -1 = diag V -1 2 θ (2.31)
where V θ is the estimation of the element-wise variance of the gradient and it can be estimated by exponential moving average as follow,

∆ V θ = -τ -1 V θ + τ -1 ∇(U(θ)) 2
(2.32)

The averaging window is specified by the choice of τ, which-again-can be automatically determined by using a procedure similar to adaptive learning rate for s g d [START_REF] Tieleman | RMSProp: Divide the gradient by a running average of its recent magnitude[END_REF],

∆τ = -g 2 θ V -1 θ τ + 1 (2.33) ∆g θ = -τ -1 g θ + τ -1 ∇ U(θ) (2.34)
where g θ is a smoothed estimation of the gradients ∇U(θ). Finally, we can also use V θ for B, specifically B = 1 2 V θ ε. By substituting these new quantities in the original dynamics we obtain

∆θ = ν (2.35) ∆ν = -ε 2 V -1 2 θ ∇ U(θ) -ε V -1 2 θ Cν + N 0, 2ε 3 V -1 2 θ C V -1 2 θ -ε 4 I (2.36) where ν is defined as ν = εM -1 ρ = ε V -1 2 θ ρ. I N I T I A L I Z A T I O N S O F V A R I A T I O N A L I N F E R E N C E F O R B A Y E S I A N N E U R A L N E T W O R K S

Stochastic variational inference is an established way to carry out approximate

Bayesian inference for deep models flexibly and at scale. While there have been effective proposals for good initializations for loss minimization in deep learning, far less attention has been devoted to the issue of initialization of stochastic variational inference. In this chapter we analyze this problem for Bayesian deep neural network and we propose a novel layer-wise initialization strategy based on Bayesian linear models. The proposed method is extensively validated on regression and classification tasks, including Bayesian deep neural networks and Bayesian convolutional neural networks, showing faster and better convergence compared to alternatives inspired by the literature on initializations for loss minimization.

.1

Deep neural networks (d n ns) and convolutional neural networks (c n ns) have become the preferred choice to tackle various learning tasks, due to their ability to model complex problems and the mature development of regularization techniques to control overfitting (LeCun, [START_REF] Lecun | Deep learning[END_REF][START_REF] Srivastava | Dropout: A Simple Way to Prevent Neural Networks from Overfitting[END_REF]. There has been a recent surge of interest in the issues associated with their over-confidence in predictions, and proposals to mitigate these [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF][START_REF] Kendall | What Uncertainties Do We Need in Bayesian Deep Learning for Computer Vision?[END_REF][START_REF] Lakshminarayanan | Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles[END_REF]. Bayesian techniques offer a natural framework to deal with such issues, but they are characterized by computational intractability (C. M. [START_REF] Bishop | Pattern recognition and machine learning[END_REF][START_REF] Ghahramani | Probabilistic Machine Learning and Artificial Intelligence[END_REF]. A popular way to recover tractability is to use variational inference [START_REF] Jordan | An Introduction to Variational Methods for Graphical Models[END_REF]. In variational inference (v i), an approximate posterior distribution is introduced and its parameters are adapted by optimizing a variational objective, which is a lower bound to the marginal likelihood. Stochastic variational inference offers a practical way to carry out stochastic optimization of the variational objective. Here stochasticity is introduced with a doubly stochastic approximation of the expectation term, which is unbiasedly approximated using Monte Carlo and by selecting a subset of the training points [START_REF] Graves | Practical Variational Inference for Neural Networks[END_REF][START_REF] Kingma | Auto-Encoding Variational Bayes[END_REF]. While v i is an attractive and practical way to perform approximate inference for d n ns, there are limitations. For example, the form of the approximating distribution can be too simple to accurately approximate complex posterior distributions [START_REF] Ha | HyperNetworks[END_REF][START_REF] Ranganath | Deep Exponential Families[END_REF][START_REF] Rezende | Variational Inference with Normalizing Flows[END_REF]. Furthermore, v i increases the number of optimization parameters compared to optimizing model parameters through, e.g., loss minimization; for example, a fully factorized Gaussian posterior over model parameters doubles the number of parameters in the optimization compared to loss minimization. This has motivated research on other ways to perform approximate Bayesian inference for dnns by establishing connections between variational inference and dropout (Gal and Ghahramani, 2016b;Gal and Ghahramani, 2016a;[START_REF] Gal | Concrete Dropout[END_REF].

.1 35

A theoretical understanding of the optimization landscape of d n ns and c n ns is still in its early stages of development [START_REF] Dziugaite | Computing Nonvacuous Generalization Bounds for Deep (Stochastic) Neural Networks with Many More Parameters than Training Data[END_REF][START_REF] Garipov | Loss Surfaces, Mode Connectivity, and Fast Ensembling of DNNs[END_REF], and most works have focused on the practical aspects characterizing the optimization of their parameters [START_REF] Duchi | Adaptive Subgradient Methods for Online Learning and Stochastic Optimization[END_REF][START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF][START_REF] Srivastava | Dropout: A Simple Way to Prevent Neural Networks from Overfitting[END_REF]. If this lack of theory is apparent for optimization of model parameters, this is even more so for the understanding of the optimization landscape of the objective in variational inference, where variational parameters enter in a nontrivial way in the objective [START_REF] Graves | Practical Variational Inference for Neural Networks[END_REF][START_REF] Rezende | Stochastic Backpropagation and Approximate Inference in Deep Generative Models[END_REF]. Initialization plays a huge role in the convergence of vi; the illustrative example in Figure 3.1 shows how a poor initialization can prevent vi to converge to good solutions in short amount of time. The problem is even more severe for complex architectures, such as the ones that we discuss in the experiments; for example, vi systematically converges to trivial solutions (posterior equal to the prior) when applied to c n ns. In this chapter, we focus on this issue affecting v i for d n ns and c n ns. While there is an established literature on ways to initialize model parameters of d n ns when minimizing its loss [START_REF] Glorot | Understanding the Difficulty of Training Deep Feedforward Neural Networks[END_REF][START_REF] Saxe | Exact Solutions to the Nonlinear Dynamics of Learning in Deep Linear Neural Networks[END_REF][START_REF] Mishkin | All You Need is a Good Init[END_REF], to the best of our knowledge, there is no study that systematically tackles this issue for v i for Bayesian d n ns and c n ns. Inspired by the literature on residual networks (K. [START_REF] He | Deep Residual Learning for Image Recognition[END_REF] and greedy initialization of d n ns [START_REF] Bengio | Greedy Layer-Wise Training of Deep Networks[END_REF][START_REF] Mishkin | All You Need is a Good Init[END_REF], we propose a novel initialization strategy for v i grounded on Bayesian linear modeling, which we call i b l m. Iterating from the first layer, i b l m initializes the posterior at layer (l) by learning Bayesian linear models which regress from the input, propagated up to layer (l), to the labels.

We show how i b l m can be applied in a scalable way and without considerable overhead to regression and classification problems, and how it can be applied to initialize vi not only for dnns but also for cnns. Through a series of experiments, we demonstrate that iblm leads to faster convergence compared to other initializations inspired by prior work on loss minimization for dnns. Furthermore, we show that iblm makes it possible for vi with a Gaussian approximation applied to cnns to compete with Monte Carlo dropout (m c d) (Gal and Ghahramani, 2016a) and noisy natural gradients (noisy-kfac; G. [START_REF] Zhang | Noisy Natural Gradient as Variational Inference[END_REF]), which are state-of-art methods to perform approximate inference for cnns.

In summary, in this chapter we make the following contributions:

• we propose a novel way to initialize v i for d n ns based on Bayesian linear models;

• we show how this can be done for regression and classification;

• we show how to apply our strategy to cnns;

• we empirically demonstrate that our proposal allows us to achieve performance superior to other initializations of vi inspired by the literature on loss minimization;

• for the first time, we are able to run Gaussian v i on large-scale c n ns, obtaining remarkable performances.

.1 .1 A review of the role of initialization in deep learning

The problem of initialization of weights and biases in d n ns for gradient-based loss minimization has been extensively tackled in the literature since early breakthroughs in the field [START_REF] Rumelhart | Learning Internal Representations by Error Propagation[END_REF][START_REF] Baldi | Neural Networks and Principal Component Analysis: Learning from Examples without Local Minima[END_REF]. et al. (1986) study the limitations of backpropagation-based gradient descend to find good local minima in the parameter space, showing that this issue arises even in small dnns; the authors address this issue by initalizing the optimization using random weights. Few years later, [START_REF] Baldi | Neural Networks and Principal Component Analysis: Learning from Examples without Local Minima[END_REF] report that for layered linear feed-forward neural networks trained using a quadratic loss function, there exists a unique global minimum corresponding to a orthogonal projection of the ordinary least square solution into a subspace spanned by eigenvalues of the covariance matrix of training samples. Nevertheless, the restrictions required to derive this closed-form solution (linear layers without bias on a single hidden layer networks with equal number of input and output features) severely limit the applicability of these results to a wider and more interesting class of d n ns. To overcome these limitations, [START_REF] Lecun | Gradient-Based Learning Applied to Document Recognition[END_REF] discusses practical tricks to achieve an efficient loss minimization through back-propagation; choosing separate learning rate for each weight and a Gaussian random initialization are some of the proposed ways to obtain an effective optimization.

Rumelhart
More recently, [START_REF] Bengio | Greedy Layer-Wise Training of Deep Networks[END_REF] propose a greedy layer-wise unsupervised pre-training, which proved to help optimization and generalization. A justification can be found in [START_REF] Erhan | Why Does Unsupervised Pre-training Help Deep Learning?[END_REF], where the authors show that pre-training can act as regularization; by initializing the parameters in a region corresponding to a better basin of attraction for the optimization procedure, the model can reach a better local minimum and increase its generalization capabilities. [START_REF] Glorot | Understanding the Difficulty of Training Deep Feedforward Neural Networks[END_REF] propose a simple way to estimate the variance for random initialization of weights, which makes it possible to avoid saturation both in forward and back-propagation steps. Another possible strategy can be found in the work by [START_REF] Saxe | Exact Solutions to the Nonlinear Dynamics of Learning in Deep Linear Neural Networks[END_REF], that investigates the dynamics of gradient descend optimization, and proposes a random orthogonal initialization of the weights based on the singular value decomposition of a Gaussian random matrix. This algorithm takes a random weight matrix filled with Gaussian noise, decomposes it to orthonormal basis using a singular value decomposition and replaces the weights with one of the components. Building on this work, [START_REF] Mishkin | All You Need is a Good Init[END_REF] propose a data-driven weight initialization by scaling the orthonormal matrix of weights to make the variance of the output as close to one as possible.

Variational inference addresses the problem of intractable Bayesian inference by reinterpreting inference as an optimization problem. The reparameterization trick [START_REF] Kingma | Auto-Encoding Variational Bayes[END_REF] allows for the stochastic optimization of the variational lower bound through automatic differentiation [START_REF] Duchi | Adaptive Subgradient Methods for Online Learning and Stochastic Optimization[END_REF][START_REF] Zeiler | ADADELTA: An Adaptive Learning Rate Method[END_REF][START_REF] Sutskever | On the importance of initialization and momentum in deep learning[END_REF][START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF]. Despite the tight connection between loss minimization and evidence lower bound (elbo) maximization [START_REF] Graves | Practical Variational Inference for Neural Networks[END_REF], to the best of our knowledge there is no study that either empirically or theoretically addresses the problem of initialization of parameters for v i; we could only find a mention of this in [START_REF] Krishnan | On the Challenges of Learning with Inference Networks on Sparse, High-Dimensional Data[END_REF] for variational autoencoders. We aim to fill this gap by proposing a novel way to initialize parameters in v i for probabilistic deep models.

.2 :

In this section, we introduce our proposed Iterative Bayesian Linear Model (iblm) initialization for vi. We first introduce iblm for regression with dnns, and we then show how this can be extended to classification and to cnns.

. 2 .1 Initialization of DNNs for Regression

In order to initialize the weights of d n ns, we proceed iteratively as follows.

Before applying the nonlinearity through the activation function, each layer in a Bayesian d n n can be seen as multivariate Bayesian linear regression model.

We use this observation as an inspiration to initialize the variational parameters as follows. Starting from the first layer, we can set the parameters of q(W (0) ) by running Bayesian linear regression with inputs X and labels y. Let's considering for a moment a particular column j of the weights matrix, which we denote

X y q(W (0) ,0 ) (a) X y q(W (0) ,1 ) (b) X y q(W (1) ) (c
with w = W (0)
•j . By choosing a Gaussian prior and a Gaussian likelihood as follows

p(w) = N(w | 0, Λ) (3.1) p(y | w, X) = N(y | Xw, L) (3.2)
the posterior p lm (w | y, X) as simple as apply the rule of Gaussian random variables (C. M. [START_REF] Bishop | Pattern recognition and machine learning[END_REF],

p lm (w | y, X) = N(w | m, S) = N(w | SX L -1 y, S) , (3.3)
where

S -1 = Λ -1 + X L -1 X.
This process is eventually repeating the process for multiple hidden units (the posterior factorizes on the columns of W (0) ). After this, we initialize the approximate posterior over the weights at the second layer q(W (1) )

by running Bayesian linear regression with inputs X = a(X W (0) ) and labels y,

where W (0) is a sample from q(W (0) ). We then proceed iteratively in the same way

Remember that a(•) denotes for us the element-wise application of the activation function to the argument up to the last layer. Figure 3.2 gives an illustration of the proposed method for a simple architecture.

The intuition behind i b l m is as follows. If one layer is enough to capture the complexity of a regression task, we expect to be able to learn an effective mapping right after the initialization of the first layer. In this case, we also expect that the mapping at the next layers implements simple transformations, close to the identity. Learning a set of weights with these characteristics starting from a random initialization is far from trivial, which also motivated the work on residual networks (K. [START_REF] He | Deep Residual Learning for Image Recognition[END_REF]. Our iblm initialization takes this observation as an intuition to initialize v i for general deep models.

From a complexity point of view, denoting by D l the number of output neurons at layer (l), this is equivalent to D l univariate Bayesian linear models. Instead of using the entire training set to learn the linear models, each one of these is inferred based on a random mini-batch of data, whose inputs are propagated through the previous layers. The complexity of iblm is linear in the batch size and cubic in the number of neurons to be initialized. Later on, we will provide an evaluation of the effect of batch size and a timing profiling of iblm.

. 2 . 2 From the Bayesian linear model posterior to the variational approximation

The proposed i b l m initialization of variational parameters can be used with any choice for the form of the approximate posterior. The exact posterior of Bayesian linear regression is not factorized, so one needs to match this with the form of the chosen approximate posterior. For simplicity of notation, let w be one of the parameters vector of interest. We can formulate this problem by minimizing the Kullback-Leibler (k l) divergence from the variational approximation q(w; ν) to the posterior obtained from the linear model p lm (w | X, y). In the case of a fully factorized approximate posterior over the weights q(w; ν) = N(w | µ, Iσ 2 ), this minimization can be done analytically.

Recall that the expression for the kl divergence between two multivariate Gaussians

p 0 = N(w | µ 0 , Σ 0 ) and p 1 = N(w | µ 1 , Σ 1 ) is KL [p 0 p 1 ] = 1 2 Tr Σ -1 1 Σ 0 + 1 2 (µ 1 -µ 0 ) Σ -1 1 (µ 1 -µ 0 ) - D 2 + 1 2 log detΣ 1 detΣ 0 (3.4)
The kl divergence is not symmetric, so the order in which we take this matters. If we consider KL [q(w) p lm (w | X, y)], the expression becomes:

KL [q(w) p(w | X, y)] = 1 2 Tr S -1 diag(σ 2 ) + 1 2 (m -µ) S -1 (m -µ) - D 2 + 1 2 log detS i σ 2 i (3.5)
It is a simple matter to show that the optimal mean µ is m as µ appears only in the quadratic form which is clearly minimized when µ = m. For the variances σ 2 , we need to take the derivative of the k l divergence and set it to zero:

∂KL [q(w) p(w | X, y)] ∂σ 2 i = 1 2 ∂ Tr S -1 diag(σ 2 ) ∂σ 2 i - 1 2 ∂ i log σ 2 i ∂σ 2 i = 0 (3.6)
Rewriting the trace term as the sum of the Hadamard product of the matrices in

the product ij (S -1 diag(σ 2 )) ij = i σ 2 i (S -1 ) ii , this yields ∂KL [q(w) p(w | X, y)] ∂σ 2 i = 1 2 ∂σ 2 i (S -1 ) ii ∂σ 2 i - 1 2 ∂ log σ 2 i ∂σ 2 i = 0 (3.7)
This results in (σ 2 i ) -1 = (S -1 ) ii . This approximation has the effect of underestimating the variance for each variable [START_REF] Murphy | Machine learning -A Probabilistic Perspective[END_REF]. In case we consider KL [p(w | X, y) q(w)], the expression for the mean would remain the same but the variances would become σ 2 i = S ii , which is the simplest way to approximate the correlated posterior over w but it is going to inflate the variance in case of strong correlations. Similar results can be also obtained for different posterior distributions, such as Gaussian posteriors with full or low-rank covariance, or matrix-variate Gaussian posteriors [START_REF] Louizos | Structured and Efficient Variational Deep Learning with Matrix Gaussian Posteriors[END_REF].

. 2 . Initialization for classi cation and convolutional layers

In this section we show how our proposal can be extended to k-class classification problems. We assume a one-hot encoding of the labels, so that Y is an n × k matrix of zeros and ones (one for each row of Y). Recently, it has been shown that it is possible to obtain an accurate modeling of the posterior over classification functions by applying regression on a transformation of the labels [START_REF] Milios | Dirichletbased Gaussian Processes for Large-scale Calibrated Classification[END_REF]. This is interesting because it allows us to apply Bayesian linear regression as before in order to initialize v i for d n ns. The transformation of the labels is based on the formalization of a simple intuition, which is the inversion of the softmax transformation. One-hot encoded labels are viewed as a set of parameters of a degenerate Dirichlet distribution. We resolve the degeneracy of the Dirichlet distribution by adding a small regularization, say α = 0.01, to the parameters. At this point, we leverage the fact that Dirichlet distributed random variables can be constructed as a ratio of Gamma random variables, that is, if

x i ∼ Gamma(a i , b), then x i j x j
∼ Dir(a). We can then approximate the Gamma random variables with log-Normals by moment matching, which become Gaussian after a logarithm transformation. By doing so, we obtain a representation of the labels which allows us to use standard regression with a Gaussian likelihood, and which retrieves an approximate Dirichlet when mapping predictions back using the softmax transformation. As a result, the latent functions obtained represent probabilities of class labels. The only small complication is that the transformation imposes a different noise level for labels that are 0 or 1, and this is due to the nonsymmetric nature of the transformation. Nevertheless, it is a simple matter to extend Bayesian linear regression to handle heteroscedasticity; see the supplementary material and [START_REF] Milios | Dirichletbased Gaussian Processes for Large-scale Calibrated Classification[END_REF] for more insights on the transformation to apply regression on classification problems. i b l m can also be applied to c n ns. Convolutional layers are commonly implemented as matrix multiplication (e.g. as a linear model) between a batched matrix of patches and a reshaped filter matrix [START_REF] Jia | Learning Semantic Image Representations at a Large Scale[END_REF]. Rather than using the outputs of the previous layer as they are, for convolutional layers each Bayesian linear model learns the mapping from spatial patches to output features. In Algorithm 2 we sum- marize a sketch of the proposed method for regression as well as for classification and convolutional layers. X IBLM ← output of previous layer; if layer is convolutional then

X IBLM ← patch extraction(X IBLM ); if likelihood is classification then var(Y IBLM ) ← log [(Y + α) -1 + 1]; mean(Y IBLM ) ← log (Y + α) -var(Y IBLM )/2; else Y IBLM ← Y; Compute p lm (w | X, Y);
Initialize q(w) with the best approx. of p lm (w | X, Y)

.3

In this section, we compare different initialization algorithms for variational inference to prove the effectiveness of i b l m. At layer (l), we choose priors p(W (l) ) =

3 . 3 43 i,j N(w (l) i,j | 0, 1 D (l-1) )
, where D (l-1) denotes the number of input features at layer (l), and focus on fully-factorized variational posteriors q(W (l) ) = i,j N(w

(l) i,j | µ (l) i,j , (σ 2 ) (l) i,j );
We propose a number of competitors inspired from the literature developed for loss minimization in d n ns and c n ns:

• Uninformative. The posteriors at each layer are initialized with zero mean and unit variance.

• Random Heuristic. An extension to commonly used heuristics with µ

(l) i,j = 0 and (σ 2 ) (l) i,j = 1 D (l-1)
. Because this is the same as for the prior, this yields an initial k l divergence in the e lbo equal to zero.

• Xavier Normal. Originally proposed by [START_REF] Glorot | Understanding the Difficulty of Training Deep Feedforward Neural Networks[END_REF], it samples all weights independently from a Gaussian distribution with zero mean and

(σ 2 ) (l) i,j = 2 D (l-1) +D (l)
. This variance-based scaling avoids issues with vanishing or exploding gradients. We extend this to vi by directly setting µ (l) i,j = 0 and

(σ 2 ) (l) i,j = 2 D (l-1) +D (l)
, given that the sampling is performed during the Monte Carlo estimate of the log-likelihood.

• Orthogonal. Starting from the analysis of learning dynamics of d n ns with linear activations, [START_REF] Saxe | Exact Solutions to the Nonlinear Dynamics of Learning in Deep Linear Neural Networks[END_REF] propose an initialization scheme with orthonormal weight matrices. The idea is to decompose a Gaussian random matrix onto an orthonormal basis, and use the resulting orthogonal matrix for initialization. We adapt this method to vi by initializing the mean matrix with the orthogonal matrix and (σ 2 )

(l) i,j = 1 D (l-1) .
In the experiments, we make use of the Pytorch QR-decomposition [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF].

• Layer-Sequential Unit-Variance (LSUV). Starting from the orthogonal initialization, [START_REF] Mishkin | All You Need is a Good Init[END_REF] propose a data-driven greedy layer-wise variance scaling of the weight matrices. We implement Layer-Sequential Unit-Variance (lsuv) for the means, while the variances are set to (σ 2 )

(l) i,j = 1 D (l-1) .
Figure 3.3 presents an illustration of the different initialization for the variational parameters in a simple 1D regression task. i b l m being data dependent already shows some visible fitting behaviour. We report an extensive validation of our proposal, with a series of experiments involving dnns and cnns for regression and classification, comparing it against other initializations for v i. For bigger models, like the cnns, we will benchmark mean-field Gaussian variational inference with mcd (Gal and Ghahramani, 2016b) and Natural Noisy Gradients or noisy-kfac run on a server equipped with two 16c/32t Intel Xeon CPU and four NVIDIA Tesla P100, with a maximum time budget of 24 hours (never reached). To better understand the effectiveness of different initializations, all learning curves are plotted w.r.t. training iteration rather than wall-clock time.

. .1 The e ect of initialization in deep variational neural networks

In this experiment we compare initialization methods for a shallow d n n architecture on two regression datasets, Powerplant (N = 9568, D = 4) and Protein (N = 45730, D = 9). The architecture used in these experiments has one single hidden layer with 100 hidden neurons and rectified linear unit (ReLU) activations.

Figure 3.4 on the top row shows the learning curves repeated over five different train/test splits. First of all we observe that initialization hugely impact the convergence properties of variational inference. On Powerplant, for instance, the impact is easily quantified: between the fastest and the slowest configurations to converge there are differences as big as 30/40 times more training interations. This is an expected result and a trivial conclusion, but the effects of poor choices of initialization schemes in the convergence of variational inference were never discussed.

Furthermore, i b l m allows for a better initialization compared to the competitors, leading to lower m n l l on the test set (curves for the root mean squared error (r m s e) are similar and available in the Appendix). Similar considerations hold when increasing the depth of the model, keeping the same experimental setup. Fig-

ure 3.4
shows what happens we we switch to a deeper architecture with five hidden layers and 100 hidden neurons per layer (R eLU activations).

. . 2 Scaling up variational inference to deep convolutional neural networks

For this experiment, we implemented a Bayesian version of the original LeNet architecture proposed by [START_REF] Lecun | Gradient-Based Learning Applied to Document Recognition[END_REF] with two convolutional layers of 6 and 16 filters, respectively and ReLU activations applied after all convolutional layers and fully-connected layers. We tested our framework on MNIST and on CIFAR10 .

Before showing and discussing the comparisons, this is a challenging setup to validate two claims we made in the previous sections. Previously we claimed that (i) small batches of data are sufficient to solve the Bayesian linear model and that (ii) our initialization does not incur significant overheads. To justify such claims, we initialize a LeNet on MNIST with an increasing number of samples per batch; Finally, Figure 3.6 reports the test mnll after initialization as a function of the time required to initialize the model (orange points correspond to Pareto-optimal points).

The timings and likelihood values are obtained by repeating the experiment eight time independently on a graphic processing unit (gpu) fully dedicated. The iblm initializer takes roughly one minute to finish, which in the full timing budget of training such models is negligible.

Now we can move to the comparisons, to which we include a further initialization based on the maximum-a-posteriori (map). For this case we optimize the map loss with the same prior and we use this solution to initialize the µ (l) i,j , while we set the variances to a small value of log[(σ 2 ) (l) i,j ] = -5.5 (for a fair comparison, we allow the m a p to train for the same amount of time required by i b l m to complete). Even for relatively small c n ns, we observe somethig alarming: the only initialization strategies that achieve convergence are the orthogonal and lsuv, along with iblm and map, which are both data dependent; the other methods systematically make the optimization to not make progress, pushing the posterior back to the prior. ?? reports the progression of the error rate and mnll. For both MNIST and CIFAR10 , i b l m places the parameters where the network can consistently deliver better performance both in terms of error rate and mnll throughout the entire learning procedure.

. . Comparison with variational inference beyond mean eld

Gaussian

Monte Carlo Dropout (Gal and Ghahramani, 2016a) offers a simple and effective way to perform approximate Bayesian c n n inference, thanks to the connection between dropout (Simonyan and Zisserman, 2014) and variational inference. In this experiment, we aim to compare and discuss benefits and disadvantages of using a Gaussian posterior approximation with respect to the Bernoulli approximation that characterizes mcds. For a fair comparison, we implemented the same lenet architecture and the same learning procedure as in Gal and Ghahramani (2016a) 1 , with the two convolutional layers having 20 and 50 filters, respectively. We also use the same learning rate policy λ × (1 + τ × iter) -p with τ = 0.0001, p = 0.75, Calibration Error represents a summary statistic of the calibration [START_REF] Naeini | Obtaining Well Calibrated Probabilities Using Bayesian Binning[END_REF] and it is defined as follows,

1
ECE = M m=1 |B m | n |acc(B m ) -conf(B m )| , (3.8)
where the test set is divided into M disjoint subsets {B 1 , . . . , B M }, each of them corresponding to a given interval of the predicted confidence. Figure 4.9 shows the reliability diagrams and the ece for AlexNet trained on CIFAR10.

Finally, we demonstrate that -provided with a sensible initialization -even simple factorized Gaussian posterior can achieve state-of-the-art performance on CIFAR10 with vgg2 , a large scale cnn [START_REF] Simonyan | Very Deep Convolutional Networks for Large-Scale Image Recognition[END_REF]. In this experiment, in addition to m c d, we also compare with n o i s yk f ac, an approximation of matrix-variate Gaussian posterior using noisy natural gradients introduced by G. initialization of this kind of model, the approximate posterior is drastically different from a spherical Gaussian prior and the variational objective is majorly dominated by the regularization term rather than the reconstruction likelihood. To deal with such issue, we propose and implement a simple policy to gradually include the kl term in the elbo. Given the generic expression for the elbo, we modify the lower bound as follow:

L e l b o (ν) = E q(θ;ν) log p(y | θ) -λKL [q(θ; ν) p(θ)] with λ = γ (1 + exp(-α(iter -β))) -1
. This way, we start the optimization of the e l b o with a low regularization, and progressively increase it throughout the optimization. Effectively, this is equivalent to start with a prior with very weak regularization strength and slowly moving towards a proper spherical Gaussian.

For the experiment on VGG16, we used α = 2 • 10 -3 , β = 2.5 • 10 4 and γ = 10 -1 .

The practice of down-playing the k l term was initially introduced for generative models like variational autoencoders (vaes) [START_REF] Bowman | Generating Sentences from a Continuous Space[END_REF][START_REF] Sønderby | Ladder Variational Autoencoders[END_REF], but now accepted by the Bayesian deep learning community as a quick remedy for poor performance.

.4

Before moving to some final remarks, let's quick summarize the empirical evidence up to this point. We prove that the role of initialization is severely under reported in the literature of variational inference for Bayesian neural network This work fills an important gap in the literature of Bayesian deep learning, that is how to effectively initialize variational parameters in vi. We proposed a novel way to do so, i b l m, which is based on an iterative layer-wise initialization based on Bayesian linear models. Through a series of experiments, including regression and classification with dnns and cnns, we demonstrated the ability of our approach to consistently initialize the optimization in a way that makes convergence faster than alternatives inspired from the state-of-the-art in loss minimization for deep learning.

E F F I C I E N T PA R A M E T E R I Z A T I O N S F O R V A R I A T I O N A L P O S T E R I O R S
Over-parameterized models, such as deep neural networks and convolutional neural networks, form a class of models that are routinely adopted in a wide variety of applications, and for which Bayesian inference is desirable but extremely challenging. Variational inference offers the tools to tackle this challenge in a scalable way and with some degree of flexibility on the approximation, but for over-parameterized models this is challenging due to the over-regularization property of the variational objective. Inspired by the literature on kernel methods, and in particular on structured approximations of distributions of random matrices, this chapter will discuss Walsh-Hadamard variational inference (w h v i), a new parameterization for variational inference which leverages Walsh-Hadamard-based factorization strategies to reduce the parameterization and accelerate computations.

Extensive theoretical and empirical analyses demonstrate that w h v i yields considerable speedups and model reductions compared to other techniques to carry out approximate inference for over-parameterized models, and ultimately show how advances in kernel methods can be translated into advances in approximate Bayesian inference for Deep Learning.

.1

Since its inception, variational inference (vi) [START_REF] Jordan | An Introduction to Variational Methods for Graphical Models[END_REF] has continuously gained popularity as a scalable and flexible approximate inference scheme for a variety of models for which exact Bayesian inference is intractable. Bayesian neural networks [START_REF] Mackay | Bayesian methods for backpropagation networks[END_REF][START_REF] Neal | Bayesian Learning for Neural Networks[END_REF] represent a good example of models for which inference is intractable, and for which v i -and approximate inference in general -is challenging due to the nontrivial form of the posterior distribution and the large dimensionality of the parameter space [START_REF] Graves | Practical Variational Inference for Neural Networks[END_REF]Gal and Ghahramani, 2016b). Recent advances in v i allow one to effectively deal with these issues in various ways. For instance, a flexible class of posterior approximations can be constructed using, e.g., normalizing flows (D. [START_REF] Rezende | Variational Inference with Normalizing Flows[END_REF],

whereas the need to operate with large parameter spaces has pushed the research in the direction of Bayesian compression [START_REF] Louizos | Bayesian Compression for Deep Learning[END_REF][START_REF] Molchanov | Variational Dropout Sparsifies Deep Neural Networks[END_REF].

We saw in the previous chapters that employing v i is notoriously challenging for Bayesian deep learning, and even more so for over-parameterized statistical models. Let's consider a supervised learning task with N input vectors and corresponding labels collected in X = {x 1 , . . . , x N } and y = {y 1 , . . . , y N }, respectively; furthermore, let's consider deep neural networks (d n ns) with weight matrices L) , likelihood p(y | θ, X), and prior p(θ). Following standard variational arguments, after introducing an approximation q(θ; ν) to the posterior p(W | y, X) it is possible to obtain a lower bound to the log-marginal likelihood log p(y | X) as follows:

θ = W (1) , . . . , W (
log p(y | X) E q(θ;ν) log p(y | θ, X) -KL [q(θ; ν) p(θ)] . (4.1)
The first term acts as a model fitting term, whereas the second one acts as a regularizer, penalizing solutions where the posterior is far away from the prior. It is easy to verify that the Kullback-Leibler (k l) term can be the dominant one in the objective for over-parameterized models. For example, a mean field posterior approximation turns the kl term into a sum of as many kl terms as the number of model parameters, say Q, which can dominate the overall objective when Q N.

As a result, the optimization focuses on keeping the approximate posterior close to the prior, disregarding the rather important model fitting term. This issue has been observed in a variety of deep models [START_REF] Bowman | Generating Sentences from a Continuous Space[END_REF], where it was proposed to gradually include the kl term throughout the optimization [START_REF] Bowman | Generating Sentences from a Continuous Space[END_REF][START_REF] Sønderby | Ladder Variational Autoencoders[END_REF] to scale up the model fitting term (A. G. Wilson and Izmailov, 2020a;[START_REF] Wenzel | How Good is the Bayes Posterior in Deep Neural Networks Really?[END_REF] or to improve the initialization of variational parameters [START_REF] Rossi | Good Initializations of Variational Bayes for Deep Models[END_REF]. Alternatively, other approximate inference methods for deep models with connections to vi have been proposed, notably Monte Carlo

.1 .1 Contributions

In this chapter we make the several contributions. We propose a novel strategy to cope with model over-parameterization when using variational inference, which is inspired by the literature on kernel methods. Our proposal is to reparameterize the variational posterior over model parameters by means of a structured decomposition based on random matrix theory [START_REF] Tropp | Improved Analysis of the subsampled Randomized Hadamard Transform[END_REF], which has inspired a number of fundamental contributions in the literature on approximations for kernel methods, such as Fastfood [START_REF] Le | Fastfood -Approximating Kernel Expansions in Loglinear Time[END_REF] and Orthogonal Random Features (ortogonal random features (o r f), (F. X. [START_REF] Yu | Orthogonal Random Features[END_REF]). The key operation within our proposal is the Walsh-Hadamard transform, and this is why we name our proposal Walsh-Hadamard Variational Inference (whvi). Finally, unlike mean field v i, w h v i induces a matrix-variate distribution to approximate the posterior over the weights, thus increasing flexibility at a log-linear cost in D instead of linear.

We can think of our proposal as a specific factorization of the weight matrix, so we can speculate that other tensor factorizations [START_REF] Novikov | Tensorizing Neural Networks[END_REF] of the weight matrix could equally yield such benefits. Our comparison against various matrix factorization alternatives, however, shows that w h v i is superior to other parameterizations that have the same complexity. Furthermore, while matrix-variate posterior approximations have been proposed in the literature of vi [START_REF] Louizos | Structured and Efficient Variational Deep Learning with Matrix Gaussian Posteriors[END_REF], this comes at the expense of increasing the complexity, while our proposal keeps the complexity to log-linear in D.

Through a wide range of experiments on dnns and convolutional neural networks (c n ns), we demonstrate that our approach enables the possibility to run variational inference on complex over-parameterized models, while being competitive with state-of-the-art alternatives. Ultimately, our proposal shows how advances in kernel methods can be instrumental in improving v i, much like previous works showed how kernel methods can improve, e.g., Markov chain Monte Carlo sampling [START_REF] Sejdinovic | Kernel Adaptive Metropolis-Hastings[END_REF][START_REF] Strathmann | Gradientfree Hamiltonian Monte Carlo with Efficient Kernel Exponential Families[END_REF] and statistical testing [START_REF] Gretton | A Kernel Statistical Test of Independence[END_REF][START_REF] Gretton | A Kernel Two-sample Test[END_REF][START_REF] Zaremba | B-test: A Non-parametric, Low Variance Kernel Two-sample Test[END_REF].

.2

whvi is inspired by a line of works that developed from random feature expansions for kernel machines [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF], which we briefly review here. A positive-definite kernel function κ(x i , x j ) induces a mapping φ(x), which can be infinite dimensional depending on the choice of κ(•, •). Among the large literature of scalable kernel machines, random feature expansion techniques aim at constructing a finite approximation to φ(•). Denoting by K the matrix with elements κ(x i , x j ) and by Φ the matrix with rows φ(x i ) , this approximation yields K ≈ Φ Φ and this allows one to turn kernel machines into linear models in this kernel-induced representation Φ. For many kernel functions [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF][START_REF] Cho | Kernel Methods for Deep Learning[END_REF], this approximation is built by applying a nonlinear transformation to a random projection XΩ, where Ω has entries N(ω ij |0, 1). If the matrix of training points X is N × D and we are aiming to construct D random features, that is Ω is D × D, this requires N times O(D 2 ) time, which can be prohibitive when D is large.

Fastfood [START_REF] Le | Fastfood -Approximating Kernel Expansions in Loglinear Time[END_REF] tackles the issue of large dimensional problems by replacing the matrix Ω with a random matrix for which the space complexity is reduced from O(D 2 ) to O(D) and time complexity of performing products with input vectors is reduced from O(D 2 ) to O(D log D). In Fastfood, the matrix Ω is replaced by

Ω ≈ SHGΠHB, (4.2)
where Π is a permutation matrix, H is the Walsh-Hadamard matrix, whereas G and B are diagonal random matrices with standard Normal and Rademacher ({±1}) distributions, respectively. The Walsh-Hadamard matrix is defined recursively starting from

H 2 = 1 1 1 -1
, and then

H 2D = H D H D H D -H D , (4.3)
which it could possibly be scaled by D -1/2 to make it orthonormal. The product Hx can be computed in O(D log D) time and O(1) space using the in-place version of the Fast Walsh-Hadamard Transform (fwht, [START_REF] Fino | Unified Matrix Treatment of the Fast Walsh-Hadamard Transform[END_REF]. S is also diagonal with i.i.d. entries, and it is chosen such that the elements of Ω obtained by this series of operations are approximately independent and follow a standard Normal (see [START_REF] Tropp | Improved Analysis of the subsampled Randomized Hadamard Transform[END_REF] for more details). Fastfood inspired a series of other works on kernel approximations , whereby Gaussian random matrices are approximated by a series of products between diagonal Rademacher and Walsh-Hadamard matrices (F. X. [START_REF] Yu | Orthogonal Random Features[END_REF][START_REF] Bojarski | Structured Adaptive and Random Spinners for Fast Machine Learning Computations[END_REF].

.3

-Fastfood and its variants yield cheap approximations to Gaussian random matrices with pseudo-independent entries, and zero mean and unit variance. The question we address here is whether we can use these types of approximations as cheap approximating distributions for v i. Let's start by considering a prior for the elements of the diagonal matrix G = diag(g) and a variational posterior q(g) =

posterior with some desirable properties as discussed next. Let W = W (l) ∈ R D×D be the weight matrix of a d nn at layer (l), and consider W ∼ q(W) s.t. W = S 1 Hdiag( g)HS 2 with g ∼ q(g).

(4.4)

. 3 .1 Statistical properties of the structure induced by

The choice of a Gaussian q(g) and the linearity of the operations in Equation (4.4) induce a parameterization of a matrix-variate Gaussian distribution for W, which is controlled by S 1 and S 2 if we assume that we can optimize these diagonal matrices.

For a generic D 1 × D 2 matrix-variate Gaussian distribution, we have

W ∼ MN(M, U, V) if and only if vect(W) ∼ N(vect(M), V ⊗ U), (4.5)
where

M ∈ R D 1 ×D 2 is the mean matrix and U ∈ R D 1 ×D 1 and V ∈ R D 2 ×D 2 are two
positive definite covariance matrices among rows and columns, and ⊗ denotes the Kronecker product. Matrix-variate Gaussian posteriors for variational inference have been introduced in Louizos and Welling (2016); however, assuming full covariance matrices U and V is memory and computationally intensive (quadratic and cubic in D, respectively). w h v i captures covariances across weights (see Figure 4.1), while keeping memory requirements linear in D and complexity loglinear in D, and we shall now see why.

In w h v i, as S 2 is diagonal, . . .

HS 2 = [v 1 , . . . , v D ] with v i = (S 2 ) i,i (
SHdiag(v d-1 )     (4.6)
This rewriting, shows that the choice of q(g) yields q(vect(W)) = N(Aµ, AΣA ), proving that w h v i assumes a matrix-variate distribution q(W), see Figure 4.1 for an illustration of this.

We now derive the parameters of the matrix-variate distribution q(W) = MN(M, U, V)

of the weight matrix W ∈ H D×D given by whvi, W = S 1 Hdiag( g)HS 2 with g ∼ N(µ, Σ).

(4.7)

The mean M = S 1 Hdiag(µ)HS 2 derives from the linearity of the expectation. The covariance matrices U and V are non-identifiable: for any scale factor s > 0, we have MN(M, U, V) equals MN(M, sU, 1 s V). Therefore, we constrain the parameters such that Tr(V) = 1. The covariance matrices verify (see e.g. Section 1 in the supplement of [START_REF] Ding | Dimension folding PCA and PFC for matrix-valued predictors[END_REF])

U = E q(W) (W -M)(W -M) V = 1 Tr(U) E q(W) (W -M) (W -M) .
Denoting by Σ 1/2 a root of Σ and considering p(ε) = N(0, I), we have

Note that the Walsh-Hadamard matrix H is symmetric, while S 1 and S 2 are diagonal.

U = E p(ε) S 1 Hdiag(Σ 1/2 ε)HS 2 2 Hdiag(Σ 1/2 ε)HS 1 . (4.8)
Let's define the matrix T 2 ∈ H D×D 2 where the i th row is the column-wise vectorization of the matrix (Σ

1/2 i,j (HS 2 ) i,j ) j,j D , T 1 =     vect Σ 1,: (HS 2 ) 1,: . . . vect Σ D,: (HS 2 ) d,:     .
(4.9)

We have that

(T 2 T 2 ) i,i = D j,j =1 Σ 1/2 i,j Σ 1/2 i ,j (HS 2 ) i,j (HS 2 ) i ,j = D j,j ,j =1 Σ 1/2 i,j (HS 2 ) i,j E p(ε) [ j j ]Σ 1/2 i ,j (HS 2 ) i ,j = D j =1 E p(ε)     D j=1 j Σ 1/2 i,j (HS 2 ) i,j     D j =1 j Σ 1/2 i ,j (HS 2 ) i ,j     = E p(ε) diag(Σ 1/2 ε)HS 2 2 Hdiag(Σ 1/2 ε) i,i
.

Using Equation (4.8), a root of U = U 1/2 U 1/2 can be found: for one choice of S 1 , S 2 and q(g). The empirical quantities are obtained with Monte Carlo samples.

U 1/2 = S 1 HT 2 . ( 4
Similarly for V, we have

V 1/2 = 1 Tr(U) S 2 HT 1 , with T 1 =     vect Σ 1,: (HS 1 ) 1,: . . . vect Σ D,: (HS 1 ) d,:     . (4.11)
The form of the moments for the matrix variate Gaussian are non-trivial and they can be confusing; to verify the correctness of these formulas we run a numerical study: we compute the mean and the two covariances analitically and empirically after sampling 15 000 times from Equation (4.4) and we report their . The matrix A in Equation (4.6) expresses the linear relationship between the weights W = S 1 Hdiag(g)HS 2 and the variational random vector g, i.e. vect(W) = Ag. Recall the definition of

A =     S 1 Hdiag(v 1 )
. . .

S 1 Hdiag(v D )     , with v i = (S 2 ) i,i (H) :,i . (4.12)
We show that a LQ-decomposition of A can be explicitly formulated as.

vect(W) = [s (2) i S 1 Hdiag(h i )] i=1,...,D g = LQg, (4.13) where h i is the i th column of H, L = diag((s (2) i s) i=1,...,D ), diag(s (1) ) = S 1 , diag(s (2) ) = S 2 , and Q = [Hdiag(h i )] i=1,.
..,D . Equation (4.13) derives directly from block matrix and vector operations. As L is clearly lower triangular (even diagonal), let us proof that Q has orthogonal columns. Defining the d × d matrix Q (i) = Hdiag(h i ), we have:

Q Q = D i=1 Q (i) Q (i) = D i=1 diag(h i )H Hdiag(h i ) = D i=1 diag(h 2 i ) = D i=1 D -1 I = I.
The mean of the structured matrix-variate posterior assumed by whvi can span a D-dimensional linear subspace within the whole D 2 -dimensional parameter space, and the orientation is controlled by the matrices S 1 and S 2 . In more details, this decomposition gives direct insight on the role of the Walsh-Hadamard transforms: with complexity D log(D), they perform fast rotations Q of vectors living in a space of dimension D into a space of dimension D 2 . Treated as parameters gathered in L, S 1 and S 2 control the orientation of the subspace by distortion of the canonical axes.

. 3 . 2 Reparameterizations in for stochastic optimization

The so-called reparameterization trick [START_REF] Kingma | Auto-Encoding Variational Bayes[END_REF]) is a standard way to make the variational lower bound in Equation (4.1) a deterministic function of the variational parameters, so as to be able to carry out gradient-based optimization despite the stochasticity of the objective. Considering input vectors h i to a given layer, an improvement over this approach is to consider the distribution of the product Wh i . This is also known as the local reparameterization trick [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF], and it reduces the variance of stochastic gradients in the optimization, thus improving convergence. The product Wh i follows the distribution N(m, AA ) (A. K. [START_REF] Gupta | Matrix variate distributions[END_REF], with m = S 1 Hdiag(µ)HS 2 h i , and

A = S 1 Hdiag(HS 2 h i )Σ 1/2 . (4.14)
A sample from this distribution can be efficiently computed thanks to the Walsh-Hadamard transform as:

W(µ)h i + W(Σ 1/2 ε)h i , with W a linear matrix-valued function W(u) = S 1 Hdiag(u)HS 2 .

.,

The choice of the parameterization of W in w h v i leaves space to several possible alternatives, which we compare in Figure 4.3. For all of them, G is learned variationally and the remaining diagonal S i (if any) are either optimized or treated variationally (Gaussian mean-field). Figure 4.3 shows the behavior of these alternatives when applied to a 2 × 64 network with rectified linear unit (ReLU) activations.

With the exception of the simple and highly constrained alternative GH, all parameterizations are converging quite easily and the comparison with m c d shows that indeed the proposed w h v i performs better both in terms of error rate and mnll. Furthermore, this result allows us to speculate that our proposal is the best.

w h v i is effectively imposing a factorization of W, where parameters are either optimized or treated variationally. Tensor decompositions for dnns and cnns have been proposed in [START_REF] Novikov | Tensorizing Neural Networks[END_REF]; here W is decomposed into k small matrices

(tensor cores), such that W = W 1 W 2 • • • W k , where each W i has dimensions r i-1 × r i (with r 1 = r k = D).
We adapt this idea to make a comparison with w h v i. In order to match the space and time complexity of w h v i, assuming {r i = R|∀i = 2, . . . , k -1}, we set: R ∝ log 2 D and K ∝ D (log 2 D) 2 . Also, to match the number of variational parameters, all internal cores (i = 2, . . . , k -1) are learned with fully factorized Gaussian posterior, while the remaining are optimized (see Table 4.1). Given the same asymptotic complexity, Figure 4.4 reports the results of this comparison again on a 2 hidden layer network. Not only whvi can reach better solutions in terms of test performance, but optimization is also faster. We speculate that this is attributed to the redundant variational parameterization induced by the tensor cores, which makes the optimization landscapes highly multi-modal, leading to slow convergence.

..1 Extensions

. For the sake of presentation, so far we have assumed W ∈ R D×D with D = 2 d , but we can easily extend whvi to handle parameters of any shape W ∈ R D out ×D in . One possibility is to use w h v i with a large D × D matrix with D = 2 d , such that a subset of its elements represent W.

Alternatively, a suitable value of d can be chosen so that W is a concatenation by M is the number of pseudo-data used to sample from a matrix Gaussian distribution, see [START_REF] Louizos | Structured and Efficient Variational Deep Learning with Matrix Gaussian Posteriors[END_REF].

Complexity Space Time

Mean field Gaussian

O(D 2 ) O(D 2 ) Gaussian matrix variate O(D 2 ) O(D 2 + M 3 ) Tensor factorization O(KR 2 ) O(R 2 ) w h v i O(D) O(D log D)
Algorithm 3: Setup dimensions for non-squared matrix row/column of square matrices of size D = 2 d , padding if necessary (Algorithm 3 shows this case). When one of the dimensions is equal to one so that the parameter matrix is a vector (W = w ∈ R D ), this latter approach is not ideal, as whvi would fall back on mean-field vi. whvi can be extended to handle these cases efficiently by reshaping the parameter vector into a matrix of size 2 d with suitable d, again by padding if necessary. Thanks to the reshaping, w h v i uses √ D parameters to model a posterior over D, and allows for computations in O( √ D log D) rather than D. This is possible by reshaping the vector that multiplies the weights in a similar way.

. Normalizing flows (D. [START_REF] Rezende | Variational Inference with Normalizing Flows[END_REF] are a family of parameterized distributions that allow for flexible approximations. In the general setting, consider a set of invertible, continuous and differentiable functions

f k : R D → R D with parameters λ k . Given θ 0 ∼ q 0 (θ 0 ), θ 0 is transformed with a chain of K flows to θ K = (f K • • • • • f 1 )(θ 0 ).
The variational lower bound slightly differs from Equation (4.1) to take into account the determinant of the Jacobian of the transformation, yielding a new variational objective as follows:

E q log p(y | θ, X) -KL [q 0 (θ 0 ) p(θ K )] + E q 0 (θ 0 ) K k=1 log det ∂f k (θ k-1 ; λ k ) ∂θ k-1 . (4.15)
Setting the initial distribution q 0 to a fully factorized Gaussian N(θ 0 | µ, σI) and assuming a Gaussian prior on the generated θ K , the k l term is still tractable.

The tranformation f is generally chosen to allow for fast computation of the determinant of the Jacobian. The parameters of the initial density q 0 as well as the flow parameters λ are optimized. In our case, we consider q K as a distribution over the elements of g. This approach increases the flexibility of the form of the variational posterior in w h v i, which is no longer Gaussian, while still capturing covariances across weights. This is obtained at the expense of losing the possibility of employing the local reparameterization trick. In the remaining part of the chapter, we will use planar flows (D. [START_REF] Rezende | Variational Inference with Normalizing Flows[END_REF]. Although this is a simple flow parameterization, a planar flow requires only O(D) parameters and thus it does not increase the time/space complexity of whvi.

.5

In this section we will provide experimental evaluations of our proposal, with experiments ranging from regression on classic benchmark datasets to image classification with large-scale convolutional neural networks. We will also comment on the computational efficiency and some potential limitation of our proposal.

.

The experiments on Bayesian d n n are run with the following setup.

For w h v i, we used a zero-mean prior over g with fully factorized covariance λI; λ = 10 -5 was chosen to obtain sensible variances in the output layer. It is possible to design a prior over g such that the prior on W has constant marginal variance and low correlations although empirical evaluations showed not to yield a significant improvement compared to the previous (simpler) choice. In the final implementation of whvi that we used in all experiments, S 1 and S 2 are optimized.

We used classic Gaussian likelihood with optimized noise variance for regression and softmax likelihood for classification. Training is performed for 500 steps with fixed noise variance and for other 50000 steps with optimized noise variance. Batch size is fixed to 64 and for estimation of the stochastic gradients we use one sample at train time and 64 samples at test time. We choose the Adam optimizer [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF] with exponential learning rate decay λ t+1 = λ 0 (1 + γt) -p , with λ 0 = 0.001, p = 0.3, γ = 0.0005 and t being the current iteration. Similar setup was also used for the Bayesian cnn experiment. The only differences are the batch size -increased to 256 -and the optimizer, which is run without learning rate decay.

4 . 5 67

-2 0 2 -1 0 1 2 G-WHVI (This work) -2 0 2 NF-WHVI (This work) -2 0 2 SGHMC -2 0 2 -1 0 1 2 MCD -2 0 2
Mean-eld Gaussian ) and with planar normalizing ow (10 ows for a total of 4141 param.), mean-eld Gussian (35k param.) and Monte Carlo dropout ( ) (17k param.). The two shaded areas represent the 95th and the 75th percentile of the predictions. As "ground truth", we also show the predictive posterior obtained by running Stochastic Gradient Hamiltonian Monte Carlo (

) on the same model (R < 1.05, [START_REF] Gelman | Bayesian Data Analysis[END_REF]).

. 5 .1 Toy example

We begin our experimental validation with a 1D-regression problem. We generated a 1D toy regression problem with 128 inputs sampled from U[-1, 2], and removed 20% inputs on a predefined interval; targets are noisy realizations of a random function (noise variance σ 2 = exp(-3)). We model these data using a d n n with 2 hidden layers of 128 features and cosine activations. We test four models: variational inference with mean-field Gaussian, Monte Carlo dropout (m c d, Gal and Ghahramani, 2016b) with dropout rate 0.4 and two variants of w h v isgw h v i with Gaussian posterior and n fw h v i with planar flows (10 planar flows). We also show the free form posterior obtained by running a Markov chain Monte Carlo (m c m c) algorithm, s g h m c in this case (T. [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF][START_REF] Springenberg | Bayesian Optimization with Robust Bayesian Neural Networks[END_REF], for several thousands steps. As Figure 4.5 shows, w h v is offers a sensible modeling of the uncertainty on the input domain, whereas mean-field and m c d seem to be slightly over-confident.

. 5 . 2 Empirical comparison on the UCI benchmark

We conduct now a series of comparisons with state-of-the-art v i schemes for Bayesian dnns; see the Supplement for the list of data sets used in the experiments. 5.12 (0.79) 8.01 (0.41) 8.21 (0.55) 4.70 (0.72) 3.20 (0.36) 3.41 (0.05) 3.56 (0.08) 3.17 (0.37) Energy 2.07 (0.11) 3.10 (0.14) 1.96 (0.28) 0.58 (0.07) 4.15 (0.15) 4.91 (0.09) 2.11 (0.12) 2.00 (0.60) Kin8nm 0.09 (0.00) 0.12 (0.00) 0.07 (0.00) 0.08 (0.00) -0.87 (0.02) -0.83 (0.02) -1.19 (0.04) -1.19 (0.04) Naval 0.30 (0.30) 0.01 (0.00) 0.00 (0.00) 0.01 (0.00) -1.00 (2.27) -6.23 (0.01) -6.52 (0.09) -6.25 (0.01) We compare w h v i with m c d and n n g (n o i s yk f ac, G. [START_REF] Zhang | Noisy Natural Gradient as Variational Inference[END_REF]. m c d draws on a formal connection between dropout and v i with Bernoulli-like posteriors, while the more recent n o i s yk f ac yields a matrix-variate Gaussian distribution using noisy natural gradients. To these baselines, we also add the comparison with mean field Gaussian (m f g). In w h v i, the last layer assumes a fully factorized Gaussian posterior.

Data is randomly divided into 90%/10% splits for training and testing eight times.

We standardize the input features x while keeping the targets y unnormalized. Differently from the experimental setup in Louizos and Welling (2016), G. Zhang et al.

(2018), and Hernandez-Lobato and R. Adams (2015), we use the same architecture regardless of the size of the dataset. Futhermore, to test the efficiency of w h v i in case of over-parameterized models, we set the network to have two hidden layers and 128 features with ReLU activations (as a reference, these models are ∼20 times bigger than the usual setup, which uses a single hidden layer with 50/100 units).

We report the test r m s e and the average predictive test negative log-likelihood (mnll) in Table 4.2. On the majority of the datasets, whvi outperforms mcd and n o i s yk f ac, while requiring less parameters to be optimized. Futhermore, we study how the test m n l l varies with the number of hidden units in a 2-layered network. As Figure 4.6 shows, w h v i and m c d both behave well while the other competitive methods struggle. Empirically, these results demonstrate the value of w h v i, which offers a competitive parameterization of a matrix-variate Gaussian posterior.

Being able to increase width and depth of a model without drastically increasing the number of variational parameters is also one of the competitive advantages of w h v i. Using classification datasets, we now study the effect of changing depth and width of the architecture. Figure 4.7 reports a summary of this comparison, where we considere two and three hidden layer networks with width from 64 up to 512. While reasoning about the role of overparameterization in deep and wide neural networks for increasing generalization is beyond the scope of this work, we observe that at test time increasing the number of hidden layers and the numbers of hidden features w h v i allows the model to deliver better performance, while avoiding catastrofic overfitting. This is not always the case for e.g. m c d, which with few exceptions has performance degradation the wider and deeper a model is.

We report the training and test curves for Mocap in Figure 4.8. While the analysis of the optimization problem is not a primary focus of this work, the behavior of the training evidence lower bound (elbo) reported shows stable convergence. Thanks to the structure of the weights matrices, widening and deepening the model is beneficial performance-wise and doesn't incur in convergence instability observed for other methods. Furthermore, with a non-mean field assumption of the weights, the e l b o is tighter than before and, in principle, it could be used as a suitable objective function for model selection.

. 5 . 3 Bayesian convolutional neural networks for image classi cation

We continue the experimental evaluation of whvi by analyzing its performance on c n ns. For this experiment, we replace all fully-connected layers in the c n n with the whvi parameterization, while the convolutional filters are treated variationally using mcd. In this setup, we fit vgg (Simonyan and Zisserman, 2014), alexnet [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] and resnet -18 (K. [START_REF] He | Deep Residual Learning for Image Recognition[END_REF] on CIFAR10 . Using whvi, we can reduce the number of parameters in the linear layers without affecting neither test performance nor calibration properties of the resulting model, as shown in Figure 4.9 and Table 4.3. For a l e x n e t and r e s n e t we also try our variant of w h v i with normalizing flows. Even though we lose the benefits of the local reparameterization, the higher flexibility of normalizing flows allows the model to obtain better test performance with respect to the Gaussian posterior. This can be improved even further using more complex families of normalizing flows (D. [START_REF] Rezende | Variational Inference with Normalizing Flows[END_REF][START_REF] Van Den Berg | Sylvester Normalizing Flows for Variational Inference[END_REF][START_REF] Kingma | Improved Variational Inference with Inverse Autoregressive Flow[END_REF]Louizos and Welling, 2017). With w h v i, a l e x n e t and its original ∼23.3m parameters is reduced to just ∼2.3m (9.9%) when using g-w h v i and to ∼2.4m (10.2%) with w h v i and 3 planar flows. . Convolution layers are linear layers applied on a series of patches extracted from the image and they can be implemented with matrix multiplication, once filters are reshaped in 2D. Given this interpretation of the convolutional layers, we also extended whvi for for this case. We observe though that in this case resulting models had too few parameters to obtain any interesting results. For a l e x n e t , we obtained a model with just 189k parameters, which corresponds to a sparsity of 99.2% with respect of the original model. As a reference, [START_REF] Wen | Learning Structured Sparsity in Deep Neural Networks[END_REF] was able to reach sparsity only up to 60% in the convolutional layers without impacting performance. To study this behavior in details, we take a simple cnn with two convolutional layers and one linear layer (Figure 4.10). We see that the combination of mcd and whvi performs very well in terms of convergence and test performance, while the use of w h v i on the convolutional filters brings an overall degradation of the performance. Interestingly, though, we also observe that m c d with the same number of parameters as for w h v i (referred to as lowrank m c d) performs even worse than the baseline: this once again confirms the parameterization of w h v i as an efficient alternative.

. 5 .Comments on computational e ciency

whvi builds his computational efficiency on the Fast Walsh-Hadamard Transform (f w h t), which allows to cut the complexity of a D-dimensional matrix-vector multiplication Hx from a naive O(D 2 ) to O(D log D), without generating and storing H. For our experimental evaluation, we implemented a batched version of the fwht in PyTorch in C++ and CUDA to leverage the full computational capabilities of modern graphic processing units (g p us). Figure 4.11a presents a timing profiling of our implementation versus the naive matmul (batch size of 512 samples and profiling repeated 1000 times). The breakeven point for the c p u implementation is in the neighborhood of 512/1024 features, while on g p u we see that f w h t is consistently faster. In Figure 4.11b we also compare the time required to sample and compute the predictions on the test set on two regression datasets as a function of the number of hidden units in a two-layer dnn. The workstation used is equipped with two Intel Xeon c p us, four NVIDIA Tesla P100 and 512 GB of RAM. Each experiment is carried out on a g p u fully dedicated to it. The n n g algorithm is implemented in Tensorflow 1 while the others are written in PyTorch. We made sure to fully exploit all parallelization opportunities in the competiting methods and ours; we believe that the timings are not severely affected by external factors other than the actual implementation of the algorithms. We also report the storage footprint (in Bytes) for several configuration of depth/width in Figure 4.12. To put things into prospective, the test m n l l is also reported. The Pareto analysis shows that for most of the cases w h v i shows good trade-off between number of parameters and modeling performance. Indeed, there are configurations for which mcd delivers better test mnll, but all of them are at the expeses of bigger models. using whvi and mcd and noisy natural gradients. The regularity of the algorithm for computing the f w h t and its reduced memory footprint result on an overall higher utilization of the gpu, 85% for whvi versus ∼ 70% for mcd. This translates into an increase of energy efficiency up to 33% w.r.t m c d, despite being 51%

Finally, in

faster. We speculate that the poor performance of n n g is due to the inversion of the approximation to the Fisher matrix, which scales cubically in the number of units.

. 5 . 5 Exploring the parameter e ciency of WHVI with Gaussian processes

We conclude this empirical evaluation with a discussion on a possible application of whvi for scalable Gaussian processs (gps). gps represent the perfect controlled environment where we can analyze the trade off between model approximation and inference approximation, with a focus on parameter efficiency. To this extent, we focus on g ps with random feature expansions [START_REF] Lázaro-Gredilla | Sparse Spectrum Gaussian Process Regression[END_REF].

Given a random feature expansion of the kernel martix, say K ≈ ΦΦ , the latent variables can be rewritten as f = Φw, with w ∼ N(0, I). The random features Φ are constructed by randomly projecting the input matrix X using a Gaussian random matrix Ω and applying a nonlinear transformation, which depends on the choice of the kernel function. For a given set of random features, we can treat the weights of the resulting model variationally, with e.g. a mean field Gaussian q(w). By reshaping the vector of parameters w into a D × D matrix, w h v i allows for the linearized g p to reduce the number of parameters to optimize or, given the same complexity, increase the number of random features (see ??). This latter case will be also compared with another approximation of Gaussian process using variational inference and inducing points [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF][START_REF] Hensman | MCMC for Variationally Sparse Gaussian Processes[END_REF]. We report the results on four datasets in Figure 4.14. This variant of w h v i that reshapes w into a matrix largely outperforms the direct variational inference with mean field Gaussian. However, it appears that this improvement of the random feature inference for gps is not enough to reach the performance of inducing points. Inducing point approximations are based on the Nyström approximation of kernel matrices, which are known to lead to lower approximation error on the elements on the kernel matrix compared to random features approximations. This is the The task is repeated 16 consecutive times and pro ling is carried out using the nvidia-smi tool. reason we attribute to the lower performance of whvi compared to inducing points approximations in this experiment.

.6

In the early sections of the chapter, we have already briefly reviewed some of the literature on vi and Bayesian dnns and cnns; here we complement the literature by including other relevant works that have connections with whvi.

Our work takes inspiration from the works on random features for kernel approximation [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF] and Fastfood [START_REF] Le | Fastfood -Approximating Kernel Expansions in Loglinear Time[END_REF]. Random feature expansions have had a wide impact on the literature on kernel methods. Such approximations have been successfully used to scale a variety of models, such as Support Vector Machines [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF], Gaussian processes [START_REF] Lázaro-Gredilla | Sparse Spectrum Gaussian Process Regression[END_REF] and Deep Gaussian processes [START_REF] Cutajar | Random Feature Expansions for Deep Gaussian Processes[END_REF]Gal and Ghahramani, 2016b). This has contributed to bridging the gap between Deep g ps and Bayesian d n ns and c n ns [START_REF] Neal | Bayesian Learning for Neural Networks[END_REF][START_REF] Duvenaud | Avoiding pathologies in very deep networks[END_REF][START_REF] Cutajar | Random Feature Expansions for Deep Gaussian Processes[END_REF]Gal and Ghahramani, 2016a), which is an active area of research which aims to gain a better understanding of deep learning models through the use of kernel methods (G. [START_REF] Matthews | Gaussian Process Behaviour in Wide Deep Neural Networks[END_REF][START_REF] Dunlop | How Deep Are Deep Gaussian Processes?[END_REF][START_REF] Garriga-Alonso | Deep Convolutional Networks as shallow Gaussian Processes[END_REF]. Structured random features [START_REF] Le | Fastfood -Approximating Kernel Expansions in Loglinear Time[END_REF][START_REF] Yu | Orthogonal Random Features[END_REF][START_REF] Bojarski | Structured Adaptive and Random Spinners for Fast Machine Learning Computations[END_REF] have been also applied to the problem of handling large dimensional convolutional features (Z. [START_REF] Yang | Deep Fried Convnets[END_REF] and Convolutional g ps (G.-L. [START_REF] Tran | Calibrating Deep Convolutional Gaussian Processes[END_REF].

Bayesian inference on dnns and cnns has been research topic of several seminar works (see e.g. [START_REF] Graves | Practical Variational Inference for Neural Networks[END_REF][START_REF] Hernandez-Lobato | Probabilistic Backpropagation for Scalable Learning of Bayesian Neural Networks[END_REF][START_REF] Blundell | Weight Uncertainty in Neural Network[END_REF]Gal and Ghahramani, 2016b;Gal and Ghahramani, 2016a). Recent advances in d n ns have investigated the effect of over-parameterization and how model compression can be used during or after training [START_REF] Hubara | Binarized Neural Networks[END_REF][START_REF] Louizos | Bayesian Compression for Deep Learning[END_REF][START_REF] Zhu | To Prune, or Not to Prune: Exploring the Efficacy of Pruning for Model Compression[END_REF]. Our current understanding shows that model performance is affected by the network size with bigger and wider neural networks being more resilient to overfit [START_REF] Neyshabur | In Search of the Real Inductive Bias: On the Role of Implicit Regularization in Deep Learning[END_REF][START_REF] Neyshabur | The role of overparametrization in generalization of neural networks[END_REF]. For variational inference, and Bayesian inference in general, over-parameterization is reflected on over-regularization of the objective, leading the optimization to converge to trivial solutions (posterior equal to prior). Several works have encountered and proposed solutions to such issue [START_REF] Higgins | beta-VAE: Learning Basic Visual Concepts with a Constrained Variational Framework[END_REF][START_REF] Burgess | Understanding disentangling in β-VAE[END_REF][START_REF] Bowman | Generating Sentences from a Continuous Space[END_REF][START_REF] Sønderby | Ladder Variational Autoencoders[END_REF][START_REF] Rossi | Good Initializations of Variational Bayes for Deep Models[END_REF]. The problem of how to run accurate Bayesian inference on over-parametrized models like Bayesian neural networks (b n ns) is still an ongoing open question (A. G. Wilson and Izmailov, 2020a;[START_REF] Wenzel | How Good is the Bayes Posterior in Deep Neural Networks Really?[END_REF]

.7

Inspired by the literature on scalable kernel methods, this chapter proposed Walsh-Hadamard Variational Inference (w h v i). w h v i offers a novel parameterization of the variational posterior, which is particularly attractive for over-parameterized models, such as modern dnns and cnns. whvi assumes a matrix-variate posterior distribution, which therefore captures covariances across weights. Crucially, unlike 4 . 7 77 previous work on matrix-variate posteriors for v i, this is achieved with a light parameterization and fast computations, bypassing the over-regularization issues of vi for over-parameterized models. The large experimental campaign, demonstrates that whvi is a strong competitor of other variational approaches for such models, while offering considerable speedups. One limitation of the parameterization induced by whvi is that its mean cannot span the whole D 2 -dimensional space. A simple and cheap remedy could be to modify the parameterization from S 1 HGHS 2 to S 1 HGHS 2 + M with E[G] = 0 so that the mean can span the whole space thanks to M, while the rest would allow to model covariances across weights. However, though a preliminary numerical study on the rmse between the weights induced by whvi and arbitrary weight matrices, we see constant behavior w.r.t. D. Possible extensions include the possibility to capture the covariance between weights across layers, by either sharing the matrix G across, or by concatenating all weights into a single matrix which is then treated using whvi, with the necessary adaptations to handle the sequential nature of computations.

T H E E F F E C T O F S E L E C T I N G T H E P R I O R F O R B A Y E S I A N D E E P L E A R N I N G

The Bayesian treatment of neural networks dictates that a prior distribution is specified over their weight and bias parameters. This poses a challenge because modern neural networks are characterized by a large number of parameters, and the choice of these priors has an uncontrolled effect on the induced functional prior, which is the distribution of the functions obtained by sampling the parameters from their prior distribution. We argue that this is a hugely limiting aspect of Bayesian deep learning, and this chapter tackles this limitation in a practical and effective way.

.1

The concept of prior distribution in Bayesian inference allows us to describe the family of solutions that we consider acceptable, before having seen any data.

While in some cases selecting an appropriate prior is easy or intuitive given the context [START_REF] O'hagan | Bayes-Hermite quadrature[END_REF][START_REF] Rasmussen | Bayesian Monte Carlo[END_REF][START_REF] Srinivas | Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design[END_REF][START_REF] Cockayne | A Bayesian Conjugate Gradient Method (with Discussion)[END_REF][START_REF] Briol | Probabilistic Integration: A Role in Statistical Computation?[END_REF], for nonlinear parametric models with thousands (or millions) of parameters, like deep neural networks (d n ns) and convolutional neural networks (c n ns), this choice is not straightforward. Despite many advances in the field [START_REF] Kendall | What Uncertainties Do We Need in Bayesian Deep Learning for Computer Vision?[END_REF][START_REF] Rossi | Good Initializations of Variational Bayes for Deep Models[END_REF][START_REF] Osawa | Practical Deep Learning with Bayesian Principles[END_REF][START_REF] Rossi | Walsh-Hadamard Variational Inference for Bayesian Deep Learning[END_REF], it is reported that in some cases the predictive posteriors are not competitive to non-Bayesian alternatives, making these models-and Bayesian deep learning, in general-less than ideal solutions for a number of applications. [START_REF] Wenzel | How Good is the Bayes Posterior in Deep Neural Networks Really?[END_REF], for example, raise concerns about the quality and the usefulness of Bayesian neural network (b n n) posteriors, where it is found that tempering the posterior distribution improves the performance of some deep models. Observations of this kind should not be really surprising.

Bayesian inference is a recipe with exactly three ingredients: the prior distribution, the likelihood and the Bayes' rule. Regarding the Bayes' rule, that is simply a consequence of the axioms of probability. The fact that the posterior might not be useful in some cases should never be attributed to the Bayesian method itself. In fact, it is very easy to construct Bayesian models with poor priors and/or likelihoods, which result in poor predictive posteriors. One should therefore turn to the other two components, which encode model assumptions.

In this chapter we will discuss the role of the priors in the predictive performance of bnns. The common practice in the literature is to define a prior distribution on the network weights and biases, which is often chosen to be Gaussian. We argue that this choice can be catastrophic if miss-placed; nonetheless even with a simple Gaussian there are optimal configurations. In Figure 5.1, for example, A. G. Wilson and Izmailov (2020a) show the effect of the prior variance for two deep cnns. Note that all weights and biases in all the layers share the same parameter, which-as we will see in a moment-might not be the best solution. From this analysis we can draw two conclusions: (i) priors have huge effect in the modeling task, (ii) there are configurations which are much better than the simple N(0, 1) priors but for which we don't grasp the intuition.

While these argument seams only logical, in practice they are not. One might expect, for instance, that provided with a sufficiently large amount of data, the effect of the prior distribution will decrease. While this is true, the problem with Bayesian deep learning stands with the huge parameters dimensionality of the models considered, which pushes this limit further away. We show a simple case to exemplify this behaviour. In Figure 5.2 we consider a very simple LeNet classifier which is trained with variational inference and mean-field Gaussian posteriors.

For different combination of prior variance for linear and convolutional layers we report the "converged" test error and test likelihood. As can be seen, it appears that the effect of the prior is more critical for linear layers rather then convolutional ones, which reflect the allocation of trainable parameters. This justify a careful analysis of the effect of the priors in deep convolutional and feedforward neural networks.

.1 .1 Pathologies of deep prior functions

A prior over the parameters induces a prior on the functions generated by the model, which also depends on the network architecture. However, due to the nonlinear nature of the model, the effect of this prior on the functional output is not obvious to characterize and control. Consider the example in Figure 5.3, where we show the functions generated by sampling the weights of bnns with tanh activation from their Gaussian prior N(0, 1). We see that as depth is increased, the samples tend to form straight horizontal lines, which is a well-known pathology stemming from increasing model's depth [START_REF] Neal | Bayesian Learning for Neural Networks[END_REF][START_REF] Duvenaud | Avoiding pathologies in very deep networks[END_REF][START_REF] Matthews | Gaussian Process Behaviour in Wide Deep Neural Networks[END_REF]. We stress that a fixed Gaussian prior on the parameters is not always problematic, but it can be, especially for deeper architectures. More recently, issues of these kind of b n n priors have been recently exposed by Wenzel et al., 2020, who suggest to consider a temperate version of the posterior-effectively reducing the strength of the prior-to increase performance. Other recent works (T. [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF][START_REF] Springenberg | Bayesian Optimization with Robust Bayesian Neural Networks[END_REF] consider a hierarchical structure for the prior, where the variance of the normally-distributed b n n weights is governed by a Gamma distribution. This setting introduces additional flexibility on the space of functions, but it still does not provide much intuition regarding the properties of the prior.

Bayesian model selection constitutes a principled approach to select an appropriate prior distribution. Model selection is based on the marginal likelihood -the normalizing constant of the posterior distribution -which may be estimated from the training data. This practice is usually used to select hyperparameters of a Gaussian process (gp) as its marginal likelihood is available in closed form [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]. However, the marginal likelihood of bnns is generally intractable, and lower bounds are difficult to obtain. Graves (2011) first and[START_REF] Blundell | Weight Uncertainty in Neural Network[END_REF] later used the variational lower bound of the marginal likelihood for optimizing the parameters of a prior, yielding in some cases worse results. The Laplace's method is another alternative to approximate the marginal likelihood (Immer, Bauer, et al., 2021;[START_REF] Mackay | Probable Networks and Plausible Predictions -a Review of Practical Bayesian Methods for Supervised Neural Networks[END_REF] which is scalable and differentiable with respect to the prior hyperparameters.

Nonetheless, this kind of generative priors on the functions is very different from shallow Bayesian models, such as gps, where the selection of an appropriate prior typically reflects certain attributes that we expect from the generated functions. A gp defines a distribution of functions that is characterized by a mean and a kernel function κ. The gp prior specification can be more interpretable than the one induced by the prior over the weights of a bnn, in the sense that the kernel effectively governs the properties of possible functions, such as shape, variability and smoothness.

For example, shift-invariant kernels may impose a certain characteristic length-scale on the prior distribution over functions.

.1 . 2 Contributions

The main research question that we investigate now is how to impose functional priors on b n ns. We seek to tune the prior distributions over b n ns parameters so that the induced functional priors exhibit interpretable properties, similar to shallow g ps. While b n n priors induce a regularization effect that penalizes large values for the network weights, a gp-adjusted prior induces regularization directly on the space of functions.

We consider the Wasserstein distance between the distribution of b n n functions induced by a prior over their parameters, and a target g p prior. We propose an algorithm that optimizes such a distance with respect to the bnn prior parameters and hyper-parameters. An attractive property of our proposal is that estimating the Wasserstein distance relies exclusively on samples from both distributions, which are easy to generate. We explore the effect of g p-induced priors on the predictive posterior distribution of bnns by carrying out fully Bayesian inference of neural network models with these priors through the use of scalable Markov chain Monte Carlo (m c m c) sampling (T. [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF]. We demonstrate systematic performance improvements over alternative choices of priors and state-of-the-art approximate Bayesian deep learning approaches, including convolutional neural networks.

.2

Many recent attempts in the literature have turned their attention towards defining priors in the space of functions, rather than the space of weights. For example, [START_REF] Nalisnick | Predictive Complexity Priors[END_REF] consider a family of priors that penalize the complexity of predictive functions. [START_REF] Hafner | Noise Contrastive Priors for Functional Uncertainty[END_REF] propose a prior that is imposed on training inputs, as well as out-of-distribution inputs. This is achieved by creating pseudodata by means of perturbing the training inputs; the posterior is approximated by a variational scheme. W. [START_REF] Yang | Output-Constrained Bayesian Neural Networks[END_REF] present a methodology to induce prior knowledge by specifying certain constraints on the network output. [START_REF] Pearce | Expressive Priors in Bayesian Neural Networks: Kernel Combinations and Periodic Functions[END_REF] explore dnn architectures that recreate the effect of certain kernel combinations for gps. This result in an expressive family of network priors that converge to gps in the infinite-width limit. A different approach is proposed by [START_REF] Karaletsos | Gaussian Process Meta-Representations For Hierarchical Neural Network Weight Priors[END_REF] and [START_REF] Karaletsos | Hierarchical Gaussian Process Priors for Bayesian Neural Network Weights[END_REF], who consider a g p model for the network parameters that can capture weight correlations.

A similar direction of research focuses not only on priors but also inference in the space of functions for b n ns. For example, [START_REF] Ma | Variational Implicit Processes[END_REF] consider a b n n as an implicit prior in function space and then use gp for inference. Conversely, Sun, G. [START_REF] Sun | Functional Variational Bayesian Neural Networks[END_REF] propose a functional variational inference which employs a gp prior to regularize bnns directly in the function space by estimating the Kullback-Leibler (k l) divergence between these two stochastic processes. However, this method relies on a gradient estimator which can be inaccurate in high dimensions. In another route, M. E. E. [START_REF] Khan | Approximate Inference Turns Deep Networks into Gaussian Processes[END_REF] derive a g p posterior approximation for neural networks by means of the Laplace and generalized Gauss-Newton (ggn) approximations, leading to an implicit linearization. [START_REF] Immer | Improving predictions of Bayesian neural nets via local linearization[END_REF] make this linearization explicit and apply it to improve the performance of b n n predictions. In general, these approaches either rely heavily on non-standard inference methods or are constrained to use a certain approximate inference algorithm such as variational inference or Laplace approximation.

A different line of work focuses on meta-learning by adjusting priors based on the performance of previous tasks [START_REF] Amit | Meta-Learning by Adjusting Priors Based on Extended PAC-Bayes Theory[END_REF]. In contrast to these approaches, we aim to define a suitable prior distribution entirely a priori. We acknowledge that our choice to impose g p (or hierarchical g p) priors on neural networks is essentially heuristic: there is no particular theory that necessarily claims superiority for this kind of prior distributions. In some applications, it could be preferable to use priors that are tailored to certain kinds of data or architectures, such the deep weight prior [START_REF] Atanov | The Deep Weight Prior[END_REF]. However, we are encouraged by the empirical success and the interpretability of gp models, and we seek to investigate their suitability as b n n priors on a wide range of regression and classification problems.

. 3 85

Our work is most closely related to a family of works that attempt to map g p priors to b n ns. [START_REF] Flam-Shepherd | Mapping Gaussian Process Priors to Bayesian Neural Networks[END_REF] propose to minimize the k l between the bnn prior and some desired gp. As there is no analytical form for this kl, the authors rely on approximations based on moment matching and projections on the observation space. This limitation was later addressed [START_REF] Flam-Shepherd | Characterizing and Warping the Function space of Bayesian Neural Networks[END_REF] by means of a hypernetwork [START_REF] Ha | HyperNetworks[END_REF], which generates the weight parameters of the original b n n; the hypernetwork parameters were trained so that b n n fit the samples of a g p. In our work, we also pursue to minimize a sample-based distance between the bnn prior and some desired gp, but we avoid the difficulties in working with the kl divergence, as its evaluation is challenging due to the empirical entropy term. To the best of our knowledge, the Wasserstein distance scheme we propose is novel, and it demonstrates satisfactory convergence for compatible classes of g ps and b nns.

.3

The equivalence between function-space view and weight-space view of linear models, like Bayesian linear regression and gps [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF], is a straightforward application of Gaussian identities, but it allows us to seamlessly switch point of view accordingly to which characteristics of the model we are willing to observe or impose. We would like to leverage this equivalence also for bnns but the nonlinear nature of such models makes it analytically intractable (or impossible, for non-invertible activation functions). We argue that for bnns-and Bayesian deep learning models, in general-starting from a prior over the weights is not ideal, given the impossibility of interpreting its effect on the family of functions that the model can represent. We therefore rely on an optimization-based procedure to impose functional priors on bnns using the Wasserstein distance as a similarity metric between such distributions, as described next.

. 3 .1 A quick introduction to the Wasserstein distance

The concept of distance between probability measures is central to this work, as we frame the problem of imposing a g p prior on a b n n as a distance minimization problem. We present some known results on the Wasserstein distance that will be used in the sections that follow. Given two Borel's probability measures π(x) and ν(y) defined on the Polish space X and Y (i.e. any complete separable metric space), the generic formulation of the p-Wasserstein distance is defined as follows:

W p (π, ν) = inf γ∈Γ (π,ν) X×Y D(x, y) p γ(x, y) dx dy 1/p , (5.1)
where D(x, y) is a proper distance metric between two points x and y in the space X × Y and Γ (π, ν) is the set of functionals of all possible joint densities γ whose marginals are π and ν.

When the spaces of x and y coincide (i.e. x, y ∈ X ⊆ R d ), with D(x, y) being the Euclidian norm distance, the Wasserstein-1 distance (also known in the literature as Earth-Mover distance) takes the following shape,

W 1 (π, ν) = inf γ∈Γ (π,ν) X×X
xy γ(x, y) dx dy .

(5.2)

With the exception of few cases where the solution is available analytically (e.g.

π and ν being Gaussians), solving Equation (5.2) directly or via optimization is intractable. On the other hand, the Wasserstein distance defined in Equation ( 5.2) admits the following dual form [START_REF] Kantorovich | On the transfer of masses[END_REF][START_REF] Kantorovich | On a problem of Monge[END_REF],

W 1 (π, ν) = sup φ L 1 φ(x)π(x)dx -φ(y)ν(y)dy = sup φ L 1 E π φ(x) -E ν φ(x) , (5.3)
where φ is a 1-Lipschitz continuous function defined on X → R. This is effectively a functional maximization over φ on the difference two expectations of φ under π and ν. A revised proof of this dual form by [START_REF] Villani | Topics in Optimal Transportation[END_REF] is available in the Appendix.

Following recent literature [START_REF] Goodfellow | Generative Adversarial Nets[END_REF][START_REF] Arjovsky | Wasserstein Generative Adversarial Networks[END_REF], we parameterize the Lipschitz function by a neural network with parameters θ. In this setup, the supremum is replaced with a maximization problem with respect to the parameters θ. In order to enforce the Lipschitz constraint on φ θ , [START_REF] Arjovsky | Wasserstein Generative Adversarial Networks[END_REF] propose to clip the weights θ to lie within a compact space [-c, c] such that all function φ θ (•) will be K-Lipschitz, with K usually unknown. This approach usually biases the resulting φ θ (•) towards simple functions. Based on the fact that a differentiable function is 1-Lipschitz if and only if the norm of its gradient is 5 . 3 at most one everywhere, [START_REF] Gulrajani | Improved Training of Wasserstein GANs[END_REF] propose to constrain the gradient norm of the output of the Lipschitz function φ θ (•) with respect to its input. More specifically, the loss of the Lipschitz function is augmented by a regularization term

max θ E π(x) φ(x; θ) -E ν(x) φ(x; θ) + λE p( x) ∇ x φ( x) 2 -1 2 (5.4) Here p( x) is the distribution of x = ε x π + (1 -ε) x ν for ε ∼ U[0, 1] and x π ∼ π(x),
x ν ∼ ν(x) being the sample functions from the two marginals and λ the penalty coefficient.

. 3 . 2 Using the Wasserstein distance to impose the GP behaviour in

BNN

Let now go back to our original problem of optimizing priors for bnns. Assume a prior distribution p(w; ψ) on the weights of a bnn, where ψ is a set of parameters that determine the prior (e.g. for a Gaussian prior, ψ = {µ, σ}; we discuss more options on the parametrization of b n n priors in the section that follows). This prior over weights induces a prior distribution over functions: (5.5) where p(f | w) is deterministically defined by the network architecture.

p nn (f; ψ) = p(f | w)p(w; ψ) dw,
We consider as our target distribution a g p defined as p gp (f | 0, K), where K is the covariance matrix obtained by computing the kernel function κ for each pair of datapoints {x i , x j }. We aim at matching these two stochastic processes at a finite number of measurement points

X M def = [x 1 , ..., x M ] .
To achieve this, we propose a samplebased approach using the 1-Wasserstein distance as objective:

min ψ max θ E p gp [φ θ (f M )] -E p nn [φ θ (f M )] + λE p f ∇ fφ( f) 2 -1 2 . (5.6)
where f M denotes the set of random variables associated with the inputs at X M , and φ θ is 1-Lipschitz function.

Regarding the optimization of the θ and ψ parameters we alternate between n Lipschitz steps of maximizing L with respect to the Lipschitz function's parameters θ and one step of minimizing the Wasserstein distance with respect to the prior's parameters ψ. We therefore use two independent optimizers for θ and ψ. Figure 5.4 offers a high-level schematic representation of the proposed procedure. Given samples from two stochastic processes, the Wasserstein distance is estimated by considering the inner maximization of Equation (5.6), resulting in an optimal φ * . This inner optimization step is repeated for every step of the outer optimization loop. Notice that the objective is fully sample-based. As a result, it is not necessary to know the closed-form of the marginal density p nn (f; ψ). One may consider any stochastic process as a target prior over functions, as long as we can draw samples from it.

. 3 . 3 Prior Parameterization for Neural Networks

In the previous section, we have treated the parameters of a b n n prior p nn (f; ψ) in a rather abstract manner. Now we explore three different parametrizations of increasing complexity. The only two requirements needed to design a new parametrization are (1) to be able to generate samples and (2) to compute the log-density at any point; the latter is required to be able to draw samples from the posterior over model parameters using most m c m c sampling methods, such as Stochastic Gradient Hamiltonian Monte Carlo (sghmc). be analytically computed by taking into account to Jacobian of the transformations as follows,

log p(T K (w l )) = log p 0 (w l ) - K k=1 log det ∂t k (w l k-1 ) ∂w l k-1 , (5.7)
where

w l k-1 = (t k-1 • ... • t 2 • t 1 )(w l ) for k > 1,
and w l 0 = w l .

We shall refer to this class of bnn priors as the gp-induced bnn prior, parametrized by normalizing flows (gpi-nf). We note that nfs are typically used differently in the literature; while previous works showed how to use this distributions for better approximation of the posterior in variational inference (D. [START_REF] Rezende | Variational Inference with Normalizing Flows[END_REF][START_REF] Kingma | Improved Variational Inference with Inverse Autoregressive Flow[END_REF]Louizos and Welling, 2017) or for parametric density estimation (e.g. [START_REF] Grover | Flow-GAN: Combining Maximum Likelihood and Adversarial Learning in Generative Models[END_REF], or for enlarging the flexibility of a prior for variational autoencoders (vaes) (e.g. X. [START_REF] Chen | Variational Lossy Autoencoder[END_REF], as far as we are aware this is the first time that nfs are used to characterize a prior distribution for b n ns.

A simple setup consists in setting the initial distribution p 0 (w l ) to a fully-factorized

Gaussian N(w l | 0, σ 2 l I) and then employing a sequence of planar flows (D. [START_REF] Rezende | Variational Inference with Normalizing Flows[END_REF], each defined as

t k (w l k-1 ) = w l k-1 + u l k h(θ l k w l k-1 + b l k ),
(5.8)

where u l k ∈ R D l , θ l k ∈ R D l , b l k ∈ R are trainable parameters, and h(•) = tanh(•). The log-determinant of the Jacobian of t k is log det ∂t k (w l k-1 ) ∂w l k-1 = log 1 + u l k θ l k h (θ l k w l k-1 + b l k ) .
(5.9)

Thus for the l-th bnn layer, the parameters to optimize are

ψ l = {σ 2 l } {u l k , θ l k , b l k } K k=1 .

.4

In the following experiments, we consider two fixed priors: (1) fixed Gaussian (f g) prior, N(0, 1);

(2) fixed hierarchical (f h) prior where the prior variance for each layer is sampled from an Inverse-Gamma distribution, Γ -1 (1, 1) [START_REF] Springenberg | Bayesian Optimization with Robust Bayesian Neural Networks[END_REF]; and two g p-induced neural network (n n) priors, namely: (3) g pinduced Gaussian (g pi-g) prior, (4) g p-induced hierarchical (g pi-h) prior. Since T; σ 2 , and ( α, β) denote the optimized parameters for the i-and i-priors, respectively. κ corresponds to optimized kernel parameters, while σ 2 LA shows that the parameters are optimized on the Laplace approximation of the marginal likelihood. References are [a] for Wenzel et al. (2020), [b] for Springenberg et al. (2016), [c] for [START_REF] Lakshminarayanan | Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles[END_REF].

Priors Inference Name p(σ 2 ) p(w | σ 2 ) p(f) Reference Fixed Gaussian ( ) prior - N(0, σ 2 I) → ? Fixed hierarchical ( ) prior Γ -1 (α, β) → N(0, σ 2 I) → ? + Gibbs [b]
Fixed Gaussian prior and TS ( + ) -

N(0, σ 2 I) → ? Tempered [a] Deep ensemble - ? ? Ensemble [c] GP-induced Gaussian ( i-) prior - N(0, σ 2 I) ← GP(0, κ) GP-induced hierarchical ( i-) prior Γ -1 ( α, β) ← N(0, σ 2 I) ← GP(0, κ) + Gibbs GP-induced norm. ow ( i-) prior - NF(T K ) ← GP(0, κ)
the computational cost of the gpi-nf prior is high, we don't consider this prior for large scale models.

We additionally compare b n ns against Deep Ensemble [START_REF] Lakshminarayanan | Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles[END_REF], arguably one of the state-of-the-art approaches for uncertainty estimation in deep learning [START_REF] Ashukha | Pitfalls of In-Domain Uncertainty Estimation and Ensembling in Deep Learning[END_REF][START_REF] Ovadia | Can you trust your model's uncertainty? Evaluating predictive uncertainty under dataset shift[END_REF]. This non-Bayesian method combines solutions that maximize the predictive log-likelihood for multiple neural networks trained with different initializations. We employ an ensemble of 5 neural networks in all experiments. Following [START_REF] Lakshminarayanan | Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles[END_REF], we use Adam optimizer [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF] to train the individual networks.

Furthermore, we compare the results obtained by sampling from the posterior obtained with g p-induced priors against "tempered" posterior [START_REF] Wenzel | How Good is the Bayes Posterior in Deep Neural Networks Really?[END_REF] that uses the f g prior and temperature scaling; we refer to this approach as f g+ t s.

Finally, Table 6.1 summaries an overview of methods considered in the experiments. ). The rst two rows illustrate the prior sample and distributions, whereas the last two rows show the corresponding posterior distributions. The means and the 95% credible intervals are represented by red lines and shaded areas, respectively. The middle row shows progressions of the prior optimization.

.4 .1 Visualization on a 1D regression synthetic dataset

The dataset used is built as follows: (1) we uniformly sample 64 input locations

x in the interval [-10, 10] ; (2) we rearrange the locations on a defined interval (3) we sample a function f from the g p prior (l = 0.6, α = 1) computed at locations x; (4) we corrupt the targets with i.i.d.

Gaussian noise (σ 2 = 0.1). In this example, we consider a three hidden-layer m l p. Figure 5.5 shows all the results. The first two rows illustrate the different choice of priors. For the Wasserstein-based functional priors (gpi-g, gpi-h, gpi-nf), the third row shows the convergence of the optimization procedure. Finally, the last two rows represent the posterior collected by running sghmc with the corresponding priors.

From the analysis of these plots, we clearly see the benefit of placing a prior on the functions rather than on the parameters. First, the Wasserstein distance plots show satisfactory convergence, with the normalizing flow prior closely matching the gp prior. Second, as expected, the posteriors exhibit similar behavior according to the possible solutions realizable from the prior: classic priors tend to yield degenerate functions resulting in overconfidence in regions without data, while our g p-based priors (g pi-g, g pi-h, g pi-n f) retain information regarding lengthscale and amplitude.

.4 . 2 Comparison for Bayesian convolutional neural networks

We proceed with the experimental campaign on the CIFAR10 benchmark [START_REF] Krizhevsky | Learning multiple layers of features from tiny images[END_REF] with a number of popular cnn architectures: LeNet [START_REF] Lecun | Gradient-Based Learning Applied to Document Recognition[END_REF], VGG [START_REF] Simonyan | Very Deep Convolutional Networks for Large-Scale Image Recognition[END_REF] and PreResNet (K. [START_REF] He | Deep Residual Learning for Image Recognition[END_REF].

Regarding posterior inference, we implement sghmc which, after a burn-in phase of 10,000 iterations, collects 200 samples with 10,000 simulation steps in between. For a fair comparison, we do not use techniques such as data augmentation or adversarial examples in any of the experiments. Regarding the target g p prior, we place a hyper-prior LogNormal(log 8, 0.3) for variance; whereas the hyper-prior for length-scale is LogNormal(log 512, 0.3). We use a mini-batch size of N s = 128 and N M = 32 measurement points sampled from the empirical distribution of the training data regarding prior optimization.

Table 5.2 summarizes the results on the CIFAR10 test set with respect to accuracy and mean negative loglikelihood (mnll). These results demonstrate the effectiveness of the g p-induced priors, as evidenced by the improvements in predictive performance when using g pi-g and g pi-h priors compared to using f g and f h priors, respectively. Noticeably, the g pi-h prior offers the best performance with 76.51%, 87.03%, and 88.20% predictive accuracy on LeNet, VGG, and PreResNet, respectively. We observe that, for complex models (e.g., PreResNet and VGG), f g prior's results are improved by a large margin by tempering the posterior. This is in line with the results showed by [START_REF] Wenzel | How Good is the Bayes Posterior in Deep Neural Networks Really?[END_REF]. By contrast, in the case of LeNet, the predictive performance dramatically degraded when using temperature scaling. This is clear evidence supporting the hypothesis that a "tempered" posterior is only useful for over-parameterized models. Instead, by using g p induced priors, we consistently obtain the best results in most cases.

.4 . 3 Optimizing priors with data

Although we advocate for functional priors over b n ns, we acknowledge that a prior of this kind is essentially heuristic. A potentially more useful prior might be discovered by traditional means such as cross-validation (c v) or by running an empirical Bayes procedure (a.k.a. type-II maximum likelihood), which maximizes the marginal likelihood p(D; ψ) = p(D | w)p(w; ψ) dw w.r.t. the prior parameters. ) and using Bayesian optimization ( ). In the plots, each corresponds to a run of a single con guration, while highlights the Pareto front of the cross validation procedure. The gure also reports the N(0, 1) prior as , while is our proposal of using functional prior ( i-).

These methods impose however significant challenges: (i) for c v, the number of hyperparameters that needs to be optimized becomes exponentially large as the complexity of the neural network grows and the more fine-grained grids are considered the higher the number of the combinations of parameters to explore becomes; (ii) for empirical Bayes, we need to compute the exact marginal likelihood, which is always intractable for b n ns, thus requiring additional approximations like variational inference (v i) or the Laplace approximation.

We consider a simple case of a b n n with one hidden layer only; by adopting the simple parameterization of § 5.3.3, we shall have four parameters to optimize in total (i.e. the weight and bias variances of the hidden and the output layer).

In Figure 5.6, we demonstrate how our scheme behaves on a simple regression task (dataset is Powerplant) in comparison with a c v strategy featuring a grid size of 9 (for a total of 6561 configurations). To get results for the cross-validation procedure and to massively exploit all possible parallelization opportunities, we allocated a cloud platform with 16 workstations, for a total of 512 computing cores and 64 maximum parallel jobs. This required a bit more than one day, although the total CPU time approached 3 months. While grid-based routines are widely adopted by practitioners for cross-validation, we acknowledge that there are more efficient alternatives. To this extent, we also include Bayesian optimization [START_REF] Močkus | On Bayesian Methods for Seeking the Extremum[END_REF][START_REF] Snoek | Practical Bayesian Optimization of Machine Learning Algorithms[END_REF][START_REF] Nogueira | Bayesian Optimization: Open source constrained global optimization tool for Python[END_REF], a classical method for blackbox optimization which uses a Gaussian process as the surrogate function to be maximized (or minimized). As expected, c v indeed found marginally better configurations, but the amount of resources and time needed, even for such a small model, is orders of magnitude larger than what required by our scheme, making this procedure computationally infeasible for larger models, like cnns. To reiterate and to put things into perspective, the Wasserstein-based functional prior could be run on a 4-cores laptop in a reasonable time, while c v required an entire rack of servers.

.'

A common way to estimate hyper-parameters (i.e., prior parameters ψ) is to rely on the Bayesian Occam's razor (a.k.a. empirical Bayes), which dictates that the marginal likelihood p ψ (y) should be optimized with respect to ψ. There are countless examples where such simple procedure succeeds in practice (see, e.g., [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]Immer, Bauer, et al., 2021). We note however that marginal likelihood maximization for a large number of hyper-parameters can suffer from overfitting [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF][START_REF] Ober | The Promises and Pitfalls of Deep Kernel Learning[END_REF]. Nevertheless, we do not expect significant overfitting issues in our setting, as we focus on data that are characterized by a high level of structure (i.e. images). As we have seen, regular choices for the prior completely fail to capture the properties of such highly-structured outputs.

The marginal likelihood is obtained by marginalizing out the outputs f and the model parameters w, (5.10) where p(y | f) and p ψ (f) are the likelihood and the prior on functions, respectively.

p ψ (y) = p(y | f)p ψ (f)df ,
Unfortunately, in our context it is impossible to carry out this optimization due to the intractability of Equation (5.10).

Classic results in the statistics literature draw parallels between maximum likelihood estimation (m l e) and kl minimization [START_REF] Akaike | Information Theory and an Extension of the Maximum Likelihood Principle[END_REF], minimization of the divergence between the true data distribution and the marginal

p ψ (y).
This alternative view allows one to obtain a viable optimization strategy that relies on an empirical estimate of the data distribution π(y). This presents additional challenges however, as the empirical evaluation and optimization of kl divergences remains a well-known challenging problem [START_REF] Flam-Shepherd | Mapping Gaussian Process Priors to Bayesian Neural Networks[END_REF]. Although it is possible to evaluate k l (or any other f-divergence) empirically by leveraging results from convex analysis [START_REF] Nguyen | Estimating Divergence Functionals and the Likelihood Ratio by Convex Risk Minimization[END_REF], we have opted to substitute kl with an alternative metric that is more convenient from a computational perspective.

We thus employ the Wasserstein distance as a surrogate for kl divergence, so that we avoid the challenges of empirical k l estimation.

To summarize: (1) we would like to do prior selection by carrying out type-II m l e;

(2) the m l e objective is analytically intractable but the connection with k l minimization allows us to (3) swap the divergence with the Wasserstein distance, yielding a practical framework for choosing priors.

. .1 The distributionally-sliced Wasserstein distance

Given the two probability measures π and p ψ , both defined on R D for simplicity, the p-Wasserstein distance between π and p ψ is given by

W p p (π, p ψ ) = inf γ∈Γ (π,p ψ )
yy p γ(y, y )dydy , (5.12)

where Γ (π, p ψ ) is the set of all possible distributions γ(y, y ) such that the marginals are π(y) and p ψ (y ) [START_REF] Villani | Optimal Transport: Old and New[END_REF]. While usually analytically unavailable or computationally intractable, for D = 1 the distance has a simple closed form solution, that can be easily estimated using samples only [START_REF] Kolouri | Generalized Sliced Wasserstein Distances[END_REF].

The distributionally-sliced Wasserstein distance (dswd) takes advantage of this result by projecting the estimation of distances for high-dimensional distributions into simpler estimation of multiple distances in one dimension. The projection is done using the Radon transform R, an operator that maps a generic density function ϕ defined in R D to the set of its integrals over hyperplanes in R D , (5.13) where S D-1 is the unit sphere in R D and δ(•) is the Dirac delta [START_REF] Helgason | Integral Geometry and Radon Transforms[END_REF].

Rϕ(t, θ) := ϕ(r)δ(t -r θ)dr , ∀t ∈ R , ∀θ ∈ S D-1 ,
Using the Radon transform, for a given direction (or slice) θ we can project the two densities π and p ψ into one dimension and we can solve the optimal transport problem in this projected space. Furthermore, to avoid unnecessary computations, instead of considering all possible directions in S D-1 , distributionally-sliced Wasserstein distance (d s w d) proposes to find the optimal probability measure of slices σ(θ) on the unit sphere S D-1 , (5.14) where, for C > 0, M C is the set of probability measures σ such that E θ,θ ∼σ θ θ C (a constraint that aims to avoid directions to lie in only one small area). The 5 . 5 ' 99 direct computation of DSW p in Equation (5.14) is still challenging but admits an equivalent dual form,

DSW p (π, p ψ ) := sup σ∈M C E σ(θ) W p p Rπ(t, θ), Rp ψ (t, θ) 1/p ,
sup h∈H E σ(θ) W p p Rπ(t, h(θ)), Rp ψ (t, h(θ)) 1/p -λ C E θ,θ ∼ σ h(θ) h(θ ) + λ C C , (5.15)
where σ is a uniform distribution in S D-1 , H is the set of functions h : S D-1 → S D-1 and λ C is a regularization hyper-parameter. The formulation in Equation (5.15) is obtained by employing the Lagrangian duality theorem and by reparameterizing σ(θ) as push-forward transformation of a uniform measure in S D-1 via h. Now, by parameterizing h using a deep neural network with parameters φ, defined as h φ , Equation (5.15) becomes an optimization problem with respect to the network parameters. The final step is to approximate the analytically intractable expectations with Monte Carlo integration, max

φ 1 K K i=1 W p p Rπ(t, h φ (θ i )), Rp ψ (t, h φ (θ i )) 1/p - λ C K 2 K i,j=1 |h φ (θ i ) h φ (θ j )| + λ C C , with θ i ∼ σ(θ)
. Finally, we can use stochastic gradient methods to update φ and then use the resulting optima for the estimation of the original distance.

. . 2 Matching the marginal distribution to the data distribution via

Wasserstein distance minimization

We aim at learning the prior parameters by optimizing the marginal p ψ (y) obtained after integrating out the weights from the joint p ψ (y, w). The connection with empirical Bayes and k l minimization suggests that we can find the optimal ψ by minimizing the k l between the true data distribution π(y) and the marginal p ψ (y) . However, matching these two distributions is non-trivial due to their high dimensionality and the unavailability of their densities. To overcome this problem, we propose a sample-based approach using the distributional sliced 2-Wasserstein distance (Equation (5.15)) as objective: π(y) . (5.16) This objective function is flexible and does not require the closed-form of either p ψ (y) or π(y). The only requirement is that we can draw samples from these two distributions. Note that we can sample from p ψ (y), by first computing f after sampling from p ψ (w) and then perturbing the generated f by sampling from the likelihood p(y | f). In the next section, we will take a challenging application (unsupervised learning with an auto-encoding architecture) where this formulation of the prior learning problem will be implemented.

ψ = arg min ψ DSW 2 p ψ (y),

.6

The problem of learning useful representations of data that facilitate the solution of downstream tasks such as clustering, generative modeling and classification, is at the crux of the success of many machine learning applications (see, e.g., Bengio, Courville, et al., 2013, and references therein). From a plethora of potential solutions to this problem, unsupervised approaches based on autoencoders [START_REF] Cottrell | Image Compression by Back Propagation: A Demonstration of Extensional Programming[END_REF] are particularly appealing as, by definition, they do not require label information and have proved effective in tasks such as dimensionality reduction and information retrieval (G. E. Hinton and R. R. Salakhutdinov,

2006).

Autoencoders are neural network models composed of two parts, usually referred to as the encoder and the decoder. The encoder maps each input y i to a set of lowerdimensional latent variables z i . The decoder maps the latent variables z i back to the observations y i . The bottleneck introduced by the low-dimensional latent space is what characterizes the compression and representation learning capabilities of autoencoders. It is not surprising that these models have connections with principal component analysis, factor analysis and density networks (D. J. MacKay and Gibbs, 1999), and latent variable models [START_REF] Lawrence | Probabilistic Non-linear Principal Component Analysis with Gaussian Process Latent Variable Models[END_REF].

In applications where quantification of uncertainty is a primary requirement or where data is scarce, it is important to carry out a Bayesian treatment of these models by specifying a prior distribution over their parameters, i.e., the weights of the encoder/decoder. However, estimating the posterior distribution over the parameters of these models, which we refer to as Bayesian auto-encoders (baes), is generally intractable and requires approximations. Furthermore, the need to specify priors for a large number of parameters, coupled with the fact that autoencoders are not generative models, has motivated the development of vaes as an alternative that can overcome these limitations [START_REF] Kingma | Auto-Encoding Variational Bayes[END_REF]. Indeed, va es 5 .6 101 have found tremendous success and have become one of the preferred methods in modern machine-learning applications (see, e.g., Kingma and Welling, 2019, and references therein).

To recap, three potential limitations of baes hinder their widespread applicability in order to achieve a similar or superior adoption to their variational counterpart:

(i) lack of generative modeling capabilities; (ii) intractability of inference and (iii) difficulty of setting sensible priors over their parameters. In this work we revisit b a es and deal with these limitations in a principled way. In particular, we address the first limitation in (i) by employing density estimation in the latent space. Furthermore, we deal with the second limitation in (ii) by exploiting recent advances in m c m c and, in particular, s g h m c (T. [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF]. Finally, we believe that the third limitation (iii), which we refer to as the difficulty of carrying out model selection, requires a more detailed treatment because choosing sensible priors for Bayesian neural networks is an extremely difficult problem, and this is the main focus of this work.

.6.1 Formalization of Bayesian Autoencoders

An auto-encoder (ae) is a neural network parameterized by a set of parameters w, which transforms an unlabelled dataset,

y def = {y n } N n=1 , into a set of reconstructions ŷ def = {f n } N n=1 , with y n , f n ∈ R D .
An a e is composed of two components: (1) an encoder f enc which maps an input sample y n to a latent code z n ∈ R K , K D; and

(2) a decoder f dec which maps the latent code to a reconstructed datapoint f n . In

short, f = f(x; w) = (f dec • f enc )(x)
, where we denote w := {w enc , w dec } the union of parameters of the encoder and decoder. The Bayesian treatment of a es dictates that a prior distribution p(w) is placed over all parameters of f enc and f dec , and that this prior knowledge is transformed into a posterior distribution by means of Note that each conditional likelihood term is determined by the model architecture, the choice of w, and the input x n , but in order to keep the notation uncluttered, we write them simply as p(y n | w). In the Bayesian scheme, the prior and likelihood are both modeling choices. Before giving an in-depth treatment on priors for b a es in the next section, we briefly discuss the likelihood, which can be chosen according to the type of data. In our experiments, we mainly investigate image datasets, where pixel values are normalized in the [0, 1] range. Therefore, we rely on the continuous Bernoulli distribution (Loaiza-Ganem and Cunningham, 2019):

Ouput with Output with

p(y n | w) = D i=1 K(λ i )λ y n,i i (1 -λ i ) 1-y n,i := p(y n | f n ),
(5.18) where K(λ i ) is a properly defined normalization constant (Loaiza-Ganem and Cunningham, 2019) and

λ i = f i (x n ; w) = f n,i ∈ [0, 1]
is the i-th output from the bae given the input x n . We note that, as f n depends deterministically on w, we will use the above expression to refer to both p(y n | w) and p(y n | f n ), where the latter term will be of crucial importance when we define the functional prior induced over the reconstruction.

.6 103

.6. 2 The pathology of standard priors for s and how to x it

The choice of the prior is important for the Bayesian treatment of any model as it characterizes the hypothesis space (D. J. C. [START_REF] Mackay | A Practical Bayesian Framework for Backpropagation Networks[END_REF][START_REF] Murray | A Note on the Evidence and Bayesian Occam's Razor[END_REF]. Specifically for baes, one should note that placing a prior on the parameters of the encoder and decoder has an implicit effect on the prior over the network output (i.e. the reconstruction). In addition, the highly nonlinear nature of these models implies that interpreting the effect of the architecture is theoretically intractable and practically challenging. Several works argue that a vague prior such as N(0, 1) is good enough for some tasks and models, like classification with cnns (A. G. Wilson and Izmailov, 2020b).

However, for baes this is not enough, as illustrated in Figure 5.8. The realizations obtained by sampling weights/biases from a N(0, 1) prior indicate that this choice provides poor inductive bias. Meanwhile, by encoding better beliefs via an optimized prior, the samples can capture main characteristics intrinsic to the data, even when the model is fed with out-of-distribution inputs.

To overcome this problem, we propose a sample-based approach using the distributional sliced 2-Wasserstein distance in Equation (5.15) as objective:

ψ = arg min ψ DSW 2 p ψ (y), π(y) .
(5.19)

We will refer to the original paper (B. [START_REF] Tran | Model Selection for Bayesian Autoencoders[END_REF] for the experimental evaluations, as mainly perfomed by the first author. As a summary, the experiments include a demonstration of the effect of our model selection strategy, by considering scenarios in the small-data regime where the prior might not be necessarily tuned on the training set. In this way we are able to impose inductive bias beyond what is available in the training data. The bae with optimized priors performs better than the competing methods (va es and the b a e with standard prior) in the inference task for all training sizes, with slightly diminishing effect for larger sets, as expected.

Also, this pattern is true for different latent dimensions, where regardless of the dimensionality of the latent space, b a es with optimized priors deliver higher performances.

.7

Being able to perform Bayesian inference of neural networks represents a much sought-after objective to equip extremely flexible models with the capability of expressing uncertainty in a sound way [START_REF] Neal | Bayesian Learning for Neural Networks[END_REF]. Recent advances in m c m c sampling enabling for efficient parameter space exploration, combined with minibatching (T. [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF], have turned this long-standing challenge into a concrete possibility. However despite these advances, there have been only

few success stories involving the use of Bayesian inference techniques for neural networks [START_REF] Osawa | Practical Deep Learning with Bayesian Principles[END_REF][START_REF] Zhang | Cyclical Stochastic Gradient MCMC for Bayesian Deep Learning[END_REF][START_REF] Izmailov | What Are Bayesian Neural Network Posteriors Really Like?[END_REF]. We attribute this to the difficulties in specifying sensible priors for thousands/millions of parameters, while being able to understand and control the effect of these choices in the behavior of their output functions [START_REF] Duvenaud | Avoiding pathologies in very deep networks[END_REF].

In this chapter we analyzed the problem of selection good priors from two different perspectives: functional priors with target stochastic processes and empirical Bayes with the Wasserstein distance as a proxy to the marginal likelihood.

The difficulty in reasoning about functional priors for neural networks made us consider, for the first work, the possibility to enforce these by minimizing their distance to tractable functional priors, effectively optimizing the priors over model parameters so as to reflect these functional specifications. We chose to consider Gaussian processes, as they are a natural and popular choice to construct functional priors, whereby the characteristics of prior functions are determined by the form and parameters of Gaussian process kernel/covariance functions. While previous works attempted this by using the k l divergence between the functional priors [START_REF] Flam-Shepherd | Mapping Gaussian Process Priors to Bayesian Neural Networks[END_REF][START_REF] Flam-Shepherd | Characterizing and Warping the Function space of Bayesian Neural Networks[END_REF], the objective proves difficult to work with due to the need to estimate an entropy term based on samples, which is notoriously difficult. We proposed a novel objective based on the Wasserstein distance, and we showed that this objective offers a tractable and stable way to optimize the priors over model parameters. The attractive property of this objective is that it does not require a closed form for the target functional prior, as long as it is possible to obtain samples from it.

For the second work, we have reconsidered as a challenging benchmark the Bayesian treatment of autoencoders (a e) in light of recent advances in Bayesian neural networks. We have found that the main limitation of b a es lies in the difficulty of Variational inference techniques based on inducing variables provide an elegant framework for scalable posterior estimation in Gaussian process (g p) models. Besides enabling scalability, one of their main advantages over sparse approximations using direct marginal likelihood maximization is that they provide a robust alternative for point estimation of the inducing inputs, i.e. the location of the inducing variables. In this work we challenge the common wisdom that optimizing the inducing inputs in the variational framework yields optimal performance. We show that, by revisiting old model approximations such as the fully-independent training conditionals endowed with powerful sampling-based inference methods, treating both inducing locations and g p hyper-parameters in a Bayesian way can improve performance significantly. Based on stochastic gradient Hamiltonian Monte Carlo, we develop a fully Bayesian approach to scalable gp and deep gp models, and demonstrate its state-of-the-art performance through an extensive experimental campaign across several regression and classification problems.

.1

Bayesian kernel machines based on gps combine the modeling flexibility of kernel methods with the ability to carry out sound quantification of uncertainty [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]. Modeling and inference in g p models have evolved considerably over the last few years with key contributions in the direction of scalability to virtually any number of datapoints and generality within automatic differentiation frameworks (A. G. [START_REF] Matthews | GPflow: A Gaussian Process Library using TensorFlow[END_REF][START_REF] Krauth | AutoGP: Exploring the Capabilities and Limitations of Gaussian Process Models[END_REF]. This has been possible thanks to the combination of stochastic variational inference techniques with representations based on inducing variables [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF][START_REF] Lazaro-Gredilla | Inter-domain Gaussian Processes for Sparse Inference using Inducing Features[END_REF][START_REF] Hensman | Gaussian Processes for Big Data[END_REF], random features [START_REF] Rahimi | Random Features for Large-Scale Kernel Machines[END_REF][START_REF] Cutajar | Random Feature Expansions for Deep Gaussian Processes[END_REF]Gal and Ghahramani, 2016b), and structured approximations (A. [START_REF] Wilson | Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP)[END_REF]A. G. Wilson, Hu, R. R. Salakhutdinov, et al., 2016). These advancements have now made gps attractive to a variety of applications and likelihoods (A. G. [START_REF] Matthews | GPflow: A Gaussian Process Library using TensorFlow[END_REF][START_REF] Wilk | Convolutional Gaussian Processes[END_REF][START_REF] Bonilla | Generic Inference in Latent Gaussian Process Models[END_REF].

In this work, we focus on the variationally sparse gp framework originally formulated by [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF] and later developed by [START_REF] Hensman | Gaussian Processes for Big Data[END_REF] and Hensman, A. Matthews, et al. (2015) to scale up to large datasets via stochastic optimization. In these formulations, the gp prior is augmented with inducing variables (drawn from the same prior) and their posterior is estimated via variational inference. In contrast, the location of the inducing variables, which we refer to as the inducing inputs, are simply optimized along with covariance hyper-parameters.

In line with earlier evidence that Bayesian treatments of g ps are beneficial [START_REF] Neal | Monte Carlo Implementation of Gaussian Process Models for Bayesian Regression and Classification[END_REF][START_REF] Barber | Gaussian Processes for Bayesian Classification via Hybrid Monte Carlo[END_REF][START_REF] Murray | Slice Sampling Covariance Hyperparameters of Latent Gaussian Models[END_REF][START_REF] Filippone | Pseudo-Marginal Bayesian Inference for Gaussian Processes[END_REF], posterior inference of the inducing variables jointly with covariance hyper-parameters has been shown to improve performance [START_REF] Hensman | MCMC for Variationally Sparse Gaussian Processes[END_REF].

Despite these significant insights with regards to the benefits of full Bayesian inference over latent variables in g p models, the common practice is to optimize the inducing inputs, even in very recent gp developments [START_REF] Havasi | Inference in Deep Gaussian Processes using Stochastic Gradient Hamiltonian Monte Carlo[END_REF] et al., 2020;[START_REF] Giraldo | A Fully Natural Gradient Scheme for Improving Inference of the Heterogeneous Multi-Output Gaussian Process Model[END_REF]. In fact, the original work of [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF] advocates for a treatment of the inducing inputs as variational parameters to avoid overfitting. Furthermore, later work concludes that point estimation of the inducing inputs through optimization of the variational objective is an 'optimal' treatment (Hensman, A. G. Matthews, et al., 2015, §3). As we will see in § 6.2.1, the justification for inducing-input optimization in [START_REF] Hensman | MCMC for Variationally Sparse Gaussian Processes[END_REF] relies on being able to optimize both the prior and the posterior, and therefore, contradicts the fundamental principles of Bayesian inference. We summarize previous works on inference methods for gps in Table 6.1, which we will use for comparison in our experiments.

Thus, we revisit the role of the inducing inputs in g p models and their treatment as variational parameters or even hyper-parameters. Given their potential high dimensionality and that the typical number of inducing variables goes beyond hundreds/thousands [START_REF] Shi | Sparse Orthogonal Variational Inference for Gaussian Processes[END_REF], we argue that they should be treated simply as model variables and, therefore, having priors and carrying out efficient posterior inference over them is an important-although challenging-problem. An illustration of the richer modeling capabilities offered by treating inducing inputs in a Bayesian fashion is given in Figure 6.1.

( )

.

Firstly, we challenge the common wisdom that optimizing the inducing inputs in the variational framework yields optimal performance. We show that, by revisiting old model approximations such as the fully independent training conditionals (f i t c; see, e.g., [START_REF] Quiñonero-Candela | A Unifying View of Sparse Approximate Gaussian Process Regression[END_REF] endowed with powerful sampling-based inference methods, treating both inducing locations and g p hyper-parameters in a Bayesian way can improve performance significantly. We describe the conceptual justification and the mathematical details of our general formulation in § 6.2 and § 6.3. We then demonstrate that our approach yields state-of-the-art performance across a wide range of competitive benchmark methods, large-scale datasets and a variety of gp and deep gp models in § 6.4.

.2

We are interested in supervised learning problems with N input-label training

pairs {X, y} def = {(x i , y i )} N i=1
, where we consider a conditional likelihood p(y | f) and f is drawn from a zero-mean g p prior with covariance function k(x, x ; θ) with hyper-parameters θ. Thus, we have that p(f) = N(0, K xx|θ ), where K xx|θ is the N × N covariance matrix obtained by evaluating k(x i , x j ; θ) over all input pairs {x i , x j }. Inference in these types of models generally involves the costly O(N 3 ) operations to compute the inverse and log-determinant of the covariance matrix

K xx|θ .
. Sparse g ps are a family of approximate models that address the scalability issue by introducing a set of M inducing variables u = (u 1 , . . . , u M ) at corresponding inducing inputs Z = {z 1 , . . . , z M } such that u i = f(z i ) (see, e.g., [START_REF] Quiñonero-Candela | A Unifying View of Sparse Approximate Gaussian Process Regression[END_REF]. These inducing variables are assumed to be drawn from the same gp as the original process, yielding the joint prior p(f, u) = p(u)p(f|u). We consider a general formulation where we place priors p ψ (θ) over covariance hyper-parameters and p ξ (Z) over inducing inputs with hyper-parameters ψ, ξ,

p(θ, Z, u, f, y | X) = p ψ (θ)p ξ (Z)p(u | Z, θ)p(f | u, X, Z, θ)p(y | f), (6.1) where p(u | Z, θ) = N(0, K zz|θ ), p(f | u, X, Z, θ) = N(K xz|θ K -1 zz|θ u, K xx|θ -K xz|θ K -1 zz|θ K xz|θ ).
The matrices K zz|θ , K xz|θ denote the covariance matrices computed between points and make no assumptions about the other distributions. In our formulation, approaches that do not consider priors over covariance hyper-parameters or inducing inputs correspond to improper uniform priors in Equation (6.1).

2 .1 On scalable inference frameworks for GP models

Let Ψ def = {u, Z, θ} be the variables whose posterior we wish to infer. Our main object of interest is the log joint marginal obtained by integrating out the latent variables f in Equation (6.1), i.e., log p(y,

Ψ | X) = log f p(y | f)p(f | Ψ, X)df + log p(Ψ).
In particular, we are interested in discussing approximations to this that decompose over observations, allowing the use of stochastic optimization techniques to scale up to large datasets. In the literature of sparse g ps (see, e.g., [START_REF] Bauer | Understanding Probabilistic Sparse Gaussian Process Approximations[END_REF][START_REF] Bui | A Unifying Framework for Gaussian Process Pseudo-Point Approximations using Power Expectation Propagation[END_REF], two of the most influential methods for carrying out inference on such models are based on the variational free energy (vfe) framework [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF] and the fully independent training conditional (fitc) framework (E. [START_REF] Snelson | Sparse Gaussian Processes using Pseudo-Inputs[END_REF].

. [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF] is the definition of the approximate posterior q(f, u) def = q(u)p(f | Ψ, X), where q(u) is the variational posterior, which yields the evidence lower bound (elbo)

The key innovation in

p(y | X, Z, θ) -KL [q(u) p(u | Z, θ)] + E q(f,u) log p(y | f) def = L elbo . (6.2)
We note that this approach does not incorporate priors over inducing inputs or hyper-parameters. Inference involves constraining q(u) to a parametric form and finding its parameters to optimize the elbo. [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF] correctly argues that in the regression setting the variational approach to inducing variable approximations should be more robust to overfitting than a direct marginal likelihood maximization approach of traditional approximate models such as those described in Quiñonero-Candela and Rasmussen (2005). Indeed, if inducing inputs Z are optimized then the resulting e l b o provides an additional regularization term (see Titsias, 2009b, §3 for details). However, as we shall see later, the benefits of being Bayesian about the inducing inputs and estimating their posterior distribution can be superior to those obtained by this regularization.

Restricting the form of q(u) is suboptimal, and Hensman, A. G. Matthews, et al. (2015) proposes to sample from the optimal posterior approximation instead. By (6.3) This is the same expression derived in [START_REF] Hensman | MCMC for Variationally Sparse Gaussian Processes[END_REF], although following a different derivation showing that p vfe indeed yields the optimal distribution under the v f e framework of Equation (6.2). However, Hensman, A. G. Matthews, et al., 2015 argues that a Bayesian treatment of inducing inputs is unnecessary and concludes that the optimal prior is p(Z) = q(Z) = δ(Z -Z), where δ(•) is Dirac's delta function and Z is the set of inducing inputs that maximizes the elbo (Hensman, A. G. Matthews, et al., 2015, §3). We find such a justification flawed as it contradicts the fundamental principles of Bayesian inference. Indeed, the derivation in [START_REF] Hensman | MCMC for Variationally Sparse Gaussian Processes[END_REF] relies on minimizing both sides of the KL term in Equation (6.2), allowing for a 'free-form' optimization of the prior, which ultimately negates the necessity of all prior choices and defeats the purpose of a Bayesian treatment.

.

As an alternative, we can approximate the log joint of Equation (6.1) by imposing independence in the conditional distribution (see [START_REF] Quiñonero-Candela | A Unifying View of Sparse Approximate Gaussian Process Regression[END_REF], for details),

p(f | Ψ, X) = N K xz|θ K -1 zz|θ u, diag K xx|θ -K xz|θ K -1 zz|θ K xz|θ (6.4) which implies, log p(y, Ψ | X) ≈ N n=1 log E p(f n | Ψ,X) [p(y n | f n )] + log p(Ψ) def =log p f i t c (y,Ψ | X)
.

(6.5)

This same formulation of the fitc objective can be also been obtain by modifying the likelihood or the prior rather then the conditional distribution (E. [START_REF] Snelson | Sparse Gaussian Processes using Pseudo-Inputs[END_REF][START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF][START_REF] Bauer | Understanding Probabilistic Sparse Gaussian Process Approximations[END_REF].

We now see that, when considering i.i.d. conditional likelihoods, both approximations, log p V F E and log p F I T C , yield objectives that decompose on the observations, enabling scalable inference methods. In particular, we aim to sample from the posterior over all the latent variables using scalable approaches such as Stochastic Gradient Hamiltonian Monte Carlo (sghmc) (T. [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF]. The main ques-

2 113

tion is what approach should be preferred and how they relate to their optimization counterparts.

2 . 2 Sampling with VFE or FITC?

We will show in § 6.4 that our proposal that samples from the posterior according to Equation (6.5) consistently outperforms that in Equation ( 6.3). To understand why the f i t c objective makes sense we need to go back to the original work of [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF] and Titsias (2009a) and the seminal work of Quiñonero-Candela and Rasmussen (2005). For this, we will consider the regression case but it can be easily generalized to classification. Indeed, Titsias, 2009a shows that, in the standard regression case with homoskedastic observation noise, vfe yields exactly the same predictive posterior as the projected process (pp) approximation [START_REF] Seeger | Fast Forward Selection to Speed Up Sparse Gaussian Process Regression[END_REF], which is referred to as the deterministic training conditional (dtc) approximation.

The optimal variational posterior distribution is given by:

q * (u | θ) = N(u; m, S), (6.6)
m = σ -2 K zz K zz + σ -2 K zx K xz -1 K zx y S = K zz K zz + σ -2 K zx K xz -1 K zz
where σ 2 is the observation-noise variance. It is easy to show that, given a Gaussian posterior over the inducing variables with mean and covariance m and S, the posterior predictive distribution at test point x is a Gaussian with mean and variance

µ y (x ) = κ(x , Z)K -1 zz m (6.7) σ 2 y (x ) = κ(x , x ) -κ(x , Z)K -1 zz κ(Z, x ) + κ(x , Z)K -1 zz SK -1 zz κ(Z, x ).
Thus, replacing Equation (6.6) in Equation (6.7) we obtain:

µ y (x ) = σ -2 κ(x , Z)ΣK zx y (6.8) σ 2 y (x ) = κ(x , x ) -κ(x , Z)K -1 zz κ(Z, x ) + κ(x , Z) K zz + σ -2 K zx K xz -1 κ(Z, x ),
which indeed corresponds to the predictive distribution of the dtc/pp approximation. Despite this equivalence, as highlighted in Titsias, 2009b, the main difference is that the v f e framework provides a more robust approach to hyper-parameter estimation as the resulting elbo corresponds to a regularized marginal likelihood

( )
of the dtc approach and hence should be more robust to overfitting. Nevertheless, the dtc/pp, and consequently the vfe, predictive distribution has been shown to be less accurate than the f i t c approximation [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF][START_REF] Quiñonero-Candela | A Unifying View of Sparse Approximate Gaussian Process Regression[END_REF][START_REF] Snelson | Flexible and efficient Gaussian process models for machine learning[END_REF]. Effectively, as described in Quiñonero-Candela and Rasmussen, 2005, the v f e's solution (which is the same as d t c's) can be understood as considering a deterministic conditional prior p(f | u), i.e. a conditional prior with zero variance.

The fitc approximation considers the following approximate conditional prior:

p(f | u) ≈ N f; K xz K -1 zz u, diag K xx -K xz K -1 zz K zx = N n=1 p(f n | u) = N n=1 N(f n ; µ n , σ 2 n ), with µ n = κ(x n , Z)K -1
zz u (6.9)

σ 2 n = κ(x n , x n ) -κ(x n , Z)K -1 zz κ(Z, x n ).
(6.10)

As we shall see later, is this factorization assumption in the conditional prior that will yield a decomposable objective amenable to stochastic gradient techniques. For now, consider the posterior predictive distribution under the fitc approximation1

µ f i t c (x ) = κ(x , Z)Σ f i t c K zx Λ -1 y σ 2 f i t c (x ) = κ(x , x ) -κ(x , Z)K -1 zz κ(Z, x ) + κ(x , Z)Σ fitc κ(Z, x ),
where

Λ = diag(K xx -K xz K -1 zz K zx + σ 2 I)
and (6.11) We now see why fitc's predictive distribution above is more accurate than vfe's in Equation (6.8), as we can obtain f i t c's by replacing σ 2 I in v f e's solution with diag(K xx -K xz K -1 zz K zx ) + σ 2 I. Effectively, as described in Quiñonero-Candela and Rasmussen, 2005, v f e's solution (which is the same as d t c's) can be understood as considering a deterministic conditional prior p(f | u), i.e. with zero variance.

Σ f i t c = (K zz + K zx Λ -1 K xz ) -1 .

2 . 3 Stochastic Updates Using the FITC Approximation

Now we can understand why the log of the expectation can provide more accurate results than the expectation of the log. Basically in the former we are using the fitc approximation while in the later we are using the v f e/d t c/p p approximation. It is easy to show that when using the f i t c approximation, one can obtain a decomposable objective function that can be implemented at large scale using stochastic gradient techniques. Here we focus only on the expectation of the conditional likelihood (which is the crucial term) and in the regression setting for simplicity but the extension to the classification case (e.g. using quadrature) is straightforward.

log p(y | u, θ, Z) = log f p(f | u, θ, Z)p(y | f) df = log f 1 ,...,f N N n=1 p(y n | f n )p(f n | u, θ, Z) df = log N n=1 f n N(y n | f n , σ 2 )N(f n | µ n , σ 2 n ) df n = log N n=1 p(y n | u, θ, Z) = N n=1 log N(y n | µ n , σ 2 n + σ 2 ), (6.12)
where µ n , σ 2 n are given by Equation (6.9) and Equation (6.10). Similar results can be derived for binary classification with Bernoulli likelihood and response function λ(f):

log p(y | uθ) = log E p(f | u,θ) [p(y | f)] = log N n=1 f n N(f n ; µ n , σ 2 n )Bern(y n ; λ(f n )) df n . (6.13)
When the response function is the cdf of a standard Normal distribution, i.e.,

λ(f n ) = Φ(f n ) def = f n -∞ N(f n ; 0, 1)
df n , which is also known as the probit regression model, the expectation above can be computed analytically to obtain:

log p(y, u | θ) = N n=1 log Bern(y n ; Φ( µ n / 1 + σ 2 n )).
(6.14) For other response functions the expectation in Equation (6.13) can be estimated using quadrature.

( )

2 .4 An heteroskedastic version of the Gaussian likelihood

As Titsias (2009a) discussed in Appendix C, the f i t c approximation corresponds to a g p regression with heteroskedastic noise variance

p(y|f) = N(y|f, σ 2 I + diag[K xx -K xz K zz K zx ]).
(6.15)

If we apply this augmented likelihood to the variational expectations term, we get

E q(f) log p(y|f, σ 2 , θ) = - 1 2 n j=1 log 2π(σ 2 + σ 2 j ) + (y j -µ j ) 2 + σ 2 j σ 2 + σ 2 j . (6.16)
Since Titsias, 2009a considers this vfe formulation, we also compare with it.

2 . 5 Concluding Remarks

We conclude this section by placing our analysis in the context of 'exact' Bayesian inference techniques that consider priors over the inducing inputs. As mentioned above, the main reason for the superior performance of v f e, despite providing a less accurate predictive posterior than f i t c's, was that inducing input estimation was more robust due to the use of the variational objective, which provided an extra regularization term. However, by placing priors over the inducing inputs as well as over covariance hyper-parameters, the problem of regularization over these parameters becomes irrelevant. It is important to highlight that a variational formulation equivalent to f i t c has also been proposed (see Titsias, 2009a, App. C). In our experiments in § 6.4 we show that, nonetheless, our fully Bayesian formulation still yields superior performance to such an approach, confirming the benefits of carrying out full posterior estimation over the inducing inputs.

A final remark is that, in the spirit of Bayesian modeling, any uncertainty in the covariance should be accounted for. Thinking of gp hyper-parameters and inducing inputs as parameters of the covariance function, a distribution over these induces a distribution over the covariance function, which enriches the modeling capabilities of these models (see, e.g., [START_REF] Jang | Scalable Lévy Process Priors for Spectral Kernel Learning[END_REF] or Figure 6.1).

6. 3 117 6 .3

In this section we describe practical considerations in our Bayesian Sparse Gaussian Process (b s g p) framework, including inference techniques, prior choices and extensions to deep Gaussian processes. Recalling that Ψ = {θ, u, Z} represents the set of variables to infer and, using Equation (6.5), their posterior can be obtained as

log p(Ψ | y, X) = log E p(f | Ψ,X) p(y | f) + log p(u | θ, Z)+ log p ξ (Z) + log p ψ (θ) -log C.
(6.17)

We use Markov chain Monte Carlo (m c m c) techniques, in particular s g h m c (T. [START_REF] Chen | Stochastic Gradient Hamiltonian Monte Carlo[END_REF][START_REF] Havasi | Inference in Deep Gaussian Processes using Stochastic Gradient Hamiltonian Monte Carlo[END_REF], to obtain samples from the intractable p(Ψ | y, X). Unlike Hamiltonian Monte Carlo (hmc), which requires computing the exact gradient ∇ log p(Ψ | y, X) and the exact unnormalized posterior to evaluate the acceptance [START_REF] Neal | MCMC Using Hamiltonian Dynamics[END_REF], sghmc obtains samples from the posterior with stochastic gradients and without evaluating the Metropolis ratio (see supplement for details).

With a factorized likelihood p(y|f) and an energy function U(Ψ) = -log p(Ψ | y, X) + log C, we sample Equation (6.17) over minibatches of data.

3 .1 Prior choices

Next, we discuss prior choices for the inducing inputs and covariance hyperparameters. The inducing inputs Z support the sparse Gaussian process interpolation, which motivates matching the inducing prior to the data distribution p(X).

We begin by proposing a simple Normal (N) prior 6.18) which matches the mean and variance of the normalized data distribution, and favors inducing inputs toward the baricenter of the data inputs.

p N (Z) = M j=1 N(z j |0, I), (
We also explore two priors based on point processes, which consider distributions over point sets [START_REF] Gonzalez | Spatio-temporal Point Process Statistics: A Review[END_REF]. Point processes can induce repulsive ef- fects penalizing configurations where inducing points are clumped together. The determinantal point process (dpp), defined through

p D (Z) ∝ det K zz|θ , (6.19)
relates the probability of inducing inputs to the volume of space spanned by the covariance [START_REF] Lavancier | Determinantal Point Process Models and Statistical Inference[END_REF]. d p p is a repulsive point process, which gives higher probabilities to input diversity, controlled by the hyper-parameters ξ ≡ θ.

We then consider the Strauss process (see e.g. [START_REF] Daley | An introduction to the theory of point processes[END_REF][START_REF] Strauss | A Model for Clustering[END_REF], (6.20) where λ > 0 is the intensity, and 0 < γ 1 is the repulsion coefficient which decays the prior as a function of the number of input pairs that are within distance r.

p S (Z) ∝ λ M γ z,z ∈Z δ(|z-z |<r) ,
The Strauss prior (S) tends to maintain the minimum distance between inducing inputs, parameterized by ξ = (λ, γ, r). We finally consider an uninformative uniform prior (U), log p U (Z) = 0, which effectively provides no contribution to the evaluation of the posterior.

To gain insights on the choice of these priors, we set up a comparative analysis on the b a na na dataset (Figure 6.2). We observe that the posterior densities on the inducing inputs are multimodal and highly non-Gaussian, further confirming the necessity of free-form inference. Both Strauss and dpp-based priors encourage configurations where the inducing inputs are evenly spread. The Normal and Uniform priors, instead, focus exclusively on aligning the inducing inputs in a way that is sensible to accurately model the intricate classification boundary between the classes. This insight is confirmed by our the extensive experimental validation in § 6.4.

-. Choosing priors on the hyper-parameters has been discussed in previous works on Bayesian inference for gps (see e.g. [START_REF] Filippone | Pseudo-Marginal Bayesian Inference for Gaussian Processes[END_REF]. Throughout this chapter, we use the RBF covariance with automatic relevance determination (a r d), marginal variance σ and independent lengthscales λ i per feature [START_REF] Mackay | Bayesian methods for backpropagation networks[END_REF]. On these two hyper-parameters we place a lognormal prior with unit variance and means equal to 1 and 0.05 for λ and σ, respectively.

3 . 2 Extension to deep Gaussian processes

bsgp can be easily extended to dgp models [START_REF] Damianou | Deep Gaussian Processes[END_REF]. Deep Gaussian processes are hierarchical compositions of L sparse gps generally defined as

DGP(x) = (f L • . . . • f 1 ) (x)
Each layer f l is a vector-valued g p:

f i (•) = f (1) i (•), . . . , f (H i ) i (•) with f (j) i (•) i.i.d ∼ GP(0, κ i (•, • | θ i )).
where H i is the width of the i th g p which are taken to be independent. Each layer is associated with a set of inducing inputs Z (l) , inducing variables u (l) ( )

and hyper-parameters θ (l) [START_REF] Salimbeni | Doubly Stochastic Variational Inference for Deep Gaussian Processes[END_REF]. In our notation

Ψ = Ψ (l) L =1 = Z (l) , u (l) , θ (l) L l=1 .
In this section, we derive the mathematical basis for a Bayesian treatment of inducing inputs in a d g p setting [START_REF] Damianou | Deep Gaussian Processes[END_REF]. We assume a deep Gaussian process prior f 1) , where each f (l) is a g p. For notational brevity, we use f (0) as the input vector x. Then we can write down the joint distribution over visible and latent variables (omitting the dependency on X for clarity) as l) . (6.21)

(L) • f (L-1) • • • • f (
p y, f (l) , Ψ (l) L l=1 = p y | f (L) L l=1 p f (l) | Ψ (l) , f (l-1) p Ψ (
Our goal is to estimate the posterior after marginalizing out f (l) , l) ). (6.22)

log p ψ (l) L l=1 | y ∝ = log p y | f (L) p f (L) | ψ (L) , f (L-1) . . . p f (1) | ψ (1),f (0) df (L) df (L-1) . . . df (1) + L l=1 log p(Ψ (
While the distribution in Equation (6.22) is not immediately computable owing to the intractable expectation term, we have obtained the form of its (un-normalized) log posterior, from which we can sample using h m c methods. More calculations reveal that we can, nevertheless, obtain estimates of this expectation term with Monte Carlo sampling where,

log p y | f (L) p f (L) | ψ (L) , f (L-1) . . . p f (1) | ψ (1) df (L) df (L-1) . . . df (1) ≈ log p y | f (L) p f (L) | ψ (L) , f (L-1) df (L) (
f (1) ∼ p f (1) | ψ (1) , f (0) f (2) ∼ p f (2) | ψ (2) , f (1) • • • f (L-1) ∼ p f (L-1) | ψ (L-1) , f (L-2)
Because of the layer-wise factorization of the joint likelihood (Equation (6.21)), each step of the approximation is unbiased. While it is possible to approximate the lastlayer expectation with a Monte Carlo sample, the expectation is tractable when the likelihood is a Gaussian or a Bernoulli distribution with a probit regression model, or is computable with one-dimensional quadrature [START_REF] Hensman | MCMC for Variationally Sparse Gaussian Processes[END_REF].

.4

In this section, we provide empirical evidence that our bsgp outperforms previous inference/optimization approaches on shallow and deep gps. We use eight of the , where the posterior is constrained to be Gaussian and the remaining parameters are point-estimated, to our proposal , where we infer a free-form posterior for all Ψ = {u, θ, Z}. We refer the reader to Table 6.1 for details on the methods (colors are matched).

classic UCI benchmark datasets with standardized features and split into eight folds with 0.8/0.2 train/test ratio. We train the competing models for 10,000 iterations with a da m [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF], step size of 0.01 and a minibatch of 1,000 samples. The sampling methods are evaluated based on 256 samples collected after optimization. Following previous works (e.g. [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF][START_REF] Havasi | Inference in Deep Gaussian Processes using Stochastic Gradient Hamiltonian Monte Carlo[END_REF][START_REF] Yu | Implicit Posterior Variational Inference for Deep Gaussian Processes[END_REF], in order to evaluate and compare the full predictive posteriors we compute the mean negative loglikelihood (mnll) on the test set (rmses are reported in the supplement for reference).

.1 Prior analysis and ablation study

We start our empirical analysis with a comparative evaluation of the priors on inducing inputs described in § 6.3.1: d p p, Normal, Strauss and Uniform. We run our inference procedure on a shallow g p with 100 inducing points and we report the results in ?? (left). The results show that the Normal prior consistently outperforms the others. The uniform and Strauss priors behave similarly, while the d p p prior is consistently among the worst. We argue that the repulsive nature of the point process priors (dpp, particularly), although grounded on the intuition of covering the input space more evenly, constrains the smoothness of the functions up to the point that they become too simple to accurately model the data. With this, we select the Gaussian prior for the remaining experiments.

We now study the benefits of a Bayesian treatment of the inducing variables, inducing inputs, and hyper-parameters with an ablation study. Using the same setup as before, we start with the baseline of s v g p (Hensman, Fusi, et al We refer the reader to Table 6.1 for a description of the methods.

Hensman, A. Matthews, et al., 2015), where the posterior on u is approximated using a Gaussian and Z, θ are optimized. We then incrementally add parameters to the list of variables that are sampled rather than optimized: only u (equivalent to s g h m cd g p, [START_REF] Havasi | Inference in Deep Gaussian Processes using Stochastic Gradient Hamiltonian Monte Carlo[END_REF]), then {u, θ} and finally, our proposal, {u, θ, Z}. This experiment is repeated for different number of inducing points (10, 50 and 100).

?? (right) reports a summary of these results (full comparison in the supplement). This plot shows that each time we carry out free-form posterior inference on a bigger set of parameters rather than optimization, performance is enhanced, and our proposal outperforms previous approaches.

. 2 Choosing the objective: VFE vs FITC

In § 6.2 we discussed the role of the marginal and the variational free energy (v f e) objective when used for optimization and for sampling. In Figure 6.6 we support the discussion with empirical results. The baseline is svgp, for which the inference is approximate (Gaussian) and performed on the variational objective.

Titsias (2009a, App. C) also considers a vfe formulation of fitc which corresponds to a gp regression with heteroskedastic noise variance. The likelihood needs to be augmented to handle heteroskedasticity, but inference can be carried out exactly on the variational objective. For these two methods, {θ, Z} are optimized. We also test mcmc-svgp, the model proposed by [START_REF] Hensman | MCMC for Variationally Sparse Gaussian Processes[END_REF], implemented in GPflow (A. G. [START_REF] Matthews | GPflow: A Gaussian Process Library using TensorFlow[END_REF] with the same suggested experimental setup. This experiment indicates that having a free-form posterior on u, θ sampled from the variational objective does not dramatically improve on the exact Gaussian approximation of the f i t c model, with both of them delivering superior performance with respect to s v g p. In the same setup of Hensman, A. G. Matthews, et al. (2015) (u, θ sampled and Z optimized), we look at the effect of swapping the expectation of log with the log of expectation (which effectively means moving from the v f e objective to f i t c); on the contrary, here we observe a significant increase in performance when using the latter, further confirming the discussion of the objectives in § 6.2.2. We finally conclude this section with an experiment where we try both objectives on our proposed bsgp and also different depths of the d g p (Figure 6 1945), we test the null hypothesis of v f e objective being better than the proposed f i t c. Figure 6.7 shows that, for the majority of the cases, this can be rejected (p < 0.05).

. 3 Deep Gaussian processes on UCI benchmarks

We now report results on d g ps. We compare against two current state-of-the-art deep g p methods, s g h m cd g p [START_REF] Havasi | Inference in Deep Gaussian Processes using Stochastic Gradient Hamiltonian Monte Carlo[END_REF] and i p v id g p (H. [START_REF] Yu | Implicit Posterior Variational Inference for Deep Gaussian Processes[END_REF], and against the shallow svgp baseline (Hensman, A. Matthews, et al., 2015).

For a faithful comparison with ipvi-dgp we follow the recommended parameter configurations 2 . Using a standard setup, all models share M = 100 inducing points, the same RBF covariance with a r d and, for d g p, the same hidden dimensions (equal to the input dimension D). and 95% CI over the different folds over the UCI datasets, and also includes rank summaries. The proposed method clearly outperforms competing deep and shallow gps. The improvements are particularly evident on naval, a dataset known to be challenging to improve upon. Furthermore, the deeper models perform consistently better or on par with the shallow version, without incurring in any measurable overfitting even on small or medium sized datasets (see boston and yacht, for example).

. Similarly to the baseline algorithms, each training iteration of bsgp involves the computation of the inverse covariance with complexity O(M 3 ). In Figure 6.10 we compare the three main competitors with b s g p trained for a fixed training time budget of one hour for a shallow g p and a 2-layer d g p. The experiment is repeated four times on the same fold and the results are then averaged. Each run is performed on an isolated instance in a cloud computing platform with 8 CPU cores and 8 GB of reserved memory [START_REF] Pace | Flexible Scheduling of Distributed Analytic Applications[END_REF]. Inference on the test set is performed every 250 iterations. This shows that b s g p converges considerably faster in wall-clock time, even though a single gradient step requires slightly more time.

Computing the predictive distribution, on the other hand, is more challenging as it requires recomputing the covariance matrices K xz , K zz for each posterior sample {Z, θ}, for an overall complexity linear in the number of posterior samples. This operation can be easily parallelized and implemented on g p us but it could question the practicality of using a more involved inference method. In particular, it is relevant to study whether s v g p could deliver superior performance with a higher number of inducing points for less computational overhead. In Figure 6.11 we study this trade-off on the biggest dataset considered (Protein): while it is evident that predictions with bsgp take more time (assuming a serial computation of the covariance matrices), it is also clear that the number of inducing points svgp requires to (even marginally) improve upon b s g p is significantly larger (up to 20 times).

. Finally, we run one last comparison with methods which exploit structure in the inputs. These models allow one to scale the number of inducing variables while maintaining computational tractability. Kernel Interpolation for Scalable Structured Gaussian Processes (kiss-gp) (A. [START_REF] Wilson | Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP)[END_REF] proposes to place the inducing inputs on a fixed and equally-spaced grid and to exploit Toeplitz/Kronecker structures with an iterative conjugate gradient method to further enhance scalability. Despite these benefits, kiss-gp is known to fall short with high-dimensional data (D > 4). This shortcoming was later addressed with Deep Kernel Learning (dkl) (A. G. Wilson, Hu, R. Salakhutdinov, et al., 2016): using a deep neural network d k l projects the data in a lower dimensional manifold by learning an useful feature representation, which is then used as input to a kiss-gp. In Figure 6.12 we have the comparison of bsgp with these two methods. kiss-gp could only run on Powerplant, with a 4-dimensional grid of size 10 (for a total of 10,000 inducing points). Here, b s g p delivers better performance despite having less inducing points. For d k l we followed the suggestion of A. G. Wilson, Hu, R. Salakhutdinov, et al. (2016) to use a fully-connected neural network with a [d -1000 -1000 -500 -50 -2] architecture as feature extractor and a grid size of 100 (for again a total of 10,000 inducing points). Training is performed by alternating optimization of the neural network weights and the k i s sg p parameters. Thanks to the flexibility of the feature extractor, this configuration is very competitive with our shallow b s g p, but it yields lower performance compared to a 2-layer d g p.

.4 Large scale classi cation

The airline dataset is a classic benchmark for large scale classification. It collects delay information of all commercial flights in USA during 2008, counting more than 5 millions data points. The goal is to predict if a flight will be delayed based on 8 features, namely month, day of month, day of week, airtime, distance, arrival time, departure time and age of the plane. We pre-process the dataset following the guidelines provided in [START_REF] Hensman | MCMC for Variationally Sparse Gaussian Processes[END_REF]A. G. Wilson, Hu, R. Salakhutdinov, et al., 2016).

After a burn-in phase of 10,000 iterations, we draw 200 samples with 1000 simulation steps in between. We test on 100,000 randomly selected held-out points. We fit three models with M = 100 inducing points. Table 6.2 shows the predictive performance of three shallow GP models. The b s g p yields the best test error, m n l l, and test area under the curve (au c). We assess the convergence of the predictive posterior by evaluating the R-statistics [START_REF] Gelman | Bayesian Data Analysis[END_REF] over four independent sghmc chains. This diagnostic yielded a R = 1.02 ± 0.045, which indicates good convergence. We report further convergence analysis in the supplement.

As a further large scale example, we use the h i g g s dataset [START_REF] Baldi | Searching for Exotic Particles in High-Energy Physics with Deep Learning[END_REF], which has 11 millions data points with 28 features. This dataset was created by Monte Carlo simulations of particle dynamics in accelerators to detect the Higgs boson. We select 90% of the these points for training, while the

.5

We have developed a fully Bayesian treatment of sparse Gaussian process models that considers the inducing inputs, along with the inducing variables and covariance hyper-parameters, as random variables, places suitable priors and carries out approximate inference over them. Our approach, based on s g h m c, investigated two conventional priors (Gaussian and uniform) for the inducing inputs as well as two point process based priors (the Determinantal and the Strauss processes).

By challenging the standard belief of most previous work on sparse g p inference that assumes the inducing inputs can be estimated point-wisely, we have developed a state-of-the-art inference method and have demonstrated its outstanding performance on both accuracy and running time on regression and classification problems. We hope this work can have an impact similar (or better) to other works in machine learning that have adopted more elaborate Bayesian machinery ( ) (e.g. [START_REF] Wallach | Rethinking LDA: Why Priors Matter[END_REF] for long-standing inference problems in commonly used probabilistic models.

Finally, we believe it is worth investigating further more structured priors similar to those presented here (e.g. exploring different hyper-parameter settings), including a full joint treatment of inducing inputs and their number, i.e. p(Z, M).

F I N A L C O N S I D E R A T I O N S

In this thesis, we discussed several challenges of using approximate Bayesian inference for deep learning models. Broadly speaking, Bayesian deep learning is hindered by problems related to scalability and quality of the approximation inference. In summary, we analyzed several issues that arises from applying variational inference (vi) to deep neural networks, we studied the role of the prior for improving modeling flexibility and we exploited the connections between neural networks and Gaussian processs (g ps) to refine the model and inference approximations of sparse g ps with inducing points.

.1

This thesis makes contributions which follow four different themes:

•

Initialization of variational inference [Chapter 3]

In this chapter we discussed how the role of initialization is severely under reported in the literature of variational inference for Bayesian neural networks (b n ns). The reasons for the failing of variational inference for b n ns in practice might depends heavily on the initialization of the variational parameters. We therefore proposed a novel way to initialize variational parameters, called iterative Bayesian linear modeling (i b l m), which is based on an iterative layer-wise initialization carried out through on Bayesian linear models. Empirical evidence is shown through a series of experiments, including regression and classification with deep neural networks (d n ns)

and convolutional neural networks (cnns) which demonstrated the ability of our approach to consistently initialize the optimization in a way that makes convergence faster than alternatives.

Open questions: recent works have exploited loss-landscape analysis and the model-connectivity of deep networks to-among other things-understand how and why simple Bayesian approximations like deep ensembles work in practice. This methodology could be similarly used for studying the loss landscape of the variational objective, and possibly deriving new initializations.

•

Structured parameterizations of variational posteriors [Chapter 4]

Inspired by the literature on scalable kernel methods, this chapter proposed Walsh-Hadamard variational inference (w h v i). w h v i offers a novel parameterization of the variational posterior, which is particularly attractive for over-parameterized models, such as modern d n ns and c n ns. w h v i assumes a matrix-variate posterior distribution, which therefore captures covariances across weights. Crucially, unlike previous works on matrix-variate posteriors for v i, this is achieved with a parsimonious parameterization and fast computations, bypassing the over-regularization issues of v i for over-parameterized models. The large experimental campaign demonstrates that whvi is a strong competitor with other variational approaches for such models, while offering considerable speedups.

Open questions: over-parameterization is tightly connected with generalization capabilities of deep models; in fact model performance is affected by the network size with bigger and wider neural networks being more resilient to overfit. At the same time, how this analysis carries over to bnn trained with variational inference is still open to debate.

• Priors in Bayesian deep learning [Chapter 5]

In this chapter we analyzed the problem of selection good priors using two different perspectives: functional priors and empirical Bayes using the Wasserstein distance as a proxy to the marginal likelihood. For the first case, we proposed a novel objective based on the Wasserstein distance, and we showed that this objective offers a tractable and stable way to optimize the priors over model parameters. Starting from samples of the neural networks prior, this is done by minimizing their distances to tractable functional priors (e.g. a Gaussian process), effectively optimizing the priors over model parameters so as to reflect these functional specifications. Alternatively, we can optimize the prior parameters akin to model selection for Gaussian process via marginal likelihood maximization. For this, we derived a practical and efficient optimization framework, based on the minimization of the distributional sliced-Wasserstein distance between the distribution induced by the model and the data generating distribution.

Open questions:

the main open question is how we can choose the target stochastic process to map b n n priors to. In principle, this can be generalized, for example, to match deep Gaussian processs (d g ps) or we could use a different prior parameterization to also enforce other behaviors, like e.g. sparsity. On the other hand, we can also debate whether there could be better alternatives to empirical Bayes with marginal likelihood maximization for model selection.

• Scalable approximations for (deep) Gaussian processes [Chapter 6]

In this chapter we proposed a fully Bayesian treatment of sparse Gaussian process models by Stochastic Gradient Hamiltonian Monte Carlo (s g h m c) sampling of the Gaussian process posterior with respect to all hyperparameters, inducing inputs and also inducing locations. In fact, despite significant insights with regards to the benefits of full Bayesian inference over latent variables in g p models, the common practice is to optimize the inducing inputs. The Bayesian treatment considers a wider family of data supporting hypotheses, especially in terms of inducing locations whose marginal support we can represent more accurately. In essence, we showed that, by revisiting old model approximations such as the fully independent training conditional (f i t c) endowed with powerful sampling-based inference methods, treating both inducing locations and g p hyper-parameters in a Bayesian way can improve performance significantly.

Open questions: it is possible to investigate further more structured priors on the inducing locations (e.g. exploring different hyper-parameter settings), including a full joint treatment of inducing inputs and their number. Alternatively, this framework could be in principle extended to inter-domain Gaussian processes with inducing variables, which include-among othersthe (deep) convolutional Gaussian process.

.2 ? ?

We conclude the thesis with a digression and speculations on the future of Bayesian deep learning in general. Looking back four years at the beginning of this thesis, we can better appreciate the rapid and prolific evolution of Bayesian deep learning (Figure 7.1). By late 2017 when this thesis started (at least as a preliminary discussion with Prof. Filippone), not many works were showing practical inference of Bayesian neural networks to succeed beyond one hidden layer and few hidden units, let alone for more complex architectures like large scale cnns. In the span of a couple of years this changed, thanks to an increased interest in scalable probabilistic deep models as well as to the democratization of flexible libraries (like Tensorflow and PyTorch) and powerful computational devices (like GPUs and TPUs). This growth is even more amazing, considering that it originated concurrently for different Bayesian approximations, whether it is variational, sampling or ensembles.

We have reached a point now which makes us wonder if we can consider Bayesian deep learning methodologically solved. While this sounds very provocative, it is important to remember that the field as a whole evolved so quickly that what could not have been done four years ago is today possible. Whether it is variational inference, or Laplace approximation, or ensembles or scalable Markov-Chain Monte Carlo sampling, we have significant empirical evidence and knowledge supporting whichever method. Practitioners have today a tool-set of countless probabilistic methods to choose from, with different models being more suitable than others for different application scenarios. Yet, with only few exceptions, Bayesian deep learning is still rarely employed for end-to-end applications, despite clear benefits for uncertainty quantification and calibration. Arguably, much of the deep learning success can be attributed to simple prototyping and simple deployment and monitoring of the models. While new trends are emerging in this direction [START_REF] Novak | Neural Tangents: Fast and Easy Infinite Neural Networks in Python[END_REF]E. Daxberger et al., 2021), this is generally still out of reach for Bayesian deep learning. With this thesis, we have given the practitioners a series of new tools and suggestions to make Bayesian deep learning work in practice, with the hope that applied research could benefit from it.

From a different point of view, we can also question how much (Bayesian) deep learning is sustainable and truly scalable in the long term. In a recent survey on the future role of deep learning for discovery of new fundamental physics [START_REF] Karagiorgi | Machine Learning in the Search for New Fundamental Physics[END_REF], the Authors project that in order to guarantee statistical significance for new discoveries, the model and the overall data-processing system should be able to keep up with data streams of at least 40 Tbit/s. This is already embarrassingly problematic for deep learning (for comparison, it is equivalent to roughly 30

ImageNet every second), let alone for Bayesian inference, to the point that we might need to rethink the entire computational pipeline. Uncertainty is omnipresent in nature [START_REF] Heisenberg | Uber den anschaulichen Inhalt der quantentheoretischen Kinematik und Mechanik[END_REF][START_REF] Kennard | Zur Quantenmechanik einfacher Bewegungstypen[END_REF]) and yet our computational devices are for the most part deterministic. On a long term horizon, quantum devices could ease the burden of probabilistic inference [START_REF] Zhao | Quantum-assisted Gaussian process regression[END_REF][START_REF] Schuld | Machine learning with quantum computers[END_REF][START_REF] Huang | Power of data in quantum machine learning[END_REF]. Also, while until now Bayesian deep learning has been chasing deep learning in the race for bigger models and new architectures, with quantum computing Bayesian inference could lead deep learning-and machine learning in general-to novel and exciting research directions. This also includes the possibility to redefine our objectives. This thesis, among countless of other examples in the literature, praised the elegance of Bayes' rule for updating beliefs with evidence, but Bayesian statistics should not be considered as an immovable dogma [START_REF] Popper | The Logic of Scientific Discovery[END_REF]. For example, even without entering in the philosophical debate of "frequestist" versus "Bayesian", it is possible to generalize Bayes' theorem in such a way that inference is carried out on bounds of probability measures [START_REF] Dempster | New Methods for Reasoning Towards Posterior Distributions Based on Sample Data[END_REF][START_REF] Dempster | A Generalization of Bayesian Inference[END_REF] or to unify and derive new methods from a common root of learning algorithms (M. E. [START_REF] Khan | The Bayesian Learning Rule[END_REF]. Consider an invertible, continuous and differentiable function f : R D → R D . Given z 0 ∼ q(z 0 ), then z 1 = f( z 0 ) follows q(z 1 ) defined as q(z 1 ) = q(z 0 ) det ∂f ∂z 0 -1

. (a.1)

As a consequence, after K transformations the log-density of the final distribution is log q(z K ) = log q(z 0 ) -

K k=1 log det ∂f k-1 ∂z k-1 . (a.2)
We shall define f k (z k-1 ; λ k ) the k th transformation which takes input from the previous flow z k-1 and has parameters λ k . The final variational objective is L(θ, φ) = -E q φ (z) [log p θ (x | z)] + KL q φ (z) p(z) = = E q φ (z) [-log p θ (x | z) -log p(z) + log q φ (z)] = = E q 0 (z 0 ) [-log p θ (x | z K ) -log p(z K ) + log q K (z K )] = = E q 0 (z 0 ) -log p θ (x | z K ) -log p(z K ) + log q 0 (z 0 ) -

K k=1 log det ∂f k (z k-1 ; λ k ) ∂z k-1 =
= -E q 0 (z 0 ) log p θ (x | z) + KL [q 0 (z 0 ) p(z K )] -E q 0 (z 0 )

K k=1 log det ∂f k (z k-1 ; λ k ) ∂z k-1 . (a.3)
Setting the initial distribution q 0 to a fully factorized Gaussian N(z 0 | µ, σI) and assuming a Gaussian prior on the generated z K , the Kullback-Leibler (kl) term is analytically tractable. A possible family of transformation is the planar flow (D. [START_REF] Rezende | Variational Inference with Normalizing Flows[END_REF]. For the planar flow, f is defined as

a.2

Given two Borel's probability measures π(x) and ν(y) defined on the Polish space X and Y (i.e. any complete separable metric space such as a subset of R d ), the p-Wasserstein distance is defined as follows where D(x, y) is a proper distance metric between two points x and y in the space X × Y and Γ (π, ν) is the set of functionals of all possible joint densities whose marginals are indeed π and ν.

When the space of x and y coincides (i.e. x, y ∈ X ⊆ R d ), the most used formulation is the 1-Wasserstein distance with Euclidian norm as distance, W(π, ν) = inf γ∈Γ (π,ν) X×X

xy γ(x, y) dx dy , (a.7) This is also known in the literature as the Earth-Mover distance. Intuitively, here γ measures how much mass must be transported from x to y in order to transform

. 2 165 the distributions π into the distribution ν. Solving the Wasserstein distance means computing the minimum mass that needs to be moved. The question "How?" is answered by looking at the optimal transport plan (not the focus of these notes).

The remaining part of these notes will be dedicated to the proof of the dual formulation for Equation (a.7). It is well known in the literature of optimization that linear programming problem with convex constrains admits a dual formulation. Kantorovich introduced the dual formulation of the Wasserstein distance in 1942.

Theorem 1 On the same setup as before, the Wasserstein distance defined as

W(π, ν) = inf γ∈Γ (π,ν) X×X
xy γ(x, y) dx dy , (a.8) admits the following dual form

W(π, ν) = sup f L 1 X f(x)π(x)dx - X f(y)ν(y)dy (a.9)
where f is a 1-Lipschitz continuous function defined on X → R.

Step 1: Kantorovich duality

First of all we start with the Kantorovich duality, which defines a dual form for the generic 1-Wasserstein. (a.11)

We define ι Γ (γ) the following quantity ι Γ (γ) = sup f,g f(x)π(x) dx + g(y)ν(y) dy -[f(x) + g(y)] γ(x, y) dx dy (a.12) and we observe that Proofs that the hypothesis used for the minimax theorem hold for this case are presented in Theorem 1.9 of "Topics in Optimal Transport" [START_REF] Villani | Topics in Optimal Transportation[END_REF]. Focusing on the infimum part, we can write ) takes a negative value at some point (x 0 , y 0 ), then by choosing γ = λδ(x 0 , y 0 ) with λ → +∞ (i.e. a Dirac delta in (x 0 , y 0 )), we see that the infimum is infinite. On the other hand, is ζ(x, y) is nonnegative, then the infimum is obtained for γ = 0. Finally, this constrains can be added to the previous conditions making thus recovering the formulation in Equation (a.9).

Step 2: D-Transforms

The next challenge is to find f and g such that we can easily recover the constrain optimization above. We approach this problem by supposing to have chosen some f(x). This means that the objective is to find a good g(y) that for all x, y satisfy the condition f(x) + g(y) D(x, y) . (a.19)

The trivial solution is g(y) D(x, y)f(x). This must be true for all x, also in the worst case (when we take the infimum), If we constrain f to be D-concave, then f = f.

Step 2.1: Euclidean distance It's worth mentioning that this formulation is valid for any nonnegative measurable function D. For the Euclidian distance this simplify even further.

Theorem 3 When D(x, y) = xy and f is 1-Lipschitz, f is D-concave if and only if f = -f

We prove the necessity condition of such result. First of all, we observe that if f is 1-Lipschitz then f is 1-Lipschitz too. This is true because for any given

x fx (y) = x -y -f(x) (a.25) is 1-Lipschitz and therefore the infimum of f(y) = inf x xyf(x) is 1-Lipschitz.

Since f is 1-Lipschitz, for all x and y we have f(y) -f(x) yx =⇒ -f(x)

xy -f(y) (a.26) where the right inequality follows by choosing y = x in the infimum. We know that f ≡ f. This means thatf(x) must be equal to f(x) for the last equation to hold.

Step 3. Putting everything together

We started our discussion by proving the Kantovich duality, which states that Finally, given D(x, y) to be the Euclidean distance, we discussed the shape of f when we restrict f to be 1-Lipschitz, showing that f = -f. Putting everything together, we obtain the dual 1-Wasserstein distance in Equation (a.9), W(π, ν) = sup after 100 epochs, vi inference initialized with iblm provides on average an increase up to 14% of accuracy at test time. Full results are reported in the supplementary material. This document was typeset with L A T E X using the typographical look-and-feel classicthesis developed by André Miede and Ivo Pletikosić.

Final Version as of December 17, 2021 (1).

5 . 1

 1 The effect of initialization in deep variational neural networks 3. 3. 2 Scaling up variational inference to deep convolutional neural networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3. 3. 3 Comparison with variational inference beyond mean field Gaussian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 . 4 Final remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The problem of overparameterization in variational inference . . . . . 4. 1. 1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 . 2 Structured Approximations for Kernel Matrices . . . . . . . . . . . . . 4 . 3 From structured kernel approximations to Walsh-Hadamard Variational Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. 3. 1 Statistical properties of the structure induced by Walsh-Hadamard variational inference (whvi) . . . . . . . . . . . . . 4. 3. 2 Reparameterizations in whvi for stochastic optimization . .4 . 4 Alternative structures, tensor factorization and extensions . . . . . . . 4. 4. 1 Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 . 5 Empirical evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. 5. 1 Toy example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. 5. 2 Empirical comparison on the UCI benchmark . . . . . . . . . . 4. 5. 3 Bayesian convolutional neural networks for image classification 4. 5. 4 Comments on computational efficiency . . . . . . . . . . . . . . 4. 5. 5 Exploring the parameter efficiency of WHVI with Gaussian processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 . 6 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 . 7 Final remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 The choice of the prior matters . . . . . . . . . . . . . . . . . . . . . . . 5. 1. 1 Pathologies of deep prior functions . . . . . . . . . . . . . . . . 5. 1. 2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 . 2 The problem of choosing priors in the literature of Bayesian deep learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 . 3 Imposing Gaussian process priors in Bayesian neural networks . . . . xiii 5. 3. 1 A quick introduction to the Wasserstein distance . . . . . . . . 5. 3. 2 Using the Wasserstein distance to impose the GP behaviour in BNN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. 3. 3 Prior Parameterization for Neural Networks . . . . . . . . . . 5 . 4 Empirical evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. 4. 1 Visualization on a 1D regression synthetic dataset . . . . . . . 5. 4. 2 Comparison for Bayesian convolutional neural networks . . . 5. 4. 3 Optimizing priors with data . . . . . . . . . . . . . . . . . . . . 5 . 5 Another route for Bayesian Occam's razor . . . . . . . . . . . . . . . . 5. 5. 1 The distributionally-sliced Wasserstein distance . . . . . . . . 5. 5. 2 Matching the marginal distribution to the data distribution via Wasserstein distance minimization . . . . . . . . . . . . . . 5 . 6 Model selection for Bayesian Autoencoders . . . . . . . . . . . . . . . . 5. 6. 1 Formalization of Bayesian Autoencoders . . . . . . . . . . . . . 5. 6. 2 The pathology of standard priors for Bayesian auto-encoders (b a es) and how to fix it . . . . . . . . . . . . . . . . . . . . . . . 5 . 7 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6. 3 . 1

 31 Figure 1.1 Evolution of accuracy and model complexity used for solving the ImageNet classification task . . . . . . . . . . . . . . . Figure 1.2 Autonomous driving agent learning to turn left while avoiding collisions, being being fast and smooth (image taken from the Tesla AI keynote 2021). . . . . . . . . . . . . . . . . .

Figure 1. 3

 3 Figure 1.3 Visualization of point-estimated solutions and predictions with a probabilistic model. . . . . . . . . . . . . . . . . . . . .

Figure 1. 4

 4 Figure 1.4 Uncertainty decomposition using an heteroskedastic model of the aleatoric uncertainty and the epistemic uncertainty for the Motorcycle dataset. . . . . . . . . . . . . . . . . . . . . . . . Figure 2.1 Graphical representation of the generalization interpretation of the marginal likelihood . . . . . . . . . . . . . . . . . . . . . Figure 2.2 Example of a generic deep feedforward neural networks with L hidden layers. For visualization purposes the bias term is included in the notation W (l) . . . . . . . . . . . . . . .

Figure 2. 3

 3 Figure 2.3Analysis of the distribution of samples from the prior for a shallow neural network with tanh activation. . . . . . . . . .

Figure 2. 4

 4 Figure 2.4 Pictorial representation of running Bayesian inference on a simple architectures. The posterior and the predictive distributions are computed with numerical integration . . . Figure 2.5 Pictorial illustration of the variational inference (vi) procedure, inspired by D. Blei et al. (2016). . . . . . . . . . . . . . . Figure 2.6 Animation of a simple variational inference procedure, using analytical gradients and Monte Carlo estimation. . . . . Figure 2.7 Graphical representation of the reparameterization trick. . .

Figure 3. 1 Figure 3 . 2

 132 Figure 3.1 Due to poor initialization v i fails to converge even after 1500+ epochs while with our iterative Bayesian linear modeling (i b l m) v i easily recovers the function after fewer iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3. 7

 7 Figure 3.7 Progression of test error and test m n l l with different initializations on LeNet for MNIST and CIFAR10. Note that the Uniformative, Heuristic and Xavier intializers did not converge, and as such they are not included in the plots. . .

Figure 3. 8

 8 Figure 3.8 Progression of test error and m n l l for two different convolutional neural networks on both MNIST and CIFAR10. Variational inference is initialized with iblm. . . . . . . . . . Figure 3.9 Reliability diagram and the expected calibration error for AlexNet trained on CIFAR10. . . . . . . . . . . . . . . . . . . . Figure 3.10 Progression of test error and mnll for a very deep convolutional neural network trained on CIFAR10. . . . . . . . . . Figure 3.11 Test likelihood after initialization (on the left) and at the end of the training procedure (on the right) for LeNet on MNIST.Figure 4.1 Normalized covariance of g and vect(W). . . . . . . . . . . .

Figure 4. 1

 1 Normalized covariance of g and vect(W). . . . . . . . . . . .

Figure 4. 2

 2 Figure 4.2 Numerical verification of first and second moments for Walsh-Hadamard variational inference (whvi) for one choice of S 1 , S 2 and q(g). The empirical quantities are obtained with Monte Carlo samples. . . . . . . . . . . . . . . . . . . . . Figure 4.3 Ablation study of different structures for the parameterization of the weights distribution. . . . . . . . . . . . . . . . . .

Figure 4. 5

 5 Figure 4.5Regression example trained using w h v i with Gaussian vector (1541 param.) and with planar normalizing flow (10 flows for a total of 4141 param.), mean-field Gussian (35k param.) and Monte Carlo dropout (mcd) (17k param.). . . .

Figure 4. 6

 6 Figure 4.6 Comparison of the test m n l l as a function of the number of hidden units. The dataset used is Powerplant. . . . . . . .

Figure 4. 7

 7 Figure 4.7 Comparison of the test m n l l as a function of the number of hidden units and hidden layers. . . . . . . . . . . . . . . .

Figure 4. 8

 8 Figure 4.8 Analysis of the training curve (on the left) and test curve (on the right) for the Mocap dataset. . . . . . . . . . . . . . .

Figure 4. 9

 9 Figure 4.9 Reliability diagram and expected calibration error of vgg , a l e x n e t and r e snet with whvi. . . . . . . . . . . . . . . Figure 4.10 Test curve of LeNet when trained with different combination of variational inference parameterization. . . . . . . . . . . . .

Figure 4. 11

 11 Figure 4.11 Analysis of the timing efficiency of whvi. (a) Inference time on the test set with 128 batch size and 64 Monte Carlo samples. Experiment repeated 100 times on an g p u fully dedicated to it. (b) Timing comparison of different implementations of the vectorized Fast Walsh-Hadamard transform, with a batch size of 512 data points. . . . . . . . . . . .

Figure 4. 12

 12 Figure 4.12 Comparison of the test mnll as a function of the memory footprint for four datasets. . . . . . . . . . . . . . . . . . . . .

Figure 4

 4 Figure 4.13Power profiling during inference on the test set of CIFAR10 with a l e x n e t The task is repeated 16 consecutive times

Figure 5. 1

 1 Test performance on CIFAR10 with different priors. Figure replicated from the arXiv version of A. G. Wilson and Izmailov (2020a). . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5. 2 Figure 5 . 3

 253 Figure 5.2 Grid search analysis on the prior variance for LeNet trained via variational inference on CIFAR10. . . . . . . . . . . . . . .

Figure 5. 4

 4 Figure 5.4 Schematic representation of the process of imposing g p priors on b nns via Wasserstein distance minimization. . . .

Figure 5. 5

 5 Figure 5.5 Visualization of one-dimensional regression example with a three hidden-layer multilayer perceptron (mlp). The first two rows illustrate the prior sample and distributions, whereas the last two rows show the corresponding posterior distributions. The means and the 95% credible intervals are represented by red lines and shaded areas, respectively. The middle row shows progressions of the prior optimization. . Figure 5.6 A timing comparison between imposing functional prior and cross-validation with grid-search and using Bayesian optimization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5. 7

 7 Figure 5.7 Different loss functions for model selection share the same optimum, for linear models. . . . . . . . . . . . . . . . . . . . .

Figure 5. 8

 8 Figure 5.8 Realizations sampled from different priors given an input image. o o d stands for out-of-distribution. . . . . . . . . . .

Figure 6. 1

 1 Figure 6.1 Representation of the induced distribution on the covariance function at location x when placing priors on different set of parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 6.2 Illustration of a binary classification task on the b a na na dataset. Left: the decision bounds of the average classifier. Right: the posterior marginals of the inducing inputs. . . . . Figure 6.3 Visual representation of a 2-layer deep Gaussian process (d g p). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 6. 4

 4 Figure 6.4 Analysis of different priors on inducing locations for Bayesian Sparse Gaussian Process (b s g p) on the UCI benchmark datasets for determinantal point process (dpp), Strauss process, uniform and normal priors on Z. . . . . . . . . . . . . .

Figure 6

 6 Figure 6.7 p-values of the hypothesis test that b s g p with variational free energy (v f e) objective is better than b s g p with fully independent training conditional (fitc) objective; depth of the d g p from 1 to 5. For models with p-values < 0.05, we reject the hypothesis . . . . . . . . . . . . . . . . . . . . . . . . Figure 6.8 bsgp in with different depths of the DGP with two different objective: fitc and vfe. The number of layers corresponds to the depth of the dgp. . . . . . . . . . . . . . . . . . . . . . .

Figure 6. 9

 9 Figure 6.9 Test m n l l on UCI regression benchmarks (the error bars represent the 95%CI). The lower m n l l (i.e. to the left), the better. The number on the right of the method's name refers to the depth of the dgp. Bottom right: Rank summary of all methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 6.10 Comparison of test mnll as function of training time. The dashed line on the right hand side plot corresponds to svgp with M = 1000 inducing points. . . . . . . . . . . . . . . . . . . Figure 6.11 Comparison of test m n l l as a function of prediction time on the largest dataset (Protein). . . . . . . . . . . . . . . . . . .

Figure 6. 12

 12 Figure 6.12 Comparison with structured inducing variables methods. k i s sg p could only run on the Powerplant dataset (hence the on Protein and Kin8NM). . . . . . . . . . . . . . . . . . . Figure 7.1 Number of unique new papers submitted on arxiv with title or abstract containing different permutations of Bayesian deep learning. Data crawled from the public API on December 8th 2021. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure b.1 Progression of test r m s eand test m n l l with different initializations on a shallow architecture. . . . . . . . . . . . . . .

For

Figure 1 . 1 :

 11 Figure 1.1: Evolution of accuracy and model complexity used for solving the ImageNet classi cation task. Data extracted in September 2021 from https://paperswithcode.com/ sota/image-classification-on-imagenet

Figure 1 . 2 :

 12 Figure 1.2: Autonomous driving agent learning to turn left while avoiding collisions, being being fast and smooth (image taken from the Tesla AI keynote 2021).

  Binary classi cation in two dimensions.

Figure 1 . 3 :

 13 Figure 1.3: Visualization of point-estimated solutions and predictions with a probabilistic model.

Figure 1 . 4 :

 14 Figure 1.4: Uncertainty decomposition using an heteroskedastic model of the aleatoric uncertainty and the epistemic uncertainty for the Motorcycle dataset.

  ,...,N . A generic formulation of the Bayes' rule takes the following form,

  Figure 2.1

Figure 2 . 1 :

 21 Figure 2.1: Graphical representation of the generalization interpretation of the marginal likelihood

Figure 2 . 2 :

 22 Figure 2.2: Example of a generic deep feedforward neural networks with L hidden layers. For visualization purposes the bias term is included in the notation W (l) .

Figure 2 . 4 :

 24 Figure 2.4: Pictorial representation of running Bayesian inference on a simple architecture y =

Figure 2 . 5 :

 25 Figure 2.5: Pictorial illustration of the variational inference ( ) procedure, inspired by D. Blei et al. (2016).

Figure 2 . 6 :

 26 Figure 2.6: Animation of a simple variational inference procedure, using analytical gradients and Monte Carlo estimation.

Figure 2 . 7 :

 27 Figure 2.7: Graphical representation of the reparameterization trick.

Figure 3 . 1 :

 31 Figure 3.1: Due to poor initialization (left) fails to converge even after 1500+ epochs (RMSE = 0.434, NLL = 257.38) while with our iterative Bayesian linear modeling ( ) (right) easily recovers the function after fewer epochs (RMSE = 0.219, NLL = -88.32). The architecture has three hidden layers with 500 neurons each, and uses the activation function.

)Figure 3 . 2 :

 32 Figure 3.2: Visual representation of the proposed method for initialization. In (a) and (b), we learn two Bayesian linear models, whose outputs are used in (c) to infer the following layer.

Figure 3 . 3 :

 33 Figure 3.3: Illustration of the di erent initialization for the variational parameters in a simple 1D regression task.

Algorithm 2 :

 2 Sketch of i b lm Inputs : Model M, Dataset D Returns : Initialized model foreach layer in M do foreach outfeature in layer do X, Y ← next batch in D; propagate X;

Figure 3 . 4 :

 34 Figure 3.4: Progression of test mean negative loglikelihood ( ) with di erent initializations with shallow and deep architectures on Powerplant and Protein. Experiment repeated ve times, only the average is shown.

Figure 3

 3 Figure3.5 (on the left) shows how test mnll just after the initialization is affected by this choice. Using the full training set leads to a better estimate of the posterior but from 64/128 samples the improvement on the test mnll is only marginal. Note also that the mini-batch size affects also the heterogeneity of the posteriors, which vanishes when using the full training set. The same experiment is also repeated comparing test mnll after initialization between vi with iblm and mcd with the common Xavier initializer (Figure3.5 on the right). Similar comments apply also for this case: i b l m allows the training to start from a marginally lower negative log-likelihood.

Figure 3 . 6 :

 36 Figure 3.6: On the left, comparison of initialization time versus test , averaged out of eight successive runs (on the right, magni cation of the small portion of the plot). Orange corresponds to the Pareto frontier. Before training, four out of ve optimal initializers are .

Figure 3 . 10 :

 310 Figure 3.10: Progression of test error and for a very deep convolutional neural network trained on CIFAR10.

  [START_REF] Zhang | Noisy Natural Gradient as Variational Inference[END_REF]. The four models are trained with a maximum time budget of 100 minutes for the entire end-to-end training (curves are shifted by the initialization time). The Results are shown in Figure3.10 and in the adjacent table. In this experiment, we have experienced the situation in which, the e l b o is completely dominated by the k l divergence. The k l regularization term in the variational objective severely penalizes training of over-parameterized model. With a sensible

Figure 3 . 11 :

 311 Figure 3.11: Test likelihood after initialization (on the left) and at the end of the training procedure (on the right) for LeNet on MNIST.

  Carlo dropout (m c d);Gal and Ghahramani, 2016b) and Noisy Natural Gradients (G.[START_REF] Zhang | Noisy Natural Gradient as Variational Inference[END_REF].

  Without loss of generality, consider Bayesian d n ns with weight matrices W(l) of size D × D. Compared with mean field v i, w h v i has a number of attractive properties. The number of parameters is reduced from O(D 2 ) to O(D), thus reducing the over-regularization effect of the kl term in the variational objective. We derive expressions for the reparameterization and the local reparameterization tricks, showing that, the computational complexity is reduced from O(D 2 ) to O(D log D).

Figure 4 . 1 :

 41 Figure 4.1: Normalized covariance of g and vect(W).

  H) :,i , so W can be rewritten in terms of A ∈ R D 2 ×D and g as follows vect(W) = Ag where A =

Figure 4 . 2 :

 42 Figure 4.2: Numerical veri cation of rst and second moments for

  4.2. For the three moments the errors are well within what we can consider acceptable and, hence we can validate the expressions we derived before.

Figure 4 . 4 :

 44 Figure 4.4: Comparison between Hadamard factorization in and tensor factorization. The number in the parenthesis is the hidden dimension. Plot is w.r.t. iterations rather then time to avoid implementation artifacts. The dataset used is Drive.

Figure 4 . 5 :

 45 Figure 4.5: Regression example trained using with Gaussian vector (1541 param.) and with planar normalizing ow (10 ows for a total of 4141 param.), mean-eld Gussian (35k param.) and Monte Carlo dropout ( ) (17k param.). The two shaded areas represent the 95th and the 75th percentile of the predictions. As "ground truth", we also show the predictive posterior obtained by running Stochastic Gradient Hamiltonian Monte Carlo (

Figure 4 . 6 :

 46 Figure 4.6: Comparison of the test as a function of the number of hidden units. The dataset used is Powerplant.

Figure 4 . 8 :

 48 Figure 4.8: Analysis of the training curve (on the left) and test curve (on the right) for the Mocap dataset.

Figure 4 . 10 :

 410 Figure 4.10: Test curve of LeNet when trained with di erent combination of variational inference parameterization.

Figure 4 . 11 :Figure 4 . 12 :

 411412 Figure 4.11: Analysis of the timing e ciency of . (a) Inference time on the test set with 128 batch size and 64 Monte Carlo samples. Experiment repeated 100 times on an fully dedicated to it. (b) Timing comparison of di erent implementations of the vectorized Fast Walsh-Hadamard transform, with a batch size of 512 data points.

Figure 4 . 13 :

 413 Figure 4.13: Power pro ling during inference on the test set of CIFAR10 withThe task is repeated 16 consecutive times and pro ling is carried out using the nvidia-smi tool.

Figure 4 . 14 :

 414 Figure 4.14: Analysis of the parameter e ciency of for scalable Gaussian process regression with variational inference.

Figure 5 . 1 :

 51 Figure 5.1: Test performance on CIFAR10 with di erent priors. Figure replicated from the arXiv version of A. G. Wilson and Izmailov (2020a).

Figure 5 . 2 :

 52 Figure 5.2: Grid search analysis on the prior variance for LeNet trained via variational inference on CIFAR10.

Figure 5

 5 Figure 5.3: (Top) Sample functions of a fully-connected with 2, 4 and 8 layers obtained by placing a Gaussian prior on the weights. (Bottom) Samples from a prior with two di erent kernels.

Figure 5 . 4 :

 54 Figure 5.4: Schematic representation of the process of imposing priors on s via Wasserstein distance minimization.

Figure 5 . 5 :

 55 Figure 5.5: Visualization of one-dimensional regression example with a three hidden-layer multilayer perceptron ( ). The rst two rows illustrate the prior sample and distributions, whereas the last two rows show the corresponding posterior distributions. The means and the 95% credible intervals are represented by red lines and shaded areas, respectively. The middle row shows progressions of the prior optimization.

  gap in the dataset;

Figure 5 . 6 :

 56 Figure 5.6: A timing comparison between imposing functional prior and cross-validation with gridsearch () and using Bayesian optimization ( ). In the plots, each corresponds to a run of a single con guration, while highlights the Pareto front of the cross

Figure 5 . 7 :

 57 Figure 5.7: Di erent loss functions for model selection share the same optimum, for linear models.

Bayes

  y | w) is the conditional likelihood that factorizes as p(y | w) = N n=1 p(y n | w).

Figure 5 . 8 :

 58 Figure 5.8: Realizations sampled from di erent priors given an input image.stands for out-of-distribution.

R

  E V I S I T I N G T H E A P P R O X I M A T I O N S F O R S C A L A B L E ( D E E P ) G A U S S I A N P R O C E S S E S

Figure 6 . 1 :

 61 Figure 6.1: Representation of the induced distribution on the covariance function at location x when placing priors on di erent set of parameters.

  {X, Z}, respectively. We assume a factorized likelihood p(y| f) = N n=1 p(y n | f n )

  s inequality to bound the log joint marginal we obtain the following formulation,log p(y, Ψ | X) E p(f | Ψ,X) log p(y | f) + log p(Ψ) def = log p vfe (y, Ψ | X).

Figure 6 . 2 :

 62 Figure 6.2: Illustration of a binary classi cation task on the dataset. Left: the decision bounds of the average classi er. Right: the posterior marginals of the inducing inputs.

Figure 6 . 3 :

 63 Figure 6.3: Visual representation of a 2-layer deep Gaussian process ().

Figure 6 . 4 :

 64 Figure 6.4: Analysis of di erent priors on inducing locations for on the UCI benchmark datasets for determinantal point process ( ), Strauss process, uniform and normal priors on Z.

Figure 6 . 6 :

 66 Figure 6.6: Analysis of di erent choices of objectives when used for optimization and sampling.We refer the reader to Table6.1 for a description of the methods.

Figure 6 . 9 :

 69 Figure 6.9: Test on UCI regression benchmarks (the error bars represent the 95%CI). The lower (i.e. to the left), the better. The number on the right of the method's name refers to the depth of the . Bottom right: Rank summary of all methods.

Figure 6 . 10 :

 610 Figure 6.10: Comparison of test as function of training time. The dashed line on the right hand side plot corresponds to with M = 1000 inducing points.

Figure 6 . 11 :

 611 Figure 6.11: Comparison of test as a function of prediction time on the largest dataset (Protein).

Figure 7 . 1 :

 71 Figure 7.1: Number of unique new papers submitted on arxiv with title or abstract containing di erent permutations of Bayesian deep learning. Data crawled from the public API on December 8th 2021.

  f(z) = z + uh(w z + b) , (a.4)where λ = [u ∈ R D , w ∈ R D , b ∈ R] and h(•) = tanh(•). This is equivalent to a residual layer with single neuron m l p -as argued by[START_REF] Kingma | Improved Variational Inference with Inverse Autoregressive Flow[END_REF]. The log-determinant of the Jacobian of f is log det∂f ∂z = det I + u[h (w z + b)w] = 1 + u wh (w z + b) . (a.5)Alternatives can be found in recent works by D.[START_REF] Rezende | Variational Inference with Normalizing Flows[END_REF],[START_REF] Van Den Berg | Sylvester Normalizing Flows for Variational Inference[END_REF],[START_REF] Kingma | Improved Variational Inference with Inverse Autoregressive Flow[END_REF], andLouizos and Welling (2017).

W

  p (π, ν) = inf γ∈Γ (π,ν) X×Y D(x, y) p γ(x, y) dx dy

Theorem 2

 2 Given a nonnegative measurable function D : X × X → R, the 1-Wasserstein is computed as follows, f(x)+g(y) D(x,y) f(x)π(x) dx + g(y)ν(y) dy .

  because given the definition of Γ , if γ ∈ Γ (π, ν) then π(x) = γ(x, y) dy and ν(y) = γ(x, y) dx. By substituiting these quantities, it follows that f(x)π(x) dx + g(y)ν(y) dy = f(x) γ(x, y) dy dx + g(y) γ(x, y) dx dy = [f(x) + g(y)] γ(x, y) dx dy .(a.14)In other cases, f and g can be chosen such that the supremum becomes +∞. Given this property and the constrain on γ, we can add ι Γ (γ) to the formulation of the Wasserstein distance in Equation (a.8),W(π, ν) = inf γ∈Γ (π,ν) D(x, y)γ(x, y) dx dy + ι Γ (γ) = = inf γ D(x, y)γ(x, y) dx dy + sup f,g f(x)π(x) dx + g(y)ν(y) dy-[f(x) + g(y)] γ(x, y) dx dy , (a.15)Now, the original integral of the Wasserstein distance does not depend on f and g;therefore the supremum can be moved in front, x, y)γ(x, y) dx dy + f(x)π(x) dx + g(y)ν(y) dy-[f(x) + g(y)] γ(x, y) dx dy (a.16)Under certain conditions stated by the minimax theorem, i.e. Υ(γ, (f, g)) is convexconcave function (Υ is concave for fixed (f, g) while convex for fixed γ ), we can swap the infinum and the supremum and rewrite the definition as follows, x, y)f(x)g(y)] γ(x, y) dx dy + f(x)π(x) dx + g(y)ν(y) dy (a.17)

  x, y)f(x)g(y)] γ(x, y) dx dy =    0 if f(x) + g(y) D(x, y) , -∞ otherwise . (a.18) If the function ζ(x, y) = D(x, y) -(f(x) + g(y)

  , we observe that for a given f, if we want the supremum in Eq. 5 we cannot get a better g then taking the equality, f(y) := inf x [D(x, y)f(x)] .(a.21)We therefore have the following formulation of the Wasserstein distance,W(π, ν) = sup f f(x)π(x) dx + f(y)ν(y) dy (a.22)If now we suppose to choose g, by following the same reasoning the best f that we can get is definedf(x) = ḡ(x) := inf y [D(x, y)g(y)] . (a.23)If we replace g(y) with Eq. 17 we have yet another recursive definition of the Wasserstein distance, W(π, ν) = sup f f(x)π(x) dx + f(y)ν(y) dy (a.24)

  y)γ(x, y) dx dy = sup f,g f(x)+g(y) D(x,y) f(x)π(x) dx + g(y)ν(y) dy , +g(y) D(x,y) f(x)π(x) dx + g(y)ν(y) dy = = sup f f=inf x D-f f(x)π(x) dx + f(y)ν(y) dy ,

  we compare initialization methods for a shallow deep neural network (d n n) architecture on two datasets. The architecture used in these experiments has one single hidden layer with 100 hidden neurons and rectified linear unit (R e L U) activations. We impose that the approximate posterior has fully factorized covariance. Figure b.1 shows the learning curves on the Powerplant (n = 9568, d = 4) and Protein (n = 45730, d = 9) datasets, repeated over five different train/test splits. iterative Bayesian linear modeling (i b l m) allows for a better initialization compared to the competitors, leading to a lower root mean square error (r m s e) and lower mean negative log-likelihood (mean negative loglikelihood (mnll)) on the test for a given computational budget. We refer the reader to the supplementary material for a more detailed analysis of the results. Using the same deep dnn architecture as in the last experiment (five hidden layers with 100 neurons), we tested i b l m with classification problems on MNIST (n = 70000, d = 784), EEG (n = 14980, d = 14), Credit (n = 1000, d = 24) and Spam (n = 4601, d = 57). Interestingly, with this architecture, some initialization strategies struggled to converge, e.g., uninformative on MNIST and lsuv on EEG . The gains offered by iblm achieves are most apparent on MNIST . After less than 1000 training steps (less than an epoch), i b l m makes variational inference (v i) reach a test accuracy greater than 95%; other initializations reach such performance much later during training. Even

  Figure b.1: Progression of test and test with di erent initializations on a shallow architecture.

  Figure b.2: Progression of test and test with di erent initializations on a deep architecture.

  Practical considerations and extensions to deep GPs . . . . . . . . . .

6 . 1 Sparse Gaussian processes . . . . . . . . . . . . . . . . . . . . . . . . . . 6 . 2 Bayesian Sparse Gaussian Processes . . . . . . . . . . . . . . . . . . . . 6. 2. 1 On scalable inference frameworks for GP models . . . . . . . . 6. 2. 2 Sampling with VFE or FITC? . . . . . . . . . . . . . . . . . . . . 6. 2. 3 Stochastic Updates Using the FITC Approximation . . . . . . 6. 2. 4 An heteroskedastic version of the Gaussian likelihood . . . . 6. 2. 5 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . 6 . 3

Table 4 .

 4 1 Time and space complexity of various approaches to variational inference (v i). . . . . . . . . . . . . . . . . . . . . . . . . Table 4.2 Test root mean squared error (rmse) and test mean negative

loglikelihood (mnll) for regression datasets. Results in the format "mean (std)" . . . . . . . . . . . . . . . . . . . . . . . . .

Table 4.3 Test performance of different Bayesian convolutional neural networks (c n ns). . . . . . . . . . . . . . . . . . . . . . . . . . . Table 5.1 Glossary of methods used in the experimental campaign. . . Table 5.2 Results for different convolutional neural networks on CI-FAR10. Experimental comparison performed in B.-H. Tran et al. (2020) by the first author (Ba-Hien Tran). . . . . . . . . Table 6.1 A summary of previous works on inference methods for Gaussian processs (g ps). θ, u, Z refer to the g p hyperparameters, inducing variables and inducing inputs, respectively. () indicates that variables are optimized. . . . . . Table 6.2 a i r l i n e dataset predictive test performance. . . . . . . . . . Table 6.3 h i g g s dataset predictive test performance. . . . . . . . . . .

  (mcmc) can be efficiently implemented in Bayesian neural networks, specifically by using the Hamiltonian Monte Carlo (hmc) framework introduced by Duane et al.

	1 . 3

discuss how Markov chain Monte Carlo

  N × D out . Each part of this equation has a specific role. The

	2 .1
	will see several examples of this in the rest of this thesis. The quantity p(y | f(X; θ))
	fixes our assumptions on the observed data (e.g. continuous and noisy, binary,
	count, etc.) and it measures how likely it is for the model f with parameters θ to
	have generated the observed values. This term-the likelihood-is the probability
	density function of the model given the parameters computed at the observed
	data. Examples of likelihoods are the Gaussian likelihood N(y | f(x, θ), σ 2 ) useful for
	regression tasks with i.i.d. noisy data, or the Bernoulli likelihood Bern (y | λ(f(x, θ))),
	for binary classification (with λ : R → [0, 1]).
	prior p(θ) specifies our belief on the parameters before observing any data. This
	should capture any relevant information on the problem (if any) or being totally
	uninformative by letting "the data speak for itself". Such kind of priors (a. k. a. un-
	informative priors) are usually also improper as they don't integrate to one over their
	domain; examples are the Haldane prior (Zhu and A. Y. Lu, 2004) and the Jeffreys
	prior (Jeffreys, 1946). Note something very important: the Bayesian philosophy of
	reasoning about priors doesn't explain how to choose them (i.e., any prior can work
	with any model). In practice, some prior choices perform better than others, and we

  Progression of test error and test with di erent initializations on LeNet for MNIST and CIFAR10. Note that the Uniformative, Heuristic and Xavier intializers did not converge, and as such they are not included in the plots.

							3 . 3		47
				MNIST		0.8		CIFAR10
	Test error	0.05 0.1 0.15			Test error	0.4 0.6
		0	10 2	10 3	10 4		10 2	10 3	10 4
						2
		0.4				1.8
	Test MNLL	0.2			Test MNLL	1.4 1.6
						1.2
		0	10 2	10 3	10 4	1	10 2	10 3	10 4
				Iteration				Iteration
					Orthogonal
	Figure 3.7:

  Reliability diagram and the expected calibration error for AlexNet trained on CIFAR10.

							3 . 3	49
		1	LeNet (MNIST) Calibration on CIFAR10		AlexNet (CIFAR10)
		0 0.2 0.8 0.4 Test error rate 0 0 0.2 0.4 0.6 Accuracy	0.2	0.4 Con dence 0.6	0.8	1	0.8 0.4 0.6 2 Mean-eld Gaussian (with IBLM) ECE = 0.064 MCD ECE = 0.210
	1 Figure 3.8: Progression of test error and 10 1 10 2 10 3 10 4 10 5 0 0.2 0.4 Iteration Test Mean-eld Gaussian (with IBLM) 10 6 10 1 1 1.2 1.4 1.6 1.8 for two di erent convolutional neural networks 10 2 10 3 10 4 10 5 Step MCD on both MNIST and CIFAR10. Variational inference is initialized with . VGG16 (CIFAR10) Figure 3.9: 0.2 0.8 0.4 0.6 Test Error 2 Test MNLL Test error Gaussian VI with IBLM 0.637 0.167 Gaussian VI with MAP init. 0.750 0.201 MCD 0.821 0.215 NNG (Noisy-KFAC) 0.750 * 0.164
	Test MNLL	1 1.5					Note: Due to high variance in the test natural gradients, we report the minimum value of the last 5 for the noisy epochs.
		2000		4000	6000
			Seconds
		Diagrams are a visualization tool where sample accuracy is plotted as function of
		confidence (DeGroot and Fienberg, 1983; Niculescu-Mizil and Caruana, 2005). For
		a perfectly calibrated model, the diagram follows the identity function. Expected

https://github.com/yaringal/DropoutUncertaintyCaffeModels λ = 0.01 and weight decay of 0.0005. Figure

3

.8 on the left shows the learning curves: m c d achieves lower error rate but its rough approximation imposed by the Bernoulli-like posterior is reflected on an higher mnll compared to vi with a Gaussian posterior. Provided with a good initialization, Gaussian vi can better fit the model and deliver a better quantification of uncertainty.

Similar comments apply also for CIFAR10 on a l e x n e t

[START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF]

, a c n n composed by a stack of five convolutional layers and three fully-connected layers for a total of more than 1M parameters (2M for v i) (results in Figure

3

.8 on the right). Note that in this case, the test error seems to favor mcd but the mean-field Gaussian approximation deliver better test likelihood (1.15 vs 1.38 for m c d). Calibration of uncertainty is an important performance metric that one should take into account for comparing classification models

[START_REF] Flach | Classifier Calibration[END_REF][START_REF] Guo | On Calibration of Modern Neural Networks[END_REF]

. Reliability Diagrams and the Expected Calibration Error (e c e) are standard methods to empirically estimate the calibration uncertainty. Reliability

  Ablation study of di erent structures for the parameterization of the weights distribution. Metric: test error rate and test mean negative loglikelihood ( ) with di erent structures for the weights. Benchmark on Drive with a 2 × 64 network. Colors are coded to match the ones used in the Figure

				4 .4				,			63
	0.25			Test Error
	0.20
	0.15
	0.05 0.10							Model		Test Error Test MNLL
	0.00 0.4 0.6 0.8 1.0	0	5000	10000 Test MNLL 15000	20000	25000	GH S var HGH S 1,var HGHS 2,var H S opt HGH S 1,opt HGHS 2,opt H S 1,opt HGHS 2,opt (	)	0.097 0.226 0.043 0.061 0.054 0.031 0.026	0.249 0.773 0.159 0.190 0.199 0.146 0.094
	0.2
	0.0	0	5000	10000	15000	20000	25000
	Figure 4.3:

Table 4 . 1 :

 41 Time and space complexity of various approaches to . Note: D is the dimensionality of the feature map, K is the number of tensor cores, R is the rank of tensor cores and

  Function SetupDimensions(D in , D out ): next power ← 2 ceillog 2 D in ; if next power == 2D in then padding ← 0; else padding = next power -D in ; D in ← next power; stack, remainder = divmod(D out , D in ); if remainder != 0 then stack ← stack + 1; D out ← D in × stack; return D in , D out , padding, stack

  More complex alternatives are proposed by[START_REF] Van Den Berg | Sylvester Normalizing Flows for Variational Inference[END_REF][START_REF] Kingma | Improved Variational Inference with Inverse Autoregressive Flow[END_REF] Louizos and Welling, 2017).

Table 4 .

 4

	2: Test root mean squared error ( in the format "mean (std)"	) and test	for regression datasets. Results
	Model Dataset	MCD	MFG	NNG	Test error WHVI	MCD	MFG	NNG	Test MNLL WHVI
	Boston Concrete	3.91 (0.86) 4.47 (0.85) 3.56 (0.43) 3.14 (0.71) 6.90 (2.93)	2.99 (0.41) 2.72 (0.09) 4.33 (1.80)

Table 4

 4

	.3: Test performance of di erent Bayesian s. CIFAR10 Test error Test MNLL		0.8 Calibration of the predictions on CIFAR10 1	6000
	VGG16 MFG MCD Noisy-KFAC WHVI AlexNet MFG MCD Noisy-KFAC WHVI NF-WHVI	16.82% 21.47% 15.21% 12.85% -13.30% 20.36% 13.56% 12.72%	0.6443 0.8213 0.6374 0.6995 -0.9590 1.1990 0.6164 0.6596	Accuracy	0.2 0.4 0.6 0	0 VGG16 (ECE = 0.0229) 0.2 0.4 Con dence 0.6 AlexNet (ECE = 0.0290)	4000 2000 0 ResNet18 (ECE = 0.0088) 0.8 1 Perfect calibration
	ResNet18 MFG MCD Noisy-KFAC WHVI NF-WHVI	-10.71% -11.46% 11.42%	-0.8468 -0.4908 0.5513	Figure 4.9: Reliability diagram and expected cal-ibration error of , and with .

Table 5 . 1 :

 51 Glossary of methods used in the experimental campaign. Here, p(f) = p(f | w) dp(w) denotes the induced prior over functions; Γ -1 (α, β) denotes the Inverse-Gamma distribution with shape α, and rate β; NF(T K ) indicates a normalizing ow distribution constructed from a sequence of K invertible transformations

Table 5 . 2 :

 52 Results for di erent convolutional neural networks on CIFAR10. Experimental comparison performed in B.-H. Tran et al. (2020) by the rst author (Ba-Hien Tran).

	Architecture Method	Accuracy -% (↑)	NLL (↓)
	LeNet	Deep Ensemble	71.13 ± 0.10	0.8548 ± 0.0010
		prior + i-prior (ours)	74.65 ± 0.25 74.08 ± 0.24 75.15 ± 0.24	0.7482 ± 0.0025 0.7558 ± 0.0024 0.7360 ± 0.0024
		prior i-prior (ours)	75.22 ± 0.40 76.51 ± 0.21	0.7209 ± 0.0040 0.6952 ± 0.0021
	PreResNet Deep Ensemble	87.77 ± 0.03	0.3927 ± 0.0003
		prior + i-prior (ours)	85.34 ± 0.13 87.70 ± 0.11 86.86 ± 0.27	0.4975 ± 0.0013 0.3956 ± 0.0011 0.4286 ± 0.0027
		prior i-prior (ours)	87.26 ± 0.09 88.20 ± 0.07	0.4086 ± 0.0009 0.3808 ± 0.0007
	VGG	Deep Ensemble	81.96 ± 0.33	0.7759 ± 0.0033
		prior + i-prior (ours)	81.47 ± 0.33 82.25 ± 0.15 83.34 ± 0.53	0.5808 ± 0.0033 0.5398 ± 0.0015 0.5176 ± 0.0053
		prior i-prior (ours)	86.03 ± 0.20 87.03 ± 0.07	0.4345 ± 0.0020 0.4127 ± 0.0007

  in the context of highly-structured data, which is ubiquitous in modern machine learning applications. Consequently, we have proposed to specify priors over the autoencoder weights by means of a novel optimization of prior hyper-parameters. Inspired by connections with marginal likelihood optimization, we derived a practical and efficient optimization framework, based on the minimization of the distributional sliced-Wasserstein distance between the distribution induced by the bae and the data generating distribution. The resulting hyper-parameter optimization strategy leads to a novel way to perform model selection for baes. Note that, even if theoretically justified and empirically verified with extensive experimentation, our proposal for model selection still remains a proxy to the true marginal likelihood maximization. The d s w d formulation has nice properties of asymptotic convergence and computational tractability, but it may represent only one of the possible solutions. At the same time, we stress that the current literature does not cover this problem of b a es at all, and we believe our approach is a considerable step towards the development of practical Bayesian methods for representation learning in modern applications characterized by largescale structured data (including tabular and graph data, which are currently not covered).

	5 . 7	105
	specifying meaningful priors

Table 6 .

 6 1: A summary of previous works on inference methods for s. θ, u, Z refer to the hyper-parameters, inducing variables and inducing inputs, respectively. () indicates that variables are optimized.

					Inference
		Model		θ	u	Z Reference
	( ) ( ) ( ) ( ) ( )	--	--	-	-(heterosk.)	-	Neal (1997) and Barber and Williams (1997) Hensman, A. Matthews, et al. (2015) Titsias (2009b) Havasi et al. (2018) H. Yu et al. (2019) Hensman, A. G. Matthews, et al. (2015)
	( )						This work

;

  Shi

		6.1	109
	RBF kernel k(x, z i ) with xed λ, σ	Distribution on k(x, z i ) when sampling {λ, σ}
	z i	z i
	Distribution on k(x, z i ) when sampling {Z}	Distribution on k(x, z i ) when sampling {Z, λ, σ}
	z i	z i

  ., 2013; on the variational objective (SVGP -Hensman, A.Matthews, et al., 2015) Inference of u with heteroskedastic likelihood (equivalent to FITC) (FITC-SVGP -Titsias, 2009b) MCMC inference of u, θ on the variational objective (MCMC-SVGP -Hensman, A. G.Matthews, et al., 2015) MCMC inference of u, θ on the marginal likelihood (log of expectation) MCMC inference of u, θ, Z on the marginal likelihood[This work]

							6.4	123
		Boston	3.2	Concrete	1.8	Energy	Kin8NM
	Test MNLL	2.4 2.6	3.05 3.1 3.15		1 1.2 1.4 1.6	-1.1 -1.05
		Naval		Powerplant		Protein	Yacht
		-7	2.8		2.9	1
	-7.5	2.75		2.85	0.5
		-8	2.7		2.8	0
	Inference of u

. 2

We implemented the same architecture as in G.[START_REF] Zhang | Noisy Natural Gradient as Variational Inference[END_REF]

N(µ, Σ) of Equation (4.2). From here we can actually obtain a class of approximate

Which is, in fact, the same as in the sparse Gaussian process (s p g p) framework of E. L.[START_REF] Snelson | Flexible and efficient Gaussian process models for machine learning[END_REF].

.

We consider a layer-wise factorization with two independent zero-mean Gaussian distributions for weights and biases. The parameters to adjust

, where σ 2 l w is the prior variance shared across all weights in layer l, and σ 2 l b is the respective variance for the bias parameters. For any weight and bias entries w l , b l ∈ w l of the l-th layer, the prior is:

We refer to this case as the g p-induced b n n prior with Gaussian weights (g pig). Although this simple approach assumes a Gaussian prior on the parameters, in many cases it is sufficient to capture the target g p-based functional priors.

. A more flexible family of priors for b n ns considers a hierarchical structure where the network parameters follow a conditionally Gaussian distribution, and the prior variance for each layer follows an Inverse-Gamma distribution. For the weight and bias variances we have:

In this case, we have ψ = {α l w , β l w , α l b , β l b } L l=1 , where α l w , β l w , α l b , β l b denote the shape and rate parameters of the Inverse-Gamma distribution for the weight and biases correspondingly for layer l. The conditionally Gaussian prior over the network parameters is given as in the previous section. We refer to this parametrization as the g p-induced b n n prior with Hierarchically-distributed weights (g pi-h).

. Finally, we also consider normalizing flows (n fs) as a family of much more flexible distributions. By considering an invertible, continuous and differentiable function t : R D l → R D l , where D l is the number of parameters for l-th layer, a n f is constructed as a sequence of K of such transformations T K = {t 1 , . . . , t K } of a simple known distribution (e.g. Gaussian). Sampling from such distribution is as simple as sampling from the initial distribution and then apply the set of transformation T K . Given an initial distribution p 0 (w l ), by denoting p(T K (w l )) the final distribution, its log-density can