Interacting Contour Stochastic Gradient Langevin Dynamics
Authors: Wei Deng
Siqi Liang
Botao Hao
Guang Lin
Faming Liang
Venue: ICLR 2022 Poster (2022)
Keywords: stochastic gradient Langevin dynamics
MCMC
importance sampling
Wang-Landau algorithm
Parallel MCMC Methods
stochastic approximation

Abstract: We propose an interacting contour stochastic gradient Langevin dynamics (ICSGLD) sampler, an embarrassingly parallel multiple-chain contour stochastic gradient Langevin dynamics (CSGLD) sampler with efficient interactions. We show that ICSGLD can be theoretically more efficient than a single-chain CSGLD with an equivalent computational budget. We also present a novel random-field function, which facilitates the estimation of self-adapting parameters in big data and obtains free mode explorations. Empirically, we compare the proposed algorithm with popular benchmark methods for posterior sampling. The numerical results show a great potential of ICSGLD for large-scale uncertainty estimation tasks.

Published as a conference paper at ICLR 2022
INTERACTING CONTOUR STOCHASTIC GRADIENT LANGEVIN DYNAMICS
Wei 2 Siqi Liang', Botao Hao?, Guang Lin', Faming Liang' Purdue University MMorgan Stanley DDeppMind fml mnllang@purdee..do edu; weidengC56@gma11 com
ABSTRACT
We propose interacting contour stochastic gradient Langevin dynamics (IC- SGLD) sampler, an embarrassingly parallel multiplc-chain contour stochastic gradient Langevin dynamics (CSGLD) sampler with efficient interactions. We show that ICSGLD can be theoretically more efficient than a single-chain CSGLD with an equivalent computational budget. We also present a novel random-field function, which facilitates the estimation of self-adapting parameters in big data and obtains free mode explorations. Empirically, we compare the proposed algorithm with popular benchmark methods for posterior sampling. The numerical results show a great potential of ICSGLD for large-scale uncertainty estimation tasks.
INTRODUCTION
Stochastic gradient Langevin dynamics (SGLD) (Welling & Teh, 2011) has achieved great successes in simulations of high- nigh-dimensional systems for big data problems. It, however, yields only a fast mixing rate when the energy landscape is simple, e.g., local energy wells are shallow and not well separated. To improve its convergence for the problems with complex energy landscapes, various strategies have been proposed, such as momentum augmentation (Chen et al., 2014; Ding et al., 2014), Hessian approximation (Ahn et al., 2012: Li al., 2016), high-order numerical schemes (Chen et al., 2015; Li et al., 2019b), and cyclical learning rates (Izmailov et al., 2018; Maddox et al., 2019; Zhang al., 2020b). In spite of their asymptotic properties in Bayesian inference (Vollmer et al.. 2016) and non-convex optimization (Zhang et al., 2017), it is still difficult to achieve compelling empirical results for pathologically complex deep neural networks (DNNs).
To simulate from distributions with complex energy landscapes, e.g., those with multitude of modes well separated by high energy barriers, an emerging trend run multiple chains, where interactions between different chains can potentially accelerate the convergence of the simulation. For example, Song et al. (2014) and Futami et al. (2020) showed theoretical advantages of appropriate interactions in ensemble/popullttoon simulations. Other multiple chain methods include particle-based based nonlinear Markov ( (Vlasov) processes (Liu & Wang, 2016; Zhang et al., 2020a) and replica exchange methods (also known parallel tempering) (Deng et al., 2021a). However, the particle-based methods result in an expensive kernel matrix computation given a large number of particles (Liu & Wang, 2016); similarly, naïvely extending replica exchange methods to population chains leads to a long waiting time to swap between non-neighboring chains (Syed et al., 2021). Therefore, how to conduct interactions between different chains, while maintaining the scalability of the algorithm, is the key to the success of the parallel stochastic gradient MCMC algorithms.
In this paper, we propose an interacting contour stochastic gradient Langevin dynamics (ICSGLD) sampler, a pleasingly parallel extension of contour stochastic gradient Langevin dynamics (CSGLD) (Deng et al., 2020b) with efficient interactions. The proposed algorithm requires minimal communi- cation cost in that each chain shares with others the marginal energy likelihood estimate only. As result, the interacting mechanism improves the convergence of the simulation, while the minimal communication mode between different chains enables the proposed algorithm to be naturally adapted to distributed computing with little overhead. For the single-chain CSGLD algorithm, despite its theoretical advantages as shown in Deng et al. (2020b), estimation of the marginal energy likelihood remains challenging for big data problems with a wide energy range, jeopardizing the empirical performance of the class of importance sampling methods (Gordon al., 1993; Doucet al., 2001;

---

Published as a conference paper at ICLR 2022
Wang & Landau, 2001; Liang al., 2007; Andrieu al., 2010; Deng al., 2020b) in big data applications. To resolve this issue, we resort to a novel interacting random-fiee field function based on multiple chains for an ideal variance reduction and a more robust estimation. As such, we can greatly facilitate the estimation of the marginal energy likelihood so as to accelerate the simulations of notoriously complex distributions. To summarize, the algorithm has three main contributions:
We propose scalable interacting importance sampling method for big data problems with the min- imal communication cost. A novel random-field function is derived to tackle the incompatibility issue of the class of importance sampling methods in big data problems.
Theoretically, we study the local stability of non-linear mean-field system and justify regularity properties of the solution of Poisson's 's equation. We also prove the asymptotic normality for the stochastic approximation process in mini-batch settings and show that ICSGLD is asymptotically more efficient than the single-chain CSGLD with an equivalent computational budget.
Our proposed algorithm achieves appealing mode explorations using fixed learning rate on the MNIST dataset and obtains remarkable performance in large-scale uncertainty estimation tasks.
2
PRELIMINARIES
2.1
STOCHASTIC GRADIENT LANGEVIN DYNAMICS
standard sampling algorithm for big data problems is SGLD (Welling & Tch, 2011), which is a numerical scheme of stochastic differential equation in mini-batch settings:
N Ek
X X
(xx) V
where X X E Rd, CA is the learning rate at iteration k, N denotes the number of total data points, is the temperature, and wy is standard Gaussian vector of dimension d. In particular, is an unbiased stochastic gradient estimator based on a mini- batch data B of size n and N (x) is the unbiased energy estimator for the exact energy function Uxx). Under mild conditions on U, X is ((e)
N V (x) a
known to converge weakly to a unique invariant distribution (x) x e as E + 0.
2.2 CONTOUR STOCHASTIC GRADIENT LANGEVIN DYNAMICS
Despite its theoretical guarantees, SGLD can converge exponentially slow when U(x) is non-convex and exhibits high energy barriers. To remedy this issue, CSGLD (Deng et al., 2020b) exploits the flat histogram idea and proposes to simulate from flattened density with much lower energy barriers
W o (x) x Txx) / (U)))),
(2)
where Ç is a hyperparameter, (o( = Di= 0(i ((ii-)e(log) -lgg (() (1i 2u
In particular, {u; ) determines the partition { of X' such that X, = fx U;-1 < U(x) < u where 00 Ue < U1 < < um-1 < Um co. For practical purposes, we assume Wi+1 u, Au for i 1, ... - 2. In addition, 0 (()),, (())),..)()))) is the self-adapting
0(m))
parameter in the space + = (0(1) ... ((m)) (((),.....0(m))) 10 < ((1) ((m) < 1& Em 0(1) =
Ideally, setting < = and (()) = (1), (t), where 2x = fx. (xx)xx for i E (,, m). enables CSGLD to achieve a "random walk" in the space of energy and to penalize the over-visited partition (Wang& & Landau, 2001; Liang et al., 2007; Fort et al., 2011; 2015). However, the optimal values of O. is unknown a priori. To tackle this issue, CSGLD proposes the following procedure to adaptively estimate 0 via stochastic approximation (SA) (Robbins & Monro, 1951; Benveniste et al., 1990):
(1) Sample X++1 Xx + N V Uwo (xx) + 22T
Optimize 0.+1 02 + wh+ H( xx+1);
2

---

Published as a conference paper at ICLR 2022
Algorithm 1 Interacting contour stochastic gradient Langevin dynamics algorithm (ICSGLD).
{ 1m is pre-defined partition and Ç is hyperparameter. The update rule distributed-memory settings and discussions of hyperparameters is detailed in section B.1.1 in the supplementary material.
[1.] (Data subsampling) Draw a mini-batch data B& from D, and compute stochastic gradients
Vx VxU(xx) (P)) and energies Ü (xxx)
(x) for each x(2), where p E {,, ... ; P), [Byl = n, and [D] = N.
[2.] (Parallel simulation) Sample @ P = (x x k+11 x based on SGLD and 0g
(1) (2)
(P)
x
N + ck VRÜNO Mop x
P + v w
(4)
n
where eg is the learning rate, T is the temperature, wk denotes P independent standard Gaussian vectors, VxUv, RxU,v((8P P) (VxUv, (x(1)), Vx (x(2)).... V.Uva (x(P)))T, and
VxUvo, [+ du (log (J (x)) log (J) (x) 1) V 1)) (xÚ(xx for any E X.
[3.] (Stochastic approximation) Update the self- adaltidg parameter 0(i) for (1,2,..., m}
- (1) + wk-1 E (x())((((11)) (1 0x (i))
where 1 is an indicator function that takes value if the event appears and equals 0 otherwise.
where VxU (-) is stochastic gradient function of wv (-) to be detailed in Algorithm 1. [i((0, x) == (0, x). Hm (0,x)) is random-field function where each entry follows
m
=0 (J)(xx) (1==)g(x) -(i)) where (x = Σ il ال ال 4 (0()SS
(3)
Theoretically, CSGLD converges proaches to a fixed point dx and the samples are drawn from the flattened density Wee (x). Notably, the mean-field system is globally stable with unique stable equilibrium point in small neighborhood of Ox.. Moreover, such an appealing property holds even when U(x) is non-convex.
campling-optimizatiee equilibrium in the sense that 0g ap-
3 INTERACTING CONTOUR STOCHASTIC GRADIENT LANGEVIN DYNAMICS
The major goal of interacting CSGLD (ICSGLD) is to improve the efficiency of CSGLD. In particular, the self-adapting parameter 0 is crucial for ensuring the sampler to escape from the local traps and traverse the whole energy landscape, and how to reduce the variability of 0g, 's is the key to the success of such dynamic importance sampling algorithm. To this end, we propose an efficient variance reduction scheme via interacting parallel systems to improve the accuracy of 0 's.
3.1 INTERACTIONS IN PARALLELISM
Now we first consider naive parallel sampling scheme with P chains as follows
P P + Ck
V2Üva Vo (x (x)P) + 2 P
++1
where x8 (x(1), 22) ... PP) w 2 P denotes P independent standard Gaussian vectors, and ÜV (x = (Uno (x!) ()) Tw (x (2)))) )..),,............ Uvo (x())))
Stochastic approximation aims to find the solution O of the mean-field system h(0) such that
(())= = / ())))))=(d 0,
where we is the invariant measure simulated via SGLD that approximates wv, in (2) and H x) is the novel random-fiele function to be defined later in (8). Since h(0) is observable only up to large random perturbations (in the form of ((0,))). the optimization of 0 based isolated
3

---

Published as a conference paper at ICLR 2022
random-field functions may not be efficient enough. However, due the conditional independence of x x(2), 1... xPP) in parallel sampling, it is very natural to consider Monte Carlo average
h()) = ア Σ ir I (dx (dx(())) =
(6)
Namely, we considering the following stochastic approximation scheme
0x+1 02 wk+1 H (0,, P).
(7)
where (0g, xx+1) is an interacting random-fies faedddfieed function H (0,, xx11) = H x
Note that the Monte Carlo average is very effective to reduce the variance of the interacting random- field function H (0, x x0 P)) based on the conditionally independent random field functions. Moreover, each chain shares with others only a very short message during each iteration. Therefore, the interacting parallel system is well suited for distributed computing. where the implementations and communication costs are further detailed in section B.1.2 in the supplementary material. By contrast, each chain of the non-interacting parallel CSGLD algorithm deals with the parameter 0 and large-variance random-field function H (0, x) individually, leading coarse estimates in the end.
Formally, for the population/ensemble interaction scheme (7), we define a novel random mandom--ield function ([0, (H,(0) (0, x), (0, x). Hm (0, x)) where each component satisfies
(0. 1(1g((())) (x)) (l==J)(x)) 0(1)
As shown in Lemma 1, the corresponding mean-field function proposes to converge to different fixed point 0., s.t.
(1) x (L
dx
x (
(9)
large data set often renders the task of estimating & numerically challenging. By contrast, we resort to different solution by estimating O, instead based on large value of Ç. The proposed algorithm is summarized in Algorithm 1. For more study on the scalablity of the new scheme, we leave the discussion in section B.1.3.
3.2 RELATED WORKS
Replica exchange SGLD (Deng et al., 2020a; 2021a) has successfully extended the traditional replica exchange (Swendsen & Wang, 1986; Geyer, 1991; Earl & Deem, 2005) big data problems. However, it works with two chains only and has low swapping rate. As shown in Figure ((a), naive extension of multi-chain replica exchange SGLD yields low communication efficiency. Despite some recipe in the literature (Katzgraber et al., 2008; Bittner et al., 2008; Syed et al., 2021), how to conduct multi-chain replica exchange with low-frequency swaps is still an open question.
(a) Replica Exchange (parallel tempering)
(b) Interacting contour SGLD (ICSGLD)
Figure 1: A comparison of communication costs between replica exchange (RE) and ICSGLD. We see RE takes many iterations to swap with all the other chains; by contrast, ICSGLD possesses pleasingly parallel mechanism where the only cost comes from sharing light message.
Stein variational gradient descent (SVGD) (Liu & Wang, 2016) a popular approximate inference method to drive set of particles for posterior approximation. In particular, repulsive forces are proposed to prevent particles to collapse together into neighboring regions, which resembles our strategy of penalizing over- vesiied partition. However, SVGD tends to underestimate the uncertainty given a limited number of particles. Moreover, the quadratic cost in kernel matrix computation further raises the scalability concerns as particles are proposed.
4

---

Published as a conference paper at ICLR 2022
Admittedly, ICSGLD is not the first interacting importance sampling algorithm. For example, population stochastic approximation Monte Carlo (pop-SAMs (op-SAMC) algorithm has been proposed in Song et al. (2014), and an interacting particle Markov chain Monte Carlo (IPMCMC) algorithm has been proposed Rainforth et al. (2016). A key difference between our algorithm and others is that our algorithm is mainly devised for big data problems. The IPMCMC and pop-SAMC are gradient-free free samplers, which are hard to be adapted to high-dimensional big data problems.
Other parallel SGLD methods (Ahn al., 2014; Chen et al., 2016) aim to reduce the computational cost of gradient estimations in distributed computing, which, however, does not consider interactions for accelerating the convergence. Li et al. (2019a) proposed asynchronous protocols to reduce communication costs when the master aggregates model parameters from all workers. Instead, we don't communicate the parameter R but only share f E R and the indices, where m < d.
Our work also highly resembles the well-known Federated Averaging (FedAvg) algorithm (Li et al., 2020; Deng et al., 2021b), except that the stochastic gradient Ü (x) is replaced with the random field function H (0, x) and we only share the low-dimensional latent vector 0. Since privacy concerns and communication cost are not major bottlenecks of our problem, we leave the study of taking the Monte Carlo average (6) every K iterations for future works.
4 CONVERGENCE PROPERTIES
To study theoretical properties of ICSGLD, we first show a local stability property that is well-suited to big data problems, and then we present the asymptotic normality for the stochastic approxi- mation process in mini-batch settings, which eventually yields the desired result that ICSGLD is asymptotically more efficient than a single-chain CSGLD with an equivalent computational cost.
4.1 LOCAL STABILITY FOR NON-LINEAR MEAN-FIELD SYSTEMS IN BIG DATA
The first obstacle for the theoretical study is to approximate the components of O.e corresponding to the high energy region. To get around this issue, the random field function H x) in (8) is adopted
to estimate different target 0. 0%. As detailed in Lemma 3 in the supplementary material, the mean-field equation is now formulated as follows
h,) x (1) - (B(1)Co) + perturbations,
(10)
where Co
and 7,.00 = Dk= (-TTK) We see that (10) may not be linearly stable 27.0
20.0
Sey. =(x)dx
as in Deng al. (2020b). Although the solution of the mean-field system h()) = 0 is still unique, there may exist unstable invariant subspaces, leading us to consider the local properties. For a proper initialization of 0, which can be achieved by pre-training the model long enough time through SGLD, the mean value theorem implies a linear property in local region
ha (0) x (i) - 0(i) + perturbations. Combining the perturbation theory (Vanden-Eijndens 2001), we present the following stability result:
Lemma (Local stability, informal version of Lemma 3) Assume Assumptions A/-A4 (given in the supplementary material) hold. For any properly initialized O, we have (h(0), 0 O.) <
-6l 0 0_1,, where 0. 0. + O (sup, Var(5n + 1.), O, 0% Φ learning rate, and (x) denotes the noise in the stochastic energy estimator of batch size and Var(-) denotes the variance.
0, denotes
By justifying the drift conditions of the adaptive transition kernel and relevant smoothness properties, we can prove the existence and regularity properties of the solution of the Poisson's 's equation in Lemma 6 in the supplementary material. In what follows, we can control the fluctuations in stochastic approximation and eventually yields the 2 convergence.
Lemma 2 2L2 convergence rate, informal version of Lemma 7) Given standard Assumptions Al- A5, Og converges to Q,, where O. = 0, + O (sup, Var(s, + + 1.), such that
E 110 0.12² O (wx).
5

---

Published as a conference paper at ICLR 2022
The result differs from Theorem of Deng al. (2020b) in that the biased fixed point 0. instead of 0, is treated as the equilibrium of the continuous system, which provides us a user- user-friendly proof. Similar techniques have been adopted by Durmus Eric Moulines (2017); Xu et al. (2018). Although
the global stability (Deng et al., 2020b) may be sacrificed when ( based on Eq.(8), 0. x 0% is much easier to estimate numerically for any that yields 0x (i) based on a large > 1.
4.2 ASYMPTOTIC NORMALITY
To study the asymptotic behavior of w/ (0) 0.), where 0, is the equilibrium point s.t. 0, 0. + (Var(Sx(x)) + c+ 2 we consider fixed step size the SA step for ease of explanation. Let O, denote the solution of the mean-field system in continuous time (0g Oo), and rewrite the single-chain SA step (7) as follows
0 k+1 O(++I) = 02 Okw (H Xk+1) H (Q kw Xk+1))
+ w (H XR+1) - h(0kmw)) (O(k+1)w Onw whw (0
Further, we set Onw = - (02 - -nw) Then the stochastic approximation differs from the mean field system in that
0 (k+1)w (000 10 x++1)) tw?
(m x.+1) - h(0ww)) - remainder
martingale M,
pertorlations
[k+1)a
(k+1) R²d (O.)dWW.
Rd ho (0 (ww) (0, Oew
Mi 20
ho(0.)0.ds
Lw
where (0) = d (0) is matrix, W E R a standard Brownian motion, the last term follows from certain central limit theorem (Benveniste al., 1990) and R denotes the covariance matrix of the random-field function s.t. R(0) = Ei- Cove (H, xt). H(0, xn)).
20
We expect the weak convergence of U2: the stationary distribution of diffusion
dU2 = ha + (O.)dd.)
(11)
where U, = WT 1/2 θ.). Given that 6, converges to θ, sufficiently fast and the local linearity of he, the diffusion (11) resembles the Ornstein-Uhlenteee
process and yields the following solution
U. ~ e (h.) + (." e (ttss))oo(..) o R(O. dW,
Then we have the following theorem, whose formal proof given in section C.3.
Theorem (Asymptotic Normality) Assume Assumptions A/-A5 (given in the supplementary terial) hold. We have the following weak convergence
Wi 1/2 0.) + (0, E), where Σ -= / the, R th o. dl, ho. = ho)(..)
4.3 INTERACTING PARALLEL CHAINS ARE MORE EFFICIENT
For clarity, first denote estimate of 9 based ICSGLD with P interacting parallel chains by of and denote the estimate based on a single-long-chain CSGLD by Oxp.
Note that Theorem holds for any step size wa O(k-a), where (0.5, 1]. If simply run single-chain CSGLD algorithm with times of iterations, by Theorem 1,
WEP (0 0.) - (O.
As to ICSGLD, since the covariance Σ relies on R. which depends on the covariance of the martingale [M, the conditional independence of x(1), x²) ... (P) naturally results in an efficient variance reduction such that

---

Published as a conference paper at ICLR 2022
Corollary 1 (Asymptotic Normality for ICSGLD) Assume the same assumptions. For ICSGLD with P interacting chains, we have the following weak convergence
1/2 0 0.) N/O,
That is, under a similar computational budget, we have
[|Var(Oar- -0) waP [|Var(0." -b)]ll
2 Pl-a.
wp/p P
Corollary (Efficieney) Given a decreasing step size W/. = ( where 0.5 < a < 1, ICSGLD is asymptotically more efficient than the single-chain CSGLD with an equivalent training cost.
In practice, slowly decreasing step sizes are often preferred in stochastic algorithms for better non-asymptotic performance (Benveniste et al., 1990).
5 EXPERIMENTS
5.1 LANDSCAPE EXPLORATION ON MNIST VIA THE SCALABLE RANDOM-FIELD FUNCTION
This section shows how the novel random-fille function (8) facilitates the exploration of multiple modes on the MNIST dataset', while the standard methods, such as stochastic gradient descent (SGD) and SGLD, only gel stuck in few local modes. To simplify the experiments, we choose large batch size of 2500 and only pick the first five classes, namely digits from 0 to 4. The learning rate is fixed to le-6 and the temperature is set to 0.1 1. We see from Figure 2(a) that both SGD and SGLD lead to fast decreasing losses. By contrast, ICSGLD yields fluctuating losses that traverse freely between high energy and low energy regions. As the particles stick in local regions, the penalty of ri-vititing these zones keeps increasing until negative learning rate is injected to encourage explorations.
(a) Training Loss
(b) SGD
(c) SGLD
(d) ICSGLD
1CGGL
SGLO
80
SGD
20
1003
100
100
200 Epochs
300
400
500
10
22
Figure 2: Visualization of mode exploration on a MNIST example based on different algorithms.
We conducted a singular value decomposition (SVD) based on the first two coordinates to visualize the trajectories: We first choose a domain that includes all the coordinates, then we recover the parameter based on the grid point and truncated values in other dimensions, and finally we fine-tune the parameters and present the approximate losses of the trajectories in Figure 2(b-d). We see SGD trajectories get stuck in a local region; SGLD exploits a larger region but is still quite limited in the exploration; ICSGLD, instead, first converges to local region and then escapes it once over- visits this region. This shows the strength of ICSGLD in the simulations of complex multi- distributions. More experimental details are presented in section D. of the supplementary material.
5.2 SIMULATIONS OF MULTI-MODAL DISTRIBUTIONS
This section shows the acceleration effect of ICSGLD via a group of simulation experiments for a multi-modal distribution. The baselines include popular Monte Carlo methods such as CSGLD, SGLD, cyclical SGLD (eyeSGLD), replica exchange SGLD (reSGLD), and the particle-based SVGD.
The target multi-modee density is presented in Figure 3(a). Figure 3(b-g) displays the empirical performance of all the testing methods: the vanilla SGLD with parallel chains (xP5) undoubtedly
$ 'The random-field function (Deng et al., 2020h) requires an extra perturbation term as discussed in section D4 in the supplementary material (Deng et al., 2020b); therefore it is not practically appealing in big data. Data augmentation implicitly leads a more concentrated posterior (Wenzel al., 2020; Aitchison, 2021).
7

---

Published as a conference paper at ICLR 2022
performs the worst in this example and fails quantify the weights of each mode correctly; the single-chain cycSGLD with 5 times of iterations (xT5) improves the performance but is still not accurate enough; reSGLD (xP5) and SVGD (xP5) have good performances, while the latter is quite costly in computations: ICSGLD (xP5) does not only traverse freely over the rugged energy landscape, but also yields the most accurate approximation to the ground truth distribution. By contrast, CSGLD (xT5) performs worse than ICSGLD and overestimates the weights on the left side. For the detailed setups, the study of convergence speed, and runtime analysis, we refer interested readers to section D.2 in the supplementary material.
(a) Truth
(b) SGLD (c) eyeSGLD
(d) SVGD
(c) reSGLD ( CSGLD
(g) ICSGLD
Figure 3: Empirical behavior simulation dataset. Figure 3(c) and 3(1) show the simulation based on a single chain with 5 times of iterations (xT5) and the others run 5 parallel chains (xP5).
5.3 DEEP CONTEXTUAL BANDITS ON MUSHROOM TASKS
This section evaluates ICSGLD the contextual bandit problem based the UCI Mushroom data set as in Riquelme et al. (2018). The mushrooms are assumed to arrive sequentially and the agent needs to take an action at each time step based on past feedbacks. Our goal to minimize the cumulative regret that measures the difference between the cumulative reward obtained by the proposed policy and optimal policy. We evaluate Thompson Sampling (TS) based on a variety of approximate inference methods for posterior sampling. We choose one c-greedy policy (EpsGreedy) based on the RMSProp optimizer with a decaying learning rate (Riquelme et al., 2018) as a baseline. Two variational methods, namely stochastic gradient descent with a constant learning rate (ConstSGD) (Mandt et al., 2017) and Monte Carlo Dropout (Dropout) (Gal & Ghahramani, 2016) are compared to upproximate the posterior distribution. For the sampling algorithmss, SGLD (pSGLD) (Li et al., 2016), preconditioned CSGLD (pCSGLD) (Deng et al., 2020b), and preconditioned ICSGLD (pICSGLD). Note that all the algorithms 4 parallel chains with average outputs P4) except that PCSGLD runs a single-chain with 4 times of computational budget (xT4). For more details, we refer readers to section D.3 in the supplementary material.
include preconditioned
Figure 4 shows that EpsGreedy xP4 tends to explore too much for long horizon as expected; ConstSGD x XP4 and Dropoutx P4 perform poorly in the beginning but eventually outperform Eps- Greedy x xP4 due to the inclusion of uncertainty for exploration, whereas the uncertainty seems to be inadequate due to the nature of variational inference. By contrast, PSGLDxP4 significantly outperforms the variational methods by considering preconditioners within an exact sampling framework (SGLD). As a unique algo- rithm that runs in a single-chain manner, PCSGLDx1F leads to the worst performance due to the inefficiency adapting parameters, fortunately, PCSGLD x T4 slightly outperform PSGLDXPP in the later phase with the help of the well-estimated self-adapting parameters. Nevertheless, PICSGLD xP4 propose optimize the shared self-adapting parameters at the same time, which in greatly contributes to the simulation of the posterior. As result, PICSGLD XP4 consistently shows the lowest regret exclud- ing the very early period. This shows the great superiority of the interaction mechanism in learning the self-adapting parameters for accelerating the simulations.
9000.
learning the self-
Canuube magque
0000
Epsc sP4
ConeriGD x*4
PC5GLD
pICSGLD
0
500
1000 Staps
1500
2000
Figure 4: Cumulative regret
the mushroom task.
5.4 UNCERTAINTY ESTIMATION
This section evaluates the qualify of our algorithm in uncertainty quantification. For model architec- tures, residual networks (ResNet) (He et al., 2016) and wide ResNet (WRN) (Zagoruyko Komodakis, 2016); we choose 20, 32, and 56-layer ResNets (denoted by ResNet20, et al.) and a WRN-16-8 network, 16-layer WRN that is 8 times wider than ResNet16. We train the models
8

---

Published as a conference paper at ICLR 2022
CIFARIO0, and report the test aceuracy (ACC) and test negative log- ligeliidd (NLL) based on trials with standard error. For the out-of-distriuutio prediction performance, we test the well-trained models in Brier scores (Brier) on the Street View House Numbers dataset (SVHN).
Due to the wide adoption of momentum stochastic gradient descent (M-SGD), we use stochastic gradient Hamiltonian Monte Carlo (SGHMC) (Chen al., 2014) the baseline sampling algorithm and denote the interacting contour SGHMC by ICSHMC. In addition, we include several high performing baselines, such as SGHMC with cyclical learning rates (cycSGHMC) (Zhang al., 2020b), SWAG based on cyclic learning rates of 10 cycles (cycSWAG) (Maddox et al., 2019) and variance-reduced replica exchange SGHMC (reSGHMC) (Deng al., 2021a). For fair comparison, ICSGLD also conducts variance reduction on the energy function to alleviate the bias. Moreover, a large Ç = x 10 selected, which only induces mild gradient multipliers ranging from -1 to penalize over-visited partitions. We don't 't include SVGD (Liu & Wang. 2016) and SPOS (Zhang et al., 2020a) for scalability reasons. batch size of 256 is selected. We run parallel processes (xP4) with 500 epochs for M-SGD, reSGHMC and ICSGHMC and run cycSGHMC and cycSWAG 2000 epochs (x T4) based on a single process with 10 cycles. Refer to section D.4 of the supplementary material for the detailed settings.
TABLE 1: UNCERTAINTY ESTIMATIONS ON CIFAR100 AND SVHN.
ResNer20 NLL
ResNet32
MODEL
(C)
Brier (%e)
(CC
NLL
Brier (%c)
cycSGHMCXT4
75.4120.10
8437130
2.91=0.13
77.9310.17
7658119
3.2910.13
cycSWAGXT4
75.46:0.111
84191:26
2.78=0.12
77.9110.15
76561:22
3.1910.14
M-SGD
76.01=0.12
8175#25
2.58=0.08
78,41:0.12
7501
2.77+0.15
reSGHMCXP4
76.15=0.16
8196:127
2.73=0.10
78.57:0.07
1554=15
3.0420.09
ICSGHMCXP4
76.34-0..5
8076+31
2.54 0.14
78.7240.16
7406+22
2.76+0.15
ResNet56
WRN- 16-8
MODEL
ACC (%)
NLL
Brier (%c)
ACC (%)
NLL.
Brier (oo)
cycSGHMCxT4
81.23=0.19
6770+:59
3.18-0.08
82.98+0.033
6384+11
2.17+0.05
eycSWAGXT4
81.1410.11
67441555
3.06
83.0510.04
6359114
2.0410.07
M-SGDXP4
81.0330.114
6847+-22
82.57+0.07
6821#21
1.77+00.06
reSGHMCXP4
81.1120.16
6915140
2.92=0.12
82.7210.08
6452=19
1.92:1.00
ICSGHMCXP4
81.5120.18
6630138
2.88=0.09
83.1210.10
6338136
1.83
Table shows that the vanilla ensemble results via M-SGDxP4 surprisingly outperform cycSGHMCXT4 and cycSWAGxT4 on medium models, such ResNet20 and ResNet32, and show very good performance on the out-of-distribution samples in Brier scores. We suspect that the parallel implementation (x (xP4) provides isolated initializations with less correlated samples; by contrast, cycSGHMCXT4 T4 and eycS risk stay near the original region. reSGHMC P4 shows a remarkable performance overall, but demonstrates large variance occasionally; this indicates the insufficiency of the swaps when multiple processes are included. When it comes to testing WRN-16-8, cycSWAG> T4 shows a marvelous result and a large improvement compared to the other baselines. We conjecture that cycS is more independent of hyperparameter tuning, thus leading to better performance in larger models. We don't ton't report CSGHMCx x P4 since it becomes quite unstable during the training of ResNet56 and WRN-16-8 models and causes mediocre results. As to ICSGHMC PP,, it consistently performs remarkable in both ACC and NLL and performs comparable to M-SGD x PP in Brier scores.
explore the energy landscape contiguously, implying a
Code is available at github..cmmWWyyyDDD//tttteetttte Conne-Crontttttute Stochastic-Gratee Langgevn-- Dynnmi-sDDamiss...
6 CONCLUSION
We have proposed the ICSGLD as an efficient algorithm for sampling from distributions with a complex energy landscape, and shown theoretically that ICSGLD is indeed more efficient than the single-chain CSGLD for a slowly decreasing step size. To our best knowledge, this is the first interacting importance sampling algorithm that adapts big data problems without scalability concerns. ICSGLD has been compared with numerous state-of-the baselines for various tasks, whose remarkable results indicate its promising future in big data applications.
"The Brier score measures the mean squared error hetween the predictive and actual prohabilitics.

---

Published as a conference paper at ICLR 2022
ACKNOWLEDGMENT
Liang's research was supported in part by the grants DMS-2015498, RO1-GM117597 and R01- 26096 Lin acknowledges the support from NSF (DMSS-5550022, DMS-2053746, and DMS- 2134209), BNL Subcontract 382247, and DE-SC0021
REFERENCES
Sungjin Ahn, Anoop Korattikara, and Max Welling. Bayesian Posterior Sampling via Stochastic Gradient Fisher Scoring. In Proc. of the International Conference on Machine Learning (ICML), 2012.
Sungjin Ahn, Babak Shahbaba, and Max Welling. Distributed Stochastic Gradient MCMC. In Proc. of the International Conference Machine Leurning (ICML), 2014.
Laurence Aitchison. A Statistical Theory of Cold Posteriors in Deep Neural Networks. In Proc. of the International Conference on Learning Representation (ICLR), 2021.
C. Andrieu, E. Moulines, and P. Priouret. Stability of Stochastic Approximation under Verifiable Conditions. SIAM J. Control Optim., 44(1):283-312, 2005.
Christophe Andrieu, Arnaud Doucet, and Roman Holenstein. Particle Markov Chain Monte Carlo Methods. Journal of the Royal Statistical Society: Series (Statistical Methodology), 72(3), 2010.
Albert Benveniste, Michael Métivier, and Pierre Priouret. Adaptive Algorithms and Stochastic Approximations. Berlin: Springer, 1990.
Elmar Bittner, Andreas Nussbaumer, and Wolfhard Janke. Make Life Simple: Unleash the Full Power of the Parallel Tempering Algorithm. Physical Review Letters, 101:136603-130003, 2008
Changyou Chen, Nan Ding, and Lawrence Carin. On the Convergence of Stochastic Gradient MCMC Algorithms with High-order Integrators. In Advances in Neural Information Processing Systems (NeurIPS), pp. 22782227228 2286,22266222 2015.
Changyou Chen, Nan Ding, Chunyuan Li, Yizhe Zhang, and Lawrence Carin. Stochastic Gradient MCMC with Stale Gradients. In Advances in Neural Information Processing Systems (NeurIPS), 2016.
Tianqi Chen, Emily B. Fox, and Carlos Guestrin. Stochastic Gradient Hamiltonian Monte Carlo. In Proc. of the International Conference on Machine Learning (ICML), 2014.
Wei Deng, Qi Feng, Liyao Gao, Faming Liang, and Guang Lin. Non-Convex Learning via Replica Exchange Stochastic Gradient MCMC. In Proc. of the International Conference on Machine Learning (ICML), 2020a.
Wei Deng, Guang Lin, and Faming Liang. Contour Stochastic Gradient Langevin Dynamics Algorithm for Simulations of Multi-modal Distributions. In Advances in Neural Information Processing Systems (NeurIPS), 2020b.
Wei Deng, Qi Feng, Georgios Karagiannis, Guang Lin, and Faming Liang. Accelerating Conver- gence of Replica Exchange Stochastic Gradient MCMC via Variance Reduction. In Proc. of the International Conference Learning Representation (ICLR), 2021a.
Wei Deng, Yi-An Ma, Zhao Song, Qian Zhang, and Guang Lin. Convergence of Federated 2021b.
Averaging Langevin Dynamics.
Nan Ding, Youhan Fang, Ryan Babbush, Changyou Chen, Robert D. Skeel, and Hartmut Neven. Bayesian Sampling using Stochastic Gradient Thermostats. In Advances in Neural Information Processing Systems (NeurIPS), pp. 2203-3211, 2014.
Arnaud Doucet. Nando de Freitas, and Neil Gordon. Sequential Monte Carlo Methods in Practice. Springer Science & Business Media, 2001.
10

---

Published as a conference paper at ICLR 2022
Alain Durmus and Éric Moulines. Non-asymptotic Convergence Analysis for the Unadjusted Langevin Algorithm. Annals of Applied Probability,
Ja Earl and Michael W. Deem. Parallel Tempering: Theory, Applications, and New Perspectives.
Phys. Chem. Chem. Phys., 7:3910-3911 2005.
Murat A Erdogdu, Lester Mackey, and Ohad Shamir. Global Non-convex Optimization with Dis- cretized Diffusions. Advances in Neural Information Processing Systems (NeurIPS), 2018.
G. Fort, Moulines, and Priouret. Convergence of Adaptive and Interacting Markov Chain Monte Carlo Algorithms. Annals of Statistics, 39:3262-3288 2011.
G. Fort, B. Jourdain, E. Kuhn, T. Lelièvre, and G. Stoltz. Convergence of the Wang-Landau Algorithm. Math. Comput., (4(295):2297-2333 2015.
Futoshi Futami, Issei Sato, and Masashi Sugiyama. Accelerating the Diffusion-based Ensemble Sampling by Non-reversible Dynamics. In Proc. of the International Conference Machine Learning (ICML), 2020.
Yarin Gal and Zoubin Ghahramani. Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning. In Proc. of the International Conference Machine Learning (ICML), 2016.
Charles Geyer. Markov Chain Monte Carlo Maximum Likelihood. Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interfac, pp. 156-163, 1991.
Neil J Gordon, David J Salmond, and Adrian FM Smith. Novel Approach to Conlinear/Non-Gaussie Bayesian State Estimation. IEE Proceedings F (Radar and Signal Processing), 140(2), 1993.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep Residual Learning for Image Recognition. In The IEEE Conference Computer Vision and Pattern Recognition (CVPR), 2016.
Pavel Izmailov, Dmitry Podoprikhin, Timur Garipov, Dmitry Vetrov, and Andrew Gordon Wilson Averaging Weights Leads to Wider Optima and Better Generalization. In Proc. of the Conference on Uncertainty in Artificial Intelligence (UAI), 2018.
K. Jarrett, K. Kavukcuoglu, M. Ranzato, and Y. LeCun. What is the Best Multi-stage Architecture for Object Recognition? In Proc. of the International Conference on Computer Vision (ICCV), pp. 1146-21532 September 2009.
Helmut G Katzgraber, Simon Trebst, David A Huse, and Matthias Troyer. Feedback- -eabbaackOOtimizee Parallel Tempering Monte Carlo. Journal of Statistical Mechanics: Theory and Experiment, pp. p. P03018, 2008
Chunyuan Li, Changyou Chen, David Carlson, and Lawrence Carin. Preconditioned Stochastic Gradient Langevin Dynamics for Deep Neural Networks. In Proc. of the National Conference on Artificial Intelligence (AAAT), pp. 1788-1788 1794, 2016.
Chunyuan Li, Changyou Chen, Yunchen Pu, Ricardo Henao, and Lawrence Carin. Communication- Efficient Stochastic Gradient MCMC for Neural Networks. In Proc. of the National Conference on Artificial Intelligence (AAAI), 2019a.
Xiang Li, Kaixuan Huang, Wenhao Yang, Shusen Wang, and Zhihua Zhang, On the Convergence of FedAvg on Non-IID Data. In Proc. of the International Conference on Learning Representation (ICLR), 2020.
Xuechen Li, Denny Wu, Lester Mackey, and Murat A. Erdogdu. Stochastic Runge- Kunge-Kutta Accelerates Langevin Monte Carlo and Beyond. In Advances in Neural Information Processing Systems (NeurIPS), pp. 7746-7758, 2019b.
Faming Liang, Chuanhai Liu, and Raymond J. Carroll. Stochastic Approximation in Monte Carlo Computation. Journal of the American Statistical Association, 102:305-320, 2007.
Qiang Liu and Dilin Wang. Stein Variational Gradient Descent: A General Purpose Bayesian Inference Algorithm. In Advances in Neural Information Processing Systems (NeurIPS), 2016.

---

Published as a conference paper at ICLR 2022
Wesley Maddox, Timur Garipov, Pavel Izmailov, Dmitry Vetrov, and Andrew Gordon Wilson. A Simple Baseline for Bayesian Uncertainty in Deep Learning. Advances in Neural Information Processing Systems (NeurIPS), 2019.
Stephan Mandt, Matthew D. Hoffman, and David M. Blei. Stochastic Gradient Descent as Approxi- mate Bayesian Inference. Journal of Machine Learning Research, 18:1-35, 2017.
J.C. Mattingly, A.M. Stuartb, and DJJ. Highame. Ergodicity for SDEs and Approximations: Locally Lipschitz Vector Fields and Degenerate Noise. Stochastic Processes and their Applications, 101: 185-232, 2002.
Jonathan C. Mattingly, Andrew M. Stuart, and M.V. Tretyakov. Convergence of Numerical Time- Averaging and Stationary Measures via Poisson Equations. SIAM Journal on Numerical Analysis, 48:552-577, 2010.
Mariane Pelletier. Weak Convergence Rates for Stochastic Approximation with Application to Multiple Targets and Simulated Annealing, Annals of Applied Probability, 8:10-44, 1998.
Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky. Non-convex Learning via Stochastic Gradient Langevin Dynamics: a Nonasymptotic Analysis. In Proc. of Conference on Learning ''heory (COLT), June 2017.
Tom Rainforth, Christian A. Naesseth, Fredrik Lindsten, Brooks Paige, Hun-Willem van de Meent, Arnaud Doucet, and Frank Wood. Interacting Particle Markov Chain Monte Carlo. In Proc. of the International Conference on Machine Learning (ICML), 2016.
Carlos Riquelme, George Tucker, and Jasper Snoek. Deep Bayesian Bandits Showdown. In Proc. of the International Conference Learning Representation (ICLR), 2018.
Herbert Robbins and Sutton Monro. A Stochastic Approximation Method. Annals of Mathematical Statistics, 22:400-407, 1951.
Gareth O. Roberts and Richard L. Tweedie. Exponential Convergence of Langevin Distributions and Their Discrete Approximations. Bernoulli, 2(4):341-363 1996.
Issci Sato and Hiroshi Nakagawa. Approximation Analysis of Stochastic Gradient Langevin Dynamics by Using Fokker-Planck Equation and Ito Process. In Proc. of the International Conference Machine Learning (ICML), 2014.
Qifan Song, Minggi Wu, and Faming Liang, Weak Convergence Rates of Population versus Single- Chain Stochastic Approximation MCMC Algorithms. Advances in Applied Probability, 46: 1059-1083, 2014.
Robert H. Swendsen and Jian-Sheng Wang. Replica Monte Carlo Simulation of Spin-Glasses. Physical Review Letters, 57:2607-2609, 1986.
Saifuddin Syed, Alexandre Bouchard-Côré, George Deligiannidis, and Arnaud Doucet. Non-
Reversible Parallel Tempering: Scalable Highly Parallel MCMC scheme. Journal Royal Statistical Society, Series B, 2021.
Yee Whye Teh, Alexandre Thiéry, and Sebastian Vollmer. Consistency and Fluctuations for Stochastic Gradient Langevin Dynamics. Journal of Machine Learning Research, 17:1-33, 2016.
Eric Vanden Eijndenn Introduction to Regular Perturbation Theory. Slides, 2001. URL https: pdf..
//cims. nyu. ecuu/~eve2reegpeee
Sebastian J. Vollmer, Konstantinos C. Zygalakis, and Yee Whye Teh. Exploration of the (Non-) Asymptotic Bias and Variance of Stochastic Gradient Langevin Dynamics, Journal of Machine
Learning Research, 17(159):1-43 18,, 2016.
Fugao Wang and David P. Landau. Efficient, Multiple-range Random Walk Algorithm to Calculate the Density of States. Physical Review Letters, 86:2050-3, 2001.
T. Weinhart, A. Singh, and A.R. Thornton. Perturbation Theory & Stability Analysis. Slides, 2010.
12

---

Published as a conference paper at ICLR 2022
Max Welling and Yee Whye Teh. Bayesian Learning via Stochastic Gradient Langevin Dynamics, In Proc. of the International Conference on Machine Learning (ICML), pp. 681-688 2011.
Florian Wenzel, Kevin Roth, Bastiaan S. Veeling, Jakub Swiatkowski, Linh Tran, Stephan Mandt, Jasper Snoek, Tim Salimans, Rodolphe Jenatton, and Sebastian Nowozin. How Good is the Bayes Posterior in Deep Neural Networks Really? In Proc. of the International Conference Machine Learning (ICML), 2020.
Pan Xu, Jinghui Chen, Difan Zou, and Quanquan Gu. Global Convergence of Langevin Dynamics Based Algorithms for Nonconvex Optimization. In Advances Neural Information Processing Systems (NeurIPS), 2018.
Sergey Zagoruyko and Nikos Komodakis. Wide Residual Networks. In Proceedings of the British Machine Vision Conference (BMVC), pp. 87.1- September 2016.
Jianyi Zhang, Ruiyi Zhang, Lawrence Carin, and Changyou Chen. Stochastic Particle- Pnttize--ptimizzation Sampling and the Non-Asymptotic Convergence Theory. In Proceedings of the International Workshop on Artificial Intelligence and Statistics, 2020a.
Ruqi Zhang, Chunyuan Li, Jianyi Zhang, Changyou Chen, and Andrew Gordon Wilson Cyclical Stochastic Gradient MCMC for Bayesian Deep Learning. In Proc. of the International Conference on Learning Representation (ICLR), 2020b.
Yuchen Zhang, Percy Liang, and Moses Charikar. Hitting Time Analysis of Stochastic Gradient Langevin Dynamics. In Proc. of Conference on Learning Theory (COLT), pp. 1980-2022, 2017.
Zhun Zhong, Liang Zheng, Guoliang Kang, Shaozi Li, and Yi Yang. Random Erasing Data Augmen- tation. ArXiv e-prints, 2017
13

---

Published as a conference paper at ICLR 2022
We summarize the supplementary material follows: Section A provides the preliminary knowledge for stochastic approximation; Section B shows a local stability condition that adapts to high losses; Section C proves the main asymptotic normality for the stochastic approximation process, which naturally yields the conclusion that interacting contour stochastic gradient Langevin dynamics (ICSGLD) is more efficient than the analogous single chain based slowly decreasing step sizes; Section D details the experimental settings.
A PRELIMINARIES
A.1 STOCHASTIC APPROXIMATION
Given random-field function O x), the stochastic approximation algorithm (Benveniste et al., 1990) proposes to solve the ween-field equation h()) = 0 in the analysis of adaptive algorithms
/ (0, x x) = 0, x
h(0)
where x E X C 12, C IR (x(x) is a distribution that depends on the self-adapting parameter 0. Given the transition kernel HT (z, A) for any Borel subset A C X, the algorithm be written as follows
(1) Simulate X Ho (Xx:)) +), which yields the invariant distribution wa (-),
(2) Optimize 0 = 0 + wk+1 (0,, Xx+1).
Compared with the standard Robrins-Monro algorithm (Robbins & Monro, 1951), the algorithm proposes to simulate x from a transition kernel (= instead of the distribution wa(-) directly. In other words,, (O,, x - h(Ox) is not a Martingale but rather Markov state-dependent noise.
A.2 POISSON'S EQUATION
In the stochastic approximation algorithm, the sequence of {(xx, 0x) on the product space X' x ② is generated, which is inhomogeneous Markov chain and requires the tool of the Poisson's equation to study the convergence
( Houe (x)) = (0, x) (()).
where (-) is a function on X. The solution (x) to the Poisson's 's equation exists and is formulated in the following form when the above series converges:
mx(x)
Σ En (H((0 h(0)).
h(0)) dy). To ensure such convergence, where (H, h(0)) ( -
Benveniste et al. (1990) made the following regularity conditions on the solution (a(-) of the
Poisson's equation:
There exist a functiin V X' > [1, oo) and a positive constant > O such that 00, 0' E O, have
[To x) I CV()), (x) - Hor (x)l < C e (V(), ( S 00, (12)
where a common choice for the Lyapunov function to set V(x) = + |||| (Teh et al., 2016; Vollmer al., 2016).
A.3 GAUSSIAN DIFFUSIONS
Consider stochastic linear differential equation
dU, = Udd + (O,
do (13)
where U is m-dimensional random vector, he
d 1((0).
R(0)
2o= (H xx), H(0, xo)) is positive definite matrix that depends 0()) W m
14

---

Published as a conference paper at ICLR 2022
is standard Brownian motion. Given large enough such that 0, converges fixed point 0. sufficiently fast, we may write the diffusion associated with (q.(13) as follow
U ~ e the (0. + 1." ((-ss)hoo(O.) o R(O.)a...
(14)
Suppose that the matrix ho.) is negative definite, then U, converges in distribution Gaussian variable
E|U ee (t.)
Var(U / chholê.) R ddu
The main goal of this supplementary file is to study the Gaussian approximation of the process
We 1112 (0, 0.) the solution Eq.(14) for proper step size wi. Thereafter, the advantage of interacting mechanisms can naturally derived.
B STABILITY AND CONVERGENCE ANALYSIS
As required by the algorithm, we update contour stochastic gradient Langevin dynamics (CSGLD) simultaneously. For the notations, we denote the particle of the p-th chain at iteration by x E
(p)
X
C
T
Rd and the joint state of the P parallel particles at iteration by x
X
E
X P C RdP, We also denote the learning rate and step size at iteration & by and Wk. respectively. We denote by NO, lap) a standard d'P-dimensional Gaussian vector and denote by ( positive hyperparameter.
B. ICSGLD ALGORITHM
First, we introduce the interacting contour stochastic gradient Langevin dynamics (ICSGLD) with P parallel chains:
(1) Simulate Xk+1 X
P
8
V
x V 02) + (0,
(S))
(2) Optimize θ 0g +
Wk+ H
&P
(S2)
x+11
where V x
(x 0), VxI(x(?),0) 0). ..,(I((()()))) V.ī(x. 0) is the
stochastic adaptive gradient given by
VxI(x, 0)
1+ CT Au
(1gg
(xū( (15)
gradient multiplier
In particular, the interacting random-field function is written as
H x P
Σ
H (2).)
(16)
p=1
where each random-felk function H(0, x) (0,x).... (0,x))) follows
H,(0,x) (x)) (1,=J](x)) (()). i 1.2...., m.
(17)
Here (x(x) denotes the index E (1,2,3. m) such that u41 (x) S u, for set of
energy partitions {u and ((x) = Σien (x) where U; denotes the negative log of a posterior based on single data point and B denotes mini-batch of data of size Note that the stochastic
15

---

Published as a conference paper at ICLR 2022
energy estimator (x) results in biased estimation for the partition index (x) due to non-linear transformation. To avoid such a bias asymptotically with respect to the learning rate ek, we may consider a variance-reduced energy estimator UVR(x) following Deng et al. (2021a)
(xR(X)
Σ (Ui(x) BE
- (x (x+1)1)) (xxx1)))))
U3
X x
(18)
where the control variate x is updated every q iterations.
Compared with the na've parallelism of CSGLD, key feature of the ICSGLD algorithm lies in the joint estimation of the interacting random-field function (00 x P) in (q, for the same mean-field function h()).
B.1.1 DISCUSSIONS ON THE HYPERPARAMETERS
The most important hyperparameter is Ç. A fine-tuned Ç usually leads small slightly negative learning rate in low energy regions to avoid local-trap problems. Theoretically, affects the L convergence rate hidden in the big-O notation in Lemma 3.
The other hyperparameters can be easily tuned. For example, the ResNet models yields the full loss ranging from 10,000 to 60,000 after warm-ups, we thus partition the sample space according the energy into 200 subregions equally without tuning; since the optimization of SA is nearly convex, tuning {wx} is much easier than tuning (ex for non-convex learning.
B.1.2 DISCUSSIONS ON DISTRIBUTED COMPUTING AND COMMUNICATION COST
In shared-memory memory settings, the implementation is trivial and the details omitted. In distributed memurr settings: 0x+1 updated by the central node as follows:
The p-th worker conducts the sampling step (Si) and sends the indices J /((e()))) 's to the central
node;
The central node aggregates the indices from all worker and updates 0g based on (S2);
The central node sends 0x+1 back to each worker.
We emphasize that don't 't communicate the model parameters Rd, but rather share the self- adapting parameter E RRm where m < d. For example, WRN-16-8 has 11 M parameters (40 MB), while A can be set to dimension 200 of size 4 KB; hence, the communication cost is not big issue. Moreover, the theoretical advantage still holds if the communication frequency is slightly reduced.
B.1.3 SCALABILITY TO BIG DATA
Recall that the adaptive sampler follows that
N
log (x(()(x)) - log 0u (() - 1)
Vxū(xx) + (2TEL-1WHHH
X++1 = X2 - Ek+1
(
n
Av
gradient multiplier
The key to the success of (I)CSGLD is to generate sufficiently strong bouncy moves (negative gradient multiplier) escape local traps. To this end, can be tuned generate proper bouncy moves.
Take the CIFARI00 experiments for example:
the self-adjusting mechanism fails if the gradient multiplier uniformly "equals" to and too small value of ( could lead to this issue;
the self-adjusting mechanism works only if we choose a large enough Ç such as 3e6 to generate (desired) negative gradient multiplier in over-visited regions.
16

---

Published as a conference paper at ICLR 2022
However, when set =3e6, the original stochastic approximation (SA) update proposed in (Deng et al., 2020b) follows that
0k+1 1 (i) (i) + wh+1 (x (1==/((0--1)-a
essentially 0 for (36-
Since (()) < 1 for any E (1.... ,,), 0()) is essentially O for such a large ç, which means that the original SA fails to optimize when is large. Therefore, the limited choices of inevitably limits the scalability to big data problems. Our newly proposed SA scheme
Oa+ (i) 0x (x) Wk+1 Oz (J (1==Jg(Xe-1) -0x(1)
independent of
more independent of Ç and proposes to converge to a much smoother equilibrium 0%/ instead of Ox0. where 0z (i) fx. π (x))d x fxx e dr is the energy PDF. As such, despite the linear
stability is sacrificed, the resulting algorithm is scalable. For example, estimating -10,000x is numerically much easier than e 100,000 for a large such as 10, 000, where 10, 000 can be induced by the high losses training deep neural networks in big data.
B.2 ASSUMPTIONS
A long-standing problem for stochastic approximation is the difficulty in establishing the stability property and a practical remedy for this problem is to study e on a fixed compact set.
Assumption A1 (Compactness) The space is compact and infe 0(i) O for any i {1, 2, mL.......m..
For weaker assumptions, we refer readers to Theorem 3.2 (Fort et al., 2015), where recurrence prop- erty can be proved for the Metropolis-basee Wang-Landan algorithm, which eventually established that the estimates return to a desired compact set often enough.
Next, we lay out the smoothness assumption, which is standard in the convergence analysis of SGLD, see e.g. Mattingly et al. (2010), Raginsky et al. (2017) and Xu al. (2018).
Assumption A2 (Smoothness) (ax) is M-smooth when there exists positive constant M that satisfies vx, x' E X,
< M x -XII
(19)
In addition, we assume the dissipativity condition to ensure that the geometric ergodicity of the dynamical system holds. This assumption is also crucial for verifying the solution properties of the solution of Poisson's equation. Similar assumptions have been made in Mattingly et al. (2010); Raginsky al. (2017) and Xu et al. (2018).
Assumption A3 (Dissipativity) There exist constants m and b 2 that satisfies Vx E X' and any 0 θ
(VxL(x,0). x) < m
(20)
To further establish a bounded second moment on x E X' with respect to a proper Lyapunov function (x), we impose the following conditions the gradient noise:
Assumption A4 (Gradient noise) The stochastic gradient based on mini-batch settings is an unbi- ased estimator such that
E(VxÜ(xx) - vxU(xx)] =0: O;
furthermore, for some positive constants M and B, we have
s 12 B², E[IV2ū((x)) - (xU(xx) Vx.
where E[|] acts the distribution of the noise in the stochastic gradient (xx).
17

---

Published as a conference paper at ICLR 2022
B.3 LOCAL STABILITY VIA THE SCALABLE RANDOM- FUNCTION
Now, we are ready present our first result. Lemma 3 establishes a local stability condition for the non-linear mean-field system of ICSGLD, which implies potential convergence of 0μ to a unique fixed point that adapts to wide energy range under mild assumptions.
Lemma (Local stability, restatement of Lemma 1) Assume Assumptions A AI-A4 hold. Given any small enough learning rate e, a large enough m and batch size n,, and any 0 E O where O a small neighhorhood of 0, that contains 0,, we have (h(0). A 0-) S -- 1+1 where 0 = 0. (e)), (Sx, m(x)
(Sxm, =(x)dx)
E =0 2 (supx Var(5n + c+ + and 0,
> (Jx, x(x)dx)
DT (Ix, =(x)dx)
e info min; A (1 (()) > 0, Z,,00)) is defined below (28)) and &n (x) denotes the noise in the energy estimator Ü (x) of batch size and Var(.) denotes the variance.
Proof The aundo--ieed function (0, x) (()((x)) (1:))))111)1) (x)=J,(x)) (000)))-88)))) -8()) based on the stochastic energy estimator (x) yields a biased estimator of (0, x) (x))) (Li=J(x) (()) for any i E {1,2,... m) based on the exact energy partition function (-). By Lemma.4, we know that the bias caused by the stochastic energy is of order ((Var(Gxx)))
Now we compute the mean-field function I()) based on the measure wx(x) simulated from SGLD:
(0) -s (x)xxxxxxxxxxxx = / H: (x),)xx))x(xxxxxx++x)))))))))))) (Var(E.(x)))
= (0,
me (x)
w (x) + we (x) wa (x) (x dx O (Var(En(x))),
pappicoz-wise approtimation
In:numerical discretization
(21)
where 5y the invariant measure simulated via SGLD that approximates Wa, (x) Do (x) and wio (x) are two invariant measures that follow ww (x)
(x) and We, (x) x ((((())))
s(x)
x
(x))
(u) and (u) are piecewise continuous and constant functions, respectively
Vo(u)= = Σ
(o(k - 1)e(lox - log A(k
luz-rcusui Vo(u) = Σ 0(k) ([k)lug-i <usu
(22)
(i For the first term I, we have
[ (0.00))0,, (x)dx
1
/ ((J(x)) (1==)6x) 0(i)) T(x)
05 (x)))x)) dx
Zc+1.0
1
m
T(x)
(Lizk - -O())
()-)((k) dx
Z 15+0
k=1
Xx
(23)
1
π T(x)
1 1k=id 0(1) Σ
(S-1 dx
T(x)
Zc+i.
Nr.
(K-)
k= Nr.
Sxx nT (x)dx
1
0(1)
Z
05-1(i)
where Zct1,0
Ek=1
Ix, =xx)dx (<(k))
denotes the normalizing constant of a (x).
18

---

Published as a conference paper at ICLR 2022
Ix. =(x)ax (S-1(k)
The solution 0, that solves
O(k) 0 for any & {,,2, ... ..m satisfies 0, (k)
Ixg =(x))xx
Combining the definition of Zç.o. = -k== k=
Sx, (x]dxx (R) we have
7c.0.
1xx (x))xx
Sx. (xx)
Zç.0.
Σ
1--1((k) (k)
k=1
k=1
Ixx, m(x)dxx
25.0.
5-1 Z
Σ
SNx x(x)xx
Z5,o.
>
172
x(x))xx)
s.O
5-1
Xr
k=1
(Ixs (xx)
R=1
<x== (Ixr, x(x)) #)
which leads to Zç.o. =
In other words, the mean-field system without
(Jx, m(x))xx)
perturbations yields unique solution (z)
for any i E (1,2...., .m).
> (Ix. =(xx)dx))
(ii) For the second term 12, we have
/
(11)
1
(,, (x) Fea (x)))xx 0
(24)
m
where the result follows from the boundedness of H x) in (A1) and Lemma B4 (Deng al.,
2020b).
(iii) For the last term 13, following Theorem of Sato & Nakagawa (2014), we have for any fixed 0,
- H0) (- (x) + (x(x)) dx O
(25)
Plugging (23), (24) and (25) into (2I,, we have
Sx. π (x)dx
(0) = 7-1
(0,
0(1)
c+1,0
05
(i)
7z.0.
( 0 8C-1(1)
Sx. x(x)dx
Z5.00
= -
1+1.0 0 (i)
22.0.
+
Zg.o.
(c())
Zç.0.
(26)
Zg.o. =Z c+1.0 2
8 -
+ (i( - -0(1)CC))
Z5.
where 3,(0) is bounded term such that B
o (Var(En,x)) + H Sx, m(x)xx
O (supx (x)()((x)) + 1h) and Co = 2 when { = 1, Co = for any O E θ, which suggests that the stability condition doesn't rely the initialization of 0: however, when < # Co * when 0 * O., we see that (0) x (1)s (0())Co) + perturbations proper initialization of 0 E e.
By the definition of 7,.0 = =m (c-((b)
Zç.0.
non-linear mean-field system and requires
For any 0 E θ
2 being close enough to 0., there exists Lipschitz constant Lg = < cc. By Co. = 1, 0(i) < 1, and mean value theorem for some
ICe, -Cal
sup 0cō (1)
0(1) E [0(i), (i)), we have
(05())
- -(0()))))-11110000) -
=(((0))))--111()) - ( O(i)Ca. - O(i)Col
S ((00))))--10000)) +() - -Co
3 ((0())))-- Lg)00,00 0(1),
(27)
19

