Title: STEIN SELF-REPULSIVE DYNAMICS: BENEFITS FROM PAST SAMPLES Anonymous authors Paper under double-blind review

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

MCMC works by simulating Markov chains whose stationary distributions match the distributions of interest. Despite nice asymptotic theoretical properties, MCMC is widely criticized for its slow convergence rate in practice. In difficult problems, the samples drawn from MCMC are often found to have high auto-correlation across time, meaning that the Markov chains explore very slowly in the configuration space. When this happens, the samples returned by MCMC only approximate a small local region, and under-estimate the probability of the regions un-explored by the chain.

Stein variational gradient descent (SVGD) (Liu & Wang, 2016) is another type of approximation sampling methods designed to overcome the limitation of MCMC. Instead of drawing random samples sequentially, SVGD evolves a pre-defined number of particles (or sample points) in parallel with a special interacting particle system to match the distribution of interest by minimizing the KL divergence. In SVGD, the particles interact with each other to simultaneously move towards the high probability regions following the gradient direction, and also move away from each other due to a special repulsive force. As a result, SVGD allows us to obtain diversified samples that correctly represent the variation of the distribution of interest. SVGD has been found a promising tool for solving difficult sampling problems in which diversity promotion is critical (e.g., Feng et al., 2017; Haarnoja et al., 2017; Pu et al., 2017; Liu et al., 2017; Gong et al., 2019) . Various extensions have been developed (e.g., Han & Liu, 2018; Chen et al., 2018; Liu et al., 2019; Wang et al., 2019a) .

However, one problem of SVGD is that it theoretically requires to run an infinite number of chains in parallel in order to approximate the target distribution asymptotically (Liu, 2017) . With a finite number of particles, the fixed point of SVGD does still provide a prioritized, partial approximation to the distribution in terms of the expectation of a special case of functions (Liu & Wang, 2018) . Nevertheless, it is still desirable to develop a variant of "single-chain SVGD", which only requires to run a single chain sequentially like MCMC to achieve the correct stationary distribution asymptotically in time, with no need to take the limit of infinite number of parallel particles.

In this work, we propose an example of single-chain SVGD by integrating the special repulsive mechanism of SVGD with gradient-based MCMC such as Langevin dynamics. Our idea is to use repulsive term of SVGD to enforce the samples in MCMC away from the past samples visited at previous iterations. Such a new self-repulsive dynamics allows us to decrease the auto-correlation in MCMC and hence increase the mixing rate, but still obtain the same stationary distribution thanks to the special property of the SVGD repulsive mechanism.

We provide throughout theoretically analysis of our new method, establishing it asymptotic convergence to the target distribution. As we show in the work, the analysis is highly non-trivial, because our new self-repulsive dynamics is a non-linear, high-order Markov process. Empirically, we extensively evaluate our methods on an array of challenging sampling tasks, showing that our method yields much better uncertainty estimation and larger effective sample size.

2. BACKGROUND: LANGEVIN DYNAMICS AND SVGD:

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

Langevin dynamics Langevin dynamics is a basic gradient based MCMC method. For some target distribution on R d with density function ρ * (θ) ∝ exp(-V (θ)), where V : R d → R is the potential function, the (Euler-discrerized) Langevin dynamics simulates a Markov chain with the following rule:

θ k+1 = θ k -η∇V (θ k ) + 2ηe k , e k ∼ N (0, I),

where k denotes the number of iterations, {e k } are independent standard Gaussian noise, and η is a step size parameter. It is well known that the limiting distribution of θ k when k → ∞ approximates the target distribution when η is sufficiently small.

Because the updates in Langevin dynamics are local and incremental, new points generated by the dynamics is highly correlated to the past sample. As a result, we need to run Langevin dynamics sufficiently long in order to obtain diverse samples.

Stein Variatinal Gradient Descent (SVGD) Different from Langevin dynamics, SVGD iteratively evolves a pre-defined number of particles in parallel. Starting from an initial set of particles {θ i 0 : i = 1, ..., M }, SVGD updates the M particles in parallel by

θ i k+1 = θ i k + ηg(θ i k ; δM k ) ∀i = 1, . . . , M,

where the velocity field, which we denote by g(θ i k ; δM k ), depends the empirical distribution of the current set of particles δM k := 1 M M j=1 δ θ j k in the following way,

g(θ i k ; δM k ) = E θ∼ δM k   -K(θ, θ i k )∇V (θ) Confining Term + ∇ θ K(θ, θ i k ) Repulsive Term    .

Here δ θ is the Dirac measure centered at θ, and hence E θ∼ δM k [•] denotes averaging on the particles. The K(•, •) is a positive definite kernel, such as the Gaussian RBF kernel, specified by users. Note that g(θ i k ; δM k ) consists of a confining term and repulsive term: the confining term pushes particles to move towards high density region, and the repulsive term prevents the particles from colliding with each other. It is the balance of these two terms that allows us to asymptotically approximate the target distribution ρ * (θ) ∝ exp(-V (θ)) at the fixed point, when the number of particles goes to infinite. We refer the readers to Liu & Wang (2016) ; Liu (2017) ; Liu & Wang (2018) for throughout theoretical justifications of SVGD. But a quick, informal way to justify the SVGD update is through the Stein's identity, which shows that the velocity field g(θ; ρ) equals zero when ρ equals the true distribution ρ * , that is,

g(θ ; ρ * ) = E θ∼ρ * [-K(θ, θ )∇V (θ) + ∇ θ K(θ, θ )] = 0, ∀θ ∈ R d .

(1)

This shows that SVGD would converge if the particle distribution already forms a closed approximation to the target distribution p, meaning that the target distributions forms an (approximate) fixed point of the update.

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 sample and the past samples, hence forming a new self-avoiding dynamics with fast convergence speed.

Specifically, assume we run a single Markov chain like Langevin dynamics, where θ k denotes the particle at the k-th iteration. Denote by δM k the empirical distribution of M samples taken from the past iterations, i.e.,

δM k := 1 M M j=1 δ θ k-jcη , c η = c/η,

where c η is a thinning factor, which scales inversely with the step size η, introduced to slim the sequence of past samples. Compared with the δM k in SVGD, which is averaged over M parallel particles, δM k is averaged across time over M past samples. Given this, our Stein self-repulsive dynamics updates the sample via

EQUATION

in which the particle is updated with the typical Langevin gradient, plus a Stein repulsive force against the samples from the previous iterations. α ≥ 0 is a parameter that controls the magnitude of the Stein repulsive term. In this way, the particles are pushed away from the past samples, and hence admits lower auto-correlation and faster convergence speed. Importantly, the addition of the repulsive force does not impact the asymptotic stationary distribution, thanks to Stein's identity in (1). This is because when if the self-repulsive dynamics have converged to the target distribution ρ * , such that θ k ∼ ρ * for all k, the Stein self-repulsive term would equal to zero in expectation due to Stein's identity and hence does not introduce additional bias over Langevin dynamics. Rigorous theoretical analysis of this idea is developed in Section 4.

Practical Algorithm Because δM k 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,

EQUATION

It includes two phases. The first phase is the same as the Langevin dynamics which collects the initial M samples used in the second phase while serves as a warm start. The repulsive gradient update is introduced in the second phase to encourage the dynamics to visit the under-explored density region. We call this particular instance of our algorithm Self-Repulsive Langevin dynamics (SRLD), self-repulsive variants of more general dynamics is discussed in Section 5.

Remark Notice that, the first phase is introduced to collect the initial M samples. However, it's not really necessary to generate the initial M samples with Langevin dynamics. We can simply use some other initialization distribution and get M initial samples from that distribution. In practice, we find using Langevin dynamics to collect the initial samples is natural and it can also be viewed as the burn-in phase before sampling, so we use (3) in all of the other experiments.

Remark:

The general idea of introducing self-repulsive terms inside MCMC or other iterative algorithms is not new itself. For example, in molecular dynamics simulations, an algorithm called metadynamics (Laio & Parrinello, 2002) has been widely used, in which the particles are repelled away from the past samples in a way similar to our method, but with a typical repulsive function, such as j D(θ k , θ k-jcη ), where D(•, •) is any notation of dis-similarity. However, introducing 

Initial Sample:

Repulsive Sample 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.

2 1 0 1 2 2 1 0 1 2

Initial Sample Langevin Sample:

an arbitrary repulsive force would alter the stationary distribution of the dynamics, introducing a harmful bias into the algorithm. The key highlight of our approach, as reflected by our theoretical results in Section 4, is the unique property of the Stein repulsive term, that allows us to obtain the correct stationary distribution even with the addition of the repulsive force.

