Numerical Analysis

Motivation

Matrix approximation replaces a large or complicated matrix with a random matrix that is easier to store or compute. Averaging independent realizations improves accuracy while retaining the structure of each realization. Examples in Tropp's monograph include sparse approximations and randomized matrix products. The practical question is how the number of samples, the size of one sample, and its second moment control the spectral norm error. Tropp, §6.2

The book proves a matrix Bernstein inequality in Theorem 6.1.1, then applies it to an empirical matrix average in Corollary 6.2.1. Theorem 6.1.1 is already posed as TroppMatrixConcentration.matrix_bernstein on Prove2Me and is currently Open. This mission records the application with its own hypotheses and constants. The capstone concerns the estimator that a reader would construct, rather than a sum of abstract centered summands. Tropp, §§6.1–6.2

Setting

Let $B$ be a fixed complex $d_1\times d_2$ matrix. A sample matrix $R$ has the same size, has expectation $B$, and has spectral norm at most $L$ almost surely. The equality $\mathbb ER=B$ makes one sample an unbiased estimate of $B$. The almost-sure bound limits the size of every sample except on a probability-zero set. The dimensions $d_1,d_2$ are positive integers, and $n$ is a positive sample count.

Take independent copies $R_1,\ldots,R_n$ with the distribution of $R$. Their matrix sampling estimator is

$$\overline R_n=\frac1n\sum_{k=1}^nR_k.$$

For any rectangular matrix $A$, $A^*$ denotes its conjugate transpose and $\|A\|$ its spectral norm. The per-sample second moment is

$$m_2(R)=\max\{\|\mathbb E(RR^*)\|,\|\mathbb E(R^*R)\|\}.$$

Both terms matter for a rectangular matrix: $RR^*$ and $R^*R$ act on different coordinate spaces. All expectations of matrices are entrywise expectations, represented by Bochner integrals in Lean. The almost-sure norm bound guarantees the moments used here are finite. Tropp, equations (6.2.4)–(6.2.6)

Formalization targets

Empirical matrix approximation

Corollary 6.2.1 gives an expectation estimate and a tail estimate for the same estimator. For every real $t\ge0$,

$$\mathbb E\|\overline R_n-B\|\le \sqrt{\frac{2m_2(R)\log(d_1+d_2)}{n}}+ \frac{2L\log(d_1+d_2)}{3n},$$ $$\mathbb P\{\|\overline R_n-B\|\ge t\}\le (d_1+d_2)\exp\left(\frac{-nt^2/2}{m_2(R)+2Lt/3}\right).$$

The goal includes both bounds. The expectation estimate controls average error; the tail estimate controls the chance of an error above a chosen threshold. The dimension enters only through $\log(d_1+d_2)$ in the expectation bound and the factor $d_1+d_2$ in the tail bound. These are the book's constants, including the factors $2$ and $3$ in the linear terms. Tropp, Corollary 6.2.1, p. 81

Supporting statements

The milestone list follows Corollary 6.1.2 for uncentered summands, the bound $\|n^{-1}(R_k-B)\|\le2L/n$, and the variance estimate $v(\overline R_n-B)\le m_2(R)/n$. Equation (6.2.7) is a companion: a stated lower bound on $n$ implies $\mathbb E\|\overline R_n-B\|\le2\varepsilon$ for $\varepsilon>0$. Tropp, pp. 77, 82

Significance

The result turns a bound on one random matrix into a quantitative guarantee for an empirical approximation. It lets a user choose a sample count for a target error and evaluate both average error and failure probability. The same framework is instantiated later in the monograph for matrix sparsification and multiplication; those applications depend on controlling $L$ and $m_2(R)$ for their particular sampling laws. Tropp, §§6.2–6.4

This is a known mathematical result proved in the book. Its Lean statement is an open formalization target. The mission adds a reusable empirical-mean definition and application statements built on the already published matrix spectral norm and rectangular second-moment definitions. A proof of the general matrix Bernstein theorem remains separate; the mission can draw on it when that platform theorem is proved.

Difficulty

The estimator error is a centered sum, but the sampling assumptions are given for an uncentered random matrix $R$. Bounding the centered summands changes the uniform bound from $L$ to $2L/n$. The variance parameter must be related to the second moment of one uncentered sample, with the correct $1/n$ scaling. Dropping the common-distribution condition would leave no relation between the samples' second moments and $m_2(R)$. Applying scalar concentration to entries would not directly give the stated spectral norm bound or its dimension dependence. Tropp, proof of Corollary 6.2.1, p. 82

Formalization scope

Matrices are Matrix (Fin d₁) (Fin d₂) ℂ. The spectral norm and rectSecondMoment are imported from tc's published TroppMatrixConcentration definitions. An arbitrary probability space carries a measurable reference matrix $R_0$ and a finite family $R_k$; iIndepFun expresses independence and IdentDistrib expresses equality in law with $R_0$. The estimator is the finite sum scaled by $n^{-1}$. No particular sampling scheme or deterministic choice of $R$ is fixed: the formal goal requires only the source's unbiasedness and boundedness conditions.

The Lean statements explicitly require $d_1,d_2\ge1$, $n>0$, measurability, and the almost-sure form of $\|R\|\le L$. These make the logarithm and inverse sample count meaningful and make the random matrix laws well behaved. For the uncentered Bernstein corollary, integrability of each $S_k$ is explicit because a bounded deviation from an unspecified mean alone does not establish that the mean exists. The goal also asserts integrability of its error norm, preventing a Bochner integral from silently taking Lean's default value for a nonintegrable function. The bound on $R$ supplies that integrability. When the denominator in a tail formula is zero, the displayed total real division gives a loose bound at the degenerate case; the source's probabilistic claim remains valid.

Contributions that formalize the centered-summand calculation, the variance estimate, or the transfer from the general matrix Bernstein theorem are within scope. The imported matrix and probability definitions are reusable in other matrix concentration missions. The local empirical-mean definition is reusable whenever a matrix estimator averages finitely many random samples.

Selected references

• Joel A. Tropp, An Introduction to Matrix Concentration Inequalities, Foundations and Trends in Machine Learning 8(1–2), 2015, arXiv:1501.01571v1, especially §§6.1–6.4.

Motivation

Solving a square system of nonlinear equations $\psi_i(x) = 0$, $i = 1, \dots, n$, $x \in \mathbb{R}^n$, is one of the oldest tasks of numerical analysis. It appears whenever optimality conditions, equilibrium conditions or discretized differential equations have to be solved. N. Z. Shor's Minimization Methods for Non-Differentiable Functions (Springer 1985, translated by K. C. Kiwiel and A. Ruszczyński) treats such a system as the nonsmooth minimization problem

$$\min_{x \in E_n} f(x), \qquad f(x) = \max_{1 \le i \le n} |\psi_i(x)|,$$

whose optimal value is $0$ exactly when the system is consistent. Section 3.5 of the book applies the subgradient method with space dilation to this problem. At each iteration that method needs the gradient of one of the $\psi_i$, not all $n$ of them. Taking the dilation coefficient to its extreme value ($\beta = 0$) turns the method into a gradient orthogonalization method. The book shows that this method converges at a quadratic rate near a regular solution, like Newton's method, although it never forms or factorizes the Jacobian.

