STEIN SELF-REPULSIVE DYNAMICS: BENEFITS FROM PAST SAMPLES
Year: 2019

Abstract:
We propose a new Stein self-repulsive dynamics for obtaining diversified samples from intractable un-normalized distributions. Our idea is to introduce Stein variational gradient as a repulsive force to push the samples of Langevin dynamics away from the past trajectories. This simple idea allows us to significantly decrease the auto-correlation in Langevin dynamics and hence increase the effective sample size. Importantly, as we establish in our theoretical analysis, the asymptotic stationary distribution remains correct even with the addition of the repulsive force, thanks to the special properties of the Stein variational gradient. We perform extensive empirical studies of our new algorithm, showing that our method yields much higher sample efficiency and better uncertainty estimation than vanilla Langevin dynamics.

text: We propose a new Stein self-repulsive dynamics for obtaining diversified samples from intractable un-normalized distributions. Our idea is to introduce Stein variational gradient as a repulsive force to push the samples of Langevin dynamics away from the past trajectories. This simple idea allows us to significantly decrease the auto-correlation in Langevin dynamics and hence increase the effective sample size. Importantly, as we establish in our theoretical analysis, the asymptotic stationary distribution remains correct even with the addition of the repulsive force, thanks to the special properties of the Stein variational gradient. We perform extensive empirical studies of our new algorithm, showing that our method yields much higher sample efficiency and better uncertainty estimation than vanilla Langevin dynamics.

heading: 1 INTRODUCTION
text: Drawing samples from complex un-normalized distributions is one of the most basic problems in statistics and machine learning, with broad applications to enormous research fields that rely on probabilistic modeling. Over the past decades, large amounts of methods have been proposed for approximate sampling, including both Markov Chain Monte Carlo (MCMC) (e.g., Brooks et al., 2011) and variational inference (e.g., Wainwright et al., 2008).
MCMC works by simulating Markov chains whose stationary distributions match the distributions of interest. Despite nice asymptotic theoretical properties, MCMC is widely criticized for its slow convergence rate in practice. In difficult problems, the samples drawn from MCMC are often found to have high auto-correlation across time, meaning that the Markov chains explore very slowly in the configuration space. When this happens, the samples returned by MCMC only approximate a small local region, and under-estimate the probability of the regions un-explored by the chain.
Stein variational gradient descent (SVGD) (Liu & Wang, 2016) is another type of approximation sampling methods designed to overcome the limitation of MCMC. Instead of drawing random samples sequentially, SVGD evolves a pre-defined number of particles (or sample points) in parallel with a special interacting particle system to match the distribution of interest by minimizing the KL divergence. In SVGD, the particles interact with each other to simultaneously move towards the high probability regions following the gradient direction, and also move away from each other due to a special repulsive force. As a result, SVGD allows us to obtain diversified samples that correctly represent the variation of the distribution of interest. SVGD has been found a promising tool for solving difficult sampling problems in which diversity promotion is critical (e.g., Feng et al., 2017; Haarnoja et al., 2017; Pu et al., 2017; Liu et al., 2017; Gong et al., 2019). Various extensions have been developed (e.g., Han & Liu, 2018; Chen et al., 2018; Liu et al., 2019; Wang et al., 2019a).
However, one problem of SVGD is that it theoretically requires to run an infinite number of chains in parallel in order to approximate the target distribution asymptotically (Liu, 2017). With a finite number of particles, the fixed point of SVGD does still provide a prioritized, partial approximation to the distribution in terms of the expectation of a special case of functions (Liu & Wang, 2018). Nevertheless, it is still desirable to develop a variant of “single-chain SVGD”, which only requires to run a single chain sequentially like MCMC to achieve the correct stationary distribution asymptotically in time, with no need to take the limit of infinite number of parallel particles.
In this work, we propose an example of single-chain SVGD by integrating the special repulsive mechanism of SVGD with gradient-based MCMC such as Langevin dynamics. Our idea is to use repulsive term of SVGD to enforce the samples in MCMC away from the past samples visited at previous iterations. Such a new self-repulsive dynamics allows us to decrease the auto-correlation in MCMC and hence increase the mixing rate, but still obtain the same stationary distribution thanks to the special property of the SVGD repulsive mechanism.
We provide throughout theoretically analysis of our new method, establishing it asymptotic convergence to the target distribution. As we show in the work, the analysis is highly non-trivial, because our new self-repulsive dynamics is a non-linear, high-order Markov process. Empirically, we extensively evaluate our methods on an array of challenging sampling tasks, showing that our method yields much better uncertainty estimation and larger effective sample size.

heading: 2 BACKGROUND: LANGEVIN DYNAMICS AND SVGD
text: In this section, we give a brief introduction on Langevin dynamics (Rossky et al., 1978) and Stein Variational Gradient Descent (SVGD) (Liu & Wang, 2016), which we integrate teogehter to develop our new self-repulsive dynamics for more efficient sampling.
Langevin dynamics Langevin dynamics is a basic gradient based MCMC method. For some target distribution on Rd with density function ρ∗(θ) ∝ exp(−V (θ)), where V : Rd 7→ R is the potential function, the (Euler-discrerized) Langevin dynamics simulates a Markov chain with the following rule:
θk+1 = θk − η∇V (θk) + √
2ηek, ek ∼ N (0, I), where k denotes the number of iterations, {ek} are independent standard Gaussian noise, and η is a step size parameter. It is well known that the limiting distribution of θk when k →∞ approximates the target distribution when η is sufficiently small.
Because the updates in Langevin dynamics are local and incremental, new points generated by the dynamics is highly correlated to the past sample. As a result, we need to run Langevin dynamics sufficiently long in order to obtain diverse samples.
Stein Variatinal Gradient Descent (SVGD) Different from Langevin dynamics, SVGD iteratively evolves a pre-defined number of particles in parallel. Starting from an initial set of particles {θi0 : i = 1, ...,M}, SVGD updates the M particles in parallel by
θik+1 = θ i k + ηg(θ i k; δ̂ M k ) ∀i = 1, . . . ,M,
where the velocity field, which we denote by g(θik; δ̂ M k ), depends the empirical distribution of the current set of particles δ̂Mk := 1 M ∑M j=1 δθjk in the following way,
g(θik; δ̂ M k ) = Eθ∼δ̂Mk −K(θ,θik)∇V (θ)︸ ︷︷ ︸ Confining Term + ∇θK(θ,θik)︸ ︷︷ ︸ Repulsive Term  . Here δθ is the Dirac measure centered at θ, and hence Eθ∼δ̂Mk [·] denotes averaging on the particles. The K(·, ·) is a positive definite kernel, such as the Gaussian RBF kernel, specified by users.
Note that g(θik; δ̂ M k ) consists of a confining term and repulsive term: the confining term pushes particles to move towards high density region, and the repulsive term prevents the particles from colliding with each other. It is the balance of these two terms that allows us to asymptotically approximate the target distribution ρ∗(θ) ∝ exp(−V (θ)) at the fixed point, when the number of particles goes to infinite. We refer the readers to Liu & Wang (2016); Liu (2017); Liu & Wang (2018) for throughout theoretical justifications of SVGD. But a quick, informal way to justify the SVGD update is through the Stein’s identity, which shows that the velocity field g(θ; ρ) equals zero when ρ equals the true distribution ρ∗, that is,
g(θ′; ρ∗) = Eθ∼ρ∗ [−K(θ,θ′)∇V (θ) +∇θK(θ,θ′)] = 0, ∀θ′ ∈ Rd. (1)
This shows that SVGD would converge if the particle distribution already forms a closed approximation to the target distribution p, meaning that the target distributions forms an (approximate) fixed point of the update.