Remark Recently, (Zhang et al., 2018) proposed a different combination of SVGD and Langevin dynamics, in which the Langevin force is directly added to a set of particles that evolve in parallel with SVGD. Using our terminology, their system is

θ i k+1 = θ i k + (-ηV (θ i k ) + 2ηe i k ) + ηαg(θ i k ; δM k ), e k ∼ N (0, I) ∀i = 1, . . . , M.

This is significantly different from our method on both motivation and practical algorithm. Their algorithm still requires to simulate M chains of particles in parallel like SVGD, and was proposed to obtain easier theoretical analysis than the deterministic dynamics of SVGD. Our work is instead motivated by the practical need of decreasing the auto-correlation in Langevin dynamics, and avoiding the need of running multiple chains in SVGD, and hence must be based on self-repulsion against past samples along a single chain.

In another recent work, (Chen et al., 2018) proposed a π-SGLD method, which simulates the linear combination of the evalutionary partial differential equations of SVGD and Langevin dynamics using discrete gradient flow with blob-based method. Their method is again motivated by the theoretical interest of discovering new categories of algorithms, and does not involve self-repulsive on a single chain like our method.

An Illustrative Example Here we give an illustrative example to show the key advantage of our self-repulsive dynamics. Assume that we want to sample from a bi-variate Normal distribution shown in Figure 1 . Unlike standard settings, we assume that we have already obtained some initial samples (yellow dots in Figure 1 ) before running the dynamics. The initial samples are assumed to concentrate on the left part of the target distribution as shown in Figure 1 . In this extreme case, since the left part of the distribution is over-explored by the initial samples, it is desirable to have the subsequent new samples to concentrate more on the un-explored part of the distribution. However, standard Langevin dynamics does not take this into account, and hence yielding a bias overall representation of the true distribution (see the left panel). With our self-repulsive dynamics, the new samples are forced to explore the un-explored region more frequently, allowing us to obtain a much more accurate approximation when combining the new and initial samples.

4. THEORETICAL ANALYSIS OF STEIN SELF-REPULSIVE DYNAMICS:

In this section, we provide theoretical analysis of the self-repulsive dynamics. We establish that our self-repulsive dynamics converges to the correct target distribution asymptotically, in the limit when M approaches to infinite and the step size η approaches to 0. This is a highly non-trivial task, because the self-repulsive dynamics is a highly complex, non-linear and high order Markov stochastic process. We attack this problem by breaking the proof into the following three steps illustrated in Figure 2 :

1) At the limit of large particle sizes M → ∞ (called the mean field limit), we show that practical dynamics in (3) is closely approximated by a discrete-time mean-field dynamics characterized by (4) below.

2) By further taking the limit of small step size η → 0 + (called the continuous-time limit), the dynamics in (4) converges to a continuous-time mean-field dynamics characterized by (5) below.

3) We show that the dynamics in (5) converges to the target distribution.

Remark As we mentioned in Section 3, we introduce the first phase to collect the initial M samples for the second phase, and our theoretical analysis follows this setting to make our theory as close to the practice. However, the theoretical analysis can be generalized to the setting of drawing M initial samples from some initialization distribution with almost identical argument.

Notations We use • and •, • to represent the 2 vector norm and inner product respectively. The Lipschitz norm and bounded Lipschitz norm of a function f are defined by f Lip and f BL . The KL divergence, Wasserstein-2 distance and Bounded Lipschitz distance between distribution ρ 1 , ρ 2 are denoted as

D KL [ρ 1 ρ 2 ], W 2 [ρ 1 , ρ 2 ] and D BL [ρ 1 , ρ 2 ] respectively.

4.1. MEAN-FIELD AND CONTINUOUS-TIME LIMITS:

To fix the notation, we denote by ρ k := Law(θ k ) the distribution of θ k at time k of the practical self-repulsive dynamics (3), which we refer as the practical dynamics in the sequel, when the initial particle θ 0 is drawn from an initial continuous distribution ρ 0 supported on R d . Note that given ρ 0 , the subsequent ρ k can be recursively defined through dynamics (3). Due to the diffusion noise in Langevin dynamics, all ρ k are continuous distributions supported on R d . We now introduce the limit dynamics when we take the mean-field limit (M → +∞) and then the continuous-time limit (η → 0 + ).

Discrete-Time Mean-Field Dynamics (M → +∞) In the limit of M → ∞, we show that our practical dynamics is closely approximated by the following dynamics, in which the empirical measures are replaced by the continuous distributions of the particles themselves.

EQUATION

where ρM k = 1 M M j=1 ρk-jcη by ρk := Law( θk ) the (smooth) distribution of θk at time-step k when the dynamics is initialized with θ0 ∼ ρ0 = ρ 0 . Compared with the practical dynamics in (3), the difference is that the empirical distribution δM k is replaced by the smooth distribution ρM k . Similar to the recursive definition of ρ k following dynamics (3), ρk is also recursively defined through dynamics (4), starting from ρ0 = ρ 0 . As we show in Theorem 4.3, if the auto-correlation of θ k decays fast enough and M is sufficiently large, ρM k is well approximated by the empirical distribution δM k in (3), and further the two dynamics ((3) and ( 4)) converges to each other in the sense that W 2 [ρ k , ρk ] → 0 as M → ∞ for any k. Note that in taking the limit of M , we need to ensure that we run the dynamics for more than M c η steps. Otherwise, SRLD degenerates to Langevin dynamics as we stop the chain before we collect M samples.

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

EQUATION

where ρM t := 1 M M j=1 ρt-jc (•) and ρt = Law θt is the distribution of θt at a continuous time point t with θ 0 initialized by θ0 ∼ ρ0 = ρ 0 . We prove that ( 5) is closely approximated by ( 4) with small step size in the sense that D KL [ρ k ρkη ] → 0 as η → 0 in Theorem 4.2, and importantly, the stationary distribution of (5) equals to the target distribution ρ * (θ) ∝ exp(-V (θ)).

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(θ 1 , θ 2 ) = exp(-θ 1 -θ 2 2 /σ) for some fixed 0 < σ < ∞.

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

Assumption 4.2. (V is dissipative and smooth)

Assume that θ, -∇V (θ) ≤ b 1 -a 1 θ 2 and ∇V (θ 1 ) -∇V (θ 2 ) ≤ b 1 θ 1 -θ 2 .

We also assume that ∇V (0) ≤ b 1 . Here a 1 and b 1 are some finite positive constant.

Assumption 4.3. (Regularity Condition) Assume E θ∼ρ0 [ θ 2 ] < 0. Define ρ M k = M j=1 ρ k-jcη /M , assume there exists B < ∞ such that inf k≥M cη sup θ ≤B E g(θ; δM k ) -g(θ; ρ M k ) 2 > 0. Assumption 4.4. (Strong-convexity) Suppose that ∇V (θ 1 ) -∇V (θ 2 ), θ 1 -θ 2 ≥ L θ 1 -θ 2 2 for a positive constant L.

Remark Assumption 4.2 is standard in the existing Langevin dynamics analysis (see Dalalyan (2017); Raginsky et al. (2017) for example). Assumption 4.3 is a weak condition as it assumes that the dynamics can not degenerate into one local mode and stop moving anymore. This assumption is expected to be true when we have diffusion terms like the Gaussian noises in our self-repulsive dynamics. Assumption 4.4 is a classical assumption on the existing Langevin dynamics analysis with convex potential Dalalyan (2017); Durmus et al. (2019) . Although being a bit strong, this assumption broadly applies to posterior inference problem in the limit of big data, as the posterior distribution converges to Gaussian distributions for large training set as shown by Bernstein-von Mises theorem. It is technically possible to further generalize our results to the non-convex settings with a refined analysis, which we leave as future work. This work focuses on the classic convex setting for simplicity.

4.3. MAIN THEOREMS:

All of the proof in this section can be found in Appendix B.5.

We first prove that the limiting distribution of the continuous-time mean field dynamics ( 5) is the target distribution. This is achieved by writing dynamics (5) into the following (non-linear) partial differential equation:

∂ t ρt = ∇ • (-∇V ρt ) + ∆ρ t t ∈ [0, M c) ∇ • -∇V + αg(•, ρM t ) ρt + ∆ρ t , t ≥ M c. Theorem 4.1. (Stationary Distribution)

Given some finite M , c and α, and suppose that the limiting distribution of dynamics (5) exists. Then the limit distribution is unique, and equals to ρ * (θ) ∝ exp(-V (θ)).

We then give the upper bound on the discretization error, which can be characterized by analyzing the KL divergence between ρk and ρkη . Theorem 4.2. (Time Discretization Error)

