Stein Self-Repulsive Dynamics: Benefits from Past Samples

Mao Ye
UT Austin
[email protected]
&Tongzheng Ren ^*
UT Austin
[email protected]
&Qiang Liu
UT Austin
[email protected]
Equal Contribution

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.

1 Introduction

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; Blei et al., 2017).

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 and Wang, 2016) is a different type of approximate 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 found applications in various challenging problems (e.g., Feng et al., 2017; Haarnoja et al., 2017; Pu et al., 2017; Liu et al., 2017a; Gong et al., 2019). See Han and Liu (e.g., 2018); Chen et al. (e.g., 2018); Liu et al. (e.g., 2019); Wang et al. (e.g., 2019a) for examples of extensions.

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 and 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 thorough theoretical analysis of our new method, establishing its asymptotic convergence to the target distribution. Such result is highly non-trivial, as our new self-repulsive dynamic is a non-linear high-order Markov process. Empirically, we evaluate our methods on an array of challenging sampling tasks, showing that our method yields much better uncertainty estimation and larger effective sample size.

2 Background: Langevin dynamics & SVGD

This section gives a brief introduction on Langevin dynamics (Rossky et al., 1978) and Stein Variational Gradient Descent (SVGD) (Liu and Wang, 2016), which we integrate together to develop our new self-repulsive dynamics for more efficient sampling.

Langevin Dynamics Langevin dynamics is a basic gradient based MCMC method. For a given target distribution supported on $\mathbb{R}^{d}$ with density function $\rho^{*}(\bm{\theta})\propto\exp({-V(\bm{\theta})})$ , where $V\colon\mathbb{R}^{d}\mapsto\mathbb{R}$ is the potential function, the (Euler-discrerized) Langevin dynamics simulates a Markov chain with the following rule:

\bm{\theta}_{k+1}=\bm{\theta}_{k}-\eta\nabla V(\bm{\theta}_{k})+\sqrt{2\eta}\bm{e}_{k},~{}~{}~{}~{}~{}~{}\bm{e}_{k}\sim\mathcal{N}(0,\mathbf{I}),

where $k$ denotes the number of iterations, $\{\bm{e}_{k}\}_{k=1}^{\infty}$ are independent standard Gaussian noise, and $\eta$ is a step size parameter. It is well known that the limiting distribution of $\bm{\theta}_{k}$ when $k\to\infty$ approximates the target distribution when $\eta$ is sufficiently small.

Because the updates in Langevin dynamics are local and incremental, new points generated by the dynamics can be highly correlated to the past sample, in which case 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 $\{\bm{\theta}_{0}^{i}\}_{i=1}^{M}$ , SVGD updates the $M$ particles in parallel by

\bm{\theta}_{k+1}^{i}=\bm{\theta}_{k}^{i}+\eta{\bm{g}}(\bm{\theta}_{k}^{i};~{}\hat{\delta}_{k}^{M}),~{}~{}~{}~{}\forall i=1,\ldots,M,

where ${\bm{g}}(\bm{\theta}_{k}^{i};~{}\hat{\delta}_{k}^{M})$ is a velocity field that depends the empirical distribution of the current set of particles $\hat{\delta}_{k}^{M}:=\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k}^{j}}$ in the following way,

\displaystyle{\bm{g}}(\bm{\theta}_{k}^{i};\hat{\delta}_{k}^{M})=\mathbb{E}_{\bm{\theta}\sim\hat{\delta}_{k}^{M}}\!\!\bigg{[}\underset{\mathrm{Confining\ Term}}{\text{$\underbrace{-K(\bm{\theta},\bm{\theta}_{k}^{i})\nabla V(\bm{\theta})}$}}+\underset{\mathrm{Repulsive\ Term}}{\underbrace{\nabla_{\bm{\theta}}K(\bm{\theta},\bm{\theta}_{k}^{i})}}\!\bigg{]}\!\!.

Here $\delta_{\bm{\theta}}$ is the Dirac measure centered at $\bm{\theta}$ , and hence $\mathbb{E}_{\bm{\theta}\sim\hat{\delta}_{k}^{M}}[\cdot]$ denotes averaging on the particles. The $K(\cdot,\cdot)$ is a positive definite kernel, such as the RBF kernel, that can be specified by users.

Note that ${\bm{g}}(\bm{\theta}_{k}^{i};~{}\hat{\delta}_{k}^{M})$ 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 $\rho^{*}(\bm{\theta})\propto\exp(-V(\bm{\theta}))$ at the fixed point, when the number of particles goes to infinite. We refer the readers to Liu and Wang (2016); Liu (2017); Liu and Wang (2018) for thorough 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 ${\bm{g}}(\bm{\theta};~{}\rho)$ equals zero when $\rho$ equals the true distribution $\rho^{*}$ , that is, $\forall\bm{\theta}^{\prime}\in\mathbb{R}^{d}$ ,

\displaystyle{\bm{g}}(\bm{\theta}^{\prime};\rho^{*})=\mathbb{E}_{\bm{\theta}\sim\rho^{*}}\left[-K(\bm{\theta},\bm{\theta}^{\prime})\nabla V(\bm{\theta})+\nabla_{\bm{\theta}}K(\bm{\theta},\bm{\theta}^{\prime})\right]=0.

(1)

This equation shows that, the target distributions forms a fixed point of the update, and SVGD would converge if the particle distribution $\hat{\delta}_{k}^{M}$ gives a close approximation to the target distribution $\rho^{*}$ .

3 Stein Self-Repulsive Dynamics

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 particle and the particles from previous iterations, hence forming a new self-avoiding dynamics with fast convergence speed.

Specifically, assume we run a single Markov chain like Langevin dynamics, where $\bm{\theta}_{k}$ denotes the sample at the $k$ -th iteration. Denote by $\tilde{\delta}_{k}^{M}$ the empirical measure of $M$ samples from the past iterations:

\small\tilde{\delta}_{k}^{M}:=\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc_{\eta}}},~{}~{}~{}~{}~{}~{}~{}~{}c_{\eta}=c/\eta,

where $c_{\eta}$ is a thinning factor, which scales inversely with the step size $\eta$ , introduced to slim the sequence of past samples. Compared with the $\hat{\delta}_{k}^{M}$ in SVGD, which is averaged over $M$ parallel particles, $\tilde{\delta}_{k}^{M}$ is averaged across time over $M$ past samples. Given this, our Stein self-repulsive dynamics updates the sample via

\displaystyle\bm{\theta}_{k+1}\leftarrow\bm{\theta}_{k}~{}+~{}\underbrace{(-\eta V(\bm{\theta}_{k})+\sqrt{2\eta}\bm{e}_{k})}_{\text{Langevin}}~{}+~{}\underbrace{\eta\alpha{\bm{g}}(\bm{\theta}_{k};~{}\tilde{\delta}_{k}^{M})}_{\text{Stein Repulsive}},

(2)

in which the particle is updated with the typical Langevin gradient, plus a Stein repulsive force against the particles from the previous iterations. Here $\alpha\geq 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 if the self-repulsive dynamics has converged to the target distribution $\rho^{*}$ , such that $\theta_{k}\sim\rho^{*}$ 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 $\tilde{\delta}_{k}^{M}$ 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,

\displaystyle\bm{\theta}_{k+1}=\begin{cases}\bm{\theta}_{k}\!-\!\eta\nabla V(\bm{\theta}_{k})\!+\!\sqrt{2\eta}\bm{e}_{k},&k<Mc_{\eta},\\ \bm{\theta}_{k}\!+\!\eta\left[-\nabla V(\bm{\theta}_{k})\!+\!\alpha{\bm{g}}(\bm{\theta}_{k};\tilde{\delta}_{k}^{M})\right]\!+\!\sqrt{2\eta}\bm{e}_{k},&k\geq Mc_{\eta}.\end{cases}

