A Scalable Bayesian Sampling Method Based on Stochastic Gradient Descent Isotropization

Abstract

Stochastic gradient sg-based algorithms for Markov chain Monte Carlo sampling (sgmcmc) tackle large-scale Bayesian modeling problems by operating on mini-batches and injecting noise on sgsteps. The sampling properties of these algorithms are determined by user choices, such as the covariance of the injected noise and the learning rate, and by problem-specific factors, such as assumptions on the loss landscape and the covariance of sg noise. However, current sgmcmc algorithms applied to popular complex models such as Deep Nets cannot simultaneously satisfy the assumptions on loss landscapes and on the behavior of the covariance of the sg noise, while operating with the practical requirement of non-vanishing learning rates. In this work we propose a novel practical method, which makes the sg noise isotropic, using a fixed learning rate that we determine analytically. Extensive experimental validations indicate that our proposal is competitive with the state of the art on sgmcmc.

1. Introduction

Stochastic gradient (sg) methods have been extensively studied as a means for mcmc-based Bayesian posterior sampling algorithms to scale to large data regimes. Variants of sg-mcmc algorithms have been studied through the lens of first [1,2,3] or second-order [4,5] Langevin Dynamics, which are mathematically convenient continuous-time processes that correspond to discrete-time gradient methods with and without momentum, respectively. The common traits underlying many methods from the literature can be summarized as follows: they address large data requirements using sg and mini-batching, they inject Gaussian noise throughout the algorithm execution, and they avoid the expensive Metropolis-Hasting accept/reject tests that use the whole data [1,2,4].

Despite mathematical elegance and some promising results restricted to simple models, current approaches fall short in dealing with the complexity of the loss landscape typical of popular modern machine learning models, e.g., neural networks [6,7], for which stochastic optimization poses some serious challenges [8,9].

In general, sg-mcmc algorithms inject random noise to sg descent algorithms: the covariance of such noise and the learning rate, or step-size in the stochastic differential equation simulation community, are tightly related to the assumptions on the loss landscape, which together with the sg noise, determine the sampling properties of these methods [5]. However, current sg-mcmc algorithms applied to popular complex models such as Deep Nets, cannot simultaneously satisfy the assumptions on posterior distribution geometry and on the behavior of the covariance of the sg noise, while operating with the practical requirement of non-vanishing learning rates. In this paper, in accordance with most of the Neural Network related literature, we refer to the posterior distribution geometry as loss landscape. Some recent work [10], instead, argue for fixed step sizes, but settle for variational approximations of quadratic losses. Although we are not the first to highlight these issues, including the lack of a unified notation [5], we believe that studying the role of noise in sg-mcmc algorithms has not received enough attention, and a deeper understanding is truly desirable, as it can clarify how various methods compare. Most importantly, this endeavor can suggest novel and more practical algorithms relying on fewer parameters and less restrictive assumptions.

In this work we chose a mathematical notation that emphasizes the role of noise covariances and learning rate on the behavior of sg-mcmc algorithms (Section 2). As a result, the equivalence between learning rate annealing and extremely large injected noise covariance becomes apparent, and this allows us to propose a novel practical sg-mcmc algorithm (Section 3). We derive our proposal, by first analyzing the case where we inject the smallest complementary noise such that its combined effects with the sg noise result in an isotropic noise. Thanks to this isotropic property of the noise, it is possible to deal with intricate loss surfaces typical of deep models, and produce samples from the true posterior without learning rate annealing. This, however, comes at the expense of cubic complexity matrix operations. We address such issues through a practical variant of our scheme, which employs well-known approximations to the sg noise covariance (see, e.g., [11]). The result is an algorithm that produces approximate posterior samples with a fixed, theoretically derived, learning rate. Please note that in generic Bayesian deep learning setting, none of the existing implementations of sg-mcmc methods converge to the true posterior without learning rate annealing. In contrast, our method automatically determines an appropriate learning rate through a simple estimation procedure. Furthermore, our approach can be readily applied to pre-trained models: after a “warmup” phase to compute sg noise estimates, it can efficiently perform Bayesian posterior sampling.

We evaluate sg-mcmc algorithms (Section 4) through an extensive experimental campaign, where we compare our approach to several alternatives, including Monte Carlo Dropout (mcd) [12] and Stochastic Weighted Averaging Gaussians (swag, [9]), which have been successfully applied to the Bayesian deep learning setting. Our results indicate that our approach offers performance that are competitive to the state of the art, according to metrics that aim at assessing the predictive accuracy and uncertainty.

2. Preliminaries and Related Work

Consider a dataset of m—dimensional observations $\mathcal{D}={\left\{{\mathit{U}}_i\right\}}_{i=1}^N$. Given prior $p\left(\mathit{\theta }\right)$ for a d-dimensional set of parameters, and a likelihood model $p\left(\mathcal{D}\right|\mathit{\theta })$, the posterior is obtained by means of Bayes theorem as follows:

$$p\left(\mathit{\theta }\right|\mathcal{D})=\frac{p\left(\mathcal{D}\right|\mathit{\theta }\left)p\right(\mathit{\theta })}{p\left(\mathcal{D}\right)}$$

(1)

where $p\left(\mathcal{D}\right)$ is also known as the model evidence, defined as the integral $p\left(\mathcal{D}\right)=\int p\left(\mathcal{D}\right|\mathit{\theta }\left)p\right(\mathit{\theta })d\mathit{\theta }$. Except when the prior and the likelihood function are conjugate, Equation (1) is analytically intractable [13]. However, the joint likelihood term in the numerator is typically not hard to compute; this is a key element of many mcmc algorithms, since the normalization constant $p\left(\mathcal{D}\right)$ does not affect the shape of the distribution in any way other than scaling. The posterior distribution is necessary to obtain predictive distributions for new test observations ${\mathit{U}}_{*}$, as:

$$p\left({\mathit{U}}_{*}\right|\mathcal{D})=\int p\left({\mathit{U}}_{*}\right|\mathit{\theta })p\left(\mathit{\theta }\right|\mathcal{D})d\mathit{\theta }$$

(2)

We focus in particular on Monte Carlo methods to obtain an estimate of this predictive distribution, by averaging over $N_{\mathrm{MC}}$ samples obtained from the posterior over $\mathit{\theta }$, i.e., ${\mathit{\theta }}^{\left(i\right)}\sim p\left(\mathit{\theta }\right|\mathcal{D})$

