Bayesian learning for neural networks: an algorithmic survey

The last decade witnessed a growing interest in Bayesian learning. Yet, the technicality of the topic and the multitude of ingredients involved therein, besides the complexity of turning theory into practical implementations, limit the use of the Bayesian learning paradigm, preventing its widespread adoption across different fields and applications. This self-contained survey engages and introduces readers to the principles and algorithms of Bayesian Learning for Neural Networks. It provides an introduction to the topic from an accessible, practical-algorithmic perspective. Upon providing a general introduction to Bayesian Neural Networks, we discuss and present both standard and recent approaches for Bayesian inference, with an emphasis on solutions relying on Variational Inference and the use of Natural gradients. We also discuss the use of manifold optimization as a state-of-the-art approach to Bayesian learning. We examine the characteristic properties of all the discussed methods, and provide pseudo-codes for their implementation, paying attention to practical aspects, such as the computation of the gradients.


Introduction
Machine Learning (ML) techniques have been proven to be successful in many prediction and classification tasks across natural language processing [Young et al., 2018], computer vision [Krizhevsky et al., 2012], time-series [Längkvist et al., 2014] and finance applications [Dixon et al., 2020], among the several others.The widespread of ML methods in diverse domains is found due to their ability to scale and adapt to data, and their flexibility in addressing a variety of problems while retaining high predictive ability.Recently, Bayesian methods have gained considerable interest in ML as an attractive alternative to the classical methods providing point estimations for their inputs.Despite the numerous advantages that traditional ML methods offer, they are, broadly speaking, prone to overfitting, dimming their generalization capabilities and performance on unseen data.Furthermore, an implicit consequence of the classical point estimation and modeling setup is that it delivers models that are generally incapable of addressing uncertainties.This inability is twofold, as it includes both the estimation and prediction aspects.Indeed, as opposed to the typical practice of statistical modeling and, e.g., econometrics methods, ML methods do not directly tackle aspects related to the significance and uncertainties associated with the estimated parameters.At the same time, predictions correspond to simple point estimates without reference to the confidence levels that such estimates have.Whereas some models have been developed to, e.g., provide confidence intervals over the forecasts [e.g.Gal and Ghahramani, 2016], it has been observed that such models are generally overconfident.To estimate uncertainties implicitly embedded in ML models, Bayesian inference provides an immediate remedy and stands out as the main approach.
Bayesian methods have gained considerable interest as an attractive alternative to point estimation, especially for their ability to address uncertainty via posterior distribution, generalize while reducing overfitting [Hoeting et al., 1999], and for enabling sequential learning [Freitas et al., 2000] while re-Preprint.Under review.arXiv:2211.11865v4[stat.ML] 31 Jan 2023 taining prior and past knowledge.Although Bayesian principles have been proposed in ML decades ago [e.g.Mackay, 1992, 1995, Lampinen and Vehtari, 2001], it has been only recently that fast and feasible methods boosted a growing use of Bayesian methods in complex models, such as deep neural networks [Osawa et al., 2019, Khan et al., 2018a, Khan and Nielsen, 2018].
The most challenging task in following the Bayesian paradigm is the computation of the posterior.In the typical ML setting characterized by a high number of parameters and a considerable size of data, traditional sampling methods are prohibitive, and alternative estimation approaches such as Variational Inference (VI) have been shown to be suitable and successful [Saul et al., 1996, Wainwright and Jordan, 2008, Hoffman et al., 2013, Blei et al., 2017].Furthermore, recent research advocates the use of natural gradients for boosting the optimum search and the training [Wierstra et al., 2014], enabling fast and accurate Bayesian learning algorithms that are scalable and versatile.
Recent years witnessed enormous growth in the interest related to Bayesian ML methodologies and several contributions in the field.This survey aims at summarizing the major methodologies nowadays available, presenting them from an algorithmic, empirically-oriented perspective.With this rationale, this paper aims to provide the reader with the basic tools and concepts to understand the theory behind Bayesian Deep Learning (DL) and walk through the implementation of the several Bayesian estimation methodologies available.We should note that the focus of this paper is purely on Bayesian methods.Indeed there are a number of network architectures that can resemble a Bayesian framework by, e.g., creating a distribution for the outputs, e.g., Deep Ensembles [Osband et al., 2018], Batch Ensembles [Wen et al., 2020], Layer Ensembles [Oleksiienko and Iosifidis, 2022], or Variational Neural Networks [Oleksiienko et al., 2022].These solutions, based on particular network designs, are, however, not implicitly Bayesian and out of scope in our context.Other surveys and tutorials do exist on the general topic [e.g.Jospin et al., 2022, Heckerman, 2008, along with several lecture notes available online], yet the focus of this paper is on algorithms and mainly devoted to Variational Inference (VI) methods.In fact, despite the wide number of VI and non-VI methods published in the last decade, a comprehensive survey embracing and discussing all of them (or perhaps the major ones) is missing, and non-experts will easily find themselves lost in their pursuit to comprehend and different notions and processing steps in different methodologies.By filling this gap, we aim to promote applications and research in this area.

The Bayesian paradigm
The Bayesian paradigm in statistics is often opposed to the pure frequentist paradigm, a major area of distinction being in hypothesis testing [Etz et al., 2018].The Bayesian paradigm is based on two simple ideas.The first is that probability is a measure of belief in the occurrence of events, rather than just some limit in the frequency of occurrence when the number of samples goes toward infinity.The second is that prior beliefs influence posterior beliefs [Jospin et al., 2022].The above two are summarized in the Bayes theorem, which we now review.Let D denote the data and p(D|θ) the likelihood of the data based on a postulated model with θ ∈ Θ a k-dimensional vector of model parameters.Let p(θ) be the prior distribution on θ.The goal of Bayesian inference is the estimation of the posterior distribution [e.g., Gelman et al., 1995] where p(D) is referred to as evidence or marginal likelihood, since p(D) = Θ p(θ)p(D|θ)dθ.p(D) acts as a normalization constant for retrieving an actual probability distribution for p(θ|D).In this light, as opposed to the frequentist approach, it becomes clear that the unknown parameter θ is treated as a random variable.The prior probability p(θ), which intuitively expresses in probabilistic terms any knowledge about the parameter before the data has been collected, is updated in the posterior probability p(θ|D), mixturing prior knowledge and evidence supported by the data through the model's likelihood.Bayesian inference is generally difficult due to the fact that the marginal likelihood is often intractable and of unknown form.Indeed, only for a limited class of models, where the prior is so-called conjugate to the likelihood, the calculation of the posterior is analytically tractable.Standard examples are Normal likelihoods and prior (resulting in Normal posteriors) or Poisson likelihoods with Gamma priors (resulting in Negative Binomial posteriors).Yet, already for the simple linear regression example, Bayesian derivation is rather tedious, and already for the logistic regression, closed-form solutions are not generally available.It is clear that in complex models, such as deep neural networks typically used in ML applications, Bayesian inference can be tackled neither analytically nor numerically (consider that the integral in the marginal likelihood is multivariate, over as many dimensions as the number of parameters).
Monte-Carlo (MC) methods for sampling the posterior are certainly a possibility that has been early explored and adopted.While it still remains a valid and appropriate method for performing Bayesian inference in retractable settings, especially in high-dimensional applications, the MC approach is challenging and may become infeasible, mainly because of the need for an implicit high-dimensional sampling scheme, which is generally time-consuming and computationally demanding.As an alternative approach, VI gained much attention in recent years.VI turns the integration Bayesian problem in eq. ( 1) into an optimization problem.The idea behind VI is that of targeting an approximate form of the posterior distribution, perhaps chosen within a family of well-behaved distributions, and finding the corresponding parameter that optimizes a specific objective, i.e., that is optimal under some criterion.
In the following subsection, we review the standard non-Bayesian approach for neural network parameter estimation (Section 1.2.1), we introduce Bayesian Neural Networks (BNNs) (Section 1.4), and we provide some motivation for their use, also recalling some literature about their applications (Section 1.3).After providing the reader an introduction to standard and Bayesian neural networks, we introduce VI in Section 1.5, we describe the standard framework used in Bayesian learning, and we discuss how the standard Stochastic Gradient Descent (SGD) approach can be used for solving the optimization problem therein (Section 1.5.1).

Standard and Bayesian Neural Networks
A Bayesian Neural Network (BNN) is an Artificial Neural Network (ANN) trained with Bayesian Inference [Jospin et al., 2022].In the following, we provide a quick overview of ANNs and their typical estimation based on Backpropagation (Section 1.2.1).We then describe what a Bayesian Neural Network (BNN) is (Section 1.4), provide motivations on why to use a BNN, over a standard ANN (Section 1.3), and lastly introduce Variational Inference (VI) (Section 1.5).

Artificial Neural Networks
For completeness, we review the general ingredients, principles, ideas, and standard terminology behind Artificial Neural Network (ANN).A comprehensive and more detailed introduction to the topic is here out of scope; the interested reader can e.g., consult the accessible book [Haykin, 1998].
Neurons are elementary building blocks which can be thought of as processing units that, when combined, constitute a neural network.Each neuron processes the information presented to its input by applying a transformation to it.When affine neurons are used, the transformation corresponds to computing the weighted sum of the inputs to the neuron (received from the neurons that are connected to it or corresponding to the inputs to the neural network) and generates a value, which is further introduced to a (usually nonlinear) activation function to produce the neuron's output (input to other neurons or the neural network output).In order to account for the need of a shift to the value needed to produce an activation response, a bias is also added as an input to the activation function, which is commonly included in the weighted sum by augmenting the input to the neuron with an additional input with a constant value of 1, associated with the corresponding bias term.
While activation functions squeezing their outputs to a pre-determined range of values, like the sigmoid (with outputs in [0, 1]) or the tanh (with outputs in [−1, 1]) functions, have been widely used in the past, piece-wise linear functions, like the Rectified Linear Unit (ReLU) or the parametric ReLU functions [He et al., 2015], are nowadays widely adopted in building the hidden layers of neural networks.Linear and softmax activation functions are commonly used in the output layer for regression and classification problems, respectively.A common characteristic of activation functions used in neural networks is that they are differentiable with respect to their parameters over the range of their inputs.The transformation performed by an affine neuron is illustrated in Figure 1.
Whenever the information flow between neurons has no feedback (i.e., neurons do not process information resulting from their outputs), in the sense that information flows from the input through the neurons producing the output of the network, the network is referred to as feedforward.Neurons are arranged in layers, and a network formed by neurons in one layer is called single layer network.
When more than one layer forms a neural network, layers are generally called hidden layers since they stand between the input and the output, i.e., the "tangible" information, which consists of the •j , computing the weighted sum a l j .The socalled activation function g(•) is applied to a l j leading to the output o l j , which is sent to nodes at layer l + 1.
T is sequentially parsed to the output, from left to right, following the connections represented in grey which correspond to the weights of the network's layers.input samples and their classification targets/outputs.A feedforward neural network receiving as input a d-dimensional vector and producing a 3-dimensional output is shown in Figure 2.
The most relevant feature of a neural network is its capacity of learning.This corresponds to the ability to improve its outputs (performance in classification) by tuning the parameters (weights and biases) of its neurons.Learning algorithms of neural networks use a set of training data to iteratively update the parameters of a neural network such that some error measure is decreased or some performance measure is increased [see, e.g., Goodfellow et al., 2016].The data D consists of vectors D i = {y i , x i }, with x i representing an input and y i the corresponding target for i = 1, ..., N .Let ŷi denote the output of the network corresponding to the sample x i , that is ŷi = NN θ (x i ), with NN θ (x i ) denoting a Neural Network parametrized over θ and evaluated at x i .An error function E(D, θ) is defined at a particular parameter θ, which is used to guide the learning process.Several error functions have been used to this end, the most widely adopted ones being the mean-squared error (suitable for regression problems) and the cross-entropy (suitable for classification problems).The gradient of the error between the network's outputs ŷi and the targets y i over the entire data set (full-batch) or a subset of the data (mini-batch) is commonly used to update the network parameter values through an iterative optimization process, commonly a variant of the Backpropagation algorithm [Rumelhart et al., 1986].Widely used iterative optimization methods are the Stochastic Gradient Descent (SGD) [Robbins and Monro, 1951], Root Mean Squared Propagation (RMSProp) [Tieleman and Hinton, 2012] and Adaptive Moment Estimation (ADAM) [Kingma and Ba, 2014].
While feedforward neural networks with affine neurons have been briefly described above, a large variety of neural networks have been proposed and used for modeling different input-output data relationships.Such networks follow the main principles as those described above (i.e., they are formed by layers of neurons, which perform transformations followed by differentiable activation functions), but they are realized by using different types of neurons and/or transformations.Examples include the Radial Basis Function (RBF) networks [Broomhead and Lowe, 1988], which replace affine transformations with distance-based transformations, Convolutional Neural Networks [Homma et al., 1987], which receive a tensor input and use neurons performing convolution, Recurrent Neural Networks (e.g., Long-Short Term Memory (LSTM) [Hochreiter and Schmidhuber, 1997] and Gated Recurrent Unit (GRU) [Cho et al., 2014] networks), which model sequences of their inputs by using recurrent units, and specialized types of neural networks, such as the Temporal-Augmented Bilinear Layer (TABL) network [Tran et al., 2019] based on bilinear mapping, and the Neural Bag-of-Features network [Passalis et al., 2020], extending the classical Bag-of-Features model with a differentiable processing suitable to be used in combination with other types of neural network layers.