(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 Note 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 and 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 $\sum_{j}D(\theta_{k},\theta_{k-jc_{\eta}})$ , where $D(\cdot,\cdot)$ can be any kind of dis-similarity. However, introducing an arbitrary repulsive force would alter the stationary distribution of the dynamics, introducing a harmful bias into the algorithm. Besides, the self-adjusting mechanism in Deng et al. (2020b) can also be viewed as a repulsive force using the multiplier in gradient. 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 Recent works (Gallego and Insua, 2018; Zhang et al., 2018) also combine SVGD with Langevin dynamics, in which, however, the Langevin force is directly added to a set of particles that evolve in parallel with SVGD. Using our terminology, their system is

\displaystyle\bm{\theta}_{k+1}^{i}=\bm{\theta}_{k}^{i}+{(-\eta V(\bm{\theta}_{k}^{i})+\sqrt{2\eta}\bm{e}_{k}^{i})}+{\eta\alpha{\bm{g}}(\bm{\theta}_{k}^{i};~{}\hat{\delta}_{k}^{M})},\quad\bm{e}_{k}\sim\mathcal{N}(0,\mathbf{I}),\quad\forall i=1,\ldots,M.

This is significantly different from our method on both motivation and practical algorithm. Their algorithm still requires to simulate $M$ parallel chains of particles 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.

Refer to caption — Figure 1: Illustrating the advantage of our Self-Repulsive Langevin dynamics. With a set of initial examples locating on the left part of the target distribution (yellow dots), Self-Repulsive Langevin dynamics is forced to explore the right part more frequently, yielding an accurate approximation when combined with the initial samples. Langevin dynamics, however, does not take the past samples into account and yields a poor overall approximation.

An Illustrative Example 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 Gaussian 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 biased overall representation of the true distribution (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.

4 Theoretical Analysis of Stein Self-Repulsive Dynamics

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 particle size $M$ approaches to infinite and the step size $\eta$ approaches to $0$ . This is a highly non-trivial task, as 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:

(1)

At the limit of $M\to\infty$ (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).
(2)

By further taking the limit of $\eta\to 0^{+}$ (called the continuous-time limit), the dynamics in (4) converges to a continuous-time mean-field dynamics characterized by (5).
(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 $\left\|\cdot\right\|$ and $\left\langle\cdot,\cdot\right\rangle$ to represent the $\ell_{2}$ vector norm and inner product, respectively. The Lipschitz norm and bounded Lipschitz norm of a function $f$ are defined by $\left\|f\right\|_{\mathrm{Lip}}$ and $\left\|f\right\|_{\mathrm{BL}}$ . The KL divergence, Wasserstein-2 distance and Bounded Lipschitz distance between distribution $\rho_{1}$ , $\rho_{2}$ are denoted as $\mathbb{D}_{\mathrm{KL}}[\rho_{1}\|\rho_{2}]$ , $\mathbb{W}_{2}[\rho_{1},\rho_{2}]$ and $\mathbb{D}_{\mathrm{BL}}[\rho_{1},\rho_{2}]$ , respectively.

4.1 Mean-Field and Continuous-Time Limits

To fix the notation, we denote by $\rho_{k}:=\mathrm{Law}(\bm{\theta}_{k})$ the distribution of $\bm{\theta}_{k}$ at time $k$ of the practical self-repulsive dynamics (3), which we refer to the practical dynamics in the sequel, when the initial particle $\bm{\theta}_{0}$ is drawn from an initial continuous distribution $\rho_{0}$ . Note that given $\rho_{0}$ , the subsequent $\rho_{k}$ can be recursively defined through dynamics (3). Due to the diffusion noise in Langevin dynamics, all $\rho_{k}$ are continuous distributions supported on $\mathbb{R}^{d}$ . We now introduce the limit dynamics when we take the mean-field limit $(M\to+\infty)$ and then the continuous-time limit $(\eta\to 0^{+})$ .

Discrete-Time Mean-Field Dynamics ( $M\to+\infty$ ) In the limit of $M\to\infty$ , our practical dynamics (3) approaches to the following limit dynamics, in which the delta measures on the particles are replaced by the actual continuous distributions of the particles,

\displaystyle\tilde{\bm{\theta}}_{k+1}=\begin{cases}\tilde{\bm{\theta}}_{k}\!-\!\eta\nabla V(\tilde{\bm{\theta}}_{k})\!+\!\sqrt{2\eta}\bm{e}_{k},&k\leq Mc_{\eta},\\ \tilde{\bm{\theta}}_{k}\!+\!\eta\left[-\nabla V(\tilde{\bm{\theta}}_{k})\!+\!\alpha{\bm{g}}(\tilde{\bm{\theta}}_{k},\tilde{\rho}_{k}^{M})\right]\!+\!\sqrt{2\eta}\bm{e}_{k},&k\geq Mc_{\eta}.\end{cases}

(4)

where $\tilde{\rho}_{k}^{M}=\frac{1}{M}\sum_{j=1}^{M}\tilde{\rho}_{k-jc_{\eta}}$ and $\tilde{\rho}_{k}:=\mathrm{Law}(\tilde{\bm{\theta}}_{k})$ is the (smooth) distribution of $\tilde{\bm{\theta}}_{k}$ at time-step $k$ when the dynamics is initialized with $\tilde{\bm{\theta}}_{0}\sim\tilde{\rho}_{0}=\rho_{0}$ . Compared with the practical dynamics in (3), the difference is that the empirical distribution $\tilde{\delta}_{k}^{M}$ is replaced by the smooth distribution $\tilde{\rho}_{k}^{M}$ . Similar to the recursive definition of $\rho_{k}$ following dynamics (3), $\tilde{\rho}_{k}$ is also recursively defined through dynamics (4), starting from $\tilde{\rho}_{0}=\rho_{0}$ .

As we show in Theorem 4.3, if the auto-correlation of $\bm{\theta}_{k}$ decays fast enough and $M$ is sufficiently large, $\tilde{\rho}_{k}^{M}$ is well approximated by the empirical distribution $\tilde{\delta}_{k}^{M}$ in (3), and further the two dynamics ((3) and (4)) converges to each other in the sense that $\mathbb{W}_{2}[\rho_{k},\tilde{\rho}_{k}]\to 0$ as $M\to\infty$ for any $k$ . Note that in taking the limit of $M\to\infty$ , we need to ensure that we run the dynamics for more than $Mc_{\eta}$ steps. Otherwise, SRLD degenerates to Langevin dynamics as we stop the chain before we finish collecting the $M$ samples.

Continuous-Time Mean-Field Dynamics ( $\eta\to 0^{+}$ ) In the limit of zero step size $(\eta\to 0^{+})$ , the discrete-time mean field dynamics in (4) can be shown to converge to the following continuous-time mean-field dynamics:

\displaystyle d\bar{\bm{\theta}}_{t}

\displaystyle=\begin{cases}-\nabla V(\bar{\bm{\theta}}_{t})dt+d\mathcal{B}_{t},&t\in[0,Mc),\\ \left[-\nabla V(\bar{\bm{\theta}}_{t})+\alpha{\bm{g}}(\bar{\bm{\theta}}_{t},~{}\bar{\rho}_{t}^{M})\right]dt+d\mathcal{B}_{t},&t\geq Mc.\end{cases}

(5)

where $\bar{\rho}_{t}^{M}:=\frac{1}{M}\sum_{j=1}^{M}\bar{\rho}_{t-jc}(\cdot)$ , $\mathcal{B}_{t}$ is the Brownian motion and $\bar{\rho}_{t}=\mathrm{Law}\left(\bar{\bm{\theta}}_{t}\right)$ is the distribution of $\bar{\bm{\theta}}_{t}$ at a continuous time point $t$ with $\bm{\theta}_{0}$ initialized by $\bar{\bm{\theta}}_{0}\sim\tilde{\rho}_{0}=\rho_{0}$ . We prove that (5) is closely approximated by (4) with small step size in the sense that $\mathbb{D}_{\mathrm{KL}}[\tilde{\rho}_{k}\|\bar{\rho}_{k\eta}]\to 0$ as $\eta\to 0$ in Theorem 4.2, and importantly, the stationary distribution of (5) equals to the target distribution $\rho^{*}(\bm{\theta})\propto\exp(-V(\bm{\theta})).$

4.2 Assumptions

We first introduce the techinical assumptions used in our theoretical analysis.

Assumption 4.1 (RBF Kernel).

We use RBF kernel, i.e. $K(\bm{\theta}_{1},\bm{\theta}_{2})=\exp({-\left\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\|^{2}/\sigma})$ , for some fixed $0<\sigma<\infty$ .

We only assume the RBF kernel for the simplicity of our analysis. However, it is straightforward to generalize our theoretical result to other positive definite kernels.

Assumption 4.2 ( $V$ is dissipative and smooth).

Assume that $\left\langle\bm{\theta},-\nabla V(\bm{\theta})\right\rangle\leq b_{1}-a_{1}\left\|\bm{\theta}\right\|^{2}$ and $\left\|\nabla V(\bm{\theta}_{1})-\nabla V(\bm{\theta}_{2})\right\|\leq b_{1}\left\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\|$ . We also assume that $\left\|\nabla V(\bm{0})\right\|\leq b_{1}$ . Here $a_{1}$ and $b_{1}$ are some finite positive constant.

Assumption 4.3 (Regularity Condition).

Assume $\mathbb{E}_{\bm{\theta}\sim\rho_{0}}[\left\|\bm{\theta}\right\|^{2}]>0$ . Define $\rho_{k}^{M}=\sum_{j=1}^{M}\rho_{k-jc_{\eta}}/M$ , assume there exists $a_{2},B<\infty$ such that

\sup_{k\geq Mc_{\eta}}\frac{\mathbb{E}\left\|g(\bm{\theta}_{k};\tilde{\delta}_{k}^{M})-g(\bm{\theta}_{k};\rho_{k}^{M})\right\|^{2}}{\sup_{\left\|\bm{\theta}\right\|\leq B}\mathbb{E}\left\|g(\bm{\theta};\tilde{\delta}_{k}^{M})-g(\bm{\theta};\rho_{k}^{M})\right\|^{2}}\leq a_{2}.

Assumption 4.4 (Strong-convexity).

Suppose that $\left\langle\nabla V(\bm{\theta}_{1})-\nabla V(\bm{\theta}_{2}),\bm{\theta}_{1}-\bm{\theta}_{2}\right\rangle\geq L\left\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\|^{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; Deng et al., 2020a). 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.

4.3 Main Theorems

All of the proofs in this section can be found in Appendix E. 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:

\partial_{t}\bar{\rho}_{t}=\begin{cases}\nabla\cdot\left(-\nabla V\bar{\rho}_{t}\right)+\Delta\bar{\rho}_{t}&t\in[0,Mc)\\ \nabla\cdot\left[\left(-\nabla V+\alpha{\bm{g}}(\cdot,\bar{\rho}_{t}^{M})\right)\bar{\rho}_{t}\right]+\Delta\bar{\rho}_{t},&t\geq Mc.\end{cases}

Theorem 4.1 (Stationary Distribution).

Given some finite $M$ , $c$ and $\alpha$ , and suppose that the limit distribution of (5) exists. Then the limit distribution is unique and satisfies $\rho^{*}(\bm{\theta})\propto\exp({-V(\bm{\theta})})$ .

We then give the upper bound on the discretization error, which can be characterized by analyzing the KL divergence between $\tilde{\rho}_{k}$ and $\bar{\rho}_{k\eta}$ .

Theorem 4.2 (Time Discretization Error).

Given some sufficiently small step size $\eta$ and choose $\alpha<a_{2}/(2b_{1}+4/\sigma)$ . Under Assumption 4.1, 4.2, 4.3 and $c_{\eta}=c/\eta$ . we have for some constant $C$ ,

\displaystyle\underset{l\in\{0,...,k\}}{\max}\ \mathbb{D}_{\mathrm{KL}}\left[\bar{\rho}_{l\eta}\|\tilde{\rho}_{l}\right]\leq\begin{cases}\mathcal{O}\left(\eta+k\eta^{2}\right)&k\leq Mc_{\eta}-1\\ \mathcal{O}\left(\eta+Mc\eta+\alpha^{2}Mce^{C\alpha^{2}(k\eta-Mc)}\eta^{2}\right)&k\geq Mc_{\eta}.\end{cases}

With this theorem, we can know that if $\eta$ is small enough, then the discretization error is small and $\tilde{\rho}$ approximates $\bar{\rho}$ closely. Next we give result on the mean field limit of $M\to\infty$ .

Theorem 4.3 (Mean-Field Limit).

Under Assumption 4.1, 4.2, 4.3, and 4.4, suppose that we choose $\alpha$ and $\eta$ such that $-(a_{1}-2\alpha b_{1}/\sigma)+\eta b_{1}<0$ ; $\frac{2\alpha\eta}{\sigma}(b_{1}+1)<1$ ; $a_{2}-\alpha\left(2b_{1}+\frac{4}{\sigma}\right)>0$ ; Then there exists a constant $c_{2}$ , such that when $L/a\geq c_{2}$ and we have

\displaystyle\mathbb{W}_{2}^{2}[\rho_{k},\tilde{\rho}_{k}]=\begin{cases}\mathcal{O}\left({\alpha^{2}/{M}}+\eta^{2}\right)&\geq Mc_{\eta},\\ 0&k\leq Mc_{\eta}-1.\end{cases}

Thus, if $M$ is sufficiently large, $\rho_{k}$ can well approximate the $\tilde{\rho}_{k}$ . Combining all the above theorems, we have the following Corollary showing the convergence of the proposed practical algorithm to the target distribution.

Corollary 4.1 (Convergence to Target Distribution).

Under the assumptions of Theorem 4.1, 4.2 and 4.3, by choosing $k,\eta,M$ such that $k\eta\to\infty$ , $\exp(C\alpha^{2}k\eta)\eta^{2}=o(1)$ and $\frac{k\eta}{Mc}=\gamma\left(1+o(1)\right)$ with $\gamma>1$ , we have

\lim_{k,M\to\infty,\eta\to 0^{+}}\mathbb{D}_{\text{BL}}\left[\rho_{k},\rho^{*}\right]=0.

Remark A careful choice of parameters is needed as our system is a complicated longitudinal particle system. Also notice that if $\gamma\leq 1$ , the repulsive dynamics reduces to Langevin dynamics, as only the samples from the first phase will be collected.

5 Extension to General Dynamics

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(\bm{\theta})\propto\exp(-V(\bm{\theta}))$ :

	$\displaystyle d\bm{\theta}_{t}=-\bm{f}(\bm{\theta})dt+\sqrt{2\bm{D}(\bm{\theta})}d\mathcal{B}_{t},$
	$\displaystyle\text{with}\quad\bm{f}(\bm{\theta})=[\bm{D}(\bm{\theta})+\bm{Q}(\bm{\theta})]\nabla V(\bm{\theta})-\bm{\Gamma}(\bm{\theta}),\quad\Gamma_{i}(\bm{\theta})=\sum_{j=1}^{d}\frac{\partial}{\partial\text{$\bm{\theta}$}_{j}}\left(D_{ij}(\bm{\theta})+Q_{ij}(\bm{\theta})\right).$

where $\bm{D}$ is a positive semi-definite diffusion matrix that determines the strength of the Brownian motion and $\bm{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

\displaystyle d\bar{\bm{\theta}}_{t}

\displaystyle=\begin{cases}-\bm{f}(\bm{\theta})dt+\sqrt{2\bm{D}(\bm{\theta})}d\mathcal{B}_{t},&t\in[0,Mc)\\ -\left(\bm{f}(\bm{\theta})+\alpha{\bm{g}}(\bar{\bm{\theta}}_{t};~{}\bar{\rho}_{t}^{M})\right)dt+d\mathcal{B}_{t},&t\geq Mc\end{cases}

(6)

where $\bar{\rho}_{t}:=\mathrm{Law}(\bar{\bm{\theta}}_{t})$ is again the distribution of $\bar{\bm{\theta}}_{t}$ following (6) when initalized at $\bar{\bm{\theta}}_{0}\sim\bar{\rho}_{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 finite $M$ , $c$ and $\alpha$ , and suppose that the limiting distribution of dynamics (6) exists. Then the limiting distribution is unique and equals the target distribution $\rho^{*}(\bm{\theta})\propto\exp({-V(\bm{\theta})})$ .

6 Experiments

In this section, we evaluate the proposed method in various challenging tasks. We demonstrate the effectiveness of SRLD in high dimensions by applying it to sample the posterior of Bayesian Neural Networks. To demonstrate the superiority of the SRLD in obtaining diversified samples, we apply SRLD on contextual bandits problem, which requires the sampler efficiently explores the distribution in order to give good uncertainty estimation.

We include discussion on the parameter tuning and additional experiment on sampling high dimensional Gaussian and Gaussian mixture in Appendix B. Our code is available at https://github.com/lushleaf/Stein-Repulsive-Dynamics.

6.1 Synthetic Experiment

We first show how the repulsive gradient helps explore the whole distribution using a synthetic distribution that is easy to visualize. Following Ma et al. (2015), we compare the sampling efficiency on the following correlated 2D distribution with density

\rho^{*}([\theta_{1},\theta_{2}])\propto-\theta_{1}^{4}/10-\left(4\left(\theta_{2}+1.2\right)-\theta_{1}^{2}\right)^{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 2: Sample quality of SRLD and Langevin dynamics for sampling the correlated 2D distribution. The auto-correlation is the averaged auto-correlation of the two dimensions.

Figure 3: Sampling trajectory of the correlated 2D distribution.

Figure 2 summarizes the result with different metrics. We can see that SRLD has 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 3 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.

6.2 Bayesian Neural Network

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 and Graff, 2017). We assume the output is normal distributed, with a two-layer neural network with 50 hidden units and $\mathrm{tanh}$ activation to predict the mean of outputs. All of the datasets are randomly partitioned into $90\%$ for training and $10\%$ for testing. The results are averaged over 20 random trials. We refer readers to Appendix C for hyper-parameter tuning and other experiment details.

Dataset	Ave Test RMSE			Ave Test LL
Dataset	SVGD	LD	SRLD	SVGD	LD	SRLD
Boston	$3.300\pm 0.142$	$3.342\pm 0.187$	$\underline{{\bf 3.086\pm 0.181}}$	$-4.276\pm 0.217$	$-2.678\pm 0.092$	$\underline{{\bf-2.500\pm 0.054}}$
Concrete	$4.994\pm 0.171$	$4.908\pm 0.113$	${\bf 4.886\pm 0.108}$	$-5.500\pm 0.398$	$-3.055\pm 0.035$	$\underline{{\bf-3.034\pm 0.031}}$
Energy	$0.428\pm 0.016$	$0.412\pm 0.016$	$\underline{{\bf 0.395\pm 0.016}}$	$-0.781\pm 0.094$	$-0.543\pm 0.014$	$\underline{{\bf-0.476\pm 0.036}}$
Naval	$0.006\pm 0.000$	$0.006\pm 0.002$	$\underline{{\bf 0.003\pm 0.000}}$	$3.056\pm 0.034$	$4.041\pm 0.030$	$\underline{{\bf 4.186\pm 0.015}}$
WineRed	$0.655\pm 0.008$	$0.649\pm 0.009$	$\underline{{\bf 0.639\pm 0.009}}$	$-1.040\pm 0.018$	$-1.004\pm 0.019$	$\underline{{\bf-0.970\pm 0.016}}$
WineWhite	$\underline{{\bf 0.655\pm 0.008}}$	$0.692\pm 0.003$	$0.688\pm 0.003$	$\underline{{\bf-1.040\pm 0.019}}$	$-1.047\pm 0.004$	$-1.043\pm 0.004$
Yacht	$0.593\pm 0.071$	$0.597\pm 0.051$	${\bf 0.578\pm 0.054}$	$-1.281\pm 0.279$	$-1.187\pm 0.307$	$\underline{{\bf-0.458\pm 0.036}}$

Table 1: Averaged test RMSE and test log-likelihood on UCI datasets. Results are averaged over 20 trails. The boldface indicates the method has the best average performance and the underline marks the methods that perform the best with a significance level of

0.05

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 C that SRLD averaged on 20 particles (across time) can also achieve similar or better results as SVGD with 20 (parallel) particles.

6.3 Contextual Bandits

We consider the posterior sampling (a.k.a Thompson sampling) algorithm with Bayesian neural network as the function approximator, to demonstrate the uncertanty estimation provided by SRLD. We follow 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 compare the SRLD with the Langevin dynamics on two benchmarks from (Riquelme et al., 2018), and include SVGD as a baseline. For more detailed introduction, setup, hyper-parameter tuning and experiment details; see Appendix D.

Dataset	SVGD	LD	SRLD
Mushroom	$20.7\pm 2.0$	$4.28\pm 0.09$	$\underline{\mathbf{3.80\pm 0.16}}$
Wheel	$91.32\pm 0.17$	$38.07\pm 1.11$	$\underline{\mathbf{32.08\pm 0.75}}$

Table 2: Cumulative Regrets on two bandits problem (smaller is better). Results are averaged over 10 trails. Boldface indicates the methods with best performance and underline marks the best significant methods with significant level

0.05

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.

7 Conclusion

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.

Broader Impact Statement

This work incorporates Stein repulsive force into Langevin dynamics to improve sample efficiency. It brings a positively improvement to the community such as reinforcement learning and Bayesian neural network that needs efficient sampler. Our work do not have any negative societal impacts that we can foresee in the future.

Acknowledgement

This paper is supported in part by NSF CAREER 1846421, SenSE 2037267 and EAGER 2041327.

References

An and Huang (1996) HZ An and FC Huang. The geometrical ergodicity of nonlinear autoregressive models. Statistica Sinica, pages 943–956, 1996.
Blei et al. (2017) David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 112(518):859–877, 2017.
Blundell et al. (2015) Charles Blundell, Julien Cornebise, Koray Kavukcuoglu, and Daan Wierstra. Weight uncertainty in neural network. In Francis Bach and David Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 1613–1622, Lille, France, 07–09 Jul 2015. PMLR. URL http://proceedings.mlr.press/v37/blundell15.html.
Bradley (2005) Richard Bradley. Basic properties of strong mixing conditions. a survey and some open questions. Probability Surveys, 2:107–144, 2005.
Bradley (2007) Richard C. Bradley. Introduction to strong mixing conditions. 2007.
Brooks et al. (2011) Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng. Handbook of Markov chain Monte Carlo. CRC press, 2011.
Bubeck et al. (2012) Sébastien Bubeck, Nicolo Cesa-Bianchi, et al. Regret analysis of stochastic and nonstochastic multi-armed bandit problems. Foundations and Trends® in Machine Learning, 5(1):1–122, 2012.
Chapelle and Li (2011) Olivier Chapelle and Lihong Li. An empirical evaluation of Thompson sampling. In Advances in neural information processing systems, pages 2249–2257, 2011.
Chen et al. (2018) Changyou Chen, Ruiyi Zhang, Wenlin Wang, Bai Li, and Liqun Chen. A unified particle-optimization framework for scalable Bayesian sampling. In Proceedings of the Thirty-Fourth Conference on Uncertainty in Artificial Intelligence, UAI 2018, Monterey, California, USA, August 6-10, 2018, pages 746–755, 2018.
Dalalyan (2017) Arnak S Dalalyan. Theoretical guarantees for approximate sampling from smooth and log-concave densities. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(3):651–676, 2017.
Deng et al. (2020a) Wei Deng, Qi Feng, Liyao Gao, Faming Liang, and Guang Lin. Non-convex learning via replica exchange stochastic gradient mcmc. In International Conference on Machine Learning, pages 2474–2483. PMLR, 2020a.
Deng et al. (2020b) Wei Deng, Guang Lin, and Faming Liang. A contour stochastic gradient langevin dynamics algorithm for simulations of multi-modal distributions, 2020b.
Ding et al. (2014) Nan Ding, Youhan Fang, Ryan Babbush, Changyou Chen, Robert D. Skeel, and Hartmut Neven. Bayesian sampling using stochastic gradient thermostats. In Advances in Neural Information Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-13 2014, Montreal, Quebec, Canada, pages 3203–3211, 2014.
Dua and Graff (2017) Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL http://archive.ics.uci.edu/ml.
Durmus et al. (2019) Alain Durmus, Szymon Majewski, and Blazej Miasojedow. Analysis of langevin monte carlo via convex optimization. Journal of Machine Learning Research, 20(73):1–46, 2019.
Feng et al. (2017) Yihao Feng, Dilin Wang, and Qiang Liu. Learning to draw samples with amortized stein variational gradient descent. uncertainty in artificial intelligence, 2017.
Gallego and Insua (2018) Victor Gallego and David Rios Insua. Stochastic gradient mcmc with repulsive forces. arXiv preprint arXiv:1812.00071, 2018.
Gong et al. (2019) Chengyue Gong, Jian Peng, and Qiang Liu. Quantile Stein variational gradient descent for batch Bayesian optimization. In International Conference on Machine Learning, pages 2347–2356, 2019.
Gretton et al. (2012) Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. Journal of Machine Learning Research, 13(Mar):723–773, 2012.
Haarnoja et al. (2017) Tuomas Haarnoja, Haoran Tang, Pieter Abbeel, and Sergey Levine. Reinforcement learning with deep energy-based policies. international conference on machine learning, pages 1352–1361, 2017.
Han and Liu (2017) Jun Han and Qiang Liu. Stein variational adaptive importance sampling. Uncertainty in Artificial Intelligence, 2017.
Han and Liu (2018) Jun Han and Qiang Liu. Stein variational gradient descent without gradient. In Uncertainty in Artificial Intelligence, 2018.
Laio and Parrinello (2002) Alessandro Laio and Michele Parrinello. Escaping free-energy minima. Proceedings of the National Academy of Sciences, 99(20):12562–12566, 2002.
Liu et al. (2019) Chang Liu, Jingwei Zhuo, Pengyu Cheng, Ruiyi Zhang, and Jun Zhu. Understanding and accelerating particle-based variational inference. In International Conference on Machine Learning, pages 4082–4092, 2019.
Liu (2017) Qiang Liu. Stein variational gradient descent as gradient flow. In Advances in neural information processing systems, pages 3118–3126, 2017.
Liu and Wang (2016) Qiang Liu and Dilin Wang. Stein variational gradient descent: A general purpose bayesian inference algorithm. In Advances In Neural Information Processing Systems, pages 2378–2386, 2016.
Liu and Wang (2018) Qiang Liu and Dilin Wang. Stein variational gradient descent as moment matching. In Advances in Neural Information Processing Systems, pages 8854–8863, 2018.
Liu et al. (2017a) Yang Liu, Qiang Liu, and Jian Peng. Stein variational policy gradient. uncertainty in artificial intelligence, 2017a.
Liu et al. (2017b) Yang Liu, Prajit Ramachandran, Qiang Liu, and Jian Peng. Stein variational policy gradient. Uncertainty in Artificial Intelligence, 2017b.
Ma et al. (2015) Yi-An Ma, Tianqi Chen, and Emily Fox. A complete recipe for stochastic gradient mcmc. In Advances in Neural Information Processing Systems, pages 2917–2925, 2015.
Neal et al. (2011) Radford M Neal et al. Mcmc using hamiltonian dynamics. Handbook of markov chain monte carlo, 2(11):2, 2011.
Pu et al. (2017) Yunchen Pu, Zhe Gan, Ricardo Henao, Chunyuan Li, Shaobo Han, and Lawrence Carin. Vae learning via stein variational gradient descent. neural information processing systems, pages 4236–4245, 2017.
Raginsky et al. (2017) Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky. Non-convex learning via stochastic gradient langevin dynamics: a nonasymptotic analysis. In Proceedings of the 30th Conference on Learning Theory, COLT 2017, Amsterdam, The Netherlands, 7-10 July 2017, pages 1674–1703, 2017. URL http://proceedings.mlr.press/v65/raginsky17a.html.
Riquelme et al. (2018) Carlos Riquelme, George Tucker, and Jasper Snoek. Deep bayesian bandits showdown: An empirical comparison of bayesian deep networks for thompson sampling. In 6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings, 2018. URL https://openreview.net/forum?id=SyYe6k-CW.
Rossky et al. (1978) Peter J Rossky, JD Doll, and HL Friedman. Brownian dynamics as smart monte carlo simulation. The Journal of Chemical Physics, 69(10):4628–4633, 1978.
Schlimmer (1981) Jeff Schlimmer. Mushroom records drawn from the audubon society field guide to north american mushrooms. GH Lincoff (Pres), New York, 1981.
Thompson (1933) William R. Thompson. On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika, 25(3/4):285–294, 1933.
Wainwright et al. (2008) Martin J Wainwright, Michael I Jordan, et al. Graphical models, exponential families, and variational inference. Foundations and Trends® in Machine Learning, 1(1–2):1–305, 2008.
Wang et al. (2019a) Dilin Wang, Ziyang Tang, Chandrajit Bajaj, and Qiang Liu. Stein variational gradient descent with matrix-valued kernels. In Conference on Neural Information Processing Systems, 2019a.
Wang et al. (2019b) Ziyu Wang, Tongzheng Ren, Jun Zhu, and Bo Zhang. Function space particle optimization for bayesian neural networks. arXiv preprint arXiv:1902.09754, 2019b.
Zhang et al. (2018) Jianyi Zhang, Ruiyi Zhang, and Changyou Chen. Stochastic particle-optimization sampling and the non-asymptotic convergence theory. arXiv preprint arXiv:1809.01293, 2018.

Appendix A Discussion on Hyper-parameter Tuning

The key hyper-parameters of SRLD are 1. $\alpha$ , which balance the confining gradient and repulsive gradient; 2. $M$ the number of particles used; 3. $\sigma$ the bandwidth of kernel; 4. $\eta$ the stepsize; 5. $c_{\eta}$ the thinning factor. Among which, $\alpha$ , $M$ and $\sigma$ are the hyper-parameter introduced by the proposed repulsive gradient and thus we mainly discuss these three hyper-parameter. The number of particles $M$ and bandwidth of kennel $\sigma$ are introduced by the use the repulsive term in SVGD [Liu and Wang, 2016]. In practice, we find using a similar setting for tuning $M$ and $\sigma$ as that in SVGD [Liu and Wang, 2016] gives good performance. In specific, in order to obtain good performance, $M$ does not needs be very large, and similar to SVGD, $M=10$ already gives good enough particle approximation. A good choice of bandwidth $\sigma$ is also important for the kernel. In SVGD, instead of tuning $\sigma$ , they propose an adaptive way to adjust $\sigma$ during the running of the dynamics. Specifically, they choose $\sigma=\text{med}^{2}/\log(M)$ , where med is the median of the pairwise distance between the particles $\bm{\theta}_{k}^{i}$ , $i\in[M]$ . In this way, the bandwidth ensures that $\sum_{j=1}^{M}K(\bm{\theta}_{k}^{i},\bm{\theta}_{k}^{j})\approx 1$ . This adaptive way of choosing $\sigma$ is also widely used in current approximation inference area, e.g., Liu et al. [2017b], Han and Liu [2017], Wang et al. [2019b, a]. We also find that applying this adaptive bandwidth is able to give good empirical performance and thus we also use this method in the implementation. Now we discuss how choose $\alpha$ . Notice that $\alpha$ serves to balance the confining gradient and repulsive gradient and based on this motivation, we recommend readers to find a proper $\alpha$ using the samples at burn-in phase by setting

\alpha\approx\frac{\sum_{k=1}^{Mc_{\eta}}\left\|\nabla V(\bm{\theta}_{k})\right\|}{\sum_{k=1}^{Mc_{\eta}}\left\|g(\bm{\theta}_{k},\tilde{\delta}_{k}^{M})\right\|}.

In this way, $\alpha$ balances the two kind of gradients. And then we may further tune $\alpha$ by searching around this value. An empirical good choice of $\alpha$ is 10 for the data sets we tested and we use $\alpha=10$ for all the experiments.

The step size is important for gradient based MCMC, as too large step size gives too large discretization error while a too small step size will cause the dynamics converges very slowly. In this paper, we mainly use validation set to tune the step size. The thinning factor is also a common parameter is MCMC methods and usually MCMC methods are not sensitive to this parameter. SRLD is not sensitive to this parameter and we simply set $c_{\eta}=100$ for all experiments.

Appendix B Additional Experiment Result on Synthetic Data

In this section, we show additional experiment on synthetic data. To further visualize the role of the proposed stein repulsive gradient, we also apply our method to sample a 2D mixture of Gaussian distribution (see section B.1). To further study how different $\alpha$ influences sampling high dimension distribution, we apply SRLD to sample high dimensional Gaussian (section B.2) and high dimensional mixture of Gaussian (section B.3).

B.1 Synthetic 2D Mixture of Gaussian Experiment

Figure 4: Sample quality and autocorrelation of the mixture distribution. The auto-correlation is the averaged auto-correlation of the two dimensions.

Figure 5: Sampling trajectory of the mixture of Gaussian.

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

\displaystyle\rho^{*}(\bm{\theta})\propto 0.5\exp{\left(-\left\|\bm{\theta}-\bm{1}\right\|^{2}/2\right)}+0.5\exp{\left(-\left\|\bm{\theta}+\bm{1}\right\|^{2}/2\right)},

where $\bm{\theta}=[\theta_{1},\theta_{2}]^{\top}$ and $\mathbf{1}=[1,1]^{\top}$ . This target distribution have two mode at $-\bm{1}$ and $\bm{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 4 shows that SRLD consistently outperforms the Langevin dynamics on all of the evaluation metrics.

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.

B.2 Synthetic higher dimensional Gaussian Experiment

To show the performance of SRLD in higher dimensional case with different value of $\alpha$ , we additionally considering the problem on sampling from Gaussian distribution with $d=100$ and covariance $\bm{\Sigma}=0.5\textbf{I}$ . We run SRLD with $\alpha=100,50,20,10,0$ and the case $\alpha=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 6. In this experiment, we set one SRLD with an inappropriate $\alpha=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 $\alpha$ induces too much extra approximation error and thus its performance is not as good as these with smaller $\alpha$ (see MMD and Wasserstein distance). This phenomenon matches our theoretical finding.

Figure 6: Sample quality and autocorrelation of the higher dimensional Gaussian distribution. The auto-correlation is the averaged auto-correlation of all dimensions.

B.3 Synthetic higher dimensional Mixture of Gaussian Experiment

We also consider sampling from the mixture of Gaussian with $d=20$ . The target density is set to be

\rho^{*}(\bm{\theta})\propto\frac{1}{2}\exp\left(-0.5\left\|\bm{\theta}-\sqrt{2/d}\bm{1}\right\|^{2}\right)+\frac{1}{2}\exp\left(-0.5\left\|\bm{\theta}+\sqrt{2/d}\bm{1}\right\|^{2}\right),

where $\bm{\theta}=[\theta_{1},...,\theta_{20}]^{\top}$ and $\bm{1}=[1,...,1]^{\top}$ . And thus the mean of the two mixture component is with distance $2\sqrt{2}$ . We run SRLD with $\alpha=20,10,5,0$ (and when $\alpha=0$ , it reduces to LD). The other experiment setting is the same as the low dimensional mixture Gaussian case. Figure [7] summarizes the result. As shown in the figure, when $\alpha$ becomes larger, the repulsive forces helps the sampler better explore the density region.

Figure 7: Sample quality and autocorrelation of the higher dimensional mixture distribution. The auto-correlation is the averaged auto-correlation of all dimensions.

Appendix C BNN on UCI Datasets: Experiment Settings and Additional Results

We first give detailed experiment settings. We set a $\Gamma(1,0.1)$ prior for the inverse output variance. We set the mini-batch size to be 100. 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_{\eta}=c/\eta=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. 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 and Wang [2016]. We use $\alpha=10$ for all the data sets. For SVGD, we use the original implementation with 20 particles by Liu and Wang [2016].

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 8 and 9 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.

Dataset	Ave Test RMSE
Dataset	SRLD-20	SRLD-40	SRLD-60	SRLD-80	SRLD-100	SVGD
Boston	$\bf{3.236\pm 0.174}$	$\bf{3.173\pm 0.176}$	$\bf{3.130\pm 0.173}$	$\underline{\bf{3.101\pm 0.179}}$	$\underline{{\bf 3.086\pm 0.181}}$	$3.300\pm 0.142$
Concrete	${\bf 4.959\pm 0.109}$	${\bf 4.921\pm 0.111}$	${\bf 4.906\pm 0.109}$	${\bf 4.891\pm 0.108}$	${\bf 4.886\pm 0.108}$	$4.994\pm 0.171$
Energy	${\bf 0.422\pm 0.016}$	${\bf 0.409\pm 0.016}$	$\underline{{\bf 0.405\pm 0.016}}$	$\underline{{\bf 0.399\pm 0.016}}$	$\underline{{\bf 0.395\pm 0.016}}$	$0.428\pm 0.016$
Naval	${\bf 0.005\pm 0.001}$	$\underline{{\bf 0.004\pm 0.000}}$	$\underline{{\bf 0.003\pm 0.000}}$	$\underline{{\bf 0.003\pm 0.000}}$	$\underline{{\bf 0.003\pm 0.000}}$	$0.006\pm 0.000$
WineRed	${\bf 0.654\pm 0.009}$	$\underline{{\bf\mathbf{0.647\pm 0.009}}}$	$\underline{{\bf 0.644\pm 0.009}}$	$\underline{{\bf 0.641\pm 0.009}}$	$\underline{{\bf 0.639\pm 0.009}}$	$0.655\pm 0.008$
WineWhite	$0.695\pm 0.003$	$0.692\pm 0.003$	$0.690\pm 0.003$	$0.689\pm 0.002$	$0.688\pm 0.003$	$\underline{{\bf 0.655\pm 0.008}}$
Yacht	$0.616\pm 0.055$	$0.608\pm 0.052$	$0.597\pm 0.051$	${\bf 0.587\pm 0.054}$	${\bf 0.578\pm 0.054}$	$0.593\pm 0.071$

Table 3: Comparing SRLD with different number of samples with SVGD on test RMSE. The results are computed over 20 trials. For SRLD, the value is set to be boldface if it has better average performance than SVGD. The value if with underline if it is significantly better than SVGD with significant level 0.05 using a matched pair t-test.

Dataset	Ave Test LL
Dataset	SRLD-20	SRLD-40	SRLD-60	SRLD-80	SRLD-100	SVGD
Boston	$\underline{{\bf-2.642\pm.088}}$	$\underline{{\bf-2.582\pm 0.084}}$	$\underline{{\bf-2.527\pm 0.612}}$	$\underline{{\bf-2.516\pm 0.062}}$	$\underline{{\bf-2.500\pm 0.054}}$	$-4.276\pm 0.217$
Concrete	$\underline{{\bf-3.084\pm 0.036}}$	$\underline{{\bf-3.061\pm 0.034}}$	$\underline{{\bf-3.050\pm 0.033}}$	$\underline{{\bf-3.040\pm 0.031}}$	$\underline{{\bf-3.034\pm 0.031}}$	$-5.500\pm 0.398$
Energy	$\underline{{\bf-0.580\pm 0.053}}$	$\underline{{\bf-0.536\pm 0.048}}$	$\underline{{\bf-0.522\pm 0.046}}$	$\underline{{\bf-0.504\pm 0.044}}$	$\underline{{\bf-0.476\pm 0.036}}$	$-0.781\pm 0.094$
Naval	$\underline{{\bf 4.033\pm 0.230}}$	$\underline{{\bf 4.100\pm 0.171}}$	$\underline{{\bf 4.140\pm 0.015}}$	$\underline{{\bf 4.167\pm 0.014}}$	$\underline{{\bf 4.186\pm 0.015}}$	$3.056\pm 0.034$
WineRed	$\underline{{\bf-1.008\pm 0.019}}$	$\underline{{\bf-0.990\pm 0.017}}$	$\underline{{\bf-0.982\pm 0.016}}$	$\underline{{\bf-0.974\pm 0.016}}$	$\underline{{\bf-0.970\pm 0.016}}$	$-1.040\pm 0.018$
WineWhite	$-1.053\pm 0.004$	$-1.049\pm 0.004$	$-1.047\pm 0.004$	$-1.044\pm 0.004$	$-1.043\pm 0.004$	$\underline{{\bf-1.040\pm 0.019}}$
Yacht	$\underline{{\bf-1.160\pm 0.256}}$	$\underline{{\bf-0.650\pm 0.173}}$	$\underline{{\bf-0.556\pm 0.096}}$	$\underline{{\bf-0.465\pm 0.037}}$	$\underline{{\bf-0.458\pm 0.036}}$	$-1.281\pm 0.279$

Table 4: Comparing SRLD with different number of samples with SVGD on test log-likelihood. The results are computed over 20 trials. For SRLD, the value is set to be boldface if it has better average performance than SVGD. The value if with underline if it is significantly better than SVGD with significant level 0.05 using a matched pair t-test.

Figure 8: Comparison between SRLD and Langevin dynamics on test RMSE. The results are computed based on 20 repeats. The error bar is calculated based on RMSE of SRLD - RMSE of Langevin dynamics in 20 repeats to rule out the variance of different data splitting

Figure 9: Comparison between SRLD and Langevin dynamics on test log-likelihood. The results are computed based on 20 repeats. The error bar is calculated based on log-likelihood of SRLD - log-likelihood of Langevin dynamics in 20 repeats to rule out the variance of data splitting.

Appendix D Contextual Bandit: Experiment Settings and More Background

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 $\mathcal{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,\cdots,N$ , where $N$ is some pre-defined time horizon (and can be given to the agent), the environment provides a context $s_{t}\in\mathcal{S}$ to the agent, then the agent should choose one action $a_{t}\in\{1,2,\cdots,K\}$ based on context $s_{t}$ . The environment will return a (stochastic) reward $r(s_{t},a_{t})$ to the agent based on the context $s_{t}$ and the action $a_{t}$ 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:

\overline{R}_{N}^{\mathcal{S}}=\max_{\pi:\mathcal{S}\to\{1,2,\cdots,K\}}\mathbb{E}\left[\sum_{t=1}^{N}r(s_{t},g(s_{t}))-\sum_{t=1}^{N}r(s_{t},a_{t})\right].

(7)

where $\pi$ denotes the deterministic mapping from the context set $\mathcal{S}$ to actions $\{1,2,\cdots,K\}$ (readers can view $\pi$ as a deterministic policy in reinforcement learning). Intuitively, this pseudo-regret measures the difference of cumulative reward between the action sequence $a_{t}$ and the best action sequence $\pi(s_{t})$ . Thus, an algorithm that can minimize the pseudo-regret (7) can also find the best $\pi$ .

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 and Li, 2011]. It works by first placing an user-specified prior $\mu^{0}_{s,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 $\mu^{t}_{s,a}$ with the observation $r(s_{t-1},a_{t-1})$ at time $t-1$ where $a_{t-1}$ is selected with the posterior distribution: each time, the action is selected with the following way:

a_{t}=\mathop{\arg\max}_{a\in\{1,2,\cdots,K\}}{\hat{r}(s_{t},a)},\quad\hat{r}(s_{t},a)\sim\mu^{t}_{s,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

\mu^{0}_{s,a}

, time horizon

N

for time

t=1,2,\cdots,N

observe a new context

s_{t}\in\mathcal{S}

sample the reward of each action

\hat{r}(s_{t},a)\sim\mu^{t}_{s,a}

a\in\{1,2,\cdots,K\}

select action

a_{t}=\mathop{\arg\max}_{a\in\{1,2,\cdots,K\}}{\hat{r}(s_{t},a)}

and get the reward

r(s_{t},a_{t})

update the posterior

\mu^{t+1}_{s_{t},a_{t}}

with

r(s_{t},a_{t})

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.

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.

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” $\delta\in(0,1)$ . The context set $\mathcal{S}$ is the unit circle $\|s\|_{2}\leq 1$ in $\mathbb{R}^{2}$ , and each turn the context $s_{t}$ is uniformly sampled from $\mathcal{S}$ . $K=5$ possible actions are provided: the first action yields a constant reward $r\sim\mathcal{N}(\mu_{1},\sigma^{2})$ ; the reward corresponding to other actions is determined by the provided context $s$ :

•

For $s\in\mathcal{S}$ s.t. $\|s\|_{2}\leq\delta$ , all of the four other actions return a suboptimal reward sampled from $\mathcal{N}(\mu_{2},\sigma^{2})$ for $\mu_{2}<\mu_{1}$ .
•

For $s\in\mathcal{S}$ s.t. $\|s\|_{2}>\delta$ , according to the quarter the context $s$ is in, one of the four actions becomes optimal. This optimal action gives a reward of $\mathcal{N}(\mu_{3},\sigma^{2})$ for $\mu_{3}\gg\mu_{1}$ , and another three actions still yield the suboptimal reward $\mathcal{N}(\mu_{2},\sigma^{2})$ .

Following the setting from [Riquelme et al., 2018], we set $\mu_{1}=1.2$ , $\mu_{2}=1.0$ , and $\mu_{3}=50$ .

When $\delta$ approaches $1$ , the inner circle $\|s\|_{2}\leq\delta$ 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 $a_{1}$ 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 $\delta=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/\eta=100$ for vanilla Langevin dynamics and SRLD, and set $M=20$ , $\alpha=10$ 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 replace the warm-up of stochastic variational inference methods in Riquelme et al. [2018] with the burn-in phase of the sampling. Similar to other methods in [Riquelme et al., 2018], we keep the initial learning rate as $10^{-1}$ for fast burn-in and the step-size for sampling is tuned on the mushroom bandit and keep the same for both the mushroom and wheel bandit. 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 use 20 particles to make the comparison fair, and also tune the step-size on the mushroom bandit.

Appendix E The Detailed analysis of SRLD

E.1 Some additional notation

We use $\|\cdot\|_{\infty}$ to denote the $\ell_{\infty}$ vector norm and define the $\mathcal{L}_{\infty}$ norm of a function $f:\mathbb{R}^{d}\to\mathbb{R}^{1}$ as $\left\|f\right\|_{\mathcal{L}_{\infty}}$ . $\mathbb{D}_{\mathrm{TV}}$ denote the Total Variation distance between distribution $\rho_{1},\rho_{2}$ respectively. Also, as $K$ is $\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{1}$ , we denote $\left\|K\right\|_{\mathcal{L}_{\infty},\mathcal{L}_{\infty}}=\sup_{\bm{\theta}_{1},\bm{\theta}_{2}}K(\bm{\theta}_{1},\bm{\theta}_{2})$ . For simplicity, we may use $\left\|K\right\|_{\infty,\infty}$ as $\left\|K\right\|_{\mathcal{L}_{\infty},\mathcal{L}_{\infty}}$ . In the appendix, we also use $\phi[\rho](\bm{\theta})\coloneqq\bm{g}(\bm{\theta};\rho)$ , where $\bm{g}(\bm{\theta};\rho)$ is defined in the main text. For the clearance, we define $\pi_{M,c/\eta}*\rho_{k}\coloneqq\rho_{k}^{M}$ , $\pi_{M,c/\eta}*\tilde{\rho}_{k}\coloneqq\tilde{\rho}_{k}^{M}$ and $\pi_{M,c}*\bar{\rho}_{t}\coloneqq\bar{\rho}_{t}^{M}$ , where $\rho_{k}^{M}$ , $\tilde{\rho}_{k}^{M}$ and $\bar{\rho}_{t}^{M}$ are defined in main text.

E.2 Geometric Ergodicity of SRLD

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/\eta+1)$ -order nonlinear autoregressive model when $k\geq Mc/\eta$ :

\bm{\theta}_{k+1}=\bm{\psi}\left(\bm{\theta}_{k},...,\bm{\theta}_{k-Mc/\eta}\right)+\sqrt{2\eta}\bm{e}_{k},

where

\bm{\psi}\left(\bm{\theta}_{k},...,\bm{\theta}_{k-Mc/\eta}\right)=\bm{\theta}_{k}+\eta\left\{-\nabla V(\bm{\theta}_{k})+\alpha\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta}_{k})\right\}.

Further, if we stack the parameter by $\bm{\Theta}_{k}=\left[\bm{\theta}_{k},...,\bm{\theta}_{k-Mc/\eta}\right]^{\top}$ and define $\bm{\Psi}\left(\bm{\Theta}_{k}\right)=\left[\bm{\psi}^{\top}\left(\bm{\Theta}_{k}\right),\bm{\Theta}_{k}^{\top}\right]^{\top}$ , we have

\bm{\Theta}_{k+1}=\bm{\Psi}\left(\bm{\Theta}_{k}\right)+\sqrt{2\eta}\bm{E}_{k},

where $\bm{E}_{k}=\left[\bm{e}_{k}^{\top},\bm{0}^{\top},...,\bm{0}^{\top}\right]^{\top}.$ In this way, we formulate $\bm{\Theta}_{k}$ as a time homogeneous Markov Chain. In the following analysis, we only analyze the second phase of SRLD given some initial stacked particles $\bm{\Theta}_{Mc/\eta-1}$ .

Theorem E.1 (Geometric Ergodicity).

Under Assumption 4.1 and Assumption 4.2, suppose we choose $\eta$ and $\alpha$ such that

\max\left(1-2\eta a_{1}+\eta^{2}b_{1}+\frac{2\alpha\eta}{\sigma}b_{1},\frac{2\alpha\eta}{\sigma}(b_{1}+1)\right)<1,

then the Markov Chain of $\bm{\Theta}_{k}$ is stationary, geometrically ergodic, i.e., for any $\bm{\Theta}^{\prime}_{0}=\bm{\Theta}_{Mc/\eta-1}$ , we have

\mathbb{D}_{\mathrm{TV}}\left[P^{k}\left(\cdot,\bm{\Theta}_{0}\right),\Pi\left(\cdot\right)\right]\leq Q\left(\bm{\Theta}_{0}\right)e^{-rk},

where $r=\mathcal{O}(\eta)$ is some positive constant, $Q(\bm{\Theta}_{0})$ is constant related to $\bm{\Theta}_{0}$ , $P^{k}$ is the $k$ -step Markov transition kernel and $\Pi$ is the stationary distribution.

We defer the proof to Appendix E.5.1.

E.3 Moment Bound

Theorem E.2 (Moment Bound).

Under Assumption 4.2, suppose that we have $\mathbb{E}_{\bm{\theta}\sim\rho_{0}}\left\|\bm{\theta}\right\|^{2}<\infty$ ; and $a_{2}-\alpha\left\|K\right\|_{\infty}\left(2b_{1}+\frac{4}{\sigma}\right)>0$ , we have

		$\displaystyle\sup_{k}\mathbb{E}_{\bm{\theta}\sim\rho_{k}}\left\\|\bm{\theta}\right\\|^{2}\lor\sup_{k}\ \mathbb{E}_{\bm{\theta}\sim\tilde{\rho}_{k}}\left\\|\bm{\theta}\right\\|^{2}\lor\sup_{t}\mathbb{E}_{\bm{\theta}\sim\bar{\rho}_{t}}\left\\|\bm{\theta}\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle\mathbb{E}_{\bm{\theta}\sim\rho_{0}}\left\\|\bm{\theta}\right\\|^{2}+\frac{b_{1}+1+\eta}{a_{2}-\left\\|K\right\\|_{\mathcal{L}_{\infty},\mathcal{L}_{\infty}}\frac{2\alpha}{\sigma}-\alpha\left\\|K\right\\|_{\mathcal{L}_{\infty},\mathcal{L}_{\infty}}\left(2b_{1}+\frac{2}{\sigma}\right)}.$

And by Lemma E.1, we thus have

		$\displaystyle\sup_{k}\ \mathbb{E}_{\bm{\theta}\sim\rho_{k}}\left\\|\nabla V(\bm{\theta})\right\\|^{2}\lor\sup_{k}\ \mathbb{E}_{\bm{\theta}\sim\tilde{\rho}_{k}}\left\\|\nabla V(\bm{\theta})\right\\|^{2}\lor\sup_{t}\mathbb{E}_{\bm{\theta}\sim\bar{\rho}_{t}}\left\\|\nabla V(\bm{\theta})\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle b_{1}\mathbb{E}_{\bm{\theta}\sim\rho_{0}}\left\\|\bm{\theta}\right\\|^{2}+\frac{b_{1}(b_{1}+1+\eta)}{a_{2}-\left\\|K\right\\|_{\mathcal{L}_{\infty},\mathcal{L}_{\infty}}\frac{2\alpha}{\sigma}-\alpha\left\\|K\right\\|_{\mathcal{L}_{\infty},\mathcal{L}_{\infty}}\left(2b_{1}+\frac{2}{\sigma}\right)}+1$

The proof can be found at Appendix E.5.2.

E.4 Technical Lemma

Definition E.1 ( $\alpha$ -mixing).

For any two $\sigma$ -algebras $\mathcal{A}$ and $\mathcal{B}$ , the $\alpha$ -mixing coefficient is defined by

\alpha(\mathcal{A},\mathcal{B})=\underset{A\in\mathcal{A},B\in\mathcal{B}}{\sup}\left|\mathbb{P}\left(A\cap B\right)-\mathbb{P}\left(A\right)\mathbb{P}\left(B\right)\right|.

Let $\left(X_{k},k\geq 1\right)$ be a sequence of real random variable defined on $\left(\Omega,\mathcal{A},\mathbb{P}\right)$ . This sequence is $\alpha$ -mixing if

\alpha(n)\coloneqq\underset{k\geq 1}{\sup}\ \alpha\left(\mathcal{M}_{k},\mathcal{G}_{k+n}\right)\to 0,\ \mathrm{as}\ n\to\infty,

where $\mathcal{M}_{j}\coloneqq\sigma\left(X_{i},i\leq j\right)$ and $\mathcal{G}_{j}\coloneqq\sigma\left(X_{i},i\geq j\right)$ for $j\geq 1$ . Alternatively, as shown by Theorem 4.4 of Bradley [2007]

\alpha(n)\coloneqq\frac{1}{4}\sup\left\{\frac{\mathrm{Cov}\left(f,g\right)}{\left\|f\right\|_{\mathcal{L}_{\infty}}\left\|g\right\|_{\mathcal{L}_{\infty}}};\ f\in\mathcal{L}_{\infty}\left(\mathcal{M}_{k}\right),\ g\in\mathcal{L}_{\infty}\left(\mathcal{G}_{k+n}\right)\right\}.

Definition E.2 ( $\beta$ -mixing).

For any two $\sigma$ -algebras $\mathcal{A}$ and $\mathcal{B}$ , the $\alpha$ -mixing coefficient is defined by

\beta\left(\mathcal{A},\mathcal{B}\right)\coloneqq\sup\frac{1}{2}\sum_{i=1}^{I}\sum_{j=1}^{J}\left|\mathbb{P}\left(A_{i}\cap B_{j}\right)-\mathbb{P}\left(A_{i}\right)\mathbb{P}\left(B_{j}\right)\right|,

where the supremum is taken over all pairs of finite partitions $\left\{A_{1},...,A_{I}\right\}$ and $\left\{B_{1},...,B_{J}\right\}$ of $\Omega$ such that $A_{i}\in\mathcal{A}$ and $B_{j}\in\mathcal{B}$ for each $i$ , $j$ . Let $\left(X_{k},k\geq 1\right)$ be a sequence of real random variable defined on $\left(\Omega,\mathcal{A},\mathbb{P}\right)$ . This sequence is $\beta$ -mixing if

\beta(n)\coloneqq\underset{k\geq 1}{\sup}\ \beta\left(\mathcal{M}_{k},\mathcal{G}_{k+n}\right)\to 0,\ \mathrm{as}\ n\to\infty.

Proposition E.1 ( $\beta$ -mixing implies $\alpha$ -mixing)).

For any two $\sigma$ -algebras $\mathcal{A}$ and $\mathcal{B}$ ,

\alpha\left(\mathcal{A},\mathcal{B}\right)\leq\frac{1}{2}\beta\left(\mathcal{A},\mathcal{B}\right).

This proposition can be found in Equation 1.11 of Bradley [2005].

Proposition E.2.

A (strictly) stationary Markov Chain is geometric ergodicity if and only if $\beta(n)\to 0$ at least exponentially fast as $n\to\infty$ .

This proposition is Theorem 3.7 of Bradley [2005].

Lemma E.1 (Regularity Conditions).

By Assumption 4.2, we have $\left\|\nabla V(\bm{\theta})\right\|\leq b_{1}\left(\left\|\bm{\theta}_{1}\right\|+1\right)$ and $\left\|\bm{\theta}-\eta\nabla V(\bm{\theta})\right\|\leq\left(1-2\eta a_{1}+\eta^{2}b_{1}\right)\left\|\bm{\theta}\right\|^{2}+\eta^{2}b_{1}+2\eta b_{1}.$

Lemma E.2 (Properties of RBF Kernel).

For RBF kernel with bandwidth $\sigma$ , we have $\left\|K\right\|_{\infty,\infty}\leq 1$ and

	$\displaystyle\left\\|K(\bm{\theta}^{\prime},\bm{\theta}_{1})-K(\bm{\theta}^{\prime},\bm{\theta}_{2})\right\\|$	$\displaystyle\leq$	$\displaystyle\left\\|e^{-(\cdot)^{2}/\sigma}\right\\|_{\mathrm{Lip}}\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}$
	$\displaystyle\left\\|\nabla_{\bm{\theta}^{\prime}}K(\bm{\theta}^{\prime},\bm{\theta}_{1})-\nabla_{\bm{\theta}^{\prime}}K(\bm{\theta}^{\prime},\bm{\theta}_{2})\right\\|$	$\displaystyle\leq$	$\displaystyle\left\\|\frac{2}{\sigma}e^{-(\cdot)^{2}/\sigma}(\cdot)\right\\|_{\mathrm{Lip}}\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}.$

Lemma E.3 (Properties of Stein Operator).

For any distribution $\rho$ such that $\mathbb{E}_{\bm{\theta}\sim\rho}\left\|\nabla V(\bm{\theta})\right\|<\infty,$ we have

	$\displaystyle\left\\|\phi[\rho](\cdot)\right\\|_{\mathrm{Lip}}$	$\displaystyle\leq\left\\|e^{-(\cdot)^{2}/\sigma}\right\\|_{\mathrm{Lip}}\mathbb{E}_{\bm{\theta}\sim\rho}\left\\|\nabla V(\bm{\theta})\right\\|+\left\\|\frac{2}{\sigma}e^{-(\cdot)^{2}/\sigma}(\cdot)\right\\|_{\mathrm{Lip}},$
	$\displaystyle\left\\|\phi[\rho](\bm{\theta})\right\\|$	$\displaystyle\leq\left\\|K\right\\|_{\infty}\mathbb{E}_{\bm{\theta}^{\prime}\sim\rho}\left[\left\\|\nabla V(\bm{\theta}^{\prime})\right\\|+\frac{2}{\sigma}\left(\left\\|\bm{\theta}^{\prime}\right\\|+\left\\|\bm{\theta}\right\\|\right)\right]$
		$\displaystyle\leq\left\\|K\right\\|_{\infty}b_{1}+\mathbb{E}_{\bm{\theta}^{\prime}\sim\rho}\left[\left(\frac{2}{\sigma}+b_{1}\right)\left\\|\bm{\theta}^{\prime}\right\\|\right]+\left\\|\bm{\theta}\right\\|.$

Lemma E.4 (Bounded Lipschitz of Stein Operator).

Given $\bm{\theta}^{\prime}$ , define $\bar{\phi}_{\bm{\theta}^{\prime}}(\bm{\theta})\coloneqq\phi[\delta_{\bm{\theta}^{\prime}}](\bm{\theta})=K(\bm{\theta}^{\prime},\bm{\theta})\nabla V(\bm{\theta}^{\prime})+\nabla_{1}K(\bm{\theta}^{\prime},\bm{\theta})$ . We also denote $\bar{\phi}_{\bm{\theta}^{\prime}}(\bm{\theta})=[\bar{\phi}_{\bm{\theta}^{\prime},1}(\bm{\theta}),...,\bar{\phi}_{\bm{\theta}^{\prime},d}(\bm{\theta})]^{\top}$ . We have

	$\displaystyle\sum_{i=1}^{d}\left\\|\bar{\phi}_{\bm{\theta}^{\prime},i}(\bm{\theta})\right\\|_{\mathrm{Lip}}^{2}$	$\displaystyle\leq$	$\displaystyle 2\left\\|\nabla V(\bm{\theta}^{\prime})\right\\|^{2}\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathrm{Lip}}^{2}+2d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathrm{Lip}}^{2}$
	$\displaystyle\sum_{i=d}^{d}\left\\|\bar{\phi}_{\bm{\theta}^{\prime},i}(\bm{\theta})\right\\|_{\mathcal{L}_{\infty}}^{2}$	$\displaystyle\leq$	$\displaystyle 2d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathcal{L}_{\infty}}^{2}+2\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathcal{L}_{\infty}}^{2}\left\\|\nabla V(\bm{\theta}^{\prime})\right\\|^{2}.$