$$p\left({\mathit{U}}_{*}\right|\mathcal{D})\approx \frac{1}{N_{\mathrm{MC}}}\sum _{i=1}^{N_{\mathrm{MC}}}p\left({\mathit{U}}_{*}\right|{\mathit{\theta }}^{\left(i\right)})$$

(3)

We develop our work by working with an unnormalized version of the logarithm of the posterior density, by expressing the negative logarithm of the joint distribution of the dataset $\mathcal{D}$ and parameters $\mathit{\theta }$ as:

$$-f\left(\mathit{\theta }\right)=\sum _{i=1}^{N}logp\left({\mathit{U}}_i\right|\mathit{\theta })+logp\left(\mathit{\theta }\right).$$

(4)

For computational efficiency, we use a minibatch stochastic gradient $\mathit{g}\left(\mathit{\theta }\right)$, which guarantees that the estimated gradient is an unbiased estimate of the true gradient $\nabla f\left(\mathit{\theta }\right)$, and we assume that the randomness due to the minibatch introduces a Gaussian noise:

$$\mathit{g}\left(\mathit{\theta }\right)\sim N(\nabla f(\mathit{\theta }),2\mathit{B}(\mathit{\theta }\left)\right),$$

(5)

where the matrix $\mathit{B}\left(\mathit{\theta }\right)$ denotes the sg noise covariance, which depends on the parametric model, the data distribution and the minibatch size.

A survey of algorithms to sample from the posterior using sg methods can be found in Ma et al. [5]. In Appendix A we report some well-known facts which are relevant for the derivations in our paper. As shown in the literature [10,14], there are structural similarities between sg-mcmc algorithms and stochastic optimization methods, and both can be used to draw samples from posterior distributions. Notice that the original goal of stochastic optimization is to find the minimum of a given cost function, and the stochasticity is introduced by sub-sampling the dataset to scale. sg-mcmc methods instead aim at sampling from a given distribution, i.e., collecting multiple values, and the stochasticity is necessary explore the whole landscape. In what follows, we use a unified notation to compare many existing algorithms in light of the role played by their noise components.

It is well-known [15,16,17] that stochastic gradient descent (sgd), with and without momentum, can be studied through the following stochastic differential equation (sde), when the learning rate $\eta$ is small enough (In this work we do not consider discretization errors. The reader can refer to classical sde texts such as [18] to investigate the topic in greater depth.):

$$d{\mathit{z}}_t=\mathit{s}\left({\mathit{z}}_t\right)dt+\sqrt{2\eta \mathit{D}\left({\mathit{z}}_t\right)}d{\mathit{W}}_t.$$

(6)

where $\mathit{s}$ is usually referred to as driving force and $\mathit{D}$ as diffusion matrix We use a generic form of the sde, with variable $\mathit{z}$ instead of $\mathit{\theta }$, which accommodates sgd variants, with and without momentum. By doing this, we will be able to easily cast the expression for the two cases in what follows (The operator ${\nabla }^{\top }$ applied to matrix $\mathit{D}\left(\mathit{z}\right)$ produces a row vector whose elements are the divergences of the $\mathit{D}\left(\mathit{z}\right)$ columns. Our notation is aligned with Chen et al. [4]).

Definition 1.

A distribution $\rho \left(\mathit{z}\right)\propto exp(-\varphi (\mathit{z}\left)\right)$ is said to be astationarydistribution for thesdeof the form (6), if and only if it satisfies the following Fokker-Planck equation (fpe):

$$\begin{array}{c}\hfill 0=\mathrm{Tr}\left\{\nabla \left[-\mathit{s}{\left(\mathit{z}\right)}^{\top }\rho \left(\mathit{z}\right)+{\nabla }^{\top }\left(\mathit{D}\left(\mathit{z}\right)\rho \left(\mathit{z}\right)\right)\right]\right\}.\end{array}$$

(7)

Please note that in general, the stationary distribution does not converge to the desired posterior distribution, i.e., $\varphi \left(\mathit{z}\right)\ne f\left(\mathit{z}\right)$, as shown by Chaudhari and Soatto [8]. Additionally, given an initial condition for ${\mathit{z}}_t$, its distribution is going to converge to $\rho \left(\mathit{z}\right)$ only for $t\to \infty$. In practice, we observe the sde dynamics for a finite amount of time: then, we declare that the process is approximately in the stationary regime once the potential has reached low and stable values.

Next, we briefly overview known approaches to Bayesian posterior sampling, and interpret them as variants of an sgd process, using the fpe formalism.

2.1. Gradient Methods without Momentum

The generalized updated rule of sgd, described as a discrete-time stochastic process, writes as:

$$\begin{array}{c}\hfill \delta {\mathit{\theta }}_n=-\eta \mathit{P}\left({\mathit{\theta }}_{n-1}\right)(\mathit{g}\left({\mathit{\theta }}_{n-1}\right)+{\mathit{w}}_n),\end{array}$$

(8)

where $\mathit{P}\left({\mathit{\theta }}_{n-1}\right)$ is a user-defined preconditioning matrix, and ${\mathit{w}}_n$ is a noise term, distributed as ${\mathit{w}}_n\sim N(\mathbf{0},2\mathit{C}\left({\mathit{\theta }}_n\right))$, with a user-defined covariance matrix $\mathit{C}\left({\mathit{\theta }}_n\right)$. Then, the corresponding continuous-time sde is [15]:

$$d{\mathit{\theta }}_t=-\mathit{P}\left({\mathit{\theta }}_t\right)\nabla f\left({\mathit{\theta }}_t\right)dt+\sqrt{2\eta \mathit{P}{\left({\mathit{\theta }}_t\right)}^2\sum \left({\mathit{\theta }}_t\right)}d{\mathit{W}}_t.$$

(9)

In this paper we use the symbol n to indicate discrete time, while t for continuous time. We denote by $\mathit{C}\left(\mathit{\theta }\right)$ the covariance of the injected noise and $\sum \left(\mathit{\theta }\right)$ the composite noise covariance. Please note that $\sum \left({\mathit{\theta }}_t\right)=\mathit{B}\left({\mathit{\theta }}_t\right)+\mathit{C}\left({\mathit{\theta }}_t\right)$ combines the sg and the injected noise. Notice that our choice of notation is different from the standard one, in which the starting discrete-time process is in the form $\delta {\mathit{\theta }}_n=-\eta \mathit{P}\left({\mathit{\theta }}_{n-1}\right)\left(\mathit{g}\left({\mathit{\theta }}_{n-1}\right)\right)+{\mathit{w}}_n$. By directly grouping the injected noise with the stochastic gradient we can better appreciate the relationship between annealing the learning rate and extremely large injected noise. Moreover, as will be explained in Section 3, this allows derivation of a new sampling algorithm.