Motivation for adopting Bayesian Neural Networks
Bayesian Neural Networks are interesting tools under three perspectives: (i) theoretical, (ii) methodological, and (iii) practical.In the following, we shall briefly discuss what we mean by the above three interconnected perspectives.
From a theoretical perspective, BNNs allow for differentiating and quantifying two different sources of uncertainty, namely epistemic uncertainty, and aleatoric uncertainty [see, e.g.Der Kiureghian and Ditlevsen, 2009, from a ML perspective].Epistemic uncertainty is the one referring to the lack of knowledge, and it is captured by p(θ|D).In light of the Bayes theorem, epistemic uncertainty can be reduced with the use of additional data so that the lack of knowledge is addressed as more data are collected.After the data is collected, this results in the update of the prior belief (before the experiment is conducted) to the posterior.Thus, the Bayesian perspective allows the mixing of expert knowledge with experimental evidence.This is quite relevant in small-sample applications where the amount of collected data is inappropriate for classical statistical tools and results to apply (e.g., inference based on asymptotic theory), yet it nevertheless allows the update of the a-priori belief on the parameters, p(θ), into the posterior.On the other hand, the likelihood term captures the aleatoric uncertainty, that is the intrinsic uncertainty naturally embedded in the data, i.e., p(y|θ), in the Bayesian framework is clearly distinguished and separated from the aleatoric one.
Methodologically, is remarkable the ability of Bayesian methods to learn from small data and eventually converge to, e.g., non-Bayesian maximum likelihood estimates or, more generally, to agree with alternative frequentist methods.When the amount of the collected data overwhelms the role of the prior in the likelihood-prior mixture, Bayesian methods can be clearly seen as generalizations of standard non-Bayesian approaches.Within the Bayesian methods family, certain research areas such as PAC-Bayes [Alquier, 2021], Empirical Bayes [Casella, 1985] and Approximate Bayes Computations [Csilléry et al., 2010] deal with such connections very tightly.In this regard, there are many examples in the statistics literature; we focus on the ML perspective.For instance, regularization, ensemble, meta-learning, Monte-Carlo dropout, etc., can all be understood as Bayesian methods, and, e.g., Variational Bayes can be seen as standard linear regression [Salimans and Knowles, 2013].More in general, many ML methods can be seen as approximate Bayesian methods, whose approximate nature makes them simpler and of practical use.Furthermore, as the learned posterior can be reused and re-updated once new data become available, Bayesian learning methods are well-suited for online learning [Opper and Winther, 1999].In this regard, also the explicit use of the prior in Bayesian formulations is aligned with the No-Free-Lunch Theorem [Wolpert, 1996] whose philosophical interpretation, among the others, is that any supervised algorithm implicitly embeds and encodes some form of prior, establishing a tight connection with Bayesian theory [Serafino, 2013, Guedj andPujol, 2021].
From a practical perspective, the Bayesian approach implicitly allows for dealing with uncertainties, both in the estimated parameters and in the predictions.For a practitioner, this is by far the most relevant aspect in shifting from a standard ANN approach to BNNs.Thus, with little surprise, Bayesian methods have been well-received in high-risk application domains where quantifying uncertainties is of high importance.Examples can be found across different fields, such as industrial applications [Vehtari and Lampinen, 1999], medical applications [e.g.Chakraborty and Ghosh, 2012, Kwon et al., 2020, Lisboa et al., 2003], finance [e.g.Jang and Lee, 2017, Sariev and Germano, 2020, Magris et al., 2022a,b], fraud detection [e.g.Viaene et al., 2005], engineering [e.g.Cai et al., 2018, Du et al., 2020, Goh et al., 2005], and genetics [e.g.Ma and Wang, 1999, Liang and Kelemen, 2004, Waldmann, 2018].
As widely recognized, the estimation of BNN is not a simple task due to the generally non-conjugacy between the prior and the likelihood and the non-trivial computation of the integral involved in the marginal likelihood.For this reason, application of BNNs is relatively infrequent, and their use is not widespread across the different domains.As of now, applying Bayesian principles in a plug-and-play fashion is challenging for the general practitioner.On top of that, several estimation approaches have been developed, and navigating through them can indeed be confusing.In this survey, we collect and present parameter estimation and inference methods for Bayesian DL at an accessible level to promote the use of the Bayesian framework.

Bayesian Neural Networks
From the description in Section 1.2.1, it can be seen that the goal of approximating a function relating the input to the output in classical ANNs is treated under an entirely deterministic perspective.Switching towards a Bayesian perspective in mathematical terms is rather straightforward.In place of estimating the parameter vector θ, BNNs target the estimation of the posterior distribution p(θ|D x , D y ), that is [Jospin et al., 2022]: which stands as a simple application of the Bayes theorem.Here we assume, as it is usually the case, that the data D is composed of an input set D x and the corresponding set of outputs D y .In general, D x is a matrix of regressors, and D y is either the vector or matrix (depending on whether the nature of the output is univariate or not) of the variables that the networks aim at modeling based on D x .Alternatively but analogously, D can be thought as the collection of all input-output pairs , where N denotes the sample size, and x i and y i are the input and output vectors of observations for the i-th sample, respectively.Using this notation, While eq. ( 2) provides a theoretical prescription for obtaining the posterior distribution, in practice solving for the form of the posterior distribution and retrieving its parameters is a very challenging task.The estimation of a BNN with MC techniques and VI is discussed in the remainder of the review, here we continue the discussion towards different aspects.
Equation (2) involves all the ingredients required for performing Bayesian inference on ML models, and specifically neural networks.In the first place, eq.(2) involves a likelihood function for the data D y conditional on the observed sample D x and the parameter vector θ.The forward pass parses the input into predictions via some parameter values, such outputs (conditional on the data and the parameters) follow a prescribed likelihood function.Intuitive examples are the Gaussian likelihood (for regression) and the Binomial one (for classification).An underlying neural network is implicit in the likelihood term p(D y |D x , θ), that links the inputs to the outputs.In other words, as is the case for ANNs, the first step in designing a BNN is that of identifying a suitable neural network architecture (e.g., how many layers and of which kind and size) followed by a reasonable assumption for the likelihood function.
A major difference between ANNs and BNNs is that the latter requires the introduction of the prior distribution over the model parameters.After all, a prior must be in place for Bayesian inference to be performed; thus, priors are required in the BNN setup [Jospin et al., 2022].This means that the practitioner needs to decide on the parametric form of the prior over the parameters.
Example 1.Consider a BNN to model the variables D y = {y i } N i=1 where y i ∈ {0, 1}, based on the matrix of covariates D x .The likelihood is of a certain form and parametrized over a neural network whose weights are denoted by θ, i.e., NN θ (•).
We can approach the above problem as a 2-class classification with y i ∈ [0, 1], and derive the likelihood from the Bernoulli distribution where pi = NN θ (x i ) denotes the output of the network for the i-th sample, that is the probability that sample i belongs to class 1.The prior (on the network parameters) can be a diagonal Gaussian p(θ) = N (θ|0, τ I), where τ > 0 is a scalar and I the identity matrix.We can also approach the above problem as a regression to y i ∈ R d and derive the likelihood from the Multivariate Normal distribution where ŷi = NN θ (x i ).Assuming that the covariance matrix Σ −1 is known, the prior on θ could be as well a diagonal Gaussian.If Σ is unknown, the prior could be the product of the above Gaussian prior with, e.g., an Inverse Wishart prior distribution on Σ.In this case, the goal of the Bayesian inference is that of estimating the joint posterior of (θ, Σ).
The inference goal is the posterior distribution.(i) If the problem has a form for which the posterior can be solved analytically, we find p(θ|D x D y ) to be of a known parametric form and identify the parameters characterizing it (standard Bayesian setting, so-called conjugacy between the prior and the likelihood, [e.g., Gelman et al., 1995]).(ii) In general, we may proceed via MC sampling, in which case the estimation leads to a sample, of arbitrary size, approximating the true posterior.
The true posterior remains unknown in its exact form, yet MC enables sampling from it and thus estimating an approximate representation [e.g., Gamerman and Lopes, 2006], see Section 2. (iii) Alternatively, following VI, one sets a certain chosen parametric form for the posterior and optimizes its parameters for a certain objective function [e.g., Nakajima et al., 2019], see Section 1.5.While the actual posterior remains unknown, in VI one seeks an approximation that is optimal in some sense of optimization of a certain objective on the provided data.
Figure 3 provides an analogous representation of Figure 2, now for a BNN.Opposed to traditional ANNs, weights in BNNs are stochastic and represented with distributions.A probability distribution over the weights is learned by updating the prior with the evidence supported by the data.Even though Figure 3 might give the opposite impression, the posterior over the weights is, in general, a truly multivariate distribution where independence among its dimensions generally does not hold.
While the above clarifies that the estimation goal is a distribution whose, e.g., variance can be indicative of the level of confidence in the estimated parameters, the uncertainty associated with the outputs and the generation of the model outputs themselves remains unaddressed.The predictive distribution is defined as [e.g., Gelman et al., 1995] As the posterior (eq.( 2)) is solved, the predictive distribution can also be recovered.Yet, in practice, it is indirectly sampled.Indeed, an intuitive MC-related approach for approximating the predictive distribution is that of sampling N s values from the posterior to create N s realizations of the neural network, each based on a different parameter sample, which are used to provide predictions.This results in a collection of predictions that approximate the actual predictive distribution.In this way, it is relatively simple to recover (approximations of) the predictive distribution from which, e.g., confidence intervals can be constructed.A way to reduce the sample forecast to single values conveying relevant information is by, e.g., using common (sampling) moment estimators [e.g., Casella and Berger, 2021, Ch. 7.2.1].One may evaluate to approximate the posterior mean through model averaging (across the different realizations θ j , j = 1, . . ., N s and thus different outputs) or compute with ε j,i = NN θj (x i ) − ŷi , (8) to approximate the covariance matrix, which is indicative of the uncertainty associated with the prediction.N s corresponds to the number of samples generated from the posterior and used to generate the prediction of the network NN θj (•) receiving as input x i .In classification, one may analogously approximate predictive densities for the joint probability of the different classes and average such probabilities to summarize the average probabilities of each class and implicitly the uncertainties associated with a certain class decision, which is typically determined by the predicted class of maximum probability [e.g., Osawa et al., 2019, Magris et al., 2022a]: with C being the total number of classes and pi,c the predicted probability of class c for the sample i.

