INTERACTING CONTOUR STOCHASTIC GRADIENT
Year: 2022

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.

heading: 1 INTRODUCTION
text: Stochastic gradient Langevin dynamics (SGLD) (Welling & Teh, 2011) has achieved great successes in simulations of high-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 et 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 et 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 a multitude of modes well separated by high energy barriers, an emerging trend is to 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/population simulations. Other multiple chain methods include particle-based nonlinear Markov (Vlasov) processes (Liu & Wang, 2016; Zhang et al., 2020a) and replica exchange methods (also known as 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 communication cost in that each chain shares with others the marginal energy likelihood estimate only. As a 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 et al., 1993; Doucet et al., 2001;
Wang & Landau, 2001; Liang et al., 2007; Andrieu et al., 2010; Deng et al., 2020b) in big data applications. To resolve this issue, we resort to a novel interacting random-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 a scalable interacting importance sampling method for big data problems with the minimal 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 a non-linear mean-field system and justify regularity properties of the solution of Poisson’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 a fixed learning rate on the MNIST dataset and obtains remarkable performance in large-scale uncertainty estimation tasks.

heading: 2 PRELIMINARIES

heading: 2.1 STOCHASTIC GRADIENT LANGEVIN DYNAMICS
text: A standard sampling algorithm for big data problems is SGLD (Welling & Teh, 2011), which is a numerical scheme of a stochastic differential equation in mini-batch settings:
xk+1 = xk − k N
n ∇xŨ(xk) +
√ 2τ kwk, (1)
where xk ∈ X ∈ Rd, k is the learning rate at iteration k, N denotes the number of total data points, τ is the temperature, andwk is a standard Gaussian vector of dimension d. In particular, Nn∇xŨ(x) is an unbiased stochastic gradient estimator based on a mini-batch data B of size n and Nn Ũ(x) is the unbiased energy estimator for the exact energy function U(x). Under mild conditions on U , xk+1 is known to converge weakly to a unique invariant distribution π(x) ∝ e− U(x) τ as k → 0.