We define the stationary distribution of the sde in Equation (9) as $\rho \left(\mathit{\theta }\right)\propto exp(-\varphi (\mathit{\theta }\left)\right)$. Please note that when $\mathit{C}=\mathbf{0}$, the potential $\varphi \left(\mathit{\theta }\right)$ differs from the desired posterior $f\left(\mathit{\theta }\right)$ [8]. The following theorem, which is an adaptation of known results in light of our formalism, states the conditions for which the noisy sgd converges to the true posterior distribution (proof in Appendix A).

Theorem 1.

Consider dynamics of the form (9) and define the stationary distribution $\rho \left(\mathit{\theta }\right)\propto exp(-\varphi (\mathit{\theta }\left)\right)$. If

$${\nabla }^{\top }\left(\sum {\left(\mathit{\theta }\right)}^{-1}\right)={\mathbf{0}}^{\top }\mathrm{and}\eta \mathit{P}\left(\mathit{\theta }\right)=\sum {\left(\mathit{\theta }\right)}^{-1},$$

(10)

then $\varphi \left(\mathit{\theta }\right)=f\left(\mathit{\theta }\right)$.

Stochastic Gradient Langevin Dynamics (sgld) [1] is a simple approach to satisfy Equation (10); it uses no preconditioning, $\mathit{P}\left(\mathit{\theta }\right)=\mathit{I}$, and sets the injected noise covariance to $\mathit{C}\left(\mathit{\theta }\right)={\eta }^{-1}\mathit{I}$. In the limit for $\eta \to 0$, it holds that $\sum \left(\mathit{\theta }\right)=\mathit{B}\left(\mathit{\theta }\right)+{\eta }^{-1}\mathit{I}\simeq {\eta }^{-1}\mathit{I}$. Then, ${\nabla }^{\top }\left(\sum {\left(\mathit{\theta }\right)}^{-1}\right)=\eta {\nabla }^{\top }\mathit{I}={\mathbf{0}}^{\top }$, and $\eta \mathit{P}\left(\mathit{\theta }\right)=\sum {\left(\mathit{\theta }\right)}^{-1}$. Although sgld succeeds in (asymptotically) generating samples from the true posterior, its mixing rate is unnecessarily slow, due to the extremely small learning rate [2].

An extension to sgld is Stochastic Gradient Fisher Scoring (sgfs) [2], which can be tuned to switch between sampling from an approximate posterior, using a non-vanishing learning rate, and the true posterior, by annealing the learning rate to zero. sgfs uses preconditioning, $\mathit{P}\left(\mathit{\theta }\right)\propto \mathit{B}{\left(\mathit{\theta }\right)}^{-1}$. In practice, however, $\mathit{B}\left(\mathit{\theta }\right)$ is ill conditioned for complex models such as deep neural networks. Then, many of its eigenvalues are almost zero [8], and computing $\mathit{B}{\left(\mathit{\theta }\right)}^{-1}$ is problematic. An in-depth analysis of sgfs reveals that conditions (10) would be met with a non-vanishing learning rate only if, at convergence, ${\nabla }^{\top }\left(\mathit{B}{\left(\mathit{\theta }\right)}^{-1}\right)={\mathbf{0}}^{\top }$, which would be trivially true if $\mathit{B}\left(\mathit{\theta }\right)$ was constant. However, recent work [6,7] suggest that this condition is difficult to justify for deep neural networks.

The Stochastic Gradient Riemannian Langevin Dynamics (sgrld) algorithm [3] extends sgfs to the setting in which ${\nabla }^{\top }\left(\mathit{B}{\left(\mathit{\theta }\right)}^{-1}\right)\ne {\mathbf{0}}^{\top }$. The process dynamic is adjusted by adding the term ${\nabla }^{\top }\left(\mathit{B}{\left(\mathit{\theta }\right)}^{-1}\right)$. However, the term ${\nabla }^{\top }\left(\mathit{B}{\left(\mathit{\theta }\right)}^{-1}\right)$ has not a clear estimation procedure, restricting sgrld to cases where it can be computed analytically.

The work by [10] investigates constant-rate sgd (with no injected noise), and determines analytically the learning rate and preconditioning that minimize the Kullback–Leibler (kl) divergence between an approximation and the true posterior. Moreover, it shows that the preconditioning used in sgfs is optimal, in the sense that it converges to the true posterior, when $\mathit{B}\left(\mathit{\theta }\right)$ is constant and the true posterior has a quadratic form.

In summary, to claim convergence to the true posterior distribution, existing approaches require either vanishing learning rates or assumptions on the sg noise covariance that are difficult to verify in practice, especially when considering deep models. We instead propose a novel practical method that induces isotropic sg noise and thus satisfies Theorem 1. We determine analytically a fixed learning rate, and we require weaker assumptions on the loss shape.

2.2. Gradient Methods with Momentum

Momentum-corrected methods emerge as a natural extension to sgd approaches. The general set of update equations for (discrete-time) momentum-based algorithms is:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\delta {\mathit{\theta }}_n=\eta \mathit{P}\left({\mathit{\theta }}_{n-1}\right){\mathit{M}}^{-1}{\mathit{r}}_{n-1}\hfill \\ \delta {\mathit{r}}_n=-\eta \mathit{A}\left({\mathit{\theta }}_{n-1}\right){\mathit{M}}^{-1}{\mathit{r}}_{n-1}-\eta \mathit{P}\left({\mathit{\theta }}_{n-1}\right)(\mathit{g}\left({\mathit{\theta }}_{n-1}\right)+{\mathit{w}}_n),\hfill \end{array}\right.\end{array}$$

where $\mathit{P}\left({\mathit{\theta }}_{n-1}\right)$ is a preconditioning matrix, M is the mass matrix and $\mathit{A}\left({\mathit{\theta }}_{n-1}\right)$ is the friction matrix, as shown by [4,19]. As with the first order counterpart, the noise term is distributed as ${\mathit{w}}_n\sim N(\mathbf{0},2\mathit{C}\left({\mathit{\theta }}_n\right))$. Then, the sde to describe continuous-time system dynamics is:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}d{\mathit{\theta }}_t=\mathit{P}\left({\mathit{\theta }}_t\right){\mathit{M}}^{-1}{\mathit{r}}_{t}dt\hfill \\ d{\mathit{r}}_t=-(A\left({\mathit{\theta }}_t\right){\mathit{M}}^{-1}{r}_t+\mathit{P}\left({\mathit{\theta }}_t\right)\nabla f\left({\mathit{\theta }}_t\right))dt+\sqrt{2\eta \mathit{P}{\left({\mathit{\theta }}_t\right)}^2\sum \left({\mathit{\theta }}_t\right)}d{\mathit{W}}_t.\hfill \end{array}\right.\end{array}$$