Variational Inference (VI)
Let D denote the data and p(D|θ) the likelihood of the data based on a postulated model with θ ∈ Θ a k-dimensional vector of model parameters.Let p(θ) be the prior distribution on θ.The goal of Bayesian inference is the posterior distribution where p(D) is referred to as evidence or marginal likelihood, since p(D) = Θ (D|θ)p(θ)dθ.p(θ) acts as a normalization constant for retrieving an actual probability distribution for p(θ|D).Bayesian inference is generally difficult due to the fact that the evidence is often intractable and of unknown form.In high-dimensional applications, Monte-Carlo methods for sampling the posterior turn challenging and infeasible, and Variational Inference (VI) is an attractive alternative.
VI consists in an approximate method where the posterior distribution is approximated by the socalled variational distribution [e.g., Blei et al., 2017, Nakajima et al., 2019, Tran et al., 2021b].
The variational distribution is a probability density q(θ), belonging to some tractable class of distributions Q such as, e.g., the Exponential family.VI thus turns the Bayesian inference problem in eq. ( 10) into that of finding the best approximation q (θ) ∈ Q to p(θ|D) by minimizing the Kullback-Leibler (KL) divergence from q(θ) to p(θ|D) [Kullback and Leibler, 1951], By simple manipulations of the KL divergence definition, it can be shown that Since log p(D) is a constant not depending on the model parameters, the KL minimization problem is equivalent to the maximization problem of the so-called Lower Bound (LB) on log p(D) [e.g., Nakajima et al., 2019], For any random vector θ and a function g(θ) we denote by E f [g(θ)] the expectation of g(θ) where θ follows a probability distribution with density To make explicit the dependence of the LB on some vector of parameters ζ parametrizing the variational posterior we write We operate within the Fixed-Form Variational Inference (FFVI) framework, where the parametric form of the variational posterior is set [e.g., Tran et al., 2021b].I.e., FFVI seeks at finding the best q ≡ q ζ in the class Q of distributions indexed by a vector parameter ζ that minimizes the LB L(ζ).In this context, ζ is called variational parameter.A common choice for Q is the Exponential family, and ζ is the corresponding natural parameter.

Estimation with Stochastic Gradient Descent (SGD)
A straightforward approach to maximize L(ζ) is that of using a gradient-based method such as Stochastic Gradient Descent (SGD), Adaptive Moment Estimation (ADAM) [Kingma and Ba, 2014], or Root Mean Squared Propagation (RMSProp) [Tieleman and Hinton, 2012].The form of the basic SGD update is where t denotes the iteration, β t a (possibly adaptive) step size, and ∇ζ Under a pure Gaussian variational assumption, it is instinctive to optimize the LB for the mean vector ζ 1 = µ and variance-covariance matrix ζ 2 = Σ.In the wider FFVI setting with Q being the Exponential family, the LB is often optimized in terms of the natural parameter λ [Wainwright and Jordan, 2008].The application of the SGD update based on the standard gradient is problematic because it ignores the information geometry of the distribution q ζ [Amari, 1998], as it implicitly relies on the Euclidean distance to capture the dissimilarity between two distributions in terms of the euclidean norm ||ζ t − ζ|| 2 , which can be a quite poor and misleading measure of dissimilarity [Khan and Nielsen, 2018].By replacing the Euclidean norm with the KL divergence, the SGD update results in the following natural gradient update: The natural gradient update results in better step directions toward the optimum when optimizing the distribution parameter.The natural gradient of L(λ) is obtained by rescaling the gradient ∇ λ L(λ) by the inverse of the Fisher Information Matrix (FIM) where subscript in I −1 λ remarks that the FIM is expressed in terms of the natural parameter λ.By replacing in the above ∇ λ L(λ) with a stochastic estimate ∇λ L(λ) one obtains a stochastic natural gradient update.Example 2 (Problem statement for Variational Inference (VI)).Consider a BNN to model the targets y i , based on the covariates x i .The likelihood, of a certain form, is parametrized over a neural network, whose weights are denoted by θ.The prior could be a Gaussian distribution with, e.g., zero-mean, diagonal p(θ) = N (θ|0, I/τ ) or not p(θ) = N (θ|0, Σ 0 ).Q is the set of multivariate Gaussian distributions, specified, e.g., in terms of the natural parameter λ.
The objective is that of finding the corresponding variational parameter such that the LB E q λ [log p(θ) − log q λ (θ) + p(D|θ)] is maximized.The update of the variational parameter λ follows a gradient-based method with natural gradients.The training terminates after the LB L(λ) does not improve for a certain number of iterations: the terminal λ provides the natural parameter of the variational posterior approximation minimizing the KL divergence to the true posterior p(θ|D).

Monte Carlo Markov Chain (MCMC)
Monte Carlo Markov Chain (MCMC) is a set of methods for sampling from a probability distribution.MCMCs have numerous applications, and especially in Bayesian statistics are a fundamental tool.The foundation of MCMC methods are Markov Chains, stochastic models describing a sequence of events in which the probability of each event depends only on the state of the previous one [Gagniuc, 2017].By constructing a Markov Chain that has the desired distribution as its stationary distribution, towards which the sequence eventually converges, one can obtain samples from it, i.e., one can sample any generic probability distribution, including, e.g., a complex, perhaps multi-modal, Bayesian posterior.Early samples may be autocorrelated and not representative of the target distribution, so that MCMC methods generally require a burnout period before attaining the so-called stationary distribution.In fact, while the construction of a Markov Chain converging to the desired distribution is relatively simple, determining the number of steps to achieve such convergence with an acceptable error is much more challenging and strongly dependent on the initial setup and starting values.With burnout, the large collection of samples is practically subsampled by discarding an initial fraction of draws (e.g., 20%) to obtain a collection of approximately independent samples from the desired distribution.An accessible introduction to Markov Chains can be found in [Gagniuc, 2017], for a dedicated monograph on MCMC methods oriented toward Bayesian statistics and applications see, e.g., [Gamerman and Lopes, 2006].
Within the class of MCMC methods, some popular ones are not effective in large Bayesian problems such as BNNs.For example, the plain Gibbs sampler [Geman and Geman, 1984], despite its simplicity and desirable properties [Casella and George, 1992], suffers from residual autocorrelation between successive samples and becomes increasingly difficult as the dimensionality increases in multivariate distributions [e.g.Lynch, 2007, Ch. 4].We review the most widespread MC approaches in the context of performing Bayesian learning for Neural networks.

Metropolis-Hastings (MH)
The Metropolis-Hastings (MH) algorithm [Metropolis et al., 1953, Hastings, 1970] is particularly helpful in Bayesian inference as it allows drawing samples from any probability distribution p, given that a function f proportional up to a constant to p can be computed.This is particularly convenient as it allows to sample a Bayesian posterior by only evaluating f (θ) = p(θ|y)p(θ), completely excluding the normalization factor from the computations.The values of the Markov Chain are sampled iteratively, with each value depending solely on the preceding one: at each iteration, based on the current value, the algorithm picks a candidate value θ (proposed value), which is either accepted or rejected randomly with a probability that depends on the current and earlier values.Upon acceptance, the proposed value is used for the next iteration, otherwise is discarded, and the current value is used in the next iteration.As the algorithm proceeds and more sample values are generated, the sampled-value distribution more and more closely approximates the target distribution p.
A key ingredient in MH is the proposal density determining the drawing of the proposed value at each iteration.This is formalized by an arbitrary probability density g(θ |•), upon which depends the probability of drawing θ given the previous value θ. g is usually assumed symmetric, and a common choice is provided by a Gaussian distribution centered on θ.Algorithm 1 summarizes the above steps.

Algorithm 1 Metropolis-Hastings
Initialize The acceptance ratio α is representative of the likelihood of the proposed sample θ over the current one θ t according to p. Indeed, α = f (θ )/f (θ t ) = p(θ )/p(θ t ) as f is proportional to p.A proposed sample value θ that is more probable than θ t (α > 1) is always accepted; otherwise, it may be rejected with probability α.The algorithm thus moves around the sample space, tending to stay in regions where p is of high density and, seldomly, in regions of low density.The final collection of samples follows the distribution p.As the Markov chain eventually converges to the target distribution p, initial samples may be quite incompatible with p, especially if the algorithm is initialized at a low-density region.Thus, it is customary to discard a number B of samples and retain only the subsample {θ t } N t=B .Note that by construction, successive samples of the Markov chain are correlated.Even though the chain eventually converges to p nearby samples are correlated, causing a reduction of the effective sample size (e.g., for E θ∼p θ [θ] the central limit theorem applies but, e.g., the limiting variance is inflated by the non-zero autocorrelation in the chain).
An important feature of the MH algorithm is that it is applicable to high dimensions as it does not suffer from the course of dimensionality problem, causing an increasing rejection rate as the number of dimensions increases.This makes MH suitable for large Bayesian inference problems such as training BNNs.