Given some sufficiently small step size η and choose α < a 2 /(2b 1 + 4/σ). Under assumption 4.1, 4.2, 4.3 and c η = c/η. we have for some constant C, max l∈{0,...,k}

D KL [ρ lη ρl ] ≤ O η + kη 2 k ≤ M c η -1 O η + M cη + α 2 M ce Cα 2 (kη-M c) η 2 k ≥ M c η .

With this theorem, we can know that if η is small enough, then the discretization error is small and ρ approximates ρ closely. Next we give result on the mean field limit of M → ∞. Theorem 4.3. (Mean-Field Limit)

Under assumption 4.1, 4.2, 4.3, and 4.4, suppose that we choose α and η such that

-(a 1 -2αb 1 /σ)+ ηb 1 < 0; 2αη σ (b 1 + 1) < 1; a 2 -α 2b 1 + 4 σ > 0;

Then there exists a constant c 2 , such that when L/a ≥ c 2 and we have

W 2 2 [ρ k , ρk ] = O α 2 /M + η 2 for k ≥ M c η ,

and

W 2 2 [ρ k , ρk ] = 0 if k ≤ M c η -1.

Remark Combine all the result, we have lim k,M →∞,η→0 + D BL [ρ k , ρ * ] = 0, where the joint limit requires kη → ∞, exp(Cα 2 kη)η 2 = o(1) and kη M c = γ (1 + o(1)), with γ > 1. Notice that if γ ≤ 1, the dynamics is degenerated to Langevin dynamics.

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(θ) ∝ exp(-V (θ)):

dθ t = -f (θ)dt + 2D(θ)dB t , with f (θ) = [D(θ) + Q(θ)]∇V (θ) -Γ(θ), Γ i (θ) = d j=1 ∂ ∂θ j (D ij (θ) + Q ij (θ)) .

where D is a positive semi-definite diffusion matrix that determines the strength of the Brownian motion and Q is a skew-symmetric curl matrix that can represent the traverse effect (e.g. in Neal et al., 2011; Ding et al., 2014) . By adding the Stein repulsive force, we obtain the following general self-repulsive dynamics where ρt := Law( θt ) is again the distribution of θt following ( 6) when initalized at θ0 ∼ ρ0 . Similar to the case of Langevin dynamics, this process also converges to the correct target distribution, and can be simulated by practical dynamics similar to (3). Theorem 5.1. (Stationary Distribution)

EQUATION

Given some finite M , c and α, and suppose that the limiting distribution of dynamics (6) exists. Then the limiting distribution is unique and equals the target distribution ρ * (θ) ∝ exp(-V (θ)).

6. EXPERIMENTS:

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

6.1. SYNTHETIC EXPERIMENT:

A Correlated 2D Distribution This experiment aims to show how the repulsive gradient helps explore the whole distribution. Following Ma et al. (2015) , we compare the sampling efficiency on the following correlated 2D distribution with density

ρ * ([θ 1 , θ 2 ]) ∝ -θ 4 1 /10 -4 (θ 2 + 1.2) -θ 2 1 2 /2.

We compare the SRLD with vanilla Langevin dynamics, and evaluate the sample quality by Maximum Mean Discrepancy (MMD) (Gretton et al., 2012) , Wasserstein-1 Distance and effective sample size (ESS). Notice that the finite sample quality of gradient based MCMC method is highly related to the step size. Compared with Langevin dynamics, we have an extra repulsive gradient and thus we implicitly have larger step size. To rule out this effect, we set different step sizes of the two dynamics so that the gradient of the two dynamics has the same magnitude. In addition, to decrease the influence of random noise, the two dynamics are set to have the same initialization and use the same sequence of Gaussian noise. We collect the sample of every iteration. We repeat the experiment 20 times with different initialization and sequence of Gaussian noise.

Figure 3 summarizes the result with different metrics. We can see that SRLD have a significantly smaller MMD and Wasserstein-1 Distance as well as a larger ESS compared with the vanilla Langevin dynamics. Moreover, the introduced repulsive gradient creates a negative auto-correlation between samples. Figure 4 shows a typical trajectory of the two sampling dynamics. We can see that SRLD have a faster mixing rate than vanilla Langevin dynamics. Note that since we use the same sequence of Gaussian noise for both algorithms, the difference is mainly due to the use of repulsive gradient rather than the randomness.

Mixture of Gaussian Distribution:

We aim to show how the repulsive gradient helps the particle escape from the local high density region by sampling the 2D mixture of Gaussian distribution using SRLD and Langevin dynamics. The target density is set to be where 1 = [1, 1] . This target distribution have two mode at -1 and 1, and vanilla Langevin dynamics can stuck in one mode while keeps the another mode under-explored (as the gradient of energy function can dominate the update of samples). We use the same evaluation method, step sizes, initialization and Gaussian noise as the previous experiment. We collect one sample every 100 iterations and the experiment is repeated for 20 times. Figure 5 shows that SRLD again consistently outperforms the Langevin dynamics on all of the evaluation metrics.

ρ * (θ) ∝ 1 2 exp -θ -1 2 /2 + 1 2 exp -θ + 1 2 /2 , θ = [θ 1 , θ 2 ] ,

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 & Graff, 2017) . We assume the output is normal distributed, with a two-layer neural network with 50 hidden units and tanh activation to predict the mean of outputs. We set a Γ(1, 0.1) prior for the inverse output variance. All of the datasets are randomly partitioned into 90% for training and 10% for testing. The results are averaged over 20 random trials. We set the mini-batch size to be 100 and the number of past samples M to be 10.

In all experiments, we use RBF kernel with bandwidth set by the median trick as suggested in Liu & Wang (2016) . For SVGD, we use the original implementation with 20 particles by Liu & Wang (2016) . We run 50000 iterations for each methods, and for LD and SRLD, the first 40000 iteration is discarded as burn-in. We use a thinning factor of c η = c/η = 100 and in total we collect 100 samples from the posterior distribution. For each dataset, we generate 3 extra data splits for tuning the step size for each method.

Table 1 shows the average test RMSE and test log-likelihood and their standard deviation. The method that has the best average performance is marked as boldface. We observe that a large portion of the variance is due to the random partition of the dataset. Therefore, to show the statistical significance, we use the matched pair t-test to test the statistical significance, mark the methods that perform the best with a significance level of 0.05 with underlines. Note that the results of SRLD/LD and SVGD is not very comparable, because SRLD/LD are single chain methods which averages across time, and SVGD is a multi-chain method that only use the results of the last iteration. We provide additional results in appendix that SRLD averaged on 20 particles (across time) can also achieve similar or better results as SVGD with 20 (parallel) particles.

6.3. CONTEXTUAL BANDITS:

We evaluate the quality of uncertainty estimation provided by our methods on several contextual bandits problems. Uncertainty estimation is a key component of contextual bandits. If the agent makes decisions with a poorly estimated uncertainty, the decisions will finally lead to catastrophic failure through the feedback loops (Riquelme et al., 2018) .

Though in principle all of the MCMC methods return the samples follow the true posterior if we can run infinite MCMC steps, in practice we can only obtain finite samples as we only have finite time to run the MCMC sampler. In this case, the auto-correlation issue can lead to the under-estimate the uncertainty, which will cause the failure on all of the reinforcement learning problems that need exploration.

We consider the posterior sampling (a.k.a Thompson sampling) algorithm with Bayesian neural network as the function approximator. We follows the experimental setting from Riquelme et al. (2018) . The only difference is that we change the optimization of the objective (e.g. evidence lower bound (ELBO) in variational inference methods) into running MCMC samplers. We set the step of samplers equal to the optimization step, and use a thinning factor of 100. We compare the SRLD with the Langevin dynamics on the Mushroom and Wheel bandits from (Riquelme et al., 2018) , and include SVGD as a baseline. For more detailed introduction and setup, see Appendix A.5.

The cumulative regret is shown in Table 2 . SVGD is known to have the under-estimated uncertainty for Bayesian neural network if particle number is limited (Wang et al., 2019b) , and as a result, has the worst performance among the three methods. SRLD is slightly better than vanilla Langevin dynamics on the simple Mushroom bandits. On the much more harder Wheel bandits, SRLD is significantly better than the vanilla Langevin dynamics, which shows the improving uncertainty estimation of our methods within finite number of samples. 

Dataset:



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.

A EXPERIMENT DETAILS AND ADDITIONAL EXPERIMENT RESULT:

A.1 SYNTHETIC 2D MIXTURE OF GAUSSIAN EXPERIMENT

3 2 1 0 1 2 3 2 1 0 1 2 3

Repulsive Dynamics

