Cyclical Stochastic Gradient MCMC for Bayesian Deep Learning

Introduction

Deep neural networks are often trained with stochastic optimization methods such as stochastic gradient decent (SGD) and its variants. Bayesian methods provide a principled alternative, accounting for model uncertainty in weight space (MacKay, 1992; Neal, 1996; Wilson, 2020), and achieve an automatic balance between model complexity and data fitting. Indeed, Bayesian methods have been shown to improve the generalization performance of DNNs (Hernández-Lobato & Adams, 2015; Blundell et al., 2015; Li et al., 2016a; Maddox et al., 2019; Wilson & Izmailov, 2020), while providing a principled representation of uncertainty on predictions, which is crucial for decision making.

Approximate inference for Bayesian deep learning has typically focused on deterministic approaches, such as variational methods (Hernández-Lobato & Adams, 2015; Blundell et al., 2015). By contrast, MCMC methods are now essentially unused for inference with modern deep neural networks, despite previously providing the gold standard of performance with smaller neural networks (Neal, 1996). Stochastic gradient Markov Chain Monte Carlo (SG-MCMC) methods (Welling & Teh, 2011; Chen et al., 2014; Ding et al., 2014; Li et al., 2016a) provide a promising direction for a sampling based approach to inference in Bayesian deep learning. Indeed, it has been shown that stochastic methods, which use mini-batches of data, are crucial for finding weight parameters that provide good generalization in modern deep neural networks (Keskar et al., 2016).

However, SG-MCMC algorithms for inference with modern neural networks face several challenges: $(\textup{\it i})$ In theory, SG-MCMC asymptotically converges to target distributions via a decreasing stepsize scheme, but suffers from a bounded estimation error in limited time (Teh et al., 2016; Chen et al., 2015). $(\textup{\it ii})$ In practice, empirical successes have been reported by training DNNs in relatively short time (Li et al., 2016b; Chen et al., 2014; Gan et al., 2016; Neelakantan et al., 2016; Saatchi & Wilson, 2017). For example, Saatchi & Wilson (2017) apply SG-MCMC to generative adversarial networks (GANs) to solve the mode collapse problem and capture diverse generation styles. However, the loss surface for DNNs is highly multimodal (Auer et al., 1996; Choromanska et al., 2015). In order for MCMC to be effective for posterior inference in modern neural networks, a crucial question remains: how do we make SG-MCMC efficiently explore a highly multimodal parameter space given a practical computational budget?

Several attempts have been made to improve the sampling efficiency of SG-MCMC. Stochastic Gradient Hamiltonian Monte Carlo (SGHMC) (Chen et al., 2014) introduces momentum to Langevin dynamics. Preconditioned stochastic gradient Langevin dynamics (pSGLD) (Li et al., 2016a) adaptively adjusts the sampler’s step size according to the local geometry of parameter space. Though simple and promising, these methods are still inefficient at exploring multimodal distributions in practice. It is our contention that this limitation arises from difficulties escaping local modes when using the small stepsizes that SG-MCMC methods typically require. Note that the stepsize in SG-MCMC controls the sampler’s behavior in two ways: the magnitude to deterministically drift towards high density regions wrt. the current stochastic gradient, and the level of injecting noise to randomly explore the parameter space. Therefore, a small stepsize reduces both abilities, resulting in a large numbers of iterations for the sampler to move across the modes.

In this paper, we propose to replace the traditional decreasing stepsize schedule in SG-MCMC with a cyclical variant. To note the distinction from traditional SG-MCMC, we refer to this method as Cyclical SG-MCMC (cSG-MCMC). The comparison is illustrated in Figure 1. The blue curve is the traditional decay, while the red curve shows the proposed cyclical schedule. Cyclical SG-MCMC operates in two stages: $(\textup{\it i})$ Exploration: when the stepsize is large (dashed red curves), we consider this stage as an effective burn-in mechanism, encouraging the sampler to take large moves and leave the local mode using the stochastic gradient. $(\textup{\it ii})$ Sampling: when the stepsize is small (solid red curves), the sampler explores one local mode. We collect samples for local distribution estimation during this stage. Further, we propose two practical techniques to improve estimation efficiency: (1) a system temperature for exploration and exploitation; (2) A weighted combination scheme for samples collected in different cycles to reflect their relative importance.

This procedure can be viewed as SG-MCMC with warm restarts: the exploration stage provides the warm restarts for its following sampling stage. cSG-MCMC combines the advantages from (1) the traditional SG-MCMC to characterize the fine-scale local density of a distribution and (2) the cyclical schedule in optimization to efficiently explore multimodal posterior distributions of the parameter space. In limited time, cSG-MCMC is a practical tool to provide significantly better mixing than the traditional SG-MCMC for complex distributions. cSG-MCMC can also be considered as an efficient approximation to parallel MCMC; cSG-MCMC can achieve similar performance to parallel MCMC with only a fraction of cost (reciprocal to the number of chains) that parallel MCMC requires.

To support our proposal, we also prove the non-asymptotic convergence for the cyclical schedule. We note that this is the first convergence analysis of a cyclical stepsize algorithm (including work in optimization). Moreover, we provide extensive experimental results to demonstrate the advantages of cSG-MCMC in sampling from multimodal distributions, including Bayesian neural networks and uncertainty estimation on several large and challenging datasets such as ImageNet.

In short, cSG-MCMC provides a simple and automatic approach to inference in modern Bayesian deep learning, with promising results, and theoretical support. This work is a step towards enabling MCMC approaches in Bayesian deep learning. We release code at https://github.com/ruqizhang/csgmcmc.

Preliminaries: SG-MCMC with a Decreasing Stepsize Schedule

SG-MCMC is a family of scalable sampling methods that enables inference with mini-batches of data. For a dataset $\mathcal{D}=\{d_{i}\}_{i=1}^{N}$ and a $\theta$-parameterized model, we have the likelihood $p(\mathcal{D}|\theta)$ and prior $p(\theta)$. The posterior distribution is $p(\theta|\mathcal{D})\propto\exp(-U(\theta))~{},$ where $U(\theta)$ is the potential energy given by $U(\theta)=-\log p(\mathcal{D}|\theta)-\log p(\theta)~{}.$

To guarantee asymptotic consistency with the true distribution, SG-MCMC requires that the step sizes satisfy the following assumption:

The step sizes $\{\alpha_{k}\}$ are decreasing, i.e., $0<\alpha_{k+1}<\alpha_{k}$, with 1) $\sum_{k=1}^{\infty}\alpha_{k}=\infty$; and 2) $\sum_{k=1}^{\infty}\alpha_{k}^{2}<\infty$.

Without a decreasing step-size, the estimation error from numerical approximations is asymptotically biased. One typical decaying step-size schedule is $\alpha_{k}=a(b+k)^{-\gamma}$, with $\gamma\in(0.5,1]$ and $(a,b)$ some positive constants (Welling & Teh, 2011).

Cyclical SG-MCMC

We now introduce our cyclical SG-MCMC (cSG-MCMC) algorithm. cSG-MCMC consists of two stages: exploration and sampling. In the following, we first introduce the cyclical step-size schedule, and then describe the exploration stage in Section 3.1 and the sampling stage in Section 3.2. We propose an approach to combining samples for testing in Section F.

Assumption 1 guarantees the consistency of our estimation with the true distribution in the asymptotic time. The approximation error in limited time is characterized as the risk of an estimator $R=B^{2}+V$, where $B$ is the bias and $V$ is the variance. In the case of infinite computation time, the traditional SG-MCMC setting can reduce the bias and variance to zero. However, the time budget is often limited in practice, and there is always a trade-off between bias and variance. We therefore decrease the overall approximation error $R$ by reducing the variance through obtaining more effective samples. The effective sample size can be increased if fewer correlated samples from different distribution modes are collected.

For deep neural networks, the parameter space is highly multimodal. In practice, SG-MCMC with the traditional decreasing stepsize schedule becomes trapped in a local mode, though injecting noise may help the sampler to escape in the asymptotic regime (Zhang et al., 2017). Inspired to improve the exploration of the multimodal posteriors for deep neural networks, with a simple and automatic approach, we propose the cyclical cosine stepsize schedule for SG-MCMC. The stepsize at iteration $k$ is defined as:

where $\alpha_{0}$ is the initial stepsize, $M$ is the number of cycles and $K$ is the number of total iterations (Loshchilov & Hutter, 2016; Huang et al., 2017).

The stepsize $\alpha_{k}$ varies periodically with $k$. In each period, $\alpha_{k}$ starts at $\alpha_{0}$, and gradually decreases to . Within one period, SG-MCMC starts with a large stepsize, resulting in aggressive exploration in the parameter space; as the stepsize is decreasing, SG-MCMC explores local regions. In the next period, the Markov chain restarts with a large stepsize, encouraging the sampler to escape from the current mode and explore a new area of the posterior.

In optimization, the cyclical cosine annealing stepsize schedule has been demonstrated to be able to find diverse solutions in multimodal objectives, though not specifically different modes, using stochastic gradient methods (Loshchilov & Hutter, 2016; Huang et al., 2017; Garipov et al., 2018; Fu et al., 2019). Alternatively, we adopt the technique to SG-MCMC as an effective scheme for sampling from multimodal distributions.

1 Exploration

The first stage of cyclical SG-MCMC, exploration, discovers parameters near local modes of an objective function. Unfortunately, it is undesirable to directly apply the cyclical schedule in optimization to SG-MCMC for collecting samples at every step. SG-MCMC often requires a small stepsize in order to control the error induced by the noise from using a minibatch approximation. If the stepsize is too large, the stationary distribution of SG-MCMC might be far away from the true posterior distribution. To correct this error, it is possible to do stochastic Metropolis-Hastings (MH) (Korattikara et al., 2014; Bardenet et al., 2014; Chen et al., 2016b). However, stochastic MH correction is still computationally too expensive. Further, it is easy to get rejected with an aggressive large stepsize, and every rejection is a waste of gradient computations.

To alleviate this problem, we propose to introduce a system temperature $T$ to control the sampler’s behaviour: $p(\theta|\mathcal{D})\propto\exp(-U(\theta)/T)$. Note that the setting $T=1$ corresponds to sampling from the untempered posterior. When $T\rightarrow 0$, the posterior distribution becomes a point mass. Sampling from $\lim_{T\rightarrow 0}\exp(-U(\theta)/T)$ is equivalent to minimizing $U(\theta)$; in this context, SG-MCMC methods become stochastic gradient optimization methods.

One may increase the temperature $T$ from 0 to 1 when the step-size is decreasing. We simply consider $T=0$ and perform optimization as the burn-in stage, when the completed proportion of a cycle $r(k)=\frac{\mod(k-1,\lceil K/M\rceil)}{\lceil K/M\rceil}$ is smaller than a given threshold: $r(k)<\beta$. Note that $\beta\in(0,1)$ balances the proportion of the exploration and sampling stages in cSG-MCMC.

2 Sampling

The sampling stage corresponds to $T=1$ of the exploration stage. When $r(k)>\beta$ or step-sizes are sufficiently small, we initiate SG-MCMC updates and collect samples until this cycle ends.

One may consider the exploration stage as automatically providing warm restarts for the sampling stage. Exploration alleviates the inefficient mixing and inability to traverse the multimodal distributions of the traditional SG-MCMC methods. SG-MCMC with warm restarts explores different parts of the posterior distribution and captures multiple modes in a single training procedure.

In summary, the proposed cyclical SG-MCMC repeats the two stages, with three key advantages: $(\textup{\it i})$ It restarts with a large stepsize at the beginning of a cycle which provides enough perturbation and encourages the model to escape from the current mode. $(\textup{\it ii})$ The stepsize decreases more quickly inside one cycle than a traditional schedule, making the sampler better characterize the density of the local regions. $(\textup{\it iii})$ This cyclical stepsize shares the advantage of the “super-convergence” property discussed in Smith & Topin (2017): cSG-MCMC can accelerate convergence for DNNs by up to an order of magnitude.

Connection to the Santa algorithm.

It is interesting to note that our approach inverts steps of the Santa algorithm (Chen et al., 2016a) for optimization. Santa is a simulated-annealing-based optimization algorithm with an exploration stage when $T=1$, then gradually anneals $T\rightarrow 0$ in a refinement stage for global optimization. In contrast, our goal is to draw samples for multimodal distributions, thus we explore with $T=0$ and sample with $T=1$. Another fundamental difference is that Santa adopts the traditional stepsize decay, while we use the cyclical schedule.

