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

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

Motivation

Many optimization models minimize a smooth transformation followed by a nonsmooth penalty. A penalty may encourage sparse coefficients, restrict a feasible region, or measure deviation from observed data. Even when the penalty is convex, its composition with a nonlinear map need not be convex. This makes a standard convex optimality certificate unavailable, while numerical methods still need a way to judge whether a small step means that an iterate is close to a stationary solution. Drusvyatskiy and Lewis study this question for the prox-linear method, which repeatedly solves a convex approximation of the composite objective (Drusvyatskiy and Lewis, 2016, §5).

The paper establishes several links among step lengths, near-stationarity, local error bounds, and rates of convergence. The result selected here is Theorem 5.10: a local error bound expressed through the prox-linear step is equivalent, with stated constants, to subregularity of the composite subdifferential. This comparison lets analysts use whichever of the two properties is easier to verify for a given model (Theorem 5.10, p. 20).

Setting

Work in Euclidean spaces $\mathbb R^n$ and $\mathbb R^m$. A convex-composite objective has the form $\varphi(x)=g(x)+h(c(x))$. Here $g:\mathbb R^n\to\overline{\mathbb R}$ is proper, closed, and convex; $h:\mathbb R^m\to\mathbb R$ is finite-valued and convex; and $c:\mathbb R^n\to\mathbb R^m$ is continuously differentiable. The value $+\infty$ of $g$ can exclude points from its domain. The paper defines the relevant composite subdifferential by

$$\partial\varphi(x):=\partial g(x)+\nabla c(x)^*\partial h(c(x)),$$

and calls $x$ stationary when $0\in\partial\varphi(x)$. These are the paper's definitions on p. 13; the subdifferential is not replaced by a different generalized derivative (§5, p. 13).

At a current point $x$, the method substitutes the affine approximation $c(x)+\nabla c(x)(y-x)$ for $c(y)$. It then minimizes the prox-linear subproblem

$$\varphi_t(x;y)=g(y)+h(c(x)+\nabla c(x)(y-x))+\frac{\|x-y\|^2}{2t},\qquad t>0.$$

Write its unique minimizer as $x^t$ and the prox-gradient as $\mathcal G_t(x)=t^{-1}(x-x^t)$. Under the stated assumptions, a zero prox-gradient is equivalent to stationarity. The quantitative estimates additionally assume that $h$ is $L$-Lipschitz and that the Jacobian $\nabla c$ is $\beta$-Lipschitz, with $L,\beta\ge0$ (§5, p. 14).

For a set-valued mapping $F$, subregularity at a graph point $(\bar x,\bar y)$ with constant $l>0$ means that some neighborhood $\mathcal X$ of $\bar x$ satisfies

$$\operatorname{dist}(x,F^{-1}(\bar y))\le l\operatorname{dist}(\bar y,F(x)) \qquad(x\in\mathcal X).$$

The error bound in this mission applies this definition to $F=\mathcal G_t$ and $\bar y=0$. The other instance is $F=\partial\varphi$ (Definition 5.7, p. 19).

Formalization targets

The goal is the two-way quantitative statement of Theorem 5.10. If $\bar x$ is stationary and $t>0$, then

$$\mathcal G_t\text{ subregular at }(\bar x,0)\text{ with }\hat l \quad\Longrightarrow\quad \partial\varphi\text{ subregular there with }2\hat l,$$

and

$$\partial\varphi\text{ subregular there with }l \quad\Longrightarrow\quad \mathcal G_t\text{ subregular there with }(3L\beta t+2)l+2t.$$

The milestone list follows the paper's quantitative statements: the quadratic comparison (5.2), the gradient inequalities (5.3)–(5.5) of Lemma 5.1, Ekeland's variational principle as Theorem 5.2, the zero-step stationarity claim, the nearby near-stationary point of Theorem 5.3, and the step-length comparison (5.7) from Theorem 5.6. A final companion milestone is the per-step decrease (5.6), $\varphi(x_{k+1})\le\varphi(x_k)-\frac{1}{2L\beta}\|\mathcal G_t(x_k)\|^2$, stated at $t=1/(L\beta)$. The Ekeland result is an existing proved platform theorem, included by reference. The other targets are the paper's claims in their original order of dependence (§5, pp. 14–20).