heading: 3 STEIN SELF-REPULSIVE DYNAMICS
text: In this work, we propose to integrate Langevin dynamics and SVGD to simultaneously decrease the auto-correlation of Langevin dynamics and eliminate the need for running parallel chains in SVGD. The idea is to use Stein repulsive force between the the current sample and the past samples, hence forming a new self-avoiding dynamics with fast convergence speed.
Specifically, assume we run a single Markov chain like Langevin dynamics, where θk denotes the particle at the k-th iteration. Denote by δ̃Mk the empirical distribution of M samples taken from the past iterations, i.e.,
δ̃Mk := 1
M M∑ j=1 δθk−jcη , cη = c/η,
where cη is a thinning factor, which scales inversely with the step size η, introduced to slim the sequence of past samples. Compared with the δ̂Mk in SVGD, which is averaged over M parallel particles, δ̃Mk is averaged across time over M past samples. Given this, our Stein self-repulsive dynamics updates the sample via
θk+1 ← θk + (−ηV (θk) + √
2ηek)︸ ︷︷ ︸ Langevin + ηαg(θk; δ̃ M k )︸ ︷︷ ︸ Stein Repulsive , (2)
in which the particle is updated with the typical Langevin gradient, plus a Stein repulsive force against the samples from the previous iterations. α ≥ 0 is a parameter that controls the magnitude of the Stein repulsive term. In this way, the particles are pushed away from the past samples, and hence admits lower auto-correlation and faster convergence speed. Importantly, the addition of the repulsive force does not impact the asymptotic stationary distribution, thanks to Stein’s identity in (1). This is because when if the self-repulsive dynamics have converged to the target distribution ρ∗, such that θk ∼ ρ∗ for all k, the Stein self-repulsive term would equal to zero in expectation due to Stein’s identity and hence does not introduce additional bias over Langevin dynamics. Rigorous theoretical analysis of this idea is developed in Section 4.
Practical Algorithm Because δ̃Mk is averaged across the past samples, it is necessary to introduce a burn-in phase with the repulsive dynamics. Therefore, our overall procedure works as follows,
θk+1 =
{ θk − η∇V (θk) + √ 2ηek, k < Mcη
θk + η [ −∇V (θk) + αg(θk; δ̃Mk ) ] + √ 2ηek, k ≥Mcη. (3)
It includes two phases. The first phase is the same as the Langevin dynamics which collects the initial M samples used in the second phase while serves as a warm start. The repulsive gradient update is introduced in the second phase to encourage the dynamics to visit the under-explored density region. We call this particular instance of our algorithm Self-Repulsive Langevin dynamics (SRLD), self-repulsive variants of more general dynamics is discussed in Section 5.
Remark Notice that, the first phase is introduced to collect the initial M samples. However, it’s not really necessary to generate the initial M samples with Langevin dynamics. We can simply use some other initialization distribution and get M initial samples from that distribution. In practice, we find using Langevin dynamics to collect the initial samples is natural and it can also be viewed as the burn-in phase before sampling, so we use (3) in all of the other experiments.
Remark The general idea of introducing self-repulsive terms inside MCMC or other iterative algorithms is not new itself. For example, in molecular dynamics simulations, an algorithm called metadynamics (Laio & Parrinello, 2002) has been widely used, in which the particles are repelled away from the past samples in a way similar to our method, but with a typical repulsive function, such as ∑ j D(θk, θk−jcη ), where D(·, ·) is any notation of dis-similarity. However, introducing
an arbitrary repulsive force would alter the stationary distribution of the dynamics, introducing a harmful bias into the algorithm. The key highlight of our approach, as reflected by our theoretical results in Section 4, is the unique property of the Stein repulsive term, that allows us to obtain the correct stationary distribution even with the addition of the repulsive force.
Remark Recently, (Zhang et al., 2018) proposed a different combination of SVGD and Langevin dynamics, in which the Langevin force is directly added to a set of particles that evolve in parallel with SVGD. Using our terminology, their system is
θik+1 = θ i k + (−ηV (θik) + √ 2ηeik) + ηαg(θ i k; δ̂ M k ), ek ∼ N (0, I) ∀i = 1, . . . ,M.
This is significantly different from our method on both motivation and practical algorithm. Their algorithm still requires to simulate M chains of particles in parallel like SVGD, and was proposed to obtain easier theoretical analysis than the deterministic dynamics of SVGD. Our work is instead motivated by the practical need of decreasing the auto-correlation in Langevin dynamics, and avoiding the need of running multiple chains in SVGD, and hence must be based on self-repulsion against past samples along a single chain.
In another recent work, (Chen et al., 2018) proposed a π-SGLD method, which simulates the linear combination of the evalutionary partial differential equations of SVGD and Langevin dynamics using discrete gradient flow with blob-based method. Their method is again motivated by the theoretical interest of discovering new categories of algorithms, and does not involve self-repulsive on a single chain like our method.
An Illustrative Example Here we give an illustrative example to show the key advantage of our self-repulsive dynamics. Assume that we want to sample from a bi-variate Normal distribution shown in Figure 1. Unlike standard settings, we assume that we have already obtained some initial samples (yellow dots in Figure 1) before running the dynamics. The initial samples are assumed to concentrate on the left part of the target distribution as shown in Figure 1. In this extreme case, since the left part of the distribution is over-explored by the initial samples, it is desirable to have the subsequent new samples to concentrate more on the un-explored part of the distribution. However, standard Langevin dynamics does not take this into account, and hence yielding a bias overall representation of the true distribution (see the left panel). With our self-repulsive dynamics, the new
samples are forced to explore the un-explored region more frequently, allowing us to obtain a much more accurate approximation when combining the new and initial samples.

heading: 4 THEORETICAL ANALYSIS OF STEIN SELF-REPULSIVE DYNAMICS
text: In this section, we provide theoretical analysis of the self-repulsive dynamics. We establish that our self-repulsive dynamics converges to the correct target distribution asymptotically, in the limit when M approaches to infinite and the step size η approaches to 0. This is a highly non-trivial task, because the self-repulsive dynamics is a highly complex, non-linear and high order Markov stochastic process. We attack this problem by breaking the proof into the following three steps illustrated in Figure 2:
1) At the limit of large particle sizes M → ∞ (called the mean field limit), we show that practical dynamics in (3) is closely approximated by a discrete-time mean-field dynamics characterized by (4) below.
2) By further taking the limit of small step size η → 0+ (called the continuous-time limit), the dynamics in (4) converges to a continuous-time mean-field dynamics characterized by (5) below.
3) We show that the dynamics in (5) converges to the target distribution.
Remark As we mentioned in Section 3, we introduce the first phase to collect the initial M samples for the second phase, and our theoretical analysis follows this setting to make our theory as close to the practice. However, the theoretical analysis can be generalized to the setting of drawing M initial samples from some initialization distribution with almost identical argument.
Notations We use ‖·‖ and 〈·, ·〉 to represent the `2 vector norm and inner product respectively. The Lipschitz norm and bounded Lipschitz norm of a function f are defined by ‖f‖Lip and ‖f‖BL. The KL divergence, Wasserstein-2 distance and Bounded Lipschitz distance between distribution ρ1, ρ2 are denoted as DKL[ρ1 ‖ ρ2], W2[ρ1, ρ2] and DBL[ρ1, ρ2] respectively.

heading: 4.1 MEAN-FIELD AND CONTINUOUS-TIME LIMITS
text: To fix the notation, we denote by ρk := Law(θk) the distribution of θk at time k of the practical self-repulsive dynamics (3), which we refer as the practical dynamics in the sequel, when the initial particle θ0 is drawn from an initial continuous distribution ρ0 supported on Rd. Note that given ρ0, the subsequent ρk can be recursively defined through dynamics (3). Due to the diffusion noise in Langevin dynamics, all ρk are continuous distributions supported on Rd. We now introduce the limit dynamics when we take the mean-field limit (M → +∞) and then the continuous-time limit (η → 0+).
Discrete-Time Mean-Field Dynamics (M → +∞) In the limit of M → ∞, we show that our practical dynamics is closely approximated by the following dynamics, in which the empirical measures are replaced by the continuous distributions of the particles themselves.
θ̃k+1 =
{ θ̃k − η∇V (θ̃k) + √ 2ηek, k ≤Mcη
θ̃k + η [ −∇V (θ̃k) + αg(θk, ρ̃Mk ) ] + √ 2ηek, k ≥Mcη (4)
where ρ̃Mk = 1 M ∑M j=1 ρ̃k−jcη by ρ̃k := Law(θ̃k) the (smooth) distribution of θ̃k at time-step k when the dynamics is initialized with θ̃0 ∼ ρ̃0 = ρ0. Compared with the practical dynamics
in (3), the difference is that the empirical distribution δ̃Mk is replaced by the smooth distribution ρ̃Mk . Similar to the recursive definition of ρk following dynamics (3), ρ̃k is also recursively defined through dynamics (4), starting from ρ̃0 = ρ0. As we show in Theorem 4.3, if the auto-correlation of θk decays fast enough and M is sufficiently large, ρ̃Mk is well approximated by the empirical distribution δ̃Mk in (3), and further the two dynamics ((3) and (4)) converges to each other in the sense that W2[ρk, ρ̃k] → 0 as M → ∞ for any k. Note that in taking the limit of M , we need to ensure that we run the dynamics for more than Mcη steps. Otherwise, SRLD degenerates to Langevin dynamics as we stop the chain before we collect M samples.
Continuous-Time Mean-Field Dynamics (η → 0+) In the limit of zero step size (η → 0+), the discrete-time mean field dynamics in (4) can be shown to converge to the following continuous-time mean-field dynamics:
dθ̄t = { −∇V (θ̄t)dt+ dBt, t ∈ [0,Mc)[ −∇V (θ̄t) + αg(θk, ρ̄Mt ) ] dt+ dBt, t ≥Mc,
(5)
where ρ̄Mt := 1 M ∑M j=1 ρ̄t−jc(·) and ρ̄t = Law ( θ̄t )
is the distribution of θ̄t at a continuous time point t with θ0 initialized by θ̄0 ∼ ρ̃0 = ρ0. We prove that (5) is closely approximated by (4) with small step size in the sense that DKL[ρ̃k ‖ ρ̄kη]→ 0 as η → 0 in Theorem 4.2, and importantly, the stationary distribution of (5) equals to the target distribution ρ∗(θ) ∝ exp(−V (θ)).

heading: 4.2 ASSUMPTIONS
text: We first introduce the techinical assumptions used in our theoretical analysis.
Assumption 4.1. (RBF Kernel)
We use RBF kernel i.e. K(θ1,θ2) = exp(−‖θ1 − θ2‖2 /σ) for some fixed 0 < σ <∞.
We only assume the RBF kernel for the simplicity of our analysis. However, it is straightforward to generalize our theory to other positive definite kernels.
Assumption 4.2. (V is dissipative and smooth)
Assume that 〈θ,−∇V (θ)〉 ≤ b1 − a1 ‖θ‖2 and ‖∇V (θ1)−∇V (θ2)‖ ≤ b1 ‖θ1 − θ2‖. We also assume that ‖∇V (0)‖ ≤ b1. Here a1 and b1 are some finite positive constant. Assumption 4.3. (Regularity Condition)
Assume Eθ∼ρ0 [‖θ‖ 2 ] < 0. Define ρMk = ∑M j=1 ρk−jcη/M , assume there exists B <∞ such that
inf k≥Mcη ( sup ‖θ‖≤B E ∥∥∥g(θ; δ̃Mk )− g(θ; ρMk )∥∥∥2 ) > 0.
Assumption 4.4. (Strong-convexity)
Suppose that 〈∇V (θ1)−∇V (θ2),θ1 − θ2〉 ≥ L ‖θ1 − θ2‖2 for a positive constant L.
Remark Assumption 4.2 is standard in the existing Langevin dynamics analysis (see Dalalyan (2017); Raginsky et al. (2017) for example). Assumption 4.3 is a weak condition as it assumes that the dynamics can not degenerate into one local mode and stop moving anymore. This assumption is expected to be true when we have diffusion terms like the Gaussian noises in our self-repulsive dynamics. Assumption 4.4 is a classical assumption on the existing Langevin dynamics analysis with convex potential Dalalyan (2017); Durmus et al. (2019). Although being a bit strong, this assumption broadly applies to posterior inference problem in the limit of big data, as the posterior distribution converges to Gaussian distributions for large training set as shown by Bernstein-von Mises theorem. It is technically possible to further generalize our results to the non-convex settings with a refined analysis, which we leave as future work. This work focuses on the classic convex setting for simplicity.

heading: 4.3 MAIN THEOREMS
text: All of the proof in this section can be found in Appendix B.5.
We first prove that the limiting distribution of the continuous-time mean field dynamics (5) is the target distribution. This is achieved by writing dynamics (5) into the following (non-linear) partial differential equation:
∂tρ̄t = { ∇ · (−∇V ρ̄t) + ∆ρ̄t t ∈ [0,Mc) ∇ · [( −∇V + αg(·, ρ̄Mt ) ) ρ̄t ] + ∆ρ̄t, t ≥Mc.
Theorem 4.1. (Stationary Distribution)
Given some finiteM , c and α, and suppose that the limiting distribution of dynamics (5) exists. Then the limit distribution is unique, and equals to ρ∗(θ) ∝ exp(−V (θ)).
We then give the upper bound on the discretization error, which can be characterized by analyzing the KL divergence between ρ̃k and ρ̄kη . Theorem 4.2. (Time Discretization Error)
Given some sufficiently small step size η and choose α < a2/(2b1 + 4/σ). Under assumption 4.1, 4.2, 4.3 and cη = c/η. we have for some constant C,
max l∈{0,...,k}
DKL [ρ̄lη ‖ ρ̃l] ≤
{ O ( η + kη2 ) k ≤Mcη − 1
O ( η +Mcη + α2MceCα 2(kη−Mc)η2 ) k ≥Mcη.
With this theorem, we can know that if η is small enough, then the discretization error is small and ρ̃ approximates ρ̄ closely. Next we give result on the mean field limit of M →∞. Theorem 4.3. (Mean-Field Limit)
Under assumption 4.1, 4.2, 4.3, and 4.4, suppose that we choose α and η such that−(a1−2αb1/σ)+ ηb1 < 0; 2αησ (b1 + 1) < 1; a2 − α ( 2b1 + 4 σ ) > 0; Then there exists a constant c2, such that when L/a ≥ c2 and we have
W22[ρk, ρ̃k] = O ( α2/M + η2 ) for k ≥Mcη,
and W22[ρk, ρ̃k] = 0 if k ≤Mcη − 1.
Remark Combine all the result, we have limk,M→∞,η→0+ DBL [ρk, ρ∗] = 0, where the joint limit requires kη → ∞, exp(Cα2kη)η2 = o(1) and kηMc = γ (1 + o(1)), with γ > 1. Notice that if γ ≤ 1, the dynamics is degenerated to Langevin dynamics.

heading: 5 EXTENSION TO GENERAL DYNAMICS
text: Although we have focused on self-repulsive Langevin dynamics, our Stein self-repulsive idea can be broadly combined with general gradient-based MCMC. Following Ma et al. (2015), we consider the following general class of sampling dynamics for drawing samples from p(θ) ∝ exp(−V (θ)):
dθt = −f(θ)dt+ √ 2D(θ)dBt,
with f(θ) = [D(θ) +Q(θ)]∇V (θ)− Γ(θ), Γi(θ) = d∑ j=1 ∂ ∂θj (Dij(θ) +Qij(θ)) .
where D is a positive semi-definite diffusion matrix that determines the strength of the Brownian motion and Q is a skew-symmetric curl matrix that can represent the traverse effect (e.g. in Neal et al., 2011; Ding et al., 2014). By adding the Stein repulsive force, we obtain the following general self-repulsive dynamics
dθ̄t =
{ −f(θ)dt+ √ 2D(θ)dBt, t ∈ [0,Mc)
− ( f(θ) + αg(θ̄t; ρ̄ M t ) ) dt+ dBt, t ≥Mc
(6)
where ρ̄t := Law(θ̄t) is again the distribution of θ̄t following (6) when initalized at θ̄0 ∼ ρ̄0. Similar to the case of Langevin dynamics, this process also converges to the correct target distribution, and can be simulated by practical dynamics similar to (3). Theorem 5.1. (Stationary Distribution)
Given some finiteM , c and α, and suppose that the limiting distribution of dynamics (6) exists. Then the limiting distribution is unique and equals the target distribution ρ∗(θ) ∝ exp(−V (θ)).

heading: 6 EXPERIMENTS
text: In this section, we evaluate the proposed method in various challenging problems, including sampling the posteriors of Bayesian Neural Networks, and uncertainty estimation in Reinforcement Learning. Our results show that our Self-Repulsive Langevin dynamics (SRLD) yields much higher sample efficiency than vanilla Langevin dynamics. Our code is available along the submission.

heading: 6.1 SYNTHETIC EXPERIMENT
text: A Correlated 2D Distribution This experiment aims to show how the repulsive gradient helps explore the whole distribution. Following Ma et al. (2015), we compare the sampling efficiency on the following correlated 2D distribution with density
ρ∗([θ1, θ2]) ∝ −θ41/10− ( 4 (θ2 + 1.2)− θ21 )2 /2.
We compare the SRLD with vanilla Langevin dynamics, and evaluate the sample quality by Maximum Mean Discrepancy (MMD) (Gretton et al., 2012), Wasserstein-1 Distance and effective sample size (ESS). Notice that the finite sample quality of gradient based MCMC method is highly related to the step size. Compared with Langevin dynamics, we have an extra repulsive gradient and thus we implicitly have larger step size. To rule out this effect, we set different step sizes of the two dynamics so that the gradient of the two dynamics has the same magnitude. In addition, to decrease the influence of random noise, the two dynamics are set to have the same initialization and use the same sequence of Gaussian noise. We collect the sample of every iteration. We repeat the experiment 20 times with different initialization and sequence of Gaussian noise.
Figure 3 summarizes the result with different metrics. We can see that SRLD have a significantly smaller MMD and Wasserstein-1 Distance as well as a larger ESS compared with the vanilla Langevin dynamics. Moreover, the introduced repulsive gradient creates a negative auto-correlation between samples. Figure 4 shows a typical trajectory of the two sampling dynamics. We can see that SRLD have a faster mixing rate than vanilla Langevin dynamics. Note that since we use the same sequence of Gaussian noise for both algorithms, the difference is mainly due to the use of repulsive gradient rather than the randomness.
Mixture of Gaussian Distribution We aim to show how the repulsive gradient helps the particle escape from the local high density region by sampling the 2D mixture of Gaussian distribution using SRLD and Langevin dynamics. The target density is set to be
ρ∗(θ) ∝ 1 2
exp ( −‖θ − 1‖2 /2 ) + 1
2 exp
( −‖θ + 1‖2 /2 ) , θ = [θ1, θ2] >,
where 1 = [1, 1]>. This target distribution have two mode at −1 and 1, and vanilla Langevin dynamics can stuck in one mode while keeps the another mode under-explored (as the gradient of energy function can dominate the update of samples). We use the same evaluation method, step sizes, initialization and Gaussian noise as the previous experiment. We collect one sample every 100 iterations and the experiment is repeated for 20 times. Figure 5 shows that SRLD again consistently outperforms the Langevin dynamics on all of the evaluation metrics.

heading: 6.2 BAYESIAN NEURAL NETWORK
text: Bayesian Neural Network is one of the most important methods in Bayesian Deep Learning with wide application in practice. Here we test the performance of SRLD on sampling the posterior of Bayesian Neural Network on the UCI datasets (Dua & Graff, 2017). We assume the output is normal distributed, with a two-layer neural network with 50 hidden units and tanh activation to predict the mean of outputs. We set a Γ(1, 0.1) prior for the inverse output variance. All of the datasets are randomly partitioned into 90% for training and 10% for testing. The results are averaged over 20 random trials. We set the mini-batch size to be 100 and the number of past samples M to be 10. In all experiments, we use RBF kernel with bandwidth set by the median trick as suggested in Liu & Wang (2016). For SVGD, we use the original implementation with 20 particles by Liu & Wang (2016). We run 50000 iterations for each methods, and for LD and SRLD, the first 40000 iteration is discarded as burn-in. We use a thinning factor of cη = c/η = 100 and in total we collect 100 samples from the posterior distribution. For each dataset, we generate 3 extra data splits for tuning the step size for each method.
Table 1 shows the average test RMSE and test log-likelihood and their standard deviation. The method that has the best average performance is marked as boldface. We observe that a large portion of the variance is due to the random partition of the dataset. Therefore, to show the statistical significance, we use the matched pair t-test to test the statistical significance, mark the methods that perform the best with a significance level of 0.05 with underlines. Note that the results of SRLD/LD and SVGD is not very comparable, because SRLD/LD are single chain methods which averages across time, and SVGD is a multi-chain method that only use the results of the last iteration. We
provide additional results in appendix that SRLD averaged on 20 particles (across time) can also achieve similar or better results as SVGD with 20 (parallel) particles.

heading: 6.3 CONTEXTUAL BANDITS
text: We evaluate the quality of uncertainty estimation provided by our methods on several contextual bandits problems. Uncertainty estimation is a key component of contextual bandits. If the agent makes decisions with a poorly estimated uncertainty, the decisions will finally lead to catastrophic failure through the feedback loops (Riquelme et al., 2018).
Though in principle all of the MCMC methods return the samples follow the true posterior if we can run infinite MCMC steps, in practice we can only obtain finite samples as we only have finite time to run the MCMC sampler. In this case, the auto-correlation issue can lead to the under-estimate the uncertainty, which will cause the failure on all of the reinforcement learning problems that need exploration.
We consider the posterior sampling (a.k.a Thompson sampling) algorithm with Bayesian neural network as the function approximator. We follows the experimental setting from Riquelme et al. (2018). The only difference is that we change the optimization of the objective (e.g. evidence lower bound (ELBO) in variational inference methods) into running MCMC samplers. We set the step of samplers equal to the optimization step, and use a thinning factor of 100. We compare the SRLD with the Langevin dynamics on the Mushroom and Wheel bandits from (Riquelme et al., 2018), and include SVGD as a baseline. For more detailed introduction and setup, see Appendix A.5.
The cumulative regret is shown in Table 2. SVGD is known to have the under-estimated uncertainty for Bayesian neural network if particle number is limited (Wang et al., 2019b), and as a result, has the worst performance among the three methods. SRLD is slightly better than vanilla Langevin dynamics on the simple Mushroom bandits. On the much more harder Wheel bandits, SRLD is significantly better than the vanilla Langevin dynamics, which shows the improving uncertainty estimation of our methods within finite number of samples.

heading: 7 CONCLUSION
text: We propose a Stein self-repulsive dynamics which applies Stein variational gradient to push samples from MCMC dynamics away from its past trajectories. This allows us to significantly decrease the auto-correlation of MCMC, increasing the sample efficiency for better estimation. The advantages of our method are extensive studied both theoretical and empirical analysis in our work. In future work, we plan to investigate the combination of our Stein self-repulsive idea with more general MCMC procedures, and explore broader applications.

heading: A EXPERIMENT DETAILS AND ADDITIONAL EXPERIMENT RESULT

heading: A.1 SYNTHETIC 2D MIXTURE OF GAUSSIAN EXPERIMENT
text: To provide more evidence on the effectiveness of SRLD on escaping from local high density region, we plot the sampling trajectory of SRLD and vanilla Langevin dynamics on the mixture of Gaussian mentioned in Section 6.1. We can find that, when both of the methods obtain 200 samples, SRLD have started to explore the second mode, while vanilla Langevin dynamics still stuck in the original mode. When both of the methods have 250 examples, the vanilla Langevin dynamics just start to explore the second mode, while our SRLD have already obtained several samples from the second mode, which shows our methods effectiveness on escaping the local mode.

heading: A.2 SYNTHETIC HIGHER DIMENSIONAL GAUSSIAN EXPERIMENT
text: To show the performance of SRLD in higher dimensional case with different value of α, we additionally considering the problem on sampling from Gaussian distribution with d = 100 and covariance Σ = 0.5I. We run SRLD with α = 100, 50, 20, 10, 0 and the case α = 0 reduces to Langevin. We collect 1 sample every 10 iterations. The other experiment setting is the same as the toy examples in the main text. The results are summarized at Figure [7]. In this experiment, we set one SRLD with an inappropriate α = 100. For this chain, the repulsive gradient gives strong repulsive force and thus has the largest ESS and the fastest decay of autocorrelation. While the inappropriate value α induces too much extra approximation error and thus its performance is not as good as these with smaller α (see MMD and Wasserstein distance). This phenomenon matches our theoretical finding.

heading: A.3 SYNTHETIC HIGHER DIMENSIONAL MIXTURE OF GAUSSIAN EXPERIMENT
text: We also consider sampling from the mixture of Gaussian with d = 20. The target density is set to be
ρ∗(θ) ∝ 1 2 exp
( −0.5 ∥∥∥θ −√2/d1∥∥∥2)+ 1 2 exp ( −0.5 ∥∥∥θ +√2/d1∥∥∥2) , where θ = [θ1, ..., θ20]> and 1 = [1, ..., 1]>. And thus the mean of the two mixture component is with distance 2 √ 2. We run SRLD with α = 20, 10, 5, 0 (and when α = 0, it reduces to LD). The other experiment setting is the same as the low dimensional mixture Gaussian case. Figure [8] summarizes the result. As shown in the figure, when α becomes larger, the repulsive forces helps the sampler better explore the density region.

heading: A.4 UCI DATASETS
text: We show some additional experiment result on posterior inference on UCI datasets. As mentioned in Section 6.2, the comparison between SVGD and SRLD is not direct as SVGD is a multiple-chain method with fewer particles and SRLD is a single chain method with more samples. To show more detailed comparison, we compare the SVGD with SRLD using the first 20, 40, 60, 80 and 100 samples, denoted as SRLD-n where n is the number of samples used. Table 3 shows the result of averaged test RMSE and table 4 shows the result of averaged test loglikelihood. For SRLD with different number of samples, the value is set to be boldface if it has better average performance than SVGD. If it is statistical significant with significant level 0.05 using a matched pair t-test, we add an underline on it.
Figure 9 and 10 give some visualized result on the comparison with Langevin dynamics and SRLD. To rule out the variance of different splitting on the dataset, the errorbar is calculated based on the difference between RMSE of SRLD and RMSE of Langevin dynamcis in 20 repeats (And similarily for test log-likelihood). And we only applied the error bar on Langevin dynamics.

heading: A.5 CONTEXTUAL BANDIT
text: Contextual bandit is a class of online learning problems that can be viewed as a simple reinforcement learning problem without transition. For a completely understanding of contextual bandit problems, we refer the readers to the Chapter 4 of (Bubeck et al., 2012). Here we include the main idea for completeness. In contextual bandit problems, the agent needs to find out the best action given some observed context (a.k.a the optimal policy in reinforcement learning). Formally, we define S as the context set and K as the number of action. Then we can concretely describe the contextual bandit problems as follows: for each time-step t = 1, 2, · · · , N , where N is some pre-defined time horizon (and can be given to the agent), the environment provides a context st ∈ S to the agent, then the agent should choose one action at ∈ {1, 2, · · · ,K} based on context st. The environment will return a (stochastic) reward r(st, at) to the agent based on the context st and the action at that similar to the reinforcement learning setting. And notice that, the agent can adjust the strategy at each time-step, so that this kinds of problems are called “online” learning problem.
Solving the contextual bandit problems is equivalent to find some algorithms that can minimize the pseudo-regret (Bubeck et al., 2012), which is defined as:
R S N = max
π:S→{1,2,··· ,K} E [ N∑ t=1 r(st, g(st))− N∑ t=1 r(st, at) ] . (7)
where π denotes the deterministic mapping from the context set S to actions {1, 2, · · · ,K} (readers can view π as a deterministic policy in reinforcement learning). Intuitively, this pseudo-regret measures the difference of cumulative reward between the action sequence at and the best action sequence π(st). Thus, an algorithm that can minimize the pseudo-regret (7) can also find the best π.
Posterior sampling (a.k.a. Thompson sampling; Thompson, 1933) is one of the classical yet successful algorithms that can achieve the state-of-the-art performance in practice (Chapelle & Li, 2011). It works by first placing an user-specified prior µ0s,a on the reward r(s, a), and each turn make decision based on the posterior distribution and update it, i.e. update the posterior distribution µts,a with the observation r(st−1, at−1) at time t − 1 where at−1 is selected with the posterior distribution: each time, the action is selected with the following way:
at = arg max a∈{1,2,··· ,K}
r̂(st, a), r̂(st, a) ∼ µts,a.
i.e., greedy select the action based on the sampled reward from the posterior, thus called “Posterior Sampling”. Algorithm 1 summarizes the whole procedure of Posterior Sampling.
Algorithm 1 Posterior sampling for contextual bandits
Input: Prior distribution µ0s,a, time horizon N for time t = 1, 2, · · · , N do
observe a new context st ∈ S, sample the reward of each action r̂(st, a) ∼ µts,a, a ∈ {1, 2, · · · ,K}, select action at = arg maxa∈{1,2,··· ,K} r̂(st, a) and get the reward r(st, at), update the posterior µt+1st,at with r(st, at).
end for
Notice that all of the reinforcement learning problems face the exploration-exploitation dilemma, so as the contextual bandit problem. Posterior sampling trade off the exploration and exploitation with the uncertainty provided by the posterior distribution. So if the posterior uncertainty is not estimated properly, posterior sampling will perform poorly. To see this, if we over-estimate the uncertainty, we can explore too-much sub-optimal actions, while if we under-estimate the uncertainty, we can
fail to find the optimal actions. Thus, it is a good benchmark for evaluating the uncertainty provided by different inference methods.
Here, we test the uncertainty provided by vanilla Langevin dynamics and Self-repulsive Langevin dynamics on two of the benchmark contextual bandit problems suggested by (Riquelme et al., 2018), called mushroom and wheel. One can read (Riquelme et al., 2018) to find the detail introduction of this two contextual bandit problems. For completeness, we include it as follows:
Mushroom Mushroom bandit utilizes the data from Mushroom dataset (Schlimmer, 1981), which includes different kinds of poisonous mushroom and safe mushroom with 22 attributes that can indicate whether the mushroom is poisonous or not. Blundell et al. (2015) first introduced the mushroom bandit by designing the following reward function: eating a safe mushroom will give a +5 reward, while eating a poisonous mushroom will return a reward +5 and −35 with equal chances. The agent can also choose not to eat the mushroom, which always yield a 0 reward. Same to (Riquelme et al., 2018), we use 50000 instances in this problem.
Wheel To highlight the need for exploration, (Riquelme et al., 2018) designs the wheel bandit, that can control the need of exploration with some “exploration parameter” δ ∈ (0, 1). The context set S is the unit circle ‖s‖2 ≤ 1 in R2, and each turn the context st is uniformly sampled from S. K = 5 possible actions are provided: the first action yields a constant reward r ∼ N (µ1, σ2); the reward corresponding to other actions is determined by the provided context s:
• For s ∈ S s.t. ‖s‖2 ≤ δ, all of the four other actions return a suboptimal reward sampled from N (µ2, σ2) for µ2 < µ1.
• For s ∈ S s.t. ‖s‖2 > δ, according to the quarter the context s is in, one of the four actions becomes optimal. This optimal action gives a reward of N (µ3, σ2) for µ3 µ1, and another three actions still yield the suboptimal reward N (µ2, σ2).
Following the setting from (Riquelme et al., 2018), we set µ1 = 1.2, µ2 = 1.0, and µ3 = 50.
When δ approaches 1, the inner circle ‖s‖2 ≤ δ will dominate the unit circle and the first action becomes the optimal for most of the context. Thus, inference methods with poorly estimated uncertainty will continuously choose the suboptimal action a1 for all of the contexts without exploration. This phenomenon have been confirmed in (Riquelme et al., 2018). In our experiments, as we want to evaluate the quality of uncertainty provided by different methods, we set δ = 0.95, which is pretty hard for existing inference methods as shown in (Riquelme et al., 2018), and use 50000 contexts for evaluation.
Experiment Setup Following (Riquelme et al., 2018), we use a feed-forward network with two hidden layer of 100 units and ReLU activation. We use the same step-size and thinning factor c/η = 100 for vanilla Langevin dynamics and SRLD, and set M = 20, α = 1 on both of the mushroom and wheel bandits. The update schedule is similar to (Riquelme et al., 2018), and we just change the optimization step in stochastic variational inference methods into MCMC sampler step and use a higher (2000) initial steps for burn-in like the warm-up of stochastic variational inference methods in Riquelme et al. (2018). As this is an online posterior inference problem, we only use the
last 20 samples to give the prediction. Notice that, in the original implementation of Riquelme et al. (2018), the authors only update a few steps with new observation after observing enough data, as the posterior will gradually converge to the true reward distribution and little update is needed after observing sufficient data. Similar to their implementation, after observing enough data, we only collect one new sample with the new observation each time. For SVGD, we also use 20 particles to make the comparison fair.

heading: B THE DETAILED ANALYSIS OF SRLD

heading: B.1 SOME ADDITIONAL NOTATION
text: We use ‖ · ‖∞ to denote the `∞ vector norm and define the L∞ norm of a function f : Rd → R1 as ‖f‖L∞ . DTV denote the Total Variation distance between distribution ρ1, ρ2 respectively. Also, as K is Rd × Rd → R1, we denote ‖K‖L∞,L∞ = supθ1,θ2 K(θ1,θ2). For simplicity, we may use ‖K‖∞,∞ as ‖K‖L∞,L∞ . In the appendix, we also use φ[ρ](θ) := g(θ; ρ), where g(θ; ρ) is defined in the main text. For the clearance, we define πM,c/η ∗ ρk := ρMk , πM,c/η ∗ ρ̃k := ρ̃Mk and πM,c ∗ ρ̄t := ρ̄Mt , where ρMk , ρ̃Mk and ρ̄Mt are defined in main text.