E.5 Proof of Main Theorems

E.5.1 Proof of Theorem E.1

The proof of this theorem is by verifying the condition of Theorem 3.2 of An and Huang [1996]. Suppose $\bm{\Theta}=\left[\bm{\theta}_{1},...,\bm{\theta}_{MC+1}\right]$ , where $C=c/\eta$ , we have

	$\displaystyle\left\\|\bm{\psi}\left(\bm{\Theta}\right)\right\\|$	$\displaystyle=\left\\|\bm{\theta}_{1}+\eta\left\{-\nabla V(\bm{\theta}_{1})+\alpha\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{1+jC}}](\bm{\theta}_{k})\right\}\right\\|$
		$\displaystyle=\left\\|\bm{\theta}_{1}-\eta\nabla V(\bm{\theta}_{1})+\frac{\eta\alpha}{M}\sum_{j=1}^{M}\left[e^{-\left\\|\bm{\theta}_{1+jC}-\bm{\theta}_{1}\right\\|^{2}/\sigma}\frac{2}{\sigma}\left(\bm{\theta}_{1}-\bm{\theta}_{1+jC}\right)-e^{-\left\\|\bm{\theta}_{1+jC}-\bm{\theta}_{1}\right\\|^{2}/\sigma}\nabla V(\bm{\theta}_{1+jC})\right]\right\\|$
		$\displaystyle\leq\left\\|\bm{\theta}_{1}-\eta\nabla V(\bm{\theta}_{1})+\frac{2}{\sigma}\frac{\eta\alpha}{M}\sum_{j=1}^{M}e^{-\left\\|\bm{\theta}_{1+jC}-\bm{\theta}_{1}\right\\|^{2}/\sigma}\bm{\theta}_{1}\right\\|$
		$\displaystyle+\left\\|\frac{\eta\alpha}{M}\sum_{j=1}^{M}e^{-\left\\|\bm{\theta}_{1+jC}-\bm{\theta}_{1}\right\\|^{2}/\sigma}\frac{2}{\sigma}\left(-\nabla V(\bm{\theta}_{1+jC})-\bm{\theta}_{1+jC}\right)\right\\|$
		$\displaystyle\leq\left\\|\bm{\theta}_{1}-\eta\nabla V(\bm{\theta}_{1})\right\\|+\frac{2\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty}b_{1}(1+\left\\|\bm{\theta}_{1}\right\\|)$
		$\displaystyle+\frac{2\alpha\eta}{M\sigma}\sum_{j=1}^{M}\left\\|K\right\\|_{\infty,\infty}b_{1}\left(1+(1+\frac{1}{b_{1}})\left\\|\bm{\theta}_{1+jC}\right\\|\right)$
		$\displaystyle\overset{(1)}{\leq}b_{1}(1+\frac{4\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty})+\eta^{2}b_{1}+2\eta b_{1}$
		$\displaystyle+\left(1-2\eta a_{1}+\eta^{2}b_{1}+\frac{2\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty}b_{1}\right)\left\\|\bm{\theta}_{1}\right\\|+\frac{2\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty}(b_{1}+1)\max_{i\in[MC+1]-\{1\}}\left\\|\bm{\theta}_{1+jC}\right\\|$
		$\displaystyle\leq b_{1}(1+\frac{4\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty})+\eta^{2}b_{1}+2\eta b_{1}$
		$\displaystyle+\max\left(1-2\eta a_{1}+\eta^{2}b_{1}+\frac{2\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty}b_{1},\frac{2\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty}(b_{1}+1)\right)\max_{i\in[MC+1]}\left\\|\bm{\theta}_{1+jC}\right\\|,$