heading: 2.2 CONTOUR STOCHASTIC GRADIENT LANGEVIN DYNAMICS
text: 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 a flattened density with much lower energy barriers
$Ψθ (x) ∝ π(x)/Ψ ζ θ(U(x)), (2)
where ζ is a hyperparameter, Ψθ(u) = ∑m i=1 ( θ(i− 1)e(log θ(i)−log θ(i−1)) u−ui−1 ∆u ) 1ui−1<u≤ui .
In particular, {ui}mi=0 determines the partition {Xi}mi=1 of X such that Xi = {x : ui−1 < U(x) ≤ ui}, where −∞ = u0 < u1 < · · · < um−1 < um = ∞. For practical purposes, we assume ui+1 − ui = ∆u for i = 1, · · · ,m− 2. In addition, θ = (θ(1), θ(2), . . . , θ(m)) is the self-adapting
parameter in the space Θ = { (θ(1), · · · , θ(m)) ∣∣0 < θ(1), · · · , θ(m) < 1& ∑mi=1 θ(i) = 1}.
Ideally, setting ζ = 1 and θ(i) = θ∞(i), where θ∞(i) = ∫ Xi π(x)dx for i ∈ {1, 2, · · · ,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 θ∞ is unknown a priori. To tackle this issue, CSGLD proposes the following procedure to adaptively estimate θ via stochastic approximation (SA) (Robbins & Monro, 1951; Benveniste et al., 1990):
(1) Sample xk+1 = xk + k Nn∇xŨΨθk (xk) + √ 2τ kwk,
(2) Optimize θk+1 = θk + ωk+1H̃(θk,xk+1),
Algorithm 1 Interacting contour stochastic gradient Langevin dynamics algorithm (ICSGLD). {Xi}mi=1 is pre-defined partition and ζ is a hyperparameter. The update rule in 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 Bk from D, and compute stochastic gradients ∇xŨ(x(p)k ) and energies Ũ(x (p) k ) for each x
(p), where p ∈ {1, 2, · · · , P}, |Bk| = n, and |D| = N . [2.] (Parallel simulation) Sample x ⊗ P
k+1 := (x (1) k+1,x (2) k+1, · · · ,x (P ) k+1) > based on SGLD and θk
x ⊗ P k+1 = x ⊗ P k + k N
n ∇xŨΨθk (x
⊗ P
k ) + √ 2τ kw ⊗ P k , (4)
where k is the learning rate, τ is the temperature, w ⊗ P
k denotes P independent standard Gaussian vectors, ∇xŨΨθ (x ⊗ P ) = (∇xŨΨθ (x(1)),∇xŨΨθ (x(2)), · · · ,∇xŨΨθ (x(P )))>, and
∇xŨΨθ (x) = [ 1 + ζτ∆u ( log θ(JŨ (x))− log θ((JŨ (x)− 1) ∨ 1) )] ∇xŨ(x) for any x ∈ X . [3.] (Stochastic approximation) Update the self-adapting parameter θ(i) for i ∈ {1, 2, · · · ,m}
θk+1(i) = θk(i) + ωk+1 1
P P∑ p=1 θk(JŨ (x (p) k+1)) ( 1 i=JŨ (x (p) k+1) − θk(i) ) , (5)
where 1A is an indicator function that takes value 1 if the event A appears and equals 0 otherwise.
where∇xŨΨθ (·) is a stochastic gradient function of$Ψθ (·) to be detailed in Algorithm 1. H̃(θ,x) :=( H̃1(θ,x), · · · , H̃m(θ,x) ) is random-field function where each entry follows
H̃i(θ,x) = θζ(JŨ (x)) ( 1i=JŨ (x) − θ(i) ) ,where JŨ (x) = m∑ i=1 i1ui−1<Nn Ũ(x)≤ui . (3)
Theoretically, CSGLD converges to a sampling-optimization equilibrium in the sense that θk approaches to a fixed point θ∞ and the samples are drawn from the flattened density $Ψθ∞ (x). Notably, the mean-field system is globally stable with a unique stable equilibrium point in a small neighborhood of θ∞. Moreover, such an appealing property holds even when U(x) is non-convex.

heading: 3 INTERACTING CONTOUR STOCHASTIC GRADIENT LANGEVIN DYNAMICS
text: The major goal of interacting CSGLD (ICSGLD) is to improve the efficiency of CSGLD. In particular, the self-adapting parameter θ 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 θk’s is the key to the success of such a dynamic importance sampling algorithm. To this end, we propose an efficient variance reduction scheme via interacting parallel systems to improve the accuracy of θk’s.

heading: 3.1 INTERACTIONS IN PARALLELISM
text: Now we first consider a naı̈ve parallel sampling scheme with P chains as follows
x ⊗ P k+1 = x ⊗ P k + k N
n ∇xŨΨθk (x
⊗ P
k ) + √ 2τ kw ⊗ P k ,
where x ⊗ P = (x(1),x(2), · · · ,x(P ))>, w ⊗ P
k denotes P independent standard Gaussian vectors, and ŨΨθ (x ⊗ P ) = (ŨΨθ (x (1)), ŨΨθ (x (2)), · · · , ŨΨθ (x(P )))>.
Stochastic approximation aims to find the solution θ of the mean-field system h(θ) such that
h(θ) = ∫ X H̃(θ,x)$θ(dx) = 0,
where $θ is the invariant measure simulated via SGLD that approximates $Ψθ in (2) and H̃(θ,x) is the novel random-field function to be defined later in (8). Since h(θ) is observable only up to large random perturbations (in the form of H̃(θ,x)), the optimization of θ based on isolated
random-field functions may not be efficient enough. However, due to the conditional independence of x(1),x(2), · · · ,x(P ) in parallel sampling, it is very natural to consider a Monte Carlo average
h(θ) = 1
P P∑ p=1 ∫ X H̃(θ,x(p))$θ(dx (p)) = 0. (6)
Namely, we are considering the following stochastic approximation scheme
θk+1 = θk + ωk+1H̃(θk,x ⊗ P k+1 ), (7)
where H̃(θk,x ⊗ P k+1 ) is an interacting random-field function H̃(θk,x ⊗ P k+1 ) = 1 P ∑P p=1 H̃(θk,x (p) k+1). Note that the Monte Carlo average is very effective to reduce the variance of the interacting randomfield function H̃(θ,x ⊗ 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 θ and a large-variance random-field function H̃(θ,x) individually, leading to coarse estimates in the end.
Formally, for the population/ensemble interaction scheme (7), we define a novel random-field function H̃(θ,x) = (H̃1(θ,x), H̃2(θ,x), · · · , H̃m(θ,x)), where each component satisfies
H̃i(θ,x) = θ(JŨ (x)) ( 1i=JŨ (x) − θ(i) ) . (8)
As shown in Lemma 1, the corresponding mean-field function proposes to converge to a different fixed point θ?, s.t.
θ?(i) ∝ (∫ Xi e− U(x) τ dx ) 1 ζ ∝ θ 1 ζ ∞(i). (9)
A large data set often renders the task of estimating θ∞ numerically challenging. By contrast, we resort to a different solution by estimating θ? instead based on a 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.

heading: 3.2 RELATED WORKS
text: Replica exchange SGLD (Deng et al., 2020a; 2021a) has successfully extended the traditional replica exchange (Swendsen & Wang, 1986; Geyer, 1991; Earl & Deem, 2005) to big data problems. However, it works with two chains only and has a low swapping rate. As shown in Figure 1(a), a naı̈ve 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.
Stein variational gradient descent (SVGD) (Liu & Wang, 2016) is a popular approximate inference method to drive a 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-visited 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 more particles are proposed.
Admittedly, ICSGLD is not the first interacting importance sampling algorithm. For example, a population stochastic approximation Monte Carlo (pop-SAMC) algorithm has been proposed in Song et al. (2014), and an interacting particle Markov chain Monte Carlo (IPMCMC) algorithm has been proposed in 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 samplers, which are hard to be adapted to high-dimensional big data problems.
Other parallel SGLD methods (Ahn et 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 x ∈ Rd but only share θ ∈ Rm 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 Ũ(x) is replaced with the random field function H̃(θ,x) and we only share the low-dimensional latent vector θ. Since privacy concerns and communication cost are not major bottlenecks of our problem, we leave the study of taking the Monte Carlo average in Eq.(6) every K > 1 iterations for future works.

heading: 4 CONVERGENCE PROPERTIES
text: 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 approximation 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.

heading: 4.1 LOCAL STABILITY FOR NON-LINEAR MEAN-FIELD SYSTEMS IN BIG DATA
text: The first obstacle for the theoretical study is to approximate the components of θ∞ corresponding to the high energy region. To get around this issue, the random field function H̃(θ,x) in (8) is adopted to estimate a different target θ? ∝ θ 1 ζ ∞. As detailed in Lemma 3 in the supplementary material, the mean-field equation is now formulated as follows
hi(θ) ∝ θζ?(i)− (θ(i)Cθ) ζ + perturbations, (10) where Cθ = ( Z̃ζ,θ
Z̃ζζ,θ?
) 1 ζ
and Z̃ζ,θ = ∑m k=1 ∫ Xk π(x)dx θζ−1(k) . We see that (10) may not be linearly stable
as in Deng et 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 θ, which can be achieved by pre-training the model long enough time through SGLD, the mean value theorem implies a linear property in a local region
hi(θ) ∝ θ?(i)− θ(i) + perturbations. Combining the perturbation theory (Vanden-Eijnden, 2001), we present the following stability result:
Lemma 1 (Local stability, informal version of Lemma 3) Assume Assumptions A1-A4 (given in the supplementary material) hold. For any properly initialized θ, we have 〈h(θ),θ − θ̂?〉 ≤ −φ‖θ − θ̂?‖2, where θ̂? = θ? + O ( supx Var(ξn(x)) + + 1 m ) , θ? ∝ θ 1 ζ ∞ , φ > 0, denotes a learning rate, and ξn(x) denotes the noise in the stochastic energy estimator of batch size n and Var(·) denotes the variance.
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 equation in Lemma 6 in the supplementary material. In what follows, we can control the fluctuations in stochastic approximation and eventually yields the L2 convergence.
Lemma 2 (L2 convergence rate, informal version of Lemma 7) Given standard Assumptions A1A5. θk converges to θ̂?, where θ̂? = θ? +O ( supx Var(ξn(x)) + + 1 m ) , such that
E [ ‖θk − θ̂?‖2 ] = O (ωk) .
The result differs from Theorem 1 of Deng et al. (2020b) in that the biased fixed point θ̂? instead of θ? is treated as the equilibrium of the continuous system, which provides us a user-friendly proof. Similar techniques have been adopted by Durmus & Éric Moulines (2017); Xu et al. (2018). Although the global stability (Deng et al., 2020b) may be sacrificed when ζ 6= 1 based on Eq.(8), θ? ∝ θ 1 ζ ∞ is much easier to estimate numerically for any i that yields 0 < θ∞(i) 1 based on a large ζ > 1.

heading: 4.2 ASYMPTOTIC NORMALITY
text: To study the asymptotic behavior of ω− 1 2
k (θk − θ̂?), where θ̂? is the equilibrium point s.t. θ̂? = θ?+O ( Var(ξn(x)) + + 1 m ) , we consider a fixed step size ω in the SA step for ease of explanation. Let θ̄t denote the solution of the mean-field system in continuous time (θ̄0 = θ0), and rewrite the single-chain SA step (7) as follows
θk+1 − θ̄(k+1)ω = θk − θ̄kω + ω ( H(θk,xk+1)−H(θ̄kω,xk+1) ) + ω ( H(θ̄kω,xk+1)− h(θ̄kω) ) − ( θ̄(k+1)ω − θ̄kω − ωh(θ̄kω) ) .
Further, we set θ̃kω := ω− 1 2 (θk − θ̄kω). Then the stochastic approximation differs from the mean field system in that
θ̃(k+1)ω = ω 1 2 k∑ i=0 ( H(θi,xi+1)−H(θ̄iω,xi+1) ) ︸ ︷︷ ︸
I: perturbations
+ω 1 2 k∑ i=0 ( H(θ̄iω,xi+1)− h(θ̄iω) )︸ ︷︷ ︸ II: martingale Mi −ω 1 2 · remainder
≈ ω 1 2 k∑ i=0
hθ(θiω) (θi − θ̄iω)︸ ︷︷ ︸ ≈ω 1 2 θ̃iω +ω 1 2
k∑ i=0 Mi ≈ ∫ (k+1)ω 0 hθ(θ̄s)θ̃sds+ ∫ (k+1)ω 0 R 1 2 (θ̄s)dWs,
where hθ(θ) := ddθh(θ) is a matrix,W ∈ R m is a standard Brownian motion, the last term follows from a certain central limit theorem (Benveniste et al., 1990) andR denotes the covariance matrix of the random-field function s.t. R(θ) := ∑∞ k=−∞ Covθ(H(θ,xk), H(θ,x0)).
We expect the weak convergence of Uk to the stationary distribution of a diffusion
dUt = hθ(θt)Utdt+R 1/2(θt)dWt, (11)
where Ut = ω −1/2 t (θt − θ̂?). Given that θt converges to θ̂? sufficiently fast and the local linearity of hθ, the diffusion (11) resembles the Ornstein–Uhlenbeck process and yields the following solution
Ut ≈ e−thθ(θ̂?)U0 + ∫ t
0
e−(t−s)hθ(θ̂?) ◦R(θ̂?)dWs.
Then we have the following theorem, whose formal proof is given in section C.3.
Theorem 1 (Asymptotic Normality) Assume Assumptions A1-A5 (given in the supplementary material) hold. We have the following weak convergence
ω −1/2 k (θk − θ̂?)⇒ N (0,Σ), where Σ = ∫ ∞ 0 ethθ? ◦R ◦ eth > θ?dt, hθ? = hθ(θ̂?).

heading: 4.3 INTERACTING PARALLEL CHAINS ARE MORE EFFICIENT
text: For clarity, we first denote an estimate of θ based on ICSGLD with P interacting parallel chains by θPk and denote the estimate based on a single-long-chain CSGLD by θkP .
Note that Theorem 1 holds for any step size ωk = O(k−α), where α ∈ (0.5, 1]. If we simply run a single-chain CSGLD algorithm with P times of iterations, by Theorem 1,
ω −1/2 kP (θkP − θ̂?)⇒ N (0,Σ).
As to ICSGLD, since the covariance Σ relies on R, which depends on the covariance of the martingale {Mi}i≥1, the conditional independence of x(1),x(2), · · · ,x(P ) naturally results in an efficient variance reduction such that
Corollary 1 (Asymptotic Normality for ICSGLD) Assume the same assumptions. For ICSGLD with P interacting chains, we have the following weak convergence
ω −1/2 k (θ P k − θ̂?)⇒ N (0,Σ/P ).
That is, under a similar computational budget, we have ‖Var(θkP−θ̂?)‖F ‖Var(θPk −θ̂∗)‖F = wkPwk/P ≈ P 1−α.
Corollary 2 (Efficiency) Given a decreasing step size ωk = O(k−α), where 0.5 < α < 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 a better non-asymptotic performance (Benveniste et al., 1990).

heading: 5 EXPERIMENTS

heading: 5.1 LANDSCAPE EXPLORATION ON MNIST VIA THE SCALABLE RANDOM-FIELD FUNCTION
text: This section shows how the novel random-field 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 get stuck in few local modes. To simplify the experiments, we choose a large batch size of 2500 and only pick the first five classes, namely digits from 0 to 4. The learning rate is fixed to 1e-6 and the temperature is set to 0.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 re-visiting these zones keeps increasing until a negative learning rate is injected to encourage explorations.
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 a local region and then escapes it once it overvisits this region. This shows the strength of ICSGLD in the simulations of complex multi-modal distributions. More experimental details are presented in section D.1 of the supplementary material.

heading: 5.2 SIMULATIONS OF MULTI-MODAL DISTRIBUTIONS
text: 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 (cycSGLD), replica exchange SGLD (reSGLD), and the particle-based SVGD.
The target multi-modal density is presented in Figure 3(a). Figure 3(b-g) displays the empirical performance of all the testing methods: the vanilla SGLD with 5 parallel chains (×P5) undoubtedly
§The random-field function (Deng et al., 2020b) 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 to a more concentrated posterior (Wenzel et al., 2020; Aitchison, 2021).
performs the worst in this example and fails to quantify the weights of each mode correctly; the single-chain cycSGLD with 5 times of iterations (×T5) improves the performance but is still not accurate enough; reSGLD (×P5) and SVGD (×P5) have good performances, while the latter is quite costly in computations; ICSGLD (×P5) 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 (×T5) 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.

heading: 5.3 DEEP CONTEXTUAL BANDITS ON MUSHROOM TASKS
text: This section evaluates ICSGLD on the contextual bandit problem based on 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 is 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 -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 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 4 parallel chains with average outputs (×P4) except that pCSGLD runs a single-chain with 4 times of computational budget (×T4). For more details, we refer readers to section D.3 in the supplementary material.
Figure 4 shows that EpsGreedy ×P4 tends to explore too much for a long horizon as expected; ConstSGD×P4 and Dropout×P4 perform poorly in the beginning but eventually outperform EpsGreedy ×P4 due to the inclusion of uncertainty for exploration, whereas the uncertainty seems to be inadequate due to the nature of variational inference. By contrast, pSGLD×P4 significantly
outperforms the variational methods by considering preconditioners within an exact sampling framework (SGLD). As a unique algorithm that runs in a single-chain manner, pCSGLD×T4 leads to the worst performance due to the inefficiency in learning the selfadapting parameters, fortunately, pCSGLD×T4 slightly outperform pSGLD×P4 in the later phase with the help of the well-estimated self-adapting parameters. Nevertheless, pICSGLD×P4 propose to optimize the shared self-adapting parameters at the same time, which in turn greatly contributes to the simulation of the posterior. As a result, pICSGLD×P4 consistently shows the lowest regret excluding the very early period. This shows the great superiority of the interaction mechanism in learning the self-adapting parameters for accelerating the simulations.

heading: 5.4 UNCERTAINTY ESTIMATION
text: This section evaluates the qualify of our algorithm in uncertainty quantification. For model architectures, we use residual networks (ResNet) (He et al., 2016) and a 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, a 16-layer WRN that is 8 times wider than ResNet16. We train the models on
CIFAR100, and report the test accuracy (ACC) and test negative log-likelihood (NLL) based on 5 trials with standard error. For the out-of-distribution 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 et al., 2014) as 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 et al., 2020b), SWAG based on cyclic learning rates of 10 cycles (cycSWAG) (Maddox et al., 2019) and variance-reduced replica exchange SGHMC (reSGHMC) (Deng et al., 2021a). For a fair comparison, ICSGLD also conducts variance reduction on the energy function to alleviate the bias. Moreover, a large ζ = 3× 106 is selected, which only induces mild gradient multipliers ranging from −1 to 2 to penalize over-visited partitions. We don’t include SVGD (Liu & Wang, 2016) and SPOS (Zhang et al., 2020a) for scalability reasons. A batch size of 256 is selected. We run 4 parallel processes (×P4) with 500 epochs for M-SGD, reSGHMC and ICSGHMC and run cycSGHMC and cycSWAG 2000 epochs (×T4) based on a single process with 10 cycles. Refer to section D.4 of the supplementary material for the detailed settings.
Table 1 shows that the vanilla ensemble results via M-SGD×P4 surprisingly outperform cycSGHMC×T4 and cycSWAG×T4 on medium models, such as ResNet20 and ResNet32, and show very good performance on the out-of-distribution samples in Brier scores. We suspect that the parallel implementation (×P4) provides isolated initializations with less correlated samples; by contrast, cycSGHMC×T4 and cycSWAG×T4 explore the energy landscape contiguously, implying a risk to stay near the original region. reSGHMC×P4 shows a remarkable performance overall, but demonstrates a 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 cycSWAG is more independent of hyperparameter tuning, thus leading to better performance in larger models. We don’t report CSGHMC×P4 since it becomes quite unstable during the training of ResNet56 and WRN-16-8 models and causes mediocre results. As to ICSGHMC×P4, it consistently performs remarkable in both ACC and NLL and performs comparable to M-SGD×P4 in Brier scores. Code is available at github.com/WayneDW/Interacting-Contour-Stochastic-Gradient-Langevin-Dynamics.

heading: 6 CONCLUSION
text: 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 to big data problems without scalability concerns. ICSGLD has been compared with numerous state-of-the-art baselines for various tasks, whose remarkable results indicate its promising future in big data applications.
*The Brier score measures the mean squared error between the predictive and actual probabilities.

heading: ACKNOWLEDGMENT
text: Liang’s research was supported in part by the grants DMS-2015498, R01-GM117597 and R01GM126089. Lin acknowledges the support from NSF (DMS-1555072, DMS-2053746, and DMS2134209), BNL Subcontract 382247, and DE-SC0021142.

heading: A PRELIMINARIES

heading: A.1 STOCHASTIC APPROXIMATION
text: Given a random-field function H̃(θ,x), the stochastic approximation algorithm (Benveniste et al., 1990) proposes to solve the mean-field equation h(θ) = 0 in the analysis of adaptive algorithms
h(θ) = ∫ X H̃(θ,x)$θ(dx) = 0,
where x ∈ X ⊂ Rd, θ ∈ Θ ⊂ Rm, $θ(x) is a distribution that depends on the self-adapting parameter θ. Given the transition kernel Πθ(x, A) for any Borel subset A ⊂ X , the algorithm can be written as follows
(1) Simulate xk+1 ∼ Πθk(xk, ·), which yields the invariant distribution $θk(·),
(2) Optimize θk+1 = θk + ωk+1H̃(θk,xk+1).
Compared with the standard Robbins–Monro algorithm (Robbins & Monro, 1951), the algorithm proposes to simulate x from a transition kernel Πθ(·, ·) instead of the distribution $θ(·) directly. In other words, , H̃(θk,xk+1)− h(θk) is not a Martingale but rather a Markov state-dependent noise.

heading: A.2 POISSON’S EQUATION
text: In the stochastic approximation algorithm, the sequence of {(xk,θk)}∞k=1 on the product space X ×Θ is generated, which is an inhomogeneous Markov chain and requires the tool of the Poisson’s equation to study the convergence
µθ(x)−Πθµθ(x) = H̃(θ,x)− h(θ),
where µθ(·) is a function on X . The solution µθ(x) to the Poisson’s equation exists and is formulated in the form µθ(x) := ∑ k≥0 Πkθ(H̃(θ,x)− h(θ)),
when the above series converges. To ensure such a convergence, Benveniste et al. (1990) made the following regularity conditions on the solution µθ(·) of the Poisson’s equation: There exist a Lyapunov function V : X → [1,∞) and a positive constantC > 0 such that ∀θ,θ′ ∈ Θ, we have
‖Πθµθ(x)‖ ≤ CV (x), ‖Πθµθ(x)−Πθ′µθ′(x)‖ ≤ C‖θ − θ′‖V (x), E[V (x)] ≤ ∞, (12)
where a common choice for the Lyapunov function is to set V (x) = 1 + ‖x‖2 (Teh et al., 2016; Vollmer et al., 2016).

heading: A.3 GAUSSIAN DIFFUSIONS
text: Consider a stochastic linear differential equation
dUt = hθ(θt)Utdt+R 1/2(θt)dWt, (13)
where U is a m-dimensional random vector, hθ := ddθh(θ), R(θ) :=∑∞ k=−∞ Covθ(H(θ,xk), H(θ,x0)) is a positive definite matrix that depends on θ(·), W ∈ Rm
is a standard Brownian motion. Given a large enough t such that θt converges to a fixed point θ̂? sufficiently fast, we may write the diffusion associated with Eq.(13) as follow
Ut ≈ e−thθ(θ̂?)U0 + ∫ t
0
e−(t−s)hθ(θ̂?) ◦R(θ̂?)dWs, (14)
Suppose that the matrix hθ(θ̂?) is negative definite, then Ut converges in distribution to a Gaussian variable
E[Ut] = e−thθ(θ̂?)U0 Var(Ut) = ∫ t
0
ethθ(θ̂?) ◦R ◦ ethθ(θ̂?)du.
The main goal of this supplementary file is to study the Gaussian approximation of the process ω −1/2 k (θk − θ̂?) to the solution Eq.(14) for a proper step size ωk. Thereafter, the advantage of interacting mechanisms can be naturally derived.

heading: B STABILITY AND CONVERGENCE ANALYSIS
text: As required by the algorithm, we update P contour stochastic gradient Langevin dynamics (CSGLD) simultaneously. For the notations, we denote the particle of the p-th chain at iteration k by x(p)k ∈ X ⊂
Rd and the joint state of the P parallel particles at iteration k by x ⊗ P k := ( x (1) k ,x (2) k , · · · ,x (P ) k )> ∈
X ⊗ P ⊂ RdP . We also denote the learning rate and step size at iteration k by k and ωk, respectively. We denote by N (0, IdP ) a standard dP -dimensional Gaussian vector and denote by ζ a positive hyperparameter.
B.1 ICSGLD ALGORITHM
First, we introduce the interacting contour stochastic gradient Langevin dynamics (ICSGLD) with P parallel chains:
(1) Simulate x ⊗ P k+1 = x ⊗ P k − k∇xL̃(x ⊗ P k ,θk) +N (0, 2 kτIdP ), (S1)
(2) Optimize θk+1 = θk + ωk+1H̃(θk,x ⊗ P k+1 ), (S2)
where ∇xL̃(x ⊗ P ,θ) := ( ∇xL̃(x(1),θ),∇xL̃(x(2),θ), · · · ,∇xL̃(x(P ),θ) )> , ∇xL̃(x,θ) is the
stochastic adaptive gradient given by
∇xL̃(x,θ) = N
n
[ 1 + ζτ
∆u
( log θ(JŨ (x))− log θ((JŨ (x)− 1) ∨ 1) )] ︸ ︷︷ ︸
gradient multiplier
∇xŨ(x). (15)
In particular, the interacting random-field function is written as
H̃(θk,x ⊗ P k+1 ) = 1
P P∑ p=1 H̃(θk,x (p) k+1), (16)
where each random-field function H̃(θ,x) = (H̃1(θ,x), . . . , H̃m(θ,x)) follows H̃i(θ,x) = θ(JŨ (x)) ( 1i=JŨ (x) − θ(i) ) , i = 1, 2, . . . ,m. (17)
Here JŨ (x) denotes the index i ∈ {1, 2, 3, · · · ,m} such that ui−1 < N n Ũ(x) ≤ ui for a set of energy partitions {ui}mi=0 and Ũ(x) = ∑ i∈B Ui(x) where Ui denotes the negative log of a posterior based on a single data point i and B denotes a mini-batch of data of size n. Note that the stochastic
energy estimator Ũ(x) results in a biased estimation for the partition index JŨ (x) due to a non-linear transformation. To avoid such a bias asymptotically with respect to the learning rate k, we may consider a variance-reduced energy estimator ŨVR(x) following Deng et al. (2021a)
N n ŨVR(x) = N n ∑ i∈Bk ( Ui(x)− Ui ( xqb kq c )) + N∑ i=1 Ui ( xqb kq c ) , (18)
where the control variate xqb kq c is updated every q iterations.
Compared with the naı̈ve parallelism of CSGLD, a key feature of the ICSGLD algorithm lies in the joint estimation of the interacting random-field function H̃(θ,x ⊗ P ) in Eq.(16) for the same mean-field function h(θ).

heading: B.1.1 DISCUSSIONS ON THE HYPERPARAMETERS
text: The most important hyperparameter is ζ. A fine-tuned ζ usually leads to a small or even slightly negative learning rate in low energy regions to avoid local-trap problems. Theoretically, ζ affects the L2 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 to the energy into 200 subregions equally without tuning; since the optimization of SA is nearly convex, tuning {ωk} is much easier than tuning { k} for non-convex learning.

heading: B.1.2 DISCUSSIONS ON DISTRIBUTED COMPUTING AND COMMUNICATION COST
text: In shared-memory settings, the implementation is trivial and the details are omitted.
In distributed-memory settings: θk+1 is updated by the central node as follows:
• The p-th worker conducts the sampling step (S1) and sends the indices JŨ(x(p)k+1) ’s to the central
node;
• The central node aggregates the indices from all worker and updates θk based on (S2);
• The central node sends θk+1 back to each worker.
We emphasize that we don’t communicate the model parameters x ∈ Rd, but rather share the selfadapting parameter θ ∈ Rm, where m d. For example, WRN-16-8 has 11 M parameters (40 MB), while θ can be set to dimension 200 of size 4 KB; hence, the communication cost is not a big issue. Moreover, the theoretical advantage still holds if the communication frequency is slightly reduced.

heading: B.1.3 SCALABILITY TO BIG DATA
text: Recall that the adaptive sampler follows that
xk+1 = xk − k+1 N
n
[ 1 + ζτ
log θk(JŨ (xk))− log θk((JŨ (xk)− 1) ∨ 1) ∆u ] ︸ ︷︷ ︸
gradient multiplier
∇xŨ(xk) + √ 2τ k+1wk+1,
The key to the success of (I)CSGLD is to generate sufficiently strong bouncy moves (negative gradient multiplier) to escape local traps. To this end, ζ can be tuned to generate proper bouncy moves.
Take the CIFAR100 experiments for example:
• the self-adjusting mechanism fails if the gradient multiplier uniformly “equals” to 1 and a too small value of ζ = 1 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.
However, when we set ζ =3e6, the original stochastic approximation (SA) update proposed in (Deng et al., 2020b) follows that
θk+1(i) = θk(i) + ωk+1 θ ζ k(JŨ (xk+1))︸ ︷︷ ︸
essentially 0 for ζ 1
( 1i=JŨ (xk+1) − θk(i) ) .
Since θ(i) < 1 for any i ∈ {1, · · · ,m}, θ(i)ζ is essentially 0 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
θk+1(i) = θk(i) + ωk+1 θk(JŨ (xk+1))︸ ︷︷ ︸ independent of ζ
( 1i=JŨ (xk+1) − θk(i) )
is more independent of ζ and proposes to converge to a much smoother equilibrium θ1/ζ∞ instead of θ∞, where θ∞(i) = ∫ χi π(x)dx ∝ ∫ χi e− U(x) τ dx is the energy PDF. As such, despite the linear stability is sacrificed, the resulting algorithm is more scalable. For example, estimating e−10,000× 1 ζ is numerically much easier than e−10,000 for a large ζ such as 10, 000, where 10, 000 can be induced by the high losses in training deep neural networks in big data.

heading: B.2 ASSUMPTIONS
text: 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 Θ on a fixed compact set.
Assumption A1 (Compactness) The space Θ is compact and for any i ∈ {1, 2, . . . ,m} we have infΘ θ(i) > 0. In addition, there exists a positive constant Q > 0 that satisfies ∀θ ∈ Θ and x ∈ X ,
‖θ‖ ≤ Q, ‖H̃(θ,x)‖ ≤ Q. (19)
For weaker assumptions, we refer readers to Theorem 3.2 (Fort et al., 2015), where a recurrence property can be proved for the Metropolis-based Wang-Landau 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 et al. (2018).
Assumption A2 (Smoothness) U(x) is M -smooth when there exists a positive constant M that satisfies ∀x,x′ ∈ X ,
‖∇xU(x)−∇xU(x′)‖ ≤M‖x− x′‖. (20)
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 et al. (2017) and Xu et al. (2018).
Assumption A3 (Dissipativity) There exist constants m̃ > 0 and b̃ ≥ 0 that satisfies ∀x ∈ X and any θ ∈ Θ,
〈∇xL(x,θ),x〉 ≤ b̃− m̃‖x‖2. (21)
To further establish a bounded second moment on x ∈ X with respect to a proper Lyapunov function V (x), we impose the following conditions on the gradient noise:
Assumption A4 (Gradient noise) The stochastic gradient based on mini-batch settings is an unbiased estimator such that E[∇xŨ(xk)−∇xU(xk)] = 0; furthermore, for some positive constants M and B, we have
E[‖∇xŨ(xk)−∇xU(xk)‖2] ≤M2‖x‖2 +B2,
where E[·] acts on the distribution of the noise in the stochastic gradient∇xŨ(xk).

heading: B.3 LOCAL STABILITY VIA THE SCALABLE RANDOM-FIELD FUNCTION
text: Now, we are ready to present our first result. Lemma 3 establishes a local stability condition for the non-linear mean-field system of ICSGLD, which implies a potential convergence of θk to a unique fixed point that adapts to a wide energy range under mild assumptions.
Lemma 3 (Local stability, restatement of Lemma 1) Assume Assumptions A1-A4 hold. Given any small enough learning rate , a large enough m and batch size n, and any θ ∈ Θ̃, where Θ̃ is a small neighborhood of θ? that contains θ̂?, we have 〈h(θ),θ−θ̂?〉 ≤ −φ‖θ−θ̂?‖2, where θ̂? = θ?+O(ε),
ε = O ( supx Var(ξn(x)) + + 1 m ) and θ? =
( (∫ X1 π(x)dx ) 1 ζ
∑m k=1 (∫ Xk π(x)dx ) 1 ζ , . . . ,
( ∫ Xm
π(x)dx) 1 ζ∑m
k=1 (∫ Xk π(x)dx ) 1 ζ
) ,
φ = infθ mini Ẑ −1 ζ,θ(i)
( 1−O(ε) ) > 0, Ẑζ,θ(i) is defined below Eq.(29), and ξn(x) denotes the noise
in the energy estimator Ũ(x) of batch size n and Var(·) denotes the variance. Proof The random-field function H̃i(θ,x) = θ(JŨ (x)) ( 1i=JŨ (x) − θ(i) ) based on the stochastic
energy estimator Ũ(x) yields a biased estimator of Hi(θ,x) = θ(J(x)) ( 1i=J(x) − θ(i) ) for any i ∈ {1, 2, . . . ,m} based on the exact energy partition function J(·). By Lemma.4, we know that the bias caused by the stochastic energy is of order O(Var(ξn(x))). Now we compute the mean-field function h(θ) based on the measure $θ(x) simulated from SGLD:
hi(θ) = ∫ X H̃i(θ,x)$θ(x)dx = ∫ X Hi(θ,x)$θ(x)dx +O (Var(ξn(x)))
= ∫ X Hi(θ,x) $Ψ̃θ (x)︸ ︷︷ ︸ I1 −$Ψ̃θ (x) +$Ψθ (x)︸ ︷︷ ︸ I2:piece-wise approximation −$Ψθ (x) +$θ(x)︸ ︷︷ ︸ I3:numerical discretization  dx +O (Var(ξn(x))) , (22)
where $θ is the invariant measure simulated via SGLD that approximates $Ψθ (x). $Ψθ (x) and $Ψ̃θ (x) are two invariant measures that follow $Ψθ (x) ∝ π(x)
Ψζθ(U(x)) and $Ψ̃θ (x) ∝
π(x)
Ψ̃ζθ(U(x)) ;
Ψθ(u) and Ψ̃θ(u) are piecewise continuous and constant functions, respectively Ψθ(u) = m∑ k=1 ( θ(k − 1)e(log θ(k)−log θ(k−1)) u−uk−1 ∆u ) 1uk−1<u≤uk ; Ψ̃θ(u) = m∑ k=1 θ(k)1uk−1<u≤uk .
(23)
(i) For the first term I1, we have∫ X Hi(θ,x)$Ψ̃θ (x)dx = 1 Z̃ζ+1,θ ∫ X θ(J(x)) ( 1i=J(x) − θ(i) ) π(x) θζ(J(x)) dx
= 1
Z̃ζ+1,θ m∑ k=1 ∫ Xk (1i=k − θ(i)) π(x) θζ−1(k) dx
= 1
Z̃ζ+1,θ [ m∑ k=1 ∫ Xk π(x) θζ−1(k) 1k=idx− θ(i) m∑ k=1 ∫ Xk π(x) θζ−1(k) dx ]
= 1
Z̃ζ+1,θ
[∫ Xi π(x)dx
θζ−1(i) − θ(i)Z̃ζ,θ
] ,
(24)
where Z̃ζ+1,θ = ∑m k=1 ∫ Xk π(x)dx θζ(k) denotes the normalizing constant of $Ψ̃θ (x).
The solution θ? that solves ∫ Xk π(x)dx
θζ−1(k) − θ(k)Z̃ζ,θ = 0 for any k ∈ {1, 2, · · · ,m} satisfies θ?(k) =( ∫
Xk π(x)dx
Z̃ζ,θ?
) 1 ζ
. Combining the definition of Z̃ζ,θ? = ∑m k=1 ∫ Xk π(x)dx
θζ−1? (k) , we have
Z̃ζ,θ? = m∑ k=1
∫ Xk π(x)dx
θζ−1? (k) = m∑ k=1 ∫ Xk π(x)dx( ∫ Xk π(x)dx
Z̃ζ,θ?
) ζ−1 ζ
= Z̃ ζ−1 ζ
ζ,θ? m∑ k=1 ∫ Xk π(x)dx(∫ Xk π(x)dx ) ζ−1 ζ = Z̃ ζ−1 ζ ζ,θ? m∑ k=1 (∫ Xk π(x)dx ) 1 ζ ,
which leads to Z̃ζ,θ? = (∑m k=1 (∫ Xk π(x)dx ) 1 ζ )ζ . In other words, the mean-field system without
perturbations yields a unique solution θ?(i) =
(∫ Xi π(x)dx ) 1 ζ
∑m k=1 (∫ Xk π(x)dx ) 1 ζ for any i ∈ {1, 2, · · · ,m}.
(ii) For the second term I2, we have∫ X Hi(θ,x)(−$Ψ̃θ (x) +$Ψθ (x))dx = O ( 1 m ) , (25)
where the result follows from the boundedness of H(θ,x) in (A1) and Lemma B4 (Deng et al., 2020b). (iii) For the last term I3, following Theorem 6 of Sato & Nakagawa (2014), we have for any fixed θ,∫ X Hi(θ,x) (−$Ψθ (x) +$θ(x)) dx = O( ). (26)
Plugging Eq.(24), Eq.(25) and Eq.(26) into Eq.(22), we have
hi(θ) = Z̃ −1 ζ+1,θ [ εβ̃i(θ) + ∫ Xi π(x)dx
θζ−1(i) − θ(i)Z̃ζ,θ
]
= Z̃−1ζ+1,θ Z̃ζ,θ? θζ−1(i)
[ εβ̃i(θ) θζ−1(i)
Z̃ζ,θ? +
∫ Xi π(x)dx
Z̃ζ,θ? − θζ(i) Z̃ζ,θ Z̃ζ,θ?
]
= Z̃−1ζ+1,θ Z̃ζ,θ? θζ−1(i)
[ εβ̃i(θ) θζ−1(i)
Z̃ζ,θ? + θζ?(i)− (θ(i)Cθ)
ζ ] ,
(27)
where β̃i(θ) is a bounded term such that Z̃−1ζ+1,θεβ̃i(θ) = O ( Var(ξn(x)) + + 1 m ) , ε =
O ( supx Var(ξn(x)) + + 1 m ) and Cθ = ( Z̃ζ,θ
Z̃ζ,θ?
) 1 ζ . By the definition of Z̃ζ,θ = ∑m k=1 ∫ Xk π(x)dx θζ−1(k) ,
when ζ = 1, Cθ ≡ 1 for any θ ∈ Θ, which suggests that the stability condition doesn’t rely on the initialization of θ; however, when ζ 6= 1, Cθ 6= 1 when θ 6= θ?, we see that hi(θ) ∝ θ?(i)ζ − (θ(i)Cθ)ζ + perturbations is a non-linear mean-field system and requires a proper initialization of θ ∈ Θ̃.
For any θ ∈ Θ̃ ⊂ Θ being close enough to θ?, there exists a Lipschitz constant Lθ̃ = supi≤m,θ∈Θ̃ |Cθ?−Cθ| |θ?(i)−θ(i)| < ∞. By Cθ? = 1, θ(i) ≤ 1, and mean value theorem for some θ̃(i) ∈ [θ(i), θ?(i)], we have
|θζ?(i)− (θ(i)Cθ) ζ | = ζ(θ̃(i)Cθ̃) ζ−1|θ?(i)− θ(i)Cθ|
= ζ(θ̃(i)Cθ̃) ζ−1|θ?(i)− θ(i) + θ(i)Cθ? − θ(i)Cθ| ≤ ζ(θ̃(i)Cθ̃) ζ−1|θ?(i)− θ(i)|+ θ(i)|Cθ? − Cθ| ≤ ζ(θ̃(i)Cθ̃) ζ−1(1 + Lθ̃)|θ?(i)− θ(i)|, (28)
Combining Eq.(27) and Eq.(28), we have
hi(θ) = Z̃ −1 ζ+1,θ Z̃ζ,θ? θζ−1(i)
[ εβ̃i(θ) θζ−1(i)
Z̃ζ,θ? + θζ?(i)− (θ(i)Cθ)
ζ ] = Ẑ−1ζ,θ(i) [εβi(θ) + θ?(i)− θ(i)] ,
(29)
where Ẑ−1ζ,θ(i) = Z̃−1ζ+1,θZ̃ζ,θ?
ζ(θ̃(i)C θ̃ )ζ−1(1+L θ̃ )θζ−1(i)
; βi(θ) is some bounded term such that βi(θ) ≤
β̃i(θ)θ ζ−1(i)
ζ(θ̃(i)C θ̃ )ζ−1(1+L θ̃ )Z̃ζ,θ?
; Cθ̃ = ( Z̃ ζ,θ̃
Z̃ζ,θ?
) 1 ζ
; Lθ̃ = supi≤m,θ∈Θ̃ |Cθ?−Cθ| |θ?(i)−θ(i)| <∞.
Next, we apply the perturbation theory to solve the ODE system with small disturbances (Weinhart et al., 2010) and obtain the equilibrium θ̂?,
where εβ(θ̂?) + θ? − θ̂? = 0, to the mean-field equation hi(θ) such that
hi(θ) = Ẑ −1 ζ,θ(i) [εβi(θ) + θ?(i)− θ(i)]
= Ẑ−1ζ,θ(i) [ εβi(θ)− εβi(θ̂?) + εβi(θ̂?) + θ?(i)− θ(i) ] = Ẑ−1ζ,θ(i) [ O(ε)(θ(i)− θ̂?(i)) + θ̂?(i)− θ(i)
] = Ẑ−1ζ,θ(i) ( 1−O(ε) ) ( θ̂?(i)− θ(i) ) ,
(30)
where a smoothness condition clearly holds for the β(·) function. Given a positive definite Lyapunov function V(θ) = 12‖θ̂?−θ‖ 2, the mean-field system h(θ) = Ẑ−1ζ,θ(i)(εβ(θ)+θ?−θ) = Ẑ −1 ζ,θ(i)(1−
O(ε))(θ̂? − θ) for i ∈ {1, 2, · · · ,m} enjoys the following property
〈h(θ),∇V(θ)〉 = 〈h(θ),θ − θ̂?〉
≤ −min i Ẑ−1ζ,θ(i)
( 1−O(ε) ) ‖θ − θ̂?‖2
≤ −φ‖θ − θ̂?‖2,
where φ = infθ mini Ẑ−1ζ,θ(i) ( 1 − O(ε) ) > 0 given the compactness assumption A1 and a small
enough ε = O ( supx Var(ξn(x)) + + 1 m ) .
Remark 1 The newly proposed random-field function Eq.(17) may sacrifice the global stability by including an approximately linear mean-field system Eq.(29) instead of a linear stable system (see formula (15) in Deng et al. (2020b)). The advantage, however, is that such a mechanism facilitates the estimation of θ?. We emphasize that the original energy probability in each partition{ ∫ Xk π(x)dx }m k=1 (Deng et al., 2020b) may be very difficult to estimate for big data problems. By
contrast, the estimation of {( ∫ Xk π(x)dx ) 1 ζ }m k=1 becomes much easier given a proper ζ > 0.

heading: Technical lemmas
text: Lemma 4 The stochastic energy estimator Ũ(x) leads to a controllable bias in the random-field function.
|E[H̃i(θ,x)]−Hi(θ,x)| = O (Var(ξn(x))) , where the expectation E[·] is taken with respect to 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̃i(θ,x) = θ(JŨ (x)) ( 1i=JŨ (x) − θ(i) ) and JŨ (x) ∈ {1, 2, · · · ,m} satisfies
uJŨ (x)−1 < N n Ũ(x) ≤ uJŨ (x) for a set of energy partitions {ui} m i=0. We can interpret H̃i(θ,x) as a
non-linear transformation Φ that maps Ũ(x) to (0, 1). Similarly, Hi(θ,x) maps U(x) to (0, 1). In what follows, the bias of random-field function is upper bounded as follows
|E[H̃i(θ,x)]−Hi(θ,x)| = ∣∣∣∣∫ Φ(U(x) + ξ(x))− Φ(U(x))dµ(ξ(x))∣∣∣∣
= ∣∣∣∣∫ ξ(x)Φ′(U(x)) + ξ(x)22 Φ′′(u)dµ(ξ(x)) ∣∣∣∣
= O (Var(ξn(x))) ,
where the second equality follows from Taylor expansion for some energy u and the third equality follows because the stochastic energy estimator is unbiased; Φ′(U(x)) = O( θ(J(x))−θ(J(x)−1)∆u ) is clearly bounded due to the definition of θ; a similar conclusion also applies to Φ′′(·).

heading: B.4 CONVERGENCE OF THE SELF-ADAPTING PARAMETERS
text: The following is a restatement of Lemma 3.2 of Raginsky et al. (2017), which holds for any θ in the compact space Θ.
Lemma 5 (Uniform L2 bounds) Assume Assumptions A1, A3 and A4 hold. We have a bounded second moment supk≥1 E[‖xk‖2] <∞ given a small enough learning rate.
The following lemma justifies the regularity properties of Poisson’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 a more detailed proof by utilizing a Lyapunov function V (x) = 1 + x2 and Lemma 5.
Lemma 6 (Solution of Poisson’s equation) Assume that Assumptions A1-A4 hold. There is a solution µθ(·) on X to the Poisson’s equation
µθ(x)−Πθµθ(x) = H̃(θ,x)− h(θ). (31)
Furthermore, there exists a constant C such that for all θ,θ′ ∈ Θ
E[‖Πθµθ(x)‖] ≤ C, E[‖Πθµθ(x)−Πθ′µθ′(x)‖] ≤ C‖θ − θ′‖.
(32)
Proof 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).
(DRI) By the smoothness assumption A2, we have that U(x) is continuously differentiable almost everywhere. By the dissipative assumption A3 and Theorem 2.1 (Roberts & Tweedie, 1996), we can show that the discrete dynamics system is irreducible and aperiodic. Now consider a Lyapunov function V = 1 + ‖x‖2 and any compact subset K ⊂ Θ, the drift conditions are verified as follows: (DRI1) Given small enough learning rates { k}k≥1, the smoothness assumption A2, and the dissipative assumption A3, applying Corollary 7.5 (Mattingly et al., 2002) yields the minorization condition for the CSGLD algorithm, i.e. there exists η > 0, a measure ν, and a set C such that ν(C) = 1. Moreover, we have
Pθ∈K(x,A) ≥ ην(A) ∀A ∈ X ,x ∈ C. (I)
where Pθ(x,y) := 1 2 √ (4π )d/2 E [ e− ‖y−x+ ∇xL̃(x,θ)‖2 4 |x ] denotes the transition kernel based on CSGLD with the parameter θ ∈ K and a learning rate , in addition, the expectation is taken over the adaptive gradient ∇xL̃(x,θ) in Eq.(15). Using Assumption A1-A4, we 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 α̃ ∈ (0, 1) and β̃ ≥ 0 such that
Pθ∈KV (x) ≤ α̃V (x) + β̃. (II)
Consider a Lyapunov function V = 1 + ‖x‖2 and a constant κ = α̃+ β̃, it yields that
Pθ∈KV (x) ≤ κV (x). (III)
Now we have verified the first condition (DRI1) by checking conditions (I),(II), and (III),
(DRI2) In what follows, we check the boundedness and Lipshitz conditions on the random-field function H̃(θ,x), where each subcomponent is defiend as H̃i(θ,x) = θ(JŨ (x)) ( 1i=JŨ (x) − θ(i) ) . Recall that V = 1 + ‖x‖2, the compactness assumption A1 directly leads to
sup θ∈K⊂[0,1]m
‖H(θ,x)‖ ≤ mV (x). (IV)
For any θ1,θ2 ∈ K and a fixed x ∈ X , it suffices for us to solely verify the i-th index, which is the index that maximizes |θ1(i)− θ2(i)|, then
|H̃i(θ1,x)− H̃i(θ2,x)| = θ1(JŨ (x)) ( 1i=JŨ (x) − θ1(i) ) − θ2(JŨ (x)) ( 1i=JŨ (x) − θ2(i) ) ≤ |θ1(JŨ (x))− θ2(JŨ (x))|+ |θ1(JŨ (x))θ1(i)− θ2(JŨ (x))θ2(i)|
≤ max j
( |θ1(j)− θ2(j)|+ θ1(j)|θ1(i)− θ2(i)|+ |θ1(j)− θ2(j)|θ2(i) ) ≤ 3|θ1(i)− θ2(i)|,
where the last inequality holds since θ(i) ∈ (0, 1] for any i ≤ m. (DRI3) We proceed to verify the smoothness of the transitional kernel Pθ(x,y) with respect to θ. For any θ1,θ2 ∈ K and fixed x and y, we have
|Pθ1(x,y)− Pθ2(x,y)|
= 1 2 √ (4π )d/2 E [ e− ‖y−x+ ∇xL̃(x,θ1)‖ 2 4 |x ] − 1 2 √ (4π )d/2 E [ e− ‖y−x+ ∇xL̃(x,θ2)‖ 2 4 |x ]
. |‖y − x + ∇xL̃(x,θ1)‖2 − ‖y − x + ∇xL̃(x,θ2)‖2|
. ‖∇xL̃(x,θ1)−∇xL̃(x,θ2)‖
. ‖θ1 − θ2‖,
where the first inequality (up to a finite constant) follows by ‖ex − ey‖ . ‖x− y‖ for any x, y in a compact space; the last inequality follows by the definition of the adaptive gradient in Eq.(15) and ‖ log(x)− log(y)‖ . ‖x− y‖ by the compactness assumption A1.
For f : X → Rd, define the norm ‖f‖V = supx∈X |f(x)| V (x) . Following the same technique proposed in Liang et al. (2007) (page 319), we can verify the last drift condition
‖Pθ1f − Pθ2f‖V ≤ C‖f‖V ‖θ1 − θ2‖, ∀f ∈ LV := {f : X → Rd, ‖f‖V <∞}. (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 L2 convergence of θk, we make some extra assumptions on the step size.
Assumption A5 (Learning rate and step size) The learning rate { k}k∈N is a positive nonincreasing sequence of real numbers satisfying the conditions
lim k k = 0, ∞∑ k=1 k =∞.
The step size {ωk}k∈N is a positive non-increasing sequence of real numbers such that
lim k→∞ ωk = 0, ∞∑ k=1 ωk = +∞, ∞∑ k=1 ω2k < +∞. (33)
A practical strategy is to set ωk := O(k−α) to satisfy the above conditions for any α ∈ (0.5, 1].
The following is an application of Theorem 24 (page 246) (Benveniste et al., 1990) given stability conditions (Lemma 3).
Lemma 7 (L2 convergence rate, restatement of Lemma 2) Assume Assumptions A1-A5 hold. For any θ0 ∈ Θ̃ ⊂ Θ, a large m, small learning rates { k}∞k=1, and step sizes {ωk}∞k=1, {θk}∞k=0 converges to θ̂?, where θ̂? = θ? +O ( supx Var(ξn(x)) + supk≥k0 k + 1 m ) for some k0, such that
E [ ‖θk − θ̂?‖2 ] = O (ωk) .
The theoretical novelty is that we treat the biased θ̂? as the equilibrium of the continuous system instead of analyzing how far we are away from θ? in all aspects as in Theorem 1 (Deng et al., 2020b). This enables us to directly apply Theorem 24 (page 246). Nevertheless, it can be interpreted as a special case of Theorem 1 (Deng et al., 2020b) except that there are no perturbation terms and the equilibrium is θ̂? instead of θ?.

heading: C GAUSSIAN APPROXIMATION

heading: C.1 PRELIMINARY: SUFFICIENT CONDITIONS FOR WEAK CONVERGENCE
text: To formally prove the asymptotic normality of the stochastic approximation process ω−1/2k (θk − θ̂?), we first lay out a preliminary result (Theorem 1 of Pelletier (1998)) that provides sufficient conditions to guarantee the weak convergence.
Lemma 8 (Sufficient Conditions) Consider a stochastic algorithm as follows θk+1 = θk + ωk+1h(θk) + ωk+1ν̃k+1 + ωk+1ek+1,
where ν̃k+1 denotes a perturbation and ek+1 is a random noise. Given three conditions (C1), (C2), and (C3) defined below, we have the desired weak convergence result
ω− 1 2 (θk − θ̂?)⇒ N (0,Σ), (34) where Σ = ∫∞
0 ethθ? ◦ R ◦ eth
> θ?dt, R denotes the limiting covariance of the martingale
limk→∞ E[ek+1e>k+1|Fk] andFk is the σ-algebra of the events up to iteration k, hθ? = hθ(θ̂?)+ξ̂I , ξ̂ = limk→∞ ω0.5k −ω 0.5 k+1
ω1.5k . †
(C1) There exists an equilibrium point θ̂? and a stable matrix hθ? := hθ(θ̂?) ∈ Rm×m such that for any θ ∈ {θ : ‖θ − θ̂?‖ ≤ M̃} for some M̃ > 0, the mean-field function h : Rm → Rm satisfies
h(θ̂?) = 0
‖h(θ)− hθ?(θ − θ̂?)‖ . ‖θ − θ̂?‖2,
(C2) The step size ωk decays with an order α ∈ (0, 1] such that ωk = O(k−α).
(C3) Assumptions on the disturbances . There exists constants M̃ > 0 and α̃ > 2 such that
E [ek+1|Fk] 1{‖θ−θ̂?‖≤M̃} = 0, (I1)
sup k
E [ ‖ek+1‖α̃|Fk ] 1{‖θ−θ̂?‖≤M̃} <∞, (I2)
E [ ω−1k ‖ν̃k+1‖ 2 ] 1{‖θ−θ̂?‖≤M̃} → 0, (II)
E [ ek+1e > k+1|Fk ] 1{‖θ−θ̂?‖≤M̃} → R. (III)
Remark 2 By the definition of the mean-field function h(θ) in Eq.(27), it is easy to verify the condition C1. Moreover, Assumption A5 also fulfills the condition C2. Then, the proof hinges on the verification of the condition C3.
†For example, ξ̂ = 0 if ωk = O(k−α), where α ∈ (0.5, 1] and ξ̂ = k02 if ωk = k0 k .

heading: C.2 PRELIMINARY: CONVERGENCE OF THE COVARIANCE ESTIMATORS
text: In particular, to verify the condition E [ ek+1e > k+1|Fk ] 1{‖θ−θ̂?‖≤M̃} → R, , we study the convergence of the empirical sample mean E[f(xk)] for a test function f to the posterior expectation f̄ = ∫ X f(x)$θ̂?(x)(dx). Poisson’s equation is often used to characterize the fluctuation between f(x) and f̄ :
Lg(x) = f(x)− f̄ , (35)
where L refers to an infinitesimal generator and g(x) denotes the solution of the Poisson’s equation. Similar to the proof of Lemma 6, the existence of the solution of the Poisson’s equation has been established in (Mattingly et al., 2002; Vollmer et al., 2016). Moreover, the perturbations of E[f(xk)]− f̄ are properly bounded given regularity properties for g(x), where the 0-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 et al., 2015) with decreasing learning rates { k}k≥1. The gradient biases from Theorem 2 (Chen et al., 2015) are also included to handle the adaptive biases.
Lemma 9 (Convergence of the Covariance Estimators) Suppose Assumptions A1-A5 hold. For any θ0 ∈ Θ̃ ⊂ Θ, a large m, small learning rates { k}∞k=1, step sizes {ωk}∞k=1 and any bounded function f , we have ∣∣∣∣E [f(xk)]− ∫
X f(x)$θ̂?(x)dx ∣∣∣∣→ 0, where $θ̂?(x) is the invariant measure simulated via SGLD that approximates $Ψ̃θ? (x) ∝ π(x) θζ?(J(x)) .
Proof We study the single-chain CSGLD and reformulate the adaptive algorithm as follows:
xk+1 = xk − k∇xL̃(xk,θk) +N (0, 2 kτI) = xk − k ( ∇xL̃(xk, θ̂?) + Υ(xk,θk) ) +N (0, 2 kτI),
where ∇xL̃(x,θ) = Nn [ 1 + ζτ∆u (log θ(J(x))− log θ((J(x)− 1) ∨ 1)) ] ∇xŨ(x) ‡, ∇xL̃(x,θ) is
defined in Section B.1 and the bias term is given by Υ(xk,θk) = ∇xL̃(xk,θk)−∇xL̃(xk, θ̂?). Then, by Jensen’s inequality and Lemma 7, we have
‖E[Υ(xk,θk)]‖ ≤ E[‖∇xL̃(xk,θk)−∇xL̃(xk, θ̂?)‖] . E[‖θk − θ̂?‖] ≤ √ E[‖θk − θ̂?‖2] ≤ O ( √ ωk) .
(36)
Combining Eq.(36) and Theorem 5 (Chen et al., 2015), we have
∣∣∣∣E [f(xk)]− ∫ X f(x)$θ̂?(x)dx ∣∣∣∣ = O ( 1∑k i i + ∑k i=1 ωi‖E[Υ(xi,θi)]‖∑k i ωi + ∑k i 2 i∑k i i ) → 0, as k →∞,
where the last argument directly follows from the conditions on learning rates and step sizes in Assumption A5.
‡J(x) = ∑m i=1 i1ui−1<U(x)≤ui , where the exact energy function U(x) is selected.