Hamiltonian Monte Carlo (HMC)
Hamiltonian Monte Carlo (HMC) generates efficient transitions by using the derivatives of the density function being sampled by using approximate Hamiltonian dynamics, later corrected for performing an MH-like acceptance step [Neal et al., 2011].
HMC augments the target probability density p(θ) by introducing an auxiliary momentum variable ρ and generating draws from Typically the auxiliary density is taken as a multivariate Gaussian distribution, independent of θ: Σ can be conveniently set to the identity matrix, restricted to a diagonal matrix, or estimated from warm-up draws.The Hamiltonian is defined upon the joint density p(ρ, θ): The term T (ρ|θ) = − log p(ρ|θ) is usually called kinetic energy and V (θ) = − log p(θ) is called potential energy.To generate transitions to a new state, first, a value for the momentum is drawn independently from the current θ; then, Hamilton's equations are adopted to describe the evolution of the joint system (ρ, θ), i.e.: By having the momentum density being independent of the target density, p(ρ|θ) = p(ρ), ∂T /∂θ = 0, the transitions are governed by the derivatives Note that −∂V /∂θ is simply the gradient of the negative loglikelihood, which can be computed using automatic differentiation.The main difficulty is the simulation of the Hamiltonian dynamics, for which there is a variety of approaches [see, e.g.Leimkuhler and Reich, 2005, Berry et al., 2015, Hoffman et al., 2014].Yet, to solve the above system of differential equations, a leapfrog integrator is generally used due to its simplicity and volume-preservation and reversibility properties [Neal et al., 2011].The leapfrog integrator is a numerically stable integration algorithm specific to Hamiltonian systems.It discretizes time using a step size ε and alternates half-step momentum updates and full-step parameter updates: By repeating the above steps L times, a total of Lε time is simulated, and the resulting state is (ρ , θ ).Note that both L and ε are hyperparameters, and their tuning is often difficult in practice.
In this regard, see the Generalized HMC approach of [Horowitz, 1991] and developments aimed at resolving the tuning of the leapfrog iterator [Fichtner et al., 2020, Hoffman andSountsov, 2022], Instead of generating a random momentum vector right away and sampling a new state (ρ , θ ), to account for numerical errors in the leapfrog integrator (an analysis in this regard is found in [Leimkuhler and Reich, 2005]), a Metropolis-hasting step is used.The probability of accepting the proposal (ρ , θ ) by transitioning from (ρ, θ) is If the proposal (ρ , θ ) is accepted, the leapfrog integrator is initialized with a new momentum draw and θ ; otherwise, the same (ρ, θ) parameters are returned to start the next iteration.The HMC procedure is summarized in Algorithm 2.Besides the difficulty of calibrating the hyperparameters L and ε, HMC suffers from multimodality, yet the Hamiltonian boosts the local exploration for unimodal targets.
Algorithm 2 Hamiltonian Monte Carlo Initialize θ 0 , ρ 0 , set ε for t = 0 to N max do Simulate ρ ∼ N (0, Σ) Monte Carlo Dropout (MCD) is an indirect method for Bayesian inference.Dropout has been earlier proposed as a regularization method for avoiding overfitting and improving neural networks' predictive performance [Srivastava et al., 2014].This is achieved by applying a multiplicative Bernoulli noise on the neurons constituting the layers of the network.This corresponds to randomly switching off some neurons at each training step.The dropout rate sets the probability p i of a neuron i being switched off.Though Bernoulli noise is the most common choice, note that other types of noise can be as well adopted [e.g.Shen et al., 2018].Neurons are randomly switched off only in the training phase, and the very same network configuration in terms of the activated and disabled neurons is used during backpropagation for computing gradients for weights' calibration.On the other hand, all the neurons are left activated for predictions.Though it is intuitive that the above procedure implicitly connects to model averaging across different randomly pruned architectures obtained from a certain DL network, the exact connection between MC dropout and Bayesian inference follows a quite elaborated theory.[Gal and Ghahramani, 2016] shows that a neural network of arbitrary depth and non-linearity with dropout applied before every single layer is mathematically equivalent to an approximation to the probabilistic deep Gaussian Process (GP) model [Damianou and Lawrence, 2013], and [Jakkala, 2021] for a recent survey.That is, the dropout objective minimizes the KL divergence between a certain approximate variational model and the deep GP.A treatment limited to multi-layer perceptron networks is provided in [Gal and Ghahramani, 2015].
With ŷ being the output of a Neural Network with L layers whose loss function is E, for each layer i = 1, . . ., L let W i denote the corresponding weight matrix of dimension K i × K i−1 , and b i the bias vector of dimension K i .Be y n the target for the input x n for n = 1, . . ., N and denote the input and output sets respectively with D x and D y .A typical optimization objective includes a regularization term weighted by some decay parameter λ, that is Now consider a deep Gaussian process for modeling distributions over functions corresponding to different network architectures.Assume its covariance is of the form where σ(•) is an element-wise non-linearity, and p(w), p(b) distributions.Now let W i be a random matrix of size . The predictive distribution of the deep GP model can be expressed as where p(ω||D x , D y ) is the posterior distribution.p(y n |x n , ω) is determined by the likelihood, while ŷn is a function of x n and ω: m i are vectors of size K i for each GP layer.For the intractable posterior p(ω|D x , D y ), [Gal and Ghahramani, 2016] uses the variational approximation q(ω) defined as where the collection of probabilities p i and matrices M i , i = 1, . . ., L constitute the variational parameter.Thus, q stands as a distribution over (non-random) matrices whose columns are randomly set to zero, and z i,j = 0 implies that the unit j in layer i − 1 is dropped as an input to layer i.For minimizing the KL divergence form q to p(ω|D x , D y ), the objective corresponds to By use of Monte Carlo integration and some further approximations [see Gal and Ghahramani, 2016, for details], the objective reads which, up to the constant 1 τ N , is a feasible and unbiased MC estimator of eq. ( 37) where ω denotes a single MC draw from the posterior ωn ∼ q(ω).By taking E(y n , ŷn ) = − log p(y n |x n , ωn )/τ eq. ( 38) and eq.( 29) are equivalent for an appropriate choice of the hyperparameters τ and l.This shows that the minimization of the loss in eq. ( 29) with dropout is equivalent to minimizing the KL divergence from q to p(ω|D x , D y ), thus performing VI on the deep Gaussian process.
With an SGD approach, one can maximize the above LB and estimate the variational parameters from which one can simply obtain samples from the predictive distribution q(y |x ), and approximate its mean by the naive MC estimator: x denotes a new observation, not in D x , for which the corresponding prediction is ŷ .I.e., the predictive mean is obtained by performing N s forward passes through the network with Bernoulli realizations {z s 1 , . . ., . Such average predictions are generally referred to as MC dropout estimates.Similarly, by simple moment-matching, one can estimate the predictive variance and higher-order statistics synthesizing the properties of q(y |x ).
The predictive distribution is, in general, a multi-modal distribution resulting from superposing bimodal distributions on each weight matrix column.This constitutes a drawback of MCD, as well the implicit VI on a GP.Furthermore, the VI approximation in eqs.( 34)-( 36) may be adequate or not.It is clear that even though MCD is a possibility for VI in deep-learning models, it is constrained by the very specific form in eq. ( 34) of the variational posterior that implicitly corresponds to performing VI on a deep GP.Furthermore, there is evidence that MCD does not fully capture uncertainty associated with model predictions [Chan et al., 2020], and there are issues related to the use of improper priors and singularity of the approximate posterior.The latter ones are addressed and explored in [Hron et al., 2018], suggesting the use of the so-called Quasi-KL divergence as a remedy.Clearly, high dropout rates drive the convergence rate slow, expand the network training time, and can cause important training data to be missed or given little relative importance.However, compared to the traditional approach for neural networks, applying dropout places no additional effort and is often of faster training than other VI methods.Furthermore, if a network has been trained with dropout, only by including an additional form of regularization acting as a prior turns the ANN into a BNN, without requiring re-estimation [Jospin et al., 2022].

Bayes-By-Backprop (BBB)
A common approach for estimating the variational posterior over the networks' weights is the Bayes-By-Backprop (BBB) method of Blundell et al. [2015], perhaps a breakthrough in probabilistic deeplearning as a practical solution for Bayesian inference.
The key argument in [Blundell et al., 2015] is the use of the local reparametrization trick under which the derivative of an expectation can be expressed as the expectation of a derivative.It introduces a random variable ε having a probability density given by q(ε) and a deterministic transform t(θ, ε) such that w = t(θ, ε).The main idea is that the random variable ε is a source of noise that does not depend on the variational distribution, and the weights w are sampled indirectly as a deterministic transformation of ε, leading to a training algorithm that is analogous to that used in training regular networks.Indeed, by writing w as w = t(θ, ε), in place of evaluating which can be complex and rather tedious, under the assumption q(ε)dε = q(w|θ)dw, Blundell et al. [2015] prove that With f (w, θ) = log q(w|θ) − log p(w)p(y|w), the right side of eq. ( 41) and eq.( 42) provide an alternative approach for the estimation of the gradients of the cost function with respect to the model parameters.
The sampled value ε ∼ q(ε), resampled at each iteration, is independent of the variational parameters, while w is not directly sampled but here it is a deterministic function of ε.Given ε, all the quantities in the square bracket of eq. ( 42) are non-stochastic, enabling the use of backpropagation.
A single draw for ε approximates the right side of eq. ( 41), and suffices for providing an unbiased stochastic gradient estimation of the relevant gradient on the left side.Equation ( 41) makes explicit the possibility of using automatic differentiation to compute the gradient of f with respect to the parameter θ.By using a single sampled draw θ for approximating the expectation on the right side of eq. ( 41), the only parameter in the loss is θ, and the use of backpropagation for evaluating the gradients is straightforward.Equation ( 42) instead employs backpropagation in the "usual" sense, involving gradients of the cost with respect to the network parameters w, further rescaled by ∂w/∂θ and shifted by ∂f (w, θ)/∂θ.Equation ( 42) concerns the usual backpropagation computations in terms of the network's weights, the specific form of the partial derivative with respect to θ that the choice of t implies, while the last term depends on the chosen form of the variational posterior only (w is here not seen as a function of θ, as the form of eq. ( 42) results from applying the multi-variable chain rule).This results in a general framework for learning the posterior distribution over the network's weights.The following algorithm 3 summarizes the BBB approach.
Algorithm 3 Bayes-By-Backprop Algorithm 3 is initialized by preliminary setting the initial values of the variational parameter θ and, of course, by specifying the form of the prior and the posterior along with the form of the likelihood involving the outputs of the forward pass obtained from the specified underlying network structure.The update is very similar to the one employed in standard non-Bayesian settings, where standard optimizers such as ADAM are applicable.It is the applicability of standard optimization algorithms and the use of classic backpropagation that constitute the major breakthrough element in BBB, making it a feasible approach for Bayesian learning.
To make the description more explicit and aligned with the following sections, we present the case where the variational posterior is a diagonal Gaussian with mean µ and covariance matrix σ 2 I.In this case, the transform t takes the simple and convenient form As σ is required to be always non-negative, [Blundell et al., 2015] adopts the reparametrization σ = log(1 + exp(ρ)) and the variational posterior parameter θ = (µ, ρ).In this case, Algorithm 4 summarizes the BBB approach.
As for Algorithm 4, one may backpropagate the gradients of f w.r.t.µ and ρ directly.Alternatively, as for Algorithm 3, one may use backpropagation for computing the gradients ∂f (w, θ)/∂w, which are furthermore shared across the updates for µ and ρ, or, if preferred, adopt a general automatic differentiation setup, if, e.g., the form of the variational likelihood does not allow for a simple analytic form of the gradient.

Exponential family and Natural gradients
Assume q λ (θ) belongs to an exponential family distribution.Its probability density function is parametrized as where λ ∈ Ω is the natural parameter, φ(θ) the sufficient statistic.A(λ) = log h(θ) exp(φ(θ) λ)dν is the log-partition function, determined upon the measure ν, φ and the function h.The natural parameter space is defined as When Ω is a non-empty open set, the exponential family is referred to as regular.Furthermore, if there are no linear constraints among the components of λ and φ(θ), the exponential family in eq. ( 44) is said of minimal representation.Non-minimal families can always be reduced to minimal families through a suitable transformation and reparametrization, leading to a unique parameter vector λ associated with each distribution [Wainwright and Jordan, 2008].The mean (or expectation) parameter m ∈ M is defined as a function of λ, m(λ Under minimal representation, A(λ) is convex, thus the mapping ∇ λ A = m : Ω → M is one-to-one, and I λ is positive definite and invertible [Nielsen and Garcia, 2009].M denotes the set of realizable mean parameters.Therefore, under minimal representation we can express λ in terms of m and thus L(λ) in terms of L(m) and vice versa [Khan and Nielsen, 2018].Example 3 (The Gaussian distribution as an exponential-family member).The multivariate Gaussian distribution N (µ, Σ) with k-dimensional mean vector µ and covariance matrix Σ can be seen as a member of the exponential family (eq.( 44)).Its density reads where and constitutes the common parametrization of the multivariate Gaussian distribution in terms of its mean and variance-covariance matrix.

By applying the chain rule, ∇
The quantity ∇λ L is referred to as the natural gradient of L with respect to λ and it is obtained by pre-multiplying the Euclidean gradient by the inverse of the FIM (parametrized in terms of λ).In general, L can be a generic function whose derivative with respect to a parameter λ (not necessarily the natural parameter) exists.The standard reference for natural gradients computation is the seminal work of Amari [1998].Within a SGD context, the application of simple Euclidean gradients is problematic as it ignores the information geometry of the distribution q λ .Euclidean gradients implicitly rely on the Euclidean norm to capture the dissimilarity between two distributions which can be a quite poor dissimilarity measure [Khan and Nielsen, 2018].In fact, the SGD update can be obtained by writing and setting to zero its derivative.Although the above implies that λ moves in the direction of the gradient, it remains close to the previous λ t in terms of Euclidean distance.As λ is a parameter of a distribution, the adoption of the Euclidean measure is misleading.An Exponential family distribution induces a Riemannian manifold with a metric defined by the FIM [Khan and Nielsen, 2018].By replacing the Euclidean metric with the Riemannian one, the resulting update is indeed expressed in terms of the natural parameter: generally referred to as natural gradient update.More in general, one could replace the Euclidean distance with a proximity function such as the Bregman divergence and obtain richer classes of SGD-like updates, like mirror descent (which can be interpreted as natural gradient descent), see, e.g., [Nielsen, 2020].A very interesting point on the limitations of plain gradient search is made in [Wierstra et al., 2014] concerning the impossibility of locating, even in a one-dimensional case, a quadratic optimum.The example provided therein involves the Gaussian distribution, pivotal in VI.
For an one-dimensional Gaussian distribution with mean µ and standard deviation σ, the gradient of L with respect to the parameters µ and σ lead to the following SGD updates: For the updates to converge and the optimum to be precisely located, σ must decrease (i.e., the distribution shrinks around µ).The fact that σ appears in the denominator of both the updates is problematic: as it decreases, the variance of the updates increases as ∆ µ ∝ 1 σ and ∆ σ ∝ 1 σ .The updates become increasingly unstable, and a large overshooting update makes the search start all over again rather than converging.Increased population size and small learning rates cannot avoid the problem.The choice of the starting value is problematic, too: starting with σ >> 1 makes the updates minuscule; conversely, σ << 1 makes them huge and unstable.[Wierstra et al., 2014] discusses how the use of natural gradients fixes this issue that, e.g., may arise with BBVI.
Algorithm 5 summarizes the generic scheme upon the implementation of a natural gradient update.In Algorithm 5, ζ denotes a generic variational parameter, e.g., the natural parameter or not, while methods for evaluating ∇ θ L, I, and efficiently computing its inverse I −1 are discussed in the following sections.
Algorithm 5 Stochastic Gradient Descent with natural gradients repeat Compute: A major issue in VI is that it heavily relies upon model-specific computations, on which a generalized, ready-to-use, and plug-and-play optimizer is difficult to design.Black-box methods aim at providing solutions that can be immediately applied to a wide class of models with little effort.In the first instance, the ubiquitous use of model's gradients that traditional ML and VI approaches rely upon struggles with this principle.As [Ranganath et al., 2014] describes, for a specific class of models, where conditional distributions have a convenient form and a suitable variational family exists, VI optimization can be carried out using closed-form coordinate ascent methods [Ghahramani and Beal, 2000].In general, there is no close-form solution resulting in model-specific algorithms [Jaakkola and Jordan, 1997, Blei and Lafferty, 2007, Braun and McAuliffe, 2010] or generic algorithms that involve model-specific computations [Knowles andMinka, 2011, Paisley et al., 2012].As a consequence model assumptions and model-specific functional forms play a central role, making VI practical.The general idea of Black-box VI is that of rewriting the gradient of the LB objective as the expectation of an easy-to-compute function of the latent and observed variables.The expectation is taken with respect to the variational distribution, and the gradient is estimated by using stochastic samples from it in a MC fashion.Such stochastic gradients are used to update the variational parameters following an SGD optimization approach.Within this framework, the end-user is required to develop functions only for evaluating the model log-likelihood, while the remaining calculations are easily implemented in libraries of general use applicable to several classes of models.Black-box VI falls within stochastic optimization where the optimization objective is the maximization of the LB using noisy, unbiased, estimates of its gradient.As such, variance reduction methods have a major impact on stability and convergence, among them control variates are the most effective and of immediate implementation.

Black-Box Variational Inference (BBVI)
Black-Box Variational Inference (BBVI) optimizes the LB with stochastic optimization, through an unbiased estimator of its gradients obtained from samples from the variational posterior [Ranganath et al., 2014].By using the LB definition and the log-derivative trick on the gradient of the LB with respect to the variational parameter, ∇ ζ L can be expressed as where ζ denotes the parameter of the variational distribution q ζ .The above expression rewrites the gradient as an expectation of a quantity that does not involve the model's gradients but only those of log q(w|ζ).A naive noisy unbiased estimate of the gradient of the LB is immediate to obtain with N s samples obtained from the variational distribution, where θ s ∼ q(θ|ζ).The above MC estimator enables the immediate and feasible computation of the LB gradients as, given a sample θ s , log q(θ s |θ) is a quantity that solely depends on the form of the variational posterior and can be of simple form.On the other hand, log p(D, θ) − log q(θ|θ) is immediate to compute as it only requires evaluating the logarithm of the joint p(D, θ s ) and the density of the variational distribution in θ s .This process is summarized in Algorithm 6.If sensible, one may assume that log p(D, θ) = log p(D|θ)p(θ) but this is not explicitly required as of [Ranganath et al., 2014]: there are no assumptions on the form of the model; the approach only requires the gradient of the variational likelihood with respect to the variational parameters to be feasible to compute.
In [Ranganath et al., 2014], the authors employ an adaptive learning rate satisfying the Robbins Monroe conditions t β t = ∞ and t β 2 t < ∞, and for controlling the variance of the stochastic gradient estimator adopt Rao-blackewllization [Rao, 1945, Blackwell, 1947, Robert and Roberts, 2021] and use the of control variates [e.g.Lemieux, 2014, Robert et al., 1999, Ch. 3] within Algorithm 6.

Natural-Gradient Black-Box Variational Inference (NG-BBVI)
We shall review the approach of Trusheim et al. [2018] boosting BBVI with natural gradients, referred to as Natural-Gradient Black-Box Variational Inference (NG-BBVI).The FIM corresponds to the outer product of the score function with itself (see Section 5) and is furthermore equal to the second derivative of the KL divergence to the approximate posterior q(x|ζ): For the practical implementation, [Trusheim et al., 2018] uses a mean-field restriction on the variational model, i.e. the joint is factorized into the product of K independent terms, where each term is in general a multivariate distribution: The above restriction is also suggested by Ranganath et al. [2014] in order to allow for Rao-Blackwellization [Robert and Roberts, 2021] as a tool to be used in conjunction with control variates [e.g.Lemieux, 2014, Ch. 3] for reducing the variance of the stochastic gradient estimator.Under the above assumption, the FIM simplifies to: which significantly simplifies the general form q ζ (θ) while implicitly enabling Rao-Blackwellization with the variable-wise local expectations and thus reducing the variance of the FIM, estimated via a Monte-Carlo approach.In fact, besides a few variational models it is difficult to compute the above expectations analytically so Trusheim et al. [2018] adopts the following naive MC estimator: ∼ q ζi (θ i ) denoting a sample from the i-th factor of the posterior mean-field approximation.Note that the above does not introduce additional computations as the score of the samples θ (s) i is anyway required in the computation of the LB gradient.Furthermore, instead of using a plain SGD-like update, [Trusheim et al., 2018] adopts an ADAM-like version, boosted with natural gradient computations.Algorithm 7 summarizes the NG-BBVI approach.
The NG-BBVI implementation is slightly more complex than the original BBVI, see Algorithm 7. The MC computation involves both the black-box stochastic gradient estimation and the estimation of the optimal control variate coefficient a .Thus the posterior samples are split into two subsets.The first one X aimed at estimating a , and the second one Y at implementing the MC estimators, independently from X, and with the control variate correction term a earlier computed.The computation of the FIM follows immediately from eq. ( 55), and the computation of ∇ ζ k L is analogous to BBVI.The last four lines of Algorithm 7 correspond to the implementation of the ADAM update, operators are intended to be applied element-wise, β 1 , β 2 (exponential decay rates) are typical ADAM hyperparameters, ε > 0 is a small offset preventing divisions by zero.[Trusheim et al., 2018] differs from BBVI by the use of natural gradients (and the adoption of the ADAM-like update, though applicable to BBVI as well).On the other hand, the use of control variates and Rao-Blackwellization for variance reduction is found in both BBVI and NG-BBVI.As the natural gradient approach is preferable for the reasons discussed in Section 5, NG-BBVI is favored over BBVI.
The use of the black-box framework for computing the gradients of the LB along the MC estimator for the FIM renders NG-BBVI of general applicability and not constrained to a certain form of the variational posterior.Yet the MC-computations of the FIM are implicitly approximate, whereas for certain distributions the FIM computation can be carried out analytically and in an exact form.NG-BBVI furthermore requires the inversion of the FIM, which is a computational bottleneck.The following VON Khan and Nielsen [2018], VADAM [Khan et al., 2018a] and VOGN [Khan et al., 2018a, Osawa et al., 2019] methods indeed fix this issue: assuming a variational posterior within the exponential distribution family, natural gradients are enabled without the direct computation of the FIM and its inverse.
7 Natural gradient methods for exponential-family variational distributions In the following subsections, we review methods based on Natural gradients and Exponential-family variational approximations.The following techniques are built on natural parameter updates in the natural parameter space and rely on simplified but exact FIM computations based on the natural/expectation parameter duality (eq.( 47)).

Exact gradient computations for the exponential family
The computation of the FIM required in the natural gradient computation is, in general, not trivial.In a generic perspective, not bound to a specific variational form, the sampling approach for the FIM estimation of [Trusheim et al., 2018] is feasible.Yet for certain distributions, namely for those in the Exponential family class, natural gradients can be computed in an exact form with an analytical solution which furthermore does not involve the computation of the FIM.
The theoretic foundation of such a viable approach is provided in [Khan and Nielsen, 2018] and traces back to eq. ( 47).For an Exponential family of minimal representation, the natural gradient with respect to the natural parameter λ is equal to the gradient with respect to the expectation parameter m.This is a powerful result that allows the computation of the natural gradient as an Euclidean gradient, avoiding the computation of the FIM and its inversion.
This section presents some baseline methods using the above duality for the natural gradient computation.Differently from BBB, BBVI, and NG-BBVI the following approaches explicitly deal with variational distributions members of the Exponential family with a focus on updating their natural parameter: From the above updates on natural parameters, update rules for alternative and perhaps a more usual parametrization can often be obtained, see e.g.Section 7.2.

Natural-Gradient Variational Inference (NGVI)
Natural-Gradient Variational Inference (NGVI) constitutes a baseline methodology for natural gradient computation under a Gaussian variational distribution, upon which several other approaches have been developed.
The gradients in the natural parameter space can be computed under the expectation parametrization as Euclidean gradients.[Khan and Nielsen, 2018] shows that such gradients are of simple form and correspond to By using the definition of natural gradients in terms of µ and Σ, the update in eq. ( 60) for the natural parameter The above two constitute the NGVI update rules for updating the mean µ and covariance matrix Σ of the variational posterior with a natural gradient update, that however does not involve the computation of the FIM as it relies on Euclidean gradients.
For a diagonal covariance matrix Σ = diag σ 2 , the corresponding NGVI updates read With respect to the NGVI update two points are important to stress out.First, at each iteration, the update for µ implicitly requires Σ t+1 .This means that the update for µ follows that for Σ −1 and that µ readily uses the one-step ahead updated information on Σ.Though it may appear counterintuitive, [Lyu andTsang, 2021, Magris et al., 2022b] show that this update is not optimal (in the terms therein discussed), while an update of the form µ t+1 = µ t + βΣ t [∇ µ L] would be.Also, note that the update for µ involves Σ t+1 and not Σ −1 t+1 , meaning that in the NGVI an online inversion of Σ −1 is implicitly required at each iteration.Clearly, for the diagonal case, this is trivial and effortless to obtain.Second, in the full-covariance case, there is no guarantee that the updates guarantee Σ to be a positive-definite covariance matrix.This issue is tackled in Section 8.For the diagonal case, the constraint on Σ results in guaranteeing the positivity of the entries in the diagonal.This can be achieved via a proper reparametrization, e.g.BBVI updates ρ where σ = log(1 + exp(ρ)), or [e.g.Tan, 2021] updates the Cholesky factor.Alternatively, the learning rate can be adapted to guarantee that the step size does not drive the updates σ −2 negative [e.g., Khan andNielsen, 2018, Magris et al., 2022c].

Variational Online Newton (VON)
A computational burden in NGVI is that the gradients of the LB are still required: Variational Online Newton (VON) develops on NGVI but does not require the gradients of the variational objective.Furthermore, it only involves the gradient and Hessian of the model log-likelihood which can be computed with usual backpropagation.Khan et al. [2018b] express the lower bound as where N is the sample size, and f (θ) = − 1 N N i=1 log p(D i |θ) is negative log-likelihood of the model, i.e. standard MLE objective, where D i denotes a data example, i.e.D i = (y i , x i ).VON uses the theoretical results of [Opper andArchambeau, 2009, Rezende et al., 2014] to express the gradients of the LB objective in terms of gradient and Hessian of f (θ).By linearity of the expectation, the gradients of the L consist of the sum of the gradients of three expectation terms, in particular: (71) where g = ∇ θ f (θ) and H(θ) = ∇ 2 θθ f (θ) denote the gradient and Hessian of the MLE objective, respectively.With these relations the gradients of the LB objective write and By using these gradients in the NGVI update, one obtains where the expectations can be again evaluated via MC sampling.By using a single draw θ t ∼ N (θ|µ t , Σ t ), the feasible update reads To obtain a form for the update that resembles Newton's method where the scaling matrix is estimated online, [Khan et al., 2018b] defines S t = Σ −1 t − λI /N and conversely Σ t = (N (S t + λI/N )) −1 , and write the final form of the VON update Similarly, for a diagonal covariance matrix (thus under a mean-field assumption), with Algorithm 8 Variational Online Newton where the division is intended to be element-wise.Algorithm 8 summarizes the main elements of VON implementation.
In a mini-batch setting for estimating the stochastic gradient, with M denoting a mini-batch containing M samples, the stochastic estimates enable the practical implementation of the VON update by replacing g and H.To make this statement clear, think of f (θ) as the typical negative log-likelihood of a sample (as it is an average across samples), then g is the typical gradient for a sample in θ and, analogously, H is interpreted as the typical (average) value of the Hessian evaluated in θ, resulting when using a single data point.Stochastic gradient estimation estimates g by using a single observation D i picked at random as an unbiased estimate of the actual gradient of f , which would require the parsing of the entire sample.Analogously, one constructs a stochastic estimate of the Hessian with one or M observations (the higher M the lower the variance of the estimator, which is in any case unbiased).

Variational ADAM (VADAM)
The principle of Variational ADAM (VADAM) is that of augmenting the natural gradient update by incorporating a momentum factor, i.e., which slightly extends the form of the update in eq. ( 59).
Under a Gaussian variational q, [Khan et al., 2018b] expresses the momentum update as a VON update with momentum and recovers a variational version of an RMSProp update, to obtain the following updates where ) and β,γ 1 are learning rates.Note that in the above updates the Hessian is estimated as a squared gradient: details are provided in Section 7.5.These updates can be implemented and used in their actual form, yet they correspond to an ADAMlike update.Indeed the above update has the same form of an adaptive version of Polyak's heavy ball method.[Wilson et al., 2017] establishes a relation between the form of eq. ( 87) and the ADAM update, and in particular that the ADAM update can be written as an adaptive version of the Polyak's heavy ball method.Upon introducing the typical bias correction terms of ADAM, [Khan et al., 2018b] expresses eq. ( 87) as an ADAM update.With respect to a true ADAM update, the model weights are stochastically sampled from the posterior, resulting in a Variational version of ADAM (VADAM).For the full derivation, which is quite elaborate and extensive, refer to [Khan and Nielsen, 2018].Algorithm 9 summarizes the VADAM approach for a Gaussian variational posterior with diagonal covariance.
In the diagonal VON update, the Hessian drives the update for the scaling vector s which determines the covariance matrix diag σ 2 .The Hessian can be negative, a situation that could turn σ 2 negative, which is meaningless.Instead of indirectly tackling the issue by using a constrained optimization approach (which could be difficult to implement), such as a controlled adaptive learning rate, or model reparametrization, [Khan et al., 2018b] proposes the use of the Generalized Gauss-Newton approximation for the Hessian: This enables a minor but important difference with respect to VON: with an initial positive value for σ 2 , the above approximation will remain positive leading to valid covariance updates.This provides an algorithmic advantage over VON as constraints on σ 2 are implicitly satisfied.The above implementation of the Hessian estimation, within VON, consists in the Variational Online Gauss-Newton (VOGN) approach [Khan et al., 2018b, Osawa et al., 2019].The implementation of the above approximation is not immediate as it requires per-sample gradients.The approximation averages squared gradients evaluated on a sample-per-sample basis, as opposed to batch-gradient computation which directly computes the sum of the gradients over mini-batches [Osawa et al., 2019], which can be seen by comparing eq. ( 89) with eq. ( 90) The gradient-magnitude approximation that makes use of the mini-batch squared gradient as an approximation for the Hessian, introduces a bias in the Hessian estimation.In fact, increasing the mini-batch size is not advisable as it introduces more bias.Based on the above approximation, [Khan et al., 2018b] advances an RMSProp version of the VON update.
The practical implementation of VOGN is extensively discussed in [Osawa et al., 2019], where the efficient implementation of the per-sample gradient computation for certain network layers is discussed: the additional computations needed to access individual gradients bring the run-time within 2-5 times of that of ADAM.Algorithm 10 summarizes the implementation of the VOGN optimizer.
Algorithm 10 Variational Online Gauss-Newton repeat Randomly sample a mini-batch M of size M for s = 1, . . ., N s do Simulate θ s = µ + εσ, where σ = (1/(N (s + λ))) The form of Algorithm 10 slightly differs from that of VOGN/ADAM.The sampling of the random weights is analogous to that of Algorithm 9 and Algorithm 8, yet here posterior samples are built over standard-normal random numbers rather than directly sampling from the multivariate diagonal posterior by the use of the reparametrization trick.Note the index i referring to the individual samples in the mini-batch M. While VOGN uses a single sample for evaluating the stochastic gradients, here N s draws are averaged to reduce the approximation variance.In particular, the nested for loop computes the single-observation gradient used for the Hessian approximation, each computed in the sampled weight vector θ s .Draw-specific gradients and Hessian ĝs and ĥs are thus averaged across samples (leading to ĝ and ĥ) and used in the implementation of the ADAM-like update based on momentum (thus the hyperparameters β 1 , β 2 ).The pseudo-code in [Osawa et al., 2019] involves an additional tempering parameter and data-augmentation factor along with details for the VOGN parallel implementation, to which we refer for further insights.[Osawa et al., 2019] furthermore discusses practical implementation aspects typical in ML such as batch normalization, data augmentation, momentum, and distributed computing.The feasibility of the VOGN update for large-scale experiments with big-data sizes and deep network architectures on standard datasets promotes VOGN as a state-of-the-art method for Bayesian DL.As a remark, among its limitations, note that VOGN applies to Gaussian variational posteriors with a diagonal covariance matrix only.

Quasi Black-Box Variational Inference (QBVI)
The BBVI framework of [Ranganath et al., 2014] can benefit from the use of the natural gradients.In fact, in [Trusheim et al., 2018] natural gradients are estimated via MC sampling.On the other hand eq. ( 60) provides an exact framework for computing natural gradients without relying on sampling methods, applicable for the wide class of variational posteriors within the Exponential family, yet model-specific derivations, i.e. the computation of the gradients and Hessian, are involved.The Quasi Black-Box Variational Inference (QBVI) approach [Magris et al., 2022c] merges the BBVI setting with the exact natural gradient computation.QBVI uses eq. ( 60) to turn the computation of the natural gradients into Euclidean gradients of the LB, which are computed by the use of the score estimation, resembling the BBVI framework.On a general level, the QBVI update estimates the gradient of the LB with respect to the natural parameters as which along with a plain SGD lead to the update rule Here the exact computation of the natural gradient is carried out in terms of eq. ( 60), so that the QBVI update for a generic variational distribution and prior (both within the exponential family) reads, for the natural parameters, as: Similarly to [Khan and Nielsen, 2018], eq. ( 92) uses the properties of the Exponential family distribution for the prior p, with natural parameter η, and q to simplify the first term on the right-side of eq. ( 91).This results in the natural-parameter difference η − λ, avoiding on the first instance a sampling framework for evaluating the corresponding expectation, i.e. reducing the variance of the estimate for ∇L(λ), regardless of the estimator used for E q λ [log p(D|θ)].
[ Magris et al., 2022c] focuses on the Gaussian variational case, building the QBVI update on the NGVI update, but without using the model's gradient and Hessian as for VON.Indeed, using equations ( 61) and ( 62), [Magris et al., 2022c] recovers, for a full-covariance posterior, the following updates: where v t = Σ −1 (θ − µ t ) and µ 0 , Σ 0 denote the mean vector and covariance matrix of the prior distribution on the model parameter θ, respectively.The following naive MC estimator provides a simple approach for tacking the above expectations with θ s ∼ q λ , s = 1, . . ., N s .Algorithm 11 provides the pseudo-code for the QBVI implementation.

Variational Inference on manifolds
In this section, we review a class of methods that pursue a theoretically different approach, i.e., manifold optimization.The major challenge in VI optimization is that of guaranteeing constraints on the variational parameter.In a Gaussian, or e.g. an Inverse Wishart, setting, this corresponds to guaranteeing updates under which the covariance matrix is Symmetric and Positive Definite (SPD).
We first introduce in general terms the concept and practice of Riemann optimization.Therefore, we provide an introduction to Riemann manifolds, the concepts of tangent vectors, tangent spaces, and Riemann gradient to finally provide a more rigorous discussion of the specific problem of performing valid covariance updated for Bayesian Inference under a Gaussian variational FFVI setting.This section addresses the most crucial aspects concerning the purpose of introducing the Manifold Gaussian Variational Bayes and Exact Manifold Gaussian Variational Bayes optimizers.As the topic is itself broad and quite technical, we intentionally provide a descriptive illustration suitable for a general audience, referring to the specialized literature for additional details and a rigorous mathematical treatment at the end of the following section.

Introduction to manifold optimization
Riemann optimization is an alternative to standard SGD that well fits problems of the kind where L is a real-valued function of some parameter ζ, defined on a Riemannian manifold (M, g).
A manifold is a topological space that locally resembles Euclidean space near each point, in more detail, is a set that can locally be mapped one-to-one to R k , where k is the dimension of the manifold.g stands for a metric the manifold is equipped with.
The optimization problem aims at minimizing L by finding the parameter ζ ∈ M that lies on the "smooth surface" of the Riemannian manifold (M, g) resembling a constrained optimization problem requiring the optimum ζ * to lie on the Riemannian manifold, such as a sphere or the SPD set.As with SGD in Euclidean vector spaces, Riemann optimization is generally tackled with gradient descent on the surface of the manifold, based on the gradients of L. Yet, because of the manifold constraint, there are important differences compared to the standard SGD approach.
The Euclidean vector space R n can be interpreted as a Riemannian manifold (R n , g), with g the common Euclidean metric, where the usual SGD iteratively updates the parameter ζ as where It is clear that applying the above to a generic non-Euclidean manifold M is not trivial as there is no guarantee that ζ t+1 is a valid update, i.e. that ζ t+1 lies in M. Consider an optimization problem where the parameter ζ = (x, y, z) is required to lie on a 2-dimensional spherical manifold of radius 1, embedded in a 3-dimensional ambient space.The Riemannian manifold is

Elements of Riemannian manifolds
In R k , a steepest-ascent approach updates the current iterate ζ in the direction where the first-order increase of the objective function L is most positive.Formally, the update direction is chosen to be the unit norm vector η that minimizes the directional derivative With the domain of L being the manifold M, the argument ζ + tη does not make much sense in general as M is not necessarily a vector space.This leads to the notion of a tangent vector.A possibility for generalizing the directional derivative is to replace t is called the tangent vector to the curve γ at t = 0. Note that the above definition defines γ(0) as a mapping and not as a (e.g.time) derivative as in eq. ( 103), which would be general meaningless.
We can now formally define the notion of a tangent vector.
A tangent vector ξ ζ to a manifold M at a point ζ is a mapping from F ζ (M) to R such that there exists a curve γ on M with γ(0) = ζ satisfying Such a curve γ is said to realize the tangent vector ξ ζ .The tangent space to M at ζ is the set of all tangent vectors to M at ζ and is denoted by T ζ M. Importantly, it can be shown that T ζ M admits a vector space structure, i.e.T ζ M is a vector space: it provides a local vector space approximation of the manifold.This property is useful in defining retractions used to locally transform an optimization problem on M into an optimization problem on the more friendly vector space T ζ M. To characterize which direction of motion from ζ produces the steepest increase in L, to enable a notion of length that applies to tangent vectors, we endow the tangent space T ζ M with an inner product •, • , inducing the norm ||ξ ζ || on T ζ M, from which the direction of the steepest ascent is given by arg max that is, by the unit-norm vector ξ * ζ for which directional derivative D of L in ζ in the direction ξ * ζ is maximized.
A manifold whose tangent spaces are endowed with a smoothly varying inner product is called a Riemannian manifold, and the smoothly varying inner product is called the Riemann metric.With g being such a Riemann metric on M, the Riemannian manifold is, strictly speaking, the couple (M, g).The Euclidean space is the particular Riemannian manifold consisting of a vector space endowed with an inner product.
The gradient of L defined on a Riemannian manifold M at ζ is denoted by the unique element in (106) As in correspondence with usual Euclidean gradients, and important in the light of optimization, it can be shown that the direction of grad L(ζ) is the steepest ascent direction of and that the norm of grad L(ζ) gives the steepest slope of L at ζ.
If a manifold M e is endowed with a Riemann metric, one would expect that manifolds generated from M e inherit its Riemann metric.Let M be a manifold embedded in M e (the subscript e stands for "embedding").Since every tangent space T ζ M can be regarded as a subspace of T ζ M e , the Riemann metric g e of M e induces a Riemann metric g on M turning M into a Riemannian manifold.Endowed with this metric, M is called a Riemannian submanifold of M e .As it will appear clear in the next section, the submanifold idea is simple yet powerful as any element ξ ζ in T ζ M can be decomposed into an element of T ζ M and its corresponding orthogonal element with L e being an extended version of the differentiable function L defined on M e such that its restriction on M actually coincides with L.
Perhaps the most simple tool to tackle Riemann optimization is the Riemann Stochastic Gradient Descent (RSGD), first proposed in [Bonnabel, 2013].RSGD typically involves three steps: The last step moves the point ζ t ∈ M in the direction of the gradient along a geodesic, onto ζ t+1 , lying on the manifold.This is achieved by the so-called exponential map, mapping elements from the tangent space to M. The computation of the exponential map is however a cumber-stone in practice: often a first-order approximation is used.Such first-order approximation is called retraction R Intuitively, rather than performing an exact update following the curved geodesics of the manifold, retraction first follows a straight line in the tangent space and then orthogonally projects the point in the tangent space on the manifold.Closed-form formulae for retraction on the most common manifold are available in the literature, see e.g.[Absil et al., 2009, Hu et al., 2020], and e.g.[Hosseini and Sra, 2015] for the SPD manifold.
The main sources we have used in writing this section are the exhaustive book of Absil et al. [2009], and the articles [Hu et al., 2020] and [Tran et al., 2021a].Classical specialized books on differential geometry are those of Kobayashi andNomizu [1963], Do Carmo andFlaherty Francis [1992], Boothby and Boothby [2003], while well-suited references for readers without a background in abstract topology are e.g.[Tu, 2011, Do Carmo, 2016] and furthermore at an introductory level e.g.[Brickell andClark, 1970, Abraham et al., 2012].We suggest referring to the literature involved in the above references for further bibliographical details, e.g. the bibliographical notes in Chapter 3 of [Absil et al., 2009].An exhaustive overview of the different applications in manifold optimization in different areas can be found e.g. in [Hu et al., 2020].For the first developments on SGD on Riemannian manifolds, we refer to [Bonnabel, 2013], further developments towards an RMSprop-like adaptive version of RSGD can be found in [Kasai et al., 2019], while Riemann optimization on the lines of the popular Adam and Ada-grad are discussed in [Bécigneul and Ganea, 2018].Relevant for the SPD matrix manifold optimization are the results on vector transport and retraction in e.g.[Jeuris et al., 2012] and [Sra and Hosseini, 2015], of remarkable utility for applications.In this regard, we point to [Boumal et al., 2014] for a manifold optimization package available in multiple languages.

Variational Bayes on Riemannian manifolds with natural gradients
Variational Bayes on manifolds aims at maximizing the LB L under a fixed-form Gaussian variational posterior guaranteeing a positive-definite form of the covariance matrix Σ.
That is, the natural gradient of the extended LB L e in M e corresponds to the natural gradient of the LB on the relevant manifold M. A framework for formally associating the natural gradient with the Riemannian gradient is provided by the lemma below, see [Tran et al., 2021a] for more details.Alternatively one can derive the Riemannian gradient of L requiring M to be a so-called quotient manifold induced from a Riemannian ambient manifold.In this regard, see [Tran et al., 2021a] and the references therein.
RSGD requires a proper retraction R ζ : T ζ M → M that locally maps T ζ M onto the manifold M while preserving the first-order information of the tangent space in ζ.This means that a step of size zero stays at the same point ζ and the differential of the retraction at this origin is the identity mapping [Jeuris et al., 2012].From the geodesics between two matrices in M, [Jeuris et al., 2012] develops the popular and convenient retraction method (actually a second-order approximation of the exponential map) for the SPD matrices manifold M.This is given by where We now add a practically important element to the discussion, vector transport.In order to perform, among the others, the conjugate gradient algorithm, or implement the momentum method within the RSGD update, we need to relate a tangent vector at some point ζ ∈ M to another point η ∈ M. In differential geometry, this is achieved by a parallel translation, moving tangent vectors from one tangent space to the other, while preserving the length and angle of the original tangent vector, Figure 4 (right panel).As for the exponential map, the parallel translation is often approximated by the so-called vector transport, which is easier to compute.For ξ := ξ ζ and η := η ζ ∈ T ζ M, an effective vector transform for the manifold of interest is with where T ζ→η (ξ) denotes the vector transport of the tangent vector ξ ∈ T ζ M to η ∈ T η M. The above vector transport can be written in a compact and computationally advantageous form as (see e.g.[Sra and Hosseini, 2015] for details): with We point out that within the above SPD matrix manifold setting relevant in Gaussian VI, ζ, η, ξ are matrices and the above equations are well-defined: for homogeneity in notation, we stick with the lower-case bold symbols for indicating elements of a generic space.
The above vector transport is practically relevant and essential in implementing, e.g., a momentum method on the RSGD update, that is by using a moving average of the Riemannian gradient at the previous iteration to reduce noise in the estimated gradients and boost convergence: where ω is a momentum-weight hyper-parameter.
Manifold optimization in the context of VI is relatively new, the main reference for this paper is [Tran et al., 2021a], whose approach is reviewed in Section 8.3.Besides this, VI on manifolds is also discussed in [Zhou et al., 2021, Magris et al., 2022b] and appears in [Lin et al., 2020].Other applications, not related to the purposes of this review, are here not covered, e.g.manifold optimization for variational autoencoders [Skopek et al., 2019].Regarding the specific Bayesian inference problem for Neural Networks, at the time of writing, we are not aware of any further works or developments.

Manifold Gaussian Variational Bayes (MGVB)
We review the Manifold Gaussian Variational Bayes (MGVB) method of [Tran et al., 2021a].The variational approximation q λ to the true posterior is provided by a multivariate Gaussian distribution N (µ, Σ), µ ∈ R k .The parameter µ, Σ are jointly collected in the vector ζ = (µ, vec(Σ)), denoting the variational parameter.There are no restrictions on the structure of the variance-covariance matrix Σ which is a generic member of the manifold M of the Symmetric and Positive Definite (SPD) matrices, M = Σ ∈ R k×k : Σ = Σ , Σ 0 .
The exact form of the Fisher information matrix for the multivariate normal distribution is, e.g., provided in [Mardia and Marshall, 1984] and reads with I(Σ) being the k 2 × k 2 matrix whose generic element is The MGVB method relies on the approximation I(Σ) ≈ Σ −1 ⊗ Σ −1 , where ⊗ denotes the Kronecker product.The corresponding approximate inverse FIM reads which leads to a convenient approximate form of the natural gradients of the lower bound with respect to µ and Σ computed as ∇µ L = Σ∇ µ L, (126) ∇Σ L ≈ vec −1 (Σ ⊗ Σ)∇ vec(Σ) L = Σ∇ Σ LΣ. (127) The last equality follows from the fact that for a vector v ∈ R k×k , (Σ ⊗ Σ)v = vec Σvec −1 (v)Σ .
In virtue of the natural gradient definition, the first natural gradient for µ is exact while the second one for Σ is approximate.As pointed out in [Lin et al., 2020], the actual natural gradient for the above Gaussian distribution should read 2Σ∇ Σ LΣ, as I(Σ) = 2Σ −1 ⊗ Σ −1 , therefore the MGVB approximation.Thus, [Tran et al., 2021a] adopts the following updates for the parameters of the variational posterior: where R Σ (•) denotes a suitable retraction for Σ on the manifold M, and β is the learning rate.Momentum gradients can be used in place of natural ones.In particular [Tran et al., 2021a] uses retraction in eq. ( 116) and momentum gradients for the updating Σ.In this regard, [Tran et al., 2021a] adopts the parallel transport in eq. ( 118) for granting that at each iteration the weighted gradient remains in the tangent space of the manifold M.
The actual computation of the gradients ∇µ L and ∇Σ L boils down to computing ∇ µ L and ∇ Σ L, which in MGVB is achieved with the black-box estimator where with q ∼ N (µ, Σ), ζ = (µ, vec(Σ)), and L(ζ) ≡ L(µ, Σ).In particular, the gradient of L with respect to ζ is estimated using N s samples from the variational posterior through the unbiased estimator where θ s ∼ N (µ t , Σ t ) and the h-function is evaluated in the current value of the parameters, i.e. in ζ t = (µ t , vec(Σ t )).For a Gaussian distribution q ∼ N (µ, Σ) it can be shown that [e.g.Wierstra et al., 2014, Magris et al., 2022c]: Algorithm 12 summarizes the above process.

Exact Manifold Gaussian Variational Bayes (EMGVB)
The covariance matrix Σ is positive definite, its inverse exists and it is as well symmetric and positive definite.Therefore, Σ −1 lies within the manifold M and can be updated with a suitable retraction algorithm as for Σ in Section 8.3, Algorithm 12 Manifold Gaussian Variational Bayes 1: Simulate: θ s ∼ q µ,Σ , s = 1, . . ., N s 2: Compute: ĝµ , ĝΣ 3: m µ = ĝµ , m Σ −1 = ĝΣ 4: repeat m Σ = T Σold→Σ (m Σ ) + (1 − ω)ĝ Σ 11: until stopping condition is met Opposed to the EMGVB update, relying on the approximation I −1 (Σ) ≈ Σ −1 ⊗ Σ −1 , for tackling a positive-definite update of Σ, [Magris et al., 2022b] targets at updating Σ −1 for which its natural gradient is available in an exact form, by primarily exploiting the duality between the gradients in the natural and expectation parameter space as for eq.( 47), that circumvents the computation and the approximate form of the FIM.
where the weight 0 < ω < 1 is a hyper-parameter.As for eq.( 92), by using a Gaussian prior along with a Gaussian posterior, the natural parameter difference becomes particularly simple.With ζ = (µ, Σ), evaluating ∇ζ L accounts to practically estimating ∇ζ E q ζ [log p(D|θ)] only.Whether or not one uses the results in eq. ( 141) and eq.( 142) under a Gaussian prior assumption, or prefers to use the gradient estimator based on the h-function, h ζ (θ) = E q ζ [log p(θ) + log p(D|θ) − log q ζ (θ)], as in MGVB, a general-form for the gradients enabling the EMGVB update is provided by where if prior is Gaussian or not.
(145) Because of the computations of the constants C Σt and c µt under the Gaussian assumption for the prior p, the MC estimators in eq. ( 143) and eq.( 144) are of reduced variance.[Magris et al., 2022b] also provides analogous simplified updates under the specific assumption that the covariance matrix of q is either diagonal, block-diagonal, or full under an isotropic Gaussian prior whose mean vector is zero and prior covariance matrix Σ −1 0 equal to τ I, with τ > 0. Algorithm 13 summarizes the updating routine.

Conclusion
In this survey, we provided an algorithmic overview of standard, as well as, more recently introduced approaches for Bayesian learning for Neural Networks.We structured our description as an easily-accessible introduction to the basic concepts and related methodologies, focused on the core elements and their implementation, providing pseudo-codes and update rules to be used as references for a large number of Bayesian Neural Network implementations.
We provided a foreword introduction to Bayesian Neural Network, their peculiarities, and motivated their use with respect to standard non-Bayesian Artificial Neural Network.In the remainder, we focused on popular and feasible approaches for their estimation.Besides describing some effective Monte-Carlo methodologies, and introducing Monte Carlo Dropout as a Bayesian tool, we presented a variety of methods based on Variational Inference and natural gradients as the main methodological ingredients in modern Bayesian inference for Neural Networks.We presented the widespread Bayes-By-Backprop optimizer, followed by two common black-box methods, namely Black-Box Variational Inference and Natural-Gradient Black-Box Variational Inference.Next, we introduced natural gradients and examined the Natural-Gradient Variational Inference, Variational Online Newton, Variational Online Gauss-Newton, and Quasi Black-Box Variational Inference approaches.Lastly, by providing an introduction to manifold optimization, we provided a discussion on methods that can implicitly deal with the positive-definite constraint over Gaussian variational specifications, presenting the Manifold Gaussian Variational Bayes and Exact Manifold Gaussian Variational Bayes solutions.
We hope that our comprehensive algorithmic treatment of the above-described methodologies will contribute to a better understanding of the connections and differences between the various Bayesian methods for Neural Networks, will support the adoption of such methods in a wide range of applications, and promote further research in this field.

Figure 1 :
Figure 1: Representation of the operations within the j-th neuron at layer is l.Connections between this neuron and neurons in layer l − 1 are represented by lines corresponding to weights θ l •j .The inputs to the neuron o l−1 • interact with the weights θ l•j , computing the weighted sum a l j .The socalled activation function g(•) is applied to a l j leading to the output o l j , which is sent to nodes at layer l + 1.

Figure 2 :
Figure 2: A feedforward network with multiple layers.Dots represent neurons across different layers (colors).The d-dimensional input vectorx i = [x 1 i , . . ., x d i ]T is sequentially parsed to the output, from left to right, following the connections represented in grey which correspond to the weights of the network's layers.

Figure 3 :
Figure 3: A BNN with multiple layers.Connections correspond to random variables, and outputs here correspond to a tri-variate distribution, whose marginals are represented in the rightmost boxes.
→ ζ + tη by a smooth curve γ on M passing through ζ, i.e. γ(0) = ζ.A smooth mapping γ : R → M : t → γ(t) is termed as curve in M. Defining a derivative γ (t) as γ (t) := lim t→0 γ(ζ+t)−γ(ζ) t fails on a general manifold as it requires a vector space structure to compute the difference γ(ζ + t) − γ(ζ), however for a smooth function L on M the function L • γ : t → L(γ(t)) is a smooth and well-defined function from R to R with a well-defined classical derivative.To sum up, let ζ be a point on M, γ a curve such that γ(0) = ζ and F ζ (M) is the set of smooth real-valued functions defined in a neighborhood of ζ in M. The mapping γ(0) from F ζ (M) to R defined by

where
Proj ζ denotes the orthogonal projection onto T ζ M, and Proj ⊥ ζ denotes the orthogonal projection onto (T ζ M)⊥ .In this light, by properly defining the embedding ambient space M e , one may

Figure 4 :
Figure 4: Left: Tangent space and projection of Riemannian gradient.Center: retraction map.Right: vector transport.
(i) evaluate the gradient of L e in T ζ M e with respect to ζ at the current value ζ t , Figure 4 (left panel), (ii) project the gradient onto the tangent space of the manifold M at ζ t , and (iii) update the parameter by performing a gradient step on the surface following the direction of grad L(ζ), Figure 4 (central panel).

Lemma 1 .
The natural gradient of the function L e on the Riemannian manifold M e with the Fisher-Rao metric is the Riemannian gradient of L e .In particular, the natural gradient at ζ belongs to the tangent space to L e at ζ.This means that with respect to the embedding space M e , ∇ζ L e (ζ) is the actual Riemannian gradient, lying on the tangent space T ζ M e of L e at ζ. Yet we need to associate the Riemannian gradient in M e to the LB L in M, the actual objective of RSGD optimization.To this end, we naturally equip the submanifold M with the same Riemann metric inherited from M e .For ν ζ , ξ ζ now both in T ζ M e , ν ζ , ξ ζ = ν ζ I ζ ξ ζ , (114) and we obtain the Riemannian gradient of L in M as the projection of gradL e on T ζ M gradL(ζ) = Proj ζ gradL e (ζ).(115) In a Gaussian manifold, T ζ M ∼ = T ζ M e , thus the projection is the identity matrix I and grad L e = grad L. Indeed in Gaussian manifolds, M corresponds to the manifold of SPD matrices whereas M e = R k×k : T ζ M and T ζ M e differ by the fact that the first is the tangent space to a certain SPD matrix while the second is the tangent space of a generic k × k symmetric matrix.In terms of projection, the difference is irrelevant, thus Proj ζ grad L e = Igrad L e .Mind that, however, on a general level Proj ζ (•) can be rather difficult to compute.The above relationship between the Riemannian gradient in M e and the LB L in M, is established by treating M as a submanifold of M e .
117) and updates the current value of ζ on M by accounting for ξ on the tangent space T ζ M.

