Title: CYCLICAL STOCHASTIC GRADIENT MCMC FOR BAYESIAN DEEP LEARNING

Abstract: The posteriors over neural network weights are high dimensional and multimodal. Each mode typically characterizes a meaningfully different representation of the data. We develop Cyclical Stochastic Gradient MCMC (SG-MCMC) to automatically explore such distributions. In particular, we propose a cyclical stepsize schedule, where larger steps discover new modes, and smaller steps characterize each mode. We prove non-asymptotic convergence of our proposed algorithm. Moreover, we provide extensive experimental results, including ImageNet, to demonstrate the effectiveness of cyclical SG-MCMC in learning complex multimodal distributions, especially for fully Bayesian inference with modern deep neural networks.

1. 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, which account for model uncertainty in weight space (MacKay, 1992; Neal, 1996) , 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) , 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: (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) . (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: (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. (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.

2. 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 D = {d i } N i=1 and a θ-parameterized model, we have the likelihood p(D|θ) and prior p(θ). The posterior distribution is p(θ|D) ∝ exp(-U (θ)) , where U (θ) is the potential energy given by U (θ) = -log p(D|θ) -log p(θ) .

When D is too large, it is expensive to evaluate U (θ) for all the data points at each iteration. Instead, SG-MCMC methods use a minibatch to approximate U (θ): Ũ (θ) = -N N N i=1 log p(x i |θ)log p(θ) , where N N is the size of minibatch. We recommend Ma et al. (2015) for a general review of SG-MCMC algorithms. We describe two SG-MCMC algorithms considered in this paper. Welling & Teh (2011) proposed Stochastic Gradient Langevin Dynamics (SGLD), which uses stochastic gradients with Gaussian noise. Posterior samples are updated at the k-th step as:

SGLD & SGHMC:

θ k = θ k-1 -α k ∇ Ũ (θ k ) + √ 2α k k

, where α k is the stepsize and k has a standard Gaussian distribution.

To improve mixing over SGLD, Stochastic Gradient Hamiltonian Monte Carlo (SGHMC) (Chen et al., 2014) introduces an auxiliary momentum variable v. SGHMC is built upon HMC, with an additional friction term to counteract the noise introduced by a mini-batch. The update rule for posterior samples is:

θ k = θ k-1 + v k-1 , and v k = v k-1 -α k ∇ Ũ (θ k ) -ηv k-1 + 2(η -γ)α k k ,

where 1 -η is the momentum term and γ is the estimate of the noise.

To guarantee asymptotic consistency with the true distribution, SG-MCMC requires that the step sizes satisfy the following assumption: Assumption 1. The step sizes {α k } are decreasing, i.e., 0 < α k+1 < α k , with 1)

∞ k=1 α k = ∞; and 2) ∞ k=1 α 2 k < ∞.

Without a decreasing step-size, the estimation error from numerical approximations is asymptotically biased. One typical decaying step-size schedule is α k = a(b + k) -γ , with γ ∈ (0.5, 1] and (a, b) some positive constants (Welling & Teh, 2011) .

3. 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:

EQUATION

where α 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 α k varies periodically with k. In each period, α k starts at α 0 , and gradually decreases to 0. 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.

Related work in optimization. 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.

3.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(θ|D) ∝ exp(-U (θ)/T ). Note that the setting T = 1 corresponds to sampling from the untempered posterior. When T → 0, the posterior distribution becomes a point mass. Sampling from lim T →0 exp(-U (θ)/T ) is equivalent to minimizing U (θ); 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

) = mod (k-1, K/M ) K/M

is smaller than a given threshold: r(k) < β. Note that β ∈ (0, 1) balances the proportion of the exploration and sampling stages in cSG-MCMC.

3.2. SAMPLING:

Algorithm 1 Cyclical SG-MCMC.

Input::

The initial stepsize α 0 , number of cycles M , number of training iterations K and the proportion of exploration stage β.

for k = 1:K do α ← α k according to Eq equation 1.

if mod (k-1, K/M ) K/M < β then % Exploration stage θ ← θ -α∇ Ũk (θ) else % Sampling stage Collect samples using SG-MCMC methods Output: Samples {θ k }

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