Significance

The theorem turns an algorithmic condition, measured by a computed step, into a local property of the objective's first-order geometry. The constants specify how an error bound changes when moved between the two descriptions. In a convergence analysis, either condition can be established and then used through the other formulation without losing track of its quantitative dependence on $L$, $\beta$, and $t$. The result does not require a convex composite objective, only the convexity of $g$ and $h$ and the stated smoothness of $c$ (Theorem 5.10, p. 20).

The mathematical theorem was proved in the cited paper. The remaining work here is a machine-checked development of its composite subdifferential, the prox-linear argmin, the local subregularity estimates, and the comparison theorem. A general version of Ekeland's variational principle already has a machine-checked proof on Prove2Me; the prox-linear claims in this proposal are open proof targets. The definitions and estimates can also support later formalizations of prox-linear convergence, including the rate results in §§5–6, without building the algorithm into the statement of Theorem 5.10.

Difficulty

The prox-linear step is computed from a convex local model, while the quantity whose stationarity is sought belongs to the original nonlinear composite objective. A small displacement from $x$ to $x^t$ alone does not directly give a small subgradient at $x$: the model and the original objective have different first-order behavior away from their common base point. Theorem 5.3 addresses a nearby point $\hat x$ and supplies both a distance estimate and a subgradient estimate. Passing from those estimates to a neighborhood-wide error bound is the central issue, especially because the neighborhoods in the definition of subregularity are local to $\bar x$ (Theorem 5.3, p. 16; Theorem 5.10, p. 20).

Formalization scope

The Lean carrier is EuclideanSpace ℝ (Fin n); dimensions zero are allowed. Objective values use EReal, and properness explicitly excludes $-\infty$ and requires at least one finite point. Closedness is lower semicontinuity. The real-valued $h$ is embedded into EReal only when forming objective values. The Jacobian is fderiv, and its adjoint represents $\nabla c(x)^*$. The prox-linear point is an argmin predicate, so its definition does not depend on an arbitrary choice of minimizer.

The prox-gradient is represented as a set-valued map through those argmin points. Subregularity uses a neighborhood filter and requires the base graph point and a positive constant. Distance to an empty image on the right side is treated as $+\infty$ through a universal bound over its members. The inverse image on the left is nonempty because it contains the base point. In Theorem 5.3, the upper bound on distance to the subdifferential is represented by an attained subgradient; the subdifferential in this convex-composite setting is closed. These choices preserve the paper's claims at empty sets and extended-real values.

The goal retains every positive step size $t$. The companion (5.6) is stated at $t=1/(L\beta)$ with $L,\beta>0$: the page announces it for $t\le 1/(L\beta)$, but for $c(x)=x^2/2$, $h(r)=r$, $g=0$ the decrease is only $\frac{2t-t^2}{2}x_k^2<\frac12x_k^2$ when $t<1$. The linear-rate Theorem 5.5 and the results of §6 are not part of this mission. The first-order and step-length milestones do not require the $L$ and $\beta$ assumptions, matching the paper's first sentences. Contributions that establish the argmin properties, the distance conventions, the quantitative milestones, or the final two implications are within scope. A vacuous restriction to points already stationary would miss the neighborhood-wide error bound and is excluded.

Selected references

• D. Drusvyatskiy and A. S. Lewis, Error bounds, quadratic growth, and linear convergence of proximal methods, Mathematics of Operations Research 43(3), 2018. arXiv:1602.06661v2
• I. Ekeland, On the variational principle, Journal of Mathematical Analysis and Applications 47, 1974, 324–353. doi:10.1016/0022-247X(74)90025-0

Motivation