where $(1)$ is by Lemma E.1. Thus, given the step size $\eta$ , if we choose $\eta$ , $\alpha$ such that

\max\left(1-2\eta a_{1}+\eta^{2}b_{1}+\frac{2\alpha\eta}{\sigma}\left\|K\right\|_{\infty,\infty}b_{1},\frac{2\alpha\eta}{\sigma}\left\|K\right\|_{\infty,\infty}(b_{1}+1)\right)<1,

then our dynamics is geometric ergodic.

E.5.2 Proof of Theorem E.2

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 $\alpha=0$ , here we only analysis the second phase. Define $U_{t}=\underset{s\leq t}{\sup}\sqrt{\mathbb{E}\left\|\bar{\bm{\theta}}_{s}\right\|^{2}}$ , and thus

\frac{\partial}{\partial t}U_{t}^{2}\leq\mathbb{E}\left\langle\bar{\bm{\theta}}_{t},-V(\bar{\bm{\theta}})+\alpha\phi[\pi_{M,c}*\bar{\rho}_{t}](\bar{\bm{\theta}}_{t})\right\rangle\vee 0.

Now we bound $\mathbb{E}\left\langle\bar{\bm{\theta}}_{t},-V(\bar{\bm{\theta}})+\alpha\phi[\pi_{M,c}*\bar{\rho}_{t}](\bar{\bm{\theta}}_{t})\right\rangle$ :

		$\displaystyle\mathbb{E}\left\langle\bar{\bm{\theta}}_{t},-V(\bar{\bm{\theta}}_{t})+\alpha\phi[\pi_{M,c}*\bar{\rho}_{t}](\bar{\bm{\theta}}_{t})\right\rangle$
	$\displaystyle\leq$	$\displaystyle b_{1}-a_{2}\mathbb{E}\left\\|\bar{\bm{\theta}}_{t}\right\\|^{2}+\alpha\mathbb{E}\left\\|\bar{\bm{\theta}}_{t}\right\\|\left\\|\phi[\pi_{M,c}*\bar{\rho}_{t}](\bar{\bm{\theta}}_{t})\right\\|$
	$\displaystyle\stackrel{{\scriptstyle(1)}}{{\leq}}$	$\displaystyle b_{1}-a_{2}\mathbb{E}\left\\|\bar{\bm{\theta}}_{t}\right\\|^{2}+\alpha\left\\|K\right\\|_{\infty}\mathbb{E}\left\\|\bar{\bm{\theta}}_{t}\right\\|\mathbb{E}_{\bm{\theta}^{\prime}\sim\pi_{M,c}*\bar{\rho}_{t}}\left[\left\\|\nabla V(\bm{\theta}^{\prime})\right\\|+\frac{2}{\sigma}\left(\left\\|\bm{\theta}^{\prime}\right\\|+\left\\|\bm{\theta}_{t}\right\\|\right)\right]$
	$\displaystyle\leq$	$\displaystyle b_{1}-a_{2}\mathbb{E}\left\\|\bar{\bm{\theta}}_{t}\right\\|^{2}+\alpha\left\\|K\right\\|_{\infty}\mathbb{E}\left\\|\bar{\bm{\theta}}_{t}\right\\|\mathbb{E}_{\bm{\theta}^{\prime}\sim\pi_{M,c}*\bar{\rho}_{t}}\left[b_{1}\left(\left\\|\bm{\theta}^{\prime}\right\\|+1\right)+\frac{2}{\sigma}\left(\left\\|\bm{\theta}^{\prime}\right\\|+\left\\|\bm{\theta}_{t}\right\\|\right)\right]$
	$\displaystyle=$	$\displaystyle b_{1}-\left(a_{2}-\left\\|K\right\\|_{\infty}\frac{2\alpha}{\sigma}\right)\mathbb{E}\left\\|\bar{\bm{\theta}}_{t}\right\\|^{2}+\alpha\left\\|K\right\\|_{\infty}\mathbb{E}\left\\|\bar{\bm{\theta}}_{t}\right\\|\mathbb{E}_{\bm{\theta}^{\prime}\sim\pi_{M,c}*\bar{\rho}_{t}}\left(\left(b_{1}+\frac{2}{\sigma}\right)\left\\|\bm{\theta}^{\prime}\right\\|+b_{1}\right)$
	$\displaystyle\leq$	$\displaystyle b_{1}-\left(a_{2}-\left\\|K\right\\|_{\infty}\frac{2\alpha}{\sigma}\right)U_{t}^{2}+\alpha\left\\|K\right\\|_{\infty}\mathbb{E}\left\\|\bar{\bm{\theta}}_{t}\right\\|\mathbb{E}_{\bm{\theta}^{\prime}\sim\pi_{M,c}*\bar{\rho}_{t}}\left(\left(b_{1}+\frac{2}{\sigma}\right)\left\\|\bm{\theta}^{\prime}\right\\|+b_{1}\right)$
	$\displaystyle\leq$	$\displaystyle b_{1}-\left(a_{2}-\left\\|K\right\\|_{\infty}\frac{2\alpha}{\sigma}\right)U_{t}^{2}+\alpha\left\\|K\right\\|_{\infty}\left(b_{1}+\frac{2}{\sigma}\right)\frac{1}{M}\sum_{j=1}^{M}U_{t}U_{t-jc}+\alpha\left\\|K\right\\|_{\infty}b_{1}U_{t}$
	$\displaystyle\leq$	$\displaystyle b_{1}-\left(a_{2}-\left\\|K\right\\|_{\infty}\frac{2\alpha}{\sigma}\right)U_{t}^{2}+\alpha\left\\|K\right\\|_{\infty}\left(b_{1}+\frac{2}{\sigma}\right)U_{t}^{2}+\alpha\left\\|K\right\\|_{\infty}b_{1}(U_{t}^{2}+1)$
	$\displaystyle\leq$	$\displaystyle\left(b_{1}+1\right)-\left(a_{2}-\left\\|K\right\\|_{\infty}\frac{2\alpha}{\sigma}-\alpha\left\\|K\right\\|_{\infty}\left(2b_{1}+\frac{2}{\sigma}\right)\right)U_{t}^{2},$

