Smoothed Picard Hamiltonian Monte Carlo

1 Introduction

• Gaussian smoothing, Picard iteration, and Chebyshev–Lobatto quadrature yield a stable unadjusted Hamiltonian step with only two derivatives assumed for the original potential.The method combines smoothing with a high-order discretization of Hamiltonian dynamics.Source: §1 · pp. 2–3
• The low-accuracy smoothed Picard HMC sampler has a Wasserstein guarantee for sampling from the target.The introduction states a W2 guarantee, but the supplied theorem text does not include its full bound.Source: §1 · p. 3
• d^1/6 is the dimension dependence achieved by Theorem 1.1, compared with prior W2 approaches having d^1/3 dependence.The comparison passage also describes other prior rates, including a KL rate and an exact-simulation complexity.Source: The expected number of gradient and… · p. 3

The expected number of gradient and proximal queries is at most

• The warm start can initialize proximal BPS, which has a stated complexity of eO(κ^1/2 d^1/4 polylog(1/ε)) from a warm start.The paper identifies providing an implementable warm start for proximal BPS as an application of smoothed Picard HMC.Source: The expected number of gradient and… · p. 3
• The recursive warm start generator wraps smoothed samplers to upgrade W2, Wq, and Wψ2 guarantees into stronger divergence guarantees.The framework gives KL, mixed TV and Rényi, or genuine Rényi guarantees, respectively, and is applied to smoothed Picard HMC.Source: The expected number of gradient and… · pp. 3, 5
• The generator also removes the κ^2 term from the complexity, while Theorem 7.1 relaxes strong log-concavity to a log-Sobolev inequality and provides a KL guarantee.Source: The expected number of gradient and… · p. 3

2 Preliminaries

• Except for the generic warm start results, the paper assumes V is C2 and its Hessian lies between κ^-1 I and I after normalizing β = 1.The general Hessian bounds are αI ⪯ ∇2V ⪯ βI, with rescaling used for the normalization.Source: §2 · p. 5
• The exposition initially permits proximal queries as well as gradient queries, but the paper later removes proximal queries altogether.Needed proximal calls use step sizes at most 1/(2β), for which computing the proximal operator is described as well-conditioned.Source: §2 · p. 5
• A restricted Gaussian oracle is the sampling analogue of the proximal operator, and the paper gives an approximate implementation using eO(1 + (a^2d)^1/3) queries.The implementation is based on the exact diffusion sampler.Source: §2 · p. 5
• The paper uses eO(B) for B multiplied by a universal logarithmic power in dimension, condition number, inverse accuracy, and any explicit confidence parameter.Source: §2 · p. 6

3.1 The smoothing trick

• Gaussian smoothing lets the method sample from πη = π ∗ N(0, ηI) and recover an approximate sample from π through the RGO.The smoothed potential has higher-order derivative bounds, while its gradient is not explicitly computable from the conditional-expectation definition.Source: §3.1 · p. 6

3.2 The higher-order integrator and the W2 guarantee

• Choosing h ≍ d^-1/6 yields a non-dominant eO(d^1/9)-query terminal RGO step at η ≍ d^-1/3; polylogarithmic Chebyshev–Lobatto nodes suffice when h < √η.Picard iteration approximates the dynamics, and the overview attributes the dimensional scaling to accumulated local errors and refined quadrature control.Source: §3.2 · pp. 7–9

3.3 Higher moment control

• Smoothed Picard HMC can achieve Wq closeness to the smoothed target for any q ≥ 2 under the stated strong-convexity, smoothness, and reference-point conditions.The overview notes that the Wq extension costs polynomially in q, which is harmless when q is logarithmic in the problem parameters.Source: §3.3 · pp. 9–10

3.4 Recursive warm start generation