Obstacle problems, contact problems and control-constrained optimal control problems lead, after discretization, to complementarity systems of the form "find $y \le \psi$ and $\lambda \ge 0$ with $Ay + \lambda = f$ and $\lambda_i(y_i - \psi_i) = 0$". The primal-dual active set strategy solves such systems by guessing which constraints are active, solving one linear system, and updating the guess. Each iteration costs a single linear solve, and in practice the method stops after few iterations.

Hintermüller, Ito and Kunisch (SIAM J. Optim. 13(3), 2002; authors' version hal-01660511v1) showed that the strategy is a semismooth Newton method for a reformulation of the complementarity system through a max-function, which gives local superlinear convergence. For convergence from arbitrary initial data the paper gives three sufficient conditions on the matrix. The first is that $A$ is an M-matrix (Theorem 3.2), a class that many discretized second-order elliptic operators belong to. The M-matrix property is not stable under small perturbations: off-diagonal entries may become positive (p. 10). Theorem 3.4, the subject of this mission, shows that global convergence nevertheless survives small perturbations of an M-matrix. The paper notes that this also makes it relevant for numerical implementations.

Setting

Fix $n \ge 0$, vectors $f, \psi \in \mathbb{R}^n$ and a parameter $c > 0$. Vectors are ordered componentwise. For a matrix $A \in \mathbb{R}^{n\times n}$ the problem is

$$Ay + \lambda = f, \qquad \lambda - \max\big(0, \lambda + c(y - \psi)\big) = 0, \tag{3.1}$$

with the maximum taken componentwise. The paper writes $A$ for the matrix and calligraphic $\mathcal{A}$, $\mathcal{I}$ for the active and inactive index sets.

The primal-dual active set algorithm starts from an arbitrary pair $(y^0, \lambda^0)$. Given $(y^k, \lambda^k)$ it sets

$$\mathcal{A}_k = \{ i : \lambda^k_i + c(y^k_i - \psi_i) > 0 \}, \qquad \mathcal{I}_k = \{ i : \lambda^k_i + c(y^k_i - \psi_i) \le 0 \},$$

and determines $(y^{k+1}, \lambda^{k+1})$ from

$$Ay^{k+1} + \lambda^{k+1} = f, \qquad y^{k+1} = \psi \text{ on } \mathcal{A}_k, \qquad \lambda^{k+1} = 0 \text{ on } \mathcal{I}_k .$$

For an index set $\mathcal{I} \subseteq \{1,\dots,n\}$ with complement $\mathcal{A}$, $A_{\mathcal{I}}$ is the principal submatrix of $A$ on $\mathcal{I}$ and $A_{\mathcal{I}\mathcal{A}}$ the block with rows in $\mathcal{I}$ and columns in $\mathcal{A}$. The step above is uniquely solvable exactly when $A_{\mathcal{I}_k}$ is invertible.

A square matrix $M$ is an M-matrix if it is nonsingular, $m_{ij} \le 0$ for $i \ne j$, and $M^{-1} \ge 0$ entrywise. For a rectangular matrix $B$, $\|B\|_1$ is the matrix norm subordinate to the one-norms on both spaces, i.e. the maximal absolute column sum $\|B\|_1 = \max_j \sum_i |b_{ij}|$.

Formalization targets

Goal: Theorem 3.4

Let $M$ be an M-matrix. There is $\varepsilon > 0$ such that for every $K$ with $\|K\|_1 < \varepsilon$, every $f, \psi$ and every $c > 0$, with $A = M + K$:

1. (3.1) has exactly one solution $x^* = (y^*, \lambda^*)$;
2. the algorithm is well-defined: from every $(y^k, \lambda^k)$ the step has exactly one solution;
3. for arbitrary initial data,

$$\lim_{k\to\infty} (y^k, \lambda^k) = (y^*, \lambda^*).$$

The threshold $\varepsilon$ is left unspecified; the theorem asserts only its existence.

Milestones (the steps of the paper's proof)

• Principal submatrices $M_{\mathcal{I}}$ of an M-matrix are nonsingular M-matrices.
• With $\rho \ge \|M_{\mathcal{I}}^{-1}K_{\mathcal{I}}\|_1$ for every $\mathcal{I}$ and $\rho < \tfrac12$, every $A_{\mathcal{I}}$ is invertible and

$$A_{\mathcal{I}}^{-1} = \Big(I_{\mathcal{I}} + \sum_{i=1}^{\infty}\big(-M_{\mathcal{I}}^{-1}K_{\mathcal{I}}\big)^i\Big) M_{\mathcal{I}}^{-1}.$$

• Under the same hypothesis, for $v \le 0$, $\sum_{i\in\mathcal{I}} \big(-A_{\mathcal{I}}^{-1}v\big)_i \ge \frac{1-2\rho}{1-\rho}\,\|M_{\mathcal{I}}^{-1}v\|_1$.
• Whenever $A_{\mathcal{I}}$ is invertible, entrywise $A_{\mathcal{I}}^{-1}A_{\mathcal{I}\mathcal{A}} \le M_{\mathcal{I}}^{-1}K_{\mathcal{I}\mathcal{A}} - M_{\mathcal{I}}^{-1}K_{\mathcal{I}}(M+K)_{\mathcal{I}}^{-1}A_{\mathcal{I}\mathcal{A}}$.
• For $\|K\|_1$ small (uniformly in $f, \psi, c$) and $k \ge 1$: $\sum_i y^{k+1}_i < \sum_i y^k_i$ unless $y^{k+1} = y^k$, so $\mathcal{M}(y) = \sum_i y_i$ is a merit function.

A companion statement records the paper's remark after the theorem: for $\|K\|_1$ small, $M + K$ is a P-matrix.

Significance

The result. Theorem 3.2 makes global convergence depend on a sign pattern, which a perturbation of any size can destroy. Theorem 3.4 replaces the sign pattern by a neighbourhood: every matrix close enough to an M-matrix in $\|\cdot\|_1$ inherits unique solvability of (3.1), well-posedness of every iteration, and convergence from every starting point. As the paper observes, solvability of (3.1) for every $f$ is equivalent to $A$ being a P-matrix (Berman–Plemmons, Thm. 10.2.15), so the theorem also yields a sufficient condition for the P-matrix property.

Formalizing it. The paper's argument is a short sketch that refers back to the proof of Theorem 3.3 and asserts that "$\|K\|_1$ can be chosen sufficiently small". A formalization makes the uniformity of that choice explicit (it depends only on $M$), checks the strict decrease of the merit function including its equality case, and fixes the role of the first iterate, from which the sign facts used in (3.8) hold. To our knowledge none of the three global convergence theorems of the paper has a machine-checked proof.

Difficulty

The obvious argument is to transfer Theorem 3.2 by continuity. It fails because the M-matrix proof is monotonicity-based: it shows $y^* \le y^{k+1} \le y^k$, and that ordering uses $A_{\mathcal{I}}^{-1} \ge 0$ and $A_{\mathcal{I}\mathcal{A}} \le 0$, sign conditions that a perturbation breaks immediately. What survives is weaker: a scalar merit function $\sum_i y_i$ whose decrease must be quantified. The perturbation has to be controlled uniformly over all $2^n$ index sets, both in the multiplier term (where the decrease must remain strict whenever $\lambda^k_{\mathcal{I}} \ne 0$) and in the coupling term $A_{\mathcal{I}}^{-1}A_{\mathcal{I}\mathcal{A}}(y^k-\psi)_{\mathcal{A}}$, whose positive part must be absorbed by $-\sum_{\mathcal{A}}(y^k_i - \psi_i)$. Finally, strict decrease of a merit function does not by itself give convergence; the finiteness of the possible active sets must be combined with it.

Formalization scope

Vectors are Fin n → ℝ and matrices Matrix (Fin n) (Fin n) ℝ; index sets are Finset (Fin n) and blocks are submatrices along subtype inclusions. Conventions committed to:

• $\|B\|_1$ is the maximal absolute column sum, computed as a supremum in $\mathbb{R}_{\ge0}$, equal to $0$ for a matrix with no columns; it is not the entrywise sum or another operator norm;
• the M-matrix predicate requires nonsingularity explicitly (IsUnit of the determinant), because Lean's inverse of a singular matrix is $0$; it is stated for any finite index type so that it applies to principal blocks;
• "$\|K\|_1$ sufficiently small" is one $\varepsilon > 0$ depending only on $M$ (and $n$), quantified before $f$, $\psi$, $c$ and the initial data;
• a step is a relation and a run is an infinite sequence of steps from arbitrary initial data; step (iv)'s "Stop" is not modelled;
• the convergence conclusion is stated for every solution of (3.1), which by the first conclusion is the solution; convergence is in the product topology;
• the merit-function milestone holds from $k = 1$ on, as in the paper.

The goal does not assume that $M + K$ is a P-matrix, that its principal blocks are invertible, that $\rho < \tfrac12$, or that the merit function decreases: all of these are consequences. A version with $\varepsilon$ chosen after $f$, $\psi$, $c$ or the run, or with any of these as hypotheses, would be trivial or different and is ruled out. Well-definedness is a conclusion, so runs from every initial pair exist and the convergence conclusion is not vacuous.

In the g-estimate, the page's final $\|M^{-1}\lambda^k_{\mathcal{I}}\|_1$ is a misprint for $\|M_{\mathcal{I}}^{-1}\lambda^k_{\mathcal{I}}\|_1$; the statement uses $M_{\mathcal{I}}^{-1}$.

A complete development needs: principal submatrices of M-matrices, Neumann series for matrices with the column-sum norm, the block form of one step and the identity (3.8), the sign facts $\lambda^k \le 0$ on $\mathcal{I}_k$ and $y^k \ge \psi$ on $\mathcal{A}_k$ for $k \ge 1$, and continuity of principal minors for the P-matrix statement. The M-matrix and Neumann-series lemmas are reusable beyond this mission. Proofs of any milestone, of the companion, or of the goal through another argument are welcome.

Selected references

• M. Hintermüller, K. Ito, K. Kunisch, The primal-dual active set strategy as a semismooth Newton method, SIAM J. Optim. 13(3):865–888, 2002. https://doi.org/10.1137/S1052623401383558 ; authors' version https://hal.science/hal-01660511v1
• A. Berman, R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic Press, 1979 (SIAM Classics reprint 1994). https://doi.org/10.1137/1.9781611971262
• R. H. W. Hoppe, Multigrid algorithms for variational inequalities, SIAM J. Numer. Anal. 24:1046–1065, 1987. https://doi.org/10.1137/0724069

Motivation

Obstacle problems, contact problems and box-constrained quadratic programs share one finite-dimensional core: the linear complementarity problem "find $y \le \psi$ and $\lambda \ge 0$ with $Ay + \lambda = f$ and $\lambda_i(y_i - \psi_i) = 0$". After discretization by finite differences or piecewise linear finite elements, the matrix $A$ is typically an M-matrix, which is why this class sits at the centre of the numerical analysis of variational inequalities.

The primal-dual active set strategy solves such problems by guessing which constraints are active, solving one linear system on that guess, and updating the guess from the new primal and dual values. Hintermüller, Ito and Kunisch (SIAM J. Optim. 13 (2002) 865–888; authors' version hal-01660511) showed that this strategy is a semismooth Newton method for a nonsmooth reformulation of the problem. Local superlinear convergence then follows from Newton theory. This mission concerns the paper's global result: for M-matrices the method converges from every starting point, and its primal iterates decrease monotonically. A related global result for a purely primal active set method appears in Hoppe's 1987 analysis of multigrid methods for obstacle problems, which the paper cites as [H].

Setting

Fix $n \ge 0$, a matrix $A \in \mathbb{R}^{n\times n}$, vectors $f, \psi \in \mathbb{R}^n$ and a constant $c > 0$. Vectors are compared componentwise. The paper writes $A$ for the matrix and calligraphic $\mathcal{A}$, $\mathcal{I}$ for index sets; Lean uses A for the matrix and activeSet for $\mathcal{A}$.

Problem (3.1). Find $(y, \lambda) \in \mathbb{R}^n \times \mathbb{R}^n$ with

$$Ay + \lambda = f, \qquad \lambda - \max(0, \lambda + c(y - \psi)) = 0,$$

the maximum taken componentwise. For $c > 0$ the second equation says exactly $y \le \psi$, $\lambda \ge 0$ and $\lambda_i (y_i - \psi_i) = 0$ for every $i$.

M-matrix. $A$ is an M-matrix if it is nonsingular, its off-diagonal entries satisfy $a_{ij} \le 0$ for $i \ne j$, and its inverse is entrywise nonnegative, $A^{-1} \ge 0$.

Blocks. For an index set $\mathcal{I}$ with complement $\mathcal{A}$, $A_{\mathcal{I}}$ is the principal submatrix with rows and columns in $\mathcal{I}$, and $A_{\mathcal{I}\mathcal{A}}$ the submatrix with rows in $\mathcal{I}$ and columns in $\mathcal{A}$.

The algorithm. Start from arbitrary $(y^0, \lambda^0)$. Given $(y^k, \lambda^k)$, let

$$\mathcal{A}_k = \{ i : \lambda^k_i + c(y^k - \psi)_i > 0\}, \qquad \mathcal{I}_k = \{ i : \lambda^k_i + c(y^k - \psi)_i \le 0\},$$

and let $(y^{k+1}, \lambda^{k+1})$ solve

$$Ay^{k+1} + \lambda^{k+1} = f, \qquad y^{k+1} = \psi \text{ on } \mathcal{A}_k, \qquad \lambda^{k+1} = 0 \text{ on } \mathcal{I}_k.$$

A sequence obeying this rule for every $k$ is a run. The solution of (3.1) is written $x^* = (y^*, \lambda^*)$ and the iterates $x^k = (y^k, \lambda^k)$.

Formalization targets

Goal: Theorem 3.2 (p. 7)

If $A$ is an M-matrix, then for arbitrary initial data

$$x^k \to x^*, \qquad y^* \le y^{k+1} \le y^k \ \text{ for all } k \ge 1, \qquad y^k \le \psi \ \text{ for all } k \ge 2.$$

Milestones (Appendix A, p. 20)

1. For every partition $\mathcal{I}, \mathcal{A}$: $A_{\mathcal{I}}$ is regular, $A_{\mathcal{I}}^{-1} \ge 0$ and $A_{\mathcal{I}}^{-1} A_{\mathcal{I}\mathcal{A}} \le 0$.
2. Complementarity (A.1): $\lambda^k_i = 0$ or $y^k_i = \psi_i$, for all $i$ and $k \ge 1$, for every run and every matrix.
3. Monotonicity: $y^{k+1} \le y^k$ for $k \ge 1$.
4. Primal feasibility: $y^k \le \psi$ for $k \ge 2$.
5. Lower bound: $y^* \le y^k$ for $k \ge 1$.
6. Remark A.1 and dual feasibility: if $\lambda^{\bar k}_i < 0$ with $\bar k \ge 1$, then $\lambda^k_i = 0$ and $y^k_i \le \psi_i$ for all $k \ge \bar k + 1$; and $\lambda^k \ge 0$ for all $k$ beyond some $k_o$.

A companion item states that a run exists from every initial pair when $A$ is an M-matrix.

Significance

The result. Theorem 3.2 makes the active set strategy a globally convergent method for the discretized obstacle problem, with no line search, no merit function and no damping. The monotonicity statement says more: from the second iterate on, every primal iterate is feasible and lies above the solution, and the iterates descend to it. Combined with the paper's observation (p. 7) that convergent iterates of the method converge in finitely many steps, it gives finite termination from every starting point. The theorem is the model for later global results on the method under weaker hypotheses, two of which (a norm condition and small perturbations of M-matrices) appear in the same paper as Theorems 3.3 and 3.4.

Formalizing it. The theorem is proved in the paper; no machine-checked version is known. A formal proof requires the classical sign facts on principal blocks of M-matrices, which Mathlib does not contain, and a careful treatment of the iteration as a relation on pairs of vectors. The M-matrix block lemma is reusable well beyond this mission, for instance in comparison principles for discretized elliptic problems.

Difficulty

The obvious argument would bound a merit function or a contraction factor; neither is available, because the iteration map is discontinuous and the active sets change from step to step. The paper's argument is instead an order argument. Its central step is the monotonicity $y^{k+1} \le y^k$: it requires the complementarity (A.1) of the previous iterate, which holds only from $k = 1$ on, and the sign structure of the blocks $A_{\mathcal{I}_k}^{-1}$ and $A_{\mathcal{I}_k}^{-1} A_{\mathcal{I}_k \mathcal{A}_k}$ for every possible active set. The index bounds are sharp in the sense that the statements fail at $k = 0$ (respectively $k = 1$ for feasibility), because the initial data are arbitrary. Convergence of the multipliers is a separate step: negative multiplier coordinates must be shown to vanish permanently after one more iteration.

Formalization scope

• Vectors are Fin n → ℝ, matrices Matrix (Fin n) (Fin n) ℝ, the product is A *ᵥ y, and $\le$ on vectors is componentwise. Convergence of $x^k$ is in the product topology on $\mathbb{R}^n \times \mathbb{R}^n$.
• Nonsingularity of an M-matrix is IsUnit A.det. Lean's inverse of a singular matrix is the zero matrix, so without this clause every singular matrix with nonpositive off-diagonal entries would count as an M-matrix; the M-matrix definition is the page's, not an equivalent characterization.
• The solution is a hypothesis IsSolution A f ψ c ystar lamstar. The paper speaks of "the unique solution"; for an M-matrix it exists and is unique by Berman–Plemmons, so the hypothesis assumes nothing the paper does not have. The uniqueness fact is not restated.
• The algorithm is a relation. One step is IsStep, a run is IsRun with free initial data y 0, lam 0. The stopping option of step (iv) is not modelled. Because a universally quantified statement over runs would be empty if no run existed, the companion item run_exists states existence of a run from every initial pair.
• The constant $c > 0$ is a binder of every statement that needs it; (A.1) needs neither $c > 0$ nor any hypothesis on $A$, and has none.
• Remark A.1 is printed without "$\bar k \ge 1$"; its proof assumes it and the statement is false at $\bar k = 0$, so the bound is kept.
• A trivializing formalization is ruled out: no statement assumes monotonicity, feasibility or convergence of the iterates, and the M-matrix hypothesis cannot be met by a singular matrix.
• Needed infrastructure: principal submatrices of M-matrices (Berman–Plemmons, Ch. 6), and inverse-nonnegativity arguments for block systems. Proofs of the block lemma, of (A.1) and of the order lemmas are welcome independently of the goal.

The printed closing step of the proof writes $\lambda^k = Ay^k + f$; step (iii) gives $\lambda^k = f - Ay^k$. No item quotes or depends on that line.

Selected references

• M. Hintermüller, K. Ito, K. Kunisch, The primal-dual active set strategy as a semismooth Newton method, SIAM J. Optim. 13(3):865–888, 2002. https://doi.org/10.1137/S1052623401383558 (authors' version: https://hal.science/hal-01660511)
• A. Berman, R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic Press, 1979; reprinted as SIAM Classics in Applied Mathematics 9, 1994. https://doi.org/10.1137/1.9781611971262
• R. H. W. Hoppe, Multigrid algorithms for variational inequalities, SIAM J. Numer. Anal. 24(5):1046–1065, 1987. https://doi.org/10.1137/0724069
• L. Qi, J. Sun, A nonsmooth version of Newton's method, Math. Program. 58:353–367, 1993. https://doi.org/10.1007/BF01581275