We visually compare the difference between cyclical and traditional step size schedules (described in Section 2) in Figure 1. The cyclical SG-MCMC algorithm is presented in Algorithm 1.

Connection to Parallel MCMC.

Running parallel Markov chains is a natural and effective way to draw samples from multimodal distributions (VanDerwerken & Schmidler, 2013; Ahn et al., 2014). However, the training cost increases linearly with the number of chains. Cyclical SG-MCMC can be seen as an efficient way to approximate parallel MCMC. Each cycle effectively estimates a different region of posterior. Note cyclical SG-MCMC runs along a single training pass. Therefore, its computational cost is the same as single chain SG-MCMC while significantly less than parallel MCMC.

Combining Samples.

In cyclical SG-MCMC, we obtain samples from multiple modes of a posterior distribution by running the cyclical step size schedule for many periods. We provide a sampling combination scheme to effectively use the collected samples in Section F in the appendix.

Theoretical Analysis

Under Assumptions 2 in the appendix, for a smooth test function $\phi$, the bias and MSE of cSGLD are bounded as:

Convergence under the Wasserstein distance

Next, we consider the more general case of SGLD and characterize convergence rates in terms of a stronger metric of 2-Wasserstein distance, defined as:

where $\Gamma(\mu,\nu)$ is the set of joint distributions over $(\theta,\theta^{\prime})$ such that the two marginals equal $\mu$ and $\nu$, respectively.

Denote the distribution of $\theta_{t}$ in the SDE as $\nu_{t}$. According to Chiang & Hwang (1987), the stationary distribution $\nu_{\infty}$ matches our target distribution. Let $\mu_{K}$ be the distribution of the sample from our proposed cSGLD algorithm at the $K$-th iteration. Our goal is to derive a convergence bound on $W_{2}(\mu_{K},\nu_{\infty})$. We adopt standard assumptions as in most existing work, which are detailed in Assumption 3 in the appendix. Theorem 2 summarizes our main theoretical result.

Under Assumption 3 in the appendix, there exist constants $(C_{0},C_{1},C_{2},C_{3})$ independent of the stepsizes such that the convergence rate of our proposed cSGLD with cyclical stepsize sequence equation 1 is bounded for all $K$ satisfying ($K$ mod $M$ =0), as $W_{2}(\mu_{K},\nu_{\infty})\leq$

Particularly, if we further assume $\alpha_{0}=O(K^{-\beta})$ for $\forall\beta>1$, $W_{2}(\mu_{K},\nu_{\infty})\leq C_{3}+\left(6+\frac{C_{2}}{K^{\beta-1}}\right)^{\frac{1}{2}}[(\frac{2C_{1}}{K^{2\beta-1}}+\frac{2C_{0}}{K^{\beta-1}})^{\frac{1}{2}}+(\frac{C_{1}}{K^{2\beta-1}}+\frac{C_{0}}{K^{\beta-1}})^{\frac{1}{4}}]$.

$\textup{\it i})$ The bound is decomposed into two parts: the first part measures convergence speed of exact solution to the stationary distribution, i.e., $\nu_{\sum_{k}\alpha_{k}}$ to $\nu_{\infty}$; the second part measures the numerical error, i.e., between $\mu_{K}$ and $\nu_{\sum_{k}\alpha_{k}}$. $\textup{\it ii})$ The overall bound offers a same order of dependency on $K$ as in standard SGLD (please see the bound for SGLD in Section E of the appendix. See also Raginsky et al. (2017)). $\textup{\it iii})$ If one imposes stricter assumptions such as in the convex case, the bound can be further improved. Specific bounds are derived in the appendix. We did not consider this case due to the discrepancy from real applications.

Experiments

We demonstrate cSG-MCMC on several tasks, including a synthetic multimodal distribution (Section 5.1), image classification on Bayesian neural networks (Section 5.2) and uncertainty estimation in Section 5.3. We also demonstrate cSG-MCMC can improve the estimate efficiency for uni-modal distributions using Bayesian logistic regression in Section A.2 in the appendix. We choose SLGD and SGHMC as the representative baseline algorithms. Their cyclical counterpart are called cSGLD and cSGHMC, respectively.

We first demonstrate the ability of cSG-MCMC for sampling from a multi-modal distribution on a 2D mixture of 25 Gaussians. Specifically, we compare cSGLD with SGLD in two setting: (1) parallel running with 4 chains and (2) running with a single chain, respectively. Each chain runs for 50k iterations. The stepsize schedule of SGLD is $\alpha_{k}\propto 0.05k^{-0.55}$. In cSGLD, we set $M=30$ and the initial stepsize $\alpha_{0}=0.09$. The proportion of exploration stage $\beta=\frac{1}{4}$. Fig 2 shows the estimated density using sampling results for SGLD and cSGLD in the parallel setting. We observed that SGLD gets trapped in the local modes, depending on the initial position. In any practical time period, SGLD could only characterize partial distribution. In contrast, cSGLD is able to find and characterize all modes, regardless of the initial position. cSGLD leverages large step sizes to discover a new mode, and small step sizes to explore local modes. This result suggests cSGLD can be a significantly favourable choice in the non-asymptotic setting, for example only 50k iterations in this case. The single chain results and the quantitative results on mode coverage are reported in Section A.1 of the appendix.

2 Bayesian Neural Networks

We demonstrate the effectiveness of cSG-MCMC on Bayesian neural networks for classification on CIFAR-10 and CIFAR-100. We compare with $(\textup{\it i})$ traditional SG-MCMC; $(\textup{\it ii})$ traditional stochastic optimization methods, including stochastic gradient descent (SGD) and stochastic gradient descent with momentum (SGDM); and $(\textup{\it iii})$ Snapshot: a stochastic optimization ensemble method method with a the cyclical stepsize schedule (Huang et al., 2017). We use a ResNet-18 (He et al., 2016) and run all algorithms for 200 epochs. We report the test errors averaged over 3 runs, and the standard error ($\pm$) from the mean predictor.

We set $M=4$ and $\alpha_{0}=0.5$ for cSGLD, cSGHMC and Snapshot. The proportion hyper-parameter $\beta=$0.8 and 0.94 for CIFAR-10 and CIFAR-100, respectively. We collect 3 samples per cycle. In practice, we found that the collected samples share similarly high likelihood for DNNs, thus one may simply set the normalizing term $w_{i}$ in equation 33 to be the same for faster testing.