SG-MCMC with Warm Restarts. 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: (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. (ii) The stepsize decreases more quickly inside one cycle than a traditional schedule, making the sampler better characterize the density of the local regions. (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 → 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.

4. THEORETICAL ANALYSIS:

Our algorithm is based on the SDE characterizing the Langevin dynamics:

dθ t = -∇U (θ t )dt + √ 2dW t , where W t ∈ R d is a d-dimensional Brownian motion.

In this section, we prove nonasymptotic convergence rates for the proposed cSG-MCMC framework with a cyclical stepsize sequence {α k } defined in equation 1. For simplicity, we do not consider the exploration stage in the analysis as that corresponds to stochastic optimization. Generally, there are two different ways to describe the convergence behaviours of SG-MCMC. One characterizes the sample average over a particular test function (e.g., Chen et al. (2015) ; Vollmer et al. (2016) ); the other is in terms of the Wasserstein distance (e.g., Raginsky et al. (2017) ; Xu et al. (2017) ). We study both in the following.

Weak convergence Following Chen et al. (2015) and Vollmer et al. (2016) , we define the posterior average of an ergodic SDE as: φ X φ(θ)ρ(θ)dθ for some test function φ(θ) of interest. For the corresponding algorithm with generated samples (θ k ) K k=1 , we use the sample average φ defined as φ = 1 K K k=1 φ(θ k ) to approximate φ. We prove weak convergence of cSGLD in terms of bias and MSE, as stated in Theorem 1. Theorem 1. Under Assumptions 2 in the appendix, for a smooth test function φ, the bias and MSE of cSGLD are bounded as:

BIAS: E φ -φ = O 1 α 0 K + α 0 , MSE: E φ -φ 2 = O 1 α 0 K + α 2 0 .

(2)

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:

W 2 2 (µ, ν) := inf γ Ω×Ω θ -θ 2 2 dγ(θ, θ ) : γ ∈ Γ(µ, ν)

where Γ(µ, ν) is the set of joint distributions over (θ, θ ) such that the two marginals equal µ and ν, respectively.

Denote the distribution of θ t in the SDE as ν t . According to Chiang & Hwang (1987) , the stationary distribution ν ∞ matches our target distribution. Let µ 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 (µ K , ν ∞ ). 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. Theorem 2. 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 (µ K , ν ∞ ) ≤ C 3 exp(- Kα 0 2C 4 ) + 6 + C 2 Kα 0 2 1 2 [(C 1 3α 2 0 K 8 + σC 0 Kα 0 2 ) 1 2 + (C 1 3α 2 0 K 16 + σC 0 Kα 0 4 ) 1 4 ] .

Particularly, if we further assume Remark 1. i) The bound is decomposed into two parts: the first part measures convergence speed of exact solution to the stationary distribution, i.e., ν k α k to ν ∞ ; the second part measures the numerical error, i.e., between µ K and ν k α k . 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) ). 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.

α 0 = O(K -β ) for ∀β > 1, W 2 (µ K , ν ∞ ) ≤ C 3 + 6 + C2 K β-1 1 2 [( 2C1 K 2β-1 + 2C0 K β-1 ) 1 2 + ( C1 K 2β-1 + C0 K β-1 ) 1 4 ].

5. 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 α k ∝ 0.05k -0.55 . In cSGLD, we set M = 30 and the initial stepsize α 0 = 0.09. The proportion of exploration stage β = 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.