heading: C.3 PROOF OF THEOREM 1
text: Recall that the stochastic approximation based on a single process follows from
θk+1
= θk + ωk+1H(θk,xk+1)
= θk + ωk+1h(θk) + ωk+1 (µθk(xk+1)−Πθkµθk(xk+1)) = θk + ωk+1h(θk)
+ ωk+1 ( Πθk+1µθk+1(xk+1)−Πθkµθk(xk+1) +
ωk+2 − ωk+1 ωk+1 Πθk+1µθk+1(xk+1) ) ︸ ︷︷ ︸
νk+1
+ ωk+1
( 1
ωk+1
( ωk+1Πθkµθk(xk)− ωk+2Πθk+1µθk+1(xk+1) ) ︸ ︷︷ ︸
ςk+1
+µθk(xk+1)−Πθkµθk(xk)︸ ︷︷ ︸ ek+1
)
= θk + ωk+1h(θk) + ωk+1 (νk+1 + ςk+1)︸ ︷︷ ︸ perturbation +ωk+1 ek+1︸︷︷︸ martingale ,
(37)
where the second equality holds from the solution of Poisson’s equation in Eq.(31).
We denote θ̈k = θk + ωk+1Πθkµθk(xk). Adding ωk+2Πθk+1µθk+1(xk+1) on both sides of Eq.(37), we have θ̈k+1
= θ̈k + ωk+1h(θk) + ωk+1 (νk+1 + ek+1 + ςk+1) + ωk+2Πθk+1µθk+1(xk+1)− ωk+1Πθkµθk(xk) = θ̈k + ωk+1h(θk) + ωk+1 (νk+1 + ek+1)
= θ̈k + ωk+1h(θ̈k) + ωk+1 (ν̃k+1 + ek+1) , (38)
where ν̃k+1 = νk+1 + h(θk)− h(θ̈k). Next, we proceed to verify the conditions in C3. (I) By the martingale difference property of {ek} and the compactness assumption A1, we know that for any α̃ > 2
E[ek+1|Fk] = 0, sup k≥0 E[‖ek+1‖α̃|Fk] <∞. (I)
(II) By the definition of h(θk) in Eq.(27), we can easily check that h(θk) is Lipschitz continuous in a neighborhood of θ̂?. Combining Eq.(32), we have ‖h(θk) − h(θ̈k)‖ = O(‖θk − θ̈k‖) = O(‖ωk+1Πθkµθk(xk)‖) = O(ωk+1). Then E[‖νk+1‖] ≤ C‖θk − θ̈k‖+O(ωk+2) = O(ωk+1) by the step size condition Eq.(33). In what follows, we can verify
E [ ‖ν̃k+1‖2
ωk
] ≤ 2E [ ‖νk+1‖2
ωk
] + 2E [ ‖h(θk)− h(θ̈k)‖2
ωk
] = O(ωk)→ 0. (II)
(III) For the martingale difference noise ek+1 = µθk(xk+1)−Πθkµθk(xk) with mean 0, we have
E[ek+1e>k+1|Fk] = E[µθk(xk+1)µθk(xk+1)>|Fk]−Πθkµθk(xk)Πθkµθk(xk)>.
We denote E[ek+1e>k+1|Fk] by a function f(xk). Applying Lemma 9, we have E[ek+1e>k+1|Fk] = E[f(xk)]→ ∫ f(x)$θ̂?dx = limk→∞ E[ek+1e>k+1|Fk] := R, (III)
whereR := R(θ̂?) andR(θ) is also equivalent to ∑∞ k=−∞ Covθ(H(θ,xk), H(θ,x0)).
Having the conditions C1, C2 and C3 verified, we apply Lemma 8 and have the following weak convergence for θ̈k
ω −1/2 k (θ̈k − θ̂?)⇒ N (0,Σ),
where Σ = ∫∞
0 ethθ? ◦R ◦ eth > θ?dt and hθ? = hθ(θ̂?) + ξ̂I , ξ̂ = limk→∞ ω0.5k −ω 0.5 k+1
ω1.5k .
Considering the definition that θ̈k = θk + ωk+1Πθkµθk(xk) and E[‖Πθkµθk(xk)‖] is uniformly bounded by Eq.(32), we have
ω 1/2 k Πθkµθk(xk)→ 0 in probability.
By Slutsky’s theorem, we eventually have the desired result
ω −1/2 k (θk − θ̂?)⇒ N (0,Σ).
where the step size ωk decays with an order α ∈ (0.5, 1] such that ωk = O(k−α).