• A recursive warm start generator upgrades smoothed-sampler W2, Wp, or Wψ2 guarantees to KL, mixed TV and Rq, or Rq guarantees, respectively.The framework is generic and applies to any smoothed sampler implementing its interface.Source: §3.4 · pp. 10–11
• The recursive warm-start construction adds Gaussian noise so heat-flow regularization and reverse transport inequalities can upgrade Wasserstein guarantees to stronger distances.Instantiating the framework with smoothed Picard HMC establishes a warm-start guarantee.Source: §3.4 · p. 11

4.1 Setup and smoothing estimates

• Gaussian smoothing yields a smoothed potential with bounded Hessian, so its gradient gη is 1-Lipschitz; near-Gaussian cancellation supports higher-order regularity estimates.The cancellation is identified as stronger than a generic moment bound for a strongly log-concave law.Source: §4.1 · pp. 12–13

4.2 Local error analysis

• For h ≤ c/κ and J ≥ 2, the local-error theorem gives one-step contraction and N-step exponential contraction, with additive variance, bias, and discretization errors.The bounds are stated in a κ-dependent phase-space metric; the three errors arise from centered stochastic-gradient fluctuations, stochastic-gradient bias, and numerical-kernel discretization.Source: §4.2 · pp. 14–15

4.3 Properties of the exact-gradient kernel

• When h ≤ c/κ, the exact-gradient kernel contracts in a κ-dependent phase-space metric, and a synchronous coupling controls differences in its increments.These kernel properties supply the coupling conditions used in the local-error analysis.Source: §4.3 · pp. 16–17

4.4.1 Bias bound

• The bias analysis bounds the displacement between conditional-mean and exact-gradient kernels, then propagates those gradient errors through the endpoint update.The endpoint stochastic-gradient block is independent of the first block, and the resulting bound is uniform in the phase-space input and Gaussian variables.Source: §4.4.1 · pp. 17–18

4.4.2 Variance bound

• For h² ≤ 1, Lemma 4.10 states Lipschitz bounds that hold uniformly in the input, supporting a uniform Lipschitz conclusion under that condition.Source: §4.4.2 · p. 18
• Lemma 4.11 gives a second-order variance bound quadratic in A and B by exploiting Gaussian smoothing and conditionally centered perturbations.Unlike the basic coupling bound, it uses that Z is Gaussian; its AB term accounts for dependence of the perturbation on the smoothing variable.Source: §4.4.2 · p. 19

4.4.3 Discretization error bound

• Stationary iterated-derivative estimates control discretization error for step sizes satisfying h^2 ≤ 1.The derivative estimates use Sobolev divergence and Gaussian gradient bounds to control the Hamiltonian generator.Source: §4.4.3 · pp. 20–22

4.5 Proof of Theorem 1.1

• The initialization lemma bounds the initial law’s qth moment for every q ≥ 2.The initial law is a point mass at x_ref paired with a standard Gaussian momentum.Source: §4.5 · p. 23
• The terminal RGO preserves π after smoothing and contracts center-law W2 distances by a factor κ/(κ + η) ≤ 1.These properties transfer a center-law approximation to the target distribution.Source: §4.5 · p. 23
• The proof establishes W2(bμ_N R_η, π) ≤ √κ ε, with terminal-routine cost dominated by the Picard HMC bound.The error bookkeeping is uniform over 0 < ε ≤ 1.Source: §4.5 · p. 24

5.1 Wq guarantee

• The Wq guarantee uses a smoothing variance and a gradient-only implementation of Picard HMC, with the result targeting the smoothed distribution πη.The section motivates Wq bounds as needed for the warm-start application; unlike Theorem 1.1, this guarantee omits the terminal FORS step.Source: §5.1 · p. 25, §5 · p. 25, Its expected number of gradient… · p. 25

5.2 Wq local error analysis

• For q > 2, the local-error analysis replaces the W2 square-expansion framework with an Lq smoothness inequality.Its kernel comparison separates centered variance and bias errors, then combines them with weak and strong one-step errors.Source: §5.2 · pp. 25–26
• The Wq local-error bound follows by adapting the W2 argument with the Wq local-error framework, retaining contraction and error controls.The stated step-size condition is h ≤ c/(κq).Source: §5.2 · p. 26

