Source-linked AI summary
Accelerated High-Accuracy Sampling from a Warm Start via the Proximal Bouncy Particle Sampler
Fan Chen, Sinho Chewi, Jianfeng Lu, Matthew S Zhang
TL;DR
The paper studies sampling from a probability measure proportional to e^-V(x) on R^d. It introduces Proximal BPS and reports eO(κ^1/2 d^1/4) expected gradient complexity from a warm start with D2 divergence O(1).
Problem
The paper studies sampling from a probability measure µ(dx)=Z^-1e^-V(x)dx on R^d under condition (1.1).
Source: §1 · p. 2Method
Proximal BPS uses a conditional half-turn in harmonic motion, stopping at time π to preserve position memory while avoiding a trivial return.
Source: §3.1 · p. 10Results
eO(κ^1/2 d^1/4) expected gradient complexity when Δ=O(1), with logarithmic factors in d, κ, and 1/ε suppressed.
Source: §4.3 · p. 20Takeaways & Limitations
Proximal BPS matches BPS in condition-number dependence while improving dimension dependence, with full acceleration and d^1/4 dimension dependence.
Source: §1 · pp. 2–3Obtaining the warm start is outside this paper; a companion work is slated to address it with eO(κ^7/6 d^1/6) gradient queries.
Source: §1 · p. 3
Abstract
from arXiv · showhide
We study the problem of sampling from $\mu(\mathrm{d}x)\propto e^{-V(x)}\,\mathrm{d}x$ on $\mathbb{R}^d$, where $V$ is $\alpha$-strongly convex and $\beta$-smooth, and write $\kappa:=\beta/\alpha$. 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\'enyi divergence w.r.t. $\mu$, Proximal BPS returns a sample whose law is $\varepsilon$-close to $\mu$ in total variation distance using $\widetilde O(\sqrt\kappa\,d^{1/4} \,\mathrm{polylog}(1/\varepsilon))$ gradient queries in expectation.
1 Introduction
The paper studies sampling from strongly convex, smooth targets and introduces Proximal BPS, which combines proximal sampling with particle motion. From a warm start, its main theorem gives a gradient-query complexity of eO(κ^1/2 d^1/4).
1 Introduction
- Prior high-accuracy samplers cited in the paper have complexities eO(κd^1/2), eO(κ^2/3d^1/3), and eO(κ^1/2d) for MALA/proximal sampling, FORS-based sampling, and BPS, respectively.The paper also reports eO(κd^5/4) partial-derivative queries for zigzag.Source: §1 · p. 2
- From a warm start, Proximal BPS achieves eO(κ^1/2 d^1/4) complexity, matching BPS's condition-number dependence while improving its dimension dependence.The paper says this also improves on prior results in both parameters.Source: §1 · pp. 2–3
- This paper does not address how to obtain the required warm start; a companion work is slated to do so using eO(κ^7/6d^1/6) gradient queries.The stated main result assumes an order-2 Rényi divergence warm start with D2(µ0∥µ)=O(1).Source: §1 · pp. 2–3
- Proximal BPS combines local particle motions that retain momentum across problems with bounces correcting for target shape.Nearby reference-point gradients make corrections inexpensive because the gradient varies little locally, requiring only a small number of bounces.Source: §1 · p. 3
- The sampler reflects the auxiliary point through the current position, then runs local piecewise-deterministic particle motion for half an oscillation period.The authors relate this reflection to Gaussian overrelaxation and say it helps information persist across iterations.Source: §1 · p. 3
- Existing analyses of continuously refreshed dynamics do not directly guarantee convergence for this discrete chain, so the authors adapt modified L2 hypocoercivity.Their local motion stops after a half-period without velocity refreshment and is combined with auxiliary reflection and occasional resampling.Source: §1 · p. 3
2 Preliminaries
The preliminaries introduce the standard Bouncy Particle Sampler, the proximal augmentation used by Proximal BPS, and the divergence and notation conventions used in the paper.
2.1 The standard Bouncy Particle Sampler
- BPS moves a position-velocity pair in straight lines and reflects velocity against the target gradient at bounces; independent velocity refreshes help ensure ergodicity.The process leaves µ ⊗ φ invariant, with position marginal µ; Poisson thinning simulates its event clock exactly when a dominating rate is available.Source: §2.1 · p. 4
2.2 Proximal augmentation and the restricted Gaussian oracle
- Proximal augmentation couples X∼µ to Y=X+√ηZ, yielding a Gaussian-smoothed target marginal and conditional laws for Gibbs updates.The restricted Gaussian oracle returns a draw from the conditional law of X given Y; the exact conditional mode is used only in analysis, not by the implementable algorithm.Source: §2.2 · pp. 5–6
- The scale η trades cheaper, better-conditioned conditional sampling against more local transitions that may require more Gibbs iterations.The proximal map is the mode of the conditional law, while the conditional draw is its stochastic analogue.Source: §2.2 · pp. 5–6
2.3 Divergences and notation
- The paper defines an order-2 Rényi warm-start bound as D2(P∥Q)≤Δ and total variation as sup_A |P(A)−Q(A)|.It also gives a coupling bound on total variation and defines eO(B) as B times a fixed universal power of logarithms in displayed problem parameters.Source: §2.3 · pp. 6–7
3 The Proximal Bouncy Particle Sampler
This section introduces Proximal BPS, replacing independent Gibbs updates with reflected and conditional half-turn updates that retain dependence on the current state. It establishes properties and convergence for the ideal sampler, combining a direct L2 estimate with a half-turn bound on corrector change.
3.1 Overview of the construction
- Proximal BPS reflects the auxiliary variable and mixes conditional resampling with a half-turn update; the dominant half-turn branch retains dependence on the current position.With probability ρ, conditional resampling prevents periodicity; otherwise, the half-turn update carries information across iterations.Source: §3.1 · pp. 8–9
3.2 The ideal conditional half-turn process
- The conditional half-turn process uses a reference-point gradient in its harmonic motion and bounces on the residual gradient, reducing the expected bounce count to order βηd.Its time-π horizon removes dependence on the fresh momentum while retaining the displacement from the harmonic equilibrium with its direction reversed.Source: §3.1 · p. 10, §3.2 · p. 11
- For each fixed y, the ideal conditional half-turn kernel is reversible with respect to πX|Y=y, and the ideal Proximal BPS transition preserves the augmented law πη.The conditional phase-space process is also stated to be well-defined and stationary.Source: §3.2 · p. 12, §3.3 · p. 13
3.3 Ideal Proximal BPS
- Theorem 3.5 gives a χ2 convergence guarantee when η ≤ c0/β, for initial augmented-state laws with finite χ2 divergence from πη.The theorem also provides a bound for the X-marginal; the supplied passages do not include the bound’s formula.Source: §3.3 · pp. 13–14
- The half-turn estimate bounds its contribution to corrector change, complementing the direct L2 estimate to control both components.The model-specific input is Proposition B.2.Source: §3.3 · p. 14
- The analysis uses a discrete modified-L2 method of Dolbeault, Mouhot, and Schmeiser, with a sharpened formulation by Fan, Li, and Lu.Source: §3.3 · p. 14
- The proof constructs a modified energy, proves its contraction, selects ρ and a block length, and derives χ2 and total variation bounds.The proof roadmap also assigns separate steps to the ordinary L2 estimate and the reflection and half-turn estimates.Source: §3.3 · p. 14
- L2 analysis is partly motivated by, which gave obstructions to entropic convergence for piecewise-deterministic Markov processes such as BPS.Source: §3.3 · p. 14
4 Efficient implementation of Proximal BPS
Section 4 turns ideal Proximal BPS into a gradient-only algorithm using an approximate conditional sampler and a capped half-turn simulation. Its optimized expected gradient complexity is eO(√κ d^1/4) for an O(1) order-2 Rényi warm start.
4.1 Implementing the restricted Gaussian oracle
- The conditional draw is implemented with a first-order rejection sampler using a constructed approximate proximal point that meets the required residual bound.Gradient descent constructs the point, and a gradient query certifies the stopping condition; function values and an exact proximal oracle are unnecessary.Source: §4 · p. 15, §4.1 · p. 15
4.2 Simulating the half-turn process
- The capped simulation exactly follows the half-turn process with rate min{λt, λ̄}, and Proposition 4.2 bounds its coupling error and gradient-query count.The cap need not dominate the ideal rate.Source: §4.2 · pp. 16–17
4.3 Implementable Proximal BPS and total complexity
- Algorithm 4.2 combines the approximate conditional draw and capped half-turn process, accounting for both approximations in its total-variation and query bounds.The approximate draw is reused as the next position on the resampling branch and as the reference point for the capped process.Source: §4.3 · p. 17
- In particular, when Δ = O(1), expected gradient complexity is eO(√κ d^1/4), where eO suppresses logarithmic factors in d, κ, and 1/ε.Source: §4.3 · p. 20
A Conditional half-turn process: construction, invariance, and implementation
The appendix analyzes the conditional half-turn process through its construction, invariance, and implementation. It uses an integrated-hazard construction and path reversal, then studies the averaged kernel and capped-rate implementation.
A Conditional half-turn process: construction, invariance, and implementation
- The analysis constructs the fixed conditional process with an integrated-hazard method and proves reversibility by path reversal.It then averages over initial momentum and the reference point to obtain the ideal half-turn kernel and its augmented-space lift.Source: A Conditional half-turn process:… · p. 21
A.1 The fixed-reference phase space process
For fixed conditioning and reference points, the paper constructs a continuous-time phase-space process from harmonic motion and bounce events. It establishes well-posedness and proves invariance through a path-reversal argument.
A.1 The fixed-reference phase space process
- The fixed-reference process follows harmonic flow between bounces, with bounce updates given by a reflection map.The flow and bounce preserve shifted harmonic energy, which is used to bound the rate along a path.Source: §A.1 · p. 21
- The exponential-clock recursion has no accumulating jump times, so it defines a unique process for all t ≥ 0 with the time-homogeneous Markov property.The rate is bounded along the path by a quantity determined by the shifted harmonic energy.Source: §A.1 · pp. 21–22
- Path reversal establishes invariance and momentum-reversed detailed balance, yielding a self-adjoint L2 semigroup.The pathwise argument handles possible discontinuities of the bounce map; Gaussian momentum averaging removes momentum reversal to produce a self-adjoint position kernel.Source: §A.1 · pp. 22–24
A.2 The averaged ideal half-turn kernel
Averaging the fixed-reference half-turn kernels over the conditional reference-point law gives a conditional kernel on positions. The paper verifies joint measurability and obtains a self-adjoint Markov contraction on the augmented space.
A.2 The averaged ideal half-turn kernel
- Averaging over the reference point, drawn independently from the conditional law, gives a self-adjoint Markov contraction for each fixed conditioning value.The averaging measure does not depend on the current position.Source: §A.2 · p. 24
- The conditional half-turn kernels admit a jointly Borel-measurable version, and their augmented-space operator is a self-adjoint Markov contraction.The proof obtains measurability from the Borel dependence of the flow, event rate, reflection map, and jump-time recursion.Source: §A.2 · p. 25
A.3 Bounce complexity and capped-rate implementation
The appendix bounds the capped simulation through a coupling with the ideal process, then develops an operator and micro–macro framework for proving contraction. Reflection transfers macroscopic observables into the microscopic subspace, while coercivity and a half-turn estimate quantify that transfer.
A.3 Bounce complexity and capped-rate implementation
- The ideal and capped paths agree until an event from the excess rate (λt − λ̄ξ,∆)+; the probability they separate is bounded by the expected number of such events.Source: §A.3 · p. 26
- Each arrival of the capped proposal clock before time π costs one gradient query, and the arrival count is Poisson with mean πλ̄ξ,∆.The reference gradient is cached.Source: §A.3 · p. 26
B.1 Operator framework and micro–macro decomposition
- The operator framework splits observables into functions of y and a mean-zero conditional component; reflection can transfer mass from the first subspace into the second.The macroscopic subspace consists of functions depending only on y.Source: §B.1 · pp. 27–28
B.2 Estimates for the reflection and half-turn operators
- For macroscopic observables, the squared norm of reflection’s microscopic component equals an expected conditional variance.The representation uses the exchangeable pair (Y+, Y−).Source: §B.2 · p. 28
- Under βη ≤ 1, Lemma B.1 gives uniform macroscopic coercivity, relating centered observables to the microscopic component created by reflection.The proof uses a Poincaré inequality for πY.Source: §B.2 · pp. 28–29
- When βη is sufficiently small, the half-turn acts nearly as −I on the microscopic directions produced by reflection.Source: §B.2 · p. 29
B.3 Modified L2 contraction and proof of the main theorem
The proof uses a modified L2 Lyapunov function to turn one-step decay into contraction and then a χ2 and total-variation bound. Its corrector tracks how reflection transfers control from the macroscopic component to the microscopic one.
B.3 Modified L2 contraction and proof of the main theorem
- The corrector records reflection’s transfer of the macroscopic component into the microscopic subspace, allowing the contraction argument to control both components.Source: §B.3 · p. 29
- The modified Lyapunov function is first shown equivalent to squared L2 norm, then shown to decay in one step.Source: §B.3 · pp. 30–31
- A fixed block of steps converts one-step Lyapunov decay into L2 contraction, which is then applied to densities.Source: §B.3 · p. 33
- The proof obtains the stated χ2 estimate and uses the X-marginal μ of πη, together with contraction of total variation under projection, to reach the target-space bound.Source: §B.3 · p. 34
C Macroscopic coercivity and the half-turn bound
- The two model-specific inputs are macroscopic coercivity and a global half-turn bound.Source: C Macroscopic coercivity and the… · p. 34
C.1 Proof of Lemma B.1
The proof of Lemma B.1 analyzes reflection on macroscopic observables using conditional expectations and covariance identities. It relates the microscopic component’s norm to the residual from the macroscopic reflection operator.
C.1 Proof of Lemma B.1
- The microscopic component’s squared norm is the difference between the squared norm of f and that of its macroscopic reflection.Source: §C.1 · p. 35
- The gradient of the macroscopic reflection is a conditional expected gradient minus a conditional covariance with ∇V(X).Source: §C.1 · pp. 36–37
- The proof extends the gradient identity beyond smooth functions by establishing continuity under L2 approximation, using conditional Poincaré and Cauchy–Schwarz bounds.Source: §C.1 · p. 37
C.2 The half-turn operator on affine functions
The conditional half-turn is analyzed in normalized coordinates, where harmonic motion reverses affine observables at time π. A Duhamel representation bounds departures from that reversal, yielding a uniform estimate for centered affine functions.
C.2 The half-turn operator on affine functions
- Under βη ≤ 1/2, Proposition C.2 states a weighted half-turn estimate for every affine function ℓa,y(w)=a⊤(w−¯w_y), uniformly in y.Source: §C.2 · p. 38
- For a Gaussian conditional law, the bounce rate vanishes and harmonic motion maps each centered affine observable to its negative at time π.Source: §C.2 · p. 38
- The Duhamel formula expresses deviation from sign reversal through an inhomogeneous term that captures the mismatch of centers and residual gradient bounces.Source: §C.2 · pp. 38–39
- The proof bounds the conditional gradient covariance in terms of I − Σy, which quantifies the covariance defect relative to a standard Gaussian.Source: §C.2 · pp. 38, 40
- The conditional estimate extends to the global assertion by integrating over y, with the stated integrability condition ensuring the observable lies in L2(πη).Source: §C.2 · p. 40
C.3 Proof of Proposition B.2: From affine estimates to the half-turn bound
The proof controls the UPP spectrum near +1 by combining conditional affine approximation and energy estimates; the appendix also states Hilbert-space, Markov-operator, and strong-log-concavity tools used in the analyses.
C.3 Proof of Proposition B.2: From affine estimates to the half-turn bound
- The conditional decomposition has a q_y-mean-zero residual whose squared energy is bounded by 4η times the conditional gradient variance.The affine coefficient is chosen to center the residual gradient.Source: §C.3 · pp. 41–42
- For sufficiently small universal βη, the proof absorbs the remaining terms to establish the stated energy estimate.The covariance decomposition separates an affine contribution from a residual covariance term.Source: §C.3 · pp. 42–43
- On the UPP spectral subspace [1/2,1], the proof combines the affine part and controlled error to obtain (C.14).It invokes Proposition C.2 for the affine part and Lemma C.6 for the error estimate.Source: §C.3 · pp. 41, 44
D.1 Hilbert space and Markov operator conventions
- A Markov operator preserving ζ is an L2(ζ) contraction, while its adjoint evolves densities and both operators fix constants.The appendix also defines Hilbert-space adjoints, projections, and spectral projections.Source: §D.1 · pp. 44–45
D.2 Strong log-concavity and covariance bounds
- The appendix invokes strong-log-concavity tools and applies a two-sided covariance bound to the conditional potential to obtain the augmented-marginal Hessian estimate (2.13).The covariance-bound proof uses Brascamp–Lieb for the upper bound and Cramér–Rao for the lower bound [BGL14, Chapter 4].Source: §D.2 · p. 45
D.3 Comparison of positive operators
The appendix gives a comparison result for positive operators under a vectorwise norm condition and notes that its logarithmic dependence on the operator range is sharp.
D.3 Comparison of positive operators
- The comparison lemma assumes positive semidefinite self-adjoint A and B, bounded B with mI ⪯ B ⪯ MI, and ∥Av∥² ≤ ∥v∥∥Bv∥ for every v.The passage defines K from the bounds on B but does not show the lemma's conclusion.Source: §D.3 · p. 46
- The logarithmic dependence on K in the positive-operator comparison lemma is sharp in general.The appendix states that the best possible constant is at least max{1, (log K)/π}.Source: §D.3 · p. 46
- Remark D.5 takes K > 1 and H = L2([1,K], dx), defines A by an integral and B by Bf(x)=xf(x), and states Carlson’s inequality yields Lemma D.4’s hypotheses.For f(x) := 1/x, the passage reports ∥Af∥2 = (log K)2/π.Source: §D.3 · p. 46