heading: D MORE ON EXPERIMENTS

heading: D.1 MODE EXPLORATION ON MNIST VIA THE SCALABLE RANDOM-FIELD FUNCTION
text: For the network structure, we follow Jarrett et al. (2009) and choose a standard convolutional neural network (CNN). Such a CNN has two convolutional (conv) layers and two 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. A large batch size of 2500 is selected to 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 ∆u = 10. The step size follows ωk = min{0.01, 1k0.6+100}.

heading: D.2 SIMULATIONS OF MULTI-MODAL DISTRIBUTIONS
text: The target density function is given by π(x) ∝ exp(−U(x)), where x = (x1, x2) and U(x) follows U(x) = 0.2(x21 + x 2 2) − 2(cos(2πx1) + cos(2πx2)).
Figure 5: Target density.
We also include a regularization term L(x) = I(x21+x22)>20 × ((x21 + x 2 2) − 20). This design leads to a highly multi-modal distribution with 25 isolated modes. Figure 5 shows the contour and the 3-D plot of the target density. The ICSGLD and baseline algorithms are applied to this example. For ICSGLD, we set k = 3e−3, τ = 1, ζ = 0.75 and total number of iterations= 8e4. Besides, we partition the sample space into 100 subregions with bandwidth ∆u = 0.125 and set ωk = min(3e −3, 1k0.6+100 ).
For comparison, we run the baseline algorithms under similar settings. For CSGLD, we run a single process 5 times of the time budget and all the settings are the same as those used by ICSGLD. For reSGLD, we run five 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−3 and a temperature of 1. For cycSGLD, we run a single-chain with 5 times of the time budget. We set the initial learning rate as 1e−2 and choose 10 cycles. For the particle-based SVGD, we run five chains in parallel. For each chain, we initialize 100 particles as being drawn from a uniform distribution over a rectangle. The learning rate is set to 3e−3.
To compare the convergence rates in terms of running steps and time between ICSGLD and other algorithms, we repeat each algorithm 20 times and calculate the mean and standard error over 20 trials. Note that we run all the algorithms based on 5 parallel chains (×P5) except that cycSGLD and CSGLD are run in a single-chain with 5 times of time budget (×T5) and the steps and running
time are also scaled accordingly. Figure 6 shows that the vanilla SGLD×P5 converges the slowest among the five algorithms due to the lack of mechanism to escape local traps; cycSGLD×T5 slightly alleviates that problem by adopting cyclical learning rates; reSGLD×P5 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; ICSGLD×P5 converges faster than all the others and the noisy energy estimators only induce a bias for the latent variables and don’t affect the convergence rate significantly.
For the particle-based SVGD method, since more particles require expensive computations while fewer particles lead to a crude approximation. Therefore, we don’t show the convergence of SVGD and only compare the Monte Carlo methods.