(11)

where $\mathit{P}{\left({\mathit{\theta }}_t\right)}^2=\mathit{P}\left({\mathit{\theta }}_t\right)\mathit{P}\left({\mathit{\theta }}_t\right)$ and we assume $\mathit{P}\left({\mathit{\theta }}_t\right)$ to be symmetric. The theorem hereafter describes the conditions for which noisy sgd with momentum converges to the true posterior distribution (Appendix A).

Theorem 2.

Consider dynamics of the form (11) and define the stationary distribution for ${\mathit{\theta }}_t$ as $\rho \left(\mathit{\theta }\right)\propto exp(-\varphi (\mathit{\theta }\left)\right)$. If

$${\nabla }^{\top }\mathit{P}\left(\mathit{\theta }\right)={\mathbf{0}}^{\top }\mathrm{and}\mathit{A}\left(\mathit{\theta }\right)=\eta \mathit{P}{\left(\mathit{\theta }\right)}^2\sum \left(\mathit{\theta }\right),$$

(12)

then $\varphi \left(\mathit{\theta }\right)=f\left(\mathit{\theta }\right)$.

In the naive case, where $\mathit{P}\left(\mathit{\theta }\right)=\mathit{I},\mathit{A}\left(\mathit{\theta }\right)=\mathbf{0},\mathit{C}\left(\mathit{\theta }\right)=\mathbf{0}$, Equation (12) are not satisfied and the stationary distribution does not correspond to the true posterior [4]. To generate samples from the true posterior it is sufficient to set $\mathit{P}\left(\mathit{\theta }\right)=\mathit{I},\mathit{A}\left(\mathit{\theta }\right)=\eta \mathit{B}\left(\mathit{\theta }\right),\mathit{C}\left(\mathit{\theta }\right)=\mathbf{0}$ (as in Equation (9) in [4]).

Stochastic Gradient Hamiltonian Monte Carlo (sghmc) [4] suggests that estimating $\mathit{B}\left(\mathit{\theta }\right)$ can be costly. Hence, the injected noise $\mathit{C}\left(\mathit{\theta }\right)$ is chosen such that $\mathit{C}\left(\mathit{\theta }\right)={\eta }^{-1}\mathit{A}\left(\mathit{\theta }\right)$, where $\mathit{A}\left(\mathit{\theta }\right)$ is user-defined. When $\eta \to 0$, the following approximation holds: $\sum \left(\mathit{\theta }\right)\simeq \mathit{C}\left(\mathit{\theta }\right)$. It is then trivial to check that conditions (12) hold without the need for explicitly estimating $\mathit{B}\left(\mathit{\theta }\right)$. A further practical reason to avoid setting $\mathit{A}\left(\mathit{\theta }\right)=\eta \mathit{B}\left(\mathit{\theta }\right)$ is that the computational cost for the operation $\mathit{A}\left({\mathit{\theta }}_{n-1}\right){\mathit{M}}^{-1}{\mathit{r}}_{n-1}$ has $\mathcal{O}\left(D^2\right)$ complexity, whereas if $\mathit{C}\left(\mathit{\theta }\right)$ is diagonal, this is reduced to $\mathcal{O}\left(D\right)$. This, however, severely slows down the sampling process.

Stochastic Gradient Riemannian Hamiltonian Monte Carlo (sgrhmc) is an extension to sghmc [5]), which considers a generic, space-varying preconditioning matrix $\mathit{P}\left(\mathit{\theta }\right)$ derived from information geometric arguments [20]. sgrhmc suggests setting $\mathit{P}\left(\mathit{\theta }\right)=\mathit{G}{\left(\mathit{\theta }\right)}^{-\frac{1}{2}}$, where $\mathit{G}\left(\mathit{\theta }\right)$ is the Fisher Information matrix. To meet the requirement ${\nabla }^{\top }\mathit{P}\left(\mathit{\theta }\right)={\mathbf{0}}^{\top }$, it includes a correction term, $-{\nabla }^{\top }\mathit{P}\left(\mathit{\theta }\right)$. The injected noise is set to $\mathit{C}\left(\mathit{\theta }\right)={\eta }^{-1}\mathit{I}-\mathit{B}\left(\mathit{\theta }\right)$, consequently $\sum ={\eta }^{-1}\mathit{I}$, and the friction matrix is set to $\mathit{A}\left(\mathit{\theta }\right)=\mathit{P}{\left(\mathit{\theta }\right)}^2$. With all these choices, Theorem 2 is satisfied. Although appealing, the main drawbacks of this method are the need for an analytical expression of ${\nabla }^{\top }\mathit{P}\left(\mathit{\theta }\right)$, and the assumption for $\mathit{B}\left(\mathit{\theta }\right)$ to be known.

From a practical standpoint, momentum-based methods suffer from the requirement to tune many hyperparameters, including the learning rate, and the parameters that govern the simulation of a second-order Langevin dynamics.

The method we propose in this work can be applied to momentum-based algorithms; in this case, it could be viewed as an extension of the work in [11], albeit addressing the complex loss landscapes typical of deep neural networks. However, we leave this avenue of research for future work.

3. Sampling by Layer-Wise Isotropization

We present a simple and practical approach to inject noise to SGD iterates to perform Bayesian posterior sampling. Our goal is to sample from the true posterior distribution (or approximations thereof) using a constant learning rate, and to rely on more lenient assumptions about the shape of the loss landscape that characterize deep models, compared to previous works. In general, in modern machine learning applications, we deal with multi-layer neural networks [21]. We exploit the natural subdivision of the parameters of these architecture into different layers to propose a practical sampling scheme

Careful inspection of Theorem 1 reveals that the matrices $\mathit{P}\left(\mathit{\theta }\right),\sum \left(\mathit{\theta }\right)$ are instrumental in determining the convergence properties of sg methods to the true posterior. Therefore, we consider the constructive approach of designing $\eta \mathit{P}\left(\mathit{\theta }\right)$ to obtain a sampling scheme that meets our goals; we set $\eta \mathit{P}\left(\mathit{\theta }\right)$ to be a constant, diagonal matrix which we constrain to be layer-wise uniform:

$$\eta \mathit{P}\left(\mathit{\theta }\right)={\mathrm{\Lambda }}^{-1}=\mathrm{diag}{\left([\underset{\mathrm{layer}1}{\underbrace{{\lambda }^{\left(1\right)},\dots ,{\lambda }^{\left(1\right)}}},\dots ,\underset{\mathrm{layer}{N}_l}{\underbrace{{\lambda }^{\left(N_l\right)},\dots {\lambda }^{\left(N_l\right)}}}]\right)}^{-1}.$$

(13)

By properly selecting the set of parameters $\left\{{\lambda }^i\right\}$ we can achieve the simultaneous result of non-vanishing learning rate and well-conditioned preconditioning matrix. This implies a layer-wise learning rate ${\eta }^{\left(p\right)}=\frac{1}{{\lambda }^{\left(p\right)}}$ for the p-th layer, without further preconditioning.

We can now prove (see Appendix B), as a corollary to Theorem 1, that our design choices can guarantee convergence to the true posterior distribution.

Corollary 1.

(Theorem 1) Consider dynamics of the form (9) and define the stationary distribution $\rho \left(\mathit{\theta }\right)\propto exp(-\varphi (\mathit{\theta }\left)\right)$. If $\eta \mathit{P}\left(\mathit{\theta }\right)={\mathrm{\Lambda }}^{-1}$ as in (13), $\mathit{C}\left(\mathit{\theta }\right)=\mathrm{\Lambda }-\mathit{B}\left(\mathit{\theta }\right)$ and $\mathit{C}\left(\mathit{\theta }\right)\succ 0\forall \mathit{\theta }$, then $\varphi \left(\mathit{\theta }\right)=f\left(\mathit{\theta }\right)$.

If aforementioned conditions are satisfied, it is in fact simple to show that the relevant matrices satisfy the conditions in Equation (10). The covariance matrix of the composite noise is said to be isotropic within the layers of (deep) models. In fact, $\sum \left(\mathit{\theta }\right)=\mathit{C}\left(\mathit{\theta }\right)+\mathit{B}\left(\mathit{\theta }\right)=\mathrm{diag}\left(\left[\begin{array}{c}{\lambda }^{\left(1\right)},\dots ,{\lambda }^{\left(1\right)},\dots ,{\lambda }^{\left(N_l\right)},\dots {\lambda }^{\left(N_l\right)}\end{array}\right]\right)$. From a practical point of view, we choose $\mathrm{\Lambda }$ to be, among all valid matrices satisfying $\mathrm{\Lambda }-\mathit{B}\left(\mathit{\theta }\right)\succ 0$, the smallest (the one with the smallest $\lambda$’s). Indeed, larger $\mathrm{\Lambda }$ induce a smaller learning rate, thus unnecessarily reducing sampling speed.

Now, let us consider an ideal case, in which we assume the sg noise covariance $\mathit{B}\left(\mathit{\theta }\right)$ and $\mathrm{\Lambda }$ to be known in advance. The procedure described in Algorithm 1 illustrates a naive sg method that uses the injected noise covariance $\mathit{C}\left(\mathit{\theta }\right)$ to sample from the true posterior.

Algorithm 1 Idealized posterior sampling

{Initialization: ${\mathit{\theta }}_0$}

$\mathbf{\hspace{1em}SAMPLE} ({\mathit{\theta }}_0, \mathit{B}(\mathit{\theta }), \mathrm{\Lambda }):$

$\mathit{\hspace{1em}\theta } \leftarrow {\mathit{\theta }}_0$

loop

$g = \nabla \stackrel{\sim }{f}\left(\mathit{\theta }\right)$

$n \sim N(0,\mathit{I})$

$\mathit{C}{\left(\mathit{\theta }\right)}^{1/2} \leftarrow {(\sum - \mathit{B}(\mathit{\theta }\left)\right)}^{1/2}$

$g \leftarrow {\sum }^{-1}(g+\sqrt{2}\mathit{C}{\left(\mathit{\theta }\right)}^{1/2}n)$

$\mathit{\theta } \leftarrow \mathit{\theta }-g$

end loop

This deceivingly simple procedure generate samples from the true posterior, with a non-vanishing learning rate, as shown earlier. However, it cannot be used in practice as $\mathit{B}\left(\mathit{\theta }\right)$ and $\mathrm{\Lambda }$ are unknown. Furthermore, the algorithm requires computationally expensive operations, i.e., to compute ${(\sum -\mathit{B}\left(\mathit{\theta }\right))}^{\frac{1}{2}}$, which requires $\mathcal{O}\left(d^3\right)$ operations, and $\mathit{C}{\left(\mathit{\theta }\right)}^{\frac{1}{2}}$, which costs $\mathcal{O}\left(d^2\right)$ multiplications.

Next, we describe a practical variant of our approach, where we use approximations at the expense of generating samples from the true posterior distribution. We note that [10] suggest exploring a related preconditioning, but do not develop this path in their work. Moreover, the proposed method shares similarities with a scheme proposed in [22] although the analysis we perform here is different.

3.1. A Practical Method: Isotropic SGD

To render the idealized sampling method practical, it is necessary to consider some additional assumptions. As we explain at the end of this section, the assumptions that follow are less strict than other approaches in the literature.

Assumption 1.

Thesgnoise covariance $\mathit{B}\left(\mathit{\theta }\right)$ can be approximated with a diagonal matrix, i.e., $\mathit{B}\left(\mathit{\theta }\right)=\mathrm{diag}\left(\mathit{b}\right(\mathit{\theta }\left)\right)$.

Assumption 2.

The signal-to-noise ratio (SNR) of a gradient is small enough such that in the stationary regime, the second-order moment of the gradient is a good estimate of the true variance. Hence, combining with Assumption 1, $\mathit{b}\left(\mathit{\theta }\right)\simeq \frac{\mathrm{E}\left[\mathit{g}\right(\mathit{\theta })\odot \mathit{g}(\mathit{\theta }\left)\right]}{2}$, where ⊙ indicates the elementwise product.

Assumption 3.

The sum of the variances of noise components, layer by layer, can be assumed to constant in the stationary regime. Then, ${\beta }^{\left(p\right)}=\sum _{j\in I_p}{b}_j\left(\mathit{\theta }\right)$, where $I_p$ is the set of indices of parameters belonging to $p_{th}$ layer.