5.3.1 Variance bound

• The variance analysis extends a second-order W2 estimate to Wq using a Gaussian divergence inequality.The text says this controls the variance term in Wq and states a variance bound for h ≤ 1.Source: §5.3.1 · p. 27

5.3.2 Discretization error bound

• The discretization analysis obtains dimension-free divergence bounds for polynomial momentum fields, then uses stationary iterated-derivative bounds to control discretization error.The divergence-bound constant is stated to be independent of tensor dimensions.Source: §5.3.2 · pp. 27–28

5.4 Proof of Theorem 5.1

• The proof first analyzes exact proximal queries, then removes proximal access to obtain an expected gradient-query implementation.The proof passage gives separate query counts for direct gradient and exact proximal evaluations before the modifications that remove proximal queries.Source: §5.4 · p. 29

6.1 Main result

• Theorem 6.1 upgrades smoothed-sampler Wasserstein guarantees to stronger guarantees with negligible overhead, and in case (ii) needs W_p accuracy only for p logarithmic in 1/δ.The theorem assumes a β-log-smooth target satisfying a logarithmic Sobolev inequality and an admissible reference point; its proximal wrapper calls the smoothed sampler only on distributions with universally bounded condition number.Source: §6.1 · pp. 30–31

6.2.1 Regularization and reverse transport inequalities

• Heat-flow regularization enables reverse transport inequalities that convert transport-distance bounds into stronger information-theoretic divergence bounds.A truncation lemma first bounds displacement outside an event of probability at most δ, enabling stronger guarantees after Gaussian convolution.Source: §6.2.1 · pp. 31–32

6.2.2 RGO calculus

• An RGO applied to another RGO remains an RGO for the original distribution, with updated parameters governing its strong convexity, smoothness, and condition number.The parameter updates are used in the recursive construction.Source: §6.2.2 · p. 32

6.3 Recursive RGO sampler

• The recursive RGO sampler has guarantees under each of the three smoothed-sampler interfaces and can be analyzed as stages that smooth, update RGO parameters, and terminate with FORS.With bounded smoothing levels, the stages make geometric progress, including contraction in the well-conditioned regime.Source: §6.3 · pp. 33–34
• Inner reference-point construction costs at most Cκ0Lq expected gradient queries across the recursive stages.In the large-conditioning regime, condition numbers decrease geometrically; later stages have universally bounded condition number.Source: §6.3 · p. 36
• Terminal reference-point construction costs at most Cκ0Lq expected gradient queries, and each terminal call contributes at most ∆0²/(J + 1) in the relevant divergence.The terminal FORS call has constant expected gradient-query cost, establishing the base cases of the backward induction.Source: §6.3 · p. 36

6.4 Outer proximal sampler

• The outer proximal chain uses step size (2β)^−1 and applies the recursive RGO sampler only to uniformly well-conditioned RGO targets, avoiding the inherited κ² query cost.Each outer RGO target has curvature between βI and 3βI, so its condition number is at most 3.Source: §6.4 · pp. 36–37
• The initial Gaussian distribution has Rényi divergence at most Cd log(eκ) from π for every q ≥ 1, with R1 interpreted as KL.Source: §6.4 · p. 37
• Outer reference-point gradient descent has logarithmic expected work per step under A1 and A3, while under A2 it is capped to control failure probability and ensure TV error at most δ.Under A2, the cap-hit probabilities sum to at most δ/2, and TV errors from approximate RGO kernels also sum to at most δ/2.Source: §6.4 · p. 38

7.1 Picard HMC warm starts

• Picard HMC supplies warm-start guarantees under β-log-smoothness, a logarithmic Sobolev inequality with constant α^-1, and a reference point satisfying V(xref) − inf V ≤ d.The section also states that the KL warm-start guarantee holds under more general assumptions.Source: §7.1 · pp. 38–39