heading: D.3 DEEP CONTEXTUAL BANDITS ON MUSHROOM TASKS
text: For the UCI Mushroom data set, each mushroom is either edible or poisonous. Eating an edible mushroom yields a reward of 5, but eating a poisonous mushroom has a 50% chance to result in a reward of -35 and a reward of 5 otherwise. Eating nothing results in 0 reward. All the agents use the same architecture. In particular, we fit a two-layer neural network with 100 neurons each and ReLU activation functions. The input of the network is a feature vector with dimension 22 (context) and 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 to online scenarios, we train models after every 20 new observations.
We choose one -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 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 4 parallel chains with average outputs (×P4) except that pCSGLD runs a single-chain with 4 times of computational budget (×T4). In particular for the two contour algorithms, we set ζ = 20 and choose a constant step size for the stochastic approximation to 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 et al. (2018), the inclusion of advanced techniques may change the optimal settings of the hyperparameters. Nevertheless, we try to 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, we are given 20 random mushrooms and train the model 16 iterations every step for the parallel algorithms (×P4); we train pCSGLD×T4 64 iterations every step. EpsGreedy decays the learning rate by a 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 set the energy depth ∆u as 10. We fix the hyperrameter ζ = 20. The step sizes for pICSGLD×P4 and pCSGLD×T4 are chosen as 0.03 and 0.006, respectively. A proper regularizer is adopted for the low importance weights. See Table 2 for details.