where $(1)$ is by E.3. By the assumption that $\lambda\coloneqq a_{2}-\left\|K\right\|_{\infty}\frac{2\alpha}{\sigma}-\alpha\left\|K\right\|_{\infty}\left(2b_{1}+\frac{2}{\sigma}\right)>0$ , we have

\frac{\partial}{\partial t}U_{t}^{2}\leq\left[\left(b_{1}+1\right)-\lambda U_{t}^{2}\right]\lor 0.

By Gronwall’s inequality, we have $U_{t}^{2}\leq U_{0}^{2}+\frac{b_{1}+1}{\lambda}$ . (If $\frac{\partial}{\partial t}U_{t}^{2}=0$ , then $U_{t}$ fix and this bound still holds.) Notice that in the first phase, as $\alpha=0$ , we have $\lambda<a_{2}$ 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 $U_{k}=\underset{s\leq k}{\sup}\sqrt{\mathbb{E}\left\|\tilde{\bm{\theta}}_{s}\right\|^{2}}$ , and thus

U_{k}^{2}-U_{k-1}^{2}\leq\left[2\eta\mathbb{E}\left\langle\tilde{\bm{\theta}}_{k-1},-\nabla V(\tilde{\bm{\theta}}_{k})+\alpha\phi[\pi_{M,c/\eta}*\tilde{\rho}_{k}](\tilde{\bm{\theta}}_{k})\right\rangle+2\eta^{2}\right]\lor 0.

By a similarly analysis, we have bound

		$\displaystyle\mathbb{E}\left\langle\tilde{\bm{\theta}}_{k-1},-\nabla V(\tilde{\bm{\theta}}_{k})+\alpha\phi[\pi_{M,c/\eta}*\tilde{\rho}_{k}](\tilde{\bm{\theta}}_{k})\right\rangle$
	$\displaystyle\leq$	$\displaystyle\left(b_{1}+1\right)-\lambda U_{t}^{2},$

where $\lambda=a_{2}-\left\|K\right\|_{\infty,\infty}\frac{2\alpha}{\sigma}-\alpha\left\|K\right\|_{\infty,\infty}\left(2b_{1}+\frac{2}{\sigma}\right)>0$ . And thus we have

U_{k}^{2}-U_{k-1}^{2}\leq\left[2\eta\left[\left(b_{1}+1\right)-\lambda U_{k-1}^{2}\right]+2\eta^{2}\right]\lor 0.

It gives that

U_{k}^{2}\leq\frac{b_{1}+1+\eta}{\lambda}+U_{0}^{2}.

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.

E.5.3 Proof of Theorem 4.1 and 5.1

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

\partial_{t}\rho_{t}=\nabla\cdot\left(\left[-\left(D(\bm{\theta})+Q(\bm{\theta})\right)\nabla V(\bm{\theta})+\alpha\phi[\pi_{M,c}*\rho_{t}](\bm{\theta}_{t})\right]\rho_{t}+\left(D(\bm{\theta})+Q(\bm{\theta})\right)\nabla\rho_{t}\right).

Notice that the limiting distribution satisfies

	$\displaystyle 0$	$\displaystyle\overset{a.e.}{=}\nabla\cdot\left(\left[-\left(D(\bm{\theta})+Q(\bm{\theta})\right)\nabla V(\bm{\theta})+\alpha\phi[\pi_{M,c}*\rho_{\infty}](\bm{\theta}_{t})\right]\rho_{\infty}+\left(D(\bm{\theta})+Q(\bm{\theta})\right)\nabla\rho_{\infty}\right)$
		$\displaystyle=\nabla\cdot\left(\left[-\left(D(\bm{\theta})+Q(\bm{\theta})\right)\nabla V(\bm{\theta})+\alpha\phi[\rho_{\infty}](\bm{\theta}_{t})\right]\rho_{\infty}+\left(D(\bm{\theta})+Q(\bm{\theta})\right)\nabla\rho_{\infty}\right)$
		$\displaystyle=\nabla\cdot\left(\left[-\left(D(\bm{\theta})+Q(\bm{\theta})\right)\nabla V(\bm{\theta})\right]\rho_{\infty}+\left(D(\bm{\theta})+Q(\bm{\theta})\right)\nabla\rho_{\infty}\right)$
		$\displaystyle+\alpha\nabla\cdot\left(K*\left(\nabla\rho_{\infty}-\nabla V(\bm{\theta})\rho_{\infty}\right)\rho_{\infty}\right).$

which implies that $\rho_{\infty}\propto\exp(-V(\bm{\theta}))$ is the stationary distribution.

E.5.4 Proof of Theorem 4.2

In the later proof we use $c_{d}$ to represent the quantity

\sqrt{\mathbb{E}_{\bm{\theta}\sim\rho_{0}}\left\|\bm{\theta}\right\|^{2}+\frac{b_{1}+1+\eta}{a_{2}-\left\|K\right\|_{\infty,\infty}\frac{2\alpha}{\sigma}-\alpha\left\|K\right\|_{\infty,\infty}\left(2b_{1}+\frac{2}{\sigma}\right)}}.

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=\eta N$ , for any $0\leq t\leq T$ , define $\underline{t}=\lfloor\frac{t}{\eta}\rfloor\eta$ . We introduce an another continuous-time interpolation dynamics:

	$\displaystyle\hat{\bm{\theta}}_{t}$	$\displaystyle=\begin{cases}-\nabla V(\hat{\bm{\theta}}_{\underline{t}})+d\mathcal{B}_{t},&t\in[0,Mc)\\ -\nabla V(\hat{\bm{\theta}}_{\underline{t}})+\alpha\phi[\pi_{M,c}*\hat{\rho}_{\underline{t}}](\hat{\bm{\theta}}_{\underline{t}})+d\mathcal{B}_{t},&t\geq Mc,\end{cases}$
	$\displaystyle\hat{\rho}_{t}$	$\displaystyle=\mathrm{Law}(\hat{\bm{\theta}}_{t}),$
	$\displaystyle\hat{\bm{\theta}}_{0}$	$\displaystyle=\bar{\bm{\theta}}_{0}\sim\bar{\rho}_{0},$

Notice that here we couples this interpolation dynamics with the same Brownian motion as that of the dynamics of $\bar{\bm{\theta}}_{t}$ . By the definition of $\hat{\bm{\theta}}_{t}$ , at any $t_{k}\coloneqq k\eta$ for some integrate $k\in[N]$ , $\hat{\bm{\theta}}_{t_{k}}$ and $\tilde{\bm{\theta}}_{k}$ has the same distribution. Define $\bar{\rho}_{t}^{\bm{\theta}_{0}}=\mathrm{Law}(\bar{\bm{\theta}}_{t})$ conditioning on $\bar{\bm{\theta}}_{0}=\bm{\theta}_{0}$ and $\hat{\rho}_{t}^{\bm{\theta}_{0}}=\mathrm{Law}(\hat{\bm{\theta}}_{t})$ conditioning on $\hat{\bm{\theta}}_{0}=\bm{\theta}_{0}$ . Followed by the argument of proving Lemma 2 in Dalalyan [2017], if $k\geq\frac{Mc}{\eta}$ , we have

		$\displaystyle\mathbb{D}_{\mathrm{KL}}\left[\bar{\rho}_{t_{k}}^{\bm{\theta}_{0}}\\|\hat{\rho}_{t_{k}}^{\bm{\theta}_{0}}\right]$
	$\displaystyle=$	$\displaystyle\frac{1}{4}\int_{0}^{t_{k}}\mathbb{E}\left\\|-\nabla V(\hat{\bm{\theta}}_{\underline{s}})+\alpha\phi[\pi_{M,c}\hat{\rho}_{\underline{s}}](\hat{\bm{\theta}}_{\underline{s}})+\nabla V(\hat{\bm{\theta}}_{s})-\alpha\phi[\pi_{M,c}\bar{\rho}_{s}](\hat{\bm{\theta}}_{s})\right\\|^{2}ds$
	$\displaystyle=$	$\displaystyle\frac{1}{4}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\mathbb{E}\left\\|-\nabla V(\hat{\bm{\theta}}_{t_{j}})+\alpha\phi[\pi_{M,c}\hat{\rho}_{t_{j}}](\hat{\bm{\theta}}_{t_{j}})+\nabla V(\hat{\bm{\theta}}_{s})-\alpha\phi[\pi_{M,c}\bar{\rho}_{s}](\hat{\bm{\theta}}_{s})\right\\|^{2}ds$
	$\displaystyle\leq$	$\displaystyle\frac{3}{4}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\mathbb{E}\left\\|\nabla V(\hat{\bm{\theta}}_{t_{j}})-\nabla V(\hat{\bm{\theta}}_{s})\right\\|^{2}ds$
	$\displaystyle+$	$\displaystyle\frac{3\alpha^{2}}{4}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\mathbb{E}\left\\|\phi[\pi_{M,c}\hat{\rho}_{t_{j}}](\hat{\bm{\theta}}_{t_{j}})-\phi[\pi_{M,c}\bar{\rho}_{s}](\hat{\bm{\theta}}_{t_{j}})\right\\|^{2}ds$
	$\displaystyle\leq$	$\displaystyle\frac{3\alpha^{2}}{4}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\mathbb{E}\left\\|\phi[\pi_{M,c}\bar{\rho}_{s}](\hat{\bm{\theta}}_{t_{j}})-\phi[\pi_{M,c}\bar{\rho}_{s}](\hat{\bm{\theta}}_{s})\right\\|^{2}ds$
	$\displaystyle=$	$\displaystyle I_{1}+I_{2}+I_{3}.$

We bound $I_{1}$ , $I_{2}$ and $I_{3}$ separately.

Bounding $I_{1}$ and $I_{3}$ By the smoothness of $\nabla V$ , we have

\left\|\nabla V(\hat{\bm{\theta}}_{t_{j}})-\nabla V(\hat{\bm{\theta}}_{s})\right\|^{2}\leq b_{1}^{2}\left\|\hat{\bm{\theta}}_{t_{j}}-\hat{\bm{\theta}}_{s}\right\|^{2}.

And by Lemma E.3 (Lipschitz of Stein Operator), we know that

		$\displaystyle\left\\|\phi[\pi_{M,c}\bar{\rho}_{s}](\bm{\theta}_{1})-\phi[\pi_{M,c}\bar{\rho}_{s}](\bm{\theta}_{2})\right\\|$
	$\displaystyle\leq$	$\displaystyle\left[\left\\|e^{-(\cdot)^{2}/\sigma}\right\\|_{\mathrm{Lip}}\mathbb{E}_{\bm{\theta}\sim\pi_{M,c}*\bar{\rho}_{s}}\left\\|\nabla V(\bm{\theta})\right\\|+\left\\|\frac{2}{\sigma}e^{-(\cdot)^{2}/\sigma}(\cdot)\right\\|_{\mathrm{Lip}}\right]\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}.$

And by the Assumption 4.2 and that $\bar{\rho}_{s}$ as finite second moment, we have

		$\displaystyle\left\\|\phi[\pi_{M,c}\bar{\rho}_{s}](\bm{\theta}_{1})-\phi[\pi_{M,c}\bar{\rho}_{s}](\bm{\theta}_{2})\right\\|$
	$\displaystyle\leq$	$\displaystyle Cc_{d}\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}.$

Combine the two bounds, we have

I_{1}+I_{3}\leq\frac{3Cc_{d}^{2}}{4}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\mathbb{E}\left\|\hat{\bm{\theta}}_{t_{j}}-\hat{\bm{\theta}}_{s}\right\|^{2}ds.

Notice that $\hat{\bm{\theta}}_{t}=\hat{\bm{\theta}}_{\underline{t}}+\left[-\nabla V(\hat{\bm{\theta}}_{\underline{t}})+\alpha\phi[\pi_{M,c}*\hat{\rho}_{\underline{t}}](\hat{\bm{\theta}}_{\underline{t}})\right](t-\underline{t})+\int_{\underline{t}}^{t}d\mathcal{B}_{s}$ . By Itô’s lemma, it implies that

	$\displaystyle I_{1}+I_{3}$	$\displaystyle\leq\frac{3Cc_{d}^{2}}{4}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\mathbb{E}\left\\|\hat{\bm{\theta}}_{t_{j}}-\hat{\bm{\theta}}_{s}\right\\|^{2}ds$
		$\displaystyle\leq\frac{3Cc_{d}^{2}}{4}\int_{t_{j}}^{t_{j+1}}\left[\mathbb{E}\left\\|-\nabla V(\hat{\bm{\theta}}_{\underline{s}})+\alpha\phi[\pi_{M,c}*\hat{\rho}_{\underline{s}}](\hat{\bm{\theta}}_{\underline{s}})\right\\|^{2}(s-t_{j})^{2}+2d(s-t_{j})\right]ds$
		$\displaystyle=Cc_{d}^{2}\eta^{3}\sum_{j=0}^{k-1}\mathbb{E}\left\\|-\nabla V(\hat{\bm{\theta}}_{t_{j}})+\alpha\phi[\pi_{M,c}*\hat{\rho}_{t_{j}}](\hat{\bm{\theta}}_{t_{j}})\right\\|^{2}+Cc_{d}^{2}dkh^{2}.$

By the assumption that $\mathbb{E}\left\|\tilde{\bm{\theta}}_{t_{j}}\right\|$ is finite and $\tilde{\bm{\theta}}_{t_{j}}\overset{d}{=}\hat{\bm{\theta}}_{t_{j}}$ , $\mathbb{E}\left\|\hat{\bm{\theta}}_{t_{j}}\right\|^{2}$ is also finite, we have

		$\displaystyle\mathbb{E}\left\\|-\nabla V(\hat{\bm{\theta}}_{\underline{t}})+\alpha\phi[\pi_{M,c}*\hat{\rho}_{\underline{t}}](\hat{\bm{\theta}}_{\underline{t}})\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle 2\mathbb{E}\left\\|\nabla V(\hat{\bm{\theta}}_{\underline{t}})\right\\|^{2}+2\alpha^{2}\mathbb{E}\left\\|\phi[\pi_{M,c}*\hat{\rho}_{\underline{t}}](\hat{\bm{\theta}}_{\underline{t}})\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle 4b_{1}^{2}+4b_{1}^{2}\mathbb{E}\left\\|\hat{\bm{\theta}}_{\underline{t}}\right\\|^{2}+2\alpha^{2}\mathbb{E}\left(\left(\frac{2}{\sigma}+b_{1}\right)\mathbb{E}_{\bm{\theta}^{\prime}\sim\pi_{M,c}*\hat{\rho}_{\underline{t}}}\left\\|\bm{\theta}^{\prime}\right\\|+\left\\|\bm{\theta}\right\\|\right)^{2}$
	$\displaystyle\leq$	$\displaystyle c_{d}^{2}C.$

Thus we conclude that

I_{1}+I_{3}\leq Cc_{d}^{2}\left(c_{d}^{2}k\eta^{3}+dk\eta^{2}\right).

Bounding $I_{2}$

		$\displaystyle\mathbb{E}\left\\|\phi[\pi_{M,c}\hat{\rho}_{t_{j}}](\hat{\bm{\theta}}_{t_{j}})-\phi[\pi_{M,c}\bar{\rho}_{s}](\hat{\bm{\theta}}_{t_{j}})\right\\|^{2}$
	$\displaystyle=$	$\displaystyle\mathbb{E}\left\\|\frac{1}{M}\sum_{l=1}^{M}\left[\phi[\hat{\rho}_{t_{j}-cl}](\hat{\bm{\theta}}_{t_{j}})-\phi[\bar{\rho}_{s-cl}](\hat{\bm{\theta}}_{t_{j}})\right]\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle\frac{1}{M}\sum_{l=1}^{M}\mathbb{E}\left\\|\phi[\hat{\rho}_{t_{j}-cl}](\hat{\bm{\theta}}_{t_{j}})-\phi[\bar{\rho}_{s-cl}](\hat{\bm{\theta}}_{t_{j}})\right\\|^{2}$
	$\displaystyle=$	$\displaystyle\frac{1}{M}\sum_{l=1}^{M}\mathbb{E}\left\\|\mathbb{E}_{\bm{\theta}\sim\hat{\rho}_{t_{j}-cl}}\bar{\phi}_{\hat{\bm{\theta}}_{t_{j}}}(\bm{\theta})-\mathbb{E}_{\bm{\theta}\sim\bar{\rho}_{s-cl}}\bar{\phi}_{\hat{\bm{\theta}}_{t_{j}}}(\bm{\theta})\right\\|^{2}$
	$\displaystyle=$	$\displaystyle\frac{1}{M}\sum_{l=1}^{M}\mathbb{E}_{\hat{\bm{\theta}}_{t_{j}}}\sum_{i=1}^{d}\left\|\mathbb{E}_{\bm{\theta}\sim\hat{\rho}_{t_{j}-cl}}\bar{\phi}_{\hat{\bm{\theta}}_{t_{j}},i}(\bm{\theta})-\mathbb{E}_{\bm{\theta}\sim\bar{\rho}_{s-cl}}\bar{\phi}_{\hat{\bm{\theta}}_{t_{j}},i}(\bm{\theta})\right\|^{2}$
	$\displaystyle\leq$	$\displaystyle\frac{1}{M}\sum_{l=1}^{M}\mathbb{E}_{\hat{\bm{\theta}}_{t_{j}}}\sum_{i=1}^{d}\left(\left\\|\bar{\phi}_{\hat{\bm{\theta}}_{t_{j}},i}(\cdot)\right\\|_{\mathcal{L}_{\infty}}\lor\left\\|\bar{\phi}_{\hat{\bm{\theta}}_{t_{j}},i}(\cdot)\right\\|_{\mathrm{Lip}}\right)^{2}\mathbb{D}_{\mathrm{BL}}^{2}\left[\hat{\rho}_{t_{j}-cl},\bar{\rho}_{s-cl}\right]$