7.2 From the proxy warm start to high accuracy

• The high-accuracy construction runs proximal BPS from the implementable warm start and uses the triangle inequality to reach error ε in total variation.The warm-start phase has expected cost at most Cκ^7/6 d^1/6 L^7, with L := log(eκd/ε).Source: §7.2 · p. 40

A Terminal RGO implementation by FORS

• The appendix reviews the FORS mechanism and describes a variant tailored to sampling from RGO distributions, emphasizing differences from prior analyses.The supplied passage identifies [CCDR26] [CCRZ26] as the prior analyses.Source: A Terminal RGO implementation by FORS · p. 41

Proposal and clipping

• Given a proposal law Q and an auxiliary-variable estimator W with mean w, clipping at level B defines the implemented path law.The desired path law is specified relative to Q through its log-weight.Source: §A.1 · p. 41

Poisson exponentiation

• Poisson exponentiation samples by drawing a Poisson count, generating that many conditionally independent bounded estimator copies, and accepting according to their combined value.The lemma assumes X is drawn from Q and the conditional estimator values lie in [−B, B].Source: §A.1 · p. 41

First-order rejection sampling

• The accepted proposal has law proportional to exp(w)Q, and each attempt is accepted with probability at least exp(−2B).The supplied lemma describes this as the acceptance probability of one attempt.Source: §A.1 · p. 41

A.2 FORS for RGO distributions

• The FORS routine proposes successive blocks, accepts endpoints, restarts when required by the stopping set, and returns the state after N accepted blocks.Each attempt initializes independently and uses Poisson-distributed proposal times.Source: §A.2 · p. 42
• Exactly integrating the linear drift makes the Girsanov control depend only on the residual force.The proposal is a variant of FORS diffusion simulation for the specified RGO distributions.Source: §A.2 · pp. 42–43
• The block count satisfies N ≤ C L^(5/3){1 + (a^2d)^(1/3)}, with the stopped implementation using O(N) expected queries and O(N + log(1/p)) queries with probability at least 1 − p.The stated parameter choices satisfy both the cubic and quadratic restrictions on the block length.Source: §A.2 · p. 44

B Chebyshev–Lobatto interpolation and quadrature

• The appendix records a Banach-valued interpolation error bound for functions in C^J([0, h]; E).The interpolation facts are used by the Picard integrator.Source: B Chebyshev–Lobatto interpolation and… · p. 44

B.1 Nodes, interpolating polynomials, and integrated weights

• The Chebyshev–Lobatto interpolant is the unique polynomial of degree at most J − 1, and the Clenshaw–Curtis weights are nonnegative and sum to h.The appendix also states coefficient bounds for these interpolation and quadrature weights.Source: §B.1 · pp. 44–45

B.2 Interpolation error bound

• The Banach-valued interpolation error bound applies to f ∈ C^J([0, h]; E) for a Banach space E.The cited passage identifies it as the Banach-valued Genocchi–Hermite interpolation remainder.Source: §B.2 · p. 46

C.1 Basic inequalities and heat smoothing

• The Gaussian-analysis appendix gives Gaussian Sobolev inequalities and derivative and divergence bounds for functions and vector fields under Gaussian measure.Its setup includes a finite-dimensional Hilbert target and defines Gaussian divergence as the gradient's L2(γ)-adjoint.Source: §C.1 · p. 47
• The authors expect the results to extend to infinite dimensions but say that setting is not needed here.The stated Gaussian-analysis setup uses a finite-dimensional Hilbert space.Source: §C.1 · p. 47
• The heat semigroup satisfies ∂tPtF = 1/2 ΔPtF, and differentiation commutes with it for weakly differentiable F.Gaussian integration by parts also gives a derivative formula.Source: §C.1 · p. 48