This mission is the eighth in a series formalizing Shor's book. It covers Sections 3.5–3.6 (printed pp. 62–71).

Setting

$E_n$ is the $n$-dimensional Euclidean space with inner product $(x, y)$ and norm $\|x\|$. The data are functions $\psi_1, \dots, \psi_n : E_n \to \mathbb{R}$ with gradients $g_{\psi_i}(x)$, and the max-residual function is $f(x) = \max_i |\psi_i(x)|$ (3.24). The Jacobian is the matrix $J(x) = \{\partial\psi_i/\partial t_j\}_{i,j=1}^n$ of partial derivatives with respect to the coordinates $x = \{t_1, \dots, t_n\}$.

An almost-gradient of a function $f$ at $x_0$ is an accumulation point of a sequence of gradients $\nabla f(x_k)$, where $x_k \to x_0$ and $f$ is differentiable at every $x_k$ (p. 18).

The regular case (p. 63) is the situation in which $f^* = \min f = 0$ is attained at a point $x^*$, the $\psi_i$ are continuously differentiable near $x^*$, and $J(x^*)$ is nonsingular.

The gradient orthogonalization method (3.27) runs in stages of $n$ steps. A stage starting at $x_0$ produces $x_1, \dots, x_n$ as follows. At step $k+1$, $0 \le k < n$, the gradient $g_{\psi_{k+1}}(x_k)$ is projected on the subspace orthogonal to the gradients $g_{\psi_1}(x_0), \dots, g_{\psi_k}(x_{k-1})$ used earlier in the stage. Call the result $\varphi_{k+1}$. The step is

$$x_{k+1} = x_k - \frac{\psi_{k+1}(x_k)}{\|\varphi_{k+1}\|^2}\, \varphi_{k+1}.$$

For $k = 0$ there is nothing to project, so $\varphi_1 = g_{\psi_1}(x_0)$. The next stage starts at $x_n$.

Formalization targets

Goal: Theorem 3.9, one stage squares the error

Let $x^* = 0$ solve the system, let the $\psi_i$ be continuously differentiable with gradients Lipschitz with constant $L$ on a ball $S_\delta = \{x : \|x\| < \delta\}$ (3.28), and let $g_{\psi_1}(0), \dots, g_{\psi_n}(0)$ be linearly independent. Then there exist $\varepsilon > 0$ and $c > 0$ such that every stage with $\|x_0\| \le \varepsilon$ satisfies

$$\|x_n\| \le c\, \|x_0\|^2 .$$

Milestones inside the proof

The proof of Theorem 3.9 passes through these displayed results, stated here in attack order:

• Eq. (3.29): near the solution, the orthogonalized gradients satisfy $\|\varphi_{k+1}\| > b > 0$, with $b$ uniform over stages.
• Eq. (3.33): after step $k+1$, $|\psi_{k+1}(x_{k+1})| \le (L/b^2)\, \psi_{k+1}^2(x_k)$.
• Lemma 3.1: a step of length $h$ in a unit direction orthogonal to $g_{\psi_i}(x_1)$ changes $\psi_i$ by at most $L h (h + \|x_1 - x_2\|)$.
• Eqs. (3.36)–(3.37): within a stage, $\max_k \|x_k\| \le c_1 \|x_0\|$, and at the end point $f(x_n) \ge c_2 \|x_n\|$.

Companion result: Theorem 3.8

In the regular case, for every $\delta > 0$ there is a neighborhood of $x^*$ in which every almost-gradient satisfies

$$(1-\delta) f(x) \le (g_f(x), x - x^*) \le (1+\delta) f(x).$$

Significance

Theorem 3.9 says that a method which evaluates one equation and one gradient per step, and never solves a linear system, has the local quadratic rate of Newton's method on a regular system. The stages combine: once one stage starts close enough to the solution, the errors of later stages satisfy $\|x_0^{(r+1)}\| \le c\|x_0^{(r)}\|^2$. Theorem 3.8 describes the geometry of $f$ near a regular solution. The ratio of $(g_f(x), x - x^*)$ to $f(x)$ tends to $1$, which is what allows the dilation parameters of the book's Theorem 3.3 to approach their extreme values near the solution. Both results connect the nonsmooth space-dilation methods of Chapter 3 to classical Newton-type theory.

All results of this mission are proved in the book. None of them has, to our knowledge, a machine-checked proof: the platform has local quadratic convergence results for cubic-regularized Newton, proximal Newton and semismooth Newton, but not for this method. The work is formalizing Shor's proof. This includes building the Gram-determinant argument behind (3.29) and the second-order estimates, which are reusable for other projection-based solvers of Kaczmarz type.

Difficulty

A naive argument analyses each step as a Newton step for one equation and stops there. That shows each equation $\psi_{k+1}$ is small right after its own step, but the theorem needs all equations to be small at the end of the stage. Later steps move the point and can undo the progress on earlier equations. The step directions are orthogonal to the earlier gradients, but those gradients were taken at earlier points, not at the current one. Controlling this drift needs Lemma 3.1 together with the bound (3.36), which says that the whole stage stays within a multiple of $\|x_0\|$.

A second difficulty is that the method is only defined near the solution. The step divides by $\|\varphi_{k+1}\|^2$, and a uniform positive lower bound on these norms comes from the nonvanishing of a Gram determinant evaluated at $n$ different points. The final passage from small residuals to a small error needs the growth bound (3.37), which uses the nonsingularity of the Jacobian once more.

Formalization scope