---

Published as a conference paper at ICLR 2022
Combining (26 and (q.( we have
Žç.o.
BC (1) (1)(1) - (0(1)C0)
h₂(0) = Zc+1, 0c (i-
cB.
Zç.O.
(28)
(,,000) (E); (1) 0()).
where (.,(1) = 20,00
B:(0) is some bounded term such that (:(
S
3(0000--(00 Co
2g,
1 = supisn smp 0.(0)-001
CCo. -Ca Cc.
Z5.0.
Next, we apply the perturbation theory to solve the ODE system with small disturbances ( (Weinhart et al., 2010) and obtain the equilibrium Q.,
where eB((0.) 0. - 0. 0, to the mean-field equation (0) such that
(,))) S. (E8;(0) + (1) [())
(T,b0) (E8((0) =B. + (0.) 0.(01) 0(i)
(29)
2 [o(=)(8() - - ( 0(1)
2-1 c,0(1) (1 O(s)) (0.(1) 0(1)
where smoothness condition clearly holds for the 3--) function f. Given a positive definite
Lyapunov function V(0)= = _110_- - 0||; the mean-field system h())) (c) + 0, - 0)
Z (1 ( 0) for E {1,2, ... mm enjoys the following property
(M()), VV(0)) (N(0), 0 0-))
I min 25.b(s) -1-0(e))10-0.11 0.1²
S 0 - 02
where info min; Z,,000) (1- O())) > 0 given the compactness assumption Al and small enough O (sup, Vx)( + + ..)
Remark The newly proposed random-field function Eq.(17) may sacrifice the global stability by including an approximately linear mean-field system (q. instead of a linear stable system (see formula (15) in Deng et al. (2020b)). The advantage, however, is that such mechanism facilitates the estimation of 0.. We emphasize that the original energy probability in each partition {Sx, x(x)dxx (Deng et al., 2020b) may be very difficult to estimate for big data problems. By contrast, the estimation of (Jxx, (x))xx) $ }m
k=1 becomes much easier given proper Ç > 0.
Technical lemmas
Lemma 4 The stochastic energy estimator Úxx leads to a controllable bias in the random-field function.
[E, x)] H, O (Var(E.(())))
where the expectation B[]] is taken with respect the random noise in the stochastic energy estimator of U-).
Proof
Denote the noise in the stochastic energy estimator by {(x), such that U..) = U(-) +
(()) Recall that H,, = (x)) (1==)g(x) - 0()) and (x)) E (1,2,---, ..,,,..,,,,. m) satisfies
"A small change of 0 won't significantly affect the perturbations caused by (24), (25) and Var(E.
20

---

Published as a conference paper at ICLR 2022
U (x) N (x) S (x) for set of energy partitions {u We can interpret (0, x) as non-linear transformation & that maps (x) to (0, 1). Similarly, H; x) maps U(x) to (0, 1). In what follows, the bias of random-field function is upper bounded as follows
/ ((U(x) + {(x)) +(U(x))d((((())
(E, H, x)|
/ ((x))(((())) + {(x)2
l ("(u)d((Exx))
2
/ (x(x) (U))))((((()))d((((())
*"(u)
f (x(x)
2
S O (x))
where v(((x)) is the probability measure associated with {(x); the second equality follows from Taylor expansion for some energy and the third equality follows because the stochastic energy estimator is unbiased; &' (x)((x)) = O (0((x)))[([((xx]]) is clearly bounded due to the definition of 0: similar conclusion also applies to (").). The last inequality easily follows by applying Cauchy Schwarz inequality.
B.4 CONVERGENCE OF THE SELF-ADAPTING PARAMETERS
The following is restatement of Lemma 3.2 of Raginsky et al. (2017), which holds for any i in the compact space C.
Lemma (Uniform 1.? bounds) Assume Assumptions AI, A3 and A4 hold. have a bounded second moment supw21 [|||x []] < given small enough learning rate.
The following lemma justifies the regularity properties of Poisson's 's equation, which is crucial in controlling the perturbations through the stochastic approximation process. The first version was proposed in Lemma B2 of Deng et al. (2020b). Now we give more detailed proof by utilizing Lyapunov function (xx)
x² and Lemma 5.
Lemma 6 (Solution of Poisson's 's equation) Assume that Assumptions Al-A4 hold. There is solu- tion (a() on X to the Poisson's 's equation
(x [o) H1 ((0).
(30)
Furthermore, there exists constant C such that for all 0,0' +
E|ll II] < C, E [[ll - 11. (x] < C 0'|
(31)
The existence and the regularity property of Poisson's equation can be used to control the perturbations. The key of the proof lies in verifying drift conditions proposed in Section 6 of Andrieu et al. (2005).
Proof
(DRI) By the smoothness assumption A2, we have that Uxx is continuously differentiable almost everywhere. By the dissipative assumption A3 and Theorem 2.1 (Roberts & Tweedie, 1996), we show that the discrete dynamics system is irreducible and aperiodic. Now consider a Lyapunov function V = + ||x and any compact subset I. c O the drift conditions are verified as follows:
(DRI1) Given small enough learning rates {EN the smoothness assumption A2, and the dissipa- tive assumption A3, applying Corollary 7.5 (Mattingly et al., 2002) yields the minorization condition for the CSGLD algorithm, i.e. there exists 17 O, measure v, and a set such that v(C) 1. Moreover, we have
Poex A) (A) VAEXXXEC.
denotes the transition kernel based on CS- V [x]
where Po(xyy) 2 (1
E
GLD with the parameter 0 E K. and learning rate e, in addition, the expectation is taken the
21

---

Published as a conference paper at ICLR 2022
adaptive gradient Vxī(x, 0) in (15). Using Assumption A1-A4, can prove the uniform L2 upper bound by following Lemma 3.2 (Raginsky et al., 2017). Further, by Theorem 7.2 (Mattingly et al., 2002), there exist E (0, 1) and 3 2 0 such that
PoexV(x)) S (V(x) + B
(II)
B, it yields that Consider Lyapunov function V = + ||| and constant
Peek ((x) S (x).).
(III)
Now we have verified the first condition (DRI1) by checking conditions (I,(I), and (III), (DR12) In what follows, check the boundedness and Lipshitz conditions the random-fied field function H (00 x), where each subcomponent is defiend as H0, x) (x)( (Li-)g(0)) - 0())
Recall that V = 1 + |x| 2, the compactness assumption Al directly leads to
sup [H(0, x)]l S mV(x).
(IV)
For any 01, 02 E K. and fixed E X. suffices for us to solely verify the i-th index, which is the index that maximizes [i) (2()) then
(01, H,,,, x)1 = (1 (LimJ((0) (1(1) - (2( (10 (4m0)))))))
(0(0) S max (0)(i) - -0(()) + -10)))0(0000(())))) - (006))))+++00)00000000000) 0260)02(2))
S (Jz((( -
-016))-220111
< 3 (i) 02(i)l,
where the last inequality holds since (()) E (0, 1] for any i S m.
(DRI3) We proceed verify the smoothness of the transitional kernel Po(x, y) with respect to For any 0,, 02 E K. and fixed x and we have
(x, y) (x.
2/ (/TTe)///2 EF 1y-xte 2(00011
1
1
Eeee 1
4e
[x]
[x]
2 (4Arc)/d/2
y x + 01)1 Ily x e 02)
V Txx, 01) - Vx VxTxxx 0x)l
00 021l.
where the first inequality (up a finite constant) follows by ||e* --==ex1p S x|x yll for any x, y in compact space; the last inequality follows by the definition of the adaptive gradient in (15)) and ||x y/l by the compactness assumption Al.
[|oo(x)) loggy))
Following the same technique proposed in Liang et al. (2007) (page 319), we can verify the last drift condition
For f: Rd, define the norm [lf = supxEx V(x)
[f(x)
|Po f - fll S &CIlfllv
O2l1. e = {f X x Rdd IIfllv < <oo..
(VI)
Having conditions (I), (II), ... and (VI) verified, we are now able to prove the drift conditions
proposed in Section 6 of Andrieu et al. (2005).
Before we present the L convergence of og, we make some extra assumptions on the step size.
Assumption AS (Learning rate and step size) The learning rate {EX]NEN a positive non- increasing sequence of real numbers satisfying the conditions
lim C = 0, h
C oc.
k=
22

---

Published as a conference paper at ICLR 2022
The step size keN is positive non-increasing sequence of real numbers such that
lim Wk = O, k-00
Wk = Σ R=1
+00,
Σ we +oc.
(32)
A practical strategy is to set wk == O to satisfy the above conditions for any a E (0.5, 1].
The following is application of Theorem 24 (page 246) (Benveniste al., 1990) given stability conditions (Lemma 3).
Lemma 7 (L² convergence rate, restatement of Lemma 2) Assume Assumptions A1-A5 hold.
For any 0g θ C O, a large m, small learning rates {Ex JE and step sizes [wk] ) converges to 0,, where 0 = 0. + O (supx Var(5n + suP42ko Ek + m) for some ko, such that
E |10. 0.112] O (wx)
The theoretical novelty is that we treat the biased θ. as the equilibrium of the continuous system instead of analyzing how far we away from 0, in all aspects as in Theorem (Deng al., 2020b). This enables us to directly apply Theorem 24 (page 246). Nevertheless, it can be interpreted as a special case of Theorem i (Deng al., 2020b) except that there are perturbation terms and the equilibrium is 0, instead of 0..
C GAUSSIAN APPROXIMATION
C.I PRELIMINARY: SUFFICIENT CONDITIONS FOR WEAK CONVERGENCE
To formally prove the asymptotic normality the stochastic approximation process W/ -1/2 O.). we first lay out preliminary result (Theorem of Pelletier (1998)) that provides sufficient conditions to guarantee the weak convergence.
Lemma 8 (Sufficient Conditions) Consider stochastic algorithm follows
0x+1 = O2 + Wk+1 l((O)) + 15411k= + Wk+1
where D denotes a perturbation and ek+1 is random noise. Given three conditions (CI), (C2). and (C3) defined below, we have the desired weak convergence result
w - (Os θ.) + - (0,
(33)
where E = So the, o R e th dt, R denotes the limiting covariance of the martingale
lim Blex+1eT411 [Fl and F is the o-algebra of the events up to iteration k. he. = he (
wps -wh1
lim
(C1) There exists an equilibrium point θ. and a stable matrix he. := ho)) IRmx such that for any 0 E 10: 0.11 S M) for some M > O, the mean-field function h 1.mh - Rm satisfies
h(0.) 0
[|^(0) - ha, (0 - - S 0 -0.1?,
(C2) The step size wi decays with an order E (0, 1] such that wi = ak- (C3) Assumptions the disturbances. There exists constants M and ā > 2 such that
E [Ft 1110-0. 0. M.M) = 0,
(4)
sup |
1 Fr
1 (1 00,
(12)
E Wi || ||2|
[110 Ø. (1555 0,
(I
+ R. For example, ? = 0 we O(k ") where E (0.5, 1] and & fg if w/ kp.
E |ektieT
[i
1 (110 0, (|
(III)
23

---

Published as a conference paper at ICLR 2022
Remark 2 By the definition of the mean-field function h(0) in Eq,(26), it is easy to verify the condition Cl. Moreover, Assumption A5 also fulfills the condition C2. Then, the proof hinges on the verification of the condition C3.
C.2 PRELIMINARY: CONVERGENCE OF THE COVARIANCE ESTIMATORS
In particular, to verify the condition E [ex+1e441 FF
[110 M. R. we study the con-
vergence the empirical sample mean (xx)] for a test function f to the posterior expectation
f = fx f(x) (x)(.)). Poisson's 's equation is often used to characterize the fluctuation between f(x) and
Lg(x) = f(x)
(34)
where L refers to an infinitesimal generator and gxx) denotes the solution of the Poisson's 's equation. Similar to the proof of Lemma 6, the existence the solution of the Poisson's 's equation has been established in (Mattingly al., 2002; Vollmer al., 2016). Moreover, the perturbations of [E [x(((x))) are properly bounded given regularity properties for gx)), where the O-th, 1st, and 2nd order of the regularity properties has been established in Erdogdu et al. (2018).
The following result helps us to identify the convergence of the covariance estimators, which is adapted from Theorem 5 (Chen al., 2015) with decreasing learning rates {ex The gradient biases from Theorem 2 (Chen et al., 2015) are also included to handle the adaptive biases.
Lemma (Convergence of Covariance Estimators) Suppose Assumptions AI-A5 hold. For any Oy O O, a large m, small learning rates {< step sizes {wi ) and any bounded function f, we have
B B[((xH]] J f())=g (x)dx + O,
where wa. (x) is the invariant measure simulated via SGLD that approximates Two (x) x =(x)
(i(())))
Proof
We study the single-chain CSGLD and reformulate the adaptive algorithm follows:
Xk+1 = X - VxT(xx, 0k) + (O, 2e
= -E (VL((xx, 0.) T(xx, 0n)) + NO 2ExT(),
where VzĪ(x, 1 + &r (log (J - log [( 1) v 1)]] vxū(xx V2I(x, is
defined in Section B.1 and the bias term is given by T(x:, 0%) V2I(xe, 0x) VxI(xx, 0.). Then, by Jensen's inequality and Lemma 7, we have
(xx, 0x)1 < |[/ [xx,
(xx, 0%) - VxI(Xx, 0.)
(35)
$ E 0. @.11 S V [[ - 0_112) < O (V
Combining (3.(35) and Theorem 5 (Chen et al., 2015), we have
E [xxx]] - / fx))FF (x)dx 0
Σt wi E|T(X, 0.)11
Σ
E/C,
E.ws
Σ G
+0 as + 00,
where the last argument directly follows from the conditions on learning rates and step sizes in Assumption A5.
+x(x) Em ilag- <0())s where the exact energy function (xx) is selected.
24

---

Published as a conference paper at ICLR 2022
C.3 PROOF OF THEOREM 1
Recall that the stochastic approximation based on a single process follows from
0x+1
0x + wk+11 H0,,, xx+1)
0x + wk+1 (0b++1((p))) + #kk+1 (Jo,, (xx+1) TTo (x++1))
On + wh+ (1x)))
[00 (xx+1) (x vo, (X++1) + Wh+2 Wh+1
No, 114 (xx+1)
wk+1
wk+1
+ wk+1
Wk+illo. Ho. (xx) Wk+2 1ett2lloss (xx++1) pox (Xx+1)
II (xx)
Wk+11
cx+1
0 + wk+1h((0)) + wk+1 (vk+1 + Sk+11
twk+1
Ck+1
perturbation
martingale
(36)
where the second equality holds from the solution of Poisson's equation in Eq.(30).
We denote O₂ = 02 Wk+1 [ic4iTT,,44 No. (xx). Adding w (X++1) on both sides of (36),
we have
02+1
0
+ wh+ 11(x))(())))))))))))))))))))))) + wk+1 (Mk+1 + eb+1 + 5k+1) + Wk-2lloo Mor- (xx+1) - wk+1 Io. (xk)
0% + + 1(xx) + wk+1 (Wk+1 + Ex+1)
0 + + (0x) +4++1 (Dx+1 + (x+1),
(37)
where x*+1 = Mp1 h(0)) h(@x). Next, we proceed verify the conditions in C3.
(I) By the martingale difference property of {ek] and the compactness assumption A1, we know that for any
E[ex+1| O, E 11 F < co.
k20
(II) By the definition of (Ox)) in (26),,, we can easily check that h(0x) is Lipschitz continuous
a neighborhood 0_. Combining E31, we have [h(O+) - h(@x) (I = O - 0xll) =
(I1w++1 To, (xx) Il) = (w)+1). Then [[[[ & C/ 0z + ((wx+2) = ((wx+1) by
the step size condition (32). In what follows, we can verify
1 + 2
MM++1
[l(Ox)) h h()))
< 2E
2E
O(we) - 0.
(II)
Wk
Wk
Wi
(III) For the martingale difference noise ex+1 = nov (x6+1) 1l. (xx) with mean 0, we have
Elex+1eT411 [FR] = Elue, (xx+1) (oe+1))e, (Xx+1) [Fx] 1l Mo, (xx lee Mex (xx)
We denote Eiextre Fi] by a function f(xx)). Applying Lemma 9, we have
E Elex-1eT Fx] [x[[(x))] + / [(x)= 1 dx= = lim Tee R.
(III)
10000
where R :: R(O.) and R(0) is also equivalent to - Cova xx), (()))))).
Having the conditions C1, C2 and C3 verified, apply Lemma and have the following weak convergence for 0₃
Wx -1/2 - - = ~(O.E) E),
25

