Markov chain Monte Carlo is one of the most useful tools for obtaining samples from complex distributions, such as posterior distributions in Bayesian inference. Bouncy particle sampler (BPS) is a method that performs sampling using uniform linear motion and elastic collisions. Since BPS does not require time discretization, it is expected to be more advantageous than the widely used Hamiltonian Monte Carlo (HMC) and Langevin Monte Carlo (LMC) for distributions with large curvature changes. To sample more efficiently from multimodal distributions, we propose BPS with parallel tempering (BPS-PT) in the infinite exchange limit, based on LMC-PT with the infinite exchange rate proposed in a previous study. Using a distribution with large curvature changes, BPS-PT showed estimation accuracy comparable to HMC-PT and LMC-PT, and achieved higher ESS per computational time. Therefore, BPS-PT is considered to be effective for distributions with large curvature changes.
Figures & tables
Figure 1: We show a trajectory of a particle generated by the BPS algorithm in a potential of the normal distribution. The particle bounces at t=t1,t2,t4 and randomly changes its velocity at t=t3 .
HMC-PT
LMC-PT
BPS
BPS-PT
Optimization [sec]
6.3×102
1.18×103
-
1.09×103
Sampling [sec]
1.44×102
1.93×102
5×10
1×102
Table 1: Computation times of MGC are shown. The number of samples were 104 for BPS-PT and 105 for the other methods, respectively. To match the optimization time, the episode length for the policy gradient of BPS-PT is 1/10 of that of HMC-PT and LMC-PT. The optimization time and sampling time are the average of 5 and 10 runs, respectively.
HMC-PT
LMC-PT
BPS
BPS-PT
ESS
2×10−4
0.2
1.0
1.0
W1
0.5
0.35
0.8
0.34
KLD
4×10−2
1.3×10−1
∞
2×10−4
estimate of a
2.2×10−2
5.5×10−1
0.7
3.0×10−1
Table 2: Performance metrics of MGC are shown. To match the computation times, BPS-PT used 104 samples and the other method used 105 samples, respectively. Since no transitions between two circumferences occurred by BPS, the KLD is infinity although sample means of a was 0.7 .
Figure 2: 2D histograms of samples obtained by each method. To match the computation times, the number of samples was set to 104 for BPS-PT, and 105 for the others. Compared to the histogram obtained from the target pdf, each method had a following drawback. (a) exhibited anisotropy. (b) concentrated in the inner circumference than the outer one, which was opposite from the target pdf. (c) had no samples from the inner circumference. (d) had 1/10 samples per computation time.
HMC-PT
LMC-PT
BPS
BPS-PT
Optimization [sec]
8×10
4×102
-
2.2×102
Sampling [sec]
6×10
7.6×10
1.00×102
6.4×102
Table 3: Computation times of MNF are shown. The number of samples were 106 for all methods. To match the optimization time, the episode length for the policy gradient of BPS-PT is 1/10 of that of HMC-PT and LMC-PT. The optimization time and sampling time are the average of 5 and 10 runs, respectively.
HMC-PT
LMC-PT
BPS
BPS-PT
ESS
1×10−4
2×10−4
3×10−4
4×10−2
W1
17
5
4
5
KLD
2.7×10−3
2.2×10−5
2.7×10−6
1.5×10−4
estimate of a
3.3×10−1
3.0×10−1
3.0×10−1
2.9×10−1
Table 4: Performance metrics of MNF are shown. BPS-PT used 105 samples and the other methods used 106 samples to match computation times.
Figure 3: 2D histograms of samples obtained by each method. To match the computation times, the number of samples was set to 105 for BPS-PT, and 106 for the others, respectively. Compared to the histogram obtained from the target pdf, the horizontal spread ( x1 ) is narrower for all methods. The peaked part in the vertical direction ( x2 ) is reproduced more closely to the target pdf by BPS and BPS-PT than by HMC and LMC.
HMC-PT
LMC-PT
BPS
BPS-PT
tuning [sec]
1.09×102
1.43×102
-
2.43×102
sampling [sec]
1.36×102
1.21×102
8.1×10
1.21×103
Table 5: Computation time of GMM are shown. The number of samples were 106 for all methods. To match the optimization time, the episode length for the policy gradient of BPS-PT is 1/10 of that of HMC-PT and LMC-PT. The optimization time and sampling time are the average of 5 and 10 runs, respectively.
HMC-PT
LMC-PT
BPS
BPS-PT
ESS
8×10−4
6.4×10−2
5.8×10−1
3×10−4
W1
3.35×10
3.59×10
6.4×10
3.30×10
KLD
4.9×10−4
2.7×10−2
3.9
2.3×10−2
Table 6: Performance metrics of GMM are shown. To match calculation times, BPS-PT used 105 samples and the other method used 106 samples, respectively.
We develop a new low-accuracy sampler, called smoothed Picard Hamiltonian Monte Carlo, which combines Gaussian smoothing, Picard iteration, and higher-order discretization. For a log-concave target π∝exp(−V) in dimension d satisfying 0≺αI⪯∇2V⪯βI, with condition number κ:=β/α, smoothed Picard HMC returns a sample with αW2(⋅,π)≤ε using O(κ2+κ7/6d1/6/ε1/3) gradient queries. We also prove stronger Wq bounds, and then develop an algorithmic framework, the recursive warm start generator, to upgrade these Wq bounds to stronger divergence guarantees. This produces a warm start for the proximal bouncy particle sampler, introduced in a companion work, leading to a high-accuracy log-concave sampler with complexity O((κ7/6d1/6+κ1/2d1/4)polylog(1/ε)).
Fan Chen, Sinho Chewi, Jianfeng Lu +1
Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology. · Department of Statistics and Data Science, Yale University. · Department of Mathematics, Duke University. +1
We study the problem of sampling from μ(dx)∝e−V(x)dx on Rd, where V is α-strongly convex and β-smooth, and write κ:=β/α. We design and analyze the Proximal Bouncy Particle Sampler (Proximal BPS), a new sampler that combines ideas from the proximal sampler and the bouncy particle sampler. From a warm start initialization with O(1) Rényi divergence w.r.t. μ, Proximal BPS returns a sample whose law is ε-close to μ in total variation distance using O(κd1/4polylog(1/ε)) gradient queries in expectation.
Fan Chen, Sinho Chewi, Jianfeng Lu +1
Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology. · Department of Statistics and Data Science, Yale University. · Department of Mathematics, Duke University. +1
Sampling from multimodal distributions is a central challenge in Bayesian inference and machine learning. In light of hardness results for sampling -- classical MCMC methods, even with tempering, can suffer from exponential mixing times -- a natural question is how to leverage additional information, such as a warm start point for each mode, to enable faster mixing across modes. To address this, we introduce Reweighted ALPS (Re-ALPS), a modified version of the Annealed Leap-Point Sampler (ALPS) that dispenses with the Gaussian approximation assumption. We prove the first polynomial-time bound that works in a general setting, under a natural assumption that each component contains significant mass relative to the others when tilted towards the corresponding warm start point. Similarly to ALPS, we define distributions tilted towards a mixture centered at the warm start points, and at the coldest level, use teleportation between warm start points to enable efficient mixing across modes. In contrast to ALPS, our method does not require Hessian information at the modes, but instead estimates component partition functions via Monte Carlo. This additional estimation step is crucial in allowing the algorithm to handle target distributions with more complex geometries besides approximate Gaussian. For the proof, we show convergence results for Markov processes when only part of the stationary distribution is well-mixing and estimation for partition functions for individual components of a mixture.