By Lemma E.4 and the Assumption 4.4 that $V$ is at most quadratic growth and that $\hat{\rho}_{\underline{t}}$ has finite second moment, we have

		$\displaystyle\mathbb{E}_{\hat{\bm{\theta}}_{t_{j}}}\sum_{i=1}^{d}\left(\left\\|\bar{\phi}_{\hat{\bm{\theta}}_{t_{j}},i}(\cdot)\right\\|_{\mathcal{L}_{\infty}}\lor\left\\|\bar{\phi}_{\hat{\bm{\theta}}_{t_{j}},i}(\cdot)\right\\|_{\mathrm{Lip}}\right)^{2}$
	$\displaystyle=$	$\displaystyle\mathbb{E}_{\hat{\bm{\theta}}_{t_{j}}}\sum_{i=1}^{d}\left(\left\\|\bar{\phi}_{\hat{\bm{\theta}}_{t_{j}},i}(\cdot)\right\\|_{\mathcal{L}_{\infty}}^{2}\lor\left\\|\bar{\phi}_{\hat{\bm{\theta}}_{t_{j}},i}(\cdot)\right\\|_{\mathrm{Lip}}^{2}\right)$
	$\displaystyle\leq$	$\displaystyle\left[4d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathrm{BL}}^{2}+4\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathrm{BL}}^{2}\mathbb{E}_{\hat{\bm{\theta}}_{t_{j}}}\left\\|\nabla V(\hat{\bm{\theta}}_{t_{j}})\right\\|^{2}\right]$
	$\displaystyle\leq$	$\displaystyle C(d+c_{d}^{2}).$

Plug in the above estimation, we have

	$\displaystyle I_{2}$	$\displaystyle=\frac{3\alpha^{2}}{4}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\mathbb{E}\left\\|\phi[\pi_{M,c}\hat{\rho}_{t_{j}}](\hat{\bm{\theta}}_{t_{j}})-\phi[\pi_{M,c}\bar{\rho}_{s}](\hat{\bm{\theta}}_{t_{j}})\right\\|^{2}ds$
		$\displaystyle\leq\alpha^{2}C(d+c_{d}^{2})\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\frac{1}{M}\sum_{l=1}^{M}\mathbb{D}_{\mathrm{BL}}^{2}\left[\hat{\rho}_{t_{j}-cl},\bar{\rho}_{s-cl}\right]ds$
		$\displaystyle\leq\alpha^{2}C(d+c_{d}^{2})\sum_{j=0}^{k-1}\frac{1}{M}\sum_{l=1}^{M}\int_{t_{j}}^{t_{j+1}}\mathbb{D}_{\mathrm{KL}}\left[\hat{\rho}_{t_{j}-cl},\bar{\rho}_{s-cl}\right]ds,$

where the last inequality is due to the relation that $\mathbb{D}_{\mathrm{BL}}^{2}\overset{\mathrm{definition}}{\leq}\mathbb{D}_{\mathrm{TV}}^{2}\overset{\mathrm{Pinsker^{\prime}s}}{\leq}\mathbb{D}_{\mathrm{KL}}$ .

Overall Bound Combine all the estimation, we have

	$\displaystyle\mathbb{D}_{\mathrm{KL}}\left[\bar{\rho}_{t_{k}}^{\bm{\theta}_{0}}\\|\hat{\rho}_{t_{k}}^{\bm{\theta}_{0}}\right]$	$\displaystyle\leq\alpha^{2}C(d+c_{d}^{2})\sum_{j=0}^{k-1}\frac{1}{M}\sum_{l=1}^{M}\int_{t_{j}}^{t_{j+1}}\mathbb{D}_{\mathrm{KL}}\left[\hat{\rho}_{t_{j}-cl},\bar{\rho}_{s-cl}\right]ds+Cc_{d}^{2}\left(c_{d}^{2}k\eta^{3}+dk\eta^{2}\right)$
		$\displaystyle=\alpha^{2}C(d+c_{d}^{2})\sum_{j=0}^{k-1}\frac{1}{M}\sum_{l=1}^{M}\int_{0}^{\eta}\mathbb{D}_{\mathrm{KL}}\left[\hat{\rho}_{t_{\left(\frac{j\eta-cl}{\eta}\right)}},\bar{\rho}_{t_{\left(\frac{j\eta-cl}{\eta}\right)}+s}\right]ds+Cc_{d}^{2}\left(c_{d}^{2}k\eta^{3}+dk\eta^{2}\right)$

Similar, if $k\leq\frac{Mc}{\eta}-1$ , we have

		$\displaystyle\mathbb{D}_{\mathrm{KL}}\left[\bar{\rho}_{t_{k}}^{\bm{\theta}_{0}}\\|\hat{\rho}_{t_{k}}^{\bm{\theta}_{0}}\right]$
	$\displaystyle=$	$\displaystyle\frac{1}{4}\int_{0}^{t_{k}}\mathbb{E}\left\\|\nabla V(\hat{\bm{\theta}}_{\underline{s}})-\nabla V(\hat{\bm{\theta}}_{s})\right\\|^{2}ds$
	$\displaystyle\leq$	$\displaystyle\frac{b_{1}^{2}}{4}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\mathbb{E}\left\\|\hat{\bm{\theta}}_{t_{j}}-\hat{\bm{\theta}}_{s}\right\\|^{2}ds$
	$\leq$	$\displaystyle\frac{b_{1}^{2}\eta^{3}}{12}\sum_{j=0}^{k-1}\mathbb{E}\left\\|\nabla V(\hat{\bm{\theta}}_{t_{j}})\right\\|^{2}+\frac{dkb_{1}^{2}\eta^{2}}{4}$
	$\displaystyle\leq$	$\displaystyle\frac{b_{1}^{2}\eta^{3}kc_{d}^{2}}{12}+\frac{dkb_{1}^{2}\eta^{2}}{4}.$

Define

u_{k}=\underset{s\in\left[t_{k},t_{k+1}\right]}{\sup}\ \mathbb{D}_{\mathrm{KL}}\left[\bar{\rho}_{\underline{s}}^{\bm{\theta}_{0}}\|\hat{\rho}_{s}^{\bm{\theta}_{0}}\right],

and $U_{k}=\underset{l\in\{0,...,k\}}{\max}\ u_{l}$ . We conclude that for $k\geq\frac{Mc}{\eta}$ , for any $k^{\prime}\leq k$ ,

	$\displaystyle u_{k^{\prime}}$	$\displaystyle\leq\alpha^{2}C(d+c_{d}^{2})\sum_{j=0}^{k-1}\frac{1}{M}\sum_{l=1}^{M}\int_{0}^{h}\mathbb{D}_{\mathrm{KL}}\left[\hat{\rho}_{t_{\left(\frac{j\eta-cl}{\eta}\right)}},\bar{\rho}_{t_{\left(\frac{j\eta-cl}{\eta}\right)}+s}\right]ds+Cc_{d}^{2}\left(c_{d}^{2}k\eta^{3}+dk\eta^{2}\right)$
		$\displaystyle\leq\alpha^{2}C(d+c_{d}^{2})\sum_{j=0}^{k-1}\frac{1}{M}\sum_{l=1}^{M}\eta u_{\left(\frac{j\eta-cl}{\eta}\right)}+Cc_{d}^{2}\left(c_{d}^{2}k\eta^{3}+dk\eta^{2}\right)$
		$\displaystyle\leq\alpha^{2}C(d+c_{d}^{2})\eta\sum_{j=0}^{k-1}U_{j}+Cc_{d}^{2}\left(c_{d}^{2}k\eta^{3}+dk\eta^{2}\right).$

For $k<\frac{Mc}{\eta}$ , which is a simpler case, we have

U_{k}\leq C\left(\eta^{3}kc_{d}^{2}+dk\eta^{2}\right)<CMc\left(\eta c_{d}^{2}+d\right)\eta.

We bound the case when $k\geq\frac{Mc}{\eta}$ ,

U_{k}\leq\alpha^{2}C(d+c_{d}^{2})\eta\sum_{j=0}^{k-1}U_{j}+Cc_{d}^{2}\left(c_{d}^{2}k\eta^{3}+dk\eta^{2}\right).

If we take $\eta$ sufficiently small, such that $c_{d}^{2}k\eta^{3}\leq dk\eta^{2}$ , we have

	$\displaystyle U_{k}$	$\displaystyle\leq\alpha^{2}C(d+c_{d}^{2})\eta\sum_{j=0}^{k-1}U_{j}+2Cc_{d}^{2}dk\eta^{2}$
		$\displaystyle\leq\alpha^{2}C(d+c_{d}^{2})\eta\sum_{j=0}^{k-1}\left(U_{j}+\eta\right).$

Define $\eta^{\prime}=\alpha^{2}C(d+c_{d}^{2})\eta$ and we can choose $\eta$ small enough such that $\eta^{\prime}<1/2$ and $\eta<1/2$ . Without loss of generality, we also assume $\eta^{\prime}\geq\eta$ and thus we have

U_{k}\leq\eta^{\prime}\sum_{j=0}^{k-1}\left(U_{j}+\eta^{\prime}\right).

Also we assume $U_{k}\geq\eta^{\prime}$ , otherwise we conclude that $U_{k}<\eta^{\prime}$ . We thus have $U_{k}\leq q\sum_{j=0}^{k-1}U_{j}$ , where $q=2\eta^{\prime}$ . Suppose that $U_{\frac{Mc}{\eta}-1}=x\leq CMc\left(\eta c_{d}^{2}+d\right)\eta$ and some algebra (which reduces to Pascal’s triangle) shows that

U_{k}\leq xq(1+q)^{k-\frac{Mc}{\eta}}.

We conclude that $U_{k}\leq xq(1+q)^{k-1}$ . Notice that $q=2\alpha^{2}C(d+c_{d}^{2})\eta.$ Thus for any $k\geq Mc/\eta$ ,

	$\displaystyle U_{k}$	$\displaystyle\leq xq(1+q)^{k-\frac{Mc}{\eta}}$
		$\displaystyle=xq(1+q)^{\left(k\eta-Mc\right)/\eta}$
		$\displaystyle=xq(1+q)^{2\alpha^{2}C(d+c_{d}^{2})\left(k\eta-Mc\right)/q}$
		$\displaystyle\leq x2\alpha^{2}C(d+c_{d}^{2})e^{2\alpha^{2}C(d+c_{d}^{2})\left(k\eta-Mc\right)}\eta$
		$\displaystyle\leq CMc\alpha^{2}\left(\eta c_{d}^{2}+d\right)(d+c_{d}^{2})e^{2\alpha^{2}C(d+c_{d}^{2})\left(k\eta-Mc\right)}\eta^{2},$

for sufficiently small $\eta$ . Combine the above two estimations, we have

U_{k}\leq\begin{cases}C\left(\eta^{3}kc_{d}^{2}+dk\eta^{2}+\eta\right)&k\leq Mc/\eta-1\\ CMc\alpha^{2}\left(\eta c_{d}^{2}+d\right)(d+c_{d}^{2})e^{2\alpha^{2}C(d+c_{d}^{2})\left(k\eta-Mc\right)}\eta^{2}+C\eta&k\geq Mc/\eta\end{cases}.

Notice that now we have $U_{k}=\underset{l\in\{0,...,k\}}{\max}\underset{s\in\left[0,\eta\right]}{\sup}\ \mathbb{D}_{\mathrm{KL}}\left[\bar{\rho}_{l\eta+s}^{\bm{\theta}_{0}}\|\tilde{\rho}_{l\eta}^{\bm{\theta}_{0}}\right]$ , which is a function of $\bm{\theta}_{0}$ . We then bound $\bar{U}_{k}=\underset{l\in\{0,...,k\}}{\max}\underset{s\in\left[0,\eta\right]}{\sup}\ \mathbb{D}_{\mathrm{KL}}\left[\bar{\rho}_{l\eta+s}\|\tilde{\rho}_{l\eta}\right]$ . Notice that the KL divergence has the following variational representation:

\mathbb{D}_{\mathrm{KL}}[\rho_{1}\|\rho_{2}]=\sup_{f}\left[\mathbb{E}_{\rho_{1}}f-\mathbb{E}_{\rho_{2}}e^{f}\right],

where the $f$ is chosen in the set that $\mathbb{E}_{\rho_{1}}f$ and $\mathbb{E}_{\rho_{2}}e^{f}$ exist. And thus we have

	$\displaystyle\mathbb{D}_{\mathrm{KL}}[\bar{\rho}_{l\eta+s}\\|\tilde{\rho}_{l\eta}]$	$\displaystyle=\sup_{f}\left[\mathbb{E}_{\bm{\theta}_{0}\sim\rho_{0}}\left(\mathbb{E}_{\bar{\rho}_{l\eta+s}^{\bm{\theta}_{0}}}f-\mathbb{E}_{\tilde{\rho}_{l\eta}^{\bm{\theta}_{0}}}e^{f}\right)\right]$
		$\displaystyle\leq\mathbb{E}_{\bm{\theta}_{0}\sim\rho_{0}}\sup_{f}\left[\left(\mathbb{E}_{\bar{\rho}_{l\eta+s}^{\bm{\theta}_{0}}}f-\mathbb{E}_{\tilde{\rho}_{l\eta}^{\bm{\theta}_{0}}}e^{f}\right)\right].$

And thus $\bar{U}_{k}\leq U_{k}$ . Also the inequality that

\bar{U}_{k}=\underset{l\in\{0,...,k\}}{\max}\underset{s\in\left[0,\eta\right]}{\sup}\ \mathbb{D}_{\mathrm{KL}}\left[\bar{\rho}_{l\eta+s}\|\tilde{\rho}_{l\eta}\right]\geq\underset{l\in\{0,...,k\}}{\max}\ \mathbb{D}_{\mathrm{KL}}\left[\bar{\rho}_{l\eta}\|\tilde{\rho}_{l\eta}\right]

holds naturally by definition. We complete the proof.

E.5.5 Proof of Theorem 4.3

The constant $h_{1}$ is defined as

h_{1}=\left\|\frac{2}{\sigma}e^{-\left\|\bm{\theta}\right\|^{2}/\sigma}\theta_{1}\right\|_{\mathrm{BL}}^{2}\lor\left\|e^{-\left\|\cdot\right\|^{2}/\sigma}\right\|_{\mathrm{BL}}^{2}\lor\left\|\frac{2}{\sigma}e^{-(\cdot)^{2}/\sigma}(\cdot)\right\|_{\mathrm{Lip}}

Now we start the proof. We couple the process of $\bm{\theta}_{k}$ and $\tilde{\bm{\theta}}_{k}$ by the same gaussian noise $\bm{e}_{k}$ in every iteration and same initialization $\tilde{\bm{\theta}}_{0}=\bm{\theta}_{0}$ . For $k\leq Mc/\eta-1$ , $\mathbb{E}\left\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\|^{2}=0$ and for $k\geq Mc/\eta$ we have the following inequality,

		$\displaystyle\mathbb{E}\left\\|\bm{\theta}_{k+1}-\tilde{\bm{\theta}}_{k+1}\right\\|^{2}-\mathbb{E}\left\\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\\|^{2}$
	$\displaystyle=$	$\displaystyle 2\eta\mathbb{E}\left\langle\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k},-\nabla V(\bm{\theta}_{k})+\nabla V(\tilde{\bm{\theta}}_{k})\right\rangle$
	$\displaystyle+$	$\displaystyle 2\eta\alpha\mathbb{E}\left\langle\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k},\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}*\tilde{\rho}_{k}](\tilde{\bm{\theta}}_{k})\right\rangle$
	$\displaystyle+$	$\displaystyle\eta^{2}\mathbb{E}\left\\|-\nabla V(\bm{\theta}_{k})+\alpha\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta}_{k})+\nabla V(\tilde{\bm{\theta}}_{k})-\alpha\phi[\pi_{M,c/\eta}*\tilde{\rho}_{k}](\tilde{\bm{\theta}}_{k})\right\\|^{2}.$

By the log-concavity, we have

		$\displaystyle\mathbb{E}\left\langle\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k},-\nabla V(\bm{\theta}_{k})+\nabla V(\tilde{\bm{\theta}}_{k})\right\rangle$
	$\displaystyle\leq$	$\displaystyle-L\mathbb{E}\left\\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\\|^{2},$

for some positive constant $L$ . And also, as $\eta$ 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

		$\displaystyle\mathbb{E}\left\langle\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k},\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc}}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}*\tilde{\rho}_{k}](\tilde{\bm{\theta}}_{k})\right\rangle$
	$\displaystyle=$	$\displaystyle\mathbb{E}\left\langle\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k},\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}*\rho_{k}](\bm{\theta}_{k})\right\rangle$
	$\displaystyle+$	$\displaystyle\mathbb{E}\left\langle\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k},\phi[\pi_{M,c/\eta}\rho_{k}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}\tilde{\rho}_{k}](\bm{\theta}_{k})\right\rangle$
	$\displaystyle+$	$\displaystyle\mathbb{E}\left\langle\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k},\phi[\pi_{M,c/\eta}\tilde{\rho}_{k}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}\tilde{\rho}_{k}](\tilde{\bm{\theta}}_{k})\right\rangle$
	$\displaystyle=$	$\displaystyle I_{1}+I_{2}+I_{3}.$

We bound $I_{1}$ , $I_{2}$ and $I_{3}$ independently.

Bounding $I_{1}$

By Holder’s inequality,

	$\displaystyle I_{1}$	$\displaystyle\leq\mathbb{E}\left[\left\\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\\|\left\\|\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}*\rho_{k}](\bm{\theta}_{k})\right\\|\right]$
		$\displaystyle\leq\sqrt{\mathbb{E}\left\\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\\|^{2}}\sqrt{\mathbb{E}\left\\|\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}*\rho_{k}](\bm{\theta}_{k})\right\\|^{2}}.$

We bound the second term on the right side of the inequality. Define

a_{2}=\underset{k}{\sup}\ \frac{\mathbb{E}\left\|\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}*\rho_{k}](\bm{\theta}_{k})\right\|^{2}}{\underset{\left\|\bm{\theta}\right\|\leq B}{\sup}\mathbb{E}\left\|\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta})-\phi[\pi_{M,c/\eta}*\rho_{k}](\bm{\theta})\right\|^{2}}

and by the regularity assumption we know that $a_{2}<\infty$ . Define $\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta})-\phi[\pi_{M,c/\eta}*\rho_{k}](\bm{\theta})=\phi^{*}[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}]$ and since the stein operator is linear functional of the distribution, we have

\mathbb{E}\phi^{*}[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta})=0,

given any $\bm{\theta}$ . By Theorem E.1 that $\bm{\Theta}_{k}$ is geometric ergodicity and thus is $\beta$ -mixing with exponentially fast decay rate by Proposition E.2. And by Proposition E.1, we know that $\bm{\Theta}_{k}$ is also $\alpha$ -mixing with exponentially fast decay rate. We have the following estimation

		$\displaystyle\mathbb{E}\left\\|\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}*\rho_{k}](\bm{\theta}_{k})\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle a_{2}\underset{\left\\|\bm{\theta}\right\\|\leq B}{\sup}\mathbb{E}\left\\|\phi^{*}[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta})\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle\frac{a_{2}}{M^{2}}\underset{\left\\|\bm{\theta}\right\\|\leq B}{\sup}\mathbb{E}\sum_{k=1}^{M}\left\\|\phi^{*}[\delta_{\bm{\theta}_{t-kc/\eta}}](\bm{\theta})\right\\|^{2}$
	$\displaystyle+$	$\displaystyle\frac{a_{2}}{M^{2}}\underset{\left\\|\bm{\theta}\right\\|\leq B}{\sup}\mathbb{E}\sum_{k\neq j}\left\langle\phi^{}[\delta_{\bm{\theta}_{t-kc/\eta}}](\bm{\theta}),\phi^{}[\delta_{\bm{\theta}_{t-jc/\eta}}](\bm{\theta})\right\rangle$
	$\displaystyle\leq$	$\displaystyle\frac{Ca_{2}}{M}\left[\frac{e^{-rc}\left(1-e^{-rMc}\right)}{1-e^{rc}}+1\right],$