We found that tempering helps improve performance for Bayesian inference with neural networks. Tempering for SG-MCMC was first used by Li et al. (2016a) as a practical technique for neural network training for fast convergence in limited timehttps://github.com/ChunyuanLI/pSGLD/issues/2. We simply use the prescribed temperature of Li et al. (2016a) without tuning, but better results of the sampling methods can be achieved by tuning the temperature. More details are in Appendix J. We hypothesize that tempering helps due to the overparametrization of neural networks. Tempering enables one to leverage the inductive biases of the network, while representing the belief that the model capacity can be misspecified. In work on Safe Bayes, also known as generalized and fractional Bayesian inference, tempered posteriors are well-known to help under misspecification (e.g., Barron & Cover, 1991; de Heide et al., 2019; Grünwald et al., 2017).

For the traditional SG-MCMC methods, we found that noise injection early in training hurts convergence. To make these baselines as competitive as possible, we thus avoid noise injection for the first 150 epochs of training (corresponding to the zero temperature limit of SGLD and SGHMC), and resume SGMCMC as usual (with noise) for the last 50 epochs. This scheme is similar to the exploration and sampling stages within one cycle of cSG-MCMC. We collect 20 samples for the MCMC methods and average their predictions in testing.

We report the testing errors in Table LABEL:tab:CIFAR to compare with the non-parallel algorithms. Snapshot and traditional SG-MCMC reduce the testing errors on both datasets. Performance variance for these methods is also relatively small, due to the multiple networks in the Bayesian model average. Further, cSG-MCMC significantly outperforms Snapshot ensembles and the traditional SG-MCMC, demonstrating the importance of (1) capturing diverse modes compared to traditional SG-MCMC, and (2) capturing fine-scale characteristics of the distribution compared with Snapshot ensembles.

Diversity in Weight Space.

To further demonstrate our hypothesis that with a limited budget cSG-MCMC can find diverse modes, while traditional SG-MCMC cannot, we visualize the 12 samples we collect from cSG-MCMC and SG-MCMC on CIFAR-100 respectively using Multidimensional Scaling (MDS) in Figure 3 (a). MDS uses a Euclidean distance metric between the weight of samples. We see that the samples of cSG-MCMC form 4 clusters, which means they are from 4 different modes in weight space. However, all samples from SG-MCMC only form one cluster, which indicates traditional SG-MCMC gets trapped in one mode and only samples from that mode.

Diversity in Prediction.

To further demonstrate the samples from different cycles of cSG-MCMC provide diverse predictions we choose one sample from each cycle and linearly interpolate between two of them (Goodfellow et al., 2014; Huang et al., 2017). Specifically, let $J(\theta)$ be the test error of a sample with parameter $\theta$. We compute the test error of the convex combination of two samples $J(\lambda\theta_{1}+(1-\lambda)\theta_{2})$, where $\lambda\in$.

We linearly interpolate between two samples from neighboring chains of cSG-MCMC since they are the most likely to be similar. We randomly select 4 samples from SG-MCMC. If the samples are from the same mode, the test error of the linear interpolation of parameters will be relatively smooth, while if the samples are from different modes, the test error of the parameter interpolation will have a spike when $\lambda$ is between 0 and 1.

We show the results of interpolation for cSG-MCMC and SG-MCMC on CIFAR-100 in Figure 3 (b). We see a spike in the test error in each linear interpolation of parameters between two samples from neighboring chains in cSG-MCMC while the linear interpolation for samples of SG-MCMC is smooth. This result suggests that samples of cSG-MCMC from different chains are from different modes while samples of SG-MCMC are from the same mode.

Although the test error of a single sample of cSG-MCMC is worse than that of SG-MCMC shown in Figure 3 (c), the ensemble of these samples significantly improves the test error, indicating that samples from different modes provide different predictions and make mistakes on different data points. Thus these diverse samples can complement each other, resulting in a lower test error, and demonstrating the advantage of exploring diverse modes using cSG-MCMC.

cSG-MCMC can be viewed as an economical alternative to parallel MCMC. We verify how closely cSG-MCMC can approximate the performance of parallel MCMC, but with more convenience and less computational expense. We also note that we can improve parallel MCMC with the proposed cyclical stepsize schedule.