The diagonal covariance assumption (i.e., Assumption 1) is common in other works, such as [2,11]. The small signal-to-noise ratio as stated in Assumption 2 is in line with recent studies, such as [11,23]. Assumption 3 is similar to those appeared in earlier work, such as [24]. Please note that Assumptions 2 and 3 must hold in the stationary regime when the process reaches the bottom valley of the loss landscape. The matrix $\left(b\right(\mathit{\theta }\left)\right)$ has been associated in the literature with the empirical Fisher information matrix [2,25]. As we do not consider this matrix for preconditioning purposes, we do not further investigate this connection.

Given our assumptions, and our design choices, it is then possible to show (see Appendix B) that the optimal (i.e., the smallest possible) $\mathrm{\Lambda }=\left[\begin{array}{c}{\lambda }^{\left(1\right)},\dots ,{\lambda }^{\left(1\right)},\dots ,{\lambda }^{\left(N_l\right)},\dots {\lambda }^{\left(N_l\right)}\end{array}\right]$ satisfying Corollary 1 can be obtained as ${\lambda }^{\left(p\right)}={\beta }^{\left(p\right)}$. Please note that we do not assume $\mathit{B}\left(\mathit{\theta }\right)$ to be known, but use a simple procedure to estimate its components by computing: ${\lambda }^{\left(p\right)}={\displaystyle \sum _{j\in I_p}}{b}_j\left(\mathit{\theta }\right)=\frac{\left|\right|{\mathit{g}}^{\left(p\right)}\left(\mathit{\theta }\right){\left|\right|}^2}{2},$ where ${\mathit{g}}^{\left(p\right)}\left(\mathit{\theta }\right)$ is the portion of stochastic gradient corresponding to the p-th layer. Then, the composite noise matrix $\sum =\mathrm{\Lambda }$ is a layer-wise isotropic covariance matrix, which inspires the name of our proposed method as Isotropic SGD (i-sgd).

The practical implementation of i-sgd is shown in Algorithm 2. The advantage of i-sgd is that it can either be used to obtain posterior samples starting from a pre-trained model, or do so by training a model from scratch. In either case, the estimates of $\mathit{B}\left(\mathit{\theta }\right)$ are used to compute $\mathrm{\Lambda }$, as discussed above. An important consideration is that once all ${\lambda }^{\left(i\right)}$ have been estimated, the learning rate, layer by layer, is determined automatically. In fact, for the p-th layer, the learning rate is: ${\eta }^{\left(p\right)}={{\lambda }^{\left(p\right)}}^{-1}$. A simpler approach is to use a unique learning rate for all layers, where the equivalent $\lambda$ is the sum of all ${\lambda }^{\left(p\right)}$.

Algorithm 2i-sgd: practical posterior sampling

$\mathbf{\hspace{1em}SAMPLE} \left({\mathit{\theta }}_0\right):$

$\mathit{\hspace{1em}\theta } \leftarrow {\mathit{\theta }}_0$

loop

$g = \nabla \stackrel{\sim }{f}\left(\mathit{\theta }\right)$

for $p \leftarrow 1$ to Nl do

$n \sim N(0,\mathit{I})$

$\mathit{C}{\left(\mathit{\theta }\right)}^{1/2} \leftarrow ({\lambda }^{\left(p\right)} - (1/2) (g^{\left(p\right)}\odot g^{\left(p\right)}\left)\right)$

$g^{\left(p\right)} \leftarrow 1/{\lambda }^{\left(p\right)}(g^{\left(p\right)} + \sqrt{2}\mathit{C}{\left(\mathit{\theta }\right)}^{1/2}n)$

end for

$\mathit{\theta } \leftarrow \mathit{\theta }-g$

end loop

A Remark on Convergence

In summary, i-sgd is a practical method to perform approximate Bayesian posterior sampling, backed up by solid theoretical foundations. Our assumptions, which are at the origin of the approximate nature of i-sgd, are less strict than those used in the literature of sg-mcmc methods. More precisely, the theory behind i-sgd can explain convergence to the true posterior with a non-vanishing learning rate in the particular case when Assumption 1 holds and the estimation of $\mathit{B}\left(\mathit{\theta }\right)$ is perfect. Even with perfect estimates, this is not the case for sgfs, which requires the correction term ${\nabla }^{\top }\mathit{B}{\left(\mathit{\theta }\right)}^{-1}=0$. Additionally, both sgrld and sgrhmc are more demanding than i-sgd because they require computing ${\nabla }^{\top }\mathit{B}{\left(\mathit{\theta }\right)}^{-1}$, for which an estimation procedure is elusive. Finally, the method by Springenberg et al. [11] needs a constant, diagonal $\mathit{B}\left(\mathit{\theta }\right)$, a condition that does not necessarily hold for deep models.

3.2. Computational Cost

The computational cost of i-sgd is as follows. As with [4], we define the cost of computing a gradient minibatch as $C_{\mathit{g}}(N_b,d)$. Thanks to Assumptions 1 and 2, the computational cost for estimating the noise covariance scales as $\mathcal{O}\left(d\right)$ multiplications. The computational cost of generating random samples with the desired covariance scales as $\mathcal{O}\left(d\right)$ square roots and $\mathcal{O}\left(d\right)$ multiplications (without considering the cost of generating random numbers). The overall cost of our method is the sum of the above terms. Notice that the cost of estimating the noise covariance does not depend on the minibatch size $N_b$. We would like to stress that in many modern models, the real computational bottleneck is the backward propagation for the computation of the gradients. As all the sg-mcmc methods considered in this work require one gradient evaluation per step, the different methods have in practice the same complexity.

The space complexity of i-sgd is the same as sghmc,sgfs and variants: it scales as $\mathcal{O}\left(N_{\mathrm{sam}}d\right)$, where $N_{\mathrm{sam}}$ is the number of posterior samples.

4. Experiments

The empirical analysis of our method, and its comparison to alternative approaches from the literature, is organized as follows. First, we proceed with a validation of i-sgd using the standard uci datasets [26] and a shallow neural network. Then we move to the case of deeper models: we begin with a simple CNN used on the mnist [27] dataset, then move to the standard ResNet-18 [28] deep network using the Cifar-10 [29] dataset.