• Representation. $E_n$ is EuclideanSpace ℝ (Fin n). The equations are a family ψ : Fin n → E_n → ℝ indexed from $0$, so ψ k is the book's $\psi_{k+1}$. Gradients are Mathlib's gradient, and $S_\delta$ is the open ball Metric.ball 0 δ. Continuous differentiability is ContDiff ℝ 1; in Theorem 3.8 it is ContDiffOn ℝ 1 on a neighborhood of $x^*$.
• The stage. A stage is a predicate IsOrthStage ψ x on a sequence $x : \mathbb{N} \to E_n$, which fixes $x_1, \dots, x_n$ from $x_0$. The projection is starProjection onto the orthogonal complement of the span of the earlier gradients.
• Corrected misprint. The printed Step 1, (3.27a), divides by $\|g_{\psi_1}(x_0)\|$ rather than its square. The formalization uses the square, which is what (3.27b) and the proof's (3.31)–(3.33) require at $k = 0$.
• Division by zero. The book gives no rule for $\varphi_{k+1} = 0$; Lean's $t/0 = 0$ would leave the point unchanged. All statements apply the method only near the solution, where (3.29) bounds $\|\varphi_{k+1}\|$ away from zero.
• Quantifier order. In Theorem 3.9 and in (3.29), (3.36)–(3.37), the constants ($\varepsilon$, $c$, $\delta'$, $b$, $\varepsilon_0$, $c_1$, $c_2$) are existentials placed before the stage. They depend only on the data $\psi, \delta, L$. A statement with $c$ chosen after $x_0$ would be trivial, with $c = \|x_n\|/\|x_0\|^2$, and is ruled out.
• Nonsingularity. In Theorem 3.9 it is the linear independence of $g_{\psi_i}(0)$, the book's primary hypothesis. In Theorem 3.8 it is $\det J(x^*) \ne 0$, with $J$ defined by coordinate partial derivatives.
• Almost-gradients. They are quantified universally: Theorem 3.8 holds for every almost-gradient, not for one chosen one.

Not formalized: Theorem 3.10 (pp. 70–71), the $n$-step quadratic convergence of the $\beta = 0$ r-algorithm with resetting. Its directional minimization rule is not pinned down on the page. The algorithm says "determined by minimizing", while the proof uses "the smallest positive root" of a stationarity equation. A contribution that fixes this rule and states the theorem is welcome as a follow-up.

Contributions of any of the five milestones are welcome independently. The Gram-determinant lower bound (3.29) and Lemma 3.1 are general facts about orthogonalization and Lipschitz gradients, and are useful beyond this mission.

Selected references

• N. Z. Shor, Minimization Methods for Non-Differentiable Functions, Springer Series in Computational Mathematics 3, Springer, 1985, §§3.5–3.6, pp. 62–71. https://doi.org/10.1007/978-3-642-82118-9
• J. M. Ortega and W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, 1970. https://doi.org/10.1137/1.9780898719468

Motivation

Many proposed quantum optimization and machine-learning methods (variational quantum eigensolvers, QAOA, quantum classifiers) tune the real parameters of a circuit by gradient descent on an objective $f:\mathbb R^d\to\mathbb R$ that is only accessible as a black box. Classically, a finite-difference estimate of $\nabla f$ costs on the order of $d$ evaluations of $f$. Quantum algorithms can access $f$ coherently through a phase oracle $O_f:|x\rangle\mapsto e^{if(x)}|x\rangle$, and the question is how many such queries are needed to learn $\nabla f(x)$ to accuracy $\varepsilon$ in the sup norm.

Timeline.

• 2005: Jordan (Phys. Rev. Lett. 95, 050501) gave an algorithm that reads off the gradient of a nearly affine function from a single oracle call on a uniform superposition over a grid, followed by an inverse quantum Fourier transform. The guarantee requires $f$ to be very close to affine on the grid.
• 2017–2019: Gilyén, Arunachalam and Wiebe (arXiv:1711.00465v3; SODA 2019) ran Jordan's algorithm on a high-degree central difference of $f$ rather than on $f$. For functions whose derivatives satisfy a Gevrey-type bound, they obtained $\widetilde O(c\sqrt d/\varepsilon)$ phase queries (Theorem 25), and proved a matching $\Omega(\sqrt d/\varepsilon)$ lower bound for a class of smooth functions.
• 2024: Nannicini's lecture notes (arXiv:2408.07086) present Jordan's algorithm in this form. Their Lean formalization on the platform contains the grid, its Fourier transform and the algorithm's state, which this mission reuses.

Setting

A point of $\mathbb R^d$ is $x=(x_1,\dots,x_d)$, with $\|x\|_\infty=\max_i|x_i|$. The paper's $\mathbb N$ is $\{1,2,\dots\}$.

Grid and Jordan's algorithm. For $n$ qubits put $N=2^n$. A bit string $b$ encoding the integer $j$ labels the grid point $x^{(b)}=j/N-1/2+1/(2N)$, and $G_n$ is the set of these points (Definition 17). Algorithm 2 prepares $d$ registers of $n$ qubits in uniform superposition over $G_n^d$. It applies the phase oracle $|x\rangle\mapsto e^{2\pi i N f(x)}|x\rangle$ once, applies the inverse grid Fourier transform to each register, and measures the outcome $(k_1,\dots,k_d)\in G_n^d$.

Index sequences and analyticity. For $\alpha\in[d]^k$, $\partial_\alpha f$ is the $k$-th partial derivative along the coordinates $\alpha_1,\dots,\alpha_k$, and $y^\alpha=\prod_j y_{\alpha_j}$ (Definition 4). The function $f$ is analytic around $x$ if

$$f(x+y)=\sum_{k=0}^\infty\sum_{\alpha\in[d]^k}\frac{y^\alpha\,\partial_\alpha f(x)}{k!}\quad\text{for all }y\in\mathbb R^d$$

(Definition 5, there at $x=0$).

Central differences. The central difference coefficients are

$$a_\ell^{(2m)}=\frac{(-1)^{\ell-1}}{\ell}\frac{\binom m{|\ell|}}{\binom{m+|\ell|}{|\ell|}}\quad(0<|\ell|\le m),\qquad a_0^{(2m)}=0.$$

The degree-$2m$ central difference of $f$ is $f_{(2m)}(y)=\sum_{\ell=-m}^m a_\ell^{(2m)}f(\ell y)$, an approximation of $\nabla f(0)\cdot y$ (Definition 7).

Queries. A fractional query $O_f^s:|x\rangle\mapsto e^{isf(x)}|x\rangle$ with $|s|\le1$ counts as one query (Definitions 10–11). The algorithm of Theorem 25 has parameters $m,n_\varepsilon,n_M,T$ and $r>0$. It runs Algorithm 2 with $n=n_\varepsilon+n_M$ qubits per register on the phase $2\pi\,2^{n_\varepsilon} g_{(2m)}(rz)$, where $g(y)=f(x+y)$. It reads each outcome $k$ as $(2^{n_M}/r)\,k$, repeats $T$ times and outputs the coordinate-wise median. One run costs $\sum_\ell\lceil 2\pi 2^{n_\varepsilon}|a^{(2m)}_\ell|\rceil$ queries to $O_f$.

Formalization targets

Goal: Theorem 25

There is an absolute constant $C$ such that for every $d\ge1$, $0<\varepsilon\le c$ and $\delta>0$ some parameter choice, made before $x$ and $f$, uses at most

$$C\,\frac{c\sqrt d}{\varepsilon}\,\log\Big(2+\frac{c\sqrt d}{\varepsilon}\Big)\log\Big(2+\log\Big(2+\frac{c\sqrt d}{\varepsilon}\Big)\Big)\log\Big(2+\frac d\delta\Big)$$

queries. For every $x$ and every smooth $f$ that is analytic around $x$ with $|\partial_\alpha f(x)|\le c^k k^{k/2}$ ($k\ge1$), it outputs $\tilde\nabla f(x)$ with $\|\nabla f(x)-\tilde\nabla f(x)\|_\infty\le\varepsilon$ with probability at least $1-\delta$.

Milestones

• Lemma 34, $\sum_\ell|a^{(2m)}_\ell\ell^{k+1}|\le6e^{-7m/6}m^{k+3/2}$ for $k\ge2m$.
• Lemma 22: (23), $|f'(0)\delta-f_{(2m)}(\delta)|\le e^{-m/2}\|f^{(2m+1)}\|_\infty|\delta|^{2m+1}$; and (24), $\sum_{\ell=0}^m|a_\ell^{(2m)}|<\sum_{\ell=1}^m 1/\ell\le\ln m+1$.
• The moment identity (proof of Theorem 24): $\sum_{\ell=-m}^m a_\ell^{(2m)}\ell^k$ is $1$ for $k=1$ and $0$ for the other $k\le2m$.
• Lemma 36, a tail bound $\mathbb P[|\sum_\alpha H_\alpha x^\alpha|\ge\sqrt2(r\sqrt{dk/2})^k]\le r^{-2k}$ for i.i.d. symmetric $x_i\in[-1/2,1/2]$.
• Theorem 24: $|\nabla f(0)y-f_{(2m)}(y)|\le\sum_{k\ge2m+1}(8Rcm\sqrt d)^k$ for all but a $1/1000$ fraction of $y\in R\cdot G_n^d$.
• The choice of $r$ (proof of Theorem 25): with $r^{-1}=9cm\sqrt d\,(81\cdot8\cdot42\pi cm\sqrt d/\varepsilon)^{1/(2m)}$,

$$\sum_{k=2m+1}^\infty(8rcm\sqrt d)^k\le\frac{\varepsilon r}{8\cdot42\pi}.$$

• Lemma 20: if $|f(x)-g\cdot x-c|\le1/(42\pi N)$ off a $1/1000$ fraction of $G_n^d$ and $\|g\|_\infty\le1/3$, then $\Pr[|k_i-g_i|>4/N]\le1/3$.
• Theorem 21 (query part): $O(\log(d/\rho))$ runs and a median give $\Pr[\|\tilde g-g\|_\infty>\varepsilon]\le\rho$.
• Display (25): $\sum_\ell\lceil|a^{(2m)}_\ell|S\rceil\le2m+S(2\log m+2)$. The page prints the middle term $S\sum_\ell a^{(2m)}_\ell$ without absolute values; the formal statement uses $|a^{(2m)}_\ell|$.

Significance

The result. Theorem 25 gives a quantum gradient algorithm whose query count grows like $\sqrt d$ in the dimension and nearly linearly in $1/\varepsilon$. Jordan's original analysis gives $\sqrt d/\varepsilon^2$, and classical finite differences need on the order of $d$ evaluations. Section 6 of the paper shows that $\sqrt d/\varepsilon$ cannot be beaten for a class of smooth functions, so Theorem 25 is optimal up to logarithmic factors. Section 5.5 shows that probability oracles built from parametrized circuits satisfy the derivative bound, so the theorem applies to gradient-based training of variational circuits.

Formalizing it. The theorems are proved in the paper; this mission is not aware of a machine-checked proof of any of them. A formal proof would check the explicit constants (1/1000, $42\pi$, $8Rcm\sqrt d$, $6e^{-7m/6}$) and the accounting of fractional queries, which are spread across the main text and Appendix A. It would also produce reusable libraries: exactness and coefficient bounds for central differences, concentration of random multilinear forms, and a phase-estimation-style analysis of Jordan's algorithm on the published Lean model.

Difficulty

Jordan's algorithm reads off a gradient only if the phase is almost exactly affine on the whole grid. A truncated Taylor expansion of $f$ makes the error of order-$k$ terms grow like $d^{k}$ on the cube, which destroys the $\sqrt d$ scaling. The approximation is only good on most of the grid, not all of it. The guarantee therefore has to tolerate a small fraction of bad grid points inside the quantum algorithm (Lemma 20) and has to control random multilinear forms (Lemma 36). The analysis also has to keep explicit, $m$-uniform control of the central difference coefficients, whose sizes enter both the approximation error and the query count.

Formalization scope

• Representation. $\mathbb R^d$ is Fin d → ℝ, whose norm is $\|\cdot\|_\infty$. $\partial_\alpha f$ is iteratedFDeriv along coordinate vectors, and functions are assumed $C^\infty$ (ContDiff ℝ ((⊤ : ℕ∞) : WithTop ℕ∞)). "Analytic" is the Taylor representation (5) at the relevant point (IsTaylorRep), not Mathlib's local AnalyticAt.
• Fractions and probabilities. "All but a $1/1000$ fraction of $G_n^d$" means the number of bad grid labels is at most $N^d/1000$. Probabilities are those of the published Algorithm 2 state QAlgorithms.gradientAlgState. Independent runs multiply.
• Indices and constants. $\mathbb N$ starts at $1$: derivative bounds are for $k\ge1$, and $m\ge1$. $O(\cdot)$ and $\widetilde O(\cdot)$ are absolute constants quantified before every parameter, with the polylog factor of display (26) written out.
• Added hypotheses.
  • Theorem 24 assumes $8Rcm\sqrt d<1$, the condition under which its printed series is finite.
  • Lemma 22 uses a bound $B$ on $|f^{(2m+1)}|$ in place of the supremum.
  • Theorem 21 requires its register size $n=n_\varepsilon+n_M$ to be a natural number.
• Query count. The count is a definition: $T\sum_\ell\lceil2\pi2^{n_\varepsilon}|a_\ell^{(2m)}|\rceil$.
• Not formalized. Gate complexity (no circuit model), the probability-to-phase oracle conversions, the lower bounds and Theorem 23.
• Ruled out. None of the following counts as Theorem 25: a goal with an existential over arbitrary functions of $f$ (which could output $\nabla f(x)$ with no queries), a query count not tied to the algorithm that is run, or a statement for one fixed $d$.

Contributions are welcome at every level. Most useful are general facts about central difference coefficients (moment identities, $\sum|a_\ell|$, the binomial ratio bounds), concentration of multilinear forms in independent symmetric variables, and the single-run analysis of gradientAlgState. These are reusable beyond this mission.

Selected references

• A. Gilyén, S. Arunachalam, N. Wiebe, Optimizing quantum optimization algorithms via faster quantum gradient computation, arXiv:1711.00465v3 (2018); Proceedings of SODA 2019. https://arxiv.org/abs/1711.00465
• S. P. Jordan, Fast quantum algorithm for numerical gradient estimation, Physical Review Letters 95(5):050501, 2005. https://arxiv.org/abs/quant-ph/0405146
• J. Li, General explicit difference formulas for numerical differentiation, Journal of Computational and Applied Mathematics 183(1):29–52, 2005.
• G. Nannicini, Quantum algorithms for optimizers, arXiv:2408.07086, 2024. https://arxiv.org/abs/2408.07086
• M. A. Nielsen, I. L. Chuang, Quantum Computation and Quantum Information, Cambridge University Press, 2002 (phase estimation, §5.2).

Motivation

Many models in statistics, machine learning and engineering are nonconvex, nonsmooth and coupled by constraints: a sparsity penalty such as SCAD, MCP or capped-$\ell_1$ on one block of variables, a smooth loss coupling the blocks, and a nonlinear equality constraint tying them together. Augmented-Lagrangian methods and ADMM-type splittings are the workhorses for such problems, but most convergence theory for them assumes linear constraints or convexity. Zhu, Zhao and Zhang (arXiv:2007.12219v1, published in Mathematics of Operations Research, 2024) propose a first-order primal-dual method, NAPP-AL (nonconvex auxiliary problem principle with augmented Lagrangian), for a nonconvex problem with a nonlinear coupling constraint, and prove that it reaches approximate stationarity at a sublinear rate. This mission targets that rate result and the chain of lemmas behind it.

Setting

Fix dimensions $n,d,m$. Problem (P) is

$$\min_{u\in U,\ v\in\mathbb R^d} F(u,v):=G(u,v)+J(u)+H(v)\quad\text{s.t.}\quad\Theta(u)+Bv=0,\qquad\Theta:=\Omega+\Phi,$$

where $G$ and $H$ are differentiable with Lipschitz gradients, $J$ is lower semicontinuous and possibly nonconvex and nonsmooth, $U\subseteq\mathbb R^n$ is closed and convex, $\Theta:\mathbb R^n\to\mathbb R^m$ is Lipschitz with constant $L^0_\Theta$, each component of $\Omega$ has a Lipschitz gradient (constants $L_{\Omega_j}$, $L_\Omega=\sum_jL_{\Omega_j}$), and $B$ is an $m\times d$ matrix with $\operatorname{Im}\Theta\subseteq\operatorname{Im}B$. These are Assumption 1 (H1)–(H7); (H7) asks $F$ to be bounded below and coercive on the feasible set $\{(u,v):\Theta(u)+Bv=0\}$.

The augmented Lagrangian of (P) at $w=(u,v,p)$ is

$$\mathcal L_\gamma(w)=F(u,v)+\langle p,\Theta(u)+Bv\rangle+\frac\gamma2\|\Theta(u)+Bv\|^2 .$$

Its stationarity set $\partial\mathcal L_\gamma(w)$ (p. 6) collects the triples $(\xi_u,\xi_v,\xi_p)$ with, for $y=p+\gamma(\Theta(u)+Bv)$,

$$\xi_u\in\nabla_uG(u,v)+\partial J(u)+N_U(u)+(\nabla\Theta(u))^\top y,\quad\xi_v=\nabla_vG(u,v)+\nabla H(v)+B^\top y,\quad\xi_p=\Theta(u)+Bv,$$

where $\partial J$ is the limiting subdifferential and $N_U$ the normal cone of $U$. A stationary point is a $w$ with $0\in\partial\mathcal L_\gamma(w)$, and $\operatorname{dist}(0,S)=\inf_{\xi\in S}\|\xi\|$, with $\operatorname{dist}(0,\emptyset)=+\infty$.

NAPP-AL takes a $\beta$-strongly convex $K$ with Lipschitz gradient (Bregman distance $D(u,u')=K(u)-K(u')-\langle\nabla K(u'),u-u'\rangle$), a penalty $\gamma>\frac{\sqrt{57}+1}{2\lambda_{\min}(B^\top B)}(L_G+L_H)$ (Assumption 2) and $0<\sigma<1$. With $q^k=p^k+\gamma(\Theta(u^k)+Bv^k)$ it sets $u^{k+1}$ to a global minimizer over $U$ of the linearized model

$$\langle\nabla_uG(u^k,v^k),u\rangle+J(u)+\langle q^k,\nabla\Omega(u^k)u+\Phi(u)\rangle+\tfrac1{\epsilon^k}D(u,u^k),$$

$v^{k+1}$ to the minimizer of a strongly convex quadratic, and $p^{k+1}=p^k+\gamma(\Theta(u^{k+1})+Bv^{k+1})$. The step $\epsilon^k$ lies in $[\sigma\delta_k,\delta_k]$, with $\delta_k$ given by the explicit formula (8), which depends on $\|q^k\|$. The analysis runs through the potential

$$\Lambda^k=\mathcal L_\gamma(w^k)+c_1\|u^{k-1}-u^k\|^2+c_2\|v^{k-1}-v^k\|^2+\tfrac1{2\gamma}\|p^{k-1}-p^k\|^2 \tag{18}$$

with explicit constants $c_1,c_2$.

Formalization targets

Goal: Proposition 4.1 (p. 18)

Under Assumptions 1 and 2, with $\delta_k$ set by (8),

$$\min_{1\le j\le k}\operatorname{dist}\big(0,\partial\mathcal L_\gamma(w^j)\big)=o\big(1/\sqrt k\big).$$

Milestones

1. Proposition 3.1 (p. 9): $\operatorname{dist}(0,\partial\mathcal L_\gamma(w^{k+1}))\le\sqrt3\,h(u^k,p^k)\|w^k-w^{k+1}\|$ with the paper's explicit $h$ (see the scope section for the factor $\sqrt3$).
2. Lemma 4.1 (p. 11): the dual step $\|p^k-p^{k+1}\|$ is controlled by the primal steps.
3. Lemma 4.2 (pp. 11–12): one-step change estimates for $\mathcal L_\gamma$.
4. Lemma 4.3 (p. 15): sufficient decrease $\Lambda^{k+1}-\Lambda^k\le-c_3\|w^k-w^{k+1}\|^2$ with $c_3=\min\{\frac12,c_4,\frac1{3\gamma}\}>0$.
5. Theorem 4.1 (p. 16): $\Lambda^k$ converges, the iterates are bounded with vanishing steps, the residuals vanish, and every cluster point is stationary.

Significance

The result gives NAPP-AL an iteration complexity of $O(1/\varepsilon^2)$ to reach an $\varepsilon$-stationary point (Remark 4.1), for problems with a nonlinear, nonconvex coupling constraint and a nonsmooth nonconvex regularizer, using only gradient evaluations and one proximal-type subproblem per iteration. It is also the base layer for the paper's later results: linear convergence under an error bound (Theorem 5.1) and the sufficient conditions for that error bound (Propositions 6.1–6.2) all start from Theorem 4.1 and Lemma 4.3.

The result is proved in the paper; nothing here is formalized elsewhere. A formal development adds a machine-checked Lyapunov analysis of an augmented-Lagrangian method with nonlinear constraints, with every constant explicit. It also records a correction: Proposition 3.1 as printed is false by a factor of at most $\sqrt3$, and the corrected bound is what the paper's proof establishes.

Difficulty

The difficulty is the dual variable. The multiplier update has no step size of its own, so the augmented Lagrangian itself need not decrease along the iterates, and the obvious descent argument for $\mathcal L_\gamma$ fails; the step $\delta_k$ of (8) depends on the current multiplier estimate $\|q^k\|$, so even the step sizes are not uniformly bounded a priori. Boundedness of the iterates is not a hypothesis, and the iterates are infeasible, so coercivity of $F$ on the feasible set does not apply to them directly. The rate needs the residual bound of Proposition 3.1 with a constant uniform in $k$, while $h(u^k,p^k)$ depends on $\|q^k\|$ and $\delta_k$. Passing to cluster points involves the limiting subdifferential of a nonconvex, nonsmooth $J$, which is not closed under arbitrary convergence of its arguments.

Formalization scope

• Spaces and norms. Points are in EuclideanSpace ℝ (Fin k). The norm of a triple $(u,v,p)$ is the Euclidean norm of the stacked vector, defined explicitly; Mathlib's sup norm on products is never used. $B$ is a continuous linear map with its operator (spectral) norm and adjoint; $\lambda_{\min}(B^\top B)$ is the Rayleigh minimum over the unit sphere.
• Residual. $\operatorname{dist}(0,\cdot)$ takes values in $[0,+\infty]$ and equals $+\infty$ on the empty set, as on p. 5. Metric.infDist (zero on $\emptyset$) is not used; with it, every residual bound would hold whenever $\partial\mathcal L_\gamma$ is empty. The goal is stated as $\sqrt k\cdot\inf_{1\le j\le k}\operatorname{dist}(0,\partial\mathcal L_\gamma(w^j))\to0$ in $[0,+\infty]$, not as a big-O bound.
• Stationarity. $\partial\mathcal L_\gamma$ is the p. 6 display with the limiting subdifferential of $J$ (the published NonconvexSplitting.Shared.LimitingSubdiff) and the normal cone of $U$, which is empty off $U$, so that stationarity forces $u\in U$.
• Algorithm. NAPP-AL is a relation on given sequences: $u^{k+1}$ is any global minimizer of the $u$-subproblem over $U$, whose existence is assumed by the paper (Remark 3.2), not asserted. The starting point is free.
• Standing readings of Assumption 1, disclosed in the items: $G,H,J$ are real valued with gradients Lipschitz on the whole space; "$J$ locally Lipschitz on $\operatorname{dom}J\cap U$" is read as Lipschitz on a neighbourhood in $\mathbb R^n$ of each point of $U$ (not Mathlib's relative LocallyLipschitzOn); $\lambda_{\min}(B^\top B)>0$ (full column rank, invoked in the proof of Lemma 4.2); (H7) is coercivity on the feasible set only.
• Indices. $\Lambda^k$ and the lemmas involving $w^{k-1}$ are stated for $k\ge1$.
• Printed slips. Proposition 3.1 carries the factor $\sqrt3$ (the printed bound fails on a one-dimensional instance with ratio about $1.108$). Theorem 4.1(iv) prints "$\partial\mathcal L_\gamma(u^k)$"; it is stated at $w^{k+1}$, as Proposition 3.1 and the later sections use it.

Contributions welcome: proofs of the milestones in any order, general facts about limiting subdifferentials and normal cones of convex sets (sum rule with a locally Lipschitz function, closedness under attentive convergence), and the "non-increasing summable sequence is $o(1/k)$" lemma. These are reusable beyond this mission. The source is the preprint arXiv:2007.12219v1, not the published version; all labels and constants are those of v1.

Selected references

• D. Zhu, L. Zhao, S. Zhang, A First-Order Primal-Dual Method for Nonconvex Constrained Optimization Based on the Augmented Lagrangian, arXiv:2007.12219v1, 2020; Mathematics of Operations Research, 2024. https://arxiv.org/abs/2007.12219v1
• B. S. Mordukhovich, Variational Analysis and Generalized Differentiation I, Springer, 2006. https://doi.org/10.1007/3-540-31247-1
• R. T. Rockafellar, R. J.-B. Wets, Variational Analysis, Springer, 1998. https://doi.org/10.1007/978-3-642-02431-3
• G. Cohen, Auxiliary problem principle and decomposition of optimization problems, Journal of Optimization Theory and Applications 32 (1980) 277–305. https://doi.org/10.1007/BF00934554

Motivation

A quantum algorithm for a system of linear equations must turn operations it can implement into an approximation of the inverse of the system matrix. For a Hermitian matrix, time evolution $e^{-iAt}$ is such an operation. Childs, Kothari and Somma use a Fourier expansion of the scalar reciprocal $1/x$ to express an approximation to $A^{-1}$ as a finite linear combination of these evolutions. Their paper reports an exponential improvement in the dependence on the requested precision over earlier linear-systems methods. The scalar approximation in Lemma 11 is the mathematical input to their Fourier algorithm and to the query and gate analyses in Theorem 3. This mission isolates that approximation and the Gaussian estimates stated around it. Theorem 3's complexity claims themselves are outside the present formalization. Childs, Kothari and Somma, §3, pp. 10–14.

Setting

Fix a condition number $\kappa\ge1$. The scalar spectral domain is

$$D_\kappa=[-1,-1/\kappa]\cup[1/\kappa,1].$$

It contains the real eigenvalues permitted after the paper's normalization and excludes zero, where the reciprocal is undefined. For an accuracy $\varepsilon>0$, two functions are $\varepsilon$-close on $D_\kappa$ when their difference has magnitude at most $\varepsilon$ at every point of that domain. A complex-valued approximant is compared with $1/x$ as a complex number. The paper gives this domain and closeness convention before the Fourier construction. Childs, Kothari and Somma, pp. 9–10.

The intermediate function $g_{y_J,z_K}(x)$ is a truncated Gaussian Fourier integral. It integrates $z e^{-z^2/2}e^{-ixyz}$ over $-z_K\le z\le z_K$ and then over $0\le y\le y_J$, multiplying by $i/\sqrt{2\pi}$. The finite approximant $h_{J,K,\Delta_y,\Delta_z}(x)$ is a Gaussian Fourier double sum with the same normalization and phase. Its samples are $y_j=j\Delta_y$ for $j=0,\ldots,J-1$ and $z_k=k\Delta_z$ for $k=-K,\ldots,K$; each term carries both mesh widths. Thus $J,K$ count the retained samples, while $\Delta_y,\Delta_z$ control their spacing. The sign in the phase and the factor $i$ are part of the definitions. Childs, Kothari and Somma, (18), p. 10, and (21), p. 11.

Formalization targets

Truncation

Lemma 12 says that for cutoffs of scale

$$y_J=\Theta\!\bigl(\kappa\sqrt{\log(\kappa/\varepsilon)}\bigr), \qquad z_K=\Theta\!\bigl(\sqrt{\log(\kappa/\varepsilon)}\bigr),$$

the function $g_{y_J,z_K}$ is $\varepsilon$-close to $1/x$ on $D_\kappa$. The proposal includes the paper's explicit error bound (27) and its Gaussian tail bound (28), as well as the Gaussian Poisson-summation identity of Lemma 13. These are separate reusable targets, not replacements for the approximation theorem. Childs, Kothari and Somma, Lemmas 12–13 and (27)–(28), pp. 11–12.

Finite Fourier sum

The goal is Lemma 11:

$$\left|h_{J,K,\Delta_y,\Delta_z}(x)-\frac1x\right|\le\varepsilon \quad\text{for every }x\in D_\kappa,$$

where

$$J=\Theta\!\left(\frac\kappa\varepsilon\log\frac\kappa\varepsilon\right),\quad K=\Theta\!\left(\kappa\log\frac\kappa\varepsilon\right),\quad \Delta_y=\Theta\!\left(\frac\varepsilon{\sqrt{\log(\kappa/\varepsilon)}}\right),\quad \Delta_z=\Theta\!\left(\frac1{\kappa\sqrt{\log(\kappa/\varepsilon)}}\right).$$

The $\Theta$ statements impose upper and lower bounds with absolute constants, uniformly across $\kappa$ and $\varepsilon$. The requested approximation error is exactly $\varepsilon$, matching the lemma's claim. The four scales specify both the size of the finite expansion and its mesh. Childs, Kothari and Somma, Lemma 11, p. 10.

Significance

The finite sum provides the paper's Fourier route from scalar reciprocal approximation to a linear combination of time-evolution operators. Its sample counts and spacings are the quantities that enter the algorithm's implementation costs. A result that only asserted approximation after arbitrarily many samples would not express the bound used there. The truncated integral and Poisson identity also remain useful for other Gaussian quadrature and Fourier-sampling arguments. Childs, Kothari and Somma, §§3.1–3.2, pp. 11–14.

The paper proves these claims on paper. This proposal asks for machine-checked Lean proofs of the scalar statements, including the complex normalization, the finite and infinite sum conventions, and the uniform scale bounds. It does not assert that the quantum algorithm, its Hamiltonian simulation subroutines, or its query and gate complexity theorems have been formalized.

Difficulty

The infinite Gaussian Fourier representation alone does not specify a finite implementable expansion. Restricting both integrals introduces two truncation errors, while replacing them by samples adds discretization and periodic-aliasing effects. These errors have to be controlled simultaneously with $J$, $K$, $\Delta_y$, and $\Delta_z$ at the scales stated in Lemma 11. Merely taking the mesh arbitrarily fine gives approximation without the theorem's finite-size bounds. The proof also uses a complex Fourier phase although its target is real, and its infinite sums must have their ordinary convergent meaning. Childs, Kothari and Somma, proof of Lemma 11, pp. 13–14.

Formalization scope

The Lean development uses $\mathbb R$ for $x,\kappa,\varepsilon$ and the cutoff and mesh parameters, $\mathbb C$ for $g$ and $h$, and the complex norm for closeness. Finite sums use exactly $0\le j<J$ and $-K\le k\le K$; the integral for $g$ is iterated and has nonnegative upper cutoffs. The optional identity (20) is also an iterated integral, with the $z$ integral evaluated before the $y$ integral. Infinite sums in Lemma 13 are over $\mathbb Z$. The Gaussian integrals and series in the statements are convergent, so Lean's totalized integral and sum conventions do not supply artificial zero values.

The paper's condition-number convention is recorded as $\kappa\ge1$ and Corollary 10's operating range as $0<\varepsilon<1/2$. Those conditions keep $D_\kappa$ nonempty and $\log(\kappa/\varepsilon)$ positive. The unnumbered reciprocal-exponential bound from p. 13 has $x\ne0$ explicitly, because the printed fractions have no value at zero. The $\Theta$ bounds quantify their constants before the problem parameters. In particular, choices with unbounded sample counts or vanishing mesh widths cannot satisfy the goal solely by brute-force approximation.

A complete development needs reusable Gaussian integration and tail estimates, complex exponential identities, summability of Gaussian integer series, and Poisson summation in Mathlib's analytic setting. Contributions that establish these general facts, the stated milestones, or Lemma 11 with its original bounds all fit this scope. This mission uses the reviewed shared definition of $D_\kappa$. Childs, Kothari and Somma, pp. 9–14.

Selected references

• A. M. Childs, R. Kothari and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM Journal on Computing 46(6), 2017; arXiv:1511.02306v2.

Why convergence bounds for Hamiltonian Monte Carlo matter

Hamiltonian Monte Carlo (HMC) is one of the most widely used Markov chain Monte Carlo methods for sampling from a probability density $\mu(dx)\propto e^{-U(x)}\,dx$ on $\mathbb R^d$. It is the default sampler of probabilistic programming systems such as Stan, and it is used throughout Bayesian statistics, molecular dynamics and machine learning. In practice the chain is run for a finite number of steps, and a user needs to know how many steps make its law close to $\mu$. That requires a quantitative, dimension-explicit convergence bound.

Before the paper of Bou-Rabee, Eberle and Zimmer, such bounds were known essentially only when $U$ is strongly convex everywhere: then two copies of the chain driven by the same random velocity contract deterministically (Mangoubi and Smith, arXiv:1708.07114). Bou-Rabee, Eberle and Zimmer (arXiv:1805.00452, Ann. Appl. Probab. 30(3), 2020) remove global convexity. They require $U$ to be strongly convex only outside a ball, and they introduce a coupling of two HMC transitions under which a suitably concave distance contracts on average for every pair of starting points. This mission formalizes that contraction theorem for exact HMC, together with its consequence, exponential convergence in the $L^1$ Wasserstein distance.

Setting

Let $U:\mathbb R^d\to\mathbb R$ be four times continuously differentiable with $\int e^{-U}<\infty$. The Hamiltonian is $H(x,v)=U(x)+\tfrac12|v|^2$, and the exact Hamiltonian flow $(q_t,p_t)$ solves

$$\tfrac{d}{dt}q_t=p_t,\qquad \tfrac{d}{dt}p_t=-\nabla U(q_t),\qquad (q_0(x,v),p_0(x,v))=(x,v).$$

Fix a duration $T>0$. One step of exact HMC from $x$ draws a fresh velocity $\xi\sim N(0,I_d)$ and moves to $X'(x)=q_T(x,\xi)$; its transition kernel is $\pi(x,\cdot)=\mathrm{law}(q_T(x,\xi))$, and $\mu$ is invariant for $\pi$.

Assumption 2.1 asks that (A1) $U$ has a local minimum at $0$ with $U(0)=0$; (A2) the second, third and fourth derivatives of $U$ are bounded, with $L=\sup\|\nabla^2U\|$, $M=\sup\|\nabla^3U\|$ and $N=\sup\|\nabla^4U\|$; and (A3) there are $\mathcal R\ge0$ and $K>0$ such that

$$(x-y)\cdot(\nabla U(x)-\nabla U(y))\ge K|x-y|^2\quad\text{whenever }|x-y|\ge\mathcal R.$$

The coupling. For two starting points $x,y$ with $z=x-y$, $e=z/|z|$, both copies use the same $\xi$. If $|x-y|\ge2\mathcal R$ the second copy also uses $\xi$ (synchronous coupling). If $|x-y|<2\mathcal R$ it uses

$$\eta=\begin{cases}\xi+\gamma z&\text{if }\tilde{\mathcal U}\le\varphi_{0,1}(e\cdot\xi+\gamma|z|)/\varphi_{0,1}(e\cdot\xi),\\ \xi-2(e\cdot\xi)e&\text{otherwise,}\end{cases}$$

with $\tilde{\mathcal U}\sim\mathrm{Unif}(0,1)$ independent of $\xi$ and $\varphi_{0,1}$ the standard normal density. This is a maximal coupling of $N(0,I_d)$ with its shift by $\gamma z$, completed by a reflection. The coupled positions are $X'=q_T(x,\xi)$ and $Y'=q_T(y,\eta)$, with $R'(x,y)=|X'-Y'|$.

The distance and the parameters. With $\gamma=\min(T^{-1},\mathcal R^{-1}/4)$, $a=T^{-1}$ and $R_1=\tfrac52(\mathcal R+T)$, set

$$f(r)=\int_0^r e^{-a\min(s,R_1)}\,ds,\qquad \rho(x,y)=f(|x-y|),$$

and let $\mathcal W_\rho(\nu,\eta)=\inf_{\gamma\in C(\nu,\eta)}\int\rho\,d\gamma$ be the Kantorovich distance over couplings $C(\nu,\eta)$.

Formalization targets

Goal: Theorem 2.3

Under Assumption 2.1 and $LT^2\le\min(K/L,\tfrac14,\tfrac1{16\Lambda})$ with $\Lambda=16L\mathcal R^2$, for all $x,y\in\mathbb R^d$,

$$E\bigl[f(R'(x,y))\bigr]\le(1-c)\,f(|x-y|),\qquad c=\frac1{10}\min\Bigl(1,\frac12KT^2\Bigl(1+\frac{\mathcal R}{T}\Bigr)e^{-\mathcal R/(2T)}\Bigr)e^{-2\mathcal R/T}.$$

This is the paper's main result for exact HMC, stated with its own constants.

Milestones (in the order of the proof)

1. (112)–(116): the parameter inequalities $\gamma T\le1$, $LT\le\gamma/4$, $\gamma\mathcal R\le\tfrac14$, $aT\ge1$, $R_1\ge\tfrac52(1+\gamma T)\mathcal R$.
2. Lemma 3.4 ($h=0$): $|q_t(x,v)-q_t(y,v)|^2\le(1-Kt^2/2)|x-y|^2$ for $|x-y|\ge2\mathcal R$.
3. Theorem 2.1: the synchronous coupling contracts far-apart points, $R'^2\le(1-\tfrac12KT^2)r^2$.
4. (119)–(120): $\inf_{r>0}rf'(r)/f(r)\ge\max(1,aR_1)e^{-aR_1}$, and the resulting bound on $c_1$.
5. Lemma 3.2 ($h=0$): first-variation bounds for the flow.
6. Term I: on $\{\eta=\xi+\gamma z\}$, $R'\le(1-\tfrac34\gamma T)r$.
7. Lemma 3.7 (62): $P[\xi-\eta\ne-\gamma z]\le|\gamma z|/\sqrt{2\pi}$.
8. (123): the truncated loss on the reflection event is at most $(\gamma T/\sqrt{2\pi})\,rf'(r)$.
9. Term III: $E[(R'-R_1)^+;\,\xi-\eta\ne-\gamma z]\le\tfrac54\gamma rT$.
10. (127): $\tfrac15\gamma T\max(1,2a\mathcal R)e^{-2a\mathcal R}\ge\tfrac1{10}e^{-2\mathcal R/T}$.

Companions

The coupling property ($\eta\sim N(0,I_d)$, p. 9), Corollary 2.8 ($\mathcal W_\rho(\nu\pi^n,\eta\pi^n)\le e^{-cn}\mathcal W_\rho(\nu,\eta)$, the $\mathcal W^1$ bound with factor $M=e^{\frac52(1+\mathcal R/T)}$, and the step count (41)), and the comparison (133), $e^{-aR_1}|x-y|\le\rho(x,y)\le|x-y|$.

Significance

Theorem 2.3 gives a contraction rate for exact HMC that depends on $d$ only through the constants $L$, $K$ and $\mathcal R$ of the potential, and that holds without global convexity. Through the coupling characterization of $\mathcal W_\rho$ it yields Corollary 2.8, an explicit number of steps after which the law of the chain is within $\epsilon$ of $\mu$ in $\mathcal W^1$. The same coupling and concave-distance technique underlies the paper's results for numerical (velocity Verlet) HMC and later work on Langevin and HMC couplings.

The theorems are proved in the paper; none of them is machine-checked. Formalizing them requires the Hamiltonian ODE flow and its first-variation estimates, a Gaussian maximal/reflection coupling, and the transfer from one-step contraction to Kantorovich contraction of Markov kernels. Each of these is reusable well beyond HMC. A formal development also settles the constants. Three printed displays of the proof are not delivered by the printed argument, and Theorem 2.1 as printed is false, while the goal itself survives: see Formalization scope.

Difficulty

The obvious argument fails for nearby points. When $|x-y|<2\mathcal R$, $U$ may be non-convex between them, and the synchronous coupling can push the two trajectories apart. The coupling must therefore change the velocity of one copy, and doing so makes the contraction only hold on average, on the event $\eta=\xi+\gamma z$. On the complementary reflection event the distance can grow, and the growth must be paid for by the concavity of $f$, through the truncation at $R_1$ and the tail term. Balancing these terms, uniformly in $x$, $y$ and $d$, is the heart of the proof. For far-apart points the difficulty is different: the contraction of Lemma 3.4 needs the trajectories to stay outside the non-convex region over the whole time interval, which is where the condition $LT^2<K/L$ enters.

Formalization scope

• Space and objects. $\mathbb R^d$ is EuclideanSpace ℝ (Fin d); $\nabla U$ is Mathlib's gradient, and $\nabla^kU$ is iteratedFDeriv. $L,M,N$ are upper bounds on the derivative norms, which is equivalent for every statement here. The flow is pinned by a predicate on $(q,p)$ (initial values and the ODE for all $t\in\mathbb R$); it is never a free function. Assumption 2.1 includes the standing hypotheses of §2.1: $U\in C^4$, $\int e^{-U}<\infty$, and $T$ deterministic.
• Probability. The randomness of one step is $(\xi,\tilde{\mathcal U})$ with law $N(0,I_d)\otimes\mathrm{Unif}(0,1)$. Expectations of nonnegative quantities are Lebesgue integrals in $[0,\infty]$. The signed expectation in (123) uses a Bochner integral together with an explicit integrability conclusion, so no bound can hold because a non-integrable integral is read as $0$. Kantorovich distances take values in $[0,\infty]$ over the published set of couplings MongeKantorovichYao.transferencePlans. The $n$-step kernel is the iterated composition $\pi^{n+1}=\pi\circ\pi^n$.
• Conventions at $\mathcal R=0$. $\gamma$ is $T^{-1}$ and $LT^2\le1/(16\Lambda)$ is written $256L^2\mathcal R^2T^2\le1$, so that no division by zero decides a constant.
• Corrections. The strict condition (56) in Lemma 3.4 is written $L^2t^2<K$, which preserves its meaning at $L=0$. Theorem 2.1 is stated as its proof delivers it, $R'^2\le(1-\tfrac12KT^2)r^2$. The printed $R'\le(1-\tfrac12KT^2)r$ fails for $U(x)=\tfrac12|x|^2$, where $R'=|\cos T|\,r$. The displays (117), (122), the last inequality of (124), and (126) are not posed, because the printed argument does not deliver their constants. The goal's rate $c$ of (32) is nevertheless unaffected: under the step condition and $d\ge1$, $c=\tfrac1{20}KT^2(1+\mathcal R/T)e^{-5\mathcal R/(2T)}$, which lies below both corrected rates $c_1$ and $c_2$.
• Ruled out. The goal is not trivialized by a vacuous hypothesis: the hypotheses are met by $U(x)=\tfrac12|x|^2$, $L=K=1$, $\mathcal R=0$, $T=\tfrac12$. Nor is it trivialized by a junk value: $f$, $R'$ and the integral are genuine, and the flow is the unique solution of the ODE.
• Contributions welcome. Gronwall-type first-variation bounds for Lipschitz Hamiltonian flows, the Gaussian maximal/reflection coupling (Lemma 3.7 and the coupling property), and the general passage from a one-step average contraction to $\mathcal W_\rho$ contraction of a Markov kernel.

Selected references

• N. Bou-Rabee, A. Eberle, R. Zimmer, Coupling and convergence for Hamiltonian Monte Carlo, Ann. Appl. Probab. 30(3), 2020. arXiv:1805.00452, doi:10.1214/19-AAP1528
• O. Mangoubi, A. Smith, Rapid mixing of Hamiltonian Monte Carlo on strongly log-concave distributions, 2017. arXiv:1708.07114
• A. Eberle, Reflection couplings and contraction rates for diffusions, Probab. Theory Related Fields 166, 2016. arXiv:1305.1233
• A. Eberle, A. Guillin, R. Zimmer, Couplings and quantitative contraction rates for Langevin dynamics, Ann. Probab. 47(4), 2019. arXiv:1703.01617