---

Published as a conference paper at ICLR 2022
where Σ = Jo ethe. o R o ch dt and he. = ho(0.) + EI, E = limg- 000
-
4/
Considering the definition that O2 On + wk++ 1ok+ (xx) and B Il, uon (xx))l) is uniformly bounded by (3)(3,) we have
w 1/2 No No, (xx) + 0 in probability.
By Slutsky's s theorem, we eventually have the desired result
112 0 O.) = +0 Σ).
We
where the step size wk decays with order (0.5, such that wi O(k
26

---

Published as a conference paper at ICLR 2022
D MORE ON EXPERIMENTS
D.1 MODE EXPLORATION ON MNIST VIA THE SCALABLE RANDOM-FIEL FUNCTION
For the network structure, we follow Jarrett et al. (2009) and choose a standard convolutional neural network (CNN). Such CNN has two convolutional (conv) layers and fully-connected (FC) layers. The two conv layers has 32 and 64 feature maps, respectively. The FC layers both have 50 hidden nodes and the network has 5 outputs. large batch size 2500 selected reduce the gradient noise and reduce the stochastic approximation bias. We fix Ç = 3e4 and weight decay 25. For simplicity, we choose 100,000 partitions and Au 10. The step size follows wi = min (0.01, 2o.4
D.2 SIMULATIONS OF MULTI-MODAL DISTRIBUTIONS
The target density function is given by (x)
exp(-U(x)), where x 2(cos((TTZ1) cos(22Tz2)).
(x), 02) and (x) follows Uxx) = 0.2(21 + x3)
We also include regularization term L(x) (af++3)>200 X
((19 + +2)) 20). This design leads to highly multi-modal distribution with 25 isolated modes. Figure shows the con- tour and the 3-D plot of the target density. The ICSGLD and baseline algorithms are applied to this example. For IC- SGLD, we set E = 3e 3. -= 1, = 0.75 and total number of iterations= S Besides, partition the sample space into 100 subregions with bandwidth Au = 0.125 and set w% min g100+100
-0
40
Figure 5: Target density.
For comparison, run the baseline algorithms under similar settings. For CSGLD, we run a single process 5 times of the the same as those used by ICSGLD. For reSGLD,
time budget and all the settings parallel chains with learning rates 0.001, 0.002, ... 0.005 and temperatures 1, 2, ... .5, respectively. We estimate the correction every 100 iterations. We fix the initial correction 30 and choose the same step size for the stochastic approximation as in ICSGLD. For SGLD, we run five chains in parallel with the learning rate 3e- and a temperature of 1. For cycSGLD, we run a single-chain with times of the time budget. We set the initial learning rate as 1e-2 and choose 10 cycles. For the particle-base based SVGD, we run five chains in parallel. For each chain, we initialize 100 particles as being drawn from a uniform distribution over rectangle, The learning rate is set to 3e
five
CSGLD XT5
CSGLD XT5
2.5
cycSGLD XT5
2.5
cycSGLD XT5
SGLD XP5
XP5
reSGLD xP5
reSGLD xP5
KL divergence
2.0
ICSGLD XP5
divergence
2.0
ICSGLD XP5
1.5
1.5
1.0-
22 1.0-
0.5
0.5-
5000
10000
15000
20000
60
100
Steps
Time (s)
Figure 6: Estimation KL. divergence versus time steps for ICSGLD and baseline methods. We repeat experiments 20 times.
To compare the convergence rates in terms of rumming steps and time between ICSGLD and other algorithms, repeat each algorithm 20 times and calculate the mean and standard error over 20 trials. Note that we run all the algorithms based on parallel chains (xP5) except that cycSGLD and CSGLD in single-chain with 5 times of time budget (xT5) and the steps and running
27

---

Published as a conference paper at ICLR 2022
time are also scaled accordingly. Figure 6 shows that the vanilla SGLDxP5 converges the slowest among the five algorithms due to the lack of mechanism to escape local traps: cycSGLDx slightly alleviates that problem by adopting cyclical learning rates; reSGLDxP5 greatly accelerates the computations by utilizing high-temperature chains for exploration and low-temperature chains for exploitation, but the large correction term inevitably slows down the convergence; ICSGLDxxPS converges faster than all the others and the noisy energy estimators only induce bias for the latent variables and don't It affect the convergence rate significantly.
For the particle-based SVGD method, since more particles require expensive computations while fewer particles lead crude approximation. Therefore, we don't show the convergence of SVGD and only compare the Monte Carlo methods.
D.3 DEEP CONTEXTUAL BANDITS ON MUSHROOM TASKS
For the UCI Mushroom data set, each mushroom is either edible or poisonous. Eating an edible mushroom yields reward of 5, but eating poisonous mushroom has a 50% chance to result in reward of -35 and reward of 5 otherwise. Eating nothing results in 0 reward. All the agents use the same architecture. In particular, we fit two-layer neural network with 100 neurons each and ReLU activation functions. The input of the network there are 2 outputs, representing the predicted reward for eating or not eating a mushroom. The mean squared loss is adopted for training the models. We initialize 1024 data points and keep a data buffer of size 4096 as the training proceeds. The size of the mini-batch data is set to 512. To adapt online scenarios, we train models after every 20 new observations.
feature vector with dimension 22 (context) and
We choose one e-greedy policy (EpsGreedy) based on the RMSProp optimizer with decaying learn- ing rate (Riquelme et al., 2018) as a baseline. Two variational methods, namely stochastic gradient descent with constant learning rate (ConstSGD) (Mandt al., 2017) and Monte Carlo Dropout (Dropout) (Gal & Ghahramani, 2016) are compared to approximate the posterior distribution. For the sampling algorithms, we include preconditioned SGLD (pSGLD) (Li et al., 2016), preconditioned CSGLD (pCSGLD) (Deng et al., 2020b), and preconditioned ICSGLD (pICSGLD). Note that all the algorithms run parallel chains with average outputs ( (xP4) except that PCSGLD runs single-chain with 4 times of computational budget (x T4). In particular for the two contour algorithms, we set <= 20 and choose constant step size for the stochastic approximation fit for the time-varying posterior distributions. For more details on the experimental setups, we refer readers to section D in the supplementary material.
We report the experimental setups for each algorithm. Similar to Table 2 of Riquelme al. (2018), the inclusion of advanced techniques may change the optimal settings the hyperparameters. Nevertheless, try report the best setups for each individual algorithm. We train each algorithm 2000 steps. We initialize 1024 mushrooms and keep a data buffer of size 4096 as the training proceeds. For each step, are given 20 random mushrooms and train the model 16 iterations every step for the parallel algorithms (xP4); we train PCSGLDX T4 64 iterations every step.
EpsGreedy decays the learning rate by factor of 0.999 every step; by contrast, all the others choose a fixed learning rate. RMSprop adopts a regularizer of 0.001 and a learning rate of 0.01 to learn the preconditioners. Dropout proposes a 50% dropout rate and each subprocess simulates 5 models for predictions. For the two importance sampling (IS) algorithms, we partition the energy space into m 100 subregions and the energy depth Au as 10. We fix the hyperrameter ( 20. The step sizes for pICSGLDxPP XP4 and PCSGLD X4 are chosen as 0.03 and 0,006, respectively, A proper regularizer is adopted for the low importance weights. See Table 2 for details.
TABLE 2: DETAILS OF THE EXPERIMENTAL SETUPS.
ALGORITHM
LEARNING RATE
Temperature
RMSprop
Train
Dropout
e-Greedy
EPSGREEDY x P4
5c-7 (0.999)
YES
No
16
No
0.3%
CONSTSGD
le-6
No
No
16
No
No
DROPOUTXP4
le-6
No
No
16
YES (50%)
No
PCSGLDxT4
5e-8
YES
64
No
No
PS
3c-7
0.3
No
16
No
No
PICSGLDxP4
3c-7
0.3
YES
YES
16
No
No
28

---

Published as a conference paper at ICLR 2022
D.4 UNCERTAINTY ESTIMATION
All the algorithms, excluding M-SGD x P4, choose temperature of 0.0003 1, We run the parallel algorithms 500 epochs Px) and the single-chain algorithms 2000 epochs (xT4). The initial learning rate is 2e-6 (Bayesian settings), which corresponds to the standard 0. for averaged data likelihood.
We train cycSGHMCx T4 and MuldSWAGx based on the cosine learning rates with 10 cycles. The learning rate in the last 15% of each cycle is fixed a constant value. MultiSWAG simulates 10 random models at the end of each cycle. M-SGDxP4 follows the same cosine learning rate strategy with one cycle.
reSGHMC P4 proposes swaps between neighboring chains and requires a fixed correction of 4000 for ResNet20, 32, and and a correction of 1000 WRN-16-8. The learning rate is annealed at 250 and 375 epochs with a factor of 0.2. ICSGHMCXP4 also applies the same learning rate. We choose m = 200 and Au = 200 for ResNet20, 32, and 56 and Au 60 for WRN-16-8, Proper regularizations may be applied to the importance weights and gradient multipliers for training deep neural networks.
Variance reduction (Deng et al., 2021a) only applies to reSGHMC P4 and ICSGHMCXP4 because they are the only two algorithms that require accurate estimations the energy. We only update control variates every 2 epochs in the last 100 epochs, which maintain reasonable training time and higher reduction of variance due to a small learning rate. Other algorithms yield a worse performance when variance reduction is applied to the gradients.
D.5 EMPIRICAL VALIDATION OF REDUCED VARIANCE
To compare the 0's learned from ICSGLD and CSGLD, we try to simulate from Gaussian mixture distribution (-4.V(-6, 6, 1)
0.6 where (u, 2) denotes a Gaussian distribution with mean and standard deviation We fix 0.9 and Au = 1. We run ICSGLD with 1,000,000 iterations based on 10 interacting parallel chains and run CSGLD with 10,000,000 iterations using single chain. We refer to them as ICSGLDx P1O and CSGLD x T10, respectively. The rest of the settings follows from the experimental setup in section
PDF in Energy (dertty of statess
4.1 (Deng et al., 2020a).
To measure the variance of the estimates, repeated the experiments 10 times and present the mean and two stan- dard deviations for both CSGLDxTI0 and ICSGLDxPII in Figure 7. The results indicate that both estimates of (f (by CSGLD and ICSGLD) converge to the equilibrium that approximates the ground truth of the density of states. No- tably, ICSGLDxX P1O yields a significantly smaller variance than CSGLDxT10 but with the same computational budget. This shows the clear advantage of ICSGLD (many interact- ing short runs) over CSGLD (a single long run) in tackling the large variance issue for importance sampling.
ICSGLD xP10
CSGLD xT10
1aa
Energy
Figure 7: ICSGLD v.s. CSGLD.
WWe various data augmentation techniques, such random flipping, cropping, and random erasing (Zhong al., 2017). This leads much concentrated posterior and requires very low temperature.
29


Reviews:
[{"review_id": "ppk5eN-KALF", "reviewer": "Reviewer_aTJ6", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The idea of Monte Carlo average of random-field functions from multiple chains in this paper is straightforward and easy to implement. The theoretical analysis seem to be solid too. However, I anticipate more discussion on the connection between empirical results and theoretical implications:\n+ The analysis proves that ICSGLD improves estimation accuracy of parameter $\\theta$ with non-infinitesimal step size. However, does parameter $\\theta$ have any geometric interpretation (e.g. density of states)? Can we visualize the $\\theta$ learned from ICSGLD and CSGLD? Can we empirically verify that the estimated variance is reduced?\n+ Could the authors compare ICSGLD with different number of chains to further justify the theoretically result?\n+ What are the main factors affecting the target $\\theta_\\ast$? Is it sensitive to model structure / model size / dataset ?\n+ Is there any theoretical guide on the choice of hyper parameters e.g. $\\zeta$? The first experiment uses $\\zeta = 3e4$ and last experiment uses $\\zeta = 3\\times 10^6$. These values are significantly larger than other two experiments where $\\zeta = 0.75$ and $\\zeta = 20$.\nHow these parameters are selected in different experiments?\n+ What is the recommended $\\Delta u$ for neural network? What are the pros and cons for choosing small $\\Delta u$?\n+ What is the definition of preconditioned CSGLD and preconditioned ICSGLD? In section 5.3, authors cite [Deng et al., 2020b] for preconditioned CSGLD. However, it seems this paper doesn't contain preconditioned version of CSGLD.\n\nOthers:\n+ How KL divergence is estimated in Figure 6?\n+ It will be appreciated if the code is provided for reproduction.\n", "summary_of_the_paper": "This paper samples from distributions with complex energy landscapes by the powerful contour stochastic gradient Langevin dynamics (CSGLD) methods. In order to reduce the variance of random-field function, this paper proposes a simple but effective method, that combines  random-field functions from multiple chains. Since the parameter $\\theta$ typically has small size, the communication overhead is marginal. This paper show that compared to CSGLD which only has one chain, the variance ICSGLD with P chains and same computation budget is asymptotically more efficient. The authors evaluate ICSGLD with multiple experiments.", "summary_of_the_review": "This paper proposes an interesting extension to CSGLD whose advantages are shown both theoretically and empirically. However, I believe many questions deserve further investigation before applying ICSGLD into practical machine learning. I would like to hear authors' response on how theoretical implications are connected to practical concerns.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "E_RapCEtLhF", "reviewer": "Reviewer_9d9U", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The strength of the paper is that it proposes a new algorithm interacting contour stochastic gradient Langevin dynamics that can be shown theoretically more efficient than contour stochastic gradient Langevin dynamics. The paper is well written, and it has both theoretical justification and numerical support. \n\nThe weakness of the paper is that many results seem to be adapted from Deng et al. (2020) for the contour stochastic gradient Langevin dynamics. Even though, a new algorithm is proposed, it should be considered as an extension to the previous work Deng et al. (2020), and some similar theoretical estimates can be found in Deng et al. (2020) and the technical novelty is not clear. ", "summary_of_the_paper": "The paper proposed an interacting contour stochastic gradient Langevin dynamics that extends the contour stochastic gradient Langevin dynamics with efficient interactions in Deng et al. (2020). They obtain asymptotic normality as well as show that this proposed algorithm can be more efficient than  the contour stochastic gradient Langevin dynamics with efficient interactions. The theoretical results are complemented by the numerical experiments.", "summary_of_the_review": "(1) Page 2. For ICSGLD and CSGLD, I have one naive question. You are sampling (2) with a flattened density compared to the original Gibbs distribution $\\pi$. If you are sampling a modified distribution instead, do you need to show that the new target distribution is close to the original distribution?\n\n(2) Page 2. I am a bit confused with the notation here. In equation (2), you have $\\Psi^{\\zeta}_{\\theta}(U(x))$ and then in the bottom of page 2, you have $\\tilde{U}_{\\Psi_{\\theta_{k}}}$. What's the relation between these two different notations?\n\n(3) Page 3. In your Algorithm 1, there is a stochastic approximation step for $\\theta_{k}$. In your Assumption A2 in the appendix, you assume that the space $\\Theta$ is compact. How do you guarantee your $\\theta_{k}$ lives in a compact space since you are doing $k\\rightarrow\\infty$ analysis.\n\n(4) Your Theorem 1 and Corollary 1 for the asymptotic normality is nice. But do you have a result saying that $x_{k}$ will converge to the Gibbs distribution corresponding to $\\theta_{\\star}$ or something like that, e.g. ergodic mean etc.?\n\n(5) Lemma 2 (Lemma 7) is very similar to the result for CSGLD (Theorem 1, Deng et al. (2020).) What's the technical novelty here?", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "6aK2w4NZ8y", "reviewer": "Reviewer_mVcL", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I liked reading the paper. Overall, the work addresses a relevant problem and provides a nice introduction to previously existing solutions. \n\nIn my opinion, the contributions are significant and the proposed algorithm is presented in a clear way. The analysis is sound and backed up by numerical evidence in different settings. I believe that the paper is worth publication. \n\nHowever, I have some minor questions that require a clarification from the authors.\n\n(i) I do not understand how the parameter $\\zeta$ is tuned in practice and I would like to have more details on this procedure in addition to the comments provided in appendix B.1.3. This point seems to play an important role, e.g., in escaping from local regions of the landscape in Figure 2, while Appendix D.1 only says that $\\zeta$ is fixed to $3e4$. Since in the other cases (appendix D.2, D.3) $\\eta$ is fixed to much smaller values (namely, $0.75$ and $20$) I would like to have more intuition on this.\n\n(ii) I would like to have some clarifications on the link between the second footnote on page 7 and the sentence it refers to.", "summary_of_the_paper": "This paper addresses the problem of sampling from distributions with complex energy landscapes. The authors propose an extension of the contour Stochastic Gradient Langevin dynamics (CSGLD) sampler to efficiently simulate from big-data distributions. The extension is called “interacting” since $P$ different CSGLD are updated simultaneously and the random-field function is obtained as a Monte-Carlo average over the random field functions of these multiple chains run in parallel. This procedure allows to reduce the variance of the self-adapting parameter $\\theta$ improving the marginal energy likelihood estimate, and is well-suited for distributed computing. The different chains only share the low-dimensional vector $\\theta$ and do not share the model parameters, that are typically high dimensional. The authors prove that the proposed algorithm is asymptotically more efficient than its single-chain counterpart with the same computational budget. This theorem holds for a step size that decreases polynomially in time. Finally, they compare their algorithm to alternative methods on different tasks (e.g., mode exploration on MNIST dataset, uncertainty estimation) achieving competitive performances. ", "summary_of_the_review": "I believe that the paper is worth of publication since it presents an efficient importance sampling method that scales to big data problems and possibly leads to additional future investigations. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "3nH0FKr3clc", "reviewer": "Reviewer_CBb1", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "My main concern is the clarity of the paper. At some points, the paper looks like a riddle where you are looking for some notation and definitions. Some points of the exposition lack intuition. \n- The choice of the piecewise continuous approximation of the target density and the function H is not apparent at all. I would suggest sharing a piece of intuition with the reader to motivate the usage of these functions early on. The clarity is also hardened by the fact that the original paper by [Deng, 2020] also doesn't provide any intuition for the choice.\n- The proposed algorithm is presented in the background section together with the algorithm by [Deng, 2020]. It would be much easier for the reader if the authors would highlight the difference between the original algorithm and its modification. For instance, they could say that the main difference between algorithms is in the estimation of parameters $\\theta$ and using all of the samples instead of a single chain or even highlight it explicitly in the formulas by using a different color.\n- Pictures in Figure 1 are not informative. The pictures depict different aspects of interacting chains algorithms. Fig. 1a represents the jumps in the state space, and Fig. 1b represents the mutual update of parameters. The authors use Fig. 1b as a motivation for better mixing properties, claiming that unlike the replica-exchange method their algorithm shares \"messages\" between samples at every step. However, one can argue that the replica-exchange method also shares information between samples at every step by running the accept/reject test trying to swap the samples.\n\nAnother major issue is that the authors do not compare their algorithm against CSGLD (algorithm by [Deng, 2020]) in one of the main experiments: Bayesian networks on CIFAR-100 and SVHN. Due to the close relations between algorithms, such kind of comparison is essential for the performed empirical study.\n", "summary_of_the_paper": "The paper proposes an MCMC algorithm containing two key ideas: the approximation of the target density with a simpler function (proposed by [Deng, 2020]) and the parallel simulation of many chains that allow for better capturing the global properties of the target density. To be more precise, the main contribution of the paper is the extension of the algorithm by [Deng, 2020] to the case of multiple chains. \n\nThe proposed extension goes as follows. The original algorithm by [Deng, 2020] approximates the target density with piecewise continuous functions using the level sets (slicing the values of the energy function). These functions are parameterized by vector $\\theta$, which needs to be defined. In practice, $\\theta$ is estimated using the samples from the chain. The current paper then proposes to run many chains in parallel and use all of the available samples for the estimation.\n\nThe efficiency of the proposed technique is analyzed both theoretically and empirically. For theoretical analysis, the authors provide convergence speed and asymptotic distribution of the residuals $\\theta - \\theta^*$. Overall, the analysis favors the proposed method over the original algorithm. The empirical comparison is made on several tasks. Except for the toy task, the authors consider learning policy for contextual bandit problem and sampling from the posterior distribution of a neural network.\n\n[Deng, 2020]: Deng, Wei, Guang Lin, and Faming Liang. \"A Contour Stochastic Gradient Langevin Dynamics Algorithm for Simulations of Multi-modal Distributions.\" arXiv preprint arXiv:2010.09800 (2020).", "summary_of_the_review": "This submission is a follow-up to [Deng, 2020] that extends the previous work to the case of multiple parallel chains. Although the extension seems quite straightforward, the authors provide a theoretical analysis of the proposed scheme and perform its empirical study. The main issues are the clarity of the exposition and one inconsistency in the experimental section.\n\n[Deng, 2020]: Deng, Wei, Guang Lin, and Faming Liang. \"A Contour Stochastic Gradient Langevin Dynamics Algorithm for Simulations of Multi-modal Distributions.\" arXiv preprint arXiv:2010.09800 (2020).", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes a new variant of a stochastic gradient Langevin dynamics sampler that relies on two key ideas: approximation of the target density with a simpler function (as in [Deng, 2020]) and the parallel simulation of many chains. The authors also prove that their approach can be theoretically more efficient than a single-chain algorithm.

The reviewers see the contribution as significant although they did raise some concerns regarding the clarity of the paper. Since these concerns do not appear to be major, I recommend acceptance but I advise the authors to address the comments of the reviewers to maximize the impact of the paper.
Decision: Accept (Poster)