5.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 (i) traditional SG-MCMC; (ii) traditional stochastic optimization methods, including stochastic gradient descent (SGD) and stochastic gradient descent with momentum (SGDM); and (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 (±) from the mean predictor.

We set M = 4 and α 0 = 0.5 for cSGLD, cSGHMC and Snapshot. The proportion hyper-parameter β =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 time1 . 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. 

CIFAR:



Testing Performance for Image Classification:

We report the testing errors in Table 1 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(θ) be the test error of a sample with parameter θ. We compute the test error of the convex combination of two samples

J(λθ 1 + (1 -λ)θ 2 ), where λ ∈ [0, 1].

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 λ 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.

Comparison to Parallel 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 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 (Cycli-cal+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 Cycli-cal+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. Results on ImageNet. 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 3, 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 cy-cle schedule is an effective technique to explore the parameter space, and diversified samples can help prevent over-fitting.

5.3. UNCERTAINTY EVALUATION:

Figure 4 : Empirical CDF for the entropy of the predictive distribution on notMNIST dataset. cSGLD and cSGHMC show lower probability for the low entropy estimate than other algorithms.

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-ofdistribution 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 notM-NIST 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 notM-NIST digits, which gives the maximum entropy.

In Figure 4 , we plot the empirical CDF for the entropy of the predictive distributions on notM-NIST. 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.

6. 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.

A EXPERIMENTAL RESULTS

A.1 SYNTHETIC MULTIMODAL DISTRIBUTION:

The density of the distribution is

F (x) = 25 i=1 λN (x|µ i , Σ),

where λ = 1 25 , µ = {-4, -2, 0, 2, 4} × {-4, -2, 0, 2, 4}, Σ = 0.03 0 0 0.03 .

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 n, we consider this mode covered. On this dataset, we choose r = 0.25 and n = 100. 

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. 

Australian German Heart:



B ASSUMPTIONS B.1 ASSUMPTIONS IN WEAK CONVERGENCE ANALYSIS:

In the analysis, we define a functional ψ that solves the following Poisson Equation:

EQUATION

)

The solution functional ψ(θ k ) characterizes the difference between φ(θ k ) and the posterior average φ for every θ k , thus would typically possess a unique solution, which is at least as smooth as φ 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, ψ, of the Poisson equation equation 3. Assumption 2. ψ and its up to 3rd-order derivatives, D k ψ, are bounded by a function V, i.e.,

D k ψ ≤ H k V p k for k = (0, 1, 2, 3), H k , p k > 0.

Furthermore, the expectation of V on {θ k } is bounded: sup l EV p (θ k ) < ∞, and V is smooth such that sup s∈(0,1) V p (sθ + (1 -s) θ ) ≤ C (V p (θ) + V p (θ )), ∀θ, θ , p ≤ max{2p k } for some C > 0.

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. Assumption 3.

• There exists some constants A ≥ 0 and B ≥ 0, such that U (0) ≤ A and ∇U (0) ≤ B.

• The function U is L U -smooth : ∇U (w) -∇U (v) ≤ L U w -v . • The function U is (m u , b) -dissipative, which means for some m U > 0 and b > 0 w, ∇U (w) ≥ m U w 2 -b.

• There exists some constant δ ∈ [0, 1),

such that E[ ∇ Ũk (w) -∇U (w) 2 ] ≤ 2σ(M 2 U w 2 + B 2 ).

• We can choose µ 0 which satisfies the requirement: κ 0 := log e w 2 µ 0 (w)dw < ∞.

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. 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)

∞ k=1 h k = ∞; and 2) lim K→∞ K k=1 h N +1 k K k=1 h k = 0.

Our prove can be derived by the following results from Chen et al. (2015) . Lemma 1 (Chen et al. (2015) ). Let S K K k=1 h k . Under Assumptions 2 and 4, for a smooth test function φ, the bias and MSE of a decreasing-step-size SG-MCMC with a N th-order integrator at time S L are bounded as:

BIAS: E φ -φ = O 1 S K + K k=1 h N +1 k S K (4) MSE: E φ -φ 2 ≤ C l h 2 k S 2 K E ∆V l 2 + 1 S K + ( K k=1 h N +1 k ) 2 S 2 K .

(5)

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.

Proof of Theorem 1. 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

EQUATION

As a result, for the bias, we have

E φ -φ = O 1 S K + K k=1 h N +1 k S K = O 1 α 0 K + 3α 2 0 K/8 α 0 K = O 1 α 0 K + α 0 .

For the MSE, note the first term l h 2 k S 2 K E ∆V l 2 has a higher order than other terms, thus it is omitted in the big-O notation, i.e.,

E φ -φ 2 = O 1 α 0 K + ( 3α 2 0 K/8 α 0 K ) 2 = O 1 α 0 K + α 2 0 .

This completes the proof.

D PROOF OF THEOREM 2:

Proof of the bound for W 2 (µ K , ν ∞ ) in cSGLD. Firstly, we introduce the following SDE

EQUATION

Let ν t denote the distribution of θ t , and the stationary distribution of equation 34 be p(θ|D), which means ν ∞ = p(θ|D).

EQUATION

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

D.1 W 2 (µ K , ν K k=1 α k )

For the first part, W 2 (µ K , ν K k=1 α k ), our proof is based on the proof of Lemma 3.6 in Raginsky et al. (2017) with some modifications. We first assume E(∇ Ũ (w)) = ∇U (w), ∀w ∈ R d , which is a general assumption according to the way we choose the minibatch. And we define p(t) which will be used in the following proof:

EQUATION

Then we focus on the following continuous-time interpolation of θ k :

EQUATION

)

where

∇ Ũ ≡ ∇ Ũk for t ∈ k i=1 α i , k+1 i=1 α i .