3 2 1 0 1 2 3 2 1 0 1 2 3 3 2 1 0 1 2 3 2 1 0 1 2 3 3 2 1 0 1 2 3 2 1 0 1 2 3 3 2 1 0 1 2 3 2 1 0 1 2 3 3 2 1 0 1 2

Num of samples = 100 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.

A.2 SYNTHETIC HIGHER DIMENSIONAL GAUSSIAN EXPERIMENT:

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

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

ρ * (θ) ∝ 1 2 exp -0.5 θ -2/d1 2 + 1 2 exp -0.5 θ + 2/d1 2 ,

where θ = [θ 1 , ..., θ 20 ] and 1 = [1, ..., 1] . And thus the mean of the two mixture component is with distance 2 √ 2. We run SRLD with α = 20, 10, 5, 0 (and when α = 0, it reduces to LD). The other experiment setting is the same as the low dimensional mixture Gaussian case. Figure [8] summarizes the result. As shown in the figure, when α becomes larger, the repulsive forces helps the sampler better explore the density region.

A.4 UCI DATASETS:

We show some additional experiment result on posterior inference on UCI datasets. As mentioned in Section 6.2, the comparison between SVGD and SRLD is not direct as SVGD is a multiple-chain method with fewer particles and SRLD is a single chain method with more samples. To show more detailed comparison, we compare the SVGD with SRLD using the first 20, 40, 60, 80 and 100 samples, denoted as SRLD-n where n is the number of samples used. Table 3 shows the result of averaged test RMSE and table 4 shows the result of averaged test loglikelihood. For SRLD with different number of samples, the value is set to be boldface if it has better average performance than SVGD. If it is statistical significant with significant level 0.05 using a matched pair t-test, we add an underline on it.

Figure 9 and 10 give some visualized result on the comparison with Langevin dynamics and SRLD. To rule out the variance of different splitting on the dataset, the errorbar is calculated based on the difference between RMSE of SRLD and RMSE of Langevin dynamcis in 20 repeats (And similarily for test log-likelihood). And we only applied the error bar on Langevin dynamics. Contextual bandit is a class of online learning problems that can be viewed as a simple reinforcement learning problem without transition. For a completely understanding of contextual bandit problems, we refer the readers to the Chapter 4 of (Bubeck et al., 2012) . Here we include the main idea for completeness. In contextual bandit problems, the agent needs to find out the best action given some observed context (a.k.a the optimal policy in reinforcement learning). Formally, we define S as the context set and K as the number of action. Then we can concretely describe the contextual bandit problems as follows: for each time-step t = 1, 2, • • • , N , where N is some pre-defined time horizon (and can be given to the agent), the environment provides a context s t ∈ S to the agent, then the agent should choose one action a t ∈ {1, 2, • • • , 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: Num of samples where π denotes the deterministic mapping from the context set S to actions {1, 2, • • • , K} (readers can view π as a deterministic policy in reinforcement learning). Intuitively, this pseudo-regret measures the difference of cumulative reward between the action sequence a t and the best action sequence π(s t ). Thus, an algorithm that can minimize the pseudo-regret (7) can also find the best π.

EQUATION

Posterior sampling (a.k.a. Thompson sampling; Thompson, 1933) is one of the classical yet successful algorithms that can achieve the state-of-the-art performance in practice (Chapelle & Li, 2011) . It works by first placing an user-specified prior µ 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 µ 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 = arg max a∈{1,2,••• ,K} r(s t , a), r(s t , a) ∼ µ 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 µ 0 s,a , time horizon N for time t = 1, 2, • • • , N do observe a new context s t ∈ S,

sample the reward of each action r(s t , a) ∼ µ t s,a , a ∈ {1, 2, • • • , K}, select action a t = arg max a∈{1,2,••• ,K} r(s t , a) and get the reward r(s t , a t ), update the posterior µ t+1 st,at 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.

Here, we test the uncertainty provided by vanilla Langevin dynamics and Self-repulsive Langevin dynamics on two of the benchmark contextual bandit problems suggested by (Riquelme et al., 2018) , called mushroom and wheel. One can read (Riquelme et al., 2018) to find the detail introduction of this two contextual bandit problems. For completeness, we include it as follows:

Mushroom Mushroom bandit utilizes the data from Mushroom dataset (Schlimmer, 1981) , which includes different kinds of poisonous mushroom and safe mushroom with 22 attributes that can indicate whether the mushroom is poisonous or not. Blundell et al. (2015) first introduced the mushroom bandit by designing the following reward function: eating a safe mushroom will give a +5 reward, while eating a poisonous mushroom will return a reward +5 and -35 with equal chances. The agent can also choose not to eat the mushroom, which always yield a 0 reward. Same to (Riquelme et al., 2018) , we use 50000 instances in this problem.

Wheel To highlight the need for exploration, (Riquelme et al., 2018) designs the wheel bandit, that can control the need of exploration with some "exploration parameter" δ ∈ (0, 1). The context set S is the unit circle s 2 ≤ 1 in R 2 , and each turn the context s t is uniformly sampled from S. K = 5 possible actions are provided: the first action yields a constant reward r ∼ N (µ 1 , σ 2 ); the reward corresponding to other actions is determined by the provided context s:

• For s ∈ S s.t. s 2 ≤ δ, all of the four other actions return a suboptimal reward sampled from N (µ 2 , σ 2 ) for µ 2 < µ 1 . • For s ∈ S s.t. s 2 > δ, according to the quarter the context s is in, one of the four actions becomes optimal. This optimal action gives a reward of N (µ 3 , σ 2 ) for µ 3 µ 1 , and another three actions still yield the suboptimal reward N (µ 2 , σ 2 ).

Following the setting from (Riquelme et al., 2018) , we set µ 1 = 1.2, µ 2 = 1.0, and µ 3 = 50.

When δ approaches 1, the inner circle s 2 ≤ δ will dominate the unit circle and the first action becomes the optimal for most of the context. Thus, inference methods with poorly estimated uncertainty will continuously choose the suboptimal action 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 δ = 0.95, which is pretty hard for existing inference methods as shown in (Riquelme et al., 2018) , and use 50000 contexts for evaluation.

Figure 11 : Visualization of the wheel bandit (δ = 0.95), taken from (Riquelme et al., 2018) .

Experiment Setup Following (Riquelme et al., 2018) , we use a feed-forward network with two hidden layer of 100 units and ReLU activation. We use the same step-size and thinning factor c/η = 100 for vanilla Langevin dynamics and SRLD, and set M = 20, α = 1 on both of the mushroom and wheel bandits. The update schedule is similar to (Riquelme et al., 2018) , and we just change the optimization step in stochastic variational inference methods into MCMC sampler step and use a higher (2000) initial steps for burn-in like the warm-up of stochastic variational inference methods in Riquelme et al. (2018) . As this is an online posterior inference problem, we only use the last 20 samples to give the prediction. Notice that, in the original implementation of Riquelme et al. (2018) , the authors only update a few steps with new observation after observing enough data, as the posterior will gradually converge to the true reward distribution and little update is needed after observing sufficient data. Similar to their implementation, after observing enough data, we only collect one new sample with the new observation each time. For SVGD, we also use 20 particles to make the comparison fair.

B THE DETAILED ANALYSIS OF SRLD B.1 SOME ADDITIONAL NOTATION:

We use • ∞ to denote the ∞ vector norm and define the L ∞ norm of a function f : R d → R 1 as f L∞ . D TV denote the Total Variation distance between distribution ρ 1 , ρ 2 respectively. Also, as

K is R d × R d → R 1 , we denote K L∞,L∞ = sup θ1,θ2 K(θ 1 , θ 2 ).

For simplicity, we may use K ∞,∞ as K L∞,L∞ . In the appendix, we also use φ[ρ](θ) := g(θ; ρ), where g(θ; ρ) is defined in the main text. For the clearance, we define π M,c/η * ρ k := ρ M k , π M,c/η * ρk := ρM k and π M,c * ρt := ρM t , where ρ M k , ρM k and ρM t are defined in main text.

B.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 (M c/η + 1)-order nonlinear autoregressive model when k ≥ M c/η:

θ k+1 = ψ θ k , ..., θ k-M c/η + 2ηe k , where ψ θ k , ..., θ k-M c/η = θ k + η    -∇V (θ k ) + αφ[ 1 M M j=1 δ θ k-jc/η ](θ k )    .

Further, if we stack the parameter by

Θ k = θ k , ..., θ k-M c/η and define Ψ (Θ k ) = ψ (Θ k ) , Θ k , we have Θ k+1 = Ψ (Θ k ) + 2ηE k ,

where E k = e k , 0 , ..., 0 . In this way, we formulate Θ k as a time homogeneous Markov Chain. In the following analysis, we only analyze the second phase of SRLD given some initial stacked particles Θ M c/η-1 . Theorem B.1. (Geometric Ergodicity) Under Assumption 4.1 and Assumption 4.2, suppose we choose η and α such that

max 1 -2ηa 1 + η 2 b 1 + 2αη σ b 1 , 2αη σ (b 1 + 1) < 1,

then the Markov Chain of Θ k is stationary, geometrically ergodic, i.e., for any Θ 0 = Θ M c/η-1 , we have

D TV P k (•, Θ 0 ) , Π (•) ≤ Q (Θ 0 ) e -rk

, where r is some positive constant, Q(Θ 0 ) is constant related to Θ 0 , P k is the k-step Markov transition kernel and Π is the stationary distribution.

We defer the proof to Appendix B.5.1.

B.3 MOMENT BOUND Theorem B.2. (Moment Bound) Under assumption 4.2, suppose that we have E θ∼ρ0 θ 2 < ∞; and a 2 -α K ∞ 2b 1 + 4 σ > 0, we have sup k E θ∼ρ k θ 2 ∨ sup k E θ∼ ρk θ 2 ∨ sup t E θ∼ ρt θ 2 ≤E θ∼ρ0 θ 2 + b 1 + 1 + η a 2 -K L∞,L∞ 2α σ -α K L∞,L∞ 2b 1 + 2 σ .

And by Lemma B.1, we thus have

sup k E θ∼ρ k ∇V (θ) 2 ∨ sup k E θ∼ ρk ∇V (θ) 2 ∨ sup t E θ∼ ρt ∇V (θ) 2 ≤b 1 E θ∼ρ0 θ 2 + b 1 (b 1 + 1 + η) a 2 -K L∞,L∞ 2α σ -α K L∞,L∞ 2b 1 + 2 σ + 1

The proof can be found at Appendix B.5.2.

B.4 TECHNICAL LEMMA Definition B.1. (α-mixing):

For any two σ-algebras A and B, the α-mixing coefficient is defined by

α(A, B) = sup A∈A,B∈B |P (A ∩ B) -P (A) P (B)| .

Let (X k , k ≥ 1) be a sequence of real random variable defined on (Ω, A, P). This sequence is

α-mixing if α(n) := sup k≥1 α (M k , G k+n ) → 0, as n → ∞,

where M j := σ (X i , i ≤ j) and G j := σ (X i , i ≥ j) for j ≥ 1. Alternatively, as shown by Theorem 4.4 of Bradley ( 2007)

α(n) := 1 4 sup Cov (f, g) f L∞ g L∞ ; f ∈ L ∞ (M k ) , g ∈ L ∞ (G k+n ) . Definition B.2. (β-mixing)

For any two σ-algebras A and B, the α-mixing coefficient is defined by

β (A, B) := sup 1 2 I i=1 J j=1 |P (A i ∩ B j ) -P (A i ) P (B j )| ,

where the supremum is taken over all pairs of finite partitions {A 1 , ..., A I } and {B 1 , ..., B J } of Ω such that A i ∈ A and B j ∈ B for each i, j. Let (X k , k ≥ 1) be a sequence of real random variable defined on (Ω, A, P). This sequence is β-mixing if

β(n) := sup k≥1 β (M k , G k+n ) → 0, as n → ∞.

Proposition B.1. (β-mixing implies α-mixing)

For any two σ-algebras A and B,

α (A, B) ≤ 1 2 β (A, B) .

This proposition can be found in Equation 1.11 of Bradley (2005) . Proposition B.2. A (strictly) stationary Markov Chain is geometric ergodicity if and only if β(n) → 0 at least exponentially fast as n → ∞.

This proposition is Theorem 3.7 of Bradley (2005) . Lemma B.1. By Assumption 4.2, we have

∇V (θ) ≤ b 1 ( θ 1 + 1) and θ -η∇V (θ) ≤ 1 -2ηa 1 + η 2 b 1 θ 2 + η 2 b 1 + 2ηb 1 .

Lemma B.2. (Some property of RBF Kernel) For RBF kernel with bandwidth σ, we have K ∞,∞ ≤ 1 and

K(θ , θ 1 ) -K(θ , θ 2 ) ≤ e -(•) 2 /σ Lip θ 1 -θ 2 2 ∇ θ K(θ , θ 1 ) -∇ θ K(θ , θ 2 ) ≤ 2 σ e -(•) 2 /σ (•) Lip θ 1 -θ 2 2 .

Lemma B.3. (Some property of Stein Operator)

For any distribution ρ such that E θ∼ρ ∇V (θ) < ∞, we have θ) . We also denote φθ (θ) = [ φθ ,1 (θ), ..., φθ ,d (θ)] . We have The proof of this theorem is by verifying the condition of Theorem 3.2 of An & Huang (1996) . Suppose Θ = [θ 1 , ..., θ M C+1 ], where C = c/η, we have