We report the testing errors in Table 5.2 to compare multiple-chain results. (1) Four chains used, each runs 200 epochs (800 epochs in total), the results are shown in the first 4 columns (Cyclical+Parallel vs Decreasing+Parallel). We see that cSG-MCMC variants provide lower errors than plain SG-MCMC. (2) We reduce the number of epochs ($\#$epoch) of parallel MCMC to 100 epoch each for decreasing stepsize schedule. The total cost is 400 epochs. We compare its performance with cyclical single chain (200 epochs in total) in the last 4 columns (Decreasing+Parallel vs Cyclical+Single). We see that the cyclical schedule running on a single chain performs best even with half the computational cost! All the results indicate the importance of warm re-starts using the proposed cyclical schedule. For a given total cost budget, the proposed cSGMCMC is preferable to parallel sampling.

Comparison to Snapshot Optimization.

We carefully compared with Snapshot, as our cSG-MCMC can be viewed as the sampling counterpart of the Snapshot optimization method. We plot the test error wrt.various number of cycles $M$ in Fig. 3. As $M$ increases, cSG-MCMC and Snapshot both improve. However, given a fixed $M$, cSG-MCMC yields substantially lower test errors than Snapshot. This result is due to the ability of cSG-MCMC to better characterize the local distribution of modes: Snapshot provides a singe minimum per cycle, while cSG-MCMC fully exploits the mode with more samples, which could provide weight uncertainty estimate and avoid over-fitting.

We further study different learning algorithms on a large-scale dataset, ImageNet. ResNet-50 is used as the architecture, and 120 epochs for each run. The results on the testing set are summarized in Table LABEL:tab:imagenet, including NLL, Top1 and Top5 accuracy ($\%$), respectively. 3 cycles are considered for both cSGHMC and Snapshot, and we collect 3 samples per cycle. We see that cSGHMC yields the lowest testing NLL, indicating that the cycle schedule is an effective technique to explore the parameter space, and diversified samples can help prevent over-fitting.

3 Uncertainty Evaluation

To demonstrate how predictive uncertainty benefits from exploring multiple modes in the posterior of neural network weights, we consider the task of uncertainty estimation for out-of-distribution samples (Lakshminarayanan et al., 2017). We train a three-layer MLP model on the standard MNIST train dataset until convergence using different algorithms, and estimate the entropy of the predictive distribution on the notMNIST dataset (Bulatov, 2011). Since the samples from the notMNIST dataset belong to the unseen classes, ideally the predictive distribution of the trained model should be uniform over the notMNIST digits, which gives the maximum entropy.

In Figure 4, we plot the empirical CDF for the entropy of the predictive distributions on notMNIST. We see that the uncertainty estimates from cSGHMC and cSGLD are better than the other methods, since the probability of a low entropy prediction is overall lower. cSG-MCMC algorithms explore more modes in the weight space, each mode characterizes a meaningfully different representation of MNIST data. When testing on the out-of-distribution dataset (notMNIST), each mode can provide different predictions over the label space, leading to more reasonable uncertainty estimates. Snapshot achieves less entropy than cSG-MCMC, since it represents each mode with a single point.

The traditional SG-MCMC methods also provide better uncertainty estimation compared to their optimization counterparts, because they characterize a local region of the parameter space, rather than a single point. cSG-MCMC can be regarded as a combination of these two worlds: a wide coverage of many modes in Snapshot, and fine-scale characterization of local regions in SG-MCMC.

Discussion

We have proposed cyclical SG-MCMC methods to automatically explore complex multimodal distributions. Our approach is particularly compelling for Bayesian deep learning, which involves rich multimodal parameter posteriors corresponding to meaningfully different representations. We have also shown that our cyclical methods explore unimodal distributions more efficiently. These results are in accordance with theory we developed to show that cyclical SG-MCMC will converge faster to samples from a stationary distribution in general settings. Moreover, we show cyclical SG-MCMC methods provide more accurate uncertainty estimation, by capturing more diversity in the hypothesis space corresponding to settings of model parameters.

While MCMC was once the gold standard for inference with neural networks, it is now rarely used in modern deep learning. We hope that this paper will help renew interest in MCMC for posterior inference in deep learning. Indeed, MCMC is uniquely positioned to explore the rich multimodal posterior distributions of modern neural networks, which can lead to improved accuracy, reliability, and uncertainty representation.

AGW was supported by an Amazon Research Award, Facebook Research, NSF I-DISRE 193471, NIH R01 DA048764-01A1, NSF IIS-1563887, and NSF IIS-1910266.

References

Appendix A Experimental Results

where $\lambda=\frac{1}{25}$, $\mu=\{-4,-2,0,2,4\}^{\top}\times\{-4,-2,0,2,4\}$, $\Sigma=\begin{bmatrix}0.03&0\\ 0&0.03\end{bmatrix}$.

In Figure 5, we show the estimated density for SGLD and cSGLD in the non-parallel setting.

To quantitatively show the ability of different algorithms to explore multi-modal distributions, we define the mode-coverage metric: when the number of samples falling within the radius $r$ of a mode center is larger than a threshold $\bar{n}$, we consider this mode covered. On this dataset, we choose $r=0.25$ and $\bar{n}=100$. Table 4 shows the mode-coverage for several algorithms, based on 10 different runs.

A.2 Bayesian Logistic Regression

We consider Bayesian logistic regression (BLR) on three real-world datasets from the UCI repository: Australian (15 covariates, 690 data points), German (25 covariates, 1000 data points) and Heart (14 covariates, 270 data points). For all experiments, we collect 5000 samples with 5000 burn-in iterations. Following the settings in Li et al. (2016a), we report median effective sample size (ESS) in Table 5.

Note that BLR is unimodal in parameter space. We use this experiment as an adversarial situation for cSG-MCMC, which we primarily designed to explore multiple modes. We note that even in the unimodal setting, cSG-MCMC more effectively explores the parameter space than popular alternatives. We can also use these experiments to understand how samplers respond to varying parameter dimensionality and training set sizes.

Overall, cSG-MCMC dramatically outperforms SG-MCMC, which demonstrates the fast mixing rate due to the warm restarts. On the small dataset Heart, SGHMC and cSGHMC achieve the same results, because the posterior of BLR on this dataset is simple. However, in higher dimensional spaces (e.g., Australian and German), cSG-MCMC shows significantly higher ESS; this result means that each cycle in cSG-MCMC can characterize a different region of the posteriors, combining multiple cycles yields more accurate overall approximation.

Appendix B Assumptions

In the analysis, we define a functional $\psi$ that solves the following Poisson Equation:

The solution functional $\psi(\theta_{k})$ characterizes the difference between $\phi(\theta_{k})$ and the posterior average $\bar{\phi}$ for every $\theta_{k}$, thus would typically possess a unique solution, which is at least as smooth as $\phi$ under the elliptic or hypoelliptic settings (Mattingly et al., 2010). Following Chen et al. (2015); Vollmer et al. (2016), we make certain assumptions on the solution functional, $\psi$, of the Poisson equation equation 3.

B.2 Assumptions in convergence under the Wasserstein distance

Following existing work in Raginsky et al. (2017), we adopt the following standard assumptions summarized in Assumption 3.

There exists some constants $A\geq 0$ and $B\geq 0$, such that $U(0)\leq A$ and $\nabla U(0)\leq B$.

The function U is $L_{U}$-smooth : $\|\nabla U(w)-\nabla U(v)\|\leq L_{U}\|w-v\|$.

The function U is $(m_{u},b)-dissipative$, which means for some $m_{U}>0$ and $b>0$ $\langle w,\nabla U(w)\rangle$ $\geq$ $m_{U}\|w\|^{2}-b$.

We can choose $\mu_{0}$ which satisfies the requirement: $\kappa_{0}:=log\int e^{\|w\|^{2}}\mu_{0}(w)dw<\infty$.

Appendix C Proof of Theorem 1

To prove the theorem, we borrow tools developed by Chen et al. (2015); Vollmer et al. (2015). We first rephrase the stepsize assumptions in general SG-MCMC in Assumption 4.

The algorithm adopts an $N$-th order integrator. The step sizes $\{h_{k}\}$ are such that $0<h_{k+1}<h_{k}$, and satisfy 1) $\sum_{k=1}^{\infty}h_{k}=\infty$; and 2) $\lim_{K\rightarrow\infty}\frac{\sum_{k=1}^{K}h_{k}^{N+1}}{\sum_{k=1}^{K}h_{k}}=0$.