heading: D.4 UNCERTAINTY ESTIMATION
text: All the algorithms, excluding M-SGD×P4, choose a temperature of 0.0003 †. We run the parallel algorithms 500 epochs (×P4) and run the single-chain algorithms 2000 epochs (×T4). The initial learning rate is 2e-6 (Bayesian settings), which corresponds to the standard 0.1 for averaged data likelihood.
We train cycSGHMC×T4 and MultiSWAG×T4 based on the cosine learning rates with 10 cycles. The learning rate in the last 15% of each cycle is fixed at a constant value. MultiSWAG simulates 10 random models at the end of each cycle. M-SGD×P4 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 56 and a correction of 1000 for WRN-16-8. The learning rate is annealed at 250 and 375 epochs with a factor of 0.2. ICSGHMC×P4 also applies the same learning rate. We choose m = 200 and ∆u = 200 for ResNet20, 32, and 56 and ∆u = 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 ICSGHMC×P4 because they are the only two algorithms that require accurate estimations of the energy. We only update control variates every 2 epochs in the last 100 epochs, which maintain a reasonable training time and a higher reduction of variance due to a small learning rate. Other algorithms yield a worse performance when variance reduction is applied to the gradients.

heading: D.5 EMPIRICAL VALIDATION OF REDUCED VARIANCE
text: To compare the θ’s learned from ICSGLD and CSGLD, we try to simulate from a Gaussian mixture distribution 0.4N(−6, 1) + 0.6N(4, 1), where N(u, v) denotes a Gaussian distribution with mean u and standard deviation v. We fix ζ = 0.9 and ∆u = 1. We run ICSGLD with 1,000,000 iterations
based on 10 interacting parallel chains and run CSGLD with 10,000,000 iterations using a single chain. We refer to them as ICSGLD×P10 and CSGLD×T10, respectively. The rest of the settings follows from the experimental setup in section 4.1 (Deng et al., 2020a).
To measure the variance of the estimates, we repeated the experiments 10 times and present the mean and two standard deviations for both CSGLD×T10 and ICSGLD×P10 in Figure 7. The results indicate that both estimates of θζ (by CSGLD and ICSGLD) converge to the equilibrium that approximates the ground truth of the density of states. Notably, ICSGLD×P10 yields a significantly smaller variance than CSGLD×T10, but with the same computational budget. This shows the clear advantage of ICSGLD (many interacting short runs) over CSGLD (a single long run) in tackling the large variance issue for importance sampling.
†We use various data augmentation techniques, such as random flipping, cropping, and random erasing (Zhong et al., 2017). This leads to a much more concentrated posterior and requires a very low temperature.

Review summaries:
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.

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

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.

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.