Σ
old = Σ, Σ = R Σ (βm Σ ) 7: Simulate: θ s ∼ q µ,Σ , s = 1, . . ., N s ωm µ + (1 − ω)ĝ µ 10: with g being the Euclidean metric, and L corresponding to a custom loss function for an arbitrary point ζ on the sphere M. Though partial derivatives ∇ ζ L are straightforward to compute or evaluate, e.g. with backpropagation, at the current parameter value ζ t , there is no guarantee that the update rule for the Euclidean space ζ t+1 = ζ t + β∇ ζ Lζ t would result in an updated parameter lying on sphere M. Intuitively, on the "curved" surfaces of Riemannian manifolds the updates should follow the "curved" geodesics instead of straight lines as on familiar R n Euclidean spaces.To this end, Riemann Stochastic Gradient Descent (RSGD) constitutes a manifold generalization of the SGD.
Thus the variational parameter ζ lies on the Riemannian manifold of Symmetric and Positive Definite (SPD) matrices M = Σ ∈ R k×k : Σ = Σ , Σ 0 .The optimization problem of concern is thus the To implement the RSGD update the manifold M of SPD matrices is viewed as embedded in the Riemannian manifold M e .Let T ζ M e be the tangent space to M at ζ ∈ M e .Aligned with the discussion in Section 5, we wish to perform natural gradient updates.To this end, we equip M e with the Fisher-Rao metric, defined by the Fisher information matrix I ζ .With such a metric, the inner product between two tangent vectors ν ζ , ξ ζ ∈ T ζ M e is defined as ν ζ , ξ ζ = ν ζ I ζ ξ ζ , (111) generalizing the usual Euclidean metric ν ζ , ξ ζ = ν ζ ξ ζ .Let L e be a differentiable function defined on M e such that its restriction on M corresponds to the LB L. It can be shown that the steepest ascent direction at ζ ∈ M e for maximizing the objective L e is the natural gradient ∇L e (ζ) = I −1 ζ ∇ ζ L e (ζ), ζ ∈ M e .(112) Note that ∇ ζ L e (ζ) is the usual Euclidean gradient vector of L e (ζ), and that, importantly, for ζ ∈ M, ∇ζ L e