φ[ρ](•) Lip ≤ e -(•) 2 /σ Lip E θ∼ρ ∇V (θ) + 2 σ e -(•) 2 /σ (•) Lip , φ[ρ](θ) ≤ K ∞ E θ ∼ρ ∇V (θ ) + 2 σ ( θ + θ ) ≤ K ∞ b 1 + E θ ∼ρ 2 σ + b 1 θ + θ . Lemma B.4. (Bounded Lipschitz of Stein Operator) Given θ , define φθ (θ) := φ[δ θ ](θ) = K(θ , θ)∇V (θ ) + ∇ 1 K(θ ,

EQUATION

ψ (Θ) = θ 1 + η    -∇V (θ 1 ) + αφ[ 1 M M j=1 δ θ 1+jC ](θ k )    = θ 1 -η∇V (θ 1 ) + ηα M M j=1 e -θ 1+jC -θ1 2 /σ 2 σ (θ 1 -θ 1+jC ) -e -θ 1+jC -θ1 2 /σ ∇V (θ 1+jC ) ≤ θ 1 -η∇V (θ 1 ) + 2 σ ηα M M j=1 e -θ 1+jC -θ1 2 /σ θ 1 + ηα M M j=1 e -θ 1+jC -θ1 2 /σ 2 σ (-∇V (θ 1+jC ) -θ 1+jC ) ≤ θ 1 -η∇V (θ 1 ) + 2αη σ K ∞,∞ b 1 (1 + θ 1 ) + 2αη M σ M j=1 K ∞,∞ b 1 1 + (1 + 1 b 1 ) θ 1+jC (1) ≤ b 1 (1 + 4αη σ K ∞,∞ ) + η 2 b 1 + 2ηb 1 + 1 -2ηa 1 + η 2 b 1 + 2αη σ K ∞,∞ b 1 θ 1 + 2αη σ K ∞,∞ (b 1 + 1) max i∈[M C+1]-{1} θ 1+jC ≤ b 1 (1 + 4αη σ K ∞,∞ ) + η 2 b 1 + 2ηb 1 + max 1 -2ηa 1 + η 2 b 1 + 2αη σ K ∞,∞ b 1 , 2αη σ K ∞,∞ (b 1 + 1) max i∈[M C+1] θ 1+jC ,

where (1) is by Lemma B.1. Thus, given the step size η, if we choose η, α such that

max 1 -2ηa 1 + η 2 b 1 + 2αη σ K ∞,∞ b 1 , 2αη σ K ∞,∞ (b 1 + 1) < 1,

then our dynamics is geometric ergodic.

B.5.2 PROOF OF THEOREM B.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 α = 0, here we only analysis the second phase. Define U t = sup s≤t E θs 2 , and thus

∂ ∂t U 2 t ≤ E θt , -V ( θ) + αφ[π M,c * ρt ]( θt ) ∨ 0. Now we bound E θt , -V ( θ) + αφ[π M,c * ρt ]( θt ) : E θt , -V ( θt ) + αφ[π M,c * ρt ]( θt ) ≤b 1 -a 2 E θt 2 + αE θt φ[π M,c * ρt ]( θt ) (1) ≤ b 1 -a 2 E θt 2 + α K ∞ E θt E θ ∼π M,c * ρt ∇V (θ ) + 2 σ ( θ + θ t ) ≤b 1 -a 2 E θt 2 + α K ∞ E θt E θ ∼π M,c * ρt b 1 ( θ + 1) + 2 σ ( θ + θ t ) =b 1 -a 2 -K ∞ 2α σ E θt 2 + α K ∞ E θt E θ ∼π M,c * ρt b 1 + 2 σ θ + b 1 ≤b 1 -a 2 -K ∞ 2α σ U 2 t + α K ∞ E θt E θ ∼π M,c * ρt b 1 + 2 σ θ + b 1 ≤b 1 -a 2 -K ∞ 2α σ U 2 t + α K ∞ b 1 + 2 σ 1 M M j=1 U t U t-jc + α K ∞ b 1 U t ≤b 1 -a 2 -K ∞ 2α σ U 2 t + α K ∞ b 1 + 2 σ U 2 t + α K ∞ b 1 (U 2 t + 1) ≤ (b 1 + 1) -a 2 -K ∞ 2α σ -α K ∞ 2b 1 + 2 σ U 2 t ,

where (1) is by B.3. By the assumption that λ :

= a 2 -K ∞ 2α σ -α K ∞ 2b 1 + 2 σ > 0, we have ∂ ∂t U 2 t ≤ (b 1 + 1) -λU 2 t ∨ 0. By Gronwall's inequality, we have U 2 t ≤ U 2 0 + b1+1 λ . (If ∂ ∂t U 2 t = 0

, then U t fix and this bound still holds.) Notice that in the first phase, as α = 0, we have λ < a 2 and thus this inequality also holds. 