We compare i-sgd to other Bayesian sampling methods such as sghmc [4], sgld [2], and to alternative approaches to approximate Bayesian inference, including mcd [12], swag [9] and vsgd [10]. In general, our result indicates that: (1) i-sgd achieves similar or superior performance regarding competitors, when measuring uncertainty quantification, even with simple datasets and models; (2) i-sgd is simple to tune, when compared to alternatives; (3) i-sgd is competitive when used for deep Bayesian modeling, even when compared to standard methods used in the literature. In particular, the proposed method shares some of the strengths of vsgd, such as learning rates determined automatically and the simplicity of sgld. Appendix B includes additional implementation details on i-sgd. Appendix C presents detailed configurations of all methods we compare, and additional experimental results.

4.1. A Disclaimer on Performance Characterization

It is important to stress a detail on the analysis of the experimental campaign. The discussion is usually focused on the goodness of the various methods for representing the true posterior distribution. Different methods can or cannot claim convergence to the true posterior according to certain assumptions and the nature of the hyperparameters. In the experimental validation of the results, however, we do not have access to the form of the true posterior as it is exactly the problem we are trying to solve. The practical solution adopted is to compare the different methods in terms of proxy metrics evaluated on the test sets, such as the accuracy and uncertainty metrics. Being better in terms of these performance metrics does not imply that the sampling method is better at approximating the posterior distribution, and outperforming competitors in terms of these metric do not provide sufficient information about the intrinsic quality of the sampling scheme.

4.2. Regression Tasks, with Simple Models

We consider several regression tasks defined on the uci datasets. We use a simple neural network configuration with two fully connected layers and a relu activation function; the hidden layer includes 50 units. In this set of experiments, we use the following metrics: the root mean square error (rmse) to judge the model predictive performance and the mean negative log-likelihood (mnll) as a proxy for uncertainty quantification. We note that the task of tuning our competitors was far from trivial. We used our own version of sghmc, based on [11], to ensure a proper understanding of the implementation internals, and we proceeded with a tuning process to find appropriate values for the numerous hyperparameters. In this set of experiments, we omit results for swag, which we keep for more involved scenarios.

Table 1. rmse results for regression on uci datasets.

Method	wine	protein	naval	kin8nm	power	boston
sgld	0.759 ± 0.07	5.687 ± 0.05	0.007 ± 0.00 (F = 6.000)	0.171 ± 0.07 (F = 3.000)	11.753 ± 3.25	9.602 ± 2.06
i-sgd	0.635 ± 0.05	4.699 ± 0.03	0.001 ± 0.00	0.079 ± 0.00	4.320 ± 0.13	3.703 ± 1.19
Baseline	0.641 ± 0.05	4.733 ± 0.05	0.001 ± 0.00	0.080 ± 0.00	4.354 ± 0.12	3.705 ± 1.19
vsgd	0.635 ± 0.05	4.699 ± 0.03	0.001 ± 0.00	0.079 ± 0.00	4.325 ± 0.13	3.588 ± 1.06 (F = 1.000)
sghmc	0.628 ± 0.04	4.712 ± 0.03	0.000 ± 0.00 (F = 2.000)	0.076 ± 0.00 (F = 1.000)	4.310 ± 0.14	3.659 ± 1.24
sgld T	0.752 ± 0.07	5.673 ± 0.04	0.007 ± 0.00 (F = 6.000)	0.169 ± 0.07 (F = 3.000)	11.351 ± 3.02	9.417 ± 2.07
drop	0.637 ± 0.04	4.968 ± 0.05	0.003 ± 0.00	0.139 ± 0.01	4.531 ± 0.16	3.803 ± 1.26
sghmc T	0.628 ± 0.04	4.684 ± 0.03	0.000 ± 0.00 (F = 6.000)	0.076 ± 0.00	4.326 ± 0.13	3.692 ± 1.19

and

Table 2. mnll results for regression on uci datasets.

Method	wine	protein	naval	kin8nm	power	boston
sgld	1.546 ± 0.25	5.604 ± 0.08	−1.751 ± 0.28 (F = 6.000)	5.140 ± 7.05 (F=3.000)	8.429 ± 3.14	30.386 ± 15.77
i-sgd	1.129 ± 0.15	4.371 ± 0.03	−2.466 ± 1.12	−0.460 ± 0.65	3.122 ± 0.07	9.799 ± 5.69
Baseline	1.182 ± 0.03	3.964 ± 0.04	0.920 ± 0.00	0.924 ± 0.00	3.071 ± 0.06	5.421 ± 2.73
vsgd	1.128 ± 0.15	4.371 ± 0.03	−2.466 ± 1.12	−0.480 ± 0.65	3.088 ± 0.06	8.413 ± 5.89 (F = 1.000)
sghmc	1.041 ± 0.12	4.142 ± 0.02	−2.763 ± 1.33 (F = 2.000)	−0.798 ± 0.39 (F = 1.000)	2.924 ± 0.04	3.097 ± 0.83
sgld T	1.526 ± 0.24	5.591 ± 0.07	−1.752 ± 0.28 (F = 6.000)	5.118 ± 7.06 (F = 3.000)	8.288 ± 3.04	33.212 ± 19.69
drop	1.065 ± 0.12	4.218 ± 0.06	−2.322 ± 0.75	−0.086 ± 0.41	2.941 ± 0.04	3.989 ± 1.23
sghmc T	1.104 ± 0.14	4.191 ± 0.02	−2.966 ± 1.89 (F = 6.000)	−0.756 ± 0.42	3.116 ± 0.07	9.826 ± 5.72

report a complete overview of our results, for a selection of uci datasets. For each method and each dataset, we also included how many out of the 10 splits considered failed to converge, indicated as $F=\dots$. As explained in Appendix C we implemented a temperature scaled version of vsgd. A clear picture emerges from this first set of experiments: while for the rmse the performance is similar for different methods, for the mnll averaging over multiple samples clearly improves the uncertainty quantification capabilities. sghmc is in many cases better than alternatives, considering however the standard deviation of the results it is difficult to claim clear superiority of one method over the others.

4.3. Classification Tasks, with Deeper Models