heading: B.2 GEOMETRIC ERGODICITY OF SRLD
text: Before we start the proof of main theorems, we give the following theorem on the geometric ergodicity of SRLD. It is noticeable that under this assumption, the practical dynamics follows an (Mc/η + 1)-order nonlinear autoregressive model when k ≥Mc/η:
θk+1 = ψ ( θk, ...,θk−Mc/η ) + √ 2ηek,
where
ψ ( θk, ...,θk−Mc/η ) = θk + η −∇V (θk) + αφ[ 1M M∑ j=1 δθk−jc/η ](θk)  . Further, if we stack the parameter by Θk = [ θk, ...,θk−Mc/η ]> and define Ψ (Θk) =[
ψ> (Θk) ,Θ > k
]> , we have
Θk+1 = Ψ (Θk) + √ 2ηEk,
where Ek = [ e>k ,0 >, ...,0> ]> . In this way, we formulate Θk as a time homogeneous Markov Chain. In the following analysis, we only analyze the second phase of SRLD given some initial stacked particles ΘMc/η−1. Theorem B.1. (Geometric Ergodicity) Under Assumption 4.1 and Assumption 4.2, suppose we choose η and α such that
max ( 1− 2ηa1 + η2b1 + 2αη
σ b1,
2αη
σ (b1 + 1)
) < 1,
then the Markov Chain of Θk is stationary, geometrically ergodic, i.e., for any Θ′0 = ΘMc/η−1, we have DTV [ P k (·,Θ0) ,Π (·) ] ≤ Q (Θ0) e−rk,
where r is some positive constant, Q(Θ0) is constant related to Θ0, P k is the k-step Markov transition kernel and Π is the stationary distribution.
We defer the proof to Appendix B.5.1.