for some positive constant $r$ that characterize the decay rate of $\alpha$ mixing. Notice that here $\eta$ is canceled because the decay rate of mixing is $\mathcal{O}(\eta)$ (on the power of exponential) and $c/\eta=\mathcal{O}(\eta^{-1})$ . Combine this two estimations, we have

I_{1}\leq\sqrt{\mathbb{E}\left\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\|^{2}}\sqrt{\frac{a_{2}C}{M}\left[\frac{e^{-rc}\left(1-e^{-rMc}\right)}{1-e^{rc}}+1\right]}.

Bounding $I_{2}$ By Holder’s inequality, we have

\displaystyle I_{2}

\displaystyle\leq\sqrt{\mathbb{E}\left\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\|^{2}}\sqrt{\mathbb{E}\left\|\phi[\pi_{M,c/\eta}*\rho_{k}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}*\tilde{\rho}_{k}](\bm{\theta}_{k})\right\|^{2}}.

We bound the second term in the right side of the inequality.

		$\displaystyle\mathbb{E}\left\\|\phi[\pi_{M,c/\eta}\rho_{k}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}\tilde{\rho}_{k}](\bm{\theta}_{k})\right\\|^{2}$
	$\displaystyle=$	$\displaystyle\mathbb{E}\left\\|\frac{1}{M}\sum_{j=1}^{M}\left[\phi[\rho_{k-jc/\eta}](\bm{\theta}_{k})-\phi[\tilde{\rho}_{k-jc/\eta}](\bm{\theta}_{k})\right]\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle\frac{1}{M}\sum_{j=1}^{M}\mathbb{E}\left\\|\phi[\rho_{k-jc/\eta}](\bm{\theta}_{k})-\phi[\tilde{\rho}_{k-jc/\eta}](\bm{\theta}_{k})\right\\|^{2}$
	$\displaystyle=$	$\displaystyle\frac{1}{M}\sum_{j=1}^{M}\mathbb{E}_{\bm{\theta}_{k}}\left\\|\mathbb{E}_{\bm{\theta}\sim\rho_{k-jc/\eta}}\bar{\phi}_{\bm{\theta}_{k}}(\bm{\theta})-\mathbb{E}_{\bm{\theta}\sim\tilde{\rho}_{k-jc/\eta}}\bar{\phi}_{\bm{\theta}_{k}}(\bm{\theta})\right\\|^{2}$
	$\displaystyle=$	$\displaystyle\frac{1}{M}\sum_{j=1}^{M}\mathbb{E}_{\bm{\theta}_{k}}\sum_{i=1}^{d}\left\|\mathbb{E}_{\bm{\theta}\sim\rho_{k-jc/\eta}}\bar{\phi}_{\bm{\theta}_{k},i}(\bm{\theta})-\mathbb{E}_{\bm{\theta}\sim\tilde{\rho}_{k-jc/\eta}}\bar{\phi}_{\bm{\theta}_{k},i}(\bm{\theta})\right\|^{2}$
	$\displaystyle\leq$	$\displaystyle\frac{1}{M}\sum_{j=1}^{M}\mathbb{E}_{\bm{\theta}_{k}}\sum_{i=1}^{d}\left(\left\\|\bar{\phi}_{\bm{\theta}_{k},i}(\cdot)\right\\|_{\mathcal{L}_{\infty}}\lor\left\\|\bar{\phi}_{\bm{\theta}_{k},i}(\cdot)\right\\|_{\mathrm{Lip}}\right)^{2}\mathbb{D}_{\mathrm{BL}}^{2}\left[\rho_{k-jc/\eta},\tilde{\rho}_{k-jc/\eta}\right].$

By Lemma E.4, we have

		$\displaystyle\sum_{i=1}^{d}\left(\left\\|\bar{\phi}_{\bm{\theta}_{k},i}(\cdot)\right\\|_{\mathcal{L}_{\infty}}\lor\left\\|\bar{\phi}_{\bm{\theta}_{k},i}(\cdot)\right\\|_{\mathrm{Lip}}\right)^{2}$
	$\displaystyle=$	$\displaystyle\sum_{i=1}^{d}\left(\left\\|\bar{\phi}_{\bm{\theta}_{k},i}(\cdot)\right\\|_{\mathcal{L}_{\infty}}^{2}\lor\left\\|\bar{\phi}_{\bm{\theta}_{k},i}(\cdot)\right\\|_{\mathrm{Lip}}^{2}\right)$
	$\displaystyle\leq$	$\displaystyle\left[4d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathrm{BL}}^{2}+4\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathrm{BL}}^{2}\left\\|\nabla V(\bm{\theta}_{k})\right\\|^{2}\right].$

Plug in the above estimation and by the relation that $\mathbb{D}_{\mathrm{BL}}\leq\mathbb{W}_{1}\leq\mathbb{W}_{2}$ , we have

		$\displaystyle\mathbb{E}\left\\|\phi[\pi_{M,c}\rho_{k}](\bm{\theta}_{k})-\phi[\pi_{M,c}\tilde{\rho}_{k}](\bm{\theta}_{k})\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle\left[4d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathrm{BL}}^{2}+4\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathrm{BL}}^{2}\mathbb{E}_{\bm{\theta}_{k}}\left\\|\nabla V(\bm{\theta}_{k})\right\\|^{2}\right]\frac{1}{M}\sum_{j=1}^{M}\mathbb{D}_{\mathrm{BL}}^{2}\left[\rho_{k-cj},\tilde{\rho_{k-cj}}\right]$
	$\displaystyle\leq$	$\displaystyle\left[4d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathrm{BL}}^{2}+4\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathrm{BL}}^{2}\mathbb{E}_{\bm{\theta}_{k}}\left\\|\nabla V(\bm{\theta}_{k})\right\\|^{2}\right]\frac{1}{M}\sum_{j=1}^{M}\mathbb{W}_{2}^{2}\left[\rho_{k-cj},\tilde{\rho_{k-cj}}\right].$

And combined all the estimation and by the definition of Wasserstein-distance, we conclude that

	$\displaystyle I_{2}$	$\displaystyle\leq\sqrt{4d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathrm{BL}}^{2}+4\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathrm{BL}}^{2}\mathbb{E}_{\bm{\theta}_{k}}\left\\|\nabla V(\bm{\theta}_{k})\right\\|^{2}}\sqrt{\frac{1}{M}\sum_{j=1}^{M}\mathbb{W}_{2}^{2}\left[\rho_{k-cj},\tilde{\rho}_{k-cj}\right]}$
		$\displaystyle\leq\sqrt{4d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathrm{BL}}^{2}+4\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathrm{BL}}^{2}\mathbb{E}_{\bm{\theta}_{k}}\left\\|\nabla V(\bm{\theta}_{k})\right\\|^{2}}\sqrt{\frac{1}{M}\sum_{j=1}^{M}\mathbb{E}\left\\|\bm{\theta}_{k-cj}-\tilde{\bm{\theta}}_{k-cj}\right\\|^{2}}.$

Bounding $I_{3}$

By Holder’s inequality,

\displaystyle I_{3}

\displaystyle\leq\sqrt{\mathbb{E}\left\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\|^{2}}\sqrt{\mathbb{E}\left\|\phi[\pi_{M,c/\eta}*\tilde{\rho}_{k}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}*\tilde{\rho}_{k}](\tilde{\bm{\theta}}_{k})\right\|^{2}}.

We bound the last term on the right side of the inequality. By assumption and Lemma E.3, we have

		$\displaystyle\mathbb{E}\left\\|\phi[\pi_{M,c/\eta}\tilde{\rho}_{k}](\bm{\theta}_{k})-\phi[\pi_{M,c/\eta}\tilde{\rho}_{k}](\tilde{\bm{\theta}}_{k})\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle\left[\left\\|e^{-(\cdot)^{2}/\sigma}\right\\|_{\mathrm{Lip}}\mathbb{E}_{\bm{\theta}\sim\tilde{\rho}_{k}}\left\\|\nabla V(\bm{\theta})\right\\|+\left\\|\frac{2}{\sigma}e^{-(\cdot)^{2}/\sigma}(\cdot)\right\\|_{\mathrm{Lip}}\right]^{2}\mathbb{E}\left\\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\\|^{2}.$

And combine the estimation, we have

\displaystyle I_{3}

\displaystyle\leq\left[\left\|e^{-(\cdot)^{2}/\sigma}\right\|_{\mathrm{Lip}}\mathbb{E}_{\bm{\theta}\sim\tilde{\rho}_{k}}\left\|\nabla V(\bm{\theta})\right\|+\left\|\frac{2}{\sigma}e^{-(\cdot)^{2}/\sigma}(\cdot)\right\|_{\mathrm{Lip}}\right]\mathbb{E}\left\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\|^{2}.

Overall Bound Combine all the results, we have the following bound: for $k\geq Mc$ ,

		$\displaystyle\mathbb{E}\left\\|\bm{\theta}_{k+1}-\tilde{\bm{\theta}}_{k+1}\right\\|^{2}-\mathbb{E}\left\\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle-2\eta L\mathbb{E}\left\\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\\|^{2}$
	$\displaystyle+$	$\displaystyle 2\eta\alpha\sqrt{\mathbb{E}\left\\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\\|^{2}}\frac{c_{1}}{\sqrt{M}}$
	$\displaystyle+$	$\displaystyle 2\text{$\eta\alpha$}c_{2}\sqrt{\frac{1}{M}\sum_{j=1}^{M}\mathbb{E}\left\\|\bm{\theta}_{k-jc/\eta}-\tilde{\bm{\theta}}_{k-jc/\eta}\right\\|^{2}\mathbb{E}\left\\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\\|^{2}}$
	$\displaystyle+$	$\displaystyle 2\eta\alpha c_{3}\mathbb{E}\left\\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\\|^{2}$
	$\displaystyle+$	$\displaystyle\eta^{2}c_{4},$

where

c_{1}=\sqrt{a_{2}C\left[\frac{e^{-rc}\left(1-e^{-rMc}\right)}{1-e^{rc}}+1\right]},

c_{2}=\sqrt{4d\left\|\frac{2}{\sigma}e^{-\left\|\bm{\theta}\right\|^{2}/\sigma}\theta_{1}\right\|_{\mathrm{BL}}^{2}+4\left\|e^{-\left\|\cdot\right\|^{2}/\sigma}\right\|_{\mathrm{BL}}^{2}\underset{k}{\sup}\ \mathbb{E}_{\bm{\theta}_{k}}\left\|\nabla V(\bm{\theta}_{k})\right\|^{2}},

c_{3}=\left[\left\|e^{-(\cdot)^{2}/\sigma}\right\|_{\mathrm{Lip}}\underset{k}{\sup}\ \mathbb{E}_{\bm{\theta}\sim\tilde{\rho}_{k}}\left\|\nabla V(\bm{\theta})\right\|+\left\|\frac{2}{\sigma}e^{-(\cdot)^{2}/\sigma}(\cdot)\right\|_{\mathrm{Lip}}\right],

and

c_{4}=\underset{k\geq Mc/\eta}{\sup}\ \mathbb{E}\left\|\nabla V(\bm{\theta}_{k})+\alpha\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{k-jc/\eta}}](\bm{\theta}_{k})-\nabla V(\tilde{\bm{\theta}}_{k})-\alpha\phi[\pi_{M,c/\eta}*\tilde{\rho}_{k}](\tilde{\bm{\theta}}_{k})\right\|^{2}.

Define $u_{k}=\sqrt{\mathbb{E}\left\|\bm{\theta}_{k}-\tilde{\bm{\theta}}_{k}\right\|^{2}}$ and $U_{k}=\underset{l\in[k]}{\sup}\ u_{l}$ , we have

U_{k+1}^{2}\leq qU_{k}^{2}+\frac{2\eta\alpha c_{1}}{\sqrt{M}}U_{k}+\eta^{2}c_{4},

where $q=(1-2\eta(L-\alpha c_{2}-\alpha c_{3}))$ . By the assumption that $\alpha\leq L/(c_{2}+c_{3})$ , $q<1$ . Now we prove the bound of $U_{k}$ by induction. We take the hypothesis that $U_{k}^{2}\leq\frac{\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}+(1-q)\eta\left(c_{4}+\frac{1}{1-q}\right)\right)^{2}}{\left(1-q\right)^{2}}$ and notice that the hypothesis holds for $U_{0}=0$ . By the hypothesis, we have

	$\displaystyle U_{k+1}^{2}$	$\displaystyle\leq q\frac{\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}+(1-q)\eta\left(c_{4}+\frac{1}{1-q}\right)\right)^{2}}{\left(1-q\right)^{2}}+\frac{2\eta\alpha c_{1}}{\sqrt{M}}\frac{\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}+(1-q)\eta\left(c_{4}+\frac{1}{1-q}\right)\right)}{\left(1-q\right)}+\eta^{2}\left(c_{4}+\frac{1}{1-q}\right)$
		$\displaystyle=q\frac{\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}+(1-q)\eta\left(c_{4}+\frac{1}{1-q}\right)\right)^{2}}{\left(1-q\right)^{2}}+\frac{1}{1-q}\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}\right)^{2}+\frac{2\eta\alpha c_{1}}{\sqrt{M}}\eta\left(c_{4}+\frac{1}{1-q}\right)+\eta^{2}\left(c_{4}+\frac{1}{1-q}\right)$
		$\displaystyle=q\frac{\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}+(1-q)\eta\left(c_{4}+\frac{1}{1-q}\right)\right)^{2}}{\left(1-q\right)^{2}}$
		$\displaystyle+\frac{1-q}{\left(1-q\right)^{2}}\left[\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}\right)^{2}+(1-q)\frac{2\eta\alpha c_{1}}{\sqrt{M}}\eta\left(c_{4}+\frac{1}{1-q}\right)+(1-q)\eta^{2}\left(c_{4}+\frac{1}{1-q}\right)\right]$
		$\displaystyle\leq q\frac{\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}+(1-q)\eta\left(c_{4}+\frac{1}{1-q}\right)\right)^{2}}{\left(1-q\right)^{2}}$
		$\displaystyle+\frac{1-q}{\left(1-q\right)^{2}}\left[\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}\right)^{2}+(1-q)\frac{2\eta\alpha c_{1}}{\sqrt{M}}\eta\left(c_{4}+\frac{1}{1-q}\right)+(1-q)^{2}\eta^{2}\left(c_{4}+\frac{1}{1-q}\right)^{2}\right]$
		$\displaystyle=q\frac{\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}+(1-q)\eta\left(c_{4}+\frac{1}{1-q}\right)\right)^{2}}{\left(1-q\right)^{2}}+(1-q)\frac{\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}+(1-q)\eta\left(c_{4}+\frac{1}{1-q}\right)\right)^{2}}{\left(1-q\right)^{2}}$
		$\displaystyle=\frac{\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}+(1-q)\eta\left(c_{4}+\frac{1}{1-q}\right)\right)^{2}}{\left(1-q\right)^{2}},$

where the last second inequality holds by $(1-q)\left(c_{4}+\frac{1}{1-q}\right)\geq 1$ . Thus we complete the argument of induction and we have, for any $k$ ,

	$\displaystyle U_{k}^{2}$	$\displaystyle\leq\frac{\left(\frac{2\eta\alpha c_{1}}{\sqrt{M}}+(1-q)\eta\left(c_{4}+\frac{1}{1-q}\right)\right)^{2}}{\left(1-q\right)^{2}}$
		$\displaystyle\leq 2\frac{\frac{4\eta^{2}\alpha^{2}c_{1}^{2}}{M}+2(1-q)^{2}\eta^{2}\left(c_{4}+\frac{1}{1-q}\right)^{2}}{\left(1-q\right)^{2}}$
		$\displaystyle=\frac{2\alpha^{2}c_{1}^{2}}{(L-\alpha c_{2}-\alpha c_{3})^{2}}\frac{1}{M}+4\eta^{2}\left(c_{4}+2\eta(L-\alpha c_{2}-\alpha c_{3})\right)^{2}.$

And it implies that $\mathbb{W}_{2}^{2}[\rho_{k},\tilde{\rho}_{k}]\leq u_{k}\leq U_{k}\leq\frac{2\alpha^{2}c_{1}^{2}}{(L-\alpha c_{2}-\alpha c_{3})^{2}}\frac{1}{M}+4\eta^{2}\left(c_{4}+2\eta(L-\alpha c_{2}-\alpha c_{3})\right)^{2}.$

E.6 Proof of Technical Lemmas

E.6.1 Proof of Lemma E.1

For the first part:

		$\displaystyle\left\\|\nabla V(\bm{\theta})\right\\|$
	$\displaystyle\leq$	$\displaystyle\left\\|\nabla V(\bm{\theta})-\nabla V(\bm{0})\right\\|+\left\\|\nabla V(\bm{0})\right\\|$
	$\displaystyle\leq$	$\displaystyle b_{1}\left(\left\\|\bm{\theta}_{1}\right\\|+1\right).$

For the second part:

		$\displaystyle\left\\|\bm{\theta}-\eta\nabla V(\bm{\theta})\right\\|$
	$\displaystyle=$	$\displaystyle\left\langle\bm{\theta}-\eta\nabla V(\bm{\theta}),\bm{\theta}-\eta\nabla V(\bm{\theta})\right\rangle$
	$\displaystyle=$	$\displaystyle\left\\|\bm{\theta}\right\\|^{2}+2\eta\left\langle\bm{\theta},-\nabla V(\bm{\theta})\right\rangle+\eta^{2}\left\\|\nabla V(\bm{\theta})\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle\left\\|\bm{\theta}\right\\|^{2}+2\eta\left(-a_{1}\left\\|\bm{\theta}\right\\|^{2}+b_{1}\right)+\eta^{2}b_{1}(1+\left\\|\bm{\theta}\right\\|^{2})$
	$\displaystyle=$	$\displaystyle\left(1-2\eta a_{1}+\eta^{2}b_{1}\right)\left\\|\bm{\theta}\right\\|^{2}+\eta^{2}b_{1}+2\eta b_{1}.$

E.6.2 Proof of Lemma E.2

It is obvious that $\left\|K\right\|_{\infty,\infty}\leq 1$ .

		$\displaystyle\left\\|K(\bm{\theta}^{\prime},\bm{\theta}_{1})-K(\bm{\theta}^{\prime},\bm{\theta}_{2})\right\\|$
		$\displaystyle\left\\|e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{1}\right\\|^{2}/\sigma}-e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{2}\right\\|^{2}/\sigma}\right\\|$
	$\displaystyle\leq$	$\displaystyle\left\\|e^{-(\cdot)^{2}/\sigma}\right\\|_{\mathrm{Lip}}\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}.$

And

		$\displaystyle\left\\|\nabla_{\bm{\theta}^{\prime}}K(\bm{\theta}^{\prime},\bm{\theta}_{1})-\nabla_{\bm{\theta}^{\prime}}K(\bm{\theta}^{\prime},\bm{\theta}_{2})\right\\|$
	$\displaystyle=$	$\displaystyle\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{1}\right\\|^{2}/\sigma}\left(\bm{\theta}^{\prime}-\bm{\theta}_{1}\right)-\frac{2}{\sigma}e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{2}\right\\|^{2}/\sigma}\left(\bm{\theta}^{\prime}-\bm{\theta}_{2}\right)\right\\|$
	$\displaystyle\leq$	$\displaystyle\left\\|\frac{2}{\sigma}e^{-(\cdot)^{2}/\sigma}(\cdot)\right\\|_{\mathrm{Lip}}\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}.$

