Algorithmic Theory of ODEs and Sampling from Well-conditioned Logconcave Densities

1 Introduction

• The motivating gap is that discretized sampling algorithms still have dimension-dependent evaluation costs, despite dimension-independent bounds for continuous processes.For logistic regression, earlier work achieved dimension dependence d^1/4 under bounded incoherence, motivating the question of avoiding polynomial dependence on dimension.Source: §1 · pp. 3–4

1.1 Results

• The algorithms reduce dimension dependence to polylogarithmic function calls for the smooth, strongly logconcave class, and achieve polylogarithmic parallel depth.The paper also reports runtime and depth polynomial in log(1/ε), contrasting with prior nearly linear-time algorithms whose dependence was polynomial in 1/ε.Source: §1.1 · p. 5
• For strongly convex functions with Lipschitz gradient, the paper gives a faster sampling algorithm without requiring further smoothness assumptions.This result illustrates that the contraction bound and ODE theorem apply beyond infinitely differentiable functions.Source: §1.1 · p. 5

1.2 HMC and improved contraction rate

• The improved exact-HMC contraction rate is κ^1.5 rather than the previous κ^2 bound, and the analysis shows a larger step size is possible.A separate optimal-contraction bound cannot establish overall mixing because its time parameter depends on an unknown coupling point.Source: §1.1 · p. 6, §1.2 · p. 7, §1.3 · p. 8

1.1 Results

• For the ODEs arising in HMC, the paper gives a nearly linear-time, polylogarithmic-depth solver when solutions are close to a small known basis span.The ODE guarantee assumes a Lipschitz ODE function and a piecewise-polynomial approximation; its complexity is nearly that of representing the solution.Source: §1.1 · pp. 6–7

1.3 Techniques

• The ODE approximation strategy uses piecewise degree-two polynomials under Hessian bounds and low-degree polynomials for smooth functions.Degree bounds are obtained through the Cauchy-Kowalevsky method, majorants, and radius-of-convergence bounds for one-variable ODEs.Source: §1.3 · p. 8
• The paper focuses on polylogarithmic-depth results and omits decomposable-function algorithms that sample individual terms.The authors expect their techniques may also be improved for decomposable functions.Source: §1.3 · p. 8

2 ODE Solver for any basis

• The collocation method represents the solution using a basis and applies fixed-point iteration to an operator built from that basis.The section motivates this as a way to move beyond the low-degree-polynomial restriction of earlier collocation work.Source: §2 · p. 8, §2.2 · pp. 10–11

2.1 Picard-Lindelöf theorem

• If F is L-Lipschitz in x, the j-fold Picard operator has Lipschitz constant at most (LT)^j/j!, yielding a contraction for some j.The bound is stated for any norm on the state space.Source: §2.1 · p. 9

2.2 Intuition of collocation method

• The method assumes the solution derivative is approximable in a D-dimensional basis, which may consist of polynomials, piecewise polynomials, or Fourier functions.The basis assumption requires its integral to be bounded; the paper later establishes an O(1)-bounded basis for low-degree polynomial and piecewise-polynomial spaces on Chebyshev nodes.Source: §2.2 · pp. 10–11

2.3 Collocation method

• Collocation evaluates the ODE operator at basis nodes and applies a precomputed matrix; each iteration takes eO(dD) time for the bases considered.The representation stores the solution through its values at the nodes.Source: §2.3 · p. 11
• For the general k-th-order ODE result, the algorithm uses O(D log(C/epsilon)) evaluations of F.The theorem reduces higher-order ODEs to first-order ODEs.Source: §2.3 · pp. 12–13
• For piecewise-polynomial bases, the k-th-order result gives a piecewise-polynomial solution using O(D log(C/epsilon)) evaluations of F.The stated setting assumes known intervals and polynomial degrees on each interval.Source: §2.3 · pp. 13–14

2.4 Proof of first order ODE

• The collocation operator has Lipschitz constant at most gamma_phi LT, and the algorithm's output is obtained by iterating that operator.The bound is used in the proof of the first-order guarantee.Source: §2.4 · pp. 14, 16

2.5 A basis for piece-wise polynomials

• Every piecewise-polynomial space with fixed nodes admits a 2000-bounded basis.The construction concatenates bases for the individual intervals.Source: §2.5 · pp. 17–18
• For distinct pieces, basis-function supports are disjoint or nested relative to the interval endpoint.Source: §2.5 · p. 18
• Experiments suggest the proposed basis is 1-bounded.Source: §2.5 · p. 18
• Theorem 2.12 cites [DGR96, Section 5] for Lagrange basis polynomials on Chebyshev nodes.Source: §2.5 · p. 18

3 Improved Contraction Bound for HMC

• The section first establishes the improved contraction bound for exact HMC, then analyzes noisy HMC by coupling approximate and ideal iterates.The noisy-process argument uses the per-step approximation error when bounding the distance between coupled iterates.Source: §3 · pp. 19, 22, 23
• The ideal-HMC contraction analysis reduces the problem to a matrix ODE with a Hessian bounded between m2 I and M2 I.The proof then bounds the ODE's error dynamics, including by separating movement parallel and orthogonal to the initial direction.Source: §3 · pp. 19–20
• Theorem 3.2 gives a contraction guarantee for noisy HMC under strong convexity and Lipschitz-gradient assumptions, and states an iteration bound.The supplied theorem passages do not expose the bound's full expression.Source: §3 · p. 22

4 Strongly Convex functions with Lipschitz Gradient