Our prove can be derived by the following results from Chen et al. (2015).

Let $S_{K}\triangleq\sum_{k=1}^{K}h_{k}$. Under Assumptions 2 and 4, for a smooth test function $\phi$, the bias and MSE of a decreasing-step-size SG-MCMC with a $N$th-order integrator at time $S_{L}$ are bounded as:

Note that Assumption 4 is only required if one wants to prove the asymptotically unbias of an algorithm. Lemma 1 still applies even if Assumption 4 is not satisfied. In this case one would obtain a biased algorithm, which is the case of cSGLD.

Our results is actually a special case of Lemma 1. To see that, first note that our cSGLD adopts a first order integrator, thus $N=1$. To proceed, note that $S_{K}=\sum_{k=1}^{K}\alpha_{k}=O(\alpha_{0}K)$, and

Let $\nu_{t}$ denote the distribution of $\theta_{t}$, and the stationary distribution of equation 34 be $p(\theta|\mathcal{D})$, which means $\nu_{\infty}=p(\theta|\mathcal{D})$.

Further, let $\mu_{k}$ denote the distribution of $\theta_{k}$.

, we need to give the bounds for these two parts respectively.

Then we focus on the following continuous-time interpolation of $\theta_{k}$:

Let $\operatorname{\mathbf{P}}_{V}^{t}:=\mathcal{L}\left(V(s):0\leq s\leq t\right)$ and $\operatorname{\mathbf{P}}_{\theta}^{t}:=\mathcal{L}\left(\theta(s):0\leq s\leq t\right)$ and according to the proof of Lemma 3.6 in Raginsky et al. (2017), we can derive a similar result for the relative entropy of $\operatorname{\mathbf{P}}_{V}^{t}$ and $\operatorname{\mathbf{P}}_{\theta}^{t}$:

The last line follows the fact that $\mathcal{L}(\underline{\theta}(s))=\mathcal{L}(V(s)),~{}\forall s$. Then we will let $t=\sum_{k=1}^{K}\alpha_{k}$ and we can use the martingale property of the integral to derive:

For the first part (14), we consider some $s\in[{\sum_{k=1}^{j}\alpha_{k}},{\sum_{k=1}^{j+1}\alpha_{k}})$, for which the following holds:

Thus, we can use Lemma 3.1 and 3.2 in Raginsky et al. (2017) for the following result:

Hence we can bound the first part, (choosing $\alpha_{0}\leq 1$),

The last line (17) follows fromNote: we only focus on the case when $K\mod M=0$. equation C. The second part (15) can be bounded as follows:

Due to the data-processing inequality for the relative entropy, we have

Similar to the statement of Lemma 3.2 in Raginsky et al. (2017), we can fix $\alpha_{0}\in(0,1\wedge\frac{m_{U}}{4M_{U}^{2}})$. Then, we can know that

, where $\alpha_{min}$ is defined as $\alpha_{min}\triangleq\frac{\alpha_{0}}{2}\left[\cos\left(\frac{\pi~{}\text{mod}(\lceil K/M\rceil-1,\lceil K/M\rceil)}{\lceil K/M\rceil}\right)+1\right]$.

If $1-2\alpha_{min}m_{U}+4\alpha_{min}^{2}M_{U}^{2}\leq 0$, then from equation 18 it follows that

If $0\leq 1-2\alpha_{min}m_{U}+4\alpha_{min}^{2}M_{U}^{2}\leq 1$, then iterating equation 18 gives

Due to the expression of $\frac{\alpha_{0}}{\alpha_{min}}$, $C_{0}$ is independent of $\alpha_{0}$. Then we denote the $6L_{U}^{2}(C_{0}+d)$ as $C_{1}$ and we can derive

Then according to Proposition 3.1 in Bolley & Villani (2005) and Lemma 3.3 in Raginsky et al. (2017), if we denote $\kappa_{0}+2b+2d$ as $C_{2}$, we can derive the following result:

We can directly get the following results from (3.17) in Raginsky et al. (2017) that there exist some positive constants $(C_{3},C_{4})$,

Now combining the bounds for $W_{2}(\mu_{K},\nu_{\sum_{k=1}^{K}\alpha_{k}})$ and $W_{2}(\nu_{\sum_{k=1}^{K}\alpha_{k}},\nu_{\infty})$, substituting $\alpha_{0}=O(1/K^{\beta})$, and noting $W_{2}(\nu_{\sum_{k=1}^{K}\alpha_{k}},\nu_{\infty})$ decreases w.r.t.​ $K$, we arrive at the bound stated in the theorem.

Appendix E Relation with SGLD

For the standard polynomially-decay-stepsize SGLD, the convergence rate is bounded as

and $W_{2}(\nu_{\sum_{k=1}^{K}h_{k}},\nu_{\infty})\leq C_{3}\exp(-\frac{\sum_{k=1}^{K}h_{k}}{C_{4}})$.

Similar to the proof of $W_{2}({\mu}_{K},\nu_{\infty})$ in cSGLD, we get the following update rule for SGLD with the stepsize following a polynomial decay i.e., $h_{k}=\frac{h_{0}}{k}$,

, we need to give the bounds for these two parts respectively.

Then we focus on the following continuous-time interpolation of $\theta_{k}$:

Let $\operatorname{\mathbf{P}}_{V}^{t}:=\mathcal{L}\left(V(s):0\leq s\leq t\right)$ and $\operatorname{\mathbf{P}}_{\theta}^{t}:=\mathcal{L}\left(\theta(s):0\leq s\leq t\right)$ and according to the proof of Lemma 3.6 in Raginsky et al. (2017), we can derive the similar result for the relative entropy of $\operatorname{\mathbf{P}}_{V}^{t}$ and $\operatorname{\mathbf{P}}_{\theta}^{t}$:

The last line follows the fact that $\mathcal{L}(\underline{\theta}(s))=\mathcal{L}(V(s)),~{}\forall s$. Then we will let $t=\sum_{k=1}^{K}h_{k}$ and we can use the martingale property of integral to derive:

For the first part (29), we consider some $s\in[{\sum_{k=1}^{j}h_{k}},{\sum_{k=1}^{j+1}h_{k}})$, the following equation holds:

Thus, we can use Lemma 3.1 and 3.2 in Raginsky et al. (2017) for the following result:

Hence we can bound the first part, (choosing $h_{0}\leq 1$),

where the last line follows from the fact that

The second part (30) can be bounded as follows:

Due to the data-processing inequality for the relative entropy, we have

Then we denote the $6L_{U}^{2}(D_{0}+d)$ as $D_{1}$ and we can derive

Then according to Proposition 3.1 in Bolley & Villani (2005) and Lemma 3.3 in Raginsky et al. (2017), if we denote $\kappa_{0}+2b+2d$ as $D_{2}$, we can derive the following result,

We can directly get the following results from (3.17) in Raginsky et al. (2017) that there exist some positive constants $(C_{3},C_{4})$,

∎ Based on the convergence error bounds, we discuss an informal comparison with standard SGLD. Consider the following two cases.We must emphasize that since the term $W_{2}(\mu_{K},\nu_{\sum_{k=1}^{K}\alpha_{k}})$ in the equation 9 increases w.r.t. $K$, our $\alpha_{0}$ must be set small enough in practice. Hence, in this informal comparison, we also set $\alpha_{0}$ small enough to make $W_{2}(\mu_{K},\nu_{\sum_{k=1}^{K}\alpha_{k}})$ less important.

Appendix F Combining Samples

In cyclical SG-MCMC, we obtain samples from multiple modes of a posterior distribution by running the cyclical step size schedule for many periods. We now show how to effectively utilize the collected samples. We consider each cycle exploring different part of the target distribution $p(\theta|\mathcal{D})$ on a metric space $\Theta$. As we have $M$ cycles in total, the $m$th cycle characterizes a local region $\Theta_{m}\subset\Theta$, defining the “sub-posterior” distribution: $p_{m}(\theta|\mathcal{D})=\frac{p(\theta|\mathcal{D})\textbf{1}_{\Theta_{m}}}{w_{m}},\text{with}\ w_{m}=\int_{\Theta_{m}}p(\theta|\mathcal{D})d\theta,$ where $w_{m}$ is a normalizing constant. For a testing function $f(\theta)$, we are often interested in its true posterior expectation $\bar{f}=\int f(\theta)p(\theta|\mathcal{D})d\theta$. The sample-based estimation is

where $K_{m}$ is the number of samples from the $m$th cycle, and $\theta^{(m)}\in\Theta_{m}$.

The weight for each cycle $w_{i}$ is estimated using the harmonic mean method (Green, 1995; Raftery et al., 2006): $\hat{w}_{m}\approx[\frac{1}{K_{m}}\sum_{j=1}^{K_{m}}\frac{1}{p(\mathcal{D}|\theta_{j}^{(m)})}]^{-1}.$ This approach provides a simple and consistent estimator, where the only additional cost is to traverse the training dataset to evaluate the likelihood $p(\mathcal{D}|\theta_{j}^{(m)})$ for each sample $\theta_{j}^{(m)}$. We evaluate the likelihood once off-line and store the result for testing.

Let $\nu_{t}$ denote the distribution of $\theta_{t}$, and the stationary distribution of equation 34 be $p(\theta|\mathcal{D})$, which means $\nu_{\infty}=p(\theta|\mathcal{D})$.

However, the exact evaluation of the gradient $\nabla U$ is computationally expensive. Hence, we need to adopt noisy evaluations of $\nabla U$. For simplicity, we assume that at any point $\theta_{k}$, we can observe the value

where ${\zeta_{k}:k=0,1,2,...}$ is a sequence of random (noise) vectors. Then the algorithm is defined as:

Further, let $\mu_{k}$ denote the distribution of $\theta_{k}$.

Following the existing work in Dalalyan & Karagulyan (2019), we adopt the following standard assumptions summarized in Assumption 5,

For some positive constants m and M, it holds

(independence of updates) $\xi_{k+1}$ in equation 35 is independent of $(\zeta_{1},\zeta_{2},...,\zeta_{k})$

Under Assumption 5 in the appendix and $\alpha_{0}\in(0,\frac{1}{m}\wedge\frac{2}{M})$, if we define the $\alpha_{min}$ as $\frac{\alpha_{0}}{2}\left[\cos\left(\frac{\pi~{}\text{mod}(\lceil K/M\rceil-1,\lceil K/M\rceil)}{\lceil K/M\rceil}\right)+1\right]$, we can derive the the following bounds.

If $m\alpha_{min}+M\alpha_{0}\leq 2$, then $W_{2}(\mu_{k+1},\nu_{\infty})\leq$

If $m\alpha_{min}+M\alpha_{0}>2$, then $W_{2}(\mu_{k+1},\nu_{\infty})\leq$

where the $M,m,\delta,\sigma$ are some positive constants defined in Assumption 5

G.2 Proof

According to the equation 1, we can find that the stepsize $\alpha_{k}$ varies from $\alpha_{0}$ to $\alpha_{min}$, where $\alpha_{min}$ is defined as $\alpha_{min}\triangleq\frac{\alpha_{0}}{2}\left[\cos\left(\frac{\pi~{}\text{mod}(\lceil K/M\rceil-1,\lceil K/M\rceil)}{\lceil K/M\rceil}\right)+1\right]$. When $0<\alpha_{0}<\min(2/M,1/m)$, it is easy for us to know that $0<\alpha_{k}<\min(2/M,1/m)$ for every $k>0$. Then we can derive that all the $\rho_{k}\triangleq\max(1-m\alpha_{k},M\alpha_{k}-1)$ will satisfy $0<\rho_{k}<1$. Now according to the Proposition 2 in Dalalyan & Karagulyan (2019), we can derive the result that

Then we will use another lemma derived from Dalalyan & Karagulyan (2019).

If A,B,C are non-negative numbers such that A $\in$ (0,1) and the sequence of non-negative numbers ${y_{k}}$ satisfies the following inequality

Using Lemma 2, we can finish our proof now.

If $m\alpha_{min}+M\alpha_{0}\leq 2$, the $\rho_{k}$ will satisfy $\rho_{k}\leq 1-m\alpha_{min}$ for every $k>0$. Then the equation 38 will turn into