Discrete-Time Mean Field Dynamics:

U 2 k -U 2 k-1 ≤ 2ηE θk-1 , -∇V ( θk ) + αφ[π M,c/η * ρk ]( θk ) + 2η 2 ∨ 0.

By a similarly analysis, we have bound

E θk-1 , -∇V ( θk ) + αφ[π M,c/η * ρk ]( θk ) ≤ (b 1 + 1) -λU 2 t , where λ = a 2 -K ∞,∞ 2α σ -α K ∞,∞ 2b 1 + 2 σ > 0. And thus we have U 2 k -U 2 k-1 ≤ 2η (b 1 + 1) -λU 2 k-1 + 2η 2 ∨ 0. It gives that U 2 k ≤ b 1 + 1 + η λ + U 2 0 .

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.

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

∂ t ρ t = ∇ • ([-(D(θ) + Q(θ)) ∇V (θ) + αφ[π M,c * ρ t ](θ t )] ρ t + (D(θ) + Q(θ)) ∇ρ t ) .

Notice that the limiting distribution satisfies

0 a.e. = ∇ • ([-(D(θ) + Q(θ)) ∇V (θ) + αφ[π M,c * ρ ∞ ](θ t )] ρ ∞ + (D(θ) + Q(θ)) ∇ρ ∞ ) = ∇ • ([-(D(θ) + Q(θ)) ∇V (θ) + αφ[ρ ∞ ](θ t )] ρ ∞ + (D(θ) + Q(θ)) ∇ρ ∞ ) = ∇ • ([-(D(θ) + Q(θ)) ∇V (θ)] ρ ∞ + (D(θ) + Q(θ)) ∇ρ ∞ ) + α∇ • (K * (∇ρ ∞ -∇V (θ)ρ ∞ ) ρ ∞ ) .

which implies that ρ ∞ ∝ exp(-V (θ)) is the stationary distribution.

B.5.4 PROOF OF THEOREM 4.2

In the later proof we use c d to represent the quantity

E θ∼ρ0 θ 2 + b 1 + 1 + η a 2 -K ∞,∞ 2α σ -α K ∞,∞ 2b 1 + 2 σ .

Recall that there are two dynamics: the continuous-time mean field dynamics ( 5) and the discretized version discrete-time mean field Dynamics (4). Notice that here we couple the discrete-time mean field dynamics with the continuous-time mean field system using the same initialization. Given any T = ηN , for any 0 ≤ t ≤ T , define t = t η η. We introduce an another continuous-time interpolation dynamics:

θt = -∇V ( θt ) + dB t , t ∈ [0, M c) -∇V ( θt ) + αφ[π M,c * ρt ]( θt ) + dB t , t ≥ M c, ρt = Law( θt ), θ0 = θ0 ∼ ρ0 ,

Notice that here we couples this interpolation dynamics with the same Brownian motion as that of the dynamics of θt . By the definition of θt , at any t k := kη for some integrate k ∈ [N ], θt k and θk has the same distribution. Define ρθ0 t = Law( θt ) conditioning on θ0 = θ 0 and ρθ0 t = Law( θt ) conditioning on θ0 = θ 0 . Followed by the argument of proving Lemma 2 in Dalalyan (2017), if k ≥ M c η , we have

D KL ρθ0 t k ρθ0 t k = 1 4 t k 0 E -∇V ( θs ) + αφ[π M,c * ρs ]( θs ) + ∇V ( θs ) -αφ[π M,c * ρs ]( θs ) 2 ds = 1 4 k-1 j=0 tj+1 tj E -∇V ( θtj ) + αφ[π M,c * ρtj ]( θtj ) + ∇V ( θs ) -αφ[π M,c * ρs ]( θs ) 2 ds ≤ 3 4 k-1 j=0 tj+1 tj E ∇V ( θtj ) -∇V ( θs ) 2 ds + 3α 2 4 k-1 j=0 tj+1 tj E φ[π M,c * ρtj ]( θtj ) -φ[π M,c * ρs ]( θtj ) 2 ds ≤ 3α 2 4 k-1 j=0 tj+1 tj E φ[π M,c * ρs ]( θtj ) -φ[π M,c * ρs ]( θs ) 2 ds =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 ∇V , we have

∇V ( θtj ) -∇V ( θs ) 2 ≤ b 2 1 θtj -θs 2 .

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

[π M,c * ρs ](θ 1 ) -φ[π M,c * ρs ](θ 2 ) ≤ e -(•) 2 /σ Lip E θ∼π M,c * ρs ∇V (θ) + 2 σ e -(•) 2 /σ (•) Lip θ 1 -θ 2 2 .

And by the Assumption 4.2 and that ρs as finite second moment, we have

φ[π M,c * ρs ](θ 1 ) -φ[π M,c * ρs ](θ 2 ) ≤Cc d θ 1 -θ 2 2 .

Combine the two bounds, we have

I 1 + I 3 ≤ 3Cc 2 d 4 k-1 j=0 tj+1 tj E θtj -θs 2 ds.

Notice that θt = θt + -∇V ( θt ) + αφ[π M,c * ρt ]( θt ) (t -t) + t t dB s . By Itô's lemma, it implies that

I 1 + I 3 ≤ 3Cc 2 d 4 k-1 j=0 tj+1 tj E θtj -θs 2 ds ≤ 3Cc 2 d 4 tj+1 tj E -∇V ( θs ) + αφ[π M,c * ρs ]( θs ) 2 (s -t j ) 2 + 2d(s -t j ) ds = Cc 2 d η 3 k-1 j=0 E -∇V ( θtj ) + αφ[π M,c * ρtj ]( θtj ) 2 + Cc 2 d dkh 2 .

By the assumption that E θtj is finite and θtj

d = θtj , E θtj 2 is also finite, we have E -∇V ( θt ) + αφ[π M,c * ρt ]( θt ) 2 ≤2E ∇V ( θt ) 2 + 2α 2 E φ[π M,c * ρt ]( θt ) 2 ≤4b 2 1 + 4b 2 1 E θt 2 + 2α 2 E 2 σ + b 1 E θ ∼π M,c * ρt θ + θ 2 ≤c 2 d C. Thus we conclude that I 1 + I 3 ≤ Cc 2 d c 2 d kη 3 + dkη 2 . Bounding I 2 E φ[π M,c * ρtj ]( θtj ) -φ[π M,c * ρs ]( θtj ) 2 =E 1 M M l=1 φ[ρ tj -cl ]( θtj ) -φ[ρ s-cl ]( θtj ) 2 ≤ 1 M M l=1 E φ[ρ tj -cl ]( θtj ) -φ[ρ s-cl ]( θtj ) 2 = 1 M M l=1 E E θ∼ ρt j -cl φ θt j (θ) -E θ∼ ρs-cl φ θt j (θ) 2 = 1 M M l=1 E θt j d i=1 E θ∼ ρt j -cl φ θt j ,i (θ) -E θ∼ ρs-cl φ θt j ,i (θ) 2 ≤ 1 M M l=1 E θt j d i=1 φ θt j ,i (•) L∞ ∨ φ θt j ,i (•) Lip 2 D 2 BL ρtj-cl , ρs-cl

By Lemma B.4 and the Assumption 4.4 that V is at most quadratic growth and that ρt has finite second moment, we have