• The section gives a faster HMC-based sampling algorithm for strongly convex functions with Lipschitz gradient, while noting that specialized algorithms may improve its bound.The section's stated purpose is to show the contraction bound and ODE theorem also apply without infinite differentiability.Source: §4 · p. 23

4.1 Bounding the ODE solution

• For functions with bounded Hessian eigenvalues, an HMC trajectory can be approximated by continuously differentiable, piecewise degree-2 polynomials.The approximation result requires bounded Hessian eigenvalues; the accompanying lemmas do not require convexity.Source: §4.1 · pp. 23, 24, 26
• Starting the HMC walk at the minimum makes the gradient norm O(√(M2d)) for most iterations, a bound the paper says is tight for quadratic functions.This controls a term that the section identifies as limiting the runtime analysis.Source: §4.1 · p. 24

4.2 Sampling

• Starting at the minimum, Theorem 1.6 gives a Wasserstein-distance guarantee for sampling from a density proportional to e^-f, with gradient evaluations counted amortized over iterations.The theorem applies to strongly convex functions with Hessian eigenvalues bounded between m2 and M2.Source: §4.2 · p. 27

5 Sampling from Incoherent Logistic Loss Functions and More

• The proof represents the dynamics with s(t) = Ax(t), whose force separates into a linear term and a coordinatewise term after ignoring m2s.This structure is used to approximate the solution of the transformed equation by collocation.Source: §5 · p. 28

5.1 ℓ∞bound of the dynamic

• The analysis bounds the ℓ∞ magnitude of the transformed HMC dynamics with probability at least 1 − η under the stated assumptions.The within-step estimate includes the initial position, initial velocity, and a τM term.Source: §5.1 · pp. 28–29

5.2 Lipschitz constant of F

• For the transformed dynamics s = Ax, the force function F is (Lφ′τ + m2)-Lipschitz in ℓ∞ norm, where τ = ∥AA⊤∥∞→∞.Here Lφ′ is the Lipschitz constant of φ′.Source: §5.2 · p. 31

5.4 Main result

• Each iteration uses O(d log^3(1/δ)) evaluations of φ′ and matrix-vector multiplications by A and A⊤.The theorem’s iteration bound is also stated in terms of the condition number κ = M2/m2.Source: §5.4 · pp. 33, 36

A.1 Notation

• The appendix defines the p-to-q matrix norm, Wasserstein distance, and Cauchy estimates used in the analysis.It also records an exponential upper bound for a helper integral-equation solution.Source: §A.1 · p. 41, §A.2 · p. 42

A.2 Simple ODEs

• The appendix introduces helper lemmas for bounding solutions of simple integral equations.One stated bound is v(t) ≤ βe^γt for all t ≥ 0.Source: §A.2 · p. 42

A.3 Cauchy Estimates and Method of Majorants

• For an analytic vector field, the cited theorem gives a unique real analytic solution around t = 0 and a further bound under its stated conditions.The supplied passage does not specify the bound’s expression.Source: §A.3 · p. 43
• The majorant argument bounds derivatives of a multivariate ODE solution through a scalar ODE, then bounds that scalar solution directly.The paper notes that this makes explicit a bound often stated qualitatively, citing for an introduction to majorants.Source: §A.3 · pp. 42–43

B.1 Proof of general k-th order ODE

• The section’s proof of Theorem 2.4 derives an error bound and a bound on the number of iterations.The expressions for those bounds are not included in the supplied passages.Source: §B.1 · pp. 44, 47, 48
• The proof establishes that x*(t) is approximated by an element of Vd.The supplied passages do not state the approximation error explicitly.Source: §B.1 · p. 46

C Deferred Proof for Cauchy Estimates (Section A.3)

• Lemma A.9 compares derivatives of the multivariate solution with those of a scalar ODE under the stated analytic majorization assumptions.The proof starts from the second-order equation and differentiates it to obtain higher-derivative relations.Source: C Deferred Proof for Cauchy Estimates… · p. 48
• The comparison begins by bounding the initial first and second derivatives of the vector solution by the scalar solution’s corresponding derivatives.The proof then proceeds to higher derivatives.Source: C Deferred Proof for Cauchy Estimates… · p. 49

• Lemma D.1 assumes that f is holomorphic on U and that the closed disk of radius R centered at z0 lies inside U.Source: D Cauchy Estimates of Some Functions · p. 50
• The lemma denotes the disk boundary by γR and defines it as the circle centered at z0 with radius R.Source: D Cauchy Estimates of Some Functions · p. 50
• The lemma states its estimate for every integer n ≥ 0.Source: D Cauchy Estimates of Some Functions · p. 50

D.1 Logistic loss function

• The logistic loss derivative φ′(t) has Cauchy estimate M = 1 with radius r = 1.Source: §D.1 · p. 50
• The proof takes a real center z0 and sets the circle radius to 1, which bounds the imaginary component b between −1 and 1.Source: §D.1 · p. 50
• On that range, the proof uses cos b ∈ [0.54, 1] and the lower bound 1 + e^a + e^(2a) ≥ 1.Source: §D.1 · p. 50

D.2 Pseudo-Huber loss function

• The passage identifies Pseudo-Huber loss as a smooth alternative to Huber loss, whose applications include regression [CW15a], low-rank approximation [SWZ19], clustering [BFL16], and sparse recovery.Source: §D.2 · p. 51
• For fixed δ > 0, the Pseudo-Huber derivative φ′(x) has Cauchy estimate M = 1 with radius r = δ/2.Source: §D.2 · p. 51
• The proof applies Lemma D.1 to establish the stated estimate for all real z0.Source: §D.2 · p. 51