for every $k>0$. Then we can set $A=m\alpha_{min}$, $C=1.65M(\alpha_{0}^{3}d)^{1/2}+\alpha_{0}\delta d^{1/2}$, $B=\delta\alpha_{0}d^{1/2}$ and we can get the result.

If $m\alpha_{min}+M\alpha_{0}>2$, the $\rho_{k}$ will satisfy $\rho_{k}\leq M\alpha_{0}-1$ for every $k>0$. Then the equation 38 will turn into

for every $k>0$. Then we can set $A=2-M\alpha_{0}$, $C=1.65M(\alpha_{0}^{3}d)^{1/2}+\alpha_{0}\delta d^{1/2}$, $B=\delta\alpha_{0}d^{1/2}$ and we can get the result.

Appendix H Future Direction for the Wasserstein gradient flows

We would like to point out that the convergence theorems developed in the above several sections can be potentially applied to study the convergence of the Wasserstein gradient flows (Santambrogio, 2017), which can be regarded as a continuous-time MCMC (Chen et al., 2018; Liu et al., 2019). The theorems may shed some lights on the stepsize choice of the Wasserstein gradient flows which is less studied in the literature. We leave it as an interesting future work.

Compared to SG-MCMC, there are two additional hyperparameters in Algorithm 1: the number of cycles $M$ and the proportion of exploration stage $\beta$. We now study how sensitive they are when comparing to the parallel MCMC. With the same setup as in Section 5.2, We compare our method with $M$ cycles and $L$ epochs per cycle with running $M$ chains parallel MCMC for $L$ epochs. The training budget is 200 epochs. In Table 5.2, $M=4$ and $\beta=0.8$ on CIFAR-10. We compare cSGLD and parallel SGLD with smaller and larger values of $M$ and $\beta$. In Table 6, we see that the conclusion that cSG-MCMC is better than parallel SG-MCMC holds with different values of $M$ and $\beta$.

I.2 Hyperparameters Setting in Practice

Given the training budget, there is a trade-off between the number of cycles $M$ and the cycle length. We find that it works well in practice by setting the cycle length such that the model with optimization methods will be close to a mode after running for that length. (e.g. the cycle length for CIFAR-10 is 50 epochs. The model optimized by SGD can achieve about 5% error after 50 epochs which means the model is close but not fully converge to a mode after 50 epochs.) Once the cycle length is fixed, $M$ is fixed. $\beta$ needs tuning for different tasks by cross-validation. Generally, $\beta$ needs to be tuned so that the sampler has enough time to reach a good region before starting sampling.

Appendix J Tempering in Bayesian Neural Networks

Tempering is common in modern Bayesian deep learning, for both variational inference and MCMC approaches (e.g., Li et al., 2016a; Nguyen et al., 2017; Fortunato et al., 2017). In general, tempering reflects the belief that the model capacity is misspecified. This combination of beliefs with data is what shapes the posterior we want to use to form a good predictive distribution.

Although we use the prescribed temperature in pSGLD (Li et al., 2016a) for all neural network experiments in the main text ($T\approx 0.0045$), we here investigate the effect of temperature $T$ on performance. We show negative log-likelihood (NLL) and classification error as a function of temperature on CIFAR-10 and CIFAR-100 using cSGLD with the same setup as in Section 5.2. We consider $T\in[1,0.5,0.1,0.05,0.01,0.005,0]$. Figure 6 and 7 show the results on CIFAR-10 and CIFAR-100, respectively. On CIFAR-10, the best performance is achieved at $T=0.1$ with NLL 0.1331 and error 4.22%. On CIFAR-100, the best performance is achieved at $T=0.01$ with NLL 0.7835 and error 20.53%. We find that the optimal temperature is often less than 1. We hypothesize that this result is due to the model misspecification common to neural networks.

Indeed, modern neural networks are particularly overparametrized. Tempering enables one to use a model with similar inductive biases to a modern neural network, but with a more well calibrated capacity (which is especially important when we are doing Bayesian integration instead of optimization). Indeed, we show that by sampling from the tempered posterior, we outperform optimization. Learning the amount of tempering by cross-validation is a principled way of aligning the tempering procedure with representing a reasonable posterior. We have shown that sampling with cSGMCMC with tempering helps in terms of both NLL and accuracy, which indicates that we are finding a better predictive distribution.

For both cSGLD and cSGHMC, $M=100$, $\beta=0.01$. For cSGLD, $\alpha_{0}N=1.2,0.5,1.5$ for Austrilian, German and Hear respectively. For cSGHMC $\alpha_{0}N=0.5,0.3,1.0$ for Austrilian, German and Hear respectively. For SG-MCMC, the stepsize is $a$ for the first 5000 iterations and then switch to the decay schedule (1) with $b=0$, $\gamma=0.55$. $aN=1.2,0.5,1.5$ for Austrilian, German and Hear respectively for SGLD and $aN=0.5,0.3,1.0$ for Austrilian, German and Hear respectively for SGHMC. $\eta=0.5$ in cSGHMC and SGHMC.

Assume that we collect $\{\theta_{b}\}_{b=1}^{B}$ samples. Effective sample size (ESS) is computed by

Similar to Hoffman & Gelman (2014), $\hat{\sigma}^{2}$ and $\hat{\mu}$ are obtained by running an independent sampler. We use HMC in this paper.

K.2 Bayesian Neural Networks

For SG-MCMC, the stepsize decays from 0.1 to 0.001 for the first 150 epochs and then switch to the decay schedule (1) with $a=0.01,b=0$ and $\gamma=0.5005$. $\eta=0.9$ in cSGHMC, Snapshot-SGDM and SGHMC.

K.3 Uncertainty Evaluation

For both cSG-MCMC and Snapshot, $M=4$. $\beta=0.8$ in cSG-MCMC. $\alpha_{0}N=0.01$ and $0.008$ for cSGLD and cSGHMC respectively. For SG-MCMC, the stepsize is $a$ for the first 50 iterations and then switch to the decay schedule (1) with $b=0$, $\gamma=0.5005$. $aN=0.01$ for SGLD and $aN=0.008$ for SGHMC. $\eta=0.5$ in cSGHMC, Snapshot-SGDM and SGHMC.