Next, we compare i-sgd against competitors on image classification tasks. First, we use the mnist dataset, and a simple LeNet-5 CNN [30]. All methods are compared based on the test accuracy acc,mnll and the expected calibration error (ece, 31]). Additionally, at test time, we carry out predictions on both mnist and not-mnist; the latter is a dataset equivalent to mnist, but it represents letters rather than numbers. (http://yaroslavvb.blogspot.com/2011/09/notmnist-dataset.html, accessed on 24 October 2021) This experimental setup is often used to check whether the entropy of the predictions on not-mnist is higher than the entropy of the predictions on mnist (the entropy of the output of an $N_{cl}$ classes classifier, represented by the vector $\mathbf{p}$, is defined as $-{\displaystyle \sum _{i=1}^{N_{cl}}}{p}_{i}log{p}_i$).

Table 3. Results for classification on mnist dataset.

Method	acc	mnll	Mean ${\mathit{H}}_0$	ece	Mean ${\mathit{H}}_1$	Failed
i-sgd	9916.3333 ± 2.8674	263.5311 ± 16.3600	0.0368 ± 0.0019	0.0491 ± 0.0003	0.4558 ± 0.0591	0.0000
sghmc	9930.6667 ± 2.4944	268.2559 ± 6.8172	0.0593 ± 0.0018	0.0531 ± 0.0003	1.0369 ± 0.0346	0.0000
drop	9912.6667 ± 6.0185	362.8973 ± 24.8881	0.0960 ± 0.0090	0.0541 ± 0.0011	0.5507 ± 0.0577	0.0000
baseline	9886.6667 ± 11.0252	352.6640 ± 20.8622	0.0353 ± 0.0058	0.0468 ± 0.0001	0.0019 ± 0.0003	0.0000
baseline r	9919.0000 ± 9.4163	242.7644 ± 17.0736	0.0303 ± 0.0001	0.0482 ± 0.0006	0.0021 ± 0.0002	0.0000
swag	9917.0000 ± 2.8284	308.8182 ± 20.0979	0.0675 ± 0.0108	0.0524 ± 0.0011	0.3953 ± 0.0442	0.0000
sgld	9927.0000 ± 1.0000	279.7685 ± 16.6563	0.0556 ± 0.0034	0.0531 ± 0.0004	1.3032 ± 0.1942	1.0000
vsgd	9927.3333 ± 6.7987	225.3725 ± 16.3739	0.0274 ± 0.0008	0.0481 ± 0.0005	0.0414 ± 0.0070	0.0000
i-sgd T	9915.6667 ± 0.9428	255.9641 ± 12.8051	0.0289 ± 0.0014	0.0478 ± 0.0002	0.0284 ± 0.0122	0.0000
sghmc T	9937.0000 ± 0.0000	231.5332 ± 0.0000	0.0434 ± 0.0000	0.0518 ± 0.0000	0.4623 ± 0.0000	2.0000

indicates that all methods are essentially equivalent in terms of accuracy and mnll. We consider, together with the classical in and out of distribution entropies the regions of convergence (rocs) diagrams comparing detection of out of distribution samples and false alarms when using as test statistic the entropy. Results, reported in Figure 1, clearly shows that: (1) collecting multiple samples improve the uncertainty quantification capabilities (2) i-sgd is competitive (but not the best scheme) and importantly outperform the closest approach to ours, i.e., vsgd. The experimental results show that i-sgd improves the quality of the baseline model with respect to all metrics. To test whether the improvements are due just to “additional training” or are intrinsically due to the Bayesian averaging properties, we do consider alternative deterministic baselines (details in Appendix C). For this set of experiments the best performing is baseline R. As can be appreciated by comparing Table 3 and Figure 1, while it is possible to increase the classical metrics, i-sgd (and other methods) still outperform by a large margin the baselines in terms of detection of out of distribution samples.

We now move on to a classical image classification problem with deep convolutional networks, whereby we use the cifar10 dataset, and the ResNet-18 network architecture. For this set of experiments, we compare i-sgd, sghmc, swag, and vsgd using again test accuracy and mnll, which we report in

Table 4. Results for classification on cifar1010 dataset.

Method	acc	mnll	mean ${\mathit{H}}_0$	ece
i-sgd	8591.3333 ± 17.4611	4393.3557 ± 107.0878	0.6107 ± 0.0337	0.0731 ± 0.0075
sghmc	8634.6667 ± 5.1854	4357.8998 ± 11.2722	0.6300 ± 0.0023	0.0819 ± 0.0017
swag wd	8740.6667 ± 35.5653	3931.9900 ± 45.6605	0.4130 ± 0.0066	0.0275 ± 0.0015
swag	8061.0000 ± 11.4310	5903.2605 ± 62.8167	0.5308 ± 0.0135	0.0163 ± 0.0019
baseline	8273.3333 ± 26.7872	8050.4467 ± 109.9864	0.2250 ± 0.0005	0.0809 ± 0.0020
vsgd	8255.6667 ± 24.1155	8919.8062 ± 106.3571	0.1761 ± 0.0078	0.0905 ± 0.0020

. As usual, we compare the results against the baseline of the individual network resulting from the pre-training phase. Results are obtained averaging over three independent seeds. Notice, as expanded in Appendix C that for swag we do consider two variants: the Bayesian correct one (swag) and a second variant that has better performance (swag wd). We stress again, as highlighted in Section 4.1 that not always goodness of approximation of the posterior and performance correlate positively. Additionally in this case, we found i-sgd to be competitive with other methods and superior to the baseline. Among the competitors, we found i-sgd to the easiest to tune, given the feature of a fixed learning rate informed by theoretical considerations; we believe that this is an important aspect to consider for a wide adoption of our proposal by practitioners.

5. Conclusions

sg methods allowed Bayesian posterior sampling algorithms, such as mcmc, to regain relevance in an age when datasets have reached extremely large sizes. However, despite mathematical elegance and promising results, current approaches from the literature are restricted to simple models. Indeed, the sampling properties of these algorithms are determined by simplifying assumptions on the loss landscape, which do not hold for the kind of complex models which are popular these days, such as deep models. Meanwhile, sg-mcmc algorithms require vanishing learning rates, which force practitioners to develop creative annealing schedules that are often model specific and difficult to justify.

We have attempted to target these weaknesses by suggesting a simpler algorithm that relies on fewer parameters and less strict assumptions compared to the literature on sg-mcmc. We used a unified mathematical notation to deepen our understanding of the role of the covariance of the noise of stochastic gradients and learning rate on the behavior of sg-mcmc algorithms. We then presented a practical variant of the sgd algorithm, which uses a constant learning rate, and an additional noise to perform Bayesian posterior sampling. Our proposal is derived from the ideal method, in which it is guaranteed that samples are generated from the true posterior. When the learning rate and noise terms are empirically estimated, with no user intervention, our method offers a very good approximation to the posterior, as demonstrated by the extensive experimental campaign.

We verified empirically the quality of our approach, and compared its performance to state-of-the-art sg-mcmc and alternative methods. Results, which span a variety of settings, indicated that our method is competitive to the alternatives from the state-of-the-art, while being much simpler to use.