heading: B.3 MOMENT BOUND
text: Theorem B.2. (Moment Bound) Under assumption 4.2, suppose that we have Eθ∼ρ0 ‖θ‖ 2 < ∞; and a2 − α ‖K‖∞ ( 2b1 + 4 σ ) > 0, we have
sup k
Eθ∼ρk ‖θ‖ 2 ∨ sup k Eθ∼ρ̃k ‖θ‖ 2 ∨ sup t Eθ∼ρ̄t ‖θ‖ 2
≤Eθ∼ρ0 ‖θ‖ 2
+ b1 + 1 + η
a2 − ‖K‖L∞,L∞ 2α σ − α ‖K‖L∞,L∞
( 2b1 + 2 σ ) . And by Lemma B.1, we thus have
sup k
Eθ∼ρk ‖∇V (θ)‖ 2 ∨ sup k Eθ∼ρ̃k ‖∇V (θ)‖ 2 ∨ sup t Eθ∼ρ̄t ‖∇V (θ)‖ 2
≤b1Eθ∼ρ0 ‖θ‖ 2
+ b1(b1 + 1 + η)
a2 − ‖K‖L∞,L∞ 2α σ − α ‖K‖L∞,L∞
( 2b1 + 2 σ ) + 1
The proof can be found at Appendix B.5.2.

heading: B.4 TECHNICAL LEMMA

heading: Definition B.1. (α-mixing)
text: For any two σ-algebras A and B, the α-mixing coefficient is defined by α(A,B) = sup
A∈A,B∈B |P (A ∩B)− P (A)P (B)| .
Let (Xk, k ≥ 1) be a sequence of real random variable defined on (Ω,A,P). This sequence is α-mixing if
α(n) := sup k≥1 α (Mk,Gk+n)→ 0, as n→∞,
whereMj := σ (Xi, i ≤ j) and Gj := σ (Xi, i ≥ j) for j ≥ 1. Alternatively, as shown by Theorem 4.4 of Bradley (2007)
α(n) := 1
4 sup
{ Cov (f, g)
‖f‖L∞ ‖g‖L∞ ; f ∈ L∞ (Mk) , g ∈ L∞ (Gk+n)
} .

heading: Definition B.2. (β-mixing)

heading: For any two σ-algebras A and B, the α-mixing coefficient is defined by
text: β (A,B) := sup 1 2 I∑ i=1 J∑ j=1 |P (Ai ∩Bj)− P (Ai)P (Bj)| ,
where the supremum is taken over all pairs of finite partitions {A1, ..., AI} and {B1, ..., BJ} of Ω such that Ai ∈ A and Bj ∈ B for each i, j. Let (Xk, k ≥ 1) be a sequence of real random variable defined on (Ω,A,P). This sequence is β-mixing if
β(n) := sup k≥1 β (Mk,Gk+n)→ 0, as n→∞.
Proposition B.1. (β-mixing implies α-mixing)