E.6.3 Proof of Lemma E.3

For any distribution $\rho$ such that $\int\left\|\nabla_{\bm{\theta}}V(\bm{\theta})\right\|\rho(\bm{\theta})d\bm{\theta}<\infty$ ,

		$\displaystyle\left\\|\phi[\rho](\bm{\theta}_{1})-\phi[\rho](\bm{\theta}_{2})\right\\|$
	$\displaystyle=$	$\displaystyle\left\\|\mathbb{E}_{\bm{\theta}\sim\rho}\left\{-\left[K(\bm{\theta},\bm{\theta}_{1})-K(\bm{\theta},\bm{\theta}_{2})\right]\nabla V(\bm{\theta})+\nabla_{1}K(\bm{\theta},\bm{\theta}_{1})-\nabla_{1}K(\bm{\theta},\bm{\theta}_{2})\right\}\right\\|$
	$\displaystyle\leq$	$\displaystyle\left\\|e^{-(\cdot)^{2}/\sigma}\right\\|_{\mathrm{Lip}}\mathbb{E}_{\bm{\theta}\sim\rho}\left\\|\nabla V(\bm{\theta})\right\\|\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}$
	$\displaystyle+$	$\displaystyle\left\\|\frac{2}{\sigma}e^{-(\cdot)^{2}/\sigma}(\cdot)\right\\|_{\mathrm{Lip}}\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}.$

For proving the second result, we notice that

	$\displaystyle\left\\|\phi[\rho](\bm{\theta})\right\\|$	$\displaystyle=\mathbb{E}_{\bm{\theta}^{\prime}\sim\rho}\left[K(\bm{\theta}^{\prime},\bm{\theta})\nabla V(\bm{\theta}^{\prime})+\nabla_{1}K(\bm{\theta}^{\prime},\bm{\theta})\right]$
		$\displaystyle\leq\left\\|K\right\\|_{\infty}\mathbb{E}_{\bm{\theta}^{\prime}\sim\rho}\left[\left\\|\nabla V(\bm{\theta}^{\prime})\right\\|+\frac{2}{\sigma}\left(\left\\|\bm{\theta}^{\prime}\right\\|+\left\\|\bm{\theta}\right\\|\right)\right]$
		$\displaystyle\leq\left\\|K\right\\|_{\infty}b_{1}+\mathbb{E}_{\bm{\theta}^{\prime}\sim\rho}\left[\left(\frac{2}{\sigma}+b_{1}\right)\left\\|\bm{\theta}^{\prime}\right\\|+\left\\|\bm{\theta}\right\\|\right].$

E.6.4 Proof of Lemma E.4

Given any $\bm{\theta}^{\prime}$ ,

		$\displaystyle\sum_{i=1}^{d}\left\\|\bar{\phi}_{\bm{\theta}^{\prime},i}(\bm{\theta})\right\\|_{\mathrm{Lip}}^{2}$
	$\displaystyle=$	$\displaystyle\sum_{i=1}^{d}\left[\underset{\bm{\theta}_{1}\neq\bm{\theta}_{2}}{\sup}\frac{\left\|\bar{\phi}_{\bm{\theta}^{\prime},i}(\bm{\theta}_{1})-\bar{\phi}_{\bm{\theta}^{\prime},i}(\bm{\theta}_{2})\right\|}{\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}}\right]^{2}$
	$\displaystyle=$	$\displaystyle\sum_{i=1}^{d}\underset{\bm{\theta}_{1}\neq\bm{\theta}_{2}}{\sup}\frac{\left\|\bar{\phi}_{\bm{\theta}^{\prime},i}(\bm{\theta}_{1})-\bar{\phi}_{\bm{\theta}^{\prime},i}(\bm{\theta}_{2})\right\|^{2}}{\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}^{2}}$
	$\displaystyle\leq$	$\displaystyle 2\sum_{i=1}^{d}\underset{\bm{\theta}_{1}\neq\bm{\theta}_{2}}{\sup}\frac{\left\|\left(e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{1}\right\\|^{2}/\sigma}-e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{2}\right\\|^{2}/\sigma}\right)\frac{\partial}{\partial\bm{\theta}_{i}^{\prime}}V(\bm{\theta}^{\prime})\right\|^{2}}{\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}^{2}}$
	$\displaystyle+$	$\displaystyle 2\sum_{i=1}^{d}\underset{\bm{\theta}_{1}\neq\bm{\theta}_{2}}{\sup}\frac{\left\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{1}\right\\|^{2}/\sigma}(\bm{\theta}_{1,i}-\bm{\theta}_{i}^{\prime})-\frac{2}{\sigma}e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{2}\right\\|^{2}/\sigma}(\bm{\theta}_{2,i}-\bm{\theta}_{i}^{\prime})\right\|^{2}}{\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}^{2}}.$

For the first term on the right side of the inequality,

		$\displaystyle\sum_{i=1}^{d}\underset{\bm{\theta}_{1}\neq\bm{\theta}_{2}}{\sup}\frac{\left\|\left(e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{1}\right\\|^{2}/\sigma}-e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{2}\right\\|^{2}/\sigma}\right)\frac{\partial}{\partial\theta_{i}^{\prime}}V(\bm{\theta}^{\prime})\right\|^{2}}{\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}^{2}}$
	$\displaystyle=$	$\displaystyle\sum_{i=1}^{d}\left\|\frac{\partial}{\partial\theta_{i}}V(\bm{\theta}^{\prime})\right\|^{2}\underset{\bm{\theta}_{1}\neq\bm{\theta}_{2}}{\sup}\frac{\left\|\left(e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{1}\right\\|^{2}/\sigma}-e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{2}\right\\|^{2}/\sigma}\right)\right\|^{2}}{\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}^{2}}$
	$\displaystyle=$	$\displaystyle\left\\|\nabla V(\bm{\theta}^{\prime})\right\\|^{2}\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathrm{Lip}}^{2}.$

To bound the second term, by the symmetry of each coordinates, we have

		$\displaystyle\sum_{i=1}^{d}\underset{\bm{\theta}_{1}\neq\bm{\theta}_{2}}{\sup}\frac{\left\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{1}\right\\|^{2}/\sigma}(\theta_{1,i}-\theta^{\prime}_{i})-\frac{2}{\sigma}e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}_{2}\right\\|^{2}/\sigma}(\bm{\theta}_{1,i}-\bm{\theta}^{\prime}_{i})\right\|^{2}}{\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}^{2}}$
	$\displaystyle=$	$\displaystyle d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathrm{Lip}}^{2}.$

This finishes the first part of the lemma.

		$\displaystyle\sum_{i=d}^{d}\left\\|\bar{\phi}_{\bm{\theta}^{\prime},i}(\bm{\theta})\right\\|_{\mathcal{L}_{\infty}}^{2}$
	$\displaystyle=$	$\displaystyle\sum_{i=d}^{d}\left\\|e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}\right\\|^{2}/\sigma}\left(\frac{2}{\sigma}\bm{\theta}_{i}-\frac{2}{\sigma}\bm{\theta}_{i}^{\prime}-\frac{\partial}{\partial\bm{\theta}_{i}^{\prime}}V(\bm{\theta}^{\prime})\right)\right\\|_{\mathcal{L}_{\infty}}^{2}$
	$\displaystyle\leq$	$\displaystyle\sum_{i=d}^{d}2\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}\right\\|^{2}/\sigma}\left(\bm{\theta}_{i}-\bm{\theta}_{i}^{\prime}\right)\right\\|_{\mathcal{L}_{\infty}}^{2}+\sum_{i=d}^{d}2\left\\|e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}\right\\|^{2}/\sigma}\frac{\partial}{\partial\bm{\theta}_{i}^{\prime}}V(\bm{\theta}^{\prime})\right\\|_{\mathcal{L}_{\infty}}^{2}$
	$\displaystyle\leq$	$\displaystyle\sum_{i=d}^{d}2\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}^{\prime}-\bm{\theta}\right\\|^{2}/\sigma}\left(\bm{\theta}_{i}-\bm{\theta}_{i}^{\prime}\right)\right\\|_{\mathcal{L}_{\infty}}^{2}+2\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathcal{L}_{\infty}}^{2}\left\\|\nabla V(\bm{\theta}^{\prime})\right\\|^{2}$
	$\displaystyle\leq$	$\displaystyle 2d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathcal{L}_{\infty}}^{2}+2\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathcal{L}_{\infty}}^{2}\left\\|\nabla V(\bm{\theta}^{\prime})\right\\|^{2}.$

	$\displaystyle\left\\|K(\bm{\theta}^{\prime},\bm{\theta}_{1})-K(\bm{\theta}^{\prime},\bm{\theta}_{2})\right\\|$	$\displaystyle\leq$	$\displaystyle\left\\|e^{-(\cdot)^{2}/\sigma}\right\\|_{\mathrm{Lip}}\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}$
	$\displaystyle\left\\|\nabla_{\bm{\theta}^{\prime}}K(\bm{\theta}^{\prime},\bm{\theta}_{1})-\nabla_{\bm{\theta}^{\prime}}K(\bm{\theta}^{\prime},\bm{\theta}_{2})\right\\|$	$\displaystyle\leq$	$\displaystyle\left\\|\frac{2}{\sigma}e^{-(\cdot)^{2}/\sigma}(\cdot)\right\\|_{\mathrm{Lip}}\left\\|\bm{\theta}_{1}-\bm{\theta}_{2}\right\\|_{2}.$

	$\displaystyle\left\\|\phi[\rho](\cdot)\right\\|_{\mathrm{Lip}}$	$\displaystyle\leq\left\\|e^{-(\cdot)^{2}/\sigma}\right\\|_{\mathrm{Lip}}\mathbb{E}_{\bm{\theta}\sim\rho}\left\\|\nabla V(\bm{\theta})\right\\|+\left\\|\frac{2}{\sigma}e^{-(\cdot)^{2}/\sigma}(\cdot)\right\\|_{\mathrm{Lip}},$
	$\displaystyle\left\\|\phi[\rho](\bm{\theta})\right\\|$	$\displaystyle\leq\left\\|K\right\\|_{\infty}\mathbb{E}_{\bm{\theta}^{\prime}\sim\rho}\left[\left\\|\nabla V(\bm{\theta}^{\prime})\right\\|+\frac{2}{\sigma}\left(\left\\|\bm{\theta}^{\prime}\right\\|+\left\\|\bm{\theta}\right\\|\right)\right]$
		$\displaystyle\leq\left\\|K\right\\|_{\infty}b_{1}+\mathbb{E}_{\bm{\theta}^{\prime}\sim\rho}\left[\left(\frac{2}{\sigma}+b_{1}\right)\left\\|\bm{\theta}^{\prime}\right\\|\right]+\left\\|\bm{\theta}\right\\|.$

	$\displaystyle\sum_{i=1}^{d}\left\\|\bar{\phi}_{\bm{\theta}^{\prime},i}(\bm{\theta})\right\\|_{\mathrm{Lip}}^{2}$	$\displaystyle\leq$	$\displaystyle 2\left\\|\nabla V(\bm{\theta}^{\prime})\right\\|^{2}\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathrm{Lip}}^{2}+2d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathrm{Lip}}^{2}$
	$\displaystyle\sum_{i=d}^{d}\left\\|\bar{\phi}_{\bm{\theta}^{\prime},i}(\bm{\theta})\right\\|_{\mathcal{L}_{\infty}}^{2}$	$\displaystyle\leq$	$\displaystyle 2d\left\\|\frac{2}{\sigma}e^{-\left\\|\bm{\theta}\right\\|^{2}/\sigma}\theta_{1}\right\\|_{\mathcal{L}_{\infty}}^{2}+2\left\\|e^{-\left\\|\cdot\right\\|^{2}/\sigma}\right\\|_{\mathcal{L}_{\infty}}^{2}\left\\|\nabla V(\bm{\theta}^{\prime})\right\\|^{2}.$

	$\displaystyle\left\\|\bm{\psi}\left(\bm{\Theta}\right)\right\\|$	$\displaystyle=\left\\|\bm{\theta}_{1}+\eta\left\{-\nabla V(\bm{\theta}_{1})+\alpha\phi[\frac{1}{M}\sum_{j=1}^{M}\delta_{\bm{\theta}_{1+jC}}](\bm{\theta}_{k})\right\}\right\\|$
		$\displaystyle=\left\\|\bm{\theta}_{1}-\eta\nabla V(\bm{\theta}_{1})+\frac{\eta\alpha}{M}\sum_{j=1}^{M}\left[e^{-\left\\|\bm{\theta}_{1+jC}-\bm{\theta}_{1}\right\\|^{2}/\sigma}\frac{2}{\sigma}\left(\bm{\theta}_{1}-\bm{\theta}_{1+jC}\right)-e^{-\left\\|\bm{\theta}_{1+jC}-\bm{\theta}_{1}\right\\|^{2}/\sigma}\nabla V(\bm{\theta}_{1+jC})\right]\right\\|$
		$\displaystyle\leq\left\\|\bm{\theta}_{1}-\eta\nabla V(\bm{\theta}_{1})+\frac{2}{\sigma}\frac{\eta\alpha}{M}\sum_{j=1}^{M}e^{-\left\\|\bm{\theta}_{1+jC}-\bm{\theta}_{1}\right\\|^{2}/\sigma}\bm{\theta}_{1}\right\\|$
		$\displaystyle+\left\\|\frac{\eta\alpha}{M}\sum_{j=1}^{M}e^{-\left\\|\bm{\theta}_{1+jC}-\bm{\theta}_{1}\right\\|^{2}/\sigma}\frac{2}{\sigma}\left(-\nabla V(\bm{\theta}_{1+jC})-\bm{\theta}_{1+jC}\right)\right\\|$
		$\displaystyle\leq\left\\|\bm{\theta}_{1}-\eta\nabla V(\bm{\theta}_{1})\right\\|+\frac{2\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty}b_{1}(1+\left\\|\bm{\theta}_{1}\right\\|)$
		$\displaystyle+\frac{2\alpha\eta}{M\sigma}\sum_{j=1}^{M}\left\\|K\right\\|_{\infty,\infty}b_{1}\left(1+(1+\frac{1}{b_{1}})\left\\|\bm{\theta}_{1+jC}\right\\|\right)$
		$\displaystyle\overset{(1)}{\leq}b_{1}(1+\frac{4\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty})+\eta^{2}b_{1}+2\eta b_{1}$
		$\displaystyle+\left(1-2\eta a_{1}+\eta^{2}b_{1}+\frac{2\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty}b_{1}\right)\left\\|\bm{\theta}_{1}\right\\|+\frac{2\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty}(b_{1}+1)\max_{i\in[MC+1]-\{1\}}\left\\|\bm{\theta}_{1+jC}\right\\|$
		$\displaystyle\leq b_{1}(1+\frac{4\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty})+\eta^{2}b_{1}+2\eta b_{1}$
		$\displaystyle+\max\left(1-2\eta a_{1}+\eta^{2}b_{1}+\frac{2\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty}b_{1},\frac{2\alpha\eta}{\sigma}\left\\|K\right\\|_{\infty,\infty}(b_{1}+1)\right)\max_{i\in[MC+1]}\left\\|\bm{\theta}_{1+jC}\right\\|,$

	$\displaystyle I_{1}+I_{3}$	$\displaystyle\leq\frac{3Cc_{d}^{2}}{4}\sum_{j=0}^{k-1}\int_{t_{j}}^{t_{j+1}}\mathbb{E}\left\\|\hat{\bm{\theta}}_{t_{j}}-\hat{\bm{\theta}}_{s}\right\\|^{2}ds$
		$\displaystyle\leq\frac{3Cc_{d}^{2}}{4}\int_{t_{j}}^{t_{j+1}}\left[\mathbb{E}\left\\|-\nabla V(\hat{\bm{\theta}}_{\underline{s}})+\alpha\phi[\pi_{M,c}*\hat{\rho}_{\underline{s}}](\hat{\bm{\theta}}_{\underline{s}})\right\\|^{2}(s-t_{j})^{2}+2d(s-t_{j})\right]ds$
		$\displaystyle=Cc_{d}^{2}\eta^{3}\sum_{j=0}^{k-1}\mathbb{E}\left\\|-\nabla V(\hat{\bm{\theta}}_{t_{j}})+\alpha\phi[\pi_{M,c}*\hat{\rho}_{t_{j}}](\hat{\bm{\theta}}_{t_{j}})\right\\|^{2}+Cc_{d}^{2}dkh^{2}.$

Stein Self-Repulsive Dynamics: Benefits from Past Samples

Abstract

1 Introduction

2 Background: Langevin dynamics & SVGD

3 Stein Self-Repulsive Dynamics

4 Theoretical Analysis of Stein Self-Repulsive Dynamics

4.1 Mean-Field and Continuous-Time Limits

4.2 Assumptions

Assumption 4.1 (RBF Kernel).

Assumption 4.2 (VV is dissipative and smooth).

Assumption 4.3 (Regularity Condition).

Assumption 4.4 (Strong-convexity).

4.3 Main Theorems

Theorem 4.1 (Stationary Distribution).

Theorem 4.2 (Time Discretization Error).

Theorem 4.3 (Mean-Field Limit).

Corollary 4.1 (Convergence to Target Distribution).

5 Extension to General Dynamics

Theorem 5.1 (Stationary Distribution).

6 Experiments

6.1 Synthetic Experiment

6.2 Bayesian Neural Network

6.3 Contextual Bandits

7 Conclusion

Broader Impact Statement

Acknowledgement

References

Appendix A Discussion on Hyper-parameter Tuning

Appendix B Additional Experiment Result on Synthetic Data

B.1 Synthetic 2D Mixture of Gaussian Experiment

B.2 Synthetic higher dimensional Gaussian Experiment

B.3 Synthetic higher dimensional Mixture of Gaussian Experiment

Appendix C BNN on UCI Datasets: Experiment Settings and Additional Results

Appendix D Contextual Bandit: Experiment Settings and More Background

Appendix E The Detailed analysis of SRLD

E.1 Some additional notation

E.2 Geometric Ergodicity of SRLD

Theorem E.1 (Geometric Ergodicity).

E.3 Moment Bound

Theorem E.2 (Moment Bound).

E.4 Technical Lemma

Definition E.1 (α\alpha-mixing).

Definition E.2 (β\beta-mixing).

Proposition E.1 (β\beta-mixing implies α\alpha-mixing)).

Proposition E.2.

Lemma E.1 (Regularity Conditions).

Lemma E.2 (Properties of RBF Kernel).

Lemma E.3 (Properties of Stein Operator).

Lemma E.4 (Bounded Lipschitz of Stein Operator).

E.5 Proof of Main Theorems

E.5.1 Proof of Theorem E.1

E.5.2 Proof of Theorem E.2

E.5.3 Proof of Theorem 4.1 and 5.1

E.5.4 Proof of Theorem 4.2

E.5.5 Proof of Theorem 4.3

E.6 Proof of Technical Lemmas

E.6.1 Proof of Lemma E.1

E.6.2 Proof of Lemma E.2

E.6.3 Proof of Lemma E.3

E.6.4 Proof of Lemma E.4

Assumption 4.2 ( $V$ is dissipative and smooth).

Definition E.1 ( $\alpha$ -mixing).

Definition E.2 ( $\beta$ -mixing).

Proposition E.1 ( $\beta$ -mixing implies $\alpha$ -mixing)).