E θt j d i=1 φ θt j ,i (•) L∞ ∨ φ θt j ,i (•) Lip 2 =E θt j d i=1 φ θt j ,i (•) 2 L∞ ∨ φ θt j ,i (•) 2 Lip ≤ 4d 2 σ e -θ 2 /σ θ 1 2 BL + 4 e -• 2 /σ 2 BL E θt j ∇V ( θtj ) 2 ≤C(d + c 2 d ).

Plug in the above estimation, we have

I 2 = 3α 2 4 k-1 j=0 tj+1 tj E φ[π M,c * ρtj ]( θtj ) -φ[π M,c * ρs ]( θtj ) 2 ds ≤ α 2 C(d + c 2 d ) k-1 j=0 tj+1 tj 1 M M l=1 D 2 BL ρtj-cl , ρs-cl ds ≤ α 2 C(d + c 2 d ) k-1 j=0 1 M M l=1 tj+1 tj D KL ρtj-cl , ρs-cl ds,

where the last inequality is due to the relation that D

2 BL definition ≤ D 2 TV Pinsker s ≤ D KL .

Overall Bound Combine all the estimation, we have

D KL ρθ0 t k ρθ0 t k ≤ α 2 C(d + c 2 d ) k-1 j=0 1 M M l=1 tj+1 tj D KL ρtj-cl , ρs-cl ds + Cc 2 d c 2 d kη 3 + dkη 2 = α 2 C(d + c 2 d ) k-1 j=0 1 M M l=1 η 0 D KL ρt ( jη-cl η ) , ρt ( jη-cl η ) +s ds + Cc 2 d c 2 d kη 3 + dkη 2 Similar, if k ≤ M c η -1, we have D KL ρθ0 t k ρθ0 t k = 1 4 t k 0 E ∇V ( θs ) -∇V ( θs ) 2 ds ≤ b 2 1 4 k-1 j=0 tj+1 tj E θtj -θs 2 ds ≤ b 2 1 η 3 12 k-1 j=0 E ∇V ( θtj ) 2 + dkb 2 1 η 2 4 ≤ b 2 1 η 3 kc 2 d 12 + dkb 2 1 η 2 4 . Define u k = sup s∈[t k ,t k+1 ] D KL ρθ0 s ρθ0 s ,

and U k = max l∈{0,...,k} u l . We conclude that for k ≥ M c η , for any k ≤ k,

u k ≤ α 2 C(d + c 2 d ) k-1 j=0 1 M M l=1 h 0 D KL ρt ( jη-cl η ) , ρt ( jη-cl η ) +s ds + Cc 2 d c 2 d kη 3 + dkη 2 ≤ α 2 C(d + c 2 d ) k-1 j=0 1 M M l=1 ηu ( jη-cl η ) + Cc 2 d c 2 d kη 3 + dkη 2 ≤ α 2 C(d + c 2 d )η k-1 j=0 U j + Cc 2 d c 2 d kη 3 + dkη 2 .

For k < M c η , which is a simpler case, we have

U k ≤ C η 3 kc 2 d + dkη 2 < CM c ηc 2 d + d η. We bound the case when k ≥ M c η , U k ≤ α 2 C(d + c 2 d )η k-1 j=0 U j + Cc 2 d c 2 d kη 3 + dkη 2 .

If we take η sufficiently small, such that c 2 d kη 3 ≤ dkη 2 , we have

U k ≤ α 2 C(d + c 2 d )η k-1 j=0 U j + 2Cc 2 d dkη 2 ≤ α 2 C(d + c 2 d )η k-1 j=0 (U j + η) . Define η = α 2 C(d + c 2 d

)η and we can choose η small enough such that η < 1/2 and η < 1/2. Without loss of generality, we also assume η ≥ η and thus we have

U k ≤ η k-1 j=0 (U j + η ) .

Also we assume U k ≥ η , otherwise we conclude that U k < η . We thus have U k ≤ q k-1 j=0 U j , where q = 2η . Suppose that U M c η -1 = x ≤ CM c ηc 2 d + d η and some algebra (which reduces to Pascal's triangle) shows that for some positive constant r that characterize the decay rate of α mixing. Combine this two estimations, we have Bounding I 2 By Holder's inequality, we have

U k ≤ xq(1 + q) k-M c η . We conclude that U k ≤ xq(1 + q) k-1 . Notice that q = 2α 2 C(d + c 2 d )η. Thus for any k ≥ M c/η, U k ≤ xq(1 + q) k-M c η = xq(1 + q) (kη-M c)/η = xq(1 + q) 2α 2 C(d+c 2 d )(kη-M c)/q ≤ x2α 2 C(d + c 2 d )e 2α 2 C(d+c 2 d )(kη-M c) η ≤ CM cα 2 ηc 2 d + d (d + c 2 d )e 2α 2 C(d+c 2 d )(kη-M c) η 2 ,

I 1 ≤ E θ k -θk

I 2 ≤ E θ k -θk 2 E φ[π M,c/η * ρ k ](θ k ) -φ[π M,c/η * ρk ](θ k ) 2 .

We bound the second term in the right side of the inequality. And combined all the estimation and by the definition of Wasserstein-distance, we conclude that

I 2 ≤ 4d 2 σ e -θ 2 /σ θ 1 2 BL + 4 e -• 2 /σ 2 BL E θ k ∇V (θ k ) 2 1 M M j=1 W 2 2 [ρ k-cj , ρk-cj ] ≤ 4d 2 σ e -θ 2 /σ θ 1 2 BL + 4 e -• 2 /σ 2 BL E θ k ∇V (θ k ) 2 1 M M j=1 E θ k-cj -θk-cj 2 .

Bounding I 3 By Holder's inequality, Overall Bound Combine all the results, we have the following bound: for k ≥ M c, u l , we have

I 3 ≤ E θ k -θk

E θ k+1 -θk+1 2 -E θ k -θk 2 ≤ -2ηLE θ k -θk 2 +2ηα E θ k -θk 2 c 1 √ M +2ηαc 2 1 M M j=1 E θ k-jc/η -θk-jc/η 2 E θ k -θk 2 +2ηαc 3 E θ k -θk

U 2 k+1 ≤ qU 2 k + 2ηαc 1 √ M U k + η 2 c 4 ,

where q = (1 -2η(L -αc 2 -αc 3 )). By the assumption that α ≤ L/(c 2 + c 3 ), q < 1. Now we prove the bound of U k by induction. We take the hypothesis that

U 2 k ≤ 2ηαc 1 √ M +(1-q)η(c4+ 1 1-q ) 2 (1-q) 2

For the second part:

θ -η∇V (θ) = θ -η∇V (θ), θ -η∇V (θ) = θ 2 + 2η θ, -∇V (θ) + η 2 ∇V (θ) For any distribution ρ such that ∇ θ V (θ) ρ(θ)dθ < ∞,

2 ≤ θ 2 + 2η -a 1 θ 2 + b 1 + η 2 b 1 (1 + θ 2 ) = 1 -2ηa 1 + η 2 b 1 θ 2 + η 2 b 1 + 2ηb 1 .

φ[ρ](θ 1 ) -φ[ρ](θ 2 ) = E θ∼ρ {-[K(θ, θ 1 ) -K(θ, θ 2 )] ∇V (θ) + ∇ 1 K(θ, θ 1 ) -∇ 1 K(θ, θ 2 )} ≤ e -(•) 2 /σ Lip E θ∼ρ ∇V (θ) θ 1 -θ 2 2 + 2 σ e -(•) 2 /σ (•) Lip θ 1 -θ 2 2 .

For proving the second result, we notice that 

φ[ρ](θ) = E θ ∼ρ [K(θ , θ)∇V (θ ) + ∇ 1 K(θ , θ)] ≤ K ∞ E θ ∼ρ ∇V (θ ) + 2 σ ( θ + θ ) ≤ K ∞ b 1 + E θ ∼ρ 2 σ + b 1 θ + θ .

:



Figures and Tables:

Figure 1: An illustrative example showing the advantage of our Self-Repulsive Langevin dynamics.With a set of initial examples locating on the left part of the target distribution (show in yellow), 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.

Figure 3: 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 4: Sampling trajectory of the correlated 2D distribution.

Figure 6: Sampling trajectory of the mixture of Gaussian.

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

Figure 9: 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 10: 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 SRLDlog-likelihood of Langevin dynamics in 20 repeats to rule out the variance of data splitting.

(4) Similarly to the analysis of the continuous-time mean field dynamics (5), we only give proof of the second phase. Define U k = sup s≤k

for sufficiently small η. Combine the above two estimations, we haveU k ≤ C η 3 kc 2 d + dkη 2 + η k ≤ M c/η -1 CM cα 2 ηc 2 d + d (d + c 2 d )e 2α 2 C(d+c 2 d )(kη-M c) η 2 + Cη k ≥ M c/η .Notice that now we have U k = max l∈{0,...,k} sup s∈[0,η] D KL ρθ0 lη+s ρθ0 lη , which is a function of θ 0 . We then bound Ūk = max l∈{0,...,k} sup s∈[0,η] D KL [ρ lη+s ρlη ]. Notice that the KL divergence has the following variational representation:D KL [ρ 1 ρ 2 ] = sup f E ρ1 f -E ρ2 e f ,where the f is chosen in the set that E ρ1 f and E ρ2 e f exist. And thus we haveD KL [ρ lη+s ρlη ] = sup f KL [ρ lη+s ρlη ] ≥ max l∈{0,...,k} D KL [ρ lη ρlη ]holds naturally by definition. We complete the proof. •) 2 /σ (•) Lip Now we start the proof. We couple the process of θ k and θk by the same gaussian noise e k in every iteration and same initialization θ0 = θ 0 . For k ≤ M c/η -1, E θ k -θk ≤ Ca 2 M e -rc 1 -e -rM c 1 -e rc + 1 ,