heading: For any two σ-algebras A and B,
text: α (A,B) ≤ 1 2 β (A,B) .
This proposition can be found in Equation 1.11 of Bradley (2005). Proposition B.2. A (strictly) stationary Markov Chain is geometric ergodicity if and only if β(n)→ 0 at least exponentially fast as n→∞.
This proposition is Theorem 3.7 of Bradley (2005). Lemma B.1. By Assumption 4.2, we have ‖∇V (θ)‖ ≤ b1 (‖θ1‖+ 1) and ‖θ − η∇V (θ)‖ ≤( 1− 2ηa1 + η2b1 ) ‖θ‖2 + η2b1 + 2ηb1. Lemma B.2. (Some property of RBF Kernel) For RBF kernel with bandwidth σ, we have ‖K‖∞,∞ ≤ 1 and
‖K(θ′,θ1)−K(θ′,θ2)‖ ≤ ∥∥∥e−(·)2/σ∥∥∥
Lip ‖θ1 − θ2‖2 ‖∇θ′K(θ′,θ1)−∇θ′K(θ′,θ2)‖ ≤ ∥∥∥∥ 2σ e−(·)2/σ(·) ∥∥∥∥ Lip ‖θ1 − θ2‖2 .
Lemma B.3. (Some property of Stein Operator)
For any distribution ρ such that Eθ∼ρ ‖∇V (θ)‖ <∞, we have ‖φ[ρ](·)‖Lip ≤ ∥∥∥e−(·)2/σ∥∥∥
Lip Eθ∼ρ ‖∇V (θ)‖+ ∥∥∥∥ 2σ e−(·)2/σ(·) ∥∥∥∥
Lip
,
‖φ[ρ](θ)‖ ≤ ‖K‖∞ Eθ′∼ρ [ ‖∇V (θ′)‖+ 2
σ (‖θ′‖+ ‖θ‖) ] ≤ ‖K‖∞ b1 + Eθ′∼ρ [( 2
σ + b1
) ‖θ′‖ ] + ‖θ‖ .
Lemma B.4. (Bounded Lipschitz of Stein Operator) Given θ′, define φ̄θ′(θ) := φ[δθ′ ](θ) = K(θ′,θ)∇V (θ′) +∇1K(θ′,θ). We also denote φ̄θ′(θ) = [φ̄θ′,1(θ), ..., φ̄θ′,d(θ)]>. We have
d∑ i=1 ∥∥φ̄θ′,i(θ)∥∥2Lip ≤ 2 ‖∇V (θ′)‖2 ∥∥∥e−‖·‖2/σ∥∥∥2Lip + 2d ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2 Lip
d∑ i=d ∥∥φ̄θ′,i(θ)∥∥2L∞ ≤ 2d ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2 L∞ + 2 ∥∥∥e−‖·‖2/σ∥∥∥2 L∞ ‖∇V (θ′)‖2 .

heading: B.5 PROOF OF MAIN THEOREMS

heading: B.5.1 PROOF OF THEOREM B.1
text: The proof of this theorem is by verifying the condition of Theorem 3.2 of An & Huang (1996). Suppose Θ = [θ1, ...,θMC+1], where C = c/η, we have
‖ψ (Θ)‖ = ∥∥∥∥∥∥θ1 + η −∇V (θ1) + αφ[ 1M M∑ j=1 δθ1+jC ](θk)  ∥∥∥∥∥∥
= ∥∥∥∥∥∥θ1 − η∇V (θ1) + ηαM M∑ j=1 [ e−‖θ1+jC−θ1‖ 2/σ 2 σ (θ1 − θ1+jC)− e−‖θ1+jC−θ1‖ 2/σ∇V (θ1+jC) ]∥∥∥∥∥∥
≤ ∥∥∥∥∥∥θ1 − η∇V (θ1) + 2σ ηαM M∑ j=1 e−‖θ1+jC−θ1‖ 2/σθ1 ∥∥∥∥∥∥ +
∥∥∥∥∥∥ηαM M∑ j=1 e−‖θ1+jC−θ1‖ 2/σ 2 σ (−∇V (θ1+jC)− θ1+jC) ∥∥∥∥∥∥ ≤ ‖θ1 − η∇V (θ1)‖+ 2αη
σ ‖K‖∞,∞ b1(1 + ‖θ1‖)
+ 2αη
Mσ M∑ j=1 ‖K‖∞,∞ b1 ( 1 + (1 + 1 b1 ) ‖θ1+jC‖ ) (1) ≤ b1(1 + 4αη
σ ‖K‖∞,∞) + η 2b1 + 2ηb1
+ ( 1− 2ηa1 + η2b1 + 2αη
σ ‖K‖∞,∞ b1
) ‖θ1‖+ 2αη
σ ‖K‖∞,∞ (b1 + 1) max i∈[MC+1]−{1} ‖θ1+jC‖
≤ b1(1 + 4αη
σ ‖K‖∞,∞) + η 2b1 + 2ηb1
+ max ( 1− 2ηa1 + η2b1 + 2αη
σ ‖K‖∞,∞ b1,
2αη
σ ‖K‖∞,∞ (b1 + 1)
) max
i∈[MC+1] ‖θ1+jC‖ ,
where (1) is by Lemma B.1. Thus, given the step size η, if we choose η, α such that
max ( 1− 2ηa1 + η2b1 + 2αη
σ ‖K‖∞,∞ b1,
2αη
σ ‖K‖∞,∞ (b1 + 1)
) < 1,
then our dynamics is geometric ergodic.

heading: B.5.2 PROOF OF THEOREM B.2
text: Continuous-Time Mean Field Dynamics (5) Notice that as our dynamics has two phases and the first phase can be viewed as an special case of the second phase by setting α = 0, here we only
analysis the second phase. Define Ut = sup s≤t
√ E ∥∥θ̄s∥∥2, and thus
∂ ∂t U2t ≤ E
〈 θ̄t,−V (θ̄) + αφ[πM,c ∗ ρ̄t](θ̄t) 〉 ∨ 0.
Now we bound E 〈 θ̄t,−V (θ̄) + αφ[πM,c ∗ ρ̄t](θ̄t) 〉 :
E 〈 θ̄t,−V (θ̄t) + αφ[πM,c ∗ ρ̄t](θ̄t) 〉 ≤b1 − a2E
∥∥θ̄t∥∥2 + αE ∥∥θ̄t∥∥∥∥φ[πM,c ∗ ρ̄t](θ̄t)∥∥ (1)
≤b1 − a2E ∥∥θ̄t∥∥2 + α ‖K‖∞ E ∥∥θ̄t∥∥Eθ′∼πM,c∗ρ̄t [‖∇V (θ′)‖+ 2σ (‖θ′‖+ ‖θt‖) ] ≤b1 − a2E ∥∥θ̄t∥∥2 + α ‖K‖∞ E∥∥θ̄t∥∥Eθ′∼πM,c∗ρ̄t [b1 (‖θ′‖+ 1) + 2σ (‖θ′‖+ ‖θt‖) ]
=b1 − ( a2 − ‖K‖∞ 2α
σ ) E ∥∥θ̄t∥∥2 + α ‖K‖∞ E ∥∥θ̄t∥∥Eθ′∼πM,c∗ρ̄t ((b1 + 2σ ) ‖θ′‖+ b1 ) ≤b1 − ( a2 − ‖K‖∞ 2α
σ
) U2t + α ‖K‖∞ E ∥∥θ̄t∥∥Eθ′∼πM,c∗ρ̄t ((b1 + 2σ ) ‖θ′‖+ b1 ) ≤b1 − ( a2 − ‖K‖∞ 2α
σ
) U2t + α ‖K‖∞ ( b1 + 2
σ
) 1
M M∑ j=1 UtUt−jc + α ‖K‖∞ b1Ut
≤b1 − ( a2 − ‖K‖∞ 2α
σ
) U2t + α ‖K‖∞ ( b1 + 2
σ
) U2t + α ‖K‖∞ b1(U 2 t + 1)
≤ (b1 + 1)− ( a2 − ‖K‖∞ 2α
σ − α ‖K‖∞
( 2b1 + 2
σ
)) U2t ,
where (1) is by B.3. By the assumption that λ := a2 − ‖K‖∞ 2α σ − α ‖K‖∞ ( 2b1 + 2 σ ) > 0, we have ∂
∂t U2t ≤
[ (b1 + 1)− λU2t ] ∨ 0.
By Gronwall’s inequality, we have U2t ≤ U20 + b1+1λ . (If ∂ ∂tU 2 t = 0, then Ut fix and this bound still holds.) Notice that in the first phase, as α = 0, we have λ < a2 and thus this inequality also holds.
Discrete-Time Mean Field Dynamics (4) Similarly to the analysis of the continuous-time mean
field dynamics (5), we only give proof of the second phase. Define Uk = sup s≤k
√ E ∥∥∥θ̃s∥∥∥2, and thus
U2k − U2k−1 ≤ [ 2ηE 〈 θ̃k−1,−∇V (θ̃k) + αφ[πM,c/η ∗ ρ̃k](θ̃k) 〉 + 2η2 ] ∨ 0.
By a similarly analysis, we have bound E 〈 θ̃k−1,−∇V (θ̃k) + αφ[πM,c/η ∗ ρ̃k](θ̃k) 〉 ≤ (b1 + 1)− λU2t ,
where λ = a2 − ‖K‖∞,∞ 2α σ − α ‖K‖∞,∞ ( 2b1 + 2 σ ) > 0. And thus we have
U2k − U2k−1 ≤ [ 2η [ (b1 + 1)− λU2k−1 ] + 2η2 ] ∨ 0.
It gives that
U2k ≤ b1 + 1 + η
λ + U20 .
Practical Dynamics (3) The analysis of Practical Dynamics (3) is almost identical to that of the discrete-time mean field dynamics (4) and thus is omitted here.

heading: B.5.3 PROOF OF THEOREM 4.1 AND 5.1
text: Notice that the dynamics in Theorem 4.1 is special case of that in Theorem 5.1 and thus we only prove Theorem 5.1 here. After some algebra, we can show that the continuity equation of dynamics (6) is
∂tρt = ∇ · ([− (D(θ) +Q(θ))∇V (θ) + αφ[πM,c ∗ ρt](θt)] ρt + (D(θ) +Q(θ))∇ρt) .
Notice that the limiting distribution satisfies
0 a.e. = ∇ · ([− (D(θ) +Q(θ))∇V (θ) + αφ[πM,c ∗ ρ∞](θt)] ρ∞ + (D(θ) +Q(θ))∇ρ∞) = ∇ · ([− (D(θ) +Q(θ))∇V (θ) + αφ[ρ∞](θt)] ρ∞ + (D(θ) +Q(θ))∇ρ∞) = ∇ · ([− (D(θ) +Q(θ))∇V (θ)] ρ∞ + (D(θ) +Q(θ))∇ρ∞) + α∇ · (K ∗ (∇ρ∞ −∇V (θ)ρ∞) ρ∞) .
which implies that ρ∞ ∝ exp(−V (θ)) is the stationary distribution.

