Paper deep dive
Reflected Proximal Dynamics for Accelerated High-Accuracy Log-Concave Sampling
Bo Yuan, Wei Guo, Yongxin Chen
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We study sampling from $\pi(x)\propto\mathrm{e}^{-V(x)}$ on $\mathbb{R}^d$, where $V$ is $\alpha$-strongly convex and $\beta$-smooth with condition number $\kappa=\beta/\alpha$. For this problem, the proximal sampler preserves the target exactly and recovers overdamped Langevin dynamics in the small-step limit, but its Gaussian resampling step discards directional information. Underdamped Langevin dynamics retains such information through momentum to accelerate exploration, motivating a proximal construction that combines target preservation with directional persistence. We introduce \textit{reflected proximal dynamics}, which combines proximal sampling with an exact reflection of the auxiliary variable, conditional kinetic motion, and controlled refreshment. We develop two realizations: a flux-resampled underdamped half-turn and a persistent-momentum quarter-turn. From a warm start with constant $\chi^2$ divergence, both methods achieve total variation error $\varepsilon$ using $\widetilde O(\sqrt{\kappa}d^{1/4})$ expected oracle queries. We further remove the supplied warm-start assumption by combining either sampler with existing posterior samplers in a proximal initialization stage with an optimized coarse step size. Accounting for both initialization and sampling yields $\widetilde O(\kappa^{6/7}d^{1/4})$ expected oracle queries. In both regimes, the dependence on the target accuracy $\varepsilon$ is $\mathrm{polylog}(1/\varepsilon)$, making reflected proximal dynamics a promising framework for accelerated high-accuracy log-concave sampling.
Tags
Links
- Source: https://arxiv.org/abs/2610.16248v1
- Canonical: https://arxiv.org/abs/2610.16248v1
Trouble viewing inline? Open PDF directly →
Full Text
212,604 characters extracted from source content.
Expand or collapse full text
Reflected Proximal Dynamics for Accelerated High-Accuracy Log-Concave Sampling Bo Yuan Wei Guo Yongxin Chen Affiliation: Georgia Institute of Technology Affiliation: byuan48, wei.guo, yongchen@gatech.edu Abstract We study sampling from π(x)∝e−V(x)π(x) ^-V(x) on ℝdR^d, where V is α-strongly convex and β-smooth with condition number κ=β/ακ=β/α. For this problem, the proximal sampler preserves the target exactly and recovers overdamped Langevin dynamics in the small-step limit, but its Gaussian resampling step discards directional information. Underdamped Langevin dynamics retains such information through momentum to accelerate exploration, motivating a proximal construction that combines target preservation with directional persistence. We introduce reflected proximal dynamics, which combines proximal sampling with an exact reflection of the auxiliary variable, conditional kinetic motion, and controlled refreshment. We develop two realizations: a flux-resampled underdamped half-turn and a persistent-momentum quarter-turn. From a warm start with constant χ2χ^2 divergence, both methods achieve total variation error ε using O~(κd1/4) O( κd^1/4) expected oracle queries. We further remove the supplied warm-start assumption by combining either sampler with existing posterior samplers in a proximal initialization stage with an optimized coarse step size. Accounting for both initialization and sampling yields O~(κ6/7d1/4) O(κ^6/7d^1/4) expected oracle queries. In both regimes, the dependence on the target accuracy ε is polylog(1/ε)polylog(1/ ), making reflected proximal dynamics a promising framework for accelerated high-accuracy log-concave sampling. 1 Introduction We study sampling from a distribution π(x)∝e−V(x)π(x) ^-V(x) on ℝdR^d, where V is α-strongly convex and β-smooth, with condition number κ=β/ακ=β/α. This problem appears in fields including Bayesian inference, machine learning, and statistical physics (Gelman et al., 1995; Ghosal and van der Vaart, 2017; Newman and Barkema, 1999; Chewi, 2023), and a classical approach is Markov chain Monte Carlo (MCMC). An MCMC method constructs a Markov chain with stationary distribution π and, under suitable conditions on the transition kernel and target distribution, converges to π. Such samples allow expectations under π to be approximated even when direct integration and the normalizing constant of π are unavailable. Two workhorse continuous-time processes are (overdamped) Langevin dynamics (LD) and underdamped Langevin dynamics (ULD). LD is a gradient flow of the potential energy function V combined with Brownian noise, and by Euler discretization, it leads to the Langevin Monte Carlo (LMC, Dalalyan (2017); Durmus et al. (2019); Vempala and Wibisono (2019); Chewi et al. (2022)) algorithm. ULD introduces momentum, allowing the process to retain its direction rather than move diffusively (Cheng et al., 2018; Dalalyan and Riou-Durand, 2020; Ma et al., 2021; Betancourt and Girolami, 2019). This persistence can accelerate exploration, but discretizing the coupled position–momentum dynamics requires more sophisticated numerical schemes (Shen and Lee, 2019; Zhang et al., 2023). The proximal sampler is a discrete-time alternative to LD that preserves the target exactly (Lee et al., 2021; Chen et al., 2022). It augments the position X∼πX π with an auxiliary variable Y∼(X,ηId)Y (X,η I_d), then alternates a Gaussian draw of Y|XY X with a localized draw of X|YX Y, called a restricted Gaussian oracle (RGO) draw. In the η→0η→ 0 limit, the chain reduces to LD (Mitra and Wibisono, 2025). This construction also exposes a gap: redrawing the Gaussian displacement Y−XY-X at every iteration discards directional memory, whereas ULD retains precisely this information through momentum. We bridge this gap with reflected proximal dynamics. Instead of redrawing the auxiliary variable, we set Y←2X−Y← 2X-Y, which leaves the position fixed and reverses the current displacement. Because the Gaussian displacement is symmetric about zero, this reflection preserves the joint distribution of (X,Y)(X,Y). Reflection alone does not move X, so most iterations follow it with a conditional transition that turns the retained displacement into persistent motion; occasional independent conditional draws remove correlations that this persistent update may leave behind. The resulting framework (8) brings the directional persistence of ULD into a target-preserving construction. We instantiate this framework with two implementable realizations: a flux-resampled half-turn (Sec. 4) and a persistent quarter-turn (Sec. 5). The former uses a damped conditional half-turn and then discards its temporary momentum; the latter rotates position information into momentum and carries it across outer iterations. Both use the same reflection–refresh framework but store directional information differently. When a warm start is not supplied, we prepend a proximal warm-up with an optimized coarse step size and account for its cost together with final sampling (Sec. 6). In either regime, geometric contraction makes the ideal mixing horizon logarithmic in 1/ε1/ , and the implementable complexity depends on the target accuracy only through polylogarithmic factors. Concurrent work on the proximal bouncy particle sampler (PBPS, Chen et al. (2026b)) uses the same Gaussian augmentation and outer reflection and matches our warm-start oracle complexity. Its conditional dynamics is an undamped harmonic half-turn with deterministic momentum reflections, whereas ours are a damped flux-resampled half-turn and a persistent position–momentum quarter-turn. Unlike PBPS and our flux method, the persistent method uses an independent RGO reference only on rare refresh steps; App. A gives a detailed comparison. We summarize our main contributions as follows: • We introduce reflected proximal dynamics and prove that it preserves the target distribution if the accelerated branch preserves the conditional distribution sampled by the RGO (Sec. 3). • We provide two realizations of the framework and prove accelerated contraction: both achieve O~(κd1/4) O( κd^1/4) expected oracle complexity from a constant-χ2χ^2 warm start (Secs. 4 and 5). • We remove the warm start requirement with an optimized warm-up stage, yielding an end-to-end expected oracle complexity of O~(κ6/7d1/4) O(κ^6/7d^1/4) (Sec. 6). 2 Preliminaries Notations. We write Law(X)Law(X) for the law of a random variable X. For probability densities ν and π, write the total variation distance as TV(ν,π):=supA|ν(A)−π(A)|TV(ν,π):= _A ν(A)-π(A) . When ν≪πν π, the order-r Rényi divergence for r>1r>1 is Dr(ν∥π):=1r−1log∫(ν(x)π(x))rπ(x)dxD_r(ν\|π):= 1r-1 ( ν(x)π(x))^rπ(x)\,dx. In particular, the χ2χ^2 divergence is χ2(ν∥π):=∫(ν(x)π(x)−1)2π(x)dx=eD2(ν∥π)−1χ^2(ν\|π):= ( ν(x)π(x)-1)^2π(x)\,dx=e^D_2(ν\|π)-1, and TV(ν,π)≤12χ2(ν∥π)TV(ν,π)≤ 12 χ^2(ν\|π). We call ν a constant-D2D_2 warm start of π if D2(ν∥π)=O(1)D_2(ν\|π)=O(1). The notations O~(⋅),Θ~(⋅) O(·), (·) suppress polylogarithmic factors in the displayed problem parameters. [a]+:=maxa,0[a]_+:= \a,0\. Assumption 1 (Target distribution). We assume that the target density is π(x)∝e−V(x)π(x) ^-V(x), where the potential V:ℝd→ℝV:R^d is α-strongly convex and β-smooth, i.e., α2∥x−x′∥2≤V(x)−V(x′)−⟨∇V(x′),x−x′⟩≤β2∥x−x′∥2 α2 x-x ^2≤ V(x)-V(x )- ∇ V(x ),x-x ≤ β2 x-x ^2 for all x,x′∈ℝdx,x ^d. We denote the condition number by κ:=β/ακ:=β/α. Query conventions. We assume oracle access to V and ∇V∇ V, and let NqryN_qry denote the total number of oracle evaluations used by an algorithm. This count includes rejected proposals, event-rate evaluations, mode-finding steps, and all work internal to an approximate RGO. 2.1 Overdamped and Underdamped Langevin Dynamics For d-dimensional standard Brownian motion (Bt)t≥0(B_t)_t≥ 0, (overdamped) Langevin dynamics (LD) follows dXt=−∇V(Xt)dt+2dBt,dX_t=-∇ V(X_t)\,dt+ 2\,dB_t, whereas underdamped Langevin dynamics (ULD) introduces momentum and follows dXt=PtdtanddPt=−∇V(Xt)dt−γPtdt+2γdBt.dX_t=P_t\,dt _t=-∇ V(X_t)\,dt-γ P_t\,dt+ 2γ\,dB_t. Their stationary distributions are π and π⊗γdπ _d, respectively, where γd=(0,Id) _d=N(0,I_d). The distinction relevant to this work is how the two dynamics move through space. In LD, Brownian noise acts directly on position, so motion over short time scales is diffusive and carries no separate directional state from one interval to the next; ULD instead stores direction in the momentum: before friction substantially changes PtP_t, the position moves approximately as Xt+h≈Xt+hPtX_t+h≈ X_t+hP_t. This short-time transport can cross a region more coherently than a position-only diffusion. More precisely, let μt=Law(Xt) _t=Law(X_t) for LD and μ^t=Law(Xt,Pt) μ_t=Law(X_t,P_t) for ULD. For initial laws with finite χ2χ^2 divergence from their respective stationary distributions, under Assump. 1, there are universal constants c,C>0c,C>0 such that χ2(μt∥π)≤e−2αtχ2(μ0∥π),χ2(μ^t∥π⊗γd)≤Ce−cαtχ2(μ^0∥π⊗γd),χ^2( _t\|π) ^-2α tχ^2( _0\|π),χ^2( μ_t\|π _d)≤ Ce^-c αtχ^2( μ_0\|π _d), where the ULD bound uses the low-friction choice γ≍αγ α (Vempala and Wibisono, 2019; Chewi, 2023; Chen et al., 2026d). Thus, after normalizing β=1β=1, the two continuous-time complexities are respectively O(κ)O(κ) and O(κ)O( κ). 2.2 The Proximal Sampler The proximal sampler (Lee et al., 2021; Chen et al., 2022) augments X∼πX π with an auxiliary random variable Y under the joint distribution Πη(x,y)∝exp(−V(x)−∥x−y∥2/(2η)), _η(x,y) (-V(x)- x-y ^2/(2η)), where η>0η>0 is the step size. Its X-marginal is π, and the conditional distributions are Y∣X=x∼(x,ηId),X∣Y=y∼qy(x)∝exp(−V(x)−∥x−y∥2/(2η)). Y X=x (x,η I_d), X Y=y q_y(x) (-V(x)- x-y ^2/(2η)). (1) An exact independent draw from qyq_y is called an ideal restricted Gaussian oracle (RGO) draw. The standard proximal sampler is the Gibbs chain that draws Yk+1∼(Xk,ηId)Y_k+1 (X_k,η I_d) and then Xk+1∼qYk+1X_k+1 q_Y_k+1. qyq_y has condition number (β+1/η)/(α+1/η)(β+ 1η)/(α+ 1η), which approaches 11 as η→0η→ 0. The small-η limit recovers LD, while the ideal Gibbs chain preserves π for every η (Mitra and Wibisono, 2025). 2.3 Poisson Process and Poisson Thinning A Poisson process is a random clock that determines when events occur. Its intensity (λt)t≥0( _t)_t≥ 0 describes how frequently the clock should ring near time t. Conditional on the history up to time t, an event occurs during the next short interval [t,t+dt)[t,t+dt) with probability λtdt+o(dt) _t\,dt+o(dt); equivalently, the next event time τ after t satisfies Pr(τ>s∣(λr)r∈[t,s])=exp(−∫tsλrdr) (τ>s ( _r)_r∈[t,s])= (- _t^s _r\,dr), ∀s≥t∀ s≥ t. However, directly sampling this clock can be difficult because the integrated intensity ∫tλrr _t _r\,dr generally has no closed-form inverse. Poisson thinning avoids this inversion. Suppose λs≤λ¯ _s≤ λ for all s≥ts≥ t. One first proposes candidate times from a homogeneous Poisson process of rate λ¯ λ, by sampling the waiting time to the next event from the exponential density u↦λ¯e−λ¯u1u>0u λe^- λu1_u>0. At a candidate time s, accept it as an event with probability λs/λ¯ _s/ λ; otherwise, reject it and continue the trajectory. The accepted times thus have the desired state-dependent intensity (λt)t≥0( _t)_t≥ 0. 3 Reflected Proximal Dynamics The standard proximal sampler redraws the auxiliary variable Yk+1∼(Xk,ηId)Y_k+1 (X_k,η I_d) at every iteration. This makes each update simple, but it also discards the direction of the previous move. Reflected proximal dynamics modifies this outer update while preserving the augmented target. We first present the framework with state (X,Y)(X,Y); the momentum-augmented version is introduced in Sec. 5. 3.1 The Framework Algorithm 1: Ideal reflected proximal dynamics 1: Step size η, initial states (X0,Y0)(X_0,Y_0), total steps N, refresh probability ρ, and kernels MyM_y preserving qyq_y for any y∈ℝdy ^d 2: for k=0,…,N−1k=0,…,N-1 do 3: Yk+1←2Xk−YkY_k+1← 2X_k-Y_k ⊳ reflect 4: Draw Uk∼Unif[0,1]U_k [0,1] 5: if Uk≤ρU_k≤ρ then ⊳ RGO branch 6: Draw Xk+1∼qYk+1X_k+1 q_Y_k+1 independently 7: else⊳ accelerated branch 8: Draw Xk+1∼MYk+1(Xk,⋅)X_k+1 M_Y_k+1(X_k,·) 9: end if 10: end for 11: return XNX_N Recall the conditional distribution qy(x)q_y(x) in 1. For each fixed y, let MyM_y be a Markov kernel that preserves qyq_y. Unlike an RGO draw which depends only on Yk+1Y_k+1 when sampling Xk+1X_k+1, the kernel MYk+1M_Y_k+1 may also use the previous iterate XkX_k. The alternating updates of the ideal framework are given in Sec. 3.1. The algorithm is ideal because the RGO branch draws exactly from qyq_y. The realizations in Secs. 4 and 5 replace the required exact conditional draws by controlled finite constructions, yielding implementable algorithms. We call the second branch accelerated because it is where directional persistence is introduced; invariance of an arbitrary MyM_y alone does not imply acceleration. Proposition 1 (Exactness of the framework). If MyM_y preserves qyq_y for every y, then Sec. 3.1 preserves Πη _η. Its position marginal therefore has invariant law π. Proof. Let φη:=(0,ηId) _η:=N(0,η I_d), so that Πη(x,y)=π(x)φη(y−x) _η(x,y)=π(x) _η(y-x). For every f, ∫f(x,2x−y)π(x)φη(y−x)xy=∫f(x,y′)π(x)φη(x−y′)xdy′=∫fdΠη. f(x,2x-y)\,π(x) _η(y-x)\,dx\,dy= f(x,y )\,π(x) _η(x-y )\,dx\,dy = f\,d _η. Hence the reflection preserves Πη _η. Now fix y and let M¯y(x,x′)=(1−ρ)My(x,x′)+ρqy(x′) M_y(x,x )=(1-ρ)M_y(x,x )+ρ q_y(x ) be the mixture distribution of the accelerated and RGO branches. Since MyM_y preserves qyq_y, ∫qy(x)My(x,x′)x=qy(x′)⟹∫[∫f(x′,y)M¯y(x,x′)dx′]qy(x)x=∫f(x′,y)qy(x′)dx′. q_y(x)M_y(x,x )\,dx=q_y(x ) [ f(x ,y) M_y(x,x )\,dx ]q_y(x)\,dx= f(x ,y)q_y(x )\,dx . Thus the second update preserves the conditional law qyq_y for every y. Composing it with the reflection preserves Πη _η, whose position marginal is π. ∎ 3.2 Why Reflection and Refreshment? The standard proximal sampler redraws Y|XY X independently at every iteration. This update preserves the target, but it also replaces the current displacement Y−XY-X by fresh Gaussian noise and therefore discards the direction of the previous move. To retain this information, reflected proximal dynamics instead sets Y+=2X−Y^+=2X-Y. The reflection reverses the displacement without changing its magnitude and, by Prop. 1, preserves the augmented target Πη _η. Reflection alone leaves X unchanged, so the subsequent conditional update must turn the retained displacement into motion. The accelerated branch does this by starting from the current position and applying a qY+q_Y^+-invariant kernel. Unlike an independent RGO draw, its output can remain correlated with its input; this correlation allows the direction retained by reflection to persist across iterations. Successive position increments tend to point in the same direction, so they add rather than cancel. The two concrete kernels in Secs. 4 and 5 are designed so that this persistence produces transport over the target scale. Preserving qyq_y, however, does not by itself imply contraction: the identity kernel is invariant, and a nearly deterministic kernel can retain periodic or weakly mixed components. We therefore take an independent RGO step with a small refresh probability ρ. The accelerated branch preserves direction on most iterations, while the refresh branch removes correlations that the persistent motion alone may fail to dissipate. The one-dimensional example below illustrates how this balance changes diffusive motion into approximately ballistic transport. 3.3 A One-Dimensional Picture of Acceleration Consider the quadratic potential V(x)=αx2/2V(x)=α x^2/2. The target has standard deviation 1/α1/ α, whereas one local proximal update moves on the scale η η. In the standard proximal sampler, each Gaussian redraw forgets the previous direction, so the signed increments are centered and approximately uncorrelated. As a result, the typical net displacement after n steps is only nη nη. Reaching the target scale requires roughly 1/(αη)1/(αη) iterations. Reflection keeps the previous displacement available to the accelerated branch. If this branch makes several successive moves point in a correlated direction, their signed increments add approximately linearly (nηn η) rather than cancelling as a random walk (nη nη). Achieving net displacement 1/α1/ α then takes roughly 1/αη1/ αη iterations. The occasional RGO branch prevents the chain from remaining on a periodic trajectory without erasing its direction at every step. Later, we prove distributional contraction for the two resulting chains. For both of our proposed realizations, one iteration uses O~(1+βηd) O(1+βη d) oracle queries. Multiplying this cost by O~(1/αη) O(1/ αη) iterations and choosing an optimal η=Θ~(1/(βd))η= (1/(β d)) gives O~(κd1/4) O( κd^1/4) oracle queries. Figure 1: One outer iteration in the (x,y)(x,y)-plane. Left: Standard proximal sampler (Sec. 2.2) independently redraws Yk+1∼(Xk,ηId)Y_k+1 (X_k,η I_d), discarding the previous displacement, and then draws Xk+1∼qYk+1X_k+1 q_Y_k+1 by RGO. Right: Reflected proximal sampler (8) sets Yk+1=2Xk−YkY_k+1=2X_k-Y_k, so YkY_k and Yk+1Y_k+1 are equidistant from XkX_k, and then takes the accelerated conditional transition with probability 1−ρ1-ρ or an independent RGO refresh with probability ρ. In the illustrated accelerated branch, Xk+1−Yk+1X_k+1-Y_k+1 retain the same sign as Xk−YkX_k-Y_k, preserving directional information across the iteration; the independent RGO draw need not preserve this sign. 4 Flux-Resampled Underdamped Half-Turn We first instantiate the accelerated branch of Sec. 3.1 without retaining momentum between outer iterations. At the reflected value Yk+1=yY_k+1=y, the update starts from Xk=xX_k=x, draws a Gaussian momentum, runs an event-driven underdamped half-turn, and returns its terminal position. At a high level, it freezes the nonlinear force at an independent conditional reference, follows the resulting exactly solvable affine motion, and uses random events to correct the residual force. 4.1 The Flux-Resampled Update We present the ideal flux-resampled accelerated branch, Xk+1∼FYk+1(Xk,⋅)X_k+1 F_Y_k+1(X_k,·), in Sec. 4.1. We choose a fixed numerical constant γ0∈(0,1] _0∈(0,1], and set γ=γ0αβη21+αηγ= _0 αβη^21+αη, ω=1−γ24ω= 1- γ^24, and Tγ=πωT_γ= πω. The condition βη≤1/2βη≤ 1/2 used below (Prop. 2) ensures γ<2γ<2. (Zt,Pt)(Z_t,P_t) denotes the within-update position–momentum trajectory. At an event time t, we apply an instantaneous momentum update: ZtZ_t remains continuous, while the momentum jumps from Pt−P_t- to PtP_t.11 1 By default, the process is always assumed to be càdlàg (right-continuous with left limits). All driving random sources used by one call are fresh: conditional on (x,y)(x,y), the reference draw, initial momentum, Brownian motion, driving Poisson clock, and event marks are mutually independent. Thus, following the general framework (8), the ideal flux chain makes one exact RGO draw qyq_y in every outer iteration. Algorithm 2: Ideal flux-resampled accelerated branch FyF_y 1: Xk=xX_k=x, Yk+1=yY_k+1=y, and step size η 2: Draw X~∼qy X q_y independently of x, draw P0∼γdP_0 _d, and set Z0←xZ_0← x 3: Until t=Tγt=T_γ, evolve between events according to dZt=ηPtdtdZ_t= ηP_t\,dt and dPt=−(Zt−y)/η+η∇V(X~)dt−γPtdt+2γdBtdP_t=-\(Z_t-y)/ η+ η∇ V( X)\\,dt-γ P_t\,dt+ 2γ\,dB_t 4: Set residual gradient ht←∇V(Zt)−∇V(X~)h_t←∇ V(Z_t)-∇ V( X) and generate events at rate λt←η[Pt−Tht]+ _t← η[P_t-^Th_t]_+ 5: At each event, set nt←ht/∥ht∥n_t← h_t/ h_t , draw RtR_t from the Rayleigh density r↦re−r2/21r>0r re^-r^2/21_r>0, and update Pt←Pt−(Pt−Tnt+Rt)ntP_t← P_t--(P_t-^Tn_t+R_t)n_t 6: return Xk+1←ZTγ∼Fy(x,⋅)X_k+1← Z_T_γ F_y(x,·) 4.2 Why a Half-Turn and Why Events? We now provide intuition for the design of Sec. 4.1. At fixed y, we freeze the nonlinear force at the independent reference X~ X. The conditional force then splits as −∇xlogqy(x)=x−yη+∇V(x)=x−yη+∇V(X~)+∇V(x)−∇V(X~).- _x q_y(x)= x-yη+∇ V(x)= \ x-yη+∇ V( X) \+ \∇ V(x)-∇ V( X) \. The frozen term is a harmonic force centered at xctr=y−η∇V(X~)x_ctr=y-η∇ V( X), while the second term is the residual force. Freezing makes the between-event SDE affine and exactly solvable. Drawing X~∼qy X q_y keeps the residual small on average because the stationary position and the reference follow the same localized law. First ignore the residual events. Then Zt−xctrZ_t-x_ctr follows a damped harmonic oscillator. At its half-period Tγ=π/ωT_γ=π/ω, the coefficient of the initial momentum P0P_0 vanishes, while the centered input Z0−xctrZ_0-x_ctr is multiplied by −e−γTγ/2-e^-γ T_γ/2. Thus the affine motion approximately reflects the input across xctrx_ctr, producing a coherent move rather than a local diffusive step. The choice of γ balances the damping error with the residual-force error; Sec. B.2 makes this tradeoff precise. The affine half-turn alone is not exact for qyq_y, because freezing omits the residual force. To correct it, define the residual potential Vres(z):=V(z)−V(X~)−∇V(X~)T(z−X~)V_res(z):=V(z)-V( X)-∇ V( X)^T(z- X). Along the between-event trajectory, dtVres(Zt)=ηPtTht ddtV_res(Z_t)= ηP_t^Th_t, and therefore λt=[dtVres(Zt)]+ _t=[ ddtV_res(Z_t)]_+. Thus an event can occur only while the affine trajectory climbs the omitted potential. When the process is stationary, conditional on position, the component PtTntP_t^Tn_t is standard Gaussian. The event rate weights each positive value in proportion to its size, so the incoming value Pt−Tnt>0P_t-^Tn_t>0 observed at an event has density r↦re−r2/21r>0r re^-r^2/21_r>0. The update draws RtR_t from this Rayleigh law, sets the outgoing component to PtTnt=−RtP_t^Tn_t=-R_t, and leaves the orthogonal momentum unchanged. This rule balances the probability flux of forward and momentum-reversed paths (Vanetti et al., 2017; Michel et al., 2020). Together, the event rate and Rayleigh update pair each fixed-reference path with a momentum-reversed path of equal stationary probability. Averaging over the independent reference and discarding momentum therefore yields a reversible position kernel FyF_y for qyq_y, as formalized next. Proposition 2 (Conditional invariance and reversibility). Suppose βη≤1/2βη≤ 1/2. Conditional on the reflected auxiliary value Yk+1=yY_k+1=y, if the input position satisfies Xk|Yk+1=y∼qyX_k Y_k+1=y q_y, then the output satisfies Xk+1|Yk+1=y∼qyX_k+1 Y_k+1=y q_y as well. Moreover, FyF_y satisfies detailed balance with respect to qyq_y, i.e., qy(x)Fy(x,x′)=qy(x′)Fy(x′,x)q_y(x)F_y(x,x )=q_y(x )F_y(x ,x). Consequently, Sec. 3.1 with accelerated branch FyF_y (Sec. 4.1) preserves Πη _η, whose position marginal is π. The proof in Sec. B.1 establishes non-explosion and sums the preceding Rayleigh flux balance over all possible event times, yielding the complete path-reversal argument. 4.3 From the Ideal Kernel to an Implementable Update Sec. 4.1 is ideal in two places: it draws X~ X exactly from qyq_y, and it uses the state-dependent event rate λt _t. The implementable update replaces these operations by a certified approximate RGO and capped Poisson thinning, with a finite work limit for every auxiliary optimization. Certified approximate RGO. Gradient descent on the conditional objective z↦V(z)+∥z−y∥2/(2η)z V(z)+ z-y ^2/(2η) starts from the incoming position and stops when it finds a point x x satisfying ∥y−η∇V(x^)−x^∥≤dη y-η∇ V( x)- x ≤ dη. This is an observable approximate-proximal certificate, not yet a conditional draw; it supplies the center required by the first-order rejection sampler that produces an approximate sample from qyq_y. The solver is stopped after a deterministic query budget; if it has not certified the residual, the algorithm terminates and returns the current position. Otherwise, the approximate draw is used either as the RGO-branch output or as the accelerated branch’s reference X~ X. Capped event clock. Once ∇V(X~)∇ V( X) is stored, the affine diffusion between candidate events can be sampled exactly without further oracle queries. Candidate times are proposed at rate λ¯ λ; at a candidate time t, one query evaluates λt=η[PtT∇V(Zt)−∇V(X~)]+ _t= η[P_t^T\∇ V(Z_t)-∇ V( X)\]_+, and the proposal is accepted with probability min1,λt/λ¯ \1, _t/ λ\. This exactly simulates the truncated rate minλt,λ¯ \ _t, λ\, and the resulting bias is the probability that the truncated and ideal clocks disagree (Sec. 2.3). Coupling to the ideal chain. The implementable and ideal iterations use common randomness until one of three disagreements occurs: the center solver reaches its cap, the approximate and exact conditional draws differ, or the ideal event rate exceeds λ¯ λ. The appendix assigns each source a per-iteration error proportional to the target accuracy divided by the total number of iterations, couples the two chains until their first disagreement, and applies a union bound over the run. On the agreement event, the implementable chain follows the ideal chain exactly; Sec. B.4 gives the caps, tolerances, and complete query-complexity analysis. 4.4 Convergence and Oracle Complexity We now establish the convergence of the ideal and implementable chains. Let Kη,ρK_η,ρ be the transition kernel of one complete reflected iteration with the ideal accelerated branch FyF_y. In other words, if (X0,Y0)∼ν(X_0,Y_0) ν, then (Xm,Ym)∼νKη,ρm(X_m,Y_m) ν K_η,ρ^m. Let ρ⋆∈(0,1) _ ∈(0,1) be the fixed universal numerical constant furnished by the theorem below, and set ρ=ρ⋆αηlog(eαη)≤ρ⋆<1.ρ= _ αη ( eαη )≤ _ <1. (2) Theorem 1 (Ideal flux contraction). There exist universal constants c0,c,C>0c_0,c,C>0 and a fixed ρ⋆∈(0,1) _ ∈(0,1), independent of all problem parameters, such that, if βη≤c0βη≤ c_0, then every integer nblk≥Cαηlogeαηn_blk≥ C αη eαη satisfies, for every ν≪Πην _η and every nonnegative integer J, χ2(νKη,ρJnblk∥Πη)≤e−Jχ2(ν∥Πη).χ^2(ν K_η,ρ^Jn_blk\| _η) ^-Jχ^2(ν\| _η). (3) The proof first establishes in Sec. B.2 that one conditional half-turn nearly reverses centered affine position coordinates, with an error weighted by the same conditional covariance deficit through which reflection transfers information. At the outer-chain level, reflection moves variation in Y into the conditional position, the half-turn reverses that component, and the RGO branch damps the remainder. A modified L2L^2 energy then contracts over blocks of nblk=O~(1/αη)n_blk= O(1/ αη) iterations; the complete argument is in Sec. B.3. Theorem 2 (Implementable flux half-turn). Under Assump. 1, let 0<ε<1/40< <1/4, Δ≥1 ≥ 1, and let ν0 _0 be an initial position law satisfying D2(ν0∥π)≤ΔD_2( _0\|π)≤ . There exist universal numerical constants c,C>0c,C>0, independent of d,κ,Δ,εd,κ, , , such that the following choice is valid: L=Δ+logCdκε,η=cβ(dL+L).L= + Cdκ , η= cβ( dL+L). (4) Initialize X0∼ν0X_0 _0 and Y0|X0∼(X0,ηId)Y_0 X_0 (X_0,η I_d). Then run Sec. 3.1 with the accelerated branch in Sec. 4.1, using the approximate RGO and capped Poisson clock specified in Sec. B.4. Its output distribution ν ν satisfies TV(ν^,π)≤εTV( ν,π)≤ . Writing NqryN_qry for its total number of oracle queries, Nqry≤Cκ(d+L)1/4L3/4(Δ+log4ε)log(eκL).EN_qry≤ C κ\,(d+L)^1/4L^3/4 ( + 4 ) (eκ L). (5) In particular, a constant-D2D_2 warm start requires O~(κd1/4) O( κd^1/4) expected oracle queries. The step size in 4 simultaneously controls the approximate RGO and the capped event clock. Summing the center-solver, conditional-draw, and event-proposal costs over the ideal mixing horizon gives 5; the complete error allocation and query-complexity analysis are in Sec. B.4. 5 Persistent-Momentum Quarter-Turn The flux method uses momentum only inside one accelerated update and then discards it. The persistent construction instead augments the outer state to (X,Y,P)(X,Y,P), carries the terminal momentum into the next iteration, and requires its accelerated branch to preserve qy⊗γdq_y _d. As in the flux method, we first analyze this branch with an exact conditional reference, then replace that reference by an implementable Gaussian construction that reserves approximate RGO calls for the refresh branch. 5.1 The Persistent Update At fixed Y=yY=y, let ζy _y be a reference law depending only on y. For an incoming position–momentum pair (x,p)(x,p), Sec. 5.1 defines the phase-space kernel Gyζ((x,p),⋅)G_y^ζ((x,p),·). The ideal contraction analysis takes ζy=qy _y=q_y; the implementable reference is constructed in Sec. 5.3. As in Sec. 4.1, an event is an instantaneous momentum update Pt−↦PtP_t- P_t. Conditional on (x,p,y)(x,p,y), the reference draw, Poisson clock, and event marks are fresh and mutually independent. Algorithm 3: Ideal persistent-momentum accelerated branch GyζG_y^ζ 1: Incoming pair (x,p)(x,p), fixed auxiliary value y, step size η, and reference law ζy _y 2: Draw X~∼ζy X _y independently of (x,p)(x,p), and set (Z0,P0)←(x,p)(Z_0,P_0)←(x,p) 3: Until t=π/2t=π/2, evolve between events according to Z˙t=ηPt Z_t= ηP_t and P˙t=−(Zt−y)/η−η∇V(X~) P_t=-(Z_t-y)/ η- η∇ V( X) 4: Set residual gradient ht←∇V(Zt)−∇V(X~)h_t←∇ V(Z_t)-∇ V( X) and generate events at rate λt←η[Pt−Tht]+ _t← η[P_t-^Th_t]_+ 5: At each event, set nt←ht/∥ht∥n_t← h_t/ h_t and update Pt←Pt−2(Pt−Tnt)ntP_t← P_t--2(P_t-^Tn_t)n_t 6: return (X+,P+)←(Zπ/2,Pπ/2)∼Gyζ((x,p),⋅)(X^+,P^+)←(Z_π/2,P_π/2) G_y^ζ((x,p),·) A complete ideal persistent iteration first reflects Yk+1=2Xk−YkY_k+1=2X_k-Y_k. It then redraws both position and momentum on the RGO branch, or applies the persistent phase-space kernel on the accelerated branch: (Xk+1,Pk+1)∼qYk+1⊗γd,with probability ρ,GYk+1ζ((Xk,Pk),⋅),otherwise.(X_k+1,P_k+1) casesq_Y_k+1 _d,&with probability ρ,\\ G_Y_k+1^ζ((X_k,P_k),·),&otherwise. cases 5.2 Why a Quarter-Turn? The persistent method must transfer position information into momentum rather than discard it. Fix the reference X~ X, set xctr=y−η∇V(X~)x_ctr=y-η∇ V( X), and let Wt=(Zt−xctr)/ηW_t=(Z_t-x_ctr)/ η. Without events, W˙t=Pt W_t=P_t and P˙t=−Wt P_t=-W_t, so the motion is a rotation. After one quarter-period, (W0,P0)↦(P0,−W0)(W_0,P_0) (P_0,-W_0), so centered position information is stored in the terminal momentum and can return to position in later iterations. A half-turn instead maps (W0,P0)(W_0,P_0) to (−W0,−P0)(-W_0,-P_0), producing a sign flip without exchanging the two coordinates. Events correct the residual force through specular momentum reflection, which preserves momentum magnitude and is therefore compatible with persistence. The RGO branch supplies the occasional refreshment described in Sec. 3. Proposition 3 (Conditional invariance). For every fixed y and reference value x~ x, the accelerated-branch motion preserves qy⊗γdq_y _d. Consequently, GyζG_y^ζ preserves qy⊗γdq_y _d whenever X~∼ζy X _y is independent of the incoming position and momentum. The complete persistent iteration above therefore preserves Π^η:=Πη⊗γd _η:= _η _d. The proof in Sec. C.1 combines harmonic energy conservation with the residual event rate and specular reflection. It establishes non-explosion and pairs each stationary path with its momentum-reversed path run backward, which yields the conditional invariance stated above. 5.3 From the Ideal Kernel to an Implementable Update The ideal contraction analysis draws the accelerated-branch reference exactly from qyq_y, but the implementable update avoids this per-iteration RGO. Starting gradient descent at y, the conditional mode solve in Sec. 4.3 returns a point x^y x_y satisfying ∥y−η∇V(x^y)−x^y∥≤dη y-η∇ V( x_y)- x_y ≤ dη. The accelerated branch then draws G∼γdG _d and uses the shifted Gaussian reference X~=y−η∇V(x^y)+ηG X=y-η∇ V( x_y)+ ηG. Because this reference depends only on y and fresh randomness, Prop. 3 still gives exact conditional invariance; its approximation to qyq_y is needed only for the contraction estimate. The mode solve and event clock use the same capping principle as the flux method, while an approximate RGO is invoked only on a refresh iteration, which occurs with probability ρ. The appendix couples this implementable chain to an uncapped shifted-reference comparison chain until a mode cap, clock cap, or approximate refresh causes disagreement, and gives the complete error and query-complexity analysis in Sec. C.4. 5.4 Convergence and Oracle Complexity We now prove the convergence of the ideal and implementable chains. Let KperK_per be the transition kernel of one complete ideal persistent iteration with ζy=qy _y=q_y. In other words, if (X0,Y0,P0)∼ν(X_0,Y_0,P_0) ν, then (Xm,Ym,Pm)∼νKperm(X_m,Y_m,P_m) ν K_per^m. Set g=αηg= αη and ℓg=1+log(1/g) _g=1+ (1/g). Theorem 3 (Ideal persistent contraction). There exist universal constants c0,C>0c_0,C>0 and a fixed cρ>0c_ρ>0, independent of all problem parameters, such that the following holds. If 0<η≤c0/β0<η≤ c_0/β, set ρ=cρgℓgρ=c_ρg _g and suppose ρ≤1/4ρ≤ 1/4. Then every integer nblk≥Cℓg/gn_blk≥ C _g/g satisfies, for every ν≪Π^η=Πη⊗γdν _η= _η _d and every nonnegative integer J, χ2(νKperJnblk∥Π^η)≤e−Jχ2(ν∥Π^η).χ^2(ν K_per^Jn_blk\| _η) ^-Jχ^2(ν\| _η). (6) Consequently, initialize X0∼ν0X_0 _0, Y0|X0∼(X0,ηId)Y_0 X_0 (X_0,η I_d), and independently P0∼γdP_0 _d. Suppose D2(ν0∥π)≤ΔD_2( _0\|π)≤ . For ε≤1/2 ≤ 1/2, one obtains TV(Law(Xn),π)≤εTV(Law(X_n),π)≤ if n≥Cℓg1/2g(Δ+log4ε).n≥ C _g^1/2g ( + 4 ). (7) Since g=αηg= αη, one constant-factor reduction of χ2χ^2 divergence takes O~(1/αη) O(1/ αη) iterations, and repeating this contraction from the warm start gives 7. The proof tracks how an observable moves among the auxiliary variable Y, conditional position, and momentum: reflection creates conditional position dependence, the quarter-turn rotates it into momentum and back, and the RGO branch removes the remaining conditional dependence. The appendix first proves the covariance-weighted approximate quarter-rotation in Lem. 13, then converts this three-coordinate transfer mechanism into the L2L^2 operator estimate in Prop. 11; applying that estimate to density ratios yields 6. Theorem 4 (Implementable persistent sampler). Under Assump. 1, let ν0 _0 satisfy D2(ν0∥π)≤ΔD_2( _0\|π)≤ , and let ε≤1/2 ≤ 1/2. There exist fixed universal constants c,C>0c,C>0 so that the following holds. Define Λε=Δ+log(Cdκ/ε) _ = + (Cdκ/ ) and η=c/β(dΛε+Λε)η=c/\β( d _ + _ )\. Initialize X0∼ν0X_0 _0, Y0|X0∼(X0,ηId)Y_0 X_0 (X_0,η I_d), and independently P0∼γdP_0 _d. Run the implementable algorithm specified in Sec. C.4; its output law ν ν satisfies TV(ν^,π)≤εTV( ν,π)≤ . If NqryN_qry counts all oracle queries and NrgoN_rgo counts approximate-RGO calls, then Nqry _qry ≤O~(κ(Δ+log(8/ε))(dΛε+Λε)1/2), ≤ O( κ\,( + (8/ ))( d _ + _ )^1/2), (8) Nrgo _rgo ≤O~(Δ+log(8/ε)). ≤ O( + (8/ )). In particular, a constant-D2D_2 warm start gives O~(κd1/4) O( κd^1/4) expected oracle queries and only polylog(d,κ,1/ε)polylog(d,κ,1/ ) expected RGO calls. The ideal horizon is O~(g−1Δ+log(1/ε)) O(g^-1\ + (1/ )\), while refreshes occur with probability ρ=O~(g)ρ= O(g); their product gives the polylogarithmic RGO count. Every iteration pays for the capped mode solve; accelerated iterations then use the shifted Gaussian reference and capped event clock, whereas refresh iterations invoke the approximate RGO. At the step size in the theorem these costs are O~(1) O(1) per iteration, and summing them over the ideal horizon gives the stated total query bound. The full error and query-complexity analysis is in Sec. C.4; an optional quadratic-baseline refinement is stated and proved in Cor. 3. 6 Removing the Supplied Warm Start The rates in Secs. 4 and 5 assume a constant-D2D_2 warm start. Our cold-start pipeline has three stages: localization constructs a law with finite D2D_2 divergence, a coarse proximal chain reduces this divergence to a constant, and either reflected sampler performs the final high-accuracy sampling. When a finite bound D2(ν0∥π)≤BD_2( _0\|π)≤ B is available, we skip localization and choose between direct reflected sampling and the coarse warm-up route. The resulting bounds are end-to-end: all stages used are included in the query count. Warm-start removal preserves the high-accuracy character of the reflected final stage: the overall dependence on 1/ε1/ remains polylogarithmic. Write the coarse step size as u=r/βu=r/β, where r∈[1,κ]r∈[1,κ] interpolates between the step sizes 1/β1/β and 1/α1/α, and let the ideal coarse proximal chain be the Gibbs sampler in 1 with step size u. The Rényi contraction of Chen et al. (2022, Theorem 3) makes the ideal coarse chain constant-D2D_2-warm after O((κ/r)log(eB))O((κ/r) (eB)) iterations. The RGO has condition number (β+1/u)/(α+1/u)≤1+r(β+ 1u)/(α+ 1u)≤ 1+r. Increasing r therefore shortens the ideal coarse chain but makes each conditional draw harder. The algorithm replaces those exact draws by the cheaper of two implementable posterior samplers: a regular interface based on smoothed Picard HMC (Chen et al., 2026c) followed by reflected finalization, or a harmonic-base interface based on exact diffusion (Chen et al., 2026d; Chen et al., 2026a). For fixed r, multiplying the number of coarse iterations by the cost of the regular conditional interface gives the cost of the regular-interface route: O~(κr−1+r1/6d1/6+r−1/2d1/4log(eB)). O(κ\r^-1+r^1/6d^1/6+r^-1/2d^1/4\ (eB)). Optimizing over r∈[1,κ]r∈[1,κ] and considering the relative sizes of κ and d yields the best polynomial cost factor below, which we denote by Ψ(κ,d) (κ,d): Ψ(κ,d)≍κ2/3d1/3,when d≲κ16/7, with r≍κ,κd3/16,when κ16/7≲d≲κ8, with r≍d1/8,κd1/4,when d≳κ8, with r≍κ, (κ,d) casesκ^2/3d^1/3,&when $d κ^16/7$, with $r κ$,\\ κ d^3/16,&when $κ^16/7 d κ^8$, with $r d^1/8$,\\ κd^1/4,&when $d κ^8$, with $r κ$, cases Thus the optimized coarse route costs O~(Ψ(κ,d)log(eB)) O( (κ,d) (eB)). This bound includes the implementable warm-up and the final reflected sampler, whose lower-order cost is absorbed in every regime. It does not include the one-time localization cost needed when the input is only a point. Theorem 5 (Warm-start removal). (i) Finite-divergence input. Suppose a sample from an initial law ν0 _0 is available and the known number B≥1B≥ 1 satisfies D2(ν0∥π)≤BD_2( _0\|π)≤ B. For every 0<ε<1/40< <1/4, the wrapper using either reflected sampler as its final stage returns an output law ν ν with TV(ν^,π)≤εTV( ν,π)≤ . If NqryN_qry is the total number of oracle queries, then Nqry≤O~(minκB(dB+B)1/2,Ψ(κ,d)log(eB)).EN_qry≤ O( \ κB( dB+B)^1/2, (κ,d) (eB)\). (9) If B is polynomially bounded in d,κ,1/εd,κ,1/ , this is at most O~(κ6/7d1/4) O(κ^6/7d^1/4). (i) Supplied point. Given x0∈ℝdx_0 ^d and 0<ε<1/40< <1/4, let x∗x_* be the minimizer of V, define log+(t)=max1,logt _+(t)= \1, t\, and set Qloc=κlog+[κ2+β∥x0−x∗∥]Q_loc= κ _+[κ\2+ β x_0-x_* \]. An implementable algorithm returns an output law ν ν with TV(ν^,π)≤εTV( ν,π)≤ and Nqry≤O~(Qloc+Ψ(κ,d))≤O~(Qloc+κ6/7d1/4).EN_qry≤ O(Q_loc+ (κ,d))≤ O(Q_loc+κ^6/7d^1/4). (10) The first term in 9 skips the coarse stage and applies a reflected sampler directly with the full divergence budget B. The second term uses the coarse route just described: warm the input using step size u=r/βu=r/β, then apply a constant-warm reflected finalizer. The ideal coarse chain is used only for analysis and is not executed. Its endpoint is constant-D2D_2-warm, whereas the implementable warm-up replaces each exact conditional draw by an implementable posterior sampler. The two endpoints are close in total variation. Applying the same capped reflected finalizer to both cannot increase this discrepancy, so the convergence guarantee for the ideal coarse endpoint transfers to the algorithmic output. For a supplied point, accelerated optimization stops at the observable certificate ∥∇V(xref)∥≤α/β ∇ V(x_ref) ≤α/ β. One draw from (xref,(2β)−1Id)N(x_ref,(2β)^-1I_d) then has D2≤1+d2log(2κ)D_2≤ 1+ d2 (2κ), so the finite-divergence result applies. The algorithm does not use x∗x_*; it appears only in the localization bound. The regime optimization, error handoff, and localization proof are given in App. D. 7 Conclusion Reflected proximal dynamics accelerate the outer chain by combining exact Gaussian reflection, conditional kinetic transport, and controlled refreshment. From a constant-χ2χ^2 warm start, both realizations use O~(κd1/4) O( κd^1/4) expected oracle queries, while the persistent construction makes only a polylogarithmic expected number of approximate-RGO calls. Including localization and warm-up yields an end-to-end expected oracle complexity of O~(κ6/7d1/4) O(κ^6/7d^1/4). Both the warm-start and end-to-end guarantees have only polylogarithmic dependence on 1/ε1/ , arising from an ideal mixing horizon of order log(1/ε) (1/ ). Extending the analysis to more general settings, such as weakening strong log-concavity of the target to functional inequalities (e.g., log-Sobolev and Poincaré inequalities), is a natural direction beyond the current setting explored in this paper, and we leave this as future work. AI use statement The idea of accelerating the proximal sampler using underdamped dynamics originated with the authors. We used generative AI tools to brainstorm implementation details and algorithmic variants, assist in drafting and revising mathematical proofs. AI-generated suggestions were treated as unverified drafts rather than as authoritative sources. The authors reviewed and revised all AI-assisted material, checked the final theorem statements and proofs, and verified the citations used in the paper against the referenced works. We take responsibility for the final content of this work, including all text, claims, and artifacts roduced with the aid of generative AI. References Altschuler and Chewi (2024) J. M. Altschuler and S. Chewi Faster high-accuracy log-concave sampling via algorithmic warm starts. J. ACM 71 (3). External Links: ISSN 0004-5411, Link, Document Cited by: Table 1. Betancourt and Girolami (2019) M. Betancourt and M. Girolami Hamiltonian Monte Carlo for hierarchical models. In Current trends in Bayesian methodology with applications, p. 119–142. Cited by: §1. Bierkens et al. (2020) J. Bierkens, S. Grazzi, K. Kamatani, and G. Roberts The boomerang sampler. In Proceedings of the 37th International Conference on Machine Learning, H. D. I and A. Singh (Eds.), Proceedings of Machine Learning Research, Vol. 119, p. 908–918. External Links: Link Cited by: Appendix A, §B.1, §C.1. Bouchard-Côté et al. (2018) A. Bouchard-Côté, S. J. Vollmer, and A. Doucet The bouncy particle sampler: a nonreversible rejection-free Markov chain Monte Carlo method. Journal of the American Statistical Association 113 (522), p. 855–867. External Links: Document Cited by: Appendix A, §C.1. Chen et al. (2026a) F. Chen, S. Chewi, C. C. Daskalakis, and A. Rakhlin High-accuracy sampling for diffusion models and log-concave distributions. In Proceedings of the 43rd International Conference on Machine Learning, T. Zhang, M. Dudik, M. Jaggi, A. Agarwal, S. Li, D. Schuurmans, J. Zhu, F. Berkenkamp, H. Dong, and A. Bietti (Eds.), Proceedings of Machine Learning Research, Vol. 306, p. 13835–13869. External Links: Link Cited by: Appendix A, Table 1, §D.1, §6, Proposition 12, Proposition 9. Chen et al. (2026b) F. Chen, S. Chewi, J. Lu, and M. S. Zhang Accelerated high-accuracy sampling from a warm start via the proximal bouncy particle sampler. arXiv preprint arXiv:2609.06905. Cited by: Appendix A, Table 1, §1. Chen et al. (2026c) F. Chen, S. Chewi, J. Lu, and M. S. Zhang Smoothed Picard Hamiltonian Monte Carlo. arXiv preprint arXiv:2609.06906. Cited by: Table 1, §D.1, §D.1, §6. Chen et al. (2026d) F. Chen, S. Chewi, A. Rakhlin, and M. S. Zhang Exact simulation of diffusions and improved algorithms for log-concave sampling. arXiv preprint arXiv:2608.05022. Cited by: Table 1, §D.1, §2.1, §6. Chen et al. (2022) Y. Chen, S. Chewi, A. Salim, and A. Wibisono Improved analysis for a proximal algorithm for sampling. In Proceedings of Thirty Fifth Conference on Learning Theory, P. Loh and M. Raginsky (Eds.), Proceedings of Machine Learning Research, Vol. 178, p. 2984–3014. External Links: Link Cited by: Appendix A, Table 1, §1, §2.2, §6. Cheng et al. (2018) X. Cheng, N. S. Chatterji, P. L. Bartlett, and M. I. Jordan Underdamped Langevin MCMC: a non-asymptotic analysis. In Proceedings of the 31st Conference On Learning Theory, S. Bubeck, V. Perchet, and P. Rigollet (Eds.), Proceedings of Machine Learning Research, Vol. 75, p. 300–323. External Links: Link Cited by: Appendix A, §1. Chewi et al. (2022) S. Chewi, M. A. Erdogdu, M. Li, R. Shen, and S. Zhang Analysis of Langevin Monte Carlo from Poincaré to log-Sobolev. In Proceedings of Thirty Fifth Conference on Learning Theory, P. Loh and M. Raginsky (Eds.), Proceedings of Machine Learning Research, Vol. 178, p. 1–2. External Links: Link Cited by: §1. Chewi (2023) S. Chewi Log-concave sampling. Book draft available at https://chewisinho.github.io 9, p. 17–18. Cited by: §B.1, §C.2.1, §C.2.1, §C.4.2, §1, §2.1. Dalalyan and Riou-Durand (2020) A. S. Dalalyan and L. Riou-Durand On sampling from a log-concave density using kinetic Langevin diffusions. Bernoulli 26 (3), p. 1956 – 1988. External Links: Document, Link Cited by: Appendix A, §1. Dalalyan (2017) A. S. Dalalyan Theoretical guarantees for approximate sampling from smooth and log-concave densities. Journal of the Royal Statistical Society Series B: Statistical Methodology 79 (3), p. 651–676. External Links: ISSN 1369-7412, Document, Link, https://academic.oup.com/jrsssb/article-pdf/79/3/651/49214807/jrsssb_79_3_651.pdf Cited by: §1. Durmus et al. (2019) A. Durmus, S. Majewski, and B. Miasojedow Analysis of Langevin Monte Carlo via convex optimization. Journal of Machine Learning Research 20 (73), p. 1–46. External Links: Link Cited by: §1. Durmus et al. (2018) A. Durmus, É. Moulines, and M. Pereyra Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau. SIAM Journal on Imaging Sciences 11 (1), p. 473–506. External Links: Document, Link, https://doi.org/10.1137/16M1108340 Cited by: Appendix A. Fan et al. (2023) J. Fan, B. Yuan, and Y. Chen Improved dimension dependence of a proximal algorithm for sampling. 195, p. 1473–1521. External Links: Link Cited by: Appendix A, Table 1. Gelman et al. (1995) A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin Bayesian data analysis. chapman and hall/cRc. Cited by: §1. Ghosal and van der Vaart (2017) S. Ghosal and A. van der Vaart Fundamentals of nonparametric bayesian inference. Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press. Cited by: §1. Gopi et al. (2023) S. Gopi, Y. T. Lee, D. Liu, R. Shen, and K. Tian Algorithmic aspects of the log-Laplace transform and a non-Euclidean proximal sampler. In Proceedings of Thirty Sixth Conference on Learning Theory, G. Neu and L. Rosasco (Eds.), Proceedings of Machine Learning Research, Vol. 195, p. 2399–2439. External Links: Link Cited by: Appendix A. Gouraud et al. (2025) N. Gouraud, L. Journel, and P. Monmarché The velocity jump Langevin process and its splitting scheme: long time convergence and numerical accuracy. The SMAI Journal of Computational Mathematics 11, p. 533–584. External Links: Document Cited by: Appendix A. He et al. (2024) Y. He, A. Mousavi-Hosseini, K. Balasubramanian, and M. A. Erdogdu A separation in heavy-tailed sampling: Gaussian vs. stable oracles for proximal samplers. In Advances in Neural Information Processing Systems, A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Paquet, J. Tomczak, and C. Zhang (Eds.), Vol. 37, p. 65260–65296. External Links: Document, Link Cited by: Appendix A. Kook et al. (2024) Y. Kook, S. S. Vempala, and M. S. Zhang In-and-out: algorithmic diffusion for sampling convex bodies. In Advances in Neural Information Processing Systems, A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Paquet, J. Tomczak, and C. Zhang (Eds.), Vol. 37, p. 108354–108388. External Links: Document, Link Cited by: Appendix A. Kook and Vempala (2025) Y. Kook and S. S. Vempala Zeroth-order logconcave sampling. arXiv preprint arXiv:2507.18021. Cited by: §D.2. Lee et al. (2021) Y. T. Lee, R. Shen, and K. Tian Structured logconcave sampling with a restricted Gaussian oracle. In Proceedings of Thirty Fourth Conference on Learning Theory, M. Belkin and S. Kpotufe (Eds.), Proceedings of Machine Learning Research, Vol. 134, p. 2993–3050. External Links: Link Cited by: Appendix A, §1, §2.2. Liang and Chen (2022) J. Liang and Y. Chen A proximal algorithm for sampling from non-smooth potentials. In 2022 Winter Simulation Conference (WSC), p. 3229–3240. Cited by: Appendix A. Liang and Chen (2023) J. Liang and Y. Chen A proximal algorithm for sampling. Transactions on Machine Learning Research. Note: External Links: ISSN 2835-8856, Link Cited by: Appendix A. Liang and Chen (2026) J. Liang and Y. Chen Proximal oracles for optimization and sampling. Journal of Optimization Theory and Applications 208 (110). External Links: Document Cited by: Appendix A. Ma et al. (2021) Y. Ma, N. S. Chatterji, X. Cheng, N. Flammarion, P. L. Bartlett, and M. I. Jordan Is there an analog of Nesterov acceleration for gradient-based MCMC?. Bernoulli 27 (3), p. 1942 – 1992. External Links: Document, Link Cited by: Appendix A, §1. Michel et al. (2020) M. Michel, A. Durmus, and S. Sénécal Forward event-chain Monte Carlo: fast sampling by randomness control in irreversible Markov chains. Journal of Computational and Graphical Statistics 29 (4), p. 689–702. External Links: Document Cited by: Appendix A, §B.1, §4.2. Mitra and Wibisono (2025) S. Mitra and A. Wibisono Fast convergence of Φ -divergence along the unadjusted Langevin algorithm and proximal sampler. In Proceedings of The 36th International Conference on Algorithmic Learning Theory, G. Kamath and P. Loh (Eds.), Proceedings of Machine Learning Research, Vol. 272, p. 846–869. External Links: Link Cited by: Table 1, §1, §2.2. Monmarché et al. (2023) P. Monmarché, M. Rousset, and P. Zitt Exact targeting of Gibbs distributions using velocity-jump processes. Stochastics and Partial Differential Equations: Analysis and Computations 11 (3), p. 908–947. External Links: Document Cited by: Appendix A, §B.1, §C.1. Newman and Barkema (1999) M. E. J. Newman and G. T. Barkema Monte Carlo methods in statistical physics. Oxford University Press. External Links: ISBN 9780198517962, Document, Link Cited by: §1. Salim and Richtarik (2020) A. Salim and P. Richtarik Primal dual interpretation of the proximal stochastic gradient Langevin algorithm. In Advances in Neural Information Processing Systems, H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (Eds.), Vol. 33, p. 3786–3796. External Links: Link Cited by: Appendix A. Shen and Lee (2019) R. Shen and Y. T. Lee The randomized midpoint method for log-concave sampling. In Advances in Neural Information Processing Systems, H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Eds.), Vol. 32, p. . External Links: Link Cited by: Appendix A, §1. Vanetti et al. (2017) P. Vanetti, A. Bouchard-Côté, G. Deligiannidis, and A. Doucet Piecewise-deterministic Markov chain Monte Carlo. arXiv preprint arXiv:1707.05296. Cited by: Appendix A, §B.1, §C.1, §4.2. Vempala and Wibisono (2019) S. Vempala and A. Wibisono Rapid convergence of the unadjusted Langevin algorithm: isoperimetry suffices. In Advances in Neural Information Processing Systems, H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Eds.), Vol. 32, p. . External Links: Link Cited by: §1, §2.1. Wibisono (2026) A. Wibisono Mixing time of the proximal sampler in relative Fisher information via strong data processing inequality. IEEE Transactions on Information Theory 72 (5), p. 2667–2688. External Links: Document Cited by: Appendix A. Yuan et al. (2025a) B. Yuan, J. Fan, J. Liang, and Y. Chen Client-only distributed Markov chain Monte Carlo sampling over a network. Transactions on Machine Learning Research. Note: External Links: ISSN 2835-8856, Link Cited by: Appendix A. Yuan et al. (2025b) B. Yuan, J. Fan, J. Liang, and Y. Chen Proximal sampler with adaptive step size. In Proceedings of The 28th International Conference on Artificial Intelligence and Statistics, Y. Li, S. Mandt, S. Agrawal, and E. Khan (Eds.), Proceedings of Machine Learning Research, Vol. 258, p. 1387–1395. External Links: Link Cited by: Appendix A. Yuan et al. (2023) B. Yuan, J. Fan, J. Liang, A. Wibisono, and Y. Chen On a class of Gibbs sampling over networks. In Proceedings of the 36th Conference on Learning Theory, Proceedings of Machine Learning Research, Vol. 195, p. 5754–5780. External Links: Link Cited by: Appendix A. Zhang et al. (2023) S. Zhang, S. Chewi, M. Li, K. Balasubramanian, and M. A. Erdogdu Improved discretization analysis for underdamped Langevin Monte Carlo. In Proceedings of Thirty Sixth Conference on Learning Theory, G. Neu and L. Rosasco (Eds.), Proceedings of Machine Learning Research, Vol. 195, p. 36–71. External Links: Link Cited by: Appendix A, §1. Table of Contents Appendix A Related Work and Detailed Comparison Proximal and kinetic samplers. Earlier proximal Langevin methods use Moreau regularization or proximal mappings to discretize overdamped Langevin dynamics (Durmus et al., 2018; Salim and Richtarik, 2020). The Gaussian proximal sampler instead separates global contraction from sampling a localized conditional (Lee et al., 2021; Chen et al., 2022). Approximate-RGO algorithms improve the local implementation (Fan et al., 2023; Chen et al., 2026a); our methods instead modify the outer dynamics to retain information across conditional updates. Kinetic Langevin samplers obtain acceleration by retaining velocity, with specialized discretizations needed to control the resulting finite-step error (Cheng et al., 2018; Dalalyan and Riou-Durand, 2020; Ma et al., 2021; Shen and Lee, 2019; Zhang et al., 2023). We use the same principle inside each localized conditional problem, while an exact Gaussian reflection transports memory between successive conditionals. Event-driven kinetic updates. The bouncy particle sampler corrects deterministic motion through velocity events triggered by the positive directional derivative of the potential (Bouchard-Côté et al., 2018). Event-driven Hamiltonian samplers and the Boomerang sampler instead integrate an exactly solvable reference motion and use such events to correct the residual force (Vanetti et al., 2017; Bierkens et al., 2020). Gaussian flux-resampling of the force-aligned velocity component appears in the framework of Vanetti et al. (2017) and Forward Event-Chain Monte Carlo (Michel et al., 2020); related samplers combine Langevin motion with randomized velocity events (Monmarché et al., 2023; Gouraud et al., 2025). Our flux half-turn places these ingredients inside a reflected proximal iteration and selects the duration and damping to obtain the weighted reversal estimate used by the outer-chain analysis. End-to-end complexity. Table 1 compares representative high-accuracy results under smooth strong log-concavity. We report the error metric in which each paper states its guarantee and the initialization that its bound requires. The exact-ULD/FORS row is not a proximal method; we include it as the closest global benchmark with a different tradeoff between d and κ. Table 1: Representative high-accuracy bounds for smooth, strongly log-concave sampling. The table reports oracle-query counts in expectation unless marked otherwise. The notation O~ O suppresses polylogarithmic factors in d,κ,1/εd,κ,1/ . Paper Method Initialization Metric Complexity Chen et al. (2022) Gaussian proximal sampler with rejection RGO polynomial KL/Rényi divergence KL, χ2χ^2, Rényi O~(κd) O(κ d) Fan et al. (2023) approximate-rejection proximal sampler Gaussian cold start TV, W2W_2, χ2χ^2 O~(κd1/2) O(κ d^1/2) Altschuler and Chewi (2024) ULMC warm-up with MALA/proximal reduction known or located minimizer; no warm law TV, KL, χ2χ^2, W2W_2 O~(κd1/2) O(κ d^1/2) w.h.p. Mitra and Wibisono (2025) proximal sampler with a basic rejection compiler polynomial Φ -divergence Φ -divergence O~(κd) O(κ d) Chen et al. (2026a) proximal sampler with a FORS conditional sampler finite-Rényi Gaussian start fixed-order Rényi O~(κd1/2) O(κ d^1/2) Chen et al. (2026d) exact ULD/FORS Gaussian with D2q=O(dlogκ)D_2q=O(d κ), fixed q fixed-order Rényi O~(κ2/3d1/3) O(κ^2/3d^1/3) Chen et al. (2026c); Chen et al. (2026b) SPHMC warm-up with a PBPS finalizer no supplied warm distribution TV O~(κ7/6d1/6+κd1/4) O(κ^7/6d^1/6+ κd^1/4) This work reflected proximal dynamics with optimized initialization no supplied warm distribution TV O~(κ6/7d1/4) O(κ^6/7d^1/4) Extensions of proximal sampling. The two-variable Gaussian augmentation has also been extended to bipartite networks, yielding distributed Gibbs samplers with non-asymptotic guarantees for structured log-concave targets (Yuan et al., 2023; Yuan et al., 2025a). The proximal viewpoint has also been developed beyond smooth, strongly log-concave Euclidean targets. Alternating-sampling and proximal-oracle implementations address nonsmooth, weakly smooth, and composite potentials (Liang and Chen, 2022; Liang and Chen, 2023; Liang and Chen, 2026). Proximal sampler with adaptive stepsize is explored in Yuan et al. (2025b). Non-Euclidean variants replace the Gaussian geometry by log-Laplace regularization (Gopi et al., 2023), while related algorithmic-diffusion constructions treat uniform sampling from convex bodies (Kook et al., 2024). For heavy-tailed targets, Gaussian proximal oracles can face an intrinsic high-accuracy barrier, motivating stable-oracle analogues (He et al., 2024). Recent strong-data-processing analyses also establish convergence of the proximal sampler in relative Fisher information (Wibisono, 2026). These works broaden the target, geometry, oracle, or convergence criterion; our methods instead retain the smooth strongly log-concave setting and modify the outer dynamics to carry directional memory. The closest comparison is concurrent PBPS (Chen et al., 2026b). PBPS uses the same augmentation and outer reflection, but its conditional transition is a harmonic BPS half-turn with specular bounces and a fresh independent reference from qyq_y at every iteration. PBPS and both of our methods have O~(κd1/4) O( κd^1/4) warm-start oracle-query complexity. The flux method replaces specular reflection by Gaussian-flux resampling and adds underdamped forcing. The persistent method instead uses a quarter-turn, carries momentum between ordinary iterations, and calls a conditional-resampling subroutine only on rare refresh iterations. Table 2: Comparison of reflected proximal samplers from a constant-D2D_2 warm start. PBPS Flux half-turn Persistent quarter-turn Conditional motion harmonic BPS underdamped, flux-resampled fixed-reference bouncy Momentum across outer steps discarded discarded retained Independent qyq_y draw every step every step refresh steps only Oracle queries O~(κd1/4) O( κd^1/4) O~(κd1/4) O( κd^1/4) O~(κd1/4) O( κd^1/4) Conditional-resampling calls O~(κd1/4) O( κd^1/4) O~(κd1/4) O( κd^1/4) polylog(d,κ,1/ε)polylog(d,κ,1/ ) Appendix B Proofs for the Flux-Resampled Half-Turn This appendix proves Props. 2, 5, 1 and 2. The proof follows the algorithm: 1. establish the geometry of the localized conditional law and show that the event-driven half-turn preserves it; 2. prove that one conditional half-turn approximately reverses every affine position coordinate that is relevant to the outer chain; 3. combine this estimate with the outer reflection and occasional conditional refreshment to obtain accelerated L2L^2 contraction; 4. replace exact conditional draws and unbounded event clocks by finite implementations, then account for total-variation error and oracle queries. Standard ingredients are cited rather than reproved: Brascamp–Lieb and Poincaré inequalities for strongly log-concave measures, skew balance for velocity-jump processes, and the existing approximate-RGO implementation. The new parts are the covariance-weighted half-turn estimate and the modified-energy argument that turns it into contraction of the complete reflected proximal chain. Throughout this appendix, set a=αη,s=βη,gy(w)=η∇V(y+ηw).a=αη, s=βη, g_y(w)= η\,∇ V(y+ ηw). Thus 0<a≤s0<a≤ s. Under the change of variables w=(x−y)/ηw=(x-y)/ η, the conditional density qyq_y becomes qyw(w)∝exp−∥w∥2/2−V(y+ηw)q_y^w(w) \- w ^2/2-V(y+ ηw)\. We write I for the identity on the space at hand, use IdI_d when its matrix dimension is worth displaying, and set γd=(0,Id) _d=N(0,I_d). Unless a measure is displayed explicitly, inner products and norms below are taken in the stationary L2L^2 space of the objects under discussion; gradients of functions of y are with respect to y. Operators act on observables; probability densities therefore evolve by their adjoints. The few operators used in the contraction proof are: symbol meaning P conditional expectation given Y =I− Q=I- P remaining conditional variation U outer reflection (x,y)↦(x,2x−y)(x,y) (x,2x-y) F conditional flux half-turn K one complete ideal iteration = T= P U P reflection seen on Y-only observables Γ=(I−2)1/2 =(I- T^2)^1/2 amount transferred out of the Y-only part B.1 Conditional Geometry and Validity of the Half-Turn We first collect the conditional identities used throughout the proof and then verify that the event-driven motion is an exact conditional update. Neither step requires a classical Hessian: the identities hold in the weak sense for the strongly convex, smooth potential in Assump. 1. Let W=(X−y)/ηW=(X-y)/ η under qyq_y, and denote its mean and covariance by w¯y w_y and Σy _y. The following lemma says that the departure of this conditional law from a standard Gaussian is measured by the same matrix I−ΣyI- _y that controls the curvature of the Y-marginal. Lemma 1 (Localized conditional geometry). For every y, 11+sI⪯Σy⪯11+aI,a1+aI⪯I−Σy⪯s1+sI. 11+sI _y 11+aI, a1+aI I- _y s1+sI. If ΠηY(y)∝e−Uη(y) _η^Y(y) ^-U_η(y), then ∇Uη(y)=η−1y−[X∣Y=y],∇2Uη(y)=η−1(I−Σy),∇ U_η(y)=η^-1\y-E[X Y=y]\, ∇^2U_η(y)=η^-1(I- _y), and hence α(1+a)−1I⪯∇2Uη⪯β(1+s)−1Iα(1+a)^-1I ∇^2U_η β(1+s)^-1I. Moreover, when s<1s<1, gy(W)=−w¯y,Cov(W,gy(W))=I−Σy,Cov(gy(W))⪯s1−s(I−Σy).Eg_y(W)=- w_y, (W,g_y(W))=I- _y, (g_y(W)) s1-s(I- _y). (11) Proof. The normalized conditional potential has curvature between (1+a)I(1+a)I and (1+s)I(1+s)I. Brascamp–Lieb gives the upper covariance bound. The information inequality gives Cov(W)⪰∇2[−logqyw(W)]−1⪰(1+s)−1ICov(W) \E∇^2[- q_y^w(W)]\^-1 (1+s)^-1I, proving the remaining covariance bounds. Differentiating the conditional normalizer gives the two identities for UηU_η. For 11, use ∇logqyw(w)=−w−gy(w)∇ q_y^w(w)=-w-g_y(w). Integration by parts first gives (W+gy(W))=0E(W+g_y(W))=0 and Σy+Cov(W,gy(W))=I _y+Cov(W,g_y(W))=I. Applying it once more to gyg_y gives ∇gy(W)=I−Σy+Cov(gy(W)).E∇ g_y(W)=I- _y+Cov(g_y(W)). The conditional law has Poincaré constant at most one, while 0⪯∇gy⪯sI0 ∇ g_y sI. Therefore Cov(gy(W))⪯s∇gy(W)Cov(g_y(W)) sE∇ g_y(W); rearrangement proves the last claim. These standard Brascamp–Lieb and integration-by-parts facts are reviewed, for example, by Chewi (2023). The extension from smooth potentials and test functions to the regularity in Assump. 1 is collected in Sec. B.4.2. ∎ We next verify the conditional process. Fix y and a reference w~ w. To lighten notation write g=gyg=g_y, bref=g(w~)b_ ref=g( w), and h(w)=g(w)−brefh(w)=g(w)-b_ ref. In normalized coordinates, the generator of the diffusion between events is ℒ0=pT∇w−(w+bref)T∇p+γ(Δp−pT∇p),L_0=p^T _w-(w+b_ ref)^T _p+γ( _p-p^T _p), and events occur at rate [pTh(w)]+[p^Th(w)]_+. At an event, write n=h/∥h∥n=h/ h , decompose p=v+znp=v+zn, draw R>0R>0 with density re−r2/2re^-r^2/2, and replace p by v−Rnv-Rn. When h=0h=0, the event rate is zero and the direction n is never evaluated. More explicitly, if J denotes the event part of the generator, then (f)(w,p)=[pTh(w)]+∫0∞[f(w,p−pTn(w)+rn(w))−f(w,p)]re−r2/2dr,(Jf)(w,p)=[p^Th(w)]_+ _0^∞\! [f\! (w,p-\p^Tn(w)+r\n(w) )-f(w,p) ]re^-r^2/2\,dr, with f=0Jf=0 when h(w)=0h(w)=0. The full generator used below is ℒ=ℒ0+L=L_0+J. Proposition 4 (Fixed-reference validity). For every fixed (y,w~)(y, w), the preceding process is non-explosive and preserves qyw⊗γdq_y^w _d. If ty,w~ S_t^y, w is its transition operator and (Jf)(w,p)=f(w,−p)(Jf)(w,p)=f(w,-p), then (ty,w~)∗=Jty,w~J.( S_t^y, w)^*=J S_t^y, wJ. It is a strongly continuous contraction semigroup on L2(qyw⊗γd)L^2(q_y^w _d). Consequently, inserting an independent Gaussian momentum, running for a fixed time, and discarding momentum gives a reversible position kernel for qywq_y^w. The conclusion remains true after averaging over any reference law independent of the incoming (W,P)(W,P) conditional on Y=yY=y. Proof. The calculation follows the usual skew-balance argument for event-driven kinetic samplers (Vanetti et al., 2017; Bierkens et al., 2020; Michel et al., 2020; Monmarché et al., 2023); we give the semigroup step because our reference motion also contains Langevin noise. The proof has three steps. First, a trajectory with no event can be reversed after flipping momentum; its no-event probability supplies exactly the likelihood correction caused by the residual force. Second, the Rayleigh update balances the stationary probability flow through a single collision. Finally, we sum these identities over the number and times of the collisions. This proves reversal of the complete process as well as invariance. Non-explosion. With =1+∥w+bref∥2+∥p∥2/2V=1+\ w+b_ ref ^2+ p ^2\/2, the full generator obeys ℒ≤Cw~(1+)LV≤ C_ w(1+V), and the event rate is at most Cw~(1+)C_ w(1+V), where Cw~C_ w is a finite constant that may depend on the fixed reference. Indeed, a jump changes V by (R2−z2)/2(R^2-z^2)/2, R2=2ER^2=2, and ∥h(w)∥≤s∥w−w~∥ h(w) ≤ s w- w . A stopped Dynkin formula and Gronwall bound both (Wt,Pt)EV(W_t,P_t) and the expected event count on every finite interval. The recursive integrated-hazard construction is therefore well defined and non-explosive. Reversal between and at events. Let Zt0=(Wt0,Pt0)Z_t^0=(W_t^0,P_t^0) be the affine Langevin diffusion without events, let λ(w,p)=[pTh(w)]+λ(w,p)=[p^Th(w)]_+, and define its killed semigroup—that is, its transition weighted by the probability of seeing no event—by tf(z)=z[e−∫0tλ(Zu0)duf(Zt0)]. D_tf(z)=E_z\! [e^- _0^tλ(Z_u^0)\,duf(Z_t^0) ]. The affine diffusion is invariant and momentum-reversible for the Gaussian density μref∝exp−∥w+bref∥2/2−∥p∥2/2 _ ref \- w+b_ ref ^2/2- p ^2/2\. For the desired density μ=qyw⊗γdμ=q_y^w _d, the ratio ϱ=dμ/dμref =dμ/d _ ref satisfies ∇logϱ=−h∇ =-h. Since the position path is absolutely continuous, ϱ(Wt0)ϱ(W00)=exp−∫0tPu0Th(Wu0)du=exp−∫0t(λ−λ∘J)(Zu0)du. (W_t^0) (W_0^0)= \- _0^tP_u^0Th(W_u^0)\,du \= \- _0^t(λ-λ J)(Z_u^0)\,du \. Stationary path reversal for the affine diffusion therefore gives t∗=JtJ D_t^*=J D_tJ first for bounded nonnegative functions. Taking one test function equal to one shows μt≤μ D_t≤μ, and Jensen’s inequality extends the identity to L2(μ)L^2(μ). Let G be the gain kernel: it applies one collision, weighted by the rate at which that collision occurs. At fixed w, write the incoming momentum as p=v+znp=v+zn, where v⟂nv n and z>0z>0. Apart from a common constant, the stationary marked flux is qyw(w)∥h(w)∥φd−1(v)ze−z2/2re−r2/2dwdvdzdr,z,r>0.q_y^w(w) h(w) _d-1(v)ze^-z^2/2re^-r^2/2\,dw\,dv\,dz\,dr, z,r>0. Here φd−1 _d-1 is the standard Gaussian density on the subspace orthogonal to n. It is unchanged by (v,z,r)↦(−v,r,z)(v,z,r) (-v,r,z), which turns a marked transition into the momentum reversal of its time-reversed transition. Thus ∗=JJG^*=JGJ in the corresponding kernel sense. The contribution of paths with exactly m events is the time-simplex integral of t0t1⋯tm D_t_0G D_t_1·sG D_t_m. Reversing the simplex coordinates and using the preceding two identities shows that its adjoint is obtained by conjugating with J. Tonelli justifies this calculation for nonnegative tests, and non-explosion makes the sum over m conservative, meaning that it retains total probability one. Summing proves (ty,w~)∗=Jty,w~J( S_t^y, w)^*=J S_t^y, wJ; taking one test function equal to one proves invariance. Jensen gives L2L^2 contractivity, while stochastic continuity and density of bounded continuous compactly supported functions give strong continuity. Finally, if ιf(w,p)=f(w) f(w,p)=f(w) and u(w)=∫u(w,p)γd(p) Mu(w)= u(w,p)\, _d(dp), then =ι∗ M= ^*, Jι=ιJ = , and J= MJ= M. Thus inserting a Gaussian momentum and then discarding it are adjoint operations. Suppressing the fixed superscripts (y,w~)(y, w), (tι)∗=tι( M S_t )^*= M S_t . Conditioning on an independent reference and applying Fubini proves the last assertion. ∎ Proof of Prop. 2. Use Prop. 4 with W~∼qyw W q_y^w, independently of the incoming position. It shows that the conditional position kernel in Sec. 4.1 is reversible and invariant for qyq_y. The outer map (x,y)↦(x,2x−y)(x,y) (x,2x-y) preserves Πη _η, because it only changes the sign of y−xy-x and has unit absolute Jacobian. Proposition 1 now proves invariance of the complete row and of its position marginal π. ∎ B.2 The Covariance-Weighted Half-Turn Estimate The next proposition is the local estimate used by the outer-chain proof. It does not ask the conditional update to approximate a reflection uniformly over all functions. It asks only that centered affine position coordinates be nearly negated, with a smaller error in directions where the localized conditional is closer to the Gaussian reference. Proposition 5 (Conditional affine half-turn). There is a universal C such that, whenever s=βη≤1/2s=βη≤ 1/2, the normalized-coordinate position kernel FywF_y^w of Sec. 4.1 satisfies, uniformly in y,u∈ℝdy,u ^d, ∥(I+Fyw)ℓu,y∥L2(qyw)2≤CsuT(I−Σy)u,ℓu,y(w)=uT(w−w¯y). (I+F_y^w) _u,y _L^2(q_y^w)^2≤ Cs\,u^T(I- _y)u, _u,y(w)=u^T(w- w_y). (12) The proof has two ingredients. First, the damped harmonic flight can be solved exactly on affine observables. At the chosen half-period, dependence on the fresh initial momentum vanishes and the incoming position receives a negative coefficient. Second, the event correction appears as a forcing term whose second moment is controlled by s(I−Σy)s(I- _y) through 11. Fix y and a reference w~ w, and retain the notation g,bref,hg,b_ ref,h from Sec. B.1. Let r¯=R=π/2 r=ER= π/2, ψ(z)=z+(z+r¯)ψ(z)=z_+(z+ r), and n(w)=h(w)/∥h(w)∥n(w)=h(w)/ h(w) when h(w)≠0h(w)≠ 0. For affine position and momentum observables ℓu,y(w)=uT(w−w¯y) _u,y(w)=u^T(w- w_y) and mu(p)=uTpm_u(p)=u^Tp, define ku,w~(w,p)=−uTw¯y+bref−uTh(w)ψ(pTn(w)),k_u, w(w,p)=-u^T\ w_y+b_ ref\-u^Th(w)ψ(p^Tn(w)), where the second term is zero when h(w)=0h(w)=0. Lemma 2 (Affine evolution). The two observables belong to the L2(qyw⊗γd)L^2(q_y^w _d) generator domain and satisfy ℒℓu,y=mu,ℒmu=−ℓu,y−γmu+ku,w~.L _u,y=m_u, _u=- _u,y-γ m_u+k_u, w. Writing ω=1−γ2/4ω= 1-γ^2/4, set Aγ(t)=e−γt/2cos(ωt)+γ2ωsin(ωt),Bγ(t)=e−γt/2sin(ωt)ω.A_γ(t)=e^-γ t/2 \ (ω t)+ γ2ω (ω t) \, B_γ(t)=e^-γ t/2 (ω t)ω. Then, for the fixed-reference semigroup, ty,w~ℓu,y=Aγ(t)ℓu,y+Bγ(t)mu+∫0tBγ(t−r)ry,w~ku,w~r. S_t^y, w _u,y=A_γ(t) _u,y+B_γ(t)m_u+ _0^tB_γ(t-r) S_r^y, wk_u, w\,dr. At Tγ=π/ωT_γ=π/ω, this becomes (I+Tγy,w~)ℓu,y=(1−e−γTγ/2)ℓu,y+∫0TγBγ(Tγ−r)ry,w~ku,w~dr.(I+ S_T_γ^y, w) _u,y=(1-e^-γ T_γ/2) _u,y+ _0^T_γB_γ(T_γ-r) S_r^y, wk_u, w\,dr. (13) Proof. The deterministic part of the generator gives ℒℓu,y=muL _u,y=m_u and −ℓu,y−γmu−uT(w¯y+bref)- _u,y-γ m_u-u^T( w_y+b_ ref) in the momentum coordinate. At an event, the momentum change is −(Z+R)uTn-(Z+R)u^Tn, where Z=pTn>0Z=p^Tn>0; multiplying its conditional mean by the rate ∥h∥Z h Z gives the remaining term −uThψ(Z)-u^Th\,ψ(Z). These quantities are square integrable because qywq_y^w has Gaussian tails, h grows at most linearly, and pTnp^Tn is standard normal conditional on the position. Applying the stopped Dynkin formula and then removing the stopping time using the Lyapunov bound in Prop. 4 gives, after suppressing the fixed superscripts (y,w~)(y, w), tℓu,y−ℓu,y=∫0trmur,tmu−mu=∫0tr(−ℓu,y−γmu+ku,w~)r. S_t _u,y- _u,y= _0^t S_rm_u\,dr, S_tm_u-m_u= _0^t S_r(- _u,y-γ m_u+k_u, w)\,dr. Strong continuity and division by t prove the generator-domain statements. Solving the resulting two-dimensional linear equation by variation of constants gives the displayed formula. At TγT_γ, Bγ(Tγ)=0B_γ(T_γ)=0 and Aγ(Tγ)=−e−γTγ/2A_γ(T_γ)=-e^-γ T_γ/2, which proves 13. ∎ Proof of Prop. 5. Let W,W~W, W be independent with law qywq_y^w, and average first over the Gaussian momentum; we denote this last expectation by PE_P. Direct one-dimensional Gaussian integration gives ψ(Z)=1Eψ(Z)=1 and Varψ(Z)=5/2+π/4Varψ(Z)=5/2+π/4 for Z∼(0,1)Z (0,1). Using gy(W)=−w¯yEg_y(W)=- w_y, we obtain Pku,W~=−uTgy(W)−gy(W),E_Pk_u, W=-u^T\g_y(W)-Eg_y(W)\, while the conditional momentum variance is (5/2+π/4)[uTgy(W)−gy(W~)]2(5/2+π/4)[u^T\g_y(W)-g_y( W)\]^2. Independence of W and W~ W therefore gives ku,W~2=(6+π2)Var(uTgy(W))≤CsuT(I−Σy)u,Ek_u, W^2= (6+ π2 )Var(u^Tg_y(W))≤ Cs\,u^T(I- _y)u, (14) where the last step is 11. Insert fresh Gaussian momentum in 13, average the independent reference, and discard the terminal momentum. Jensen’s inequality, semigroup contractivity, and ∫0Tγ|Bγ(r)|r≤C _0^T_γ|B_γ(r)|\,dr≤ C yield ∥(I+Fyw)ℓu,y∥22≤Cγ2uTΣyu+ku,W~2. (I+F_y^w) _u,y _2^2≤ C \γ^2u^T _yu+Ek_u, W^2 \. Here the norm is in L2(qyw)L^2(q_y^w), and the remaining expectation includes the independent position, reference, and Gaussian momentum. The parameter choice in Sec. 4 gives γ2≤as/(1+a)γ^2≤ as/(1+a). Since I−Σy⪰a(1+a)−1I- _y a(1+a)^-1I and Σy⪯I _y I, the damping term is at most suT(I−Σy)usu^T(I- _y)u. Combining this with 14 proves 12. ∎ B.3 From the Conditional Half-Turn to Accelerated Contraction We now prove the ideal-chain theorem. The argument follows two components of an observable: its average at fixed Y, and its remaining conditional variation. Outer reflection transfers part of the first component into the second. The flux half-turn approximately reverses precisely this transferred part, while an occasional independent conditional refresh damps all conditional variation. Let ℌ=L2(Πη) H=L^2( _η), let f=[f∣Y] Pf=E[f Y], and set =I− Q=I- P. Let U be the pullback of the outer reflection and F the joint conditional half-turn. By Prop. 4, =+Q F= P+ F_Q, where Q F_Q is a self-adjoint contraction on ran()ran( Q). The first term records observables depending only on Y; the second acts on their remaining conditional variation. Thus one ideal iteration has observable operator =+(1−ρ)Q. K= U\ P+(1-ρ) F_Q\. (15) Proposition 6 (Ideal L2L^2 contraction). There are universal constants c,C>0c,C>0 such that, when s≤cs≤ c and ρ≍alog(e/a)ρ a (e/a), the ideal chain satisfies ∥n∥L02(Πη)→L02(Πη)≤3exp−cnalog(e/a). K^n _L^2_0( _η)→ L^2_0( _η)≤ 3 \- cn a (e/a) \. The rest of this subsection proves the proposition. We first quantify what reflection transfers, then lift the conditional estimate in Prop. 5 to the complete state space, and finally construct an energy that decreases in one step. B.3.1 What Outer Reflection Transfers On the centered Y-only subspace ℋY=ran()∩L02(Πη)H_Y=ran( P)∩ L^2_0( _η), define =,Γ=(I−2)1/2,=. T= P U P, =(I- T^2)^1/2, C= Q U P. Under stationarity, let X∼πX π and Z∼(0,Id)Z (0,I_d) be independent, and write Y+=X+ηZY_+=X+ ηZ, Y−=X−ηZY_-=X- ηZ. These are not two independent chain iterates: they share (X,Z)(X,Z) and are the auxiliary values before and after the same reflection, so Y−=2X−Y+Y_-=2X-Y_+. Then (f)(y)=[f(Y−)∣Y+=y]( Tf)(y)=E[f(Y_-) Y_+=y], and ∥f∥2=Var(f(Y−)∣Y+) Cf ^2=EVar(f(Y_-) Y_+). Thus T is the part retained as a function of Y+Y_+, while C is the part transferred into conditional variation. The operator Γ measures the size of this transferred part. Lemma 3 (Outer-reflection geometry). If s≤1/3s≤ 1/3, then T is self-adjoint and, on ℋYH_Y, −12I⪯1−a1+aI,Γ⪰γaI,γa=2a1+a.- 12I T 1-a1+aI, _aI, _a= 2 a1+a. (16) Moreover, 4η∥∇f∥2≤∥Γf∥24η ∇ Tf ^2≤ f ^2, and =Γ C= V for an isometry :ℋY→ran() V:H_Y ( Q). Proof. Exchangeability of (Y+,Y−)(Y_+,Y_-) makes T self-adjoint. The conditional density of Y−Y_- given Y+=yY_+=y is proportional to exp−V(y+u2)−∥y−u∥28η. \-V ( y+u2 )- y-u ^28η \. Its curvature in u is at least (α+η−1)I/4(α+η^-1)I/4, while the derivative with respect to u of its y-score has norm at most (η−1−α)/4(η^-1-α)/4. Conditional Brascamp–Lieb therefore gives η∥∇f∥2≤(1−a)24(1+a)∥f∥2≤14∥Γf∥2.η ∇ Tf ^2≤ (1-a)^24(1+a) Cf ^2≤ 14 f ^2. The Y-marginal has Poincaré constant at most α−1+ηα^-1+η by Lem. 1. Applying this inequality to the centered function f Tf and using ∥f∥2=∥f∥2+∥f∥2 f ^2= Tf ^2+ Cf ^2 gives ∥≤(1−a)/(1+a) T ≤(1-a)/(1+a). For the lower spectral edge, exchangeability and Gaussian Poincaré give ⟨h,(I−)h⟩≤2η∥∇h∥2 h,(I- T)h ≤ 2η ∇ h ^2. Apply this with h=fh= Tf and use the preceding derivative bound. Spectral calculus yields (I−)2(I+2)⪰0(I- T)^2(I+2 T) 0, hence ⪰−I/2 T -I/2. Because (1−a)/(1+a)≥1/2(1-a)/(1+a)≥ 1/2, these two edges imply the stated lower bound on Γ . Finally, unitarity of U gives ∗=I−2=Γ2 C^* C=I- T^2= ^2, so =Γ−1 V= C ^-1 is an isometry. The weak-gradient and spectral-domain justification for this calculation is given in Sec. B.4.2. ∎ Reflection creates a conditional-position signal that is almost affine. For weakly differentiable f∈ℋYf _Y, set by=[∇f(Y−)∣Y+=y],Lf=2ηbY+T(W−w¯Y+),W=X−Y+η.b_y=E[∇ f(Y_-) Y_+=y], L_f=2 η\,b_Y_+^T(W- w_Y_+), W= X-Y_+ η. Here byb_y is the average slope of f(Y−)f(Y_-) conditional on the reflected value Y+=yY_+=y, and LfL_f is the corresponding centered linear signal in the conditional position W. The next lemma separates this leading signal from a nonlinear remainder. Its quantities ℛfR_f and fD_f measure, respectively, the gradient variation not retained by conditional averaging and the curvature-weighted size of the affine signal. Lemma 4 (Affine approximation after reflection). There is an error efe_f, centered conditionally on Y+Y_+, such that f=Lf+ef Cf=L_f+e_f. If ℛf=∥∇f∥2−∥∇f∥2,f=[bY+T(I−ΣY+)bY+],R_f= ∇ f ^2- T∇ f ^2, _f=E[b_Y_+^T(I- _Y_+)b_Y_+], then, for sufficiently small universal s, ∥ef∥2≤4ηℛf,ℛf+f≤C∥∇(I−)f∥∇f∥. e_f ^2≤ 4 _f, _f+D_f≤ C ∇(I- T)f ∇ f . (17) Proof. Here ∇f T∇ f means that T is applied to every coordinate of the vector-valued gradient; in particular, (∇f)(y)=by( T∇ f)(y)=b_y. Conditionally on Y+=yY_+=y, Y−=y+2ηWY_-=y+2 ηW. Subtracting its best first-order conditional approximation gives ef=f(Y−)−f(Y+)−2ηbY+T(W−w¯Y+).e_f=f(Y_-)- Tf(Y_+)-2 η\,b_Y_+^T(W- w_Y_+). Conditional Poincaré proves ∥ef∥2≤4ηℛf e_f ^2≤ 4 _f. For the second estimate, differentiate the conditional expectation weakly. Integration by parts and 11 give ∇(I−)f=∇f−by+2(I−Σy)by+η−1/2Cov(ef,gy(W)∣Y+=y).∇(I- T)f=∇ f-b_y+2(I- _y)b_y+η^-1/2Cov(e_f,g_y(W) Y_+=y). In this identity, ∇f∇ f on the right is evaluated at y, and the conditional covariance is over W given Y+=yY_+=y. Write b(y)=byb(y)=b_y. Take its inner product with ∇f+by∇ f+b_y and average over y. The first term gives ℛfR_f, and the leading covariance term gives 4f4D_f. Also, the spectral lower bound ⪰−I/2 T -I/2 implies ∥∇f−b∥2≤3∥∇f∥2−∥b∥2=3ℛf, ∇ f-b ^2≤ 3\ ∇ f ^2- b ^2\=3R_f, The two remaining cross terms are bounded using this inequality, 11, the first inequality in 17, and I−Σy⪯sI- _y sI. They are at most Csℛff+Csℛf(sℛf+f).C sR_fD_f+C sR_f(sR_f+D_f). The leading inner product is at most 2∥∇(I−)f∥∇f∥2 ∇(I- T)f ∇ f . Young’s inequality absorbs the displayed cross terms when s is sufficiently small. The graph-norm extension to the stated weak-gradient domain is recorded in Sec. B.4.2. ∎ B.3.2 Lifting the Local Half-Turn Estimate The next result is the bridge between the conditional dynamics and the outer Markov chain. We first record the operator comparison that accounts for the logarithm in that bridge. It converts an interpolation-type estimate into a quadratic form estimate without assuming that the two operators commute. Lemma 5 (Dyadic operator comparison). Suppose mI⪯B⪯MImI B MI and a bounded operator A satisfies ∥Av∥2≤∥v∥∥Bv∥ Av ^2≤ v Bv . Then A∗A⪯C1+log(M/m)B.A^*A C\1+ (M/m)\B. Proof. Partition the spectrum of B into dyadic intervals and write v=∑j=1Jvjv= _j=1^Jv_j, where J≤C1+log(M/m)J≤ C\1+ (M/m)\. If λj _j is the lower endpoint of the j-th interval, then the hypothesis gives ∥Avj∥≤2λj∥vj∥ Av_j ≤ 2 _j v_j . Consequently, ∥Av∥2≤J∑j2λj∥vj∥2≤2J⟨v,Bv⟩, Av ^2≤ J _j2 _j v_j ^2≤ 2J v,Bv , which proves the claim. ∎ Proposition 7 (Global half-turn estimate). For sufficiently small universal s and every v∈ℋYv _Y, ∥(I+Q)v∥2≤Clog(e/γa)∥Γv∥2. (I+ F_Q) Vv ^2≤ C (e/ _a) v ^2. (18) Proof. First restrict f to the spectral subspace of T corresponding to [1/2,1][1/2,1]. Apply Prop. 5 conditionally to the affine part LfL_f of f Cf, use ∥I+Q∥≤2 I+ F_Q ≤ 2 on efe_f, and then apply Lem. 4. This gives ∥(I+Q)f∥2≤Cη∥∇(I−)f∥∇f∥. (I+ F_Q) Cf ^2≤ Cη ∇(I- T)f ∇ f . On this spectral subspace, applying the derivative estimate in Lem. 3 to −1f T^-1f and −1(I−)f T^-1(I- T)f gives η∥∇f∥≤C∥Γf∥,η∥∇(I−)f∥≤C∥Γ3f∥. η ∇ f ≤ C f , η ∇(I- T)f ≤ C ^3f . Since =Γ C= V , setting f=Γ−1vf= ^-1v, which is allowed by Γ⪰γaI _aI, shows that ht=(I+Q) A_ ht=(I+ F_Q) V obeys ∥htv∥2≤C∥v∥∥Γ2v∥ A_ htv ^2≤ C v ^2v . Applying Lem. 5 to A=htA= A_ ht, B=Γ2B= ^2, and m=γa2m= _a^2 gives the required bound without treating the generally non-self-adjoint half-turn error ht A_ ht as a positive operator. On the complementary spectral subspace, Γ2⪰3I/4 ^2 3I/4, while ∥ht∥≤2 A_ ht ≤ 2, so the result is immediate. The two output vectors need not be orthogonal; the triangle inequality and ∥r1+r2∥2≤2∥r1∥2+2∥r2∥2 r_1+r_2 ^2≤ 2 r_1 ^2+2 r_2 ^2 prove 18 for an arbitrary v. ∎ B.3.3 A Decreasing Energy For centered f∈ℌf∈ H, write fY=f_Y= Pf, fQ=f_Q= Qf, and fV=∗fQf_V= V^*f_Q. Since V is an isometry, the conditional component has the orthogonal decomposition fQ=fV+f⟂,f⟂=(I−∗)fQ.f_Q= Vf_V+f_ , f_ =(I- V V^*)f_Q. The three pieces have distinct roles. The term fYf_Y is variation across values of Y; fV Vf_V is the conditional signal created from it by reflection; and f⟂f_ is the remaining conditional variation. The cross term below tracks the exchange between the first two pieces, whereas the RGO branch controls the full fQf_Q, including f⟂f_ . Define ℰϑ(f)=∥f∥2+ϑ(fY,fV),(fY,fV)=12(∥fY∥2−∥fV∥2)−⟨Γ−1fY,fV⟩.E_ (f)= f ^2+ (f_Y,f_V), (f_Y,f_V)= 12( f_Y ^2- f_V ^2)- T ^-1f_Y,f_V . The form C is chosen so that, under an exact conditional sign reversal, the pair (fY,fV)(f_Y,f_V) undergoes the rotation displayed below and C decreases by ∥fY∥2−∥fV∥2 f_Y ^2- f_V ^2. Lemma 6 (One-step energy decay). Let Λρ=log(e/γa)+(ρ/γa)2,ϑ=c⋆ρ/Λρ _ρ= (e/ _a)+(ρ/ _a)^2, =c_ ρ/ _ρ for a sufficiently small universal c⋆c_ . If 0<ρ≤1/20<ρ≤ 1/2, then 12∥f∥2≤ℰϑ(f)≤32∥f∥2,ℰϑ(f)≤(1−cρΛρ)ℰϑ(f). 12 f ^2 _ (f)≤ 32 f ^2, _ ( Kf)≤ (1-c ρ _ρ )E_ (f). Proof. Because T commutes with Γ , the block operator associated with 22C has square diag(Γ−2,Γ−2)diag( ^-2, ^-2). Hence ||≤(2γa)−1(∥fY∥2+∥fV∥2) ≤(2 _a)^-1( f_Y ^2+ f_V ^2). The displayed choice of ϑ satisfies ϑ≤γa ≤ _a, proving energy equivalence. From 15 and contractivity of Q F_Q, ∥f∥2−∥f∥2≤−ρ∥fQ∥2. Kf ^2- f ^2≤-ρ f_Q ^2. (19) To see how the correction changes, first replace the conditional move by an exact sign reversal and set g=(−)fg= U( P- Q)f. Unitarity of U and =Γ C= V give (gYgV)=(−Γ)(fYfV),(gY,gV)−(fY,fV)=−∥fY∥2+∥fV∥2. g_Yg_V= pmatrix T&- \\ & T pmatrix f_Yf_V, (g_Y,g_V)-C(f_Y,f_V)=- f_Y ^2+ f_V ^2. The actual conditional move differs from this sign reversal by eht=∗(I+Q)fQ,eρ=ρfV+(1−ρ)eht.e_ ht= V^*(I+ F_Q)f_Q, e_ρ=ρ f_V+(1-ρ)e_ ht. The two tracked components of the actual update are (f)Y=gY+Γeρ( Kf)_Y=g_Y+ e_ρ and (f)V=gV−eρ( Kf)_V=g_V- Te_ρ. Direct substitution gives ((f)Y,(f)V)−(fY,fV)=−∥fY∥2+∥fV∥2+⟨gY,Γ−1eρ⟩+12∥eρ∥2.C(( Kf)_Y,( Kf)_V)-C(f_Y,f_V)=- f_Y ^2+ f_V ^2+ g_Y, ^-1e_ρ + 12 e_ρ ^2. By Prop. 7, ht∗ht⪯Clog(e/γa)Γ2 A_ ht^* A_ ht C (e/ _a) ^2 for ht=(I+Q) A_ ht=(I+ F_Q) V. The corresponding operator factorization gives a bounded B, with ∥B∥2≤Clog(e/γa) B ^2≤ C (e/ _a), such that ht=BΓ A_ ht=B . Since eht=ht∗fQe_ ht= A_ ht^*f_Q, this gives ∥Γ−1eht∥2≤Clog(e/γa)∥fQ∥2 ^-1e_ ht ^2≤ C (e/ _a) f_Q ^2. Therefore ∥Γ−1eρ∥2≤CΛρ∥fQ∥2 ^-1e_ρ ^2≤ C _ρ f_Q ^2, while ∥eρ∥≤2∥fQ∥ e_ρ ≤ 2 f_Q . Young’s inequality now yields ((f)Y,(f)V)−(fY,fV)≤−12∥fY∥2+CΛρ∥fQ∥2.C(( Kf)_Y,( Kf)_V)-C(f_Y,f_V)≤- 12 f_Y ^2+C _ρ f_Q ^2. Multiply by ϑ , add 19, and choose c⋆c_ small enough. Energy equivalence then proves the claimed one-step decay. ∎ Proof of Prop. 6. By 16, γa≍a _a a. With ρ≍alog(e/a)ρ a (e/a), we have Λρ≲log(e/a) _ρ (e/a) and ρ/Λρ≳a/log(e/a)ρ/ _ρ a/ (e/a). Iterating Lem. 6 and using energy equivalence proves the proposition. ∎ Proof of Thm. 1. Choose nblkn_blk as in the theorem, with a sufficiently large universal constant. Proposition 6 gives ∥nblk∥L02→L022≤e−1 K^n_blk _L^2_0→ L^2_0^2 ^-1. If ϱν=dν/dΠη _ν=dν/d _η, then the density evolves by ∗ K^* and χ2(ν∥Πη)=∥ϱν−1∥2χ^2(ν\| _η)= _ν-1 ^2. Since an operator and its adjoint have the same norm, applying the block estimate J times proves 3. ∎ B.4 Implementable Algorithm and Query-Complexity Analysis The ideal chain uses an exact conditional reference and an uncapped event clock. Its finite counterpart makes the following three replacements: source ideal operation finite operation reference exact draw from qyq_y certified approximate RGO event clock full state-dependent rate capped Poisson thinning center preparation unrestricted solve capped solver with early return The proposition is stated first; the rest of the subsection controls the three corresponding disagreement events and counts their oracle queries. Proposition 8 (Finite compiler at a fixed step size). Let D2(ν0∥π)≤ΔD_2( _0\|π)≤ , where Δ≥1 ≥ 1, and let 0<ε<1/40< <1/4. For 0<η≤1/β0<η≤ 1/β, define ℓ=KΔ+log(Kdκεαη). =K \ + ( Kdκ αη ) \. There are universal constants c,K>0c,K>0 such that, if sdℓ+ℓ≤c,s\ d + \≤ c, (20) the finite flux algorithm returns X X with TV(Law(X^),π)≤εTV(Law( X),π)≤ and Nqry≤Kalog(e/a)(Δ+log4ε)log(eκℓ)+s(dℓ+ℓ).EN_qry≤ K a (e/a) ( + 4 ) \ (eκ )+s( d + ) \. (21) B.4.1 Capped Poisson thinning For fixed (y,x~)(y, x), put xctr=y−η∇V(x~)x_ ctr=y-η∇ V( x) and U=(X−xctr)/ηU=(X-x_ ctr)/ η. Between events, (U,P)(U,P) follows a linear Gaussian diffusion. With ω=1−γ2/4ω= 1-γ^2/4, its conditional mean after time t is e−γt/2(cos(ωt)+γ2ωsin(ωt)sin(ωt)ω−sin(ωt)ωcos(ωt)−γ2ωsin(ωt))(UP),e^-γ t/2 pmatrix (ω t)+ γ2ω (ω t)& (ω t)ω\\ - (ω t)ω& (ω t)- γ2ω (ω t) pmatrix UP, coordinatewise. Its covariance is obtained by integrating the transported Gaussian noise over this same linear motion. Both mean and covariance are explicit finite-dimensional operations. Thus, once ∇V(x~)∇ V( x) has been stored, no oracle query is needed between candidate event times. Given a cap λ¯ λ, candidate times form a homogeneous Poisson process of rate λ¯ λ. At a candidate time, one oracle query computes λt=η[PtT∇V(Xt)−∇V(x~)]+, _t= η[P_t^T\∇ V(X_t)-∇ V( x)\]_+, and the candidate is accepted with probability min1,λt/λ¯ \1, _t/ λ\. This is exact thinning for the truncated rate minλt,λ¯ \ _t, λ\. Lemma 7 (Stationary rate tail and clock coupling). Conditionally on Y=yY=y, let X,X~X, X be independent draws from qyq_y, and let P∼γdP _d be independent. Then, for universal c,C>0c,C>0 and every u≥0u≥ 0, Prη[PT∇V(X)−∇V(X~)]+>Cs(du+u)≤Ce−cu. \! \ η[P^T\∇ V(X)-∇ V( X)\]_+>Cs( du+u) \≤ Ce^-cu. If, in addition, sds d is below a universal constant, then an augmented initial law ν satisfying D2(ν∥Πη)≤ΔD_2(ν\| _η)≤ admits the following consequence: for every ξ∈(0,1/2)ξ∈(0,1/2), the exact and capped half-turns can be coupled to agree with probability at least 1−ξ1-ξ, using λ¯ξ,Δ=Csdℓξ,Δ+ℓξ,Δ,ℓξ,Δ=CΔ+log(1/ξ). λ_ξ, =Cs\ d _ξ, + _ξ, \, _ξ, =C\ + (1/ξ)\. The capped half-turn uses exactly Tγλ¯ξ,ΔT_γ λ_ξ, candidate events in expectation. Proof. The conditional law is (α+η−1)(α+η^-1)-strongly log-concave. Norm concentration and its second-moment bound show that, outside an event of probability 2e−u2e^-u, both ∥X−X∥ X-EX and ∥X~−X∥ X-EX are at most Cη(d+u)C η(d+u). Smoothness then bounds the conditional standard deviation of the scalar product with P by Cs(d+u)Cs( d+ u). A Gaussian tail bound proves the displayed estimate. Integrating it above Cs(dℓ+ℓ)Cs( d + ) gives (λ−λ¯)+≤Cs(d+ℓ+1)e−cℓ.E(λ- λ)_+≤ Cs( d+ +1)e^-c . Under Πη _η, invariance from Prop. 4 makes this bound valid at every time. Couple the exact and capped clocks with the same Poisson random measure; they separate only if an excess-rate event occurs. Taking ℓ=CΔ+log(1/ξ) =C\ + (1/ξ)\ makes the stationary separation probability at most e−Δξ2e^- ξ^2. Under a D2D_2-warm initial law, the path likelihood ratio is only its initial density ratio, since all subsequent references, Gaussian noise, clocks, and marks use common kernels. Cauchy–Schwarz therefore raises the stationary failure probability by at most eΔ/2e /2, giving ξ. The final expectation is that of a homogeneous Poisson clock on an interval of length Tγ=O(1)T_γ=O(1). ∎ Approximate conditional draws. We use the following existing RGO result. Proposition 9 (Approximate RGO). There is a universal KrgoK_ rgo such that the following holds. Given y, 0<δ<1/20<δ<1/2, and a point x x satisfying ∥y−η∇V(x^)−x^∥≤dη, y-η∇ V( x)- x ≤ dη, if η−1≥Krgoβdlog(1/δ)+log(1/δ),η^-1≥ K_ rgoβ\ d (1/δ)+ (1/δ)\, the rejection sampler of Chen et al. (2026a, Theorems 3.1 and D.1) returns a law q^y q_y satisfying χ2(q^y∥qy)≤δ2,χ2(qy∥q^y)≤δ2,χ^2( q_y\|q_y)≤δ^2, χ^2(q_y\| q_y)≤δ^2, using at most KrgoK_ rgo expected oracle queries after construction of the displayed certificate. The certificate itself uses only ordinary oracle queries. The conditional potential Vy(x)=V(x)+∥x−y∥2/(2η)V_y(x)=V(x)+ x-y ^2/(2η) has condition number at most two. Gradient descent with step 2/(α+β+2η−1)2/(α+β+2η^-1), initialized at the incoming position x, reaches and certifies the required residual in at most C[1+log1+∥y−η∇V(x)−x∥dη]C [1+ \1+ y-η∇ V(x)-x dη \ ] (22) oracle queries, by standard linear convergence for strongly convex smooth optimization. The implementable algorithm and error allocation. Use ρ from 2, and set nblk=⌈Kalog(e/a)⌉,M=⌈K(Δ+log4ε)⌉nblk,δrow=ε6M.n_blk= K a (e/a) , M= K ( + 4 ) n_blk, _ row= 6M. Use the cap λ¯=Ks(dℓ+ℓ) λ=Ks( d + ) and allow at most Qctr=⌈Klog(eκℓ)⌉Q_ ctr= K (eκ ) center-solver queries per row. Initialize X0∼ν0X_0 _0, Y0|X0∼(X0,ηId)Y_0 X_0 (X_0,η I_d). Each row reflects Y, runs the capped center solver, draws an approximate reference from Prop. 9 with tolerance δrow _ row, and then uses that draw either as the independent refresh or as the reference of the capped flux half-turn. If the center solver does not certify its residual within its cap, the algorithm stops and returns the current position. This failure rule, together with the almost-sure finiteness of the rejection sampler and Poisson clock, makes the procedure almost surely finite from every initial state. Proof of Prop. 8. Augmentation preserves warmness: dLaw(X0,Y0)dΠη(x,y)=dν0dπ(x),D2(Law(X0,Y0)∥Πη)≤Δ. dLaw(X_0,Y_0)d _η(x,y)= d _0dπ(x), D_2(Law(X_0,Y_0)\| _η)≤ . Theorem 1 and the choice of M put the ideal output within ε/2 /2 of π. Because every ideal row preserves Πη _η, data processing keeps its intermediate D2D_2 divergence at most Δ . It remains to couple every finite row to its ideal counterpart. There are three possible disagreements. Center cap. Under stationarity, write Y=X+ηZY=X+ ηZ. After reflection, the normalized residual at the incoming position is ∥ηZ+η∇V(X)∥dη. ηZ+η∇ V(X) dη. Strong log-concavity, Gaussian concentration, and ∥∇V(X)∥≤β∥X−x⋆∥ ∇ V(X) ≤β X-x , where x⋆x minimizes V, show that it exceeds 2κ(1+2u/d)2 κ(1+ 2u/d) with probability at most 2e−u2e^-u. Taking u=Δ+2log(1/δrow)+log2u= +2 (1/ _ row)+ 2 and using 22 shows that the stationary cap-failure probability is at most e−Δδrow2e^- _ row^2. The warm-path Cauchy–Schwarz argument used in Lem. 7 makes it at most δrow _ row under the current ideal law. Approximate reference. The definitions imply log(1/δrow)≤Cℓ (1/ _ row)≤ C . The gate 20 therefore implies the small-step-size condition in Prop. 9. Either directed χ2χ^2 bound implies total variation at most δrow/2 _ row/2, so maximal coupling matches the approximate and exact references except with probability at most δrow _ row. Because all laws have Borel densities on ℝdR^d, these maximal couplings may be selected jointly measurably in y. On agreement, the same draw is used in both chains. Capped clock. The same parameter bound gives ℓδrow,Δ≤Cℓ _ _ row, ≤ C . After increasing K, Lem. 7, applied conditionally on the common pre-row history, couples the exact and capped half-turns with failure probability at most δrow _ row. A union bound over three events and M rows is at most ε/2 /2. Together with the ideal mixing error, this proves the total-variation claim. Early termination only reduces work. Each row uses at most QctrQ_ ctr center queries, a universal expected number of RGO queries, one query to cache ∇V(X~)∇ V( X) on the half-turn branch, and Tγλ¯T_γ λ expected event-test queries. Hence Nqry≤CMQctr+s(dℓ+ℓ),EN_qry≤ CM\Q_ ctr+s( d + )\, which is 21. ∎ Proof of Thm. 2. Take η and L as in 4, with the universal constant c sufficiently small. Then a=cκ(dL+L),s=cdL+L,ℓ≤CL,a= cκ( dL+L), s= c dL+L, ≤ CL, so 20 holds. Moreover, 1alog(e/a)≤Cκ(d+L)1/4L3/4, 1 a (e/a)≤ C κ\,(d+L)^1/4L^3/4, and the braces in 21 are at most Clog(eκL)C (eκ L). Substitution proves 5 and the total-variation conclusion. ∎ For the warm-start removal argument we also record the dependence on a nonconstant divergence budget. Corollary 1 (General warmness budget). Let B≥1B≥ 1, D2(ν0∥π)≤BD_2( _0\|π)≤ B, and set LB=B+logCdκε,AB=dLB+LB,η=cβAB.L_B=B+ Cdκ , A_B= dL_B+L_B, η= cβ A_B. Then the flux algorithm returns an ε -TV sample using Nqry≤CκAB1/2log(eκAB)(B+log4ε)log(eκLB)EN_qry≤ C κA_B^1/2 (eκ A_B) (B+ 4 ) (eκ L_B) expected oracle queries. Ignoring logarithmic factors, this is O~(κBdB+B1/2) O( κB\ dB+B\^1/2), and hence O~(κd1/4B5/4) O( κd^1/4B^5/4) when 1≤B≤d1≤ B≤ d. Proof. Apply Prop. 8 with Δ=B =B. The displayed choice gives a=c/(κAB)a=c/(κ A_B), s=c/ABs=c/A_B, and ℓ≤CLB ≤ CL_B, so the finite gate holds. Substitution into 21 gives the stated bound. ∎ B.4.2 Regularity and Domain Extensions For completeness, we collect here the regularity steps used above. The weak derivatives and integrations by parts are first justified for smooth potentials and test functions. Mollifying V preserves strong convexity and smoothness; the normalizers and conditional moments have a common Gaussian envelope, so the identities pass to the limit. The derivative estimate in Lem. 3 also begins on smooth functions. Closedness of the weak gradient extends ∇ T by continuity. On the spectral subspace [1/2,1][1/2,1], −1 T^-1 is bounded. Applying the extended estimate to −1f T^-1f and −1(I−)f T^-1(I- T)f therefore gives approximations converging both in L2L^2 and in weak-gradient norm; in particular, no nonsmooth spectral projection is applied directly to a classical derivative. The same graph-norm approximation extends Lem. 4 from smooth functions to its stated weak-gradient domain. The event-driven process itself needs no smoothing: its residual force is globally Lipschitz, its integrated hazard is continuous, and the direction h/∥h∥h/ h is evaluated only at positive event rate. The proof is complete relative to the earlier RGO result in Prop. 9 and standard strongly log-concave concentration and optimization facts. Appendix C Proofs for the Persistent Quarter-Turn This appendix proves Props. 3, 13, 3, 4 and 3. The argument is organized in the dependency order of the proof: 1. show that every fixed-reference conditional motion is exact; 2. quantify how outer reflection creates a conditional-position signal; 3. prove that a quarter-turn transfers this signal to and from momentum; 4. combine this transfer with occasional refreshment to obtain L2L^2 contraction; 5. replace the ideal reference, event clock, and refresh by finite counterparts and account for their errors; 6. absorb the known strong-convexity quadratic into the conditional motion to obtain the near-Gaussian refinement. Standard ingredients are cited rather than reproved. In particular, we use established skew-balance results for velocity-jump processes, standard Poincaré and log-Sobolev consequences of strong log-concavity, and the existing first-order RGO implementation. The new arguments are the position–momentum transfer estimate, its unit-circle contraction analysis, and the shifted-reference implementation. Throughout, π denotes the target. Uppercase X,Y,P,EX,Y,P,E denote position, auxiliary position, momentum, and reference random variables, with E corresponding to X~ X in the main text; lowercase letters denote their values. Operators act on observables, so laws evolve through their adjoints. Concretely, for a Markov transition A, (f)(z)=[f(Znext)∣Zcurrent=z]( Af)(z)=E[f(Z_ next) Z_ current=z]; thus operator products act from right to left on functions, even though the corresponding state updates occur in chronological order. The small set of operators used throughout is: symbol meaning P conditional expectation given Y =I− Q=I- P remaining conditional variation U outer reflection R momentum reversal, (x,y,p)↦(x,y,−p)(x,y,p) (x,y,-p) S conditional quarter-turn K one complete ideal iteration = T= P U P reflection seen on Y-only observables Γ=(I−2)1/2 =(I- T^2)^1/2 amount transferred out of the Y-only part C.1 Exactness of the Persistent Conditional Move We first isolate the only facts about the event-driven dynamics that are used later: it is well defined, it preserves the desired conditional law, and reversing momentum gives its stationary time reversal. These are standard skew-balance properties of bouncy and Boomerang processes (Bouchard-Côté et al., 2018; Vanetti et al., 2017; Bierkens et al., 2020; Monmarché et al., 2023). We verify that our conditional process satisfies the required identities. Fix Y=yY=y, use the normalized position w=(x−y)/ηw=(x-y)/ η, and write qyw(dw)∝exp−V(y+ηw)−∥w∥2/2dw,by(w)=η∇V(y+ηw).q_y^w(dw) \-V(y+ ηw)- w ^2/2\\,dw, b_y(w)= η\,∇ V(y+ ηw). For a nonzero vector h, let Refhp=p−2(pTh)h/∥h∥2Ref_hp=p-2(p^Th)h/ h ^2, and set Ref0p=pRef_0p=p. A physical reference E corresponds to e=(E−y)/ηe=(E-y)/ η. Conditional on (y,e)(y,e), the generator is (ℒy,eF)(w,p)= (L_y,eF)(w,p)= pT∇wF−w+by(e)T∇pF p^T _wF-\w+b_y(e)\^T _pF +[pTby(w)−by(e)]+F(w,Refby(w)−by(e)p)−F(w,p). +[p^T\b_y(w)-b_y(e)\]_+\F(w,Ref_b_y(w)-b_y(e)p)-F(w,p)\. (23) The first line is a deterministic harmonic flight centered at −by(e)-b_y(e). The last line says that a state-dependent Poisson clock rings at rate [pTby(w)−by(e)]+[p^T\b_y(w)-b_y(e)\]_+; when it rings, only the momentum is reflected across the residual-gradient direction by(w)−by(e)b_y(w)-b_y(e). Fixing e is a proof device. The accelerated branch later draws a fresh reference independently given y, and averages these fixed-reference kernels. Lemma 8 (Fixed-reference validity). For every fixed y,ey,e, 23 defines a non-explosive Markov process. If Sty,eS_t^y,e denotes its transition operator and (w,p)=(w,−p) R(w,p)=(w,-p), then qyw⊗γd is invariant,(Sty,e)∗=Sty,e.q_y^w _d is invariant, (S_t^y,e)^*= RS_t^y,e R. Both conclusions remain valid after averaging Sty,eS_t^y,e over any measurable reference law depending on y, provided the reference is conditionally independent of the incoming (W,P)(W,P). Proof. Put a=by(e)a=b_y(e), r(w)=by(w)−ar(w)=b_y(w)-a, and (w,p)=(w,Refr(w)p),s(w,p)=pTr(w). J(w,p)=(w,Ref_r(w)p), s(w,p)=p^Tr(w). Between events the process follows the harmonic flow ϕu(w,p)=(−a+(w+a)cosu+psinu,−(w+a)sinu+pcosu). _u(w,p)= (-a+(w+a) u+p u,-(w+a) u+p u ). Both this flow and J preserve H(w,p)=∥w+a∥2+∥p∥2/2H(w,p)=\ w+a ^2+ p ^2\/2. Thus every set KR=H≤RK_R=\H≤ R\ is invariant. Lipschitz continuity of byb_y bounds the event rate λ=[s]+λ=[s]_+ on KRK_R by ∥p∥∥by(w)−by(e)∥≤βη2R2R+∥e+a∥=:ΛR. p \, b_y(w)-b_y(e) ≤βη 2R\ 2R+ e+a \=: _R. Thinning a homogeneous Poisson process of rate ΛR _R therefore constructs the process from every point of KRK_R. Since every initial state lies in some KRK_R, this proves well-posedness and non-explosion. We prove skew balance directly on paths. The idea is to pair every path with its time reversal after flipping momentum: collision factors agree one by one, while the change in target density exactly cancels the change in the survival probability. Write z=(w,p)z=(w,p), set Ψ(w)=V(y+ηw)−aTw,ν=qyw⊗γd, (w)=V(y+ ηw)-a^Tw, ν=q_y^w _d, so that the Lebesgue density of ν is, up to normalization, ϱ(z)=exp−H(z)−Ψ(w) (z)= \-H(z)- (w)\. Along a flight, duΨ(wu)=s(zu) ddu (w_u)=s(z_u), while a collision changes neither H nor Ψ . Moreover, ϕu=ϕ−u,=,s∘=s∘=−s. R _u R= _-u, J R= R J, s J=s R=-s. (24) Fix t>0t>0. For =(t1,…,tn) t=(t_1,…,t_n) in the ordered simplex Δn(t)=0<t1<⋯<tn<t _n(t)=\0<t_1<·s<t_n<t\, let zu=zu(z0,)z_u=z_u(z_0, t) be the deterministic skeleton obtained by following ϕφ and inserting a collision J at every tit_i. On KRK_R, the iterated-hazard expansion of the transition operator is Sty,eF(z0)=∑n≥0∫Δn(t)exp−∫0tλ(zu)du∏i=1nλ(zti−)F(zt)d.S_t^y,eF(z_0)= _n≥ 0 _ _n(t) \! \- _0^tλ(z_u)\,du \ _i=1^nλ(z_t_i-)F(z_t)\,d t. (25) The bound by ΛR _R makes this series absolutely convergent. Pair this skeleton with the one starting at zt Rz_t and having ordered collision times (t−tn,…,t−t1)(t-t_n,…,t-t_1). By 24, it is the momentum-reversed time reversal of the original skeleton and ends at z0 Rz_0. The map (z0,t1,…,tn)⟼(zt,t−tn,…,t−t1)(z_0,t_1,…,t_n) ( Rz_t,t-t_n,…,t-t_1) preserves Lebesgue measure. Indeed, ϕu _u and R preserve volume, while J does so by applying, for every fixed w, an orthogonal change of variables in p. This argument does not require continuity of the reflection direction at r(w)=0r(w)=0, where the event rate vanishes. If z−z^- is a forward pre-collision state, the corresponding reverse pre-collision state is z− R Jz^-; hence λ(z−)=λ(z−).λ( R Jz^-)=λ(z^-). Thus the collision-rate products in 25 agree. The reverse survival hazard is ∫0t[−s(zu)]+u _0^t[-s(z_u)]_+\,du, and Ψ(wt)−Ψ(w0)=∫0ts(zu)u=∫0t[s(zu)]+−[−s(zu)]+u. (w_t)- (w_0)= _0^ts(z_u)\,du= _0^t \[s(z_u)]_+-[-s(z_u)]_+ \\,du. Since H(zt)=H(z0)H(z_t)=H(z_0) and ϱ(z)=ϱ(z) ( Rz)= (z), this is exactly the survival-weight identity ϱ(z0)exp−∫0t[s(zu)]+du=ϱ(zt)exp−∫0t[−s(zu)]+du. (z_0) \! \- _0^t[s(z_u)]_+\,du \= ( Rz_t) \! \- _0^t[-s(z_u)]_+\,du \. (26) Let νR _R be ν conditioned on KRK_R. Pairing paths term by term in 25 using 26 gives, for bounded measurable F,GF,G, ∫FSty,eGdνR=∫(G)Sty,e(F)dνR=∫(Sty,eF)GdνR. FS_t^y,eG\,d _R= ( RG)S_t^y,e( RF)\,d _R= ( RS_t^y,e RF)G\,d _R. Letting R→∞R→∞ gives the same identity under ν. Taking F=1F=1 proves invariance. Jensen’s inequality then makes Sty,eS_t^y,e a contraction on L2(ν)L^2(ν), and density of bounded functions extends the identity to L2(ν)L^2(ν), proving (Sty,e)∗=Sty,e(S_t^y,e)^*= RS_t^y,e R. Finally, the recursive hazard construction is jointly Borel in (y,e,z0)(y,e,z_0): the flow, rate, integrated hazard, and reflection are Borel, and successive event times are obtained by taking level-hitting times of the integrated hazard. Non-explosion gives the measurable pointwise limit. Consequently, conditioning on an incoming-state-independent reference and applying Fubini passes both invariance and the adjoint identity to every measurable reference mixture. ∎ We now pass from the conditional move to one complete proximal iteration. All operators below act on L2(Π^η)L^2( _η), where Π^η=Πη⊗γd _η= _η _d. Let U be the pullback of the outer reflection (x,y,p)↦(x,2x−y,p)(x,y,p) (x,2x-y,p), let f=[f∣Y] Pf=E[f Y], and set =I− Q=I- P. Finally, for the ideal chain let S be the quarter-period conditional kernel after drawing E∼qyE q_y, independently of the incoming state. Then == S P= P S= P. The conditional mixture and the observable operator of one ideal iteration are ρ=ρ+(1−ρ)=+(1−ρ),=ρ. B_ρ=ρ P+(1-ρ) S= P+(1-ρ) S Q, K= U B_ρ. (27) The first term is the independent conditional refresh: as an operator on observables, an exact draw from qy⊗γdq_y _d is precisely P. The second term is the persistent quarter-turn. Lemma 9 (Invariance and refresh dissipation). The operator K preserves Π^η _η and is a contraction. Moreover, for every f∈L2(Π^η)f∈ L^2( _η), ∥f∥2−∥f∥2≥ρ∥f∥2. f ^2- Kf ^2≥ρ Qf ^2. (28) Proof. The outer reflection is a measure-preserving involution because it sends Y−XY-X to its negative. Lemma 8 shows that the conditional mixture preserves the conditional position–momentum law; hence their composition preserves Π^η _η. Orthogonality of f Pf and f Qf, together with contractivity of S, gives ∥f∥2−∥f∥2=∥f∥2−(1−ρ)2∥f∥2≥1−(1−ρ)2∥f∥2, f ^2- Kf ^2= Qf ^2-(1-ρ)^2 S Qf ^2≥\1-(1-ρ)^2\ Qf ^2, which is stronger than 28. ∎ C.2 How Reflection Feeds the Quarter-Turn This subsection proves the sampling-specific estimate used by the contraction argument. Reflection first turns a function of Y into an almost affine function of the conditional position. The quarter-turn then sends that affine position component into momentum and sends the corresponding momentum component back with the opposite sign. At a semantic level, the signal flow is f(Y+)→f(Y−) f(Y_+) \ U\ f(Y_-) =f+f, = Tf+ Cf, :=, C:= Q U P, f Cf ≈Lf, ≈ L_f, f S Cf ≈Mf, ≈ M_f, Mf SM_f ≈−f, ≈- Cf, Mf UM_f ≈Mf. ≈ M_f. The proposition below controls the last three errors simultaneously; the intermediate affine observable LfL_f makes that control possible. Refreshment then damps whatever is not represented by these signals. Notation and the result to be proved. Let f=[f∣Y] Pf=E[f Y], =I− Q=I- P, and let U denote outer reflection. On the centered Y-only subspace ℋY=ran()∩L02H_Y=ran( P)∩ L^2_0, define =,Γ=(I−2)1/2,g=αη. T= P U P, =(I- T^2)^1/2, g= αη. Under the stationary reflection coupling, write Y+=Y_+=Y, Y−=2X−Y_-=2X-Y, and W=(X−Y+)/ηW=(X-Y_+)/ η. Conditional on Y+=yY_+=y, denote the mean and covariance of W by w¯y w_y and Σy _y, and set Dy=I−ΣyD_y=I- _y. For a differentiable f∈ℋYf _Y, define af(y)=[∇f(Y−)∣Y+=y],Lf=2ηaf(Y+)T(W−w¯Y+),Mf=2ηaf(Y+)TP.a_f(y)=E[∇ f(Y_-) Y_+=y], L_f=2 η\,a_f(Y_+)^T(W- w_Y_+), M_f=2 η\,a_f(Y_+)^TP. The letter LfL_f denotes the affine position observable; MfM_f is its momentum counterpart. Componentwise, the definition also says af=∇fa_f= T∇ f, where T acts coordinatewise on vector-valued functions. The roles of the main quantities are summarized as follows. quantity role in the transfer argument f Tf the part of the reflected observable predictable from Y+Y_+ alone Γf f the scale of the conditional-position signal created by reflection f Cf the conditional-position signal itself LfL_f the leading affine approximation to that position signal MfM_f the momentum copy of LfL_f produced by a quarter-turn These formulas are initially classical for differentiable f. The proposition below is an all–L2L^2 statement: the closed weak-gradient lemma that follows proves that every f∈ℋ+f _+ has the required weak gradient and defines af,Mfa_f,M_f by graph closure. The auxiliary observable LfL_f is needed only in the smooth-core calculation. In particular, the proposition does not assume that a spectral projection preserves smoothness. Proposition 10 (Position–momentum transfer). Let ℋ+H_+ be the spectral subspace of T corresponding to [1/2,1][1/2,1]. If βηβη is smaller than a universal constant, then every f∈ℋ+f _+ satisfies ∥(f)−Mf∥2+∥Mf+f∥2+∥Mf−Mf∥2 S( Q U Pf)-M_f ^2+ SM_f+ Q U Pf ^2+ UM_f-M_f ^2 ≤C∥Γf∥∥Γ3f∥. 99.58464pt≤ C f ^3f . (29) The proof has a technical prelude followed by three substantive steps. The prelude extends estimates from smooth observables to the closed weak-gradient domain. Step 1 analyzes the signal created by reflection, Step 2 analyzes its quarter-rotation, and Step 3 combines the two estimates and performs the single closure passage. Technical prelude: extending smooth estimates to all L2L^2. We first isolate the only closure argument needed below. A graph core here means that every function with a weak gradient can be approximated together with its gradient by smooth compactly supported functions. Let ν=ΠηYν= _η^Y, let D:Dom(D)⊂L02(ν)→L2(ν,ℝd)D:Dom(D)⊂ L^2_0(ν)→ L^2(ν;R^d) be the closed weak gradient, and put 0=Cc∞(ℝd)∩L02(ν).C_0=C_c^∞(R^d)∩ L^2_0(ν). Lemma 10 (Closed weak-gradient smoothing). The set 0C_0 is a graph core for D. Suppose that T is a self-adjoint contraction, 0⊂Dom(D) TC_0 (D), and, on 0C_0, η∥Du∥≤cD∥Γu∥. η D Tu ≤ c_D u . (30) Then DD T has an everywhere-defined bounded extension D:L02(ν)→L2(ν,ℝd) R_D:L^2_0(ν)→ L^2(ν;R^d), and in fact L02(ν)⊂Dom(D),D=D,η∥Du∥≤cD∥Γu∥. TL^2_0(ν) (D), D T= R_D, η R_Du ≤ c_D u . Consequently, if Π+=[1/2,1]() _+=1_[1/2,1]( T), then ran(Π+)⊂Dom(D)ran( _+) (D) and η∥Df∥ η Df ≤2cD∥Γf∥, ≤ 2c_D f , η∥D(I−)f∥ η D(I- T)f ≤2cD∥Γ3f∥,f∈ran(Π+). ≤ 2c_D ^3f , f ( _+). Moreover, for every f∈Dom(D)f (D), D(I−)f=Df−Df.D(I- T)f=Df- R_Df. (31) Thus convergence in the graph norm of D automatically implies convergence of D(I−)fD(I- T)f; no additional graph norm is needed. Proof. The density of ν is positive and continuous, hence is bounded above and below on every compact set. Cutoff followed by ordinary convolution therefore approximates every element of Dom(D)Dom(D) in the norm ∥f∥2+∥Df∥2 f ^2+ Df ^2. To preserve centering, subtract the mean of each approximant times a fixed ψ∈Cc∞ψ∈ C_c^∞ satisfying ∫ψν=1 ψ\,dν=1. This proves the graph-core assertion. If un∈0u_n _0 converges in L2(ν)L^2(ν), 30 makes DunD Tu_n Cauchy, while un Tu_n converges. Closedness of D shows that the limit belongs to Dom(D)Dom(D), identifies its gradient, and gives the bounded extension D=D R_D=D T. On ran(Π+)ran( _+), −1 T^-1 is bounded with norm at most two. For f in that subspace, Df=D−1f,D(I−)f=D−1(I−)f.Df= R_D T^-1f, D(I- T)f= R_D T^-1(I- T)f. The first estimate follows from γt/t≤2γt _t/t≤ 2 _t on t∈[1/2,1]t∈[1/2,1], where γt=(1−t2)1/2 _t=(1-t^2)^1/2. The second follows from 1−t=1−t2t(1+t)≤2(1−t2),t∈[1/2,1]. 1-tt= 1-t^2t(1+t)≤ 2(1-t^2), t∈[1/2,1]. Finally, f∈Dom(D) Tf (D) for every f∈L02(ν)f∈ L^2_0(ν), so 31 follows by linearity and proves the last assertion. ∎ C.2.1 Step 1: What the Outer Reflection Creates Lemma 11 (Outer-reflection geometry). Assume βη≤1/3βη≤ 1/3. On ℋYH_Y, −12I⪯1−g21+g2I,Γ⪰2g1+g2I.- 12I T 1-g^21+g^2I, 2g1+g^2I. Moreover, := C:= Q U P can be written as =Γ C= V , where V is an isometry, and f=f+Γf,f=Γf−f. U Pf= P Tf+ V f, U Vf= P f- V Tf. On ℋ+H_+, η∥∇f∥≤C∥Γf∥,η∥∇(I−)f∥≤C∥Γ3f∥. η ∇ f ≤ C f , η ∇(I- T)f ≤ C ^3f . Proof. The pair (Y+,Y−)(Y_+,Y_-) is exchangeable and (f)(y)=[f(Y−)∣Y+=y]( Tf)(y)=E[f(Y_-) Y_+=y], so T is a self-adjoint contraction. The conditional density of Y−Y_- given Y+=yY_+=y is proportional to exp−V(y+u2)−∥y−u∥28η. \-V ( y+u2 )- y-u ^28η \. Its score with respect to y is sy(u)=−∇V((y+u)/2)/2−(y−u)/(4η)s_y(u)=-∇ V((y+u)/2)/2-(y-u)/(4η). For f∈0f _0, weak differentiation under the integral gives Df(y)=Cov(f(Y−),sy(Y−)∣Y+=y).D Tf(y)=Cov(f(Y_-),s_y(Y_-) Y_+=y). Indeed, on compact sets of y, the score times the conditional density is dominated by an integrable Gaussian tail. Thus Fubini and weak integration by parts justify differentiation of both the numerator and normalizing constant. In particular, 0⊂Dom(D) TC_0 (D). Applying conditional Poincaré to this covariance gives η∥Df∥2≤(1−αη)24(1+αη)∥f∥2−∥f∥2≤14∥Γf∥2.η D Tf ^2≤ (1-αη)^24(1+αη)\ f ^2- Tf ^2\≤ 14 f ^2. The Y-marginal has Poincaré constant at most α−1+ηα^-1+η: this follows from the conditional covariance bound Cov(X∣Y)⪯η(1+αη)−1ICov(X Y) η(1+αη)^-1I and Brascamp–Lieb (Chewi, 2023). Equivalently, if WηW_η is the negative log-density of the Y-marginal, then ∇2Wη(y)=η−1I−η−2Cov(X∣Y=y)⪰α1+αηI.∇^2W_η(y)=η^-1I-η^-2Cov(X Y=y) α1+αηI. Applying the resulting Poincaré inequality to the centered function f Tf yields ∥L02≤(1−g2)/(1+g2) T _L^2_0≤(1-g^2)/(1+g^2). This inequality is first obtained on 0C_0 and extends to L02(ν)L^2_0(ν) by density. For the lower edge, exchangeability and Gaussian Poincaré imply ⟨h,(I−)h⟩=12h(Y+)−h(Y−)2≤2η∥∇h∥2. h,(I- T)h = 12E\h(Y_+)-h(Y_-)\^2≤ 2η ∇ h ^2. Take h=fh= Tf for f∈0f _0 and use the preceding derivative bound. The resulting inequality involves only bounded polynomials in T, so density extends it to all of L02(ν)L^2_0(ν). Spectral calculus gives 22(I−)⪯I−22 T^2(I- T) I- T^2, equivalently (I−)2(I+2)⪰0(I- T)^2(I+2 T) 0, and hence ⪰−I/2 T -I/2. The two spectral edges imply the lower bound on Γ . Since 1=1 T1=1, the lower bound ⪰−I/2 T -I/2 also holds on the full scalar L2(ν)L^2(ν), and hence coordinatewise on vector fields such as DfDf, whose components need not be centered. Finally, ∗=I−2=Γ2 C^* C=I- T^2= ^2. Thus =Γ−1 V= C ^-1 is an isometry, and 2=I U^2=I gives the two block identities. The core derivative estimate above is 30 with cD=1/2c_D=1/2. Applying Lem. 10 proves that all of ℋ+H_+, not merely a smooth subset of it, belongs to the weak gradient domain and gives the two final inequalities. Notice that this argument never applies the discontinuous spectral projection Π+ _+ to a smooth function. ∎ Lemma 12 (Affine approximation after reflection). For every f∈0f _0, f=Lf+ef Q U Pf=L_f+e_f, where [ef∣Y+]=0E[e_f Y_+]=0. With ℛf=∥∇f∥2−∥∇f∥2,f=[af(Y+)TDY+af(Y+)],R_f= ∇ f ^2- T∇ f ^2, _f=E[a_f(Y_+)^TD_Y_+a_f(Y_+)], one has ∥ef∥2≤4ηℛf,ℛf+f≤C∥∇(I−)f∥∇f∥. e_f ^2≤ 4 _f, _f+D_f≤ C ∇(I- T)f ∇ f . (32) Proof. Conditionally on Y+=yY_+=y, we have Y−=y+2ηWY_-=y+2 ηW. Set ef=f−Lfe_f= Q U Pf-L_f. The choice of af(y)a_f(y) makes both [ef∣Y+=y]E[e_f Y_+=y] and [∇Wef∣Y+=y]E[ _We_f Y_+=y] zero. The conditional position law is one-strongly log-concave in W, so conditional Poincaré gives the first inequality. For the second, let by(w)=η∇V(y+ηw)b_y(w)= η\,∇ V(y+ ηw). Weak integration by parts under the conditional density gives Cov(W,by(W)∣y)=Dy,Cov(by(W)∣y)⪯βη1−βηDy.Cov(W,b_y(W) y)=D_y, (b_y(W) y) βη1-βηD_y. Conditional Poincaré applied to linear functions gives Σy⪯I _y I, while the Stein identity [∇by]=Dy+Cov(by)E[∇ b_y]=D_y+Cov(b_y) and 0⪯∇by⪯βηI0 ∇ b_y βη I give 0⪯Dy⪯βηI0 D_y βη I. The second relation follows by applying conditional Poincaré to byb_y and using (∇by)2⪯(βη)∇by(∇ b_y)^2 (βη)∇ b_y. These are standard Brascamp–Lieb/Stein calculations (Chewi, 2023); for a C1,1C^1,1 potential they follow by mollification and cutoff. Differentiating f Tf weakly and inserting f=Lf+ef Q U Pf=L_f+e_f yields ∇(I−)f=∇f−af+2Dyaf+η−1/2Cov(ef,by(W)∣Y+=y).∇(I- T)f=∇ f-a_f+2D_ya_f+η^-1/2Cov(e_f,b_y(W) Y_+=y). Take its inner product with ∇f+af∇ f+a_f and average over y. Using af=∇fa_f= T∇ f, exact rearrangement gives ℛf+4f= _f+4D_f= ⟨∇(I−)f,∇f+af⟩−2[afTDY+(∇f−af)] ∇(I- T)f,∇ f+a_f -2E\! [a_f^TD_Y_+(∇ f-a_f) ] −η−1/2[Cov(ef,bY+(W)∣Y+)T(∇f+af)]. -η^-1/2E\! [Cov(e_f,b_Y_+(W) Y_+)^T(∇ f+a_f) ]. Cauchy–Schwarz, the two conditional covariance bounds above, and the first inequality control the right-hand side by 2∥∇(I−)f∥∇f∥+Cβηℛff+Cβηℛf(βηℛf+f).2 ∇(I- T)f ∇ f +C βη\,R_fD_f+C βη\,R_f(βη\,R_f+D_f). For sufficiently small βηβη, Young’s inequality absorbs the last two terms and proves 32. ∎ C.2.2 Step 2: What the Quarter-Turn does to Affine Information At fixed y, write ℓa,y(w)=aT(w−w¯y) _a,y(w)=a^T(w- w_y) and ma(p)=aTpm_a(p)=a^Tp. Let S denote the quarter-period kernel after drawing an independent exact reference from the conditional position law. All norms and expectations in the next lemma are under qyw⊗γdq_y^w _d, including the independent reference randomness. Lemma 13 (Approximate quarter rotation). If Lη:=βη≤1/2L_η:=βη≤ 1/2, then, uniformly in y, ∥ℓa,y−ma∥2+∥ma+ℓa,y∥2≤22Lη1−LηaTDya. S _a,y-m_a ^2+ Sm_a+ _a,y ^2≤ 22L_η1-L_ηa^TD_ya. Proof. For fixed reference e, the generator in 23 maps ℓa,y _a,y to mam_a and maps mam_a to −ℓa,y+ka,e- _a,y+k_a,e, where, with r=by(w)−by(e)r=b_y(w)-b_y(e) and Z=pTr/∥r∥Z=p^Tr/ r , ka,e=[aTby(W)∣y]−aTby(e)−2aTrZ2Z>0.k_a,e=E[a^Tb_y(W) y]-a^Tb_y(e)-2a^Tr\,Z^21_\Z>0\. When r=0r=0, only the final term is defined to be zero; the first two terms are unchanged. The localized Dynkin formula therefore gives, at T=π/2T=π/2, STy,eℓa,y−ma=∫0TSsy,eka,ecosss,STy,ema+ℓa,y=∫0TSsy,eka,esinss.S_T^y,e _a,y-m_a= _0^TS_s^y,ek_a,e s\,ds, S_T^y,em_a+ _a,y= _0^TS_s^y,ek_a,e s\,ds. This is the variation-of-constants identity for the two observable equations (equivalently, a localized Dynkin formula). Localization is legitimate because the harmonic energy in the proof of Lem. 8 is conserved; dominated convergence then removes the cutoff. When e is an independent conditional draw, elementary Gaussian moments give ka,e2=11Var(aTby(W)∣y)Ek_a,e^2=11Var(a^Tb_y(W) y). The covariance calculation in the proof of Lem. 12 bounds this variance by Lη(1−Lη)−1aTDyaL_η(1-L_η)^-1a^TD_ya. Contractivity of Ssy,eS_s^y,e and Cauchy–Schwarz in time prove the claimed bound. ∎ C.2.3 Step 3: Assemble the Global Transfer Estimate The single graph closure. For every centered f∈Dom(D)f (D), define af=Df,Mf=2ηaf(Y)TP,=.a_f= TDf, M_f=2 η\,a_f(Y)^TP, C= Q U P. This agrees with the preceding definition on smooth functions. Since P is conditionally standard Gaussian and T is a contraction on vector-valued L2L^2, for f,h∈Dom(D)f,h (D), ∥Mf−Mh∥2=4η∥D(f−h)∥2≤4η∥D(f−h)∥2. M_f-M_h ^2=4η TD(f-h) ^2≤ 4η D(f-h) ^2. (33) Consider the stacked error map from Dom(D)Dom(D), equipped with the graph norm ∥f∥D2=∥f∥2+η∥Df∥2 f _D^2= f ^2+η Df ^2, into three copies of the stationary L2L^2 space: f=(f−Mf,Mf+f,Mf−Mf). Ef= ( S Cf-M_f, SM_f+ Cf, UM_f-M_f ). The operators ,, C, S, U are contractions, so 33 makes E graph-norm continuous. Moreover, 31 shows that f↦D(I−)f D(I- T)f is graph-norm continuous. These two facts will pass the final three-error inequality through the core in one step; none of its intermediate conditional error quantities needs a separate closure argument. Proof of Prop. 10. First let f∈0f _0. Apply Lem. 13 conditionally with coefficient 2ηaf(y)2 η\,a_f(y), then add the remainder efe_f from Lem. 12. The first two terms on the left of 29 are bounded by Cβη2f+∥ef∥2C\βη^2D_f+ e_f ^2\. For the third term, reflection exchanges Y+Y_+ and Y−Y_- without changing momentum. Since ⪰−I/2 T -I/2, ∥Mf−Mf∥2≤4ηℛf UM_f-M_f ^2≤ 4 _f. Equation 32 now bounds the sum of all three terms by Cη∥∇(I−)f∥∇f∥.Cη ∇(I- T)f ∇ f . Equivalently, on 0C_0 we have the pre-spectral estimate ∥f∥2≤Cη∥D(I−)f∥∥Df∥. Ef ^2≤ Cη D(I- T)f Df . (34) Now take an arbitrary f∈ℋ+f _+. By Lem. 10, f∈Dom(D)f (D). Because 0C_0 is a graph core, there are fn∈0f_n _0 with fn→f_n→ f and Dfn→DfDf_n→ Df. The approximants are not required to lie in ℋ+H_+. The graph continuity proved above gives fn→f,D(I−)fn→D(I−)f. Ef_n→ Ef, D(I- T)f_n→ D(I- T)f. Passing to the limit in 34 and only then applying the two sector derivative estimates in Lem. 11 proves 29. Thus no density assertion about 0∩ℋ+C_0 _+, and no smoothness preservation by the spectral projection, is used. ∎ A bounded form used below. The transfer proposition is now complete. We next normalize its position and momentum signals to the same scale so that the contraction proof can treat them as two coordinates of one approximate rotation. The isometry V is the normalized position channel, A will be the normalized momentum channel, and F1,F2,F3F_1,F_2,F_3 record the three departures from an exact quarter-rotation. The later factorization Fj=BjΓF_j=B_j says that these errors shrink with the amount of signal being transferred. On ℋ+H_+, the definition above and Lem. 11 make f↦Mff M_f a bounded operator M: indeed ∥Mf∥2≤4η∥∇f∥2≤C∥Γf∥2 M_f ^2≤ 4η ∇ f ^2≤ C f ^2. Put =Γ−1 A= M ^-1, and recall =Γ Q U P= V . A is therefore bounded on ℋ+H_+. Define F1,F2,F3F_1,F_2,F_3 there by =+F1,=−+F2,=+F3. S V= A+F_1, S A=- V+F_2, U A= A+F_3. (35) Position observables are even in momentum and momentum observables are odd, so the ranges of V and A are orthogonal. Momentum reversal in Lem. 8 gives the adjoint versions of these relations. When an adjoint is taken in the next subsection, ,,Fj V, A,F_j mean their restrictions to ℋ+H_+, extended by zero on its orthogonal complement. Thus all subsequent adjoints are bounded, everywhere-defined Hilbert-space operators. Corollary 2 (Controlled transfer errors). Let ℓg=1+log(1/g) _g=1+ (1/g). There are bounded operators BjB_j such that Fj=BjΓ,Fj∗=ΓBj∗,∥Bj∥≤Cℓg,j=1,2,3.F_j=B_j , F_j^*= B_j^*, B_j ≤ C _g, j=1,2,3. (36) Proof. Equation 29, with f replaced by Γ−1f ^-1f, gives ∥Fjf∥2≤C∥f∥∥Γ2f∥ F_jf ^2≤ C f ^2f . Partition the spectrum of Γ2 ^2 into dyadic intervals between cg2cg^2 and 11. On each spectral component frf_r, ∥Fjfr∥2≤C⟨fr,Γ2fr⟩ F_jf_r ^2≤ C f_r, ^2f_r . There are O(ℓg)O( _g) components, so the triangle inequality followed by Cauchy–Schwarz gives Fj∗Fj⪯CℓgΓ2.F_j^*F_j C _g ^2. Since Γ is invertible on the centered subspace, setting Bj=FjΓ−1B_j=F_j ^-1 proves 36. The factor order is important; no commutation between BjB_j and Γ is used. ∎ C.3 From Position–Momentum Transfer to Contraction We now prove the ideal-chain bound. Reflection alone preserves norm, and the quarter-turn mainly rotates information between position and momentum, so a one-step Poincaré estimate cannot reveal the contraction. Instead, we rule out oscillations on the unit circle and then convert that bound into geometric decay. The resolvent strategy asks whether (z−)f=b(z- K)f=b can have a large solution when |z|=1 z =1. For each slowly varying Y-mode, the maps V and A expose its position and momentum companions; the resulting three coordinates form an approximate finite-dimensional rotation. The refresh dissipation estimate controls the orthogonal conditional variation that these coordinates do not see. A uniform bound on the unit-circle resolvent then yields geometric decay through Lem. 15. Proposition 11 (Ideal persistent contraction). Under the assumptions of Thm. 3, choosing ρ≍gℓgρ g _g, where g=αηg= αη and ℓg=1+log(1/g) _g=1+ (1/g), gives ∥n∥L02→L02≤2exp(−cngℓg). K^n _L^2_0→ L^2_0≤ 2 (- cng _g ). (37) The proof has three steps. We follow three slowly changing coordinates, invert their three-by-three update on the unit circle, and use refreshment to control all remaining conditional variation. The coordinates do not form a decomposition of the full state space; this is why the final refreshment step is essential. We allow complex-valued observables below so that z may range over the complex unit circle; this does not change any operator norm. Let k=1−ρk=1-ρ, and let K be the one-step operator in 27. Recall ,Γ,, T, , V, A from Lems. 11 and 35. Lemma 14 (Unit-circle resolvent bound). If ρ≍gℓgρ g _g, then sup|z|=1∥(z−)−1∥L02→L02≤Cℓg. _ z =1 (z- K)^-1 _L^2_0→ L^2_0≤ C _gg. Proof. Step 1: isolate the three slow coordinates. Fix |z|=1 z =1 and solve (z−)f=b(z- K)f=b. Write f⟂=f_ = Qf. Let Πslow=[1/2,1]() _ slow=1_[1/2,1]( T), and split the Y-only component into xslow=Πslowfx_ slow= _ slow Pf and xfast=(I−Πslow)fx_ fast=(I- _ slow) Pf. In the three-coordinate calculation, every occurrence of ,Γ,,,Fj T, , V, A,F_j is restricted to ℋ+=ran(Πslow)H_+=ran( _ slow); equivalently, V means Πslow V _ slow, and similarly for the other maps. The other two tested coordinates are then v=∗f⟂v= V^*f_ and m=∗f⟂m= A^*f_ . coordinate information it measures xslowx_ slow the slowly changing component stored in the auxiliary variable Y v=∗f⟂v= V^*f_ conditional variation aligned with the position signal m=∗f⟂m= A^*f_ conditional variation aligned with the momentum signal xfastx_ fast the remaining Y-modes, controlled directly by the spectral gap of z−z- T The next identities encode the approximate quarter-rotation in these coordinates; their error terms are precisely the factorized FjF_j’s from the preceding subsection. The adjoint identities used below are ∗u=−∗u+F1∗u,∗u=∗u−F2∗u,∗=−∗,∗=∗+F3∗. V^* Su=- A^*u+F_1^* Ru, A^* Su= V^*u-F_2^* Ru, V^* U Q=- T V^*, A^* U= A^*+F_3^*. They follow from 35, momentum reversal, and the even/odd momentum parity recorded there. No xfastx_ fast term is omitted below: spectral orthogonality gives ∗xfast=0 V^* U Px_ fast=0, while momentum parity gives ∗xfast=0 A^* U Px_ fast=0 and F3∗xfast=0F_3^* Px_ fast=0. Projecting the resolvent equation onto these coordinates gives (z−)xslow+kΓm (z- T)x_ slow+k m =b0+kΓr1, =b_0+k r_1, zv−Γxslow−km zv- x_ slow-k Tm =bV−kr1, =b_V-k Tr_1, (38) zm−kv zm-kv =bA+kr2, =b_A+kr_2, where b0=Πslowb_0= _ slow Pb, bV=∗b_V= V^*b, bA=∗b_A= A^*b, and r1=F1∗f⟂,r2=−F2∗f⟂+F3∗f⟂.r_1=F_1^* Rf_ , r_2=-F_2^* Rf_ +F_3^* Sf_ . Corollary 2 implies ∥Γ−1r1∥+∥Γ−1r2∥≤Cℓg∥f⟂∥. ^-1r_1 + ^-1r_2 ≤ C _g f_ . This is the only logarithmic loss. The factorization Fj∗=ΓBj∗F_j^*= B_j^* is used in this order; no commutation is assumed. Step 2: invert the three-coordinate system. All coefficients on the left of 38 are functions of the single self-adjoint operator T. At a scalar spectral value t∈[1/2,1)t∈[1/2,1), put γt=1−t2 _t= 1-t^2. The determinant of the coefficient matrix is dt,k(z)=z3−tz2−k2tz+k2.d_t,k(z)=z^3-tz^2-k^2tz+k^2. Taking the first row of the adjugate gives the exact operator identity d,k(z)xslow= d_ T,k(z)x_ slow= (z2−k2)b0−k2ΓbV−kzΓbA (z^2-k^2 T)b_0-k^2 b_V-kz b_A +kz2Γr1−k2zΓr2. +kz^2 r_1-k^2z r_2. (39) The cancellation that produces the final two terms in this form is what prevents an extra factor g−1g^-1. We use two elementary scalar estimates. If |z−t|<γt/16 z-t < _t/16, then |dt,k(z)|≥1532γt2,|z2−k2t|≤4(ρ+γt). d_t,k(z) ≥ 1532 _t^2, z^2-k^2t ≤ 4(ρ+ _t). For the first, write z=e2iuz=e^2iu. The near condition gives |u|≤γt/16 u ≤ _t/16, while the real part of e−3iudt,k(e2iu)e^-3iud_t,k(e^2iu) is (1+k2)cos(3u)−tcosu(1+k^2)\ (3u)-t u\. The identity cos(3u)−tcosu=(1−t)cosu−2sin(2u)sinu (3u)-t u=(1-t) u-2 (2u) u, together with 1−t=γt2/(1+t)1-t= _t^2/(1+t), proves the lower bound. The second estimate follows from |z2−k2t|≤2|z−1|+(1−t)+(1−k2) z^2-k^2t ≤ 2 z-1 +(1-t)+(1-k^2). On the spectral region |z−t|≥γt/16 z-t ≥ _t/16, the first row of 38 directly controls xslowx_ slow. On its complement, apply the two scalar estimates to 39. The terms involving b cost at most C(ρ/g2+g−1)∥b∥C(ρ/g^2+g^-1) b , while 36 bounds the two error terms as follows. If E denotes this near spectral region, then ∥Ed,k(z)−1ΓFj∗q∥=∥Ed,k(z)−1Γ2Bj∗q∥≤Cℓg∥q∥. 1_Ed_ T,k(z)^-1 F_j^*q = 1_Ed_ T,k(z)^-1 ^2B_j^*q ≤ C _g q . Here the lower bound on dt,kd_t,k cancels the two powers of γt _t supplied by Γ2 ^2. Since Γ⪰cgI cgI and ρ≍gℓgρ g _g, both regions yield ∥f∥≤Cℓg(∥f⟂∥+∥b∥g). Pf ≤ C _g ( f_ + b g ). (40) Here the fast component has not been dropped: on [−1/2,1/2)[-1/2,1/2), z−z- T has inverse norm at most two, so ∥xfast∥≤2(∥b∥+∥f⟂∥) x_ fast ≤ 2( b + f_ ). Step 3: control everything not tested above. Equation 28 and f=zf−b Kf=zf-b imply ρ∥f⟂∥2≤2∥f∥∥b∥ρ f_ ^2≤ 2 f b . Combine this with 40, use f=f+f⟂f= Pf+f_ , and apply Young’s inequality. The result is ∥f∥≤Cℓg∥b∥. f ≤ C _gg b . Thus z−z- K is bounded below and has closed range. It is also onto. Otherwise some nonzero u would satisfy ∗u=z¯u K^*u= zu. Because K is a contraction and |z|=1 z =1, equality of norms implies ∗u=u K K^*u=u, and hence u=zu Ku=zu. This contradicts the lower bound at z, proving the asserted resolvent estimate. ∎ The next elementary Hilbert-space lemma converts the resolvent estimate into a bound on powers. We retain its short proof because the precise linear dependence on M is needed for the final rate. Lemma 15 (Unit-circle bound implies geometric decay). Let Q be a contraction on a Hilbert space. If sup|z|=1∥(z−Q)−1∥≤M _ z =1 (z-Q)^-1 ≤ M with M≥1M≥ 1, then ∥Qn∥≤2exp(−cn/M) Q^n ≤ 2 (-cn/M). Proof. For aj=sin((j+1)π/(N+2))a_j= ((j+1)π/(N+2)), apply the resolvent bound to ∑j=0Najz−j−1Qjf _j=0^Na_jz^-j-1Q^jf and integrate over the unit circle. Parseval’s identity gives ∑j=0Naj2∥Qjf∥2≤M2∑j=0N+1(aj−aj−1)2∥Qjf∥2,a−1=aN+1=0. _j=0^Na_j^2 Q^jf ^2≤ M^2 _j=0^N+1(a_j-a_j-1)^2 Q^jf ^2, a_-1=a_N+1=0. Contractivity and ∑(aj−aj−1)2=21−cos(π/(N+2))∑aj2Σ(a_j-a_j-1)^2=2\1- (π/(N+2))\Σ a_j^2 imply ∥QN∥≤πM/(N+2) Q^N ≤π M/(N+2). Taking N0=⌈2πM⌉N_0= 2π M gives ∥QN0∥≤1/2 Q^N_0 ≤ 1/2; writing n=qN0+rn=qN_0+r proves the claim. ∎ Proof of Prop. 11. Lemmas 9, 35 and 2 provide the invariance, dissipation, row identities, and controlled errors used in Lem. 14; these ultimately follow from Lems. 8 and 10. Apply Lem. 15 with M=Cℓg/gM=C _g/g. This proves 37. ∎ Proof of Thm. 3. Choose the universal constant in nblkn_blk sufficiently large that 37 gives ∥nblk∥L02→L022≤e−1 K^n_blk _L^2_0→ L^2_0^2 ^-1. If rν=dν/dΠ^ηr_ν=dν/d _η, then the density ratio evolves under ∗ K^* and χ2(ν∥Π^η)=∥rν−1∥22χ^2(ν\| _η)= r_ν-1 _2^2. Since an operator and its adjoint have the same norm, applying the block estimate J times proves 6. For the total-variation conclusion, initialize Y0|X0∼(X0,ηId)Y_0 X_0 (X_0,η I_d) and P0∼γdP_0 _d. The initial and stationary augmented laws have the same conditional law of (Y,P)(Y,P) given X, so augmentation preserves D2D_2. Choose the smallest admissible nblkn_blk and put J=⌊n/nblk⌋J= n/n_blk . The remaining iterations cannot increase χ2χ^2 because an invariant Markov kernel is an L2L^2 contraction. Combining 6 with TV≤12χ2TV≤ 12 χ^2, and enlarging the universal constant, gives 7; marginalization proves the claim for XnX_n. ∎ C.4 Implementable Algorithm and Warm-Start Accounting We begin with the result proved in this subsection and then verify its four implementation components: the shifted Gaussian reference, a deterministic iteration cap for the mode solver, a capped event intensity simulated by Poisson thinning, and the approximate RGO used on the RGO branch. Theorem 6 (Implementation theorem). Let D2(ν0∥π)≤ΔD_2( _0\|π)≤ , 0<ε≤1/20< ≤ 1/2, and define Λε=Δ+logCdκε,η=cβ(dΛε+Λε). _ = + Cdκ , η= cβ( d _ + _ ). Put g=αηg= αη, ℓg=1+log(1/g) _g=1+ (1/g), choose ρ≍gℓgρ g _g, and set N=⌈CℓgΔ+log(8/ε)⌉.N= C _gg\ + (8/ )\ . Initialize X0∼ν0X_0 _0, Y0|X0∼(X0,ηId)Y_0 X_0 (X_0,η I_d), and P0∼γdP_0 _d, independently conditional on X0X_0. There is an implementable algorithm for the persistent iteration in Sec. 5, using Sec. 5.1 as its accelerated branch, that returns XNX_N with TV(Law(XN),π)≤εTV(Law(X_N),π)≤ and Nqry≤O~[κ(Δ+log(8/ε))(dΛε+Λε)1/2],Nrgo≤O~(Δ+log(8/ε)).EN_qry≤ O\! [ κ\,( + (8/ ))( d _ + _ )^1/2 ], _rgo≤ O( + (8/ )). C.4.1 A Comparison Principle for the Complete Path The chain induced by the implementable algorithm is not invariant because it contains truncation and approximation. We therefore keep the exact chain as a comparison process and couple the two chains until their first disagreement. Lemma 16 (Warm path comparison). Let an exact Markov kernel preserve π, and let ν≪πν π satisfy D2(ν∥π)≤ΔD_2(ν\|π)≤ . If A is an event determined by a Markov path together with all its internal randomness, then Prν(A)≤eΔ/2Prπ(A)1/2. _ν(A)≤ e /2 _π(A)^1/2. Consequently, if the stationary one-step probability that an exact and an approximate transition disagree is at most a, their N-step laws started from ν can be coupled with disagreement probability at most NeΔ/2aNe /2 a. Proof. For R=dν/dπR=dν/dπ, conditional Cauchy–Schwarz gives π[RPr(A∣X0)]≤∥R∥2Prπ(A)1/2E_π[R (A X_0)]≤ R _2 _π(A)^1/2. An invariant exact kernel cannot increase D2D_2. Apply the same bound before each exact transition and take a union bound over the first disagreement. ∎ C.4.2 Replacing the Ideal Exact-Conditional Reference on the Accelerated Branch The ideal contraction proof draws E∼qyE q_y on the accelerated branch, but this draw is used only as an independent reference for the event-driven motion; it is distinct from the actual refresh on the RGO branch. The finite algorithm replaces this ideal reference by a Gaussian centered near the conditional mode. The next comparison lemma controls the resulting normalized shift and is used only to re-establish the local quarter-turn error bound. Fix y, let x∗(y)x_*(y) minimize the conditional potential, and set z=(x−x∗(y))/ηz=(x-x_*(y))/ η. After removing a constant and a linear term, the normalized conditional density has the form q(z)∝e−∥z∥2/2−r(z),mI⪯∇2r⪯LI,q(z) e^- z ^2/2-r(z), mI ∇^2r LI, where m=αηm=αη, L=βηL=βη, and ∇r(0)=0∇ r(0)=0. Lemma 17 (Shifted-Gaussian comparison). Let γδ=(δ,Id) _δ=N(δ,I_d), and set τ=2+log(L/m)τ=2+ (L/m). If τL≤cτ L≤ c and τ(L2d+∥δ∥2)≤cτ(L^2d+ δ ^2)≤ c, then ‖dγδdq‖Lτ(q)≤eCτ(L2d+∥δ∥2),χ2(γδ∥q)≤C(L2d+∥δ∥2). \| d _δdq \|_L^τ(q)≤ e^Cτ(L^2d+ δ ^2), χ^2( _δ\|q)≤ C(L^2d+ δ ^2). Proof. Normalize r(0)=0r(0)=0. Convexity and L-smoothness give 0≤r(z)≤L∥z∥2/20≤ r(z)≤ L z ^2/2 and ∥∇r(z)∥2≤2Lr(z) ∇ r(z) ^2≤ 2Lr(z). The Gaussian log-Sobolev inequality and the standard entropy argument (Chewi, 2023) imply, for 0≤tL<10≤ tL<1, logγet(r−γr)≤t2L2d2(1−tL),logγe−t(r−γr)≤t2L2d2. _γe^t(r-E_γr)≤ t^2L^2d2(1-tL), _γe^-t(r-E_γr)≤ t^2L^2d2. For completeness, the second inequality is not obtained by simply changing the sign in the first. If A−(t)=logγe−t(r−γr)A_-(t)= _γe^-t(r-E_γr), the same entropy calculation gives tA−′(t)−A−(t)≤t2Lγr−A−′(t)tA_- (t)-A_-(t)≤ t^2L\E_γr-A_- (t)\. Convexity and A−(0)=A−′(0)=0A_-(0)=A_- (0)=0 imply A−′(t)≥0A_- (t)≥ 0, so integrating (A−(t)/t)′≤Lγr≤L2d/2(A_-(t)/t) ≤ LE_γr≤ L^2d/2 proves the displayed negative-tilt bound. Applying these bounds to dγ/dqdγ/dq gives logq(dγ/dq)τ≤Cτ2L2d _q(dγ/dq)^τ≤ Cτ^2L^2d. For the shift, dγδ/dγ=expδTz−∥δ∥2/2d _δ/dγ= \δ^Tz- δ ^2/2\. The law q is one-strongly log-concave, so its linear functionals are one-subgaussian after centering. Weak integration by parts gives qZ=−q∇r(Z)E_qZ=-E_q∇ r(Z), while q∥Z∥2≤dE_q Z ^2≤ d; hence ∥qZ∥≤Ld _qZ ≤ L d. Cauchy–Schwarz and Young’s inequality therefore add at most Cτ2(L2d+∥δ∥2)Cτ^2(L^2d+ δ ^2) to the logarithmic moment. Taking the τ-th root proves the first claim. Repeating the calculation at exponent two and using eu−1≤2ue^u-1≤ 2u in the stated smallness regime proves the χ2χ^2 bound. Mollification and cutoff justify the calculation for C1,1C^1,1 potentials. ∎ Starting from y, run gradient descent on the conditional potential until it returns a measurable point x^y x_y satisfying ∥y−η∇V(x^y)−x^y∥≤dη. y-η∇ V( x_y)- x_y ≤ dη. The accelerated branch then draws G∼(0,Id)G (0,I_d) and uses Ey=y−η∇V(x^y)+ηG.E_y=y-η∇ V( x_y)+ η\,G. (41) Strong monotonicity of x↦x+η∇V(x)−yx x+η∇ V(x)-y shows that the normalized shift of 41 from x∗(y)x_*(y) is at most Ld/(1+αη)L d/(1+αη). Thus Lem. 17 applies whenever τL≲1τ L 1 and τL2d≲1τ L^2d 1. Because the reference depends only on y and fresh Gaussian randomness, Lem. 8 gives exact conditional invariance before numerical caps are imposed. Lemma 18 (Shifted-reference transfer). If τL≲1τ L 1 and τL2d≲1τ L^2d 1, the conditional kernel using 41 satisfies Props. 10 and 2, with possibly different universal constants. Proof. Let h=∇rh=∇ r, Σ=Covq(Z) =Cov_q(Z), and D=I−ΣD=I- . For the reference force =h(G+δ)U=h(G+δ), direct expansion of the Gaussian momentum moments in the proof of Lem. 13 gives ka,E2=6Varq(aTh)+5Var(aT)+5q(aTh)−(aT)2.Ek_a,E^2=6Var_q(a^Th)+5Var(a^TU)+5\E_q(a^Th)-E(a^TU)\^2. The weak integration-by-parts calculation used in Lem. 12 yields Varq(aTh)≤CLaTDaVar_q(a^Th)≤ CL\,a^TDa. It remains to control the two Gaussian-reference terms. If H=∇2rH=∇^2r denotes the weak Hessian, then qH=D+Covq(h)⪯CDE_qH=D+Cov_q(h) CD. Hölder’s inequality, Lem. 17, and 0⪯H⪯LI0 H LI give γδ[aTHa]≤CaTDaE_ _δ[a^THa]≤ Ca^TDa. Indeed, aTHa≤L∥a∥2a^THa≤ L a ^2, q[aTHa]≥m∥a∥2E_q[a^THa]≥ m a ^2, and the choice τ=2+log(L/m)τ=2+ (L/m) makes the resulting Hölder ratio universal. Gaussian Poincaré therefore implies Varγδ(aTh)≤CLaTDaVar_ _δ(a^Th)≤ CL\,a^TDa. The χ2χ^2 part of Lem. 17, followed by Cauchy–Schwarz under q, bounds the squared difference of means by the same quantity. Hence the local quarter-turn estimate holds with the shifted reference. The affine lifting in Prop. 10, momentum parity, and the factorization proof are unchanged. ∎ C.4.3 Mode, Clock, and Refresh Subroutines Lemma 19 (Mode certificate and deterministic cap). Assume βη≤1βη≤ 1. Gradient descent with step (β+η−1)−1(β+η^-1)^-1, initialized at y, finds the certificate used in 41 using at most C[1+log+(η∥∇V(y)∥d)]C [1+ _+ ( η ∇ V(y) d ) ] oracle queries. Under the stationary augmented law, stopping after JB=C[1+log(2+κ)+log(1+B/d)]J_B=C [1+ (2+κ)+ (1+ B/d ) ] iterations fails with probability at most 2e−B2e^-B. The uncapped output is a Borel function of y. Proof. The conditional potential has Hessian between α+η−1α+η^-1 and β+η−1β+η^-1, hence condition number at most two. One gradient step therefore reduces the gradient norm by a universal factor, which gives the pointwise logarithmic bound. For the stationary tail, let x⋆x be the global minimizer of V and write Y=X+ηZY=X+ ηZ, where Z∼(0,Id)Z (0,I_d). Strong log-concavity gives ∥X−x⋆∥2≤d/αE X-x ^2≤ d/α and Gaussian concentration around that scale. Combining it with the Gaussian norm tail and ∥∇V(Y)∥≤β∥Y−x⋆∥ ∇ V(Y) ≤β Y-x gives, outside an event of probability 2e−B2e^-B, η∥∇V(Y)∥d≤2κ(1+2B/d). η ∇ V(Y) d≤ 2 κ (1+ 2B/d ). Substitution into the pointwise bound gives JBJ_B. Every finite gradient-descent iterate and the first iterate meeting the certificate are Borel functions of y. Returning the complete pre-reflection state on cap failure therefore defines a total Markov kernel without violating conditional independence. ∎ Lemma 20 (Clock cap). Suppose the normalized residual potential r has L-Lipschitz gradient with 0≤L≤10≤ L≤ 1, and the normalized reference is G+δyG+ _y, where G is standard Gaussian, y↦δy _y is measurable, and supy∥δy∥≤Ld _y _y ≤ L d. For B≥1B≥ 1, set λ¯B=CL(dB+B) λ_B=CL( dB+B). Replacing the event rate by its minimum with λ¯B λ_B changes one stationary transition with probability at most C(1+Ld)e−cBC(1+L d)e^-cB. The capped flight is simulated by homogeneous Poisson thinning and uses at most 1+(π/2)λ¯B1+(π/2) λ_B expected oracle queries, including the query at the reference. Proof. Conditionally on Y=yY=y, the normalized position is one-strongly log-concave, while G and the momentum are Gaussian. Their norm tails, the bound ∥∇r(w)−∇r(e)∥≤L∥w−e∥ ∇ r(w)-∇ r(e) ≤ L w-e , and a one-dimensional Gaussian tail give, uniformly in y, Prλt>CL(du+u)≤6e−u. \ _t>CL( du+u)\≤ 6e^-u. Invariance from Lem. 8 makes the bound valid at every time. Integrating the tail over the quarter-period bounds the expected excess hazard by C(1+Ld)e−cBC(1+L d)e^-cB. Using the same Poisson random measure for the exact and capped clocks couples the paths until the first excess event. Conditioning on Y=yY=y and then integrating proves the claim for the measurable random center. Ordinary thinning gives the query bound. ∎ Proposition 12 (Approximate RGO). There is a universal KrgoK_ rgo such that the following holds. Given y, 0<δ<1/20<δ<1/2, and a point x^y x_y satisfying ∥y−η∇V(x^y)−x^y∥≤dη, y-η∇ V( x_y)- x_y ≤ dη, if η−1≥Krgoβdlog(1/δ)+log(1/δ),η^-1≥ K_ rgoβ\ d (1/δ)+ (1/δ)\, the first-order rejection sampler of Chen et al. (2026a, Theorems 3.1 and D.1) returns q^y q_y with χ2(qy∥q^y)≤δ2χ^2(q_y\| q_y)≤δ^2 and χ2(q^y∥qy)≤δ2χ^2( q_y\|q_y)≤δ^2, using at most KrgoK_ rgo expected oracle queries after construction of the mode certificate. C.4.4 Proof of the Implementation Theorem Proof of Thm. 6. Let Aε=Δ+log(8/ε)A_ = + (8/ ) and take δrgo=ε8(1+Nρ),B=CΔ+log8Nε+log(2+d+κ). _ rgo= 8(1+Nρ), B=C \ + 8N + (2+d+κ) \. The exact comparison chain uses the measurable shifted reference in 41, without mode or clock caps, and an exact refresh. Lemma 18 lets us apply Prop. 11; the choice of N makes its distance from stationarity at most ε/2 /2. The implementable algorithm first reflects y. It then runs the y-initialized mode solver for at most JBJ_B steps. On failure it returns the complete pre-reflection state (x,y,p)(x,y,p). On the refresh branch it invokes Prop. 12 with tolerance δrgo _ rgo; otherwise it uses the shifted Gaussian reference and the clock cap in Lem. 20. Couple both chains until their first disagreement. The complete error allocation is sourcetotal TV contributionideal mixingε/2mode-solver capε/8event-clock capε/8approximate refreshNρδrgo/2≤ε/16. array[]c|csource&total TV contribution\\ mixing& /2\\ mode-solver cap& /8\\ event-clock cap& /8\\ approximate refresh&Nρ _ rgo/2≤ /16. array For the mode and clock errors, Lem. 16 converts stationary exponential tails into path probabilities under the warm exact chain. The RGO error is uniform in y, so maximal coupling and kernel telescoping apply directly. The entries sum to less than ε , and projection to X cannot increase total variation. We next check the local conditions. With τ=2+logκτ=2+ κ and L=βη=c/(dΛε+Λε)L=βη=c/( d _ + _ ), reducing c gives τL≲1τ L 1 and τL2d≲1τ L^2d 1. Also log(1/δrgo)≲Λε (1/ _ rgo) _ , so the same step size verifies the gate in Prop. 12. The expected number of refresh calls is Nρ≤CℓgAεNρ≤ C _gA_ . Each iteration uses at most CJB+1+L(dB+B)+ρKrgoC\J_B+1+L( dB+B)+ρ K_ rgo\ expected oracle queries. Here JBJ_B and ℓg _g are logarithmic, and the definition of η makes L(dB+B)=O(1)L( dB+B)=O(1). Substituting g2=αηg^2=αη proves the theorem and hence Thm. 4. ∎ C.4.5 Affine Form Used by the Quadratic Baseline Lemma 21 (Affine local compiler). Suppose a conditional potential is Φy(x)=Ry(x)+ω2∥x−y∥2,mωI⪯∇2Ry⪯LωI,0<m≤L≤1. _y(x)=R_y(x)+ ω2 x-y ^2, mω I ∇^2R_y Lω I, 0<m≤ L≤ 1. Let x∗(y)x_*(y) denote the unique minimizer of Φy _y, let G∼(0,Id)G (0,I_d), and fix 0<δ<1/20<δ<1/2. Let τ=2+log(L/m)τ=2+ (L/m), and assume τL≲1τ L 1 and τL2d≲1τ L^2d 1. If a measurable point x^y x_y satisfies ∥ω(x^y−y)+∇Ry(x^y)∥≤dω, ω( x_y-y)+∇ R_y( x_y) ≤ dω, then the reference Ey=y−ω−1∇Ry(x^y)+ω−1/2GE_y=y-ω^-1∇ R_y( x_y)+ω^-1/2G satisfies the shifted-Gaussian comparison, the local quarter-turn forcing estimate, and the clock estimate with local smoothness L. The certificate is found from y using C[1+log+∥∇Ry(y)∥dω]C [1+ _+ ∇ R_y(y) dω ] oracle queries. The RGO gate becomes 1≥KrgoLdlog(1/δ)+log(1/δ).1≥ K_ rgoL\ d (1/δ)+ (1/δ)\. Proof. In the coordinate z=ω(x−x∗(y))z= ω(x-x_*(y)), the residual has Hessian between mImI and LILI. Strong monotonicity of x↦ω(x−y)+∇Ry(x)x ω(x-y)+∇ R_y(x) bounds the normalized Gaussian shift by Ld/(1+m)L d/(1+m). Thus Lems. 17 and 20 apply after replacing βηβη by L. The proof of Lem. 18 up to its local quarter-turn estimate applies in these normalized coordinates; lifting that estimate to the global reflected chain requires the separate scaling calculation in Sec. C.5. The condition number of Φy _y is at most (1+L)/(1+m)≤2(1+L)/(1+m)≤ 2, giving the displayed logarithmic mode cost. Finally, applying Prop. 12 after the affine change of variables gives exactly the stated RGO gate. ∎ C.5 Quadratic-Baseline Refinement Corollary 3 (Quadratic-baseline refinement). For 0<ε≤1/20< ≤ 1/2, if α is known and the initial position law ν0 _0 satisfies D2(ν0∥π)=O(1)D_2( _0\|π)=O(1), a quadratic-baseline implementable algorithm under the same augmented initialization attains TV error ε with O~(1+κ−1d1/4) O(1+ κ-1\,d^1/4) expected oracle queries and polylog(d,κ,1/ε)polylog(d,κ,1/ ) expected RGO calls. We prove Cor. 3 by changing only the quadratic part integrated exactly by the conditional motion. The outer reflection and its gap are unchanged. The parameter choice below is easiest to read from its design goals. Write s=κ−1s=κ-1 for the curvature above the known α-quadratic baseline, and let θ be the desired normalized smoothness of the residual force. We choose λ=1−δ0λ=1- _0 and η so that (i) subtracting the λαλα-quadratic leaves residual strong convexity δ0α>0 _0α>0; (i) the residual smoothness divided by the new harmonic curvature is exactly θ; (i) the original small-step condition remains valid, so the outer reflection gap is still governed by g=αηg= αη. The displayed formulas are chosen to meet these requirements; C0C_0 provides fixed slack in the first one. Let s=κ−1s=κ-1 and, for a sufficiently small universal c, set Λε=Δ+logCdκε,θ=cdΛε+Λε. _ = + Cdκ , θ= c d _ + _ . Choose a sufficiently large universal C0C_0, put δ0=C0θ _0=C_0θ and λ=1−δ0λ=1- _0, and select η by αη=θs+(C0−1)θ+C0θ2.αη= θs+(C_0-1)θ+C_0θ^2. (42) Choosing C0C_0 first and then reducing c ensures 0<λ<1,βη≤c0,θ≤c0,0<λ<1, βη≤ c_0, θ≤ c_0, where c0c_0 is the universal small-step constant in Thm. 3. To see the middle inequality, split into s≥θs≥θ and s<θs<θ in βη=(s+1)θ/s+(C0−1)θ+C0θ2βη=(s+1)θ/\s+(C_0-1)θ+C_0θ^2\. For the implementable algorithm, put g=αη,ℓg=1+log(1/g),ρ≍gℓg,Aε=Δ+log(8/ε),g= αη, _g=1+ (1/g), ρ g _g, A_ = + (8/ ), and choose N=⌈CℓgAε⌉,δrgo=ε8(1+Nρ),B=CΔ+log8Nε+log(2+d+κ).N= C _ggA_ , _ rgo= 8(1+Nρ), B=C \ + 8N + (2+d+κ) \. These are the same mixing horizon and error budgets used in Thm. 6, evaluated at the baseline step size in 42. For fixed y, split the conditional potential as V(x)+∥x−y∥22η=Ry(x)+ωη2∥x−y∥2,Ry(x)=V(x)−λα2∥x−y∥2,ωη=λα+η−1.V(x)+ x-y ^22η=R_y(x)+ _η2 x-y ^2, R_y(x)=V(x)- λα2 x-y ^2, _η=λα+η^-1. The residual curvature satisfies δ0αI⪯∇2Ry⪯(β−λα)I,β−λαωη=θ,ωηη≍1. _0α I ∇^2R_y (β-λα)I, β-λα _η=θ, _η 1. Thus the harmonic motion integrates both the proximal quadratic and the known λαλα-quadratic baseline; collisions see only the remaining curvature. Local validity and implementation. In the coordinate z=ωη(x−y)z= _η(x-y), use the same fixed-reference process as in 23, with byres(z)=ωη−1/2∇Ry(y+ωη−1/2z)b_y res(z)= _η^-1/2∇ R_y(y+ _η^-1/2z) in place of byb_y. Lemma 8 applies after this affine change of variables, so the process still preserves the original conditional law. Its normalized residual curvatures are m⋆=δ0αωη,L⋆=θ.m_ = _0α _η, L_ =θ. If τ⋆=2+log(L⋆/m⋆) _ =2+ (L_ /m_ ), then L⋆m⋆=ωηC0α=1+sC0θ,τ⋆≲Λε. L_ m_ = _ηC_0α=1+ sC_0θ, _ _ . Consequently τ⋆L⋆≲1 _ L_ 1 and τ⋆L⋆2d≲1 _ L_ ^2d 1, after reducing c. The definition above also gives log(1/δrgo)≤CΛε (1/ _ rgo)≤ C _ . Hence Krgoθdlog(1/δrgo)+log(1/δrgo)≤1K_ rgoθ\ d (1/ _ rgo)+ (1/ _ rgo)\≤ 1 after another universal reduction of c. This is exactly the affine RGO gate in Lem. 21. The reference and its certificate are ∥ωη(x^y−y)+∇Ry(x^y)∥≤dωη,Ey=y−ωη−1∇Ry(x^y)+ωη−1/2G. _η( x_y-y)+∇ R_y( x_y) ≤ d _η, E_y=y- _η^-1∇ R_y( x_y)+ _η^-1/2G. In the remainder of this subsection, S denotes the rescaled conditional quarter-turn averaged over this shifted-Gaussian reference. Lemma 21 therefore supplies the shifted Gaussian comparison, event-clock bound, mode solver, and RGO gate with local smoothness L⋆=θL_ =θ. The stationary mode cap remains logarithmic: since ∇RY(Y)=∇V(Y)∇ R_Y(Y)=∇ V(Y) and ωη≥η−1 _η≥η^-1, ηd∥∇V(Y)∥2≤κ(βη)+(βη)2≤κ+1. ηdE ∇ V(Y) ^2≤κ(βη)+(βη)^2≤κ+1. The concentration argument in Lem. 19 then applies without a polynomial query factor. Transfer estimate after rescaling. It remains to verify that the ideal contraction theorem sees the same position–momentum geometry. Let Σy=η−1Cov(X∣Y=y),Dy=I−Σy,Dz(y)=I−ωηηΣy. _y=η^-1Cov(X Y=y), D_y=I- _y, D_z(y)=I- _η\, _y. Thus Σy,Dy _y,D_y are exactly the conditional covariance quantities defined in Sec. C.2, while DzD_z is their rescaled counterpart. Because ωηη≥1 _η≥ 1, Dz(y)=ωηηDy−(ωηη−1)I⪯ωηηDy.D_z(y)= _η D_y-( _η-1)I _η D_y. The correctly scaled momentum observable is Mf(α)=2ωηaf(Y)TP.M_f^(α)= 2 _ηa_f(Y)^TP. Applying the local quarter-turn estimate in the z-coordinate yields ∥Lf−Mf(α)∥2+∥Mf(α)+Lf∥2≤Cθηf. SL_f-M_f^(α) ^2+ SM_f^(α)+L_f ^2≤ Cθ _f. Its defect under the outer reflection is at most CηℛfC _f, and ∥Mf(α)∥2≤Cη∥∇f∥2 M_f^(α) ^2≤ Cη ∇ f ^2. Combining these bounds with 32 and 11 proves the same position–momentum transfer estimate 29. Momentum parity is unchanged, so Cors. 2 and 14 apply with the same universal constants. Notice that the outer gap is still g=αηg= αη; only the local collision scale changes from βηβη to θ. Complexity. From 42, 1αη=sθ+C0−1+C0θ, 1αη= sθ+C_0-1+C_0θ, and hence 1αη≤C(1+κ−1θ)=O~(1+κ−1d1/4). 1 αη≤ C (1+ κ-1 θ )= O\! (1+ κ-1\,d^1/4 ). Also θd=O(1)θ d=O(1), so the clock uses only O~(1) O(1) expected oracle queries per iteration. The ideal contraction bound and our choice of N give TV error at most ε/2 /2. With the explicit B above, Lems. 16, 19 and 20 bound the two cap errors by ε/8 /8 each. Uniform maximal coupling for the approximate refresh contributes at most Nρδrgo/2≤ε/16Nρ _ rgo/2≤ /16. Thus the complete error is below ε . The per-iteration mode and clock costs are polylogarithmic, while NρNρ is polylogarithmic, so the same accounting gives TV(Law(XN),π)≤ε,Nqry=O~(1+κ−1d1/4),Nrgo=polylog(d,κ,1/ε)TV(Law(X_N),π)≤ , _qry= O\! (1+ κ-1\,d^1/4 ), _rgo=polylog(d,κ,1/ ) when Δ=O(1) =O(1). This proves Cor. 3. Appendix D Proof of the Warm-Start Wrapper This appendix proves both parts of Thm. 5. For the finite-divergence part, we first state two implementable samplers for the localized conditionals and a capped interface for the final reflected sampler. We then show that the ideal coarse proximal chain reduces a finite D2D_2 budget to a constant, couple those steps to an algorithmic chain at adaptive observation centers, optimize the coarse step size, and transfer the final guarantee through the coupling. The last subsection proves the supplied-point part by constructing the required finite D2D_2 budget from one point. The proof uses two external posterior-sampling results, stated precisely in Prop. 13, and the standard exact-proximal hypercontractive inequality stated in Lem. 23. All remaining steps—including the random observation centers, the distinction between ideal coarse and algorithmic laws, the final error transfer, and the parameter optimization—are proved below. At a high level, the wrapper compares two routes. The direct route applies a reflected sampler immediately to the initial law. The coarse route first constructs two endpoints in parallel for the analysis: an ideal coarse endpoint νm† _m , which is proved to be constant-D2D_2-warm, and an algorithmic endpoint ν^m ν_m, which is close to νm† _m in total variation. The same capped reflected kernel is then applied to both endpoints. Only ν^m ν_m is generated by the algorithm; νm† _m is the ideal coarse endpoint used to prove convergence. D.1 Posterior Samplers Used by the Wrapper At a coarse step size u>0u>0, an observation y defines the conditional density qy(x)∝exp−V(x)−∥x−y∥22u.q_y(x) \! \-V(x)- x-y ^22u \. Its strong-convexity and smoothness parameters are mu=α+u−1m_u=α+u^-1 and Lu=β+u−1L_u=β+u^-1, so its condition number is Ku=Lu/muK_u=L_u/m_u. We use two ways to sample this conditional. The regular interface first constructs a constant-warm law and then applies one of our reflected samplers. The harmonic-base interface instead starts from a localized Gaussian and invokes exact-ULD/FORS. We abbreviate the latter by “HB” because its reference is the quadratic, or harmonic, part of the conditional potential; HB does not mean a heat-bath update. The regular interface is advantageous in the intermediate- and large-dimensional regimes, whereas the HB interface supplies the small-dimensional endpoint. We will repeatedly turn expected-work routines into total kernels with a deterministic work limit. The following elementary observation records why this does not introduce a polynomial dependence on the requested accuracy. Lemma 22 (Logarithmic work truncation). Suppose a routine has a deterministic number of stages and every random-length operation is one of the following: a geometric trial count with uniformly positive success probability; a Poisson count with known mean Λ ; or a linearly convergent optimization for which N(z)≤C1+log(1+H(z)/a),H(Z)2≤M2N(z)≤ C\1+ (1+H(z)/a)\, (Z)^2≤ M^2 under the comparison law relevant to that call. For every ξ>0ξ>0, all such operations can be stopped at deterministic limits so that the chance of any stop is at most ξ. The limits add only logarithmic factors to the original work, apart from the already charged Poisson mean Λ . Returning a fixed fallback output after a stop changes the output law by at most ξ in total variation. Proof. Geometric tails require O(log(1/ξ))O( (1/ξ)) trials. A Poisson variable of mean Λ is at most Λ+CΛlog(1/ξ)+log(1/ξ) +C\ (1/ξ)+ (1/ξ)\ except with probability ξ. Finally, Markov’s inequality gives H(Z)≤M/ξH(Z)≤ M/ ξ outside an event of probability ξ; because N depends logarithmically on H, this increases its limit by only O(log(1/ξ))O( (1/ξ)). Divide ξ by the deterministic number of operations and use a union bound. Coupling the capped and uncapped routines until the first stop proves the total-variation statement. ∎ Proposition 13 (Finite posterior interfaces). Fix a law ν with D2(ν∥π)≤BD_2(ν\|π)≤ B, draw X∼νX ν and, independently, a standard Gaussian vector G∼(0,Id)G (0,I_d), where IdI_d is the d-dimensional identity; set Y=X+uGY=X+ uG, and denote the law of Y by λ. Suppose u≤1/αu≤ 1/α. For every 0<δ<1/20<δ<1/2, each of the two interfaces above gives a jointly measurable, work-capped kernel y↦P^y P_y such that ∫TV(P^y,qy)λ(y)≤δ. ( P_y,q_y)\,λ(dy)≤δ. (43) For the regular interface, every call uses at most Creg(Ku,d):=Ku7/6d1/6+Ku1/2d1/4C_ reg(K_u,d):=K_u^7/6d^1/6+K_u^1/2d^1/4 (44) oracle queries up to factors logarithmic in d,Ku,B,1/δd,K_u,B,1/δ. For the HB interface, the corresponding bound is CHB(Ku,d):=1+Ku2/3d1/3C_ HB(K_u,d):=1+K_u^2/3d^1/3 (45) oracle queries, again up to the same type of logarithmic factors. Proof. Regular interface. After rescaling a potential whose condition number is K, Chen et al. (2026c, Theorem 7.1(i)) constructs the output law of an implementable algorithm and an analysis-only comparison law that are arbitrarily close in total variation, with the latter having constant order-two Rényi divergence from the target. At constant target divergence its query cost is O~(K7/6d1/6+K) O(K^7/6d^1/6+K); the second term is absorbed by the first for K,d≥1K,d≥ 1. Applying either reflected finalizer from this paper to the comparison law and the same capped Markov kernel to the algorithmic output law costs another O~(K1/2d1/4) O(K^1/2d^1/4). Total- variation data processing and the triangle inequality give 44. HB interface. First find a point zyz_y satisfying the observable test ∥∇Φy(zy)∥≤mu/Lu ∇ _y(z_y) ≤ m_u/ L_u, where Φy(x)=V(x)+∥x−y∥2/(2u) _y(x)=V(x)+ x-y ^2/(2u). Strong convexity then gives Lu∥zy−xy∗∥≤1 L_u z_y-x_y^* ≤ 1, where xy∗x_y^* minimizes Φy _y. The Gaussian Gy=(zy,(2Lu)−1Id)G_y=N(z_y,(2L_u)^-1I_d) then obeys Dq(Gy∥qy)=Oq(dlog(eKu))D_q(G_y\|q_y)=O_q(d (eK_u)) for every fixed q>1q>1, by the usual quadratic sandwich and completion of squares. Applying Chen et al. (2026d, Theorem 3.2(i)) with output order two, initialize phase space from GyG_y tensored with stationary Gaussian momentum. The source theorem’s initial budget is then d+D4(Gy∥qy)=O~(d)d+D_4(G_y\|q_y)= O(d). Set its Rényi tolerance to O(δ2)O(δ^2), project the output to position, and use TV2≤D2/2TV^2≤ D_2/2. This gives 45. Its high-probability query bound also includes rejected trials. Random centers and work caps. It remains to justify localization at the random center y. The conditional objective has gradient ∇V(y)∇ V(y) at y. For the regular interface, the observable test ∥∇Φy(z)∥≤2mud ∇ _y(z) ≤ 2m_ud implies Φy(z)−minΦy≤d _y(z)- _y≤ d, as required by the imported warmifier; for the HB interface we use the test in the preceding paragraph. Accelerated optimization started at y reaches either certificate in O~(Ku) O( K_u) queries, with only logarithmic dependence on 1+∥y−x∗∥1+ y-x_* . This cost is absorbed by both displayed bounds. Under the law λ, the tail needed to cap this search uniformly follows from Lem. 24 below. On cap failure the kernel returns a fixed point; the cap probability is included in 43. The cap depends only on the known parameters α,β,d,B,δα,β,d,B,δ, not on the unknown law λ. We next verify the work cap rather than inferring it from an expected-query statement. The recursive warmifier has a deterministic number of stages. Its proof controls the second moments of every center passed to its gradient-descent subroutines, while its pointwise bound for each such solve is logarithmic in the center gradient; see Chen et al. (2026c, Theorems 6.1(i) and 6.5, their proofs, and Lemmas D.3–D.4). Thus these solves satisfy the third case of Lem. 22. Its remaining rejection loops have geometric tails once their observable mode certificate holds; this is stated for the RGO component in Chen et al. (2026a, Theorems 3.1 and D.1). The appended reflected finalizer in the regular interface has only geometric rejection counts and Poisson event counts. The exact-ULD/FORS interface already gives a high-probability query bound. Applying Lem. 22 makes the displayed work bounds deterministic up to their logarithmic factors and charges every stop to 43. Finally, finite optimization iterates, Gaussian draws, acceptance tests, and capped stopping rules are Borel functions of y and the random seed. Thus P^y P_y is a genuine Markov kernel rather than a pointwise collection of samplers. ∎ We next isolate the final-stage property needed by the wrapper. Stating it as a capped Markov-kernel interface is important: after the algorithmic and comparison paths separate, the final sampler must still have a deterministic work limit. Proposition 14 (Capped reflected finalizer). Fix B≥1B≥ 1 and 0<ε,ξ<1/20< ,ξ<1/2. Either reflected sampler admits an implementable realization as a total Markov kernel K^B,ε,ξ K_B, ,ξ such that, for every ν satisfying D2(ν∥π)≤BD_2(ν\|π)≤ B, TV(νK^B,ε,ξ,π)≤ε+ξTV(ν K_B, ,ξ,π)≤ +ξ. Up to logarithmic factors, the kernel uses at most O~(κB(dB+B)1/2) O\! ( κB( dB+B)^1/2 ) (46) oracle queries. In particular, constant B costs O~(κd1/4) O( κd^1/4). The same kernel may be applied to a second, arbitrary input law: its deterministic work bound remains valid, and total-variation distance between the two inputs cannot increase. Proof. The uncapped guarantees are Cor. 1 and 8. The implementable algorithms have a fixed outer horizon. Their mode solves have logarithmic pointwise work bounds and warm-path second-moment control, their rejection counts are geometric, and their event proposal counts are Poisson with the means already included in their query-complexity bounds. Hence Lem. 22 applies. For any additional error budget ξ, each implementation can therefore be replaced by a total Markov kernel with a deterministic work bound of the same polynomial order and with at most ξ additional total-variation error from every D2D_2-warm input. For the flux method the required estimates are in Props. 8 and 7; for the persistent method it follows from Lems. 16, 19 and 20 together with the same elementary Poisson and geometric tail bounds. We use this capped version below so that the query bound remains valid even if the algorithmic and comparison paths have already separated. The final assertion is the data-processing inequality for a common Markov kernel. ∎ D.2 Ideal Coarse Proximal Contraction Let SuS_u denote one exact proximal Gibbs update: observe Y=X+uGY=X+ uG, then draw X′∼qYX q_Y. The kernel is π-reversible. The following standard hypercontractive form makes its warm-start effect explicit. Lemma 23 (Exact proximal warmification). For every density f=dν/dπf=dν/dπ, every p>1p>1, and every u>0u>0, ∥Suf∥1+(1+αu)2(p−1)≤∥f∥p. S_uf _1+(1+α u)^2(p-1)≤ f _p. Consequently, if Gm=∏j=1m(1+αuj)G_m= _j=1^m(1+α u_j), then D2(νSu1⋯Sum∥π)≤2D2(ν∥π)1+Gm2.D_2(ν S_u_1·s S_u_m\|π)≤ 2D_2(ν\|π)1+G_m^2. (47) Proof. Let πu _u be the observation marginal and define Au∗f(y)=[f(X)∣Y=y]A_u^*f(y)=E[f(X) Y=y], which maps density ratios on π to density ratios on πu _u. Its adjoint is Aug(x)=[g(Y)∣X=x]A_ug(x)=E[g(Y) X=x], and one Gibbs row has density operator Su=AuAu∗S_u=A_uA_u^*. Put cu=1+αuc_u=1+α u. The two-measure heat-flow inequality of Kook and Vempala (2025, Theorem 1.2) states that Au∗:Lp0(π)→Lp1(πu)A_u^*:L^p_0(π)→ L^p_1( _u) is a contraction whenever p1−1=cu(p0−1)p_1-1=c_u(p_0-1). Duality gives the analogous gain for AuA_u: Au:Lp1(πu)→Lp2(π)A_u:L^p_1( _u)→ L^p_2(π) is a contraction when p2−1=cu(p1−1)p_2-1=c_u(p_1-1). To see the direction of this second relation, the conjugate-exponent form of the first one is p0′−1=cu(p1′−1)p_0 -1=c_u(p_1 -1), where a prime denotes the Hölder conjugate. Thus the exponent flow through one Gibbs row is Lp0(π)→Au∗Lp1(πu)→AuLp2(π),p2−1=cu2(p0−1).L^p_0(π) \ A_u^*\ L^p_1( _u) \ A_u\ L^p_2(π), p_2-1=c_u^2(p_0-1). The two excess-exponent gains therefore multiply and prove the first claim. After iteration it gives ∥Su1⋯Sumf∥2≤∥f∥1+Gm−2 S_u_1·s S_u_mf _2≤ f _1+G_m^-2. For p=1+Gm−2p=1+G_m^-2, interpolation between L1(π)L^1(π) and L2(π)L^2(π), using ∥f∥1=1 f _1=1, gives log∥f∥p≤(p−1)D2(ν∥π)/p f _p≤(p-1)D_2(ν\|π)/p. Since the output order-two divergence is twice the logarithm of its L2L^2 density norm, 47 follows. ∎ Set r=βur=β u and restrict 1≤r≤κ1≤ r≤κ. Then αu=r/κα u=r/κ, and log(1+r/κ)≍r/κ (1+r/κ) r/κ. Therefore 47 makes the ideal coarse endpoint constant-D2D_2 warm after m=O(κrlog(eB))m=O\! ( κr (eB) ) (48) equal-step-size updates. Notice also that their conditional condition number is Ku=β+u−1α+u−1=κ(1+r)κ+r≤1+r≤2r.K_u= β+u^-1α+u^-1= κ(1+r)κ+r≤ 1+r≤ 2r. (49) D.3 Random Centers and the Algorithmic Path The posterior centers in the preceding chain are adaptive: the next center is sampled from the current position. We now verify that the conditional interfaces can be applied without assuming that these centers are fixed or uniformly bounded. Lemma 24 (Random-center and path accounting). Let ν0 _0 satisfy D2(ν0∥π)≤BD_2( _0\|π)≤ B, and run m exact proximal updates with step sizes at most 1/α1/α. Then every position law νj _j satisfies νj∥X−x∗∥2≤C(d+B)α.E_ _j X-x_* ^2≤ C(d+B)α. Moreover, all observation-dependent mode searches along the path can be capped so that the probability of any cap failure is at most δ, while their query cost changes only by a factor logarithmic in d,κ,B,m,1/δd,κ,B,m,1/δ. Proof. Every exact proximal kernel preserves π, so data processing gives D2(νj∥π)≤BD_2( _j\|π)≤ B. Hence KL(νj∥π)≤BKL( _j\|π)≤ B, and the T2T_2 inequality for an α-strongly log-concave law gives W22(νj,π)≤2B/αW_2^2( _j,π)≤ 2B/α. Combining an optimal coupling with π∥X−x∗∥2≤d/αE_π X-x_* ^2≤ d/α proves the moment bound. Since Y=X+uGY=X+ uG and u≤1/αu≤ 1/α, the same order of control holds for the observation. For the tail statement, let ℙν0(m)P_ _0^(m) be the law of the entire exact path, including its observations, and let ℙπ(m)P_π^(m) be the law of the same path started at stationarity. Because their transition kernels agree, dℙν0(m)dℙπ(m)=dν0dπ(X0),ℙν0(m)(E)≤eB/2ℙπ(m)(E)1/2 dP_ _0^(m)dP_π^(m)= d _0dπ(X_0), _ _0^(m)(E) ^B/2P_π^(m)(E)^1/2 for every path event E, where the second inequality is Cauchy–Schwarz. Under the stationary path, strong log-concavity and the independent Gaussian observation give exponential tails for ∥Y−x∗∥ Y-x_* . Because the conditional mode solver converges geometrically and its stopping test is observable, a threshold with stationary failure probability at most e−Bδ2/(4m2)e^-Bδ^2/(4m^2) increases its iteration limit only logarithmically. A union bound followed by the displayed path inequality makes the probability of any cap failure at most δ. This bound is for the whole adaptive path, so no independence between successive centers is used, and its threshold is determined entirely by the displayed problem parameters. ∎ Lemma 25 (Finite approximation of the coarse path). Let νm†=ν0Sum _m = _0S_u^m be the ideal coarse endpoint with m as in 48. Using either posterior interface in Prop. 13, there is an algorithmic endpoint ν^m ν_m such that TV(ν^m,νm†)≤δTV( ν_m, _m )≤δ. Its number of oracle queries is deterministically bounded, up to logarithmic factors, by mCreg(Ku,d)mC_ reg(K_u,d) or mCHB(Ku,d)mC_ HB(K_u,d), respectively. Proof. The row-wise coupling has the following form as long as the two paths have not yet disagreed: Xj−1†→Yj†=Xj−1†+uGjYj†→qYj†Xj†,X^j−1→Y^j=X^j−1+uGjY^j→P^j,Y^jX^j. array[]cX_j-1 & \;Y_j =X_j-1 + uG_j\;&Y_j & \;q_Y_j \;&X_j ,\\[2.5pt] X_j-1& \; Y_j= X_j-1+ uG_j\;& Y_j& \; P_j, Y_j\;& X_j. array The shared Gaussian GjG_j makes Yj†=Y^jY_j = Y_j whenever Xj−1†=X^j−1X_j-1 = X_j-1. At that common observation, we use a maximal coupling of the two conditional draws. At row j, apply Prop. 13 with total error budget δ/mδ/m; this budget includes the posterior approximation and all of its mode and work caps. Couple the exact and approximate conditional draws maximally whenever the two paths have agreed so far. To avoid conditioning on the agreement event, first generate the entire ideal coarse path and let ej(y)=TV(P^j,y,qy)e_j(y)=TV( P_j,y,q_y). The probability that the first mismatch occurs in row j is at most [1no earlier mismatchej(Yj)]E[1_no earlier mismatche_j(Y_j)], which is bounded by the unconditional exact-path average ej(Yj)≤δ/mEe_j(Y_j)≤δ/m. Summing over the rows gives total mismatch at most δ. The simultaneous choice of all random-center caps is justified by Lem. 24. This is ordinary total-variation telescoping; a Hellinger chain rule is not needed because the requested accuracy affects both posterior interfaces only logarithmically. Crucially, every call is work-capped. Thus the query bound remains valid after the first coupling disagreement, even though the later algorithmic centers no longer have the exact center law. The ideal coarse endpoint νm† _m , rather than the algorithmic endpoint, is the law to which 47 is applied. ∎ D.4 Cost Optimization and Final Error Transfer Proof of the finite-divergence part of Thm. 5. The optimization selects one of the following three coarse-stage implementations. The table states the resulting choices; the algebra below derives them. Regime Interface r Coarse-stage cost d≤κ16/7d≤κ^16/7 HB r=κr=κ O~(κ2/3d1/3log(eB)) O(κ^2/3d^1/3 (eB)) κ16/7<d<κ8κ^16/7<d<κ^8 Regular r=d1/8r=d^1/8 O~(κd3/16log(eB)) O(κ d^3/16 (eB)) d≥κ8d≥κ^8 Regular r=κr=κ O~(κd1/4log(eB)) O( κd^1/4 (eB)) Fix 1≤r≤κ1≤ r≤κ, set u=r/βu=r/β, and use the regular posterior interface. Multiplying 44 by the row count in 48, and using 49, gives the coarse-stage query complexity Qreg(r)=O~(κr−1+r1/6d1/6+r−1/2d1/4log(eB)).Q_ reg(r)= O\! (κ\r^-1+r^1/6d^1/6+r^-1/2d^1/4\ (eB) ). The r−1r^-1 term records the constant work of each row. If d≤κ8d≤κ^8, the admissible choice r=d1/8r=d^1/8 balances the last two terms and gives O~(κd3/16log(eB)) O(κ d^3/16 (eB)). If d≥κ8d≥κ^8, choose r=κr=κ. The two nonconstant terms are then κ7/6d1/6κ^7/6d^1/6 and κd1/4 κd^1/4; the first is no larger than the second precisely when d≥κ8d≥κ^8. For the remaining endpoint, again choose r=κr=κ, but use the HB interface in 45. There are only O(log(eB))O( (eB)) rows and Ku=(κ+1)/2K_u=(κ+1)/2, so this route costs O~(κ2/3d1/3log(eB)) O(κ^2/3d^1/3 (eB)). The harmonic and regular-interface bounds cross at d=κ16/7d=κ^16/7, while the two regular choices cross at d=κ8d=κ^8. Selecting the cheapest route gives exactly O~(Ψ(κ,d)log(eB)) O( (κ,d) (eB)). We now close the distributional argument. Choose the constants in 48 so that D2(νm†∥π)≤1D_2( _m \|π)≤ 1, and use Lem. 25 with δ=ε/4δ= /4. Lift both νm† _m and ν^m ν_m through the same Gaussian augmentation required by the selected finalizer, and append the same independent Gaussian momentum when needed. The lifted ideal coarse law has the same density ratio relative to stationarity, while total-variation data processing preserves the ε/4 /4 discrepancy. Invoke Prop. 14 with final-stage convergence tolerance ε/4 /4 and work-cap budget ε/4 /4, and apply the resulting single kernel to both augmented laws. Its output from the ideal coarse law is within ε/2 /2 of stationarity. Because applying the same Markov kernel cannot increase total-variation distance, the two final outputs remain within ε/4 /4 of one another. The triangle inequality therefore puts the algorithmic output within 3ε/43 /4 of stationarity. The finalizer’s O~(κd1/4) O( κd^1/4) query cost is bounded by every branch of Ψ in its selected regime. This proves the second term in the minimum of 9. Alternatively, skip the coarse stage and apply the general-budget bound 46 directly to ν0 _0; this proves the first term. Finally, every branch of Ψ is at most κ6/7d1/4κ^6/7d^1/4. For the first branch the ratio is κ−4/21d1/12≤1κ^-4/21d^1/12≤ 1 because d≤κ16/7d≤κ^16/7. For the second it is κ1/7d−1/16≤1κ^1/7d^-1/16≤ 1 because d≥κ16/7d≥κ^16/7, and the third follows from κ≥1κ≥ 1. This proves the claimed envelope. ∎ D.5 Initialization from One Point Let x∗x_* be the minimizer of V, let log+(t)=max1,logt _+(t)= \1, t\, and recall the localization budget Qloc:=κlog+[κ2+β∥x0−x∗∥].Q_ loc:= κ\, _+\! [κ\2+ β x_0-x_* \ ]. The quantity QlocQ_ loc appears only in the query bound: the algorithm does not know x∗x_* and stops using the observable gradient test stated below. Lemma 26 (Observable localization). Starting from x0x_0, accelerated optimization followed by one explicit Gaussian draw produces a law q0q_0 such that D2(q0∥π)≤1+d2log(2κ)D_2(q_0\|π)≤ 1+ d2 (2κ) using O~(Qloc) O(Q_ loc) oracle queries. Proof. Run accelerated gradient descent until the observable certificate ∥∇V(xref)∥≤α/β ∇ V(x_ ref) ≤α/ β holds. The standard strongly-convex accelerated rate and V(x0)−V(x∗)≤β∥x0−x∗∥2/2V(x_0)-V(x_*)≤β x_0-x_* ^2/2 give the stated QlocQ_ loc bound. Strong convexity then implies β∥xref−x∗∥≤1 β x_ ref-x_* ≤ 1. Draw q0=(xref,(2β)−1Id)q_0=N(x_ ref,(2β)^-1I_d). Strong convexity bounds the target normalizer by Z≤e−V(x∗)(2π/α)d/2Z ^-V(x_*)(2π/α)^d/2, while smoothness gives V(x)≤V(x∗)+β∥x−x∗∥2/2V(x)≤ V(x_*)+β x-x_* ^2/2. Completing the square yields ‖dq0dπ‖∞≤e(2κ)d/2. \| dq_0dπ \|_∞ (2κ)^d/2. An L∞L^∞ density-ratio bound by M implies D2≤logMD_2≤ M, proving the result. ∎ Proof of the supplied-point part of Thm. 5. Apply Lem. 26, followed by Thm. 5 with B=1+(d/2)log(2κ)B=1+(d/2) (2κ). Since log(eB) (eB) is polylogarithmic in d and κ, the coarse-wrapper branch costs O~(Ψ(κ,d)) O( (κ,d)). Adding the one-time localization cost proves the first inequality in 10; the envelope proved above gives the second. ∎