And for each k , θ( k i=1 α i ) and θ k have the same probability law µ k . Since θ(t) is not a Markov process, we define the following process which has the same one-time marginals as θ(t)

EQUATION

Let P t V := L (V (s) : 0 ≤ s ≤ t) and P t θ := L (θ(s) : 0 ≤ s ≤ t) 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 P t V and P t θ :

D KL (P t V P t θ ) = -d P t V log d P t V d P t θ = 1 4 t 0 E ∇U (V (s)) -G s (V (s)) 2 ds = 1 4 t 0 E ∇U (θ(s)) -G s (θ(s)) 2 ds

The last line follows the fact that L(θ(s)) = L(V (s)), ∀s. Then we will let t = K k=1 α k and we can use the martingale property of the integral to derive:

D KL (P K k=1 α k V P K k=1 α k θ ) = 1 4 K-1 j=0 j+1 k=1 α k j k=1 α k E ∇U (θ(s)) -G s (θ(s)) 2 ds ≤ 1 2 K-1 j=0 j+1 k=1 α k j k=1 α k E ∇U (θ(s)) -∇U (θ( q(s) k=1 α i ) 2 ds + 1 2 K-1 j=0 j+1 k=1 α k j k=1 α k E ∇U (θ( q(s) k=1 α i ) -G s (θ( q(s) k=1 α i ) 2 ds ≤ L 2 U 2 K-1 j=0 j+1 k=1 α k j k=1 α k E θ(s) -θ( q(s) k=1 α i ) 2 ds (14) + 1 2 K-1 j=0 j+1 k=1 α k j k=1 α k E ∇U (θ( q(s) k=1 α i ) -G s (θ( q(s) k=1 α i ) 2 ds (15)

For the first part ( 14), we consider some s ∈ [

j k=1 α k , j+1 k=1 α k ),

for which the following holds:

EQUATION

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

E θ(s) -θ( j k=1 α k ) 2 ≤ 3α 2 j+1 E ∇U (θ j ) 2 + 3α 2 j+1 E ∇U (θ j ) -∇ Ũj (θ j ) 2 + 6α j+1 d ≤ 12α 2 j+1 (L 2 U E θ j 2 + B 2 ) + 6α j+1 d

Hence we can bound the first part, (choosing α 0 ≤ 1),

L 2 U 2 K-1 j=0 j+1 k=1 α k j k=1 α k E θ(s) -θ( q(s) k=1 α i ) 2 ds ≤ L 2 U 2 K-1 j=0 12α 3 j+1 (L 2 U E θ j 2 + B 2 ) + 6α 2 j+1 d ≤ L 2 U max 0≤j≤K-1 6(L 2 U E θ j 2 + B 2 ) + 3d ( K-1 j=0 α 2 j+1 ) ≤ L 2 U max 0≤j≤K-1 6(L 2 U E θ j 2 + B 2 ) + 3d 3α 2 0 K 8 (17)

The last line (17) follows from2 equation 6. The second part (15) can be bounded as follows:

1 2 K-1 j=0 j+1 k=1 α k j k=1 α k E ∇U (θ( q(s) k=1 α i ) -G s (θ( q(s) k=1 α i ) 2 ds = 1 2 K-1 j=0 α j+1 E ∇U (θ j ) -∇ Ũ (θ j ) 2 ≤ σ max 0≤j≤K-1 (L 2 U E θ j 2 + B 2 ) K-1 j=0 α j+1 ≤ σ max 0≤j≤K-1 (L 2 U E θ j 2 + B 2 )( α 0 2 K-1 j=0 (cos( πmod(j, K/M ) K/M ) + 1)) ≤ σ max 0≤j≤K-1 (L 2 U E θ j 2 + B 2 )( Kα 0 2 )

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

D KL (µ K ν K k=1 α k ) ≤ D KL (P t V P t θ ) ≤ L 2 U 2 K-1 j=0 j+1 k=1 α k j k=1 α k E θ(s) -θ( q(s) k=1 α i ) 2 ds + 1 2 K-1 j=0 j+1 k=1 α k j k=1 α k E ∇U (θ( q(s) k=1 α i ) -G s (θ( q(s) k=1 α i ) 2 ds ≤ L 2 U max 0≤j≤K-1 6(L 2 U E θ j 2 + B 2 ) + 3d 3α 2 0 K 8 + σ max 0≤j≤K-1 (L 2 U E θ j 2 + B 2 )( Kα 0 2 )

According to the proof of Lemma 3.2 in Raginsky et al. (2017) , we can bound the term Raginsky et al. (2017) , we can fix α 0 ∈ (0, 1 ∧ m U 4M 2

E θ k 2 E θ k+1 2 ≤ (1 -2α k+1 m U + 4α 2 k+1 M 2 U )E θ k 2 + 2α k+1 b + 4α 2 k+1 B 2 + 2α k+1 d β Similar to the statement of Lemma 3.2 in

U:

). Then, we can know that

EQUATION

, where α min is defined as α min α0

2 cos π mod( K/M -1, K/M ) K/M + 1 .

There are two cases to consider.

EQUATION

Now, we have

max 0≤j≤K-1 (L 2 U E θ j 2 + B 2 ) ≤ (L 2 U (κ 0 + 2(1 ∧ α 0 m U α min )(b + 2B 2 + d)) + B 2 ) := C 0

Due to the expression of α0 αmin , C 0 is independent of α 0 . Then we denote the 6L 2 U (C 0 + d) as C 1 and we can derive

D KL (µ K ν K k=1 α k ) ≤ C 1 ( 3α 2 0 K 8 ) + σC 0 ( Kα 0 2 )

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

EQUATION

2 )

EQUATION

2 )

1 2 + ( 3C1α 2 0 K 16 + KσC0α0 4 ) 1 4 ] D.2 W 2 (ν K k=1 α k , ν ∞ )

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 ),

W 2 (ν K k=1 α k , ν ∞ ) ≤ C 3 exp(- K k=1 α k /C 4 )

Now combining the bounds for W 2 (µ K , ν K k=1 α k ) and W 2 (ν K k=1 α k , ν ∞ ), substituting α 0 = O(1/K β ), and noting W 2 (ν K k=1 α k , ν ∞ ) decreases w.r.t. K, we arrive at the bound stated in the theorem.

E RELATION WITH SGLD:

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

EQUATION

where

W 2 (μ K , ν K k=1 h k ) ≤ (6 + h 0 K k=1 1 k ) 1 2 • [(D1h 2 0 π 2 6 + σD0h0 K k=1 1 k ) 1 2 + (D1h 2 0 π 2 16 + σD0 h0 2 K k=1 1 k ) 1 4 ] and W 2 (ν K k=1 h k , ν ∞ ) ≤ C 3 exp(- K k=1 h k C4

).

Proof of the bound of W 2 (μ K , ν ∞ ) in the standard SGLD. Similar to the proof of W 2 (µ K , ν ∞ ) in cSGLD, we get the following update rule for SGLD with the stepsize following a polynomial decay i.e.,

EQUATION

Let μk denote the distribution of θ k .

Since

EQUATION

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

E.1 W 2 (μ K , ν K k=1 h k )

We first assume E(∇ Ũ (w)) = ∇U (w), ∀w ∈ R d , which is a general assumption according to the way we choose the minibatch. Following the proof in Raginsky et al. (2017) and the analysis of the SPOS method in Zhang et al. (2018) , we define the following p(t) which will be used in the following proof:

EQUATION

Then we focus on the following continuous-time interpolation of θ k :

EQUATION

where

∇ Ũ ≡ ∇ Ũk for t ∈ k i=1 h i , k+1 i=1 h i .

And for each k , θ( k i=1 h i ) and θ k have the same probability law μk . Since θ(t) is not a Markov process, we define the following process which has the same one-time marginals as θ(t)

EQUATION

Let P t V := L (V (s) 0 ≤ s ≤ t) and P t θ := L (θ(s) : 0 ≤ s ≤ t) 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 P t V and P t θ :

D KL (P t V P t θ ) = -d P t V log d P t V d P t θ = 1 4 t 0 E ∇U (V (s)) -G s (V (s)) 2 ds = 1 4 t 0 E ∇U (θ(s)) -G s (θ(s)) 2 ds

The last line follows the fact that L(θ(s)) = L(V (s)), ∀s. Then we will let t = K k=1 h k and we can use the martingale property of integral to derive:

D KL (P K k=1 h k V P K k=1 h k θ ) = 1 4 K-1 j=0 j+1 k=1 h k j k=1 h k E ∇U (θ(s)) -G s (θ(s)) 2 ds ≤ 1 2 K-1 j=0 j+1 k=1 h k j k=1 h k E ∇U (θ(s)) -∇U (θ( q(s) k=1 h i ) 2 ds + 1 2 K-1 j=0 j+1 k=1 h k j k=1 h k E ∇U (θ( q(s) k=1 h i ) -G s (θ( q(s) k=1 h i ) 2 ds ≤ L 2 U 2 K-1 j=0 j+1 k=1 h k j k=1 h k E θ(s) -θ( q(s) k=1 h i ) 2 ds (29) + 1 2 K-1 j=0 j+1 k=1 h k j k=1 h k E ∇U (θ( q(s) k=1 h i ) -G s (θ( q(s) k=1 h i ) 2 ds (30)

For the first part (29), we consider some s ∈ [ j k=1 h k , j+1 k=1 h k ), the following equation holds:

EQUATION

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

E θ(s) -θ( j k=1 h k ) 2 ≤ 3h 2 j+1 E ∇U (θ j ) 2 + 3h 2 j+1 E ∇U (θ j ) -∇ Ũj (θ j ) 2 + 6h j+1 d ≤ 12h 2 j+1 (L 2 U E θ j 2 + B 2 ) + 6h j+1 d

Hence we can bound the first part, (choosing h 0 ≤ 1),

EQUATION

)

where the last line follows from the fact that

K-1 j=0 1 (j + 1) 3 ≤ K-1 j=0 1 (j + 1) 2 ≤ ∞ j=0 1 (j + 1) 2 = π 2 6 .

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

1 2 K-1 j=0 j+1 k=1 h k j k=1 h k E ∇U (θ( q(s) k=1 h i ) -G s (θ( q(s) k=1 h i ) 2 ds = 1 2 K-1 j=0 h j+1 E ∇U (θ j ) -∇ Ũ (θ j ) 2 ≤ σ max 0≤j≤K-1 (L 2 U E θ j 2 + B 2 ) K-1 j=0 h j+1 ≤ σ max 0≤j≤K-1 (L 2 U E θ j 2 + B 2 )(h 0 K j=1 1 j )

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

D KL (μ K ν K k=1 h k ) ≤ D KL (P t V P t θ ) ≤ L 2 U 2 K-1 j=0 j+1 k=1 h k j k=1 h k E θ(s) -θ( q(s) k=1 h i ) 2 ds + 1 2 K-1 j=0 j+1 k=1 h k j k=1 h k E ∇U (θ( q(s) k=1 h i ) -G s (θ( q(s) k=1 h i ) 2 ds ≤ L 2 U max 0≤j≤K-1 6(L 2 U E θ j 2 + B 2 ) + 3d π 2 6 h 2 0 + σ max 0≤j≤K-1 (L 2 U E θ j 2 + B 2 )(h 0 K j=1 1 j )

Similar to the proof of cSGLD , we have max 0≤j≤K-1

(L 2 U E θ j 2 + B 2 ) ≤ D 0

Then we denote the 6L 2 U (D 0 + d) as D 1 and we can derive

D KL (μ K ν K k=1 h k ) ≤ D 1 h 2 0 π 2 6 + σD 0 h 0 K j=1 1 j

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

W2(μK , ν K k=1 h k ) ≤ [12 + D2( K k=1 h k )] 1/2 • [(DKL(μK ν K k=1 h k )) 1/2 + (DKL(μK ν K k=1 h k )/2) 1/4 ] = [12 + D2(h0 K j=1 1 j )] 1/2 • [(D1h 2 0 π 2 6 + σD0h0 K j=1 1 j ) 1/2 + (D1h 2 0 π 2 12 + σD0h0 K j=1 1 2j ) 1/4 ] Now we derive the bound for W 2 (μ K , ν K k=1 h k ). E.2 W 2 (ν K k=1 h k , ν ∞ )

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 ),

W 2 (ν K k=1 h k , ν ∞ ) ≤ C 3 exp(- K k=1 h k /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 (µ K , ν K k=1 α k ) in the equation 9 increases w.r.t. K, our α 0 must be set small enough in practice. Hence, in this informal comparison, we also set α 0 small enough to make W 2 (µ K , ν K k=1 α k ) less important. i) If the initial stepsizes satisfy α 0 ≥ h 0 , our algorithm cSGLD runs much faster than the standard SGLD in terms of the amount of "diffusion time" i.e., the "t" indexing θ t in the continuous-time SDE mentioned above. This result follows from

K k=1 α k = Kα0 2 and K k=1 h k = K k=1 h0 k = O(h 0 log K) Kα0 2 .

In standard SGLD, since the error described by W 2 (μ K , ν K k=1 h k ) increases w.r.t. K, h 0 needs to be set small enough in practice to reduce the error. Following the general analysis of SGLD in Raginsky et al. (2017) ; Xu et al. (2017) , the dominant term in the decomposition equation 21 will be W 2 (ν K k=1 h k , ν ∞ ) since it decreases exponentially fast with the increase of t and W 2 (μ K , ν K k=1 h k ) is small due to the setting of small h 0 . Since K k=1 α k increases much faster in our algorithm than the term K k=1 h k in standard SGLD, our algorithm thus endows less error for K iterations, i.e., W 2 (ν K k=1 α k , ν ∞ ) W 2 (ν K k=1 h k , ν ∞ ). Hence, our algorithm outperforms standard SGLD, as will be verified in our experiments.

ii) Instead of setting the h 0 small enough, one may consider increasing h 0 to make standard SGLD run as "fast" as our proposed algorithm, i.e., K k=1 h k ≈ K k=1 α k . Now the W 2 (ν K k=1 h k , ν ∞ ) in equation 21 is almost the same as the W 2 (ν K k=1 h k , ν ∞ ) in equation 9. However, in this case, it is worth noting that h 0 scales as O(α 0 K/ log K). We can notice that h 0 is much larger than the α 0 and thus the W 2 (μ K , ν K k=1 h k ) cannot be ignored. Now the h 2 0 term in W 2 (μ K , ν K k=1 h k ) would scale as O(α 2 0 K 2 / log 2 K), which makes W 2 (μ K , ν K k=1 α k ) in equation 21 much larger than our W 2 (µ K , ν K k=1 α k ) defined in equation 9 since O(α 2 0 K 2 / log 2 K) O(α 2 0 K). Again, our algorithm cSGLD achieves a faster convergence rate than standard SGLD.

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(θ|D) on a metric space Θ. As we have M cycles in total, the mth cycle characterizes a local region Θ m ⊂ Θ, defining the "sub-posterior" distribution: p m (θ|D) = ),

where K m is the number of samples from the mth cycle, and θ (m) ∈ Θ m .

The weight for each cycle w i is estimated using the harmonic mean method (Green, 1995; Raftery et al., 2006) 

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 (2) with a = 0.01, b = 0 and γ = 0.5005. η = 0.9 in cSGHMC, Snapshot-SGDM and SGHMC.

K.3 UNCERTAINTY EVALUATION:

For both cSG-MCMC and Snapshot, M = 4. β = 0.8 in cSG-MCMC. α 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 (2) with b = 0, γ = 0.5005. aN = 0.01 for SGLD and aN = 0.008 for SGHMC. η = 0.5 in cSGHMC, Snapshot-SGDM and SGHMC.

:

https://github.com/ChunyuanLI/pSGLD/issues/2

Note: we only focus on the case when K mod M = 0.

Figures and Tables:

Figure 1: Illustration of the proposed cyclical stepsize schedule (red) and the traditional decreasing stepsize schedule (blue) for SG-MCMC algorithms.

Figure 3: Results of cSG-MCMC with DNNs on the CIFAR-100 dataset. (a) MDS visualization in weight space: cSG-MCMC show larger distance than traditional schedules. (b) Testing errors (%) on the path of two samples: cSG-MCMC shows more varied performance. (c) Testing errors (%) as a function of the number of cycles M : cSGLD yields consistently lower errors.

Figure 2: Sampling from a mixture of 25 Gaussians shown in (a) for the parallel setting. With a budget of 50k × 4 = 200k samples, traditional SGLD in (b) has only discovered 4 of the 25 modes, while our cSGLD in (c) has fully explored the distribution.

Figure 5: Sampling from a mixture of 25 Gaussians in the non-parallel setting. With a budget of 50K samples, traditional SGLD has only discovered one of the 25 modes, while our proposed cSGLD has explored significantly more of the distribution.

p(θ|D)1Θ m wm, with w m = Θm p(θ|D)dθ, where w m is a normalizing constant. For a testing function f (θ), we are often interested in its true posterior expectation f = f (θ)p(θ|D)dθ. The sample-based estimation is

Figure7: NLL and error (%) as a function of temeprature on CIFAR-100 using cSGLD. The best performance of both NLL and error is achieved at T = 0.01.

Comparison of test error (%) between cSG-MCMC with parallel algorithm (M =4 chains) on CIFAR-10 and CIFAR-100. The method is reported in the format of "step-size schedule (cyclical or decreasing) + single/parallel chain". The cost is reported in the format of "#epoch per chain / #epoch used in all chains". Note that a parallel algorithm with a single chain reduces to a nonparallel algorithm. Integration of the cyclical schedule with parallel algorithms provides lower testing errors.

Table content: <table><tr><td>Method</td><td colspan="6">Cyclical+Parallel Decreasing+Parallel Decreasing+Parallel</td><td colspan="2">Cyclical+Single</td></tr><tr><td>Cost</td><td colspan="2">200/800</td><td colspan="2">200/800</td><td colspan="2">100/400</td><td colspan="2">200/200</td></tr><tr><td>Sampler</td><td colspan="3">SGLD SGHMC SGLD</td><td>SGHMC</td><td>SGLD</td><td>SGHMC</td><td colspan="2">SGLD SGHMC</td></tr><tr><td>CIFAR-10</td><td>4.09</td><td>3.95</td><td>4.15</td><td>4.09</td><td>5.11</td><td>4.52</td><td>4.29</td><td>4.27</td></tr><tr><td>CIFAR-100</td><td>19.37</td><td>19.19</td><td>20.29</td><td>19.72</td><td>21.16</td><td>20.82</td><td>20.55</td><td>20.50</td></tr></table>



Table content: <table><tr><td/><td>NLL ↓ Top1 ↑ Top5 ↑</td></tr><tr><td>SGDM</td><td>0.9595 76.046 92.776</td></tr><tr><td colspan="2">Snapshot-SGDM 0.8941 77.142 93.344</td></tr><tr><td>SGHMC</td><td>0.9308 76.274 92.994</td></tr><tr><td>cSGHMC</td><td>0.8882 77.114 93.524</td></tr></table>

Table 4 shows the mode-coverage for several algorithms, based on 10 different runs. Mode coverage over 10 different runs, ± standard error.

Table content: <table><tr><td>Algorithm</td><td>Mode coverage</td></tr><tr><td>SGLD</td><td>1.8±0.13</td></tr><tr><td>cSGLD</td><td>6.7±0.52</td></tr><tr><td>Parallel SGLD</td><td>18±0.47</td></tr><tr><td>Parallel cSGLD</td><td>24.4±0.22</td></tr></table>

Effective sample size for samples for the unimodal posteriors in Bayesian linear regression, obtained using cyclical and traditional SG-MCMC algorithms, respectively.

Table content: <table><tr><td>SGLD</td><td>1676</td><td>492</td><td>2199</td></tr><tr><td>cSGLD</td><td>2138</td><td>978</td><td>2541</td></tr><tr><td>SGHMC</td><td>1317</td><td>2007</td><td>5000</td></tr><tr><td>cSGHMC</td><td>4707</td><td>2436</td><td>5000</td></tr></table>

Peer review (iclr2020):
The paper propose a new MCMC scheme which is demonstrated to perform well for estimating Bayesian neural networks. The key idea is to not keep lowering the step sizes, but -- at pre-specified times -- go back to large step sizes. The paper is timely, the proposed algorithm is novel, and the theoretical analysis also seem quite novel. My key concern is that with MCMC sampling it is often quite difficult to tune parameters, and by introducing more parameters to tune when step sizes should increase, I fear that we end up in a *tuning nightmare*. How sensitive is the algorithm to choice of parameters? I would expect that the proposed algorithm is quite similar to just running several MCMCs in parallel. The authors does a comparison to this and show that their approach is significantly faster due to *warm restarts*. Here I wonder how sensitive this conclusion is to choice of parameters (see nightmare above) ? I would guess that opposite conclusions could be reached by tuning the algorithms differently -- is that a reasonable suspicion ? It is argued that the cyclic nature of the algorithms gives a form of *warm start* that is beneficial for MCMC. My intuition dictate that this is only true of the modes of the posterior are reasonable close to each other; otherwise I do not see how this warm starting is helpful. I would appreciate learning more about why this intuition is apparently incorrect. Minor comments: * on page 4 it is stated that the proposed algorithm *automatically* provide the warm restarts -- but is it really automatic? Isn*t this a priori determined by choice of parameters for the algorithm? * It would be good to use citet instead of cite at places, e.g. *discussed in (Smith & Topin, 2017)* should be *discussed by Smith & Topin (2017)*. This would improve readability (which is generally very good). * For the empirical studies I think it would be good to report state-of-the-art results as well. I expect that the Bayesian nets still are subpar to non-Bayesian methods, and I think the paper should report this.