(1 -e -rM c ) 1 -e rc + 1 .

φ[π M,c/η * ρ k ](θ k ) -φ[π M,c/η * ρk ](θ k ) k-jc/η ](θ k ) -φ[ρ k-jc/η ](θ k ) ρ k-jc/η ](θ k ) -φ[ρ k-jc/η ](θ k ) θ k E θ∼ρ k-jc/η φθ k (θ) -E θ∼ ρk-jc/η φθ k (θ) θ∼ρ k-jc/η φθ k ,i (θ) -E θ∼ ρk-jc/η φθ k ,i (θ) k ,i (•) L∞ ∨ φθ k ,i (•) Lip 2 D 2 BL ρ k-jc/η , ρk-jc/η .By Lemma B.4, we haved i=1 φθ k ,i (•) L∞ ∨ φθ k ,i (•)Plug in the above estimation and by the relation thatD BL ≤ W 1 ≤ W 2 , we have E φ[π M,c * ρ k ](θ k ) -φ[π M,c * ρk ](θ k ) k-cj , ρk-cj ] .

φ[π M,c/η * ρk ](θ k ) -φ[π M,c/η * ρk ]( θk ) 2We bound the last term on the right side of the inequality. By assumption and Lemma B.3, we haveE φ[π M,c/η * ρk ](θ k ) -φ[π M,c/η * ρk ]( θk ) 2 ≤ e -(•) 2 /σ Lip E θ∼ ρk ∇V (θ) + 2 σ e -(•) 2 /σ (•) Lip 2And combine the estimation, we haveI 3 ≤ e -(•) 2 /σ Lip E θ∼ ρk ∇V (θ

-jc/η ](θ k ) -∇V ( θk ) -αφ[π M,c/η * ρk ]( θk ) 2 Define u k = E θ k -θk 2 and U k = sup l∈[k]

.2 PROOF OF LEMMA B.2 It is obvious that K ∞,∞ ≤ 1. K(θ , θ 1 ) -K(θ , θ 2 ) e -θ -θ1 2 /σ -e -θ -θ2 2 /σ ≤ e -(•) 2 /σ Lip θ 1 -θ 2 2 . And ∇ θ K(θ , θ 1 ) -∇ θ K(θ , θ 2 ) = 2 σ e -θ -θ1 2 /σ (θ -θ 1 ) -2 σ e -θ -θ2 2 /σ (θ -θ 2 )

.6.4 PROOF OF LEMMA B.4 Given any θ , Lip .To bound the second term, by the symmetry of each coordinates, we have θ -θ1 2 /σ (θ1,i -θ i ) -2 σ e -θ -θ2 2 /σ (θ 1,i -θ i ) -• 2 /σ 2 L∞ ∇V (θ ) 2 .

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

Cumulative Regrets on two bandits problem. 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.

Table content: <table><tr><td/><td>SVGD</td><td>LD</td><td>SRLD</td></tr><tr><td>Mushroom</td><td>20.7 ± 2.0</td><td>4.28 ± 0.09</td><td>3.80 ± 0.16</td></tr><tr><td>Wheel</td><td colspan="3">91.32 ± 0.17 38.07 ± 1.11 32.08 ± 0.75</td></tr></table>

Comparing SRLD with different number of samples with SVGD on test RMSE. The results

Table content: <table><tr><td>Dataset</td><td>SRLD-20</td><td>SRLD-40</td><td>SRLD-60</td><td colspan="2">Ave Test LL</td><td>SRLD-80</td><td>SRLD-100</td><td>SVGD</td></tr><tr><td>Boston</td><td colspan="8">-2.642 ± .088 -2.582 ± 0.084 -2.527 ± 0.612 -2.516 ± 0.062 -2.500 ± 0.054 -4.276 ± 0.217</td></tr><tr><td colspan="9">Concrete -3.084 ± 0.036 -3.061 ± 0.034 -3.050 ± 0.033 -3.040 ± 0.031 -3.034 ± 0.031 -5.500 ± 0.398</td></tr><tr><td>Energy</td><td colspan="8">-0.580 ± 0.053 -0.536 ± 0.048 -0.522 ± 0.046 -0.504 ± 0.044 -0.476 ± 0.036 -0.781 ± 0.094</td></tr><tr><td>Naval</td><td>4.033 ± 0.230</td><td>4.100 ± 0.171</td><td colspan="2">4.140 ± 0.015</td><td colspan="2">4.167 ± 0.014</td><td>4.186 ± 0.015</td><td>3.056 ± 0.034</td></tr><tr><td colspan="9">WineRed -1.008 ± 0.019 -0.990 ± 0.017 -0.982 ± 0.016 -0.974 ± 0.016 -0.970 ± 0.016 -1.040 ± 0.018</td></tr><tr><td colspan="9">WineWhite -1.053 ± 0.004 -1.049 ± 0.004 -1.047 ± 0.004 -1.044 ± 0.004 -1.043 ± 0.004 -1.040 ± 0.019</td></tr><tr><td>Yacht</td><td colspan="8">-1.160 ± 0.256 -0.650 ± 0.173 -0.556 ± 0.096 -0.465 ± 0.037 -0.458 ± 0.036 -1.281 ± 0.279</td></tr></table>

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.

Table content: <table/>

Peer review (iclr2020):
I have the rebuttal of the authors, the paper improved indeed and some point on role of M is better clarified now although it is still a bit convoluted. The paper would be stronger if the analysis shows any theoretical advantage to the presented method. I think the author put a good effort in addressing some of my concerns and I m raising my score to 6. #### Summary of the paper: The paper proposes stein self repulsive dynamics for sampling from an unnormalized distribution. The method starts by using Langevin dynamic for up to time and then uses those pasts samples to guide the trajectories of the langevin sampling to explore new areas of the densities using the stein witness function between the current particles and the past samples( similar to Stein Variational gradient descent). The paper analyses the mixing properties under standard assumptions of the potential of the Boltzmann distribution and the kernel used in the Stein discrepancy, and shows convergence to the boltzman distribution as the number of particles goes to infinity and the step size goes to zero. Authors validate their methods and show that it indeed explores on a synthetic example new areas wrt to pure langevin dynamics. Applications in sampling from the posterior of bayesian neural networks and contextual bandits compare the performance of the proposed method to langevin dynamic and pure stein descent favorably. Clarity/presentation : The paper is well written and the intuition are well presented . The notation is not great in denoting direct measures, please using another symbol. My main concerns with the paper are the following: - the definition of averages only particles choosen for time stamps. Something is off here for the continuous approximation to work at each time stamp you need particles and then you have a past horizon . As this does not matter, but I think in your implementation you are considering at each time step particles and you average on a horizon of size . is this correct? Please clarify? - The theorem show only asymptotic behavior and don*t quantify the intuition behind the paper , that the *coverage of the samples* is higher. Can you for instance bound the wasserstein distance between the pure langevin and Stein Repulsive dynamic , and between SVGD and your method, as it was done in *The promises and pitfalls of Stochastic Gradient Langevin Dynamics*. Basically you can find a coupling between trajectories of your methods and the langevin dynamics and bound the wasserstein distance between the two methods. This will be insightful to see if one would mix faster with respect to the other one. While the appendices of the paper are lengthy I don*t think they explain the most selling point of the method, since they are asymptotic and there is a confusion between the time horizon for the past , and the number of particles. if is in this defintion , we are just running langevin dynamic, and not using Stein witness function. I see an important issue in this definition of the time horizon, I think this time horizon should be finite, and that from time steps less we sample N particles , and we let this go to infinity.