C.2 Cumulants and derivatives of smoothed potentials

• For the exponential family Rη,y, derivatives of log Zη give its cumulants of order k ≥ 2.The appendix introduces a recursive representation of cumulants through iterated Gaussian covariance.Source: §C.2 · p. 49
• The near-identity Gaussian-image cumulant bound assumes ∥DT∥op ≤ 1 and ∥DT − I∥op ≤ εT, with 0 ≤ εT ≤ 1.The cited proof bounds the relevant tensor quantities using the deviation of DT from identity.Source: §C.2 · pp. 49–50

C.3 Proofs of the second-order Wasserstein estimates

• The evolution of Law(LXt) is written as a continuity equation with an explicit velocity field.The proof constructs the interpolating variables and uses Gaussian divergence identities to obtain this representation.Source: §C.3 · pp. 50–51
• The velocity estimate contains terms proportional to B^2, AB t^-1/2, and AB (1 − t)^-1/2, with operator-norm factors.The two endpoint singularities arise from Gaussian-divergence rescaling and differentiating the heat extension.Source: §C.3 · p. 52

C.4 Gaussian polynomial bounds

• For polynomial vector fields of degree at most r, the Gaussian divergence Lq norm is bounded by C(√q + q√r) times the field's Lq norm.This is the stated Lq Gaussian polynomial estimate for q ≥ 2.Source: §C.4 · p. 53

D.1 Approximate proximal queries

• When h^2 ≤ 1, the one-step Wq error from an εprox-accurate proximal center is at most C h εprox for q ≥ 1.The proof bounds the position error by Ch^2 εprox and the momentum error by Ch εprox.Source: §D.1 · p. 53

D.2 Gradient-only implementation

• The gradient-only implementation replaces each approximate proximal query in Picard HMC with a gradient-descent routine.The routine starts at y and returns an iterate as the approximate proximal point.Source: §D.2 · p. 54
• The map Ty := y − η∇V is an η-contraction with fixed point x+, and the iteration's successive-step difference decreases geometrically.The stated bound is ∥x^(k+1) − x^k∥ ≤ η^(k+1)∥∇V(y)∥; the certificate is valid.Source: §D.2 · p. 54
• Expected complexity analysis uses the fact that centers encountered along the realized run have controlled second moments.The pointwise work bound depends on the queried center.Source: §D.2 · p. 54

• Each phase uses 2J proximal and 2J direct gradient calls.Source: Consequently the expected total number… · p. 55
• J ≤ C L_q and N ≤ CκL_q/h yield the query bound in (D.8).Source: Consequently the expected total number… · p. 55
• Replacing each proximal query with Algorithm D.2 at tolerance ε_prox gives valid ε_prox-accurate calls and the claimed W_q and query bounds of Theorem 5.1.The query bound follows from Lemma D.4 after substituting the main-text step size.Source: Consequently the expected total number… · p. 55

E Extension to arbitrary Picard depth

• The appendix records the discretization error of Picard HMC with the exact smoothed gradient at arbitrary depth.Source: E Extension to arbitrary Picard depth · p. 55
• A phase uses the smoothed gradient g_η = ∇V_η, the law Π_η = π_η ⊗ N(0, I), and independent Gaussian refresh variables.Source: E Extension to arbitrary Picard depth · p. 55
• Under κ ≥ 1, 0 < η ≤ 1, and κ^-1 I ⪯ ∇²V ⪯ I, Theorem E.2 treats every integer K ≥ 1 and q ≥ 2; K = 2 recovers Proposition 5.9.Source: E Extension to arbitrary Picard depth · p. 56
• The analysis couples inputs with the same Gaussian variables and uses a map with Lipschitz constant at most h², then bounds discretization error from the exact trajectory and Picard iteration.The quadrature rules integrate the interpolating polynomial exactly, and the momentum estimate dominates the resulting bound.Source: E Extension to arbitrary Picard depth · pp. 56–57