heading: B.5.4 PROOF OF THEOREM 4.2
text: In the later proof we use cd to represent the quantity√ Eθ∼ρ0 ‖θ‖ 2 + b1 + 1 + η
a2 − ‖K‖∞,∞ 2α σ − α ‖K‖∞,∞ ( 2b1 + 2 σ ) . Recall that there are two dynamics: the continuous-time mean field dynamics (5) and the discretized version discrete-time mean field Dynamics (4). Notice that here we couple the discrete-time mean field dynamics with the continuous-time mean field system using the same initialization. Given any T = ηN , for any 0 ≤ t ≤ T , define t = b tη cη. We introduce an another continuous-time interpolation dynamics:
θ̂t = { −∇V (θ̂t) + dBt, t ∈ [0,Mc) −∇V (θ̂t) + αφ[πM,c ∗ ρ̂t](θ̂t) + dBt, t ≥Mc,
ρ̂t = Law(θ̂t),
θ̂0 = θ̄0 ∼ ρ̄0,
Notice that here we couples this interpolation dynamics with the same Brownian motion as that of the dynamics of θ̄t. By the definition of θ̂t, at any tk := kη for some integrate k ∈ [N ], θ̂tk and θ̃k has the same distribution. Define ρ̄θ0t = Law(θ̄t) conditioning on θ̄0 = θ0 and ρ̂ θ0 t = Law(θ̂t) conditioning on θ̂0 = θ0. Followed by the argument of proving Lemma 2 in Dalalyan (2017), if k ≥ Mcη , we have
DKL [ ρ̄θ0tk ‖ ρ̂ θ0 tk ] = 1
4 ∫ tk 0 E ∥∥∥−∇V (θ̂s) + αφ[πM,c ∗ ρ̂s](θ̂s) +∇V (θ̂s)− αφ[πM,c ∗ ρ̄s](θ̂s)∥∥∥2 ds
= 1
4 k−1∑ j=0 ∫ tj+1 tj E ∥∥∥−∇V (θ̂tj ) + αφ[πM,c ∗ ρ̂tj ](θ̂tj ) +∇V (θ̂s)− αφ[πM,c ∗ ρ̄s](θ̂s)∥∥∥2 ds
≤3 4 k−1∑ j=0 ∫ tj+1 tj E ∥∥∥∇V (θ̂tj )−∇V (θ̂s)∥∥∥2 ds
+ 3α2
4 k−1∑ j=0 ∫ tj+1 tj E ∥∥∥φ[πM,c ∗ ρ̂tj ](θ̂tj )− φ[πM,c ∗ ρ̄s](θ̂tj )∥∥∥2 ds
≤3α 2
4 k−1∑ j=0 ∫ tj+1 tj E ∥∥∥φ[πM,c ∗ ρ̄s](θ̂tj )− φ[πM,c ∗ ρ̄s](θ̂s)∥∥∥2 ds
=I1 + I2 + I3.
We bound I1, I2 and I3 separately.
Bounding I1 and I3 By the smoothness of∇V , we have∥∥∥∇V (θ̂tj )−∇V (θ̂s)∥∥∥2 ≤ b21 ∥∥∥θ̂tj − θ̂s∥∥∥2 .
And by Lemma B.3 (Lipschitz of Stein Operator), we know that ‖φ[πM,c ∗ ρ̄s](θ1)− φ[πM,c ∗ ρ̄s](θ2)‖
≤ [∥∥∥e−(·)2/σ∥∥∥ Lip Eθ∼πM,c∗ρ̄s ‖∇V (θ)‖+ ∥∥∥∥ 2σ e−(·)2/σ(·) ∥∥∥∥ Lip ] ‖θ1 − θ2‖2 .
And by the Assumption 4.2 and that ρ̄s as finite second moment, we have ‖φ[πM,c ∗ ρ̄s](θ1)− φ[πM,c ∗ ρ̄s](θ2)‖ ≤Ccd ‖θ1 − θ2‖2 . Combine the two bounds, we have
I1 + I3 ≤ 3Cc2d
4 k−1∑ j=0 ∫ tj+1 tj E ∥∥∥θ̂tj − θ̂s∥∥∥2 ds.
Notice that θ̂t = θ̂t+ [ −∇V (θ̂t) + αφ[πM,c ∗ ρ̂t](θ̂t) ] (t− t)+ ∫ t t dBs. By Itô’s lemma, it implies that
I1 + I3 ≤ 3Cc2d
4 k−1∑ j=0 ∫ tj+1 tj E ∥∥∥θ̂tj − θ̂s∥∥∥2 ds
≤ 3Cc 2 d
4 ∫ tj+1 tj [ E ∥∥∥−∇V (θ̂s) + αφ[πM,c ∗ ρ̂s](θ̂s)∥∥∥2 (s− tj)2 + 2d(s− tj)] ds
= Cc2dη 3 k−1∑ j=0 E ∥∥∥−∇V (θ̂tj ) + αφ[πM,c ∗ ρ̂tj ](θ̂tj )∥∥∥2 + Cc2ddkh2.
By the assumption that E ∥∥∥θ̃tj∥∥∥ is finite and θ̃tj d= θ̂tj , E∥∥∥θ̂tj∥∥∥2 is also finite, we have
E ∥∥∥−∇V (θ̂t) + αφ[πM,c ∗ ρ̂t](θ̂t)∥∥∥2
≤2E ∥∥∥∇V (θ̂t)∥∥∥2 + 2α2E∥∥∥φ[πM,c ∗ ρ̂t](θ̂t)∥∥∥2
≤4b21 + 4b21E ∥∥∥θ̂t∥∥∥2 + 2α2E(( 2
σ + b1
) Eθ′∼πM,c∗ρ̂t ‖θ′‖+ ‖θ‖ )2 ≤c2dC.
Thus we conclude that I1 + I3 ≤ Cc2d ( c2dkη 3 + dkη2 ) .
Bounding I2 E ∥∥∥φ[πM,c ∗ ρ̂tj ](θ̂tj )− φ[πM,c ∗ ρ̄s](θ̂tj )∥∥∥2
=E ∥∥∥∥∥ 1M M∑ l=1 [ φ[ρ̂tj−cl](θ̂tj )− φ[ρ̄s−cl](θ̂tj ) ]∥∥∥∥∥ 2
≤ 1 M M∑ l=1 E ∥∥∥φ[ρ̂tj−cl](θ̂tj )− φ[ρ̄s−cl](θ̂tj )∥∥∥2
= 1
M M∑ l=1 E ∥∥∥Eθ∼ρ̂tj−cl φ̄θ̂tj (θ)− Eθ∼ρ̄s−cl φ̄θ̂tj (θ)∥∥∥2
= 1
M M∑ l=1 Eθ̂tj d∑ i=1 ∣∣∣Eθ∼ρ̂tj−cl φ̄θ̂tj ,i(θ)− Eθ∼ρ̄s−cl φ̄θ̂tj ,i(θ)∣∣∣2 ≤ 1 M M∑ l=1 Eθ̂tj d∑ i=1 (∥∥∥φ̄θ̂tj ,i(·)∥∥∥L∞ ∨ ∥∥∥φ̄θ̂tj ,i(·)∥∥∥Lip )2 D2BL [ ρ̂tj−cl, ρ̄s−cl ]
By Lemma B.4 and the Assumption 4.4 that V is at most quadratic growth and that ρ̂t has finite second moment, we have
Eθ̂tj d∑ i=1 (∥∥∥φ̄θ̂tj ,i(·)∥∥∥L∞ ∨ ∥∥∥φ̄θ̂tj ,i(·)∥∥∥Lip )2
=Eθ̂tj d∑ i=1 (∥∥∥φ̄θ̂tj ,i(·)∥∥∥2L∞ ∨ ∥∥∥φ̄θ̂tj ,i(·)∥∥∥2Lip )
≤ [ 4d ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2
BL
+ 4 ∥∥∥e−‖·‖2/σ∥∥∥2
BL Eθ̂tj
∥∥∥∇V (θ̂tj )∥∥∥2 ]
≤C(d+ c2d).
Plug in the above estimation, we have
I2 = 3α2
4 k−1∑ j=0 ∫ tj+1 tj E ∥∥∥φ[πM,c ∗ ρ̂tj ](θ̂tj )− φ[πM,c ∗ ρ̄s](θ̂tj )∥∥∥2 ds
≤ α2C(d+ c2d) k−1∑ j=0 ∫ tj+1 tj 1 M M∑ l=1 D2BL [ ρ̂tj−cl, ρ̄s−cl ] ds ≤ α2C(d+ c2d) k−1∑ j=0 1 M M∑ l=1 ∫ tj+1 tj DKL [ ρ̂tj−cl, ρ̄s−cl ] ds,
where the last inequality is due to the relation that D2BL definition ≤ D2TV Pinsker′s ≤ DKL.
Overall Bound Combine all the estimation, we have DKL [ ρ̄θ0tk ‖ ρ̂ θ0 tk ] ≤ α2C(d+ c2d)
k−1∑ j=0 1 M M∑ l=1 ∫ tj+1 tj DKL [ ρ̂tj−cl, ρ̄s−cl ] ds+ Cc2d ( c2dkη 3 + dkη2 )
= α2C(d+ c2d) k−1∑ j=0 1 M M∑ l=1 ∫ η 0 DKL [ ρ̂t ( jη−clη ) , ρ̄t ( jη−clη ) +s ] ds+ Cc2d ( c2dkη 3 + dkη2 )
Similar, if k ≤ Mcη − 1, we have
DKL [ ρ̄θ0tk ‖ ρ̂ θ0 tk ] = 1
4 ∫ tk 0 E ∥∥∥∇V (θ̂s)−∇V (θ̂s)∥∥∥2 ds
≤b 2 1
4 k−1∑ j=0 ∫ tj+1 tj E ∥∥∥θ̂tj − θ̂s∥∥∥2 ds
≤b 2 1η 3
12 k−1∑ j=0 E ∥∥∥∇V (θ̂tj )∥∥∥2 + dkb21η24
≤b 2 1η 3kc2d 12 + dkb21η 2 4 .
Define
uk = sup s∈[tk,tk+1]
DKL [ ρ̄θ0s ‖ ρ̂θ0s ] ,
and Uk = max l∈{0,...,k}
ul. We conclude that for k ≥ Mcη , for any k ′ ≤ k,
uk′ ≤ α2C(d+ c2d) k−1∑ j=0 1 M M∑ l=1 ∫ h 0 DKL [ ρ̂t ( jη−clη ) , ρ̄t ( jη−clη ) +s ] ds+ Cc2d ( c2dkη 3 + dkη2 )
≤ α2C(d+ c2d) k−1∑ j=0 1 M M∑ l=1 ηu( jη−clη ) + Cc2d ( c2dkη 3 + dkη2 ) ≤ α2C(d+ c2d)η k−1∑ j=0 Uj + Cc 2 d ( c2dkη 3 + dkη2 ) .
For k < Mcη , which is a simpler case, we have Uk ≤ C ( η3kc2d + dkη 2 ) < CMc ( ηc2d + d ) η.
We bound the case when k ≥ Mcη ,
Uk ≤ α2C(d+ c2d)η k−1∑ j=0 Uj + Cc 2 d ( c2dkη 3 + dkη2 ) .
If we take η sufficiently small, such that c2dkη 3 ≤ dkη2, we have
Uk ≤ α2C(d+ c2d)η k−1∑ j=0 Uj + 2Cc 2 ddkη 2
≤ α2C(d+ c2d)η k−1∑ j=0 (Uj + η) .
Define η′ = α2C(d + c2d)η and we can choose η small enough such that η ′ < 1/2 and η < 1/2. Without loss of generality, we also assume η′ ≥ η and thus we have
Uk ≤ η′ k−1∑ j=0 (Uj + η ′) .
Also we assume Uk ≥ η′, otherwise we conclude that Uk < η′. We thus have Uk ≤ q ∑k−1 j=0 Uj , where q = 2η′. Suppose that UMc η −1 = x ≤ CMc ( ηc2d + d ) η and some algebra (which reduces to Pascal’s triangle) shows that Uk ≤ xq(1 + q)k− Mc η . We conclude that Uk ≤ xq(1 + q)k−1. Notice that q = 2α2C(d+ c2d)η. Thus for any k ≥Mc/η,
Uk ≤ xq(1 + q)k− Mc η
= xq(1 + q)(kη−Mc)/η
= xq(1 + q)2α 2C(d+c2d)(kη−Mc)/q
≤ x2α2C(d+ c2d)e2α 2C(d+c2d)(kη−Mc)η ≤ CMcα2 ( ηc2d + d ) (d+ c2d)e 2α2C(d+c2d)(kη−Mc)η2,
for sufficiently small η. Combine the above two estimations, we have
Uk ≤
{ C ( η3kc2d + dkη 2 + η )
k ≤Mc/η − 1 CMcα2 ( ηc2d + d ) (d+ c2d)e 2α2C(d+c2d)(kη−Mc)η2 + Cη k ≥Mc/η .
Notice that now we have Uk = max l∈{0,...,k} sup s∈[0,η]
DKL [ ρ̄θ0lη+s ‖ ρ̃ θ0 lη ] , which is a function of θ0.
We then bound Ūk = max l∈{0,...,k} sup s∈[0,η] DKL [ρ̄lη+s ‖ ρ̃lη]. Notice that the KL divergence has the following variational representation:
DKL[ρ1 ‖ ρ2] = sup f
[ Eρ1f − Eρ2ef ] ,
where the f is chosen in the set that Eρ1f and Eρ2ef exist. And thus we have
DKL[ρ̄lη+s ‖ ρ̃lη] = sup f
[ Eθ0∼ρ0 ( E ρ̄ θ0 lη+s f − E ρ̃ θ0 lη ef )]
≤ Eθ0∼ρ0 sup f [( E ρ̄ θ0 lη+s f − E ρ̃ θ0 lη ef )] .
And thus Ūk ≤ Uk. Also the inequality that
Ūk = max l∈{0,...,k} sup s∈[0,η] DKL [ρ̄lη+s ‖ ρ̃lη] ≥ max l∈{0,...,k} DKL [ρ̄lη ‖ ρ̃lη]
holds naturally by definition. We complete the proof.

heading: B.5.5 PROOF OF THEOREM 4.3
text: The constant h1 is defined as
h1 = ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2
BL
∨ ∥∥∥e−‖·‖2/σ∥∥∥2 BL ∨ ∥∥∥∥ 2σ e−(·)2/σ(·) ∥∥∥∥ Lip
Now we start the proof. We couple the process of θk and θ̃k by the same gaussian noise ek in every iteration and same initialization θ̃0 = θ0. For k ≤Mc/η−1, E ∥∥∥θk − θ̃k∥∥∥2 = 0 and for k ≥Mc/η
we have the following inequality, E ∥∥∥θk+1 − θ̃k+1∥∥∥2 − E∥∥∥θk − θ̃k∥∥∥2
=2ηE 〈 θk − θ̃k,−∇V (θk) +∇V (θ̃k) 〉 +2ηαE 〈 θk − θ̃k, φ[ 1
M M∑ j=1 δθk−jc/η ](θk)− φ[πM,c/η ∗ ρ̃k](θ̃k)
〉
+η2E ∥∥∥∥∥∥−∇V (θk) + αφ[ 1M M∑ j=1 δθk−jc/η ](θk) +∇V (θ̃k)− αφ[πM,c/η ∗ ρ̃k](θ̃k) ∥∥∥∥∥∥ 2 .
By the log-concavity, we have E 〈 θk − θ̃k,−∇V (θk) +∇V (θ̃k) 〉 ≤− LE
∥∥∥θk − θ̃k∥∥∥2 , for some positive constant L. And also, as η is small, the last term on the right side of the equation is small term. Thus our main target is to bound the second term. We decompose the second term on the left side of the equation by
E 〈 θk − θ̃k, φ[ 1
M M∑ j=1 δθk−jc ](θk)− φ[πM,c/η ∗ ρ̃k](θ̃k)
〉
=E 〈 θk − θ̃k, φ[ 1
M M∑ j=1 δθk−jc/η ](θk)− φ[πM,c/η ∗ ρk](θk)
〉
+E 〈 θk − θ̃k, φ[πM,c/η ∗ ρk](θk)− φ[πM,c/η ∗ ρ̃k](θk) 〉 +E 〈 θk − θ̃k, φ[πM,c/η ∗ ρ̃k](θk)− φ[πM,c/η ∗ ρ̃k](θ̃k)
〉 =I1 + I2 + I3.
We bound I1, I2 and I3 independently.
Bounding I1 By Holder’s inequality,
I1 ≤ E ∥∥∥θk − θ̃k∥∥∥ ∥∥∥∥∥∥φ[ 1M M∑ j=1 δθk−jc/η ](θk)− φ[πM,c/η ∗ ρk](θk) ∥∥∥∥∥∥ 
≤ √ E ∥∥∥θk − θ̃k∥∥∥2 √√√√√E ∥∥∥∥∥∥φ[ 1M M∑ j=1 δθk−jc/η ](θk)− φ[πM,c/η ∗ ρk](θk) ∥∥∥∥∥∥ 2 .
We bound the second term on the right side of the inequality. Define
a2 = sup k
E ∥∥∥φ[ 1M ∑Mj=1 δθk−jc/η ](θk)− φ[πM,c/η ∗ ρk](θk)∥∥∥2
sup ‖θ‖≤B
E ∥∥∥φ[ 1M ∑Mj=1 δθk−jc/η ](θ)− φ[πM,c/η ∗ ρk](θ)∥∥∥2
and by the regularity assumption we know that
a2 ≤ sup k
E ∥∥∥φ[ 1M ∑Mj=1 δθk−jc ](θk)− φ[πM,c/η ∗ ρk](θk)∥∥∥2
inf k sup ‖θ‖≤B
E ∥∥∥φ[ 1M ∑Mj=1 δθk−jc/η ](θ)− φ[πM,c/η ∗ ρk](θ)∥∥∥2 <∞.
Define φ[ 1M ∑M j=1 δθk−jc/η ](θ) − φ[πM,c/η ∗ ρk](θ) = φ∗[ 1 M ∑M j=1 δθk−jc/η ] and since the stein operator is linear functional of the distribution, we have
Eφ∗[ 1
M M∑ j=1 δθk−jc/η ](θ) = 0,
given any θ. By Theorem B.1 that Θk is geometric ergodicity and thus is β-mixing with exponentially fast decay rate by Proposition B.2. And by Proposition B.1, we know that Θk is also α-mixing with exponentially fast decay rate. We have the following estimation
E ∥∥∥∥∥∥φ[ 1M M∑ j=1 δθk−jc/η ](θk)− φ[πM,c/η ∗ ρk](θk) ∥∥∥∥∥∥ 2
≤a2 sup ‖θ‖≤B E ∥∥∥∥∥∥φ∗[ 1M M∑ j=1 δθk−jc/η ](θ) ∥∥∥∥∥∥ 2
≤ a2 M2 sup ‖θ‖≤B E M∑ k=1 ∥∥φ∗[δθt−kc/η ](θ)∥∥2 + a2 M2 sup ‖θ‖≤B E ∑ k 6=j 〈 φ∗[δθt−kc/η ](θ), φ ∗[δθt−jc/η ](θ) 〉
≤Ca2 M
[ e−rc ( 1− e−rMc ) 1− erc + 1 ] ,
for some positive constant r that characterize the decay rate of α mixing. Combine this two estimations, we have
I1 ≤ √ E ∥∥∥θk − θ̃k∥∥∥2 √ a2C
M
[ e−rc (1− e−rMc)
1− erc + 1
] .
Bounding I2 By Holder’s inequality, we have I2 ≤ √ E ∥∥∥θk − θ̃k∥∥∥2√E∥∥φ[πM,c/η ∗ ρk](θk)− φ[πM,c/η ∗ ρ̃k](θk)∥∥2.
We bound the second term in the right side of the inequality. E ∥∥φ[πM,c/η ∗ ρk](θk)− φ[πM,c/η ∗ ρ̃k](θk)∥∥2
=E ∥∥∥∥∥∥ 1M M∑ j=1 [ φ[ρk−jc/η](θk)− φ[ρ̃k−jc/η](θk) ]∥∥∥∥∥∥ 2
≤ 1 M M∑ j=1 E ∥∥φ[ρk−jc/η](θk)− φ[ρ̃k−jc/η](θk)∥∥2
= 1
M M∑ j=1 Eθk ∥∥Eθ∼ρk−jc/η φ̄θk(θ)− Eθ∼ρ̃k−jc/η φ̄θk(θ)∥∥2
= 1
M M∑ j=1 Eθk d∑ i=1 ∣∣Eθ∼ρk−jc/η φ̄θk,i(θ)− Eθ∼ρ̃k−jc/η φ̄θk,i(θ)∣∣2 ≤ 1 M M∑ j=1 Eθk d∑ i=1
(∥∥φ̄θk,i(·)∥∥L∞ ∨ ∥∥φ̄θk,i(·)∥∥Lip)2 D2BL [ρk−jc/η, ρ̃k−jc/η] . By Lemma B.4, we have
d∑ i=1 (∥∥φ̄θk,i(·)∥∥L∞ ∨ ∥∥φ̄θk,i(·)∥∥Lip)2 =
d∑ i=1 (∥∥φ̄θk,i(·)∥∥2L∞ ∨ ∥∥φ̄θk,i(·)∥∥2Lip) ≤ [ 4d ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2
BL
+ 4 ∥∥∥e−‖·‖2/σ∥∥∥2
BL ‖∇V (θk)‖2
] .
Plug in the above estimation and by the relation that DBL ≤W1 ≤W2, we have E ‖φ[πM,c ∗ ρk](θk)− φ[πM,c ∗ ρ̃k](θk)‖2
≤ [ 4d ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2
BL
+ 4 ∥∥∥e−‖·‖2/σ∥∥∥2
BL Eθk ‖∇V (θk)‖
2
] 1
M M∑ j=1 D2BL [ρk−cj , ˜ρk−cj ]
≤ [ 4d ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2
BL
+ 4 ∥∥∥e−‖·‖2/σ∥∥∥2
BL Eθk ‖∇V (θk)‖
2
] 1
M M∑ j=1 W22 [ρk−cj , ˜ρk−cj ] .
And combined all the estimation and by the definition of Wasserstein-distance, we conclude that
I2 ≤ √ 4d ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2
BL
+ 4 ∥∥e−‖·‖2/σ∥∥2
BL Eθk ‖∇V (θk)‖
2 √√√√ 1 M M∑ j=1 W22 [ρk−cj , ρ̃k−cj ]
≤ √ 4d ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2
BL
+ 4 ∥∥e−‖·‖2/σ∥∥2
BL Eθk ‖∇V (θk)‖
2 √√√√ 1 M M∑ j=1 E ∥∥∥θk−cj − θ̃k−cj∥∥∥2.
Bounding I3 By Holder’s inequality, I3 ≤ √ E ∥∥∥θk − θ̃k∥∥∥2√E∥∥∥φ[πM,c/η ∗ ρ̃k](θk)− φ[πM,c/η ∗ ρ̃k](θ̃k)∥∥∥2.
We bound the last term on the right side of the inequality. By assumption and Lemma B.3, we have E ∥∥∥φ[πM,c/η ∗ ρ̃k](θk)− φ[πM,c/η ∗ ρ̃k](θ̃k)∥∥∥2
≤ [∥∥∥e−(·)2/σ∥∥∥ Lip Eθ∼ρ̃k ‖∇V (θ)‖+ ∥∥∥∥ 2σ e−(·)2/σ(·) ∥∥∥∥ Lip ]2 E ∥∥∥θk − θ̃k∥∥∥2 .
And combine the estimation, we have
I3 ≤ [∥∥∥e−(·)2/σ∥∥∥ Lip Eθ∼ρ̃k ‖∇V (θ)‖+ ∥∥∥∥ 2σ e−(·)2/σ(·) ∥∥∥∥ Lip ] E ∥∥∥θk − θ̃k∥∥∥2 .
Overall Bound Combine all the results, we have the following bound: for k ≥Mc,
E ∥∥∥θk+1 − θ̃k+1∥∥∥2 − E∥∥∥θk − θ̃k∥∥∥2
≤− 2ηLE ∥∥∥θk − θ̃k∥∥∥2
+2ηα √ E ∥∥∥θk − θ̃k∥∥∥2 c1√
M
+2ηαc2 √√√√ 1 M M∑ j=1 E ∥∥∥θk−jc/η − θ̃k−jc/η∥∥∥2 E∥∥∥θk − θ̃k∥∥∥2
+2ηαc3E ∥∥∥θk − θ̃k∥∥∥2
+η2c4,
where
c1 = √ a2C [ e−rc (1− e−rMc)
1− erc + 1
] ,
c2 =
√ 4d ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2
BL
+ 4 ∥∥e−‖·‖2/σ∥∥2
BL sup k
Eθk ‖∇V (θk)‖ 2 ,
c3 = [∥∥∥e−(·)2/σ∥∥∥ Lip sup k Eθ∼ρ̃k ‖∇V (θ)‖+ ∥∥∥∥ 2σ e−(·)2/σ(·) ∥∥∥∥ Lip ] ,
and
c4 = sup k≥Mc/η
E ∥∥∥∥∥∥∇V (θk) + αφ[ 1M M∑ j=1 δθk−jc/η ](θk)−∇V (θ̃k)− αφ[πM,c/η ∗ ρ̃k](θ̃k) ∥∥∥∥∥∥ 2 .
Define uk = √ E ∥∥∥θk − θ̃k∥∥∥2 and Uk = sup
l∈[k] ul, we have
U2k+1 ≤ qU2k + 2ηαc1√ M Uk + η 2c4,
where q = (1 − 2η(L − αc2 − αc3)). By the assumption that α ≤ L/(c2 + c3), q < 1. Now we prove the bound of Uk by induction. We take the hypothesis that U2k ≤ ( 2ηαc1√ M +(1−q)η(c4+ 11−q ) )2
(1−q)2
and notice that the hypothesis holds for U0 = 0. By the hypothesis, we have
U2k+1 ≤ q
( 2ηαc1√ M + (1− q)η ( c4 + 1 1−q ))2 (1− q)2 + 2ηαc1√ M ( 2ηαc1√ M + (1− q)η ( c4 + 1 1−q )) (1− q) + η2 ( c4 + 1 1− q )
= q
( 2ηαc1√ M + (1− q)η ( c4 + 1 1−q ))2 (1− q)2 + 1 1− q ( 2ηαc1√ M )2 + 2ηαc1√ M η ( c4 + 1 1− q ) + η2 ( c4 + 1 1− q )
= q
( 2ηαc1√ M + (1− q)η ( c4 + 1 1−q ))2 (1− q)2
+ 1− q
(1− q)2 [( 2ηαc1√ M )2 + (1− q)2ηαc1√ M η ( c4 + 1 1− q ) + (1− q)η2 ( c4 + 1 1− q )]
≤ q
( 2ηαc1√ M + (1− q)η ( c4 + 1 1−q ))2 (1− q)2
+ 1− q
(1− q)2 [( 2ηαc1√ M )2 + (1− q)2ηαc1√ M η ( c4 + 1 1− q ) + (1− q)2η2 ( c4 + 1 1− q )2]
= q
( 2ηαc1√ M + (1− q)η ( c4 + 1 1−q ))2 (1− q)2 + (1− q) ( 2ηαc1√ M + (1− q)η ( c4 + 1 1−q ))2 (1− q)2
=
( 2ηαc1√ M + (1− q)η ( c4 + 1 1−q ))2 (1− q)2 ,
where the last second inequality holds by (1− q) ( c4 + 1 1−q ) ≥ 1. Thus we complete the argument of induction and we have, for any k,
U2k ≤
( 2ηαc1√ M + (1− q)η ( c4 + 1 1−q ))2 (1− q)2
≤ 2 4η2α2c21 M + 2(1− q)
2η2 ( c4 + 1 1−q )2 (1− q)2
= 2α2c21 (L− αc2 − αc3)2 1 M + 4η2 (c4 + 2η(L− αc2 − αc3))2 .
And it implies that W22[ρk, ρ̃k] ≤ uk ≤ Uk ≤ 2α2c21 (L−αc2−αc3)2 1 M + 4η 2 (c4 + 2η(L− αc2 − αc3))2 .

heading: B.6 PROOF OF TECHNICAL LEMMAS

heading: B.6.1 PROOF OF LEMMA B.1
text: For the first part:
‖∇V (θ)‖ ≤‖∇V (θ)−∇V (0)‖+ ‖∇V (0)‖ ≤b1 (‖θ1‖+ 1) .
For the second part:
‖θ − η∇V (θ)‖ = 〈θ − η∇V (θ),θ − η∇V (θ)〉 = ‖θ‖2 + 2η 〈θ,−∇V (θ)〉+ η2 ‖∇V (θ)‖2
≤‖θ‖2 + 2η ( −a1 ‖θ‖2 + b1 ) + η2b1(1 + ‖θ‖2)
= ( 1− 2ηa1 + η2b1 ) ‖θ‖2 + η2b1 + 2ηb1.

heading: B.6.2 PROOF OF LEMMA B.2
text: It is obvious that ‖K‖∞,∞ ≤ 1.
‖K(θ′,θ1)−K(θ′,θ2)‖∥∥∥e−‖θ′−θ1‖2/σ − e−‖θ′−θ2‖2/σ∥∥∥ ≤ ∥∥∥e−(·)2/σ∥∥∥
Lip ‖θ1 − θ2‖2 .
And
‖∇θ′K(θ′,θ1)−∇θ′K(θ′,θ2)‖
= ∥∥∥∥ 2σ e−‖θ′−θ1‖2/σ (θ′ − θ1)− 2σ e−‖θ′−θ2‖2/σ (θ′ − θ2) ∥∥∥∥
≤ ∥∥∥∥ 2σ e−(·)2/σ(·) ∥∥∥∥ Lip ‖θ1 − θ2‖2 .

heading: B.6.3 PROOF OF LEMMA B.3
text: For any distribution ρ such that ∫ ‖∇θV (θ)‖ ρ(θ)dθ <∞,
‖φ[ρ](θ1)− φ[ρ](θ2)‖ = ‖Eθ∼ρ {− [K(θ,θ1)−K(θ,θ2)]∇V (θ) +∇1K(θ,θ1)−∇1K(θ,θ2)}‖
≤ ∥∥∥e−(·)2/σ∥∥∥
Lip Eθ∼ρ ‖∇V (θ)‖ ‖θ1 − θ2‖2
+ ∥∥∥∥ 2σ e−(·)2/σ(·) ∥∥∥∥
Lip
‖θ1 − θ2‖2 .
For proving the second result, we notice that
‖φ[ρ](θ)‖ = Eθ′∼ρ [K(θ′,θ)∇V (θ′) +∇1K(θ′,θ)] ≤ ‖K‖∞ Eθ′∼ρ [ ‖∇V (θ′)‖+ 2
σ (‖θ′‖+ ‖θ‖) ] ≤ ‖K‖∞ b1 + Eθ′∼ρ [( 2
σ + b1
) ‖θ′‖+ ‖θ‖ ] .

heading: B.6.4 PROOF OF LEMMA B.4
text: Given any θ′,
d∑ i=1 ∥∥φ̄θ′,i(θ)∥∥2Lip =
d∑ i=1
[ sup
θ1 6=θ2 ∣∣φ̄θ′,i(θ1)− φ̄θ′,i(θ2)∣∣ ‖θ1 − θ2‖2 ]2
= d∑ i=1 sup θ1 6=θ2 ∣∣φ̄θ′,i(θ1)− φ̄θ′,i(θ2)∣∣2 ‖θ1 − θ2‖22 ≤2 d∑ i=1 sup θ1 6=θ2 ∣∣∣(e−‖θ′−θ1‖2/σ − e−‖θ′−θ2‖2/σ) ∂∂θ′iV (θ′)∣∣∣2 ‖θ1 − θ2‖22
+2 d∑ i=1 sup θ1 6=θ2 ∣∣∣ 2σ e−‖θ′−θ1‖2/σ(θ1,i − θ′i)− 2σ e−‖θ′−θ2‖2/σ(θ2,i − θ′i)∣∣∣2 ‖θ1 − θ2‖22 .
For the first term on the right side of the inequality,
d∑ i=1 sup θ1 6=θ2 ∣∣∣(e−‖θ′−θ1‖2/σ − e−‖θ′−θ2‖2/σ) ∂∂θ′iV (θ′)∣∣∣2 ‖θ1 − θ2‖22
= d∑ i=1 ∣∣∣∣ ∂∂θiV (θ′) ∣∣∣∣2 sup θ1 6=θ2 ∣∣∣(e−‖θ′−θ1‖2/σ − e−‖θ′−θ2‖2/σ)∣∣∣2 ‖θ1 − θ2‖22
= ‖∇V (θ′)‖2 ∥∥∥e−‖·‖2/σ∥∥∥2
Lip .
To bound the second term, by the symmetry of each coordinates, we have
d∑ i=1 sup θ1 6=θ2 ∣∣∣ 2σ e−‖θ′−θ1‖2/σ(θ1,i − θ′i)− 2σ e−‖θ′−θ2‖2/σ(θ1,i − θ′i)∣∣∣2 ‖θ1 − θ2‖22
=d ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2
Lip
.
This finishes the first part of the lemma.
d∑ i=d ∥∥φ̄θ′,i(θ)∥∥2L∞ =
d∑ i=d ∥∥∥∥e−‖θ′−θ‖2/σ ( 2σθi − 2σθ′i − ∂∂θ′iV (θ′) )∥∥∥∥2 L∞
≤ d∑ i=d 2 ∥∥∥∥ 2σ e−‖θ′−θ‖2/σ (θi − θ′i) ∥∥∥∥2 L∞ + d∑ i=d 2 ∥∥∥∥e−‖θ′−θ‖2/σ ∂∂θ′iV (θ′) ∥∥∥∥2 L∞
≤ d∑ i=d 2 ∥∥∥∥ 2σ e−‖θ′−θ‖2/σ (θi − θ′i) ∥∥∥∥2 L∞ + 2 ∥∥∥e−‖·‖2/σ∥∥∥2 L∞ ‖∇V (θ′)‖2
≤2d ∥∥∥∥ 2σ e−‖θ‖2/σθ1 ∥∥∥∥2 L∞ + 2 ∥∥∥e−‖·‖2/σ∥∥∥2 L∞ ‖∇V (θ′)‖2 .

Review summaries:
This paper proposed another variant of Langevin dynamics, called “Stein self-repulsive dynamics,” which simultaneously decreases the auto-correlation of Langevin dynamics and eliminates the need for running parallel chains in SVGD. They combined Langevin dynamics with Stein variational gradient descent and theoretically justified that the proposed method successfully converges to the stationary distribution with only a single chain, unlike SVGD. The proposed method decreases the auto-correlation of Langevin dynamics, so the proposed method increases the sample efficiency.