Numerical Analysis

Motivation

Monotone variational inequalities cover convex minimisation, convex–concave saddle-point problems, Nash equilibria of monotone games and complementarity problems, and they are the standard model for these in optimization and operations research. First-order methods for them (extragradient, forward–backward–forward, reflected and projected gradient methods) need a stepsize below $1/L$, where $L$ is a global Lipschitz constant of the operator. That constant is often unknown, too pessimistic, or nonexistent: in composite minimisation with a locally smooth term, or in saddle-point problems with bilinear-plus-nonlinear couplings, the operator is only locally Lipschitz. The usual remedy is a linesearch, which costs extra operator or prox evaluations per iteration and complicates the complexity accounting.

Y. Malitsky, Golden Ratio Algorithms for Variational Inequalities (preprint 2018, Optimization Online 6598; published in Mathematical Programming, 2020, doi:10.1007/s10107-019-01416-w) proposes the Explicit Golden Ratio Algorithm (EGRAAL): its stepsizes are computed in closed form from the last two iterates, it uses one evaluation of $F$ and one proximal step per iteration, and it needs neither a Lipschitz constant nor a linesearch. This mission formalizes its main convergence theorem, Theorem 2 of the preprint.

Setting

Let $\mathcal E$ be a finite-dimensional real inner product space with norm $\|\cdot\|=\sqrt{\langle\cdot,\cdot\rangle}$. Let $g:\mathcal E\to(-\infty,+\infty]$ with domain $\operatorname{dom} g=\{x: g(x)<+\infty\}$, and $F:\operatorname{dom} g\to\mathcal E$. The variational inequality (1) asks for

$$z^*\in\mathcal E\quad\text{with}\quad \langle F(z^*),z-z^*\rangle+g(z)-g(z^*)\ge0\quad\forall z\in\mathcal E. \tag{1}$$

Its solution set is $S$. The standing assumptions are: (C1) $S\ne\emptyset$; (C2) $g$ is proper, convex and lower semicontinuous; (C3) $F$ is monotone on $\operatorname{dom} g$, $\langle F(u)-F(v),u-v\rangle\ge0$ for $u,v\in\operatorname{dom}g$.

The proximal operator is $\operatorname{prox}_g(w)=\operatorname{argmin}_x\{g(x)+\tfrac12\|x-w\|^2\}$. Write $\varphi=\frac{\sqrt5+1}{2}$ for the golden ratio. Algorithm 1 (EGRAAL) takes $z^0,z^1\in\mathcal E$, $\lambda_0>0$, a parameter $\phi\in(1,\varphi]$ and a cap $\bar\lambda>0$, sets $\bar z^0=z^1$, $\theta_0=1$, $\rho=\frac1\phi+\frac1{\phi^2}$, and for $k\ge1$ computes

$$\lambda_k=\min\Big\{\rho\lambda_{k-1},\ \frac{\phi\theta_{k-1}}{4\lambda_{k-1}}\frac{\|z^k-z^{k-1}\|^2}{\|F(z^k)-F(z^{k-1})\|^2},\ \bar\lambda\Big\},\qquad \bar z^k=\frac{(\phi-1)z^k+\bar z^{k-1}}{\phi},$$ $$z^{k+1}=\operatorname{prox}_{\lambda_k g}\big(\bar z^k-\lambda_kF(z^k)\big),\qquad \theta_k=\frac{\lambda_k}{\lambda_{k-1}}\phi,$$

with the convention $0/0=+\infty$ in the middle term. The paper uses the bifunction $\Psi(u,v)=\langle F(u),v-u\rangle+g(v)-g(u)$.

Formalization targets

Goal: Theorem 2

If $F$ is locally Lipschitz continuous and (C1)–(C3) hold, then for every run of Algorithm 1 there is $z^*\in S$ with

$$z^k\to z^*\qquad\text{and}\qquad \bar z^k\to z^*.$$

Nothing is fixed beyond the paper's parameter ranges: $\phi\in(1,\varphi]$, $\bar\lambda>0$, $\lambda_0>0$ and the starting points are arbitrary. The two sequences share one limit.

Milestones

1. Eq. (4), the prox-inequality: $\bar x=\operatorname{prox}_g w\iff\langle\bar x-w,x-\bar x\rangle\ge g(\bar x)-g(x)$ for all $x$.
2. Eq. (18), the estimates the step rule gives: $\lambda_k\le\rho\lambda_{k-1}$, $\theta_k\le1+\frac1\phi$, and $\lambda_k^2\|F(z^k)-F(z^{k-1})\|^2\le\frac{\theta_k\theta_{k-1}}4\|z^k-z^{k-1}\|^2$.
3. Eq. (24), an identity that follows from the averaging step: $\|z^{k+1}-z\|^2=\frac\phi{\phi-1}\|\bar z^{k+1}-z\|^2-\frac1{\phi-1}\|\bar z^k-z\|^2+\frac1\phi\|z^{k+1}-\bar z^k\|^2$.
4. Eq. (27), the energy inequality, for $z\in\operatorname{dom}g$ and $k\ge2$.
5. Lemma 2: along bounded runs, $(\lambda_k)$ and $(\theta_k)$ are bounded and bounded away from $0$.
6. Lemma 1 (Bauschke–Combettes, Theorem 5.5): a Fejér monotone sequence whose cluster points lie in a nonempty set $C$ converges to a point of $C$.

Significance

The result. Theorem 2 shows that a monotone variational inequality with a locally Lipschitz operator can be solved by a method whose stepsizes adapt to the local curvature of $F$ at no extra cost: one $F$ evaluation and one prox step per iteration, and no global constant and no backtracking. Because $F$ is only ever evaluated at the prox outputs $z^k\in\operatorname{dom}g$, the method also applies when $F$ is undefined or badly behaved outside the feasible set, where reflected-gradient methods can fail. The same analysis gives an ergodic $O(1/k)$ rate and, under an error bound, an $R$-linear rate (§2.2 of the preprint; not part of this mission). The paper also derives fixed-point algorithms for demi-contractive operators from it.

Formalizing it. The theorem has a published proof; to our knowledge no machine-checked version exists, and Mathlib has no proximal operator, no theory of monotone variational inequalities and no Fejér-monotonicity lemma. This mission produces a formal convergence proof for an adaptive first-order method with a nonsmooth convex term, together with reusable pieces: the prox-inequality for extended-real-valued convex functions, a finite-dimensional Fejér convergence lemma, and a formal model of an adaptive-step algorithm with the $0/0=+\infty$ rule.

Difficulty

The usual convergence argument for projected or extragradient methods bounds the cross term $\langle F(z^k)-F(z^{k-1}),z^k-z^{k+1}\rangle$ using a global Lipschitz constant and a fixed stepsize. Here neither exists. The stepsize at iteration $k$ depends on the iterates, and the energy that decreases changes from step to step, since it involves $\theta_{k-1}$. Local Lipschitz continuity gives a usable constant only once the iterates are known to be bounded, and boundedness has to come from the energy inequality. Stepsizes that tend to $0$ would also break the argument (Lemma 2 excludes this for bounded runs). The last step, identifying cluster points as solutions, needs lower semicontinuity of $g$ and a limit in the prox-inequality along a subsequence with convergent stepsizes.

Formalization scope

$\mathcal E$ is a real InnerProductSpace with FiniteDimensional ℝ E. $g$ is a map E → EReal; (C2) is a structure: $g$ never takes the value $\bot$, is finite somewhere, has a convex epigraph in $E\times\mathbb R$, and is LowerSemicontinuous on $E$. $F$ is a total map E → E, and every hypothesis on it (monotonicity, Lipschitz bounds) is restricted to $\operatorname{dom}g$. The variational inequality is stored as $g(z^*)\le\langle F(z^*),z-z^*\rangle+g(z)$, with $z^*\in\operatorname{dom}g$, which avoids extended-real subtraction. $\operatorname{prox}_{\lambda g}$ is an argmin predicate, so it never produces a junk value. The algorithm is a predicate on the four sequences, written for index $k+1$. The rule $0/0=+\infty$ is a case split: if $F(z^k)=F(z^{k-1})$ the step is $\min\{\rho\lambda_{k-1},\bar\lambda\}$. No condition such as $F(z^1)\ne F(z^0)$ or $\lambda_0\le\bar\lambda$ is imposed.

"Locally Lipschitz" is formalized as Lipschitz on every bounded subset of $\operatorname{dom}g$. This is the property the proof of Lemma 2 uses. It agrees with local Lipschitz continuity when $\operatorname{dom}g$ is closed (for example $g=\delta_C$ for a closed convex $C$, or $g$ finite everywhere) and is stronger otherwise. Eq. (27) is stated for $z\in\operatorname{dom}g$, where $F(z)$ and $\Psi(z,z^k)$ are defined. Lemma 1 carries the hypothesis $C\ne\emptyset$ of its cited source, without which it is false.

The statements are not vacuous: $g\equiv0$, $F\equiv0$ satisfy (C1)–(C3) and the Lipschitz hypothesis, and admit a run of Algorithm 1 with $\lambda_k=\min\{\rho\lambda_{k-1},\bar\lambda\}$. A proof of Theorem 2 must hold for every run with the paper's parameters, not only for such degenerate data.

Contributions welcome: the prox-inequality and the existence of the prox for proper convex lsc $g$ (both reusable beyond this mission), the Fejér lemma, the algebraic estimates (18) and (24), and the energy inequality (27). Once these are in place, Lemma 2 and the cluster-point argument complete Theorem 2.

Selected references

• Y. Malitsky, Golden Ratio Algorithms for Variational Inequalities, preprint, Optimization Online 6598, 2018. https://optimization-online.org/wp-content/uploads/2018/05/6598.pdf ; published in Mathematical Programming 184 (2020), 383–410. https://doi.org/10.1007/s10107-019-01416-w
• H. H. Bauschke, P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer, 2011 (Theorem 5.5). https://doi.org/10.1007/978-1-4419-9467-7
• G. M. Korpelevich, The extragradient method for finding saddle points and other problems, Ekonomika i Matematicheskie Metody 12 (1976), 747–756.
• Y. Malitsky, Projected reflected gradient methods for monotone variational inequalities, SIAM Journal on Optimization 25 (2015), 502–520. https://doi.org/10.1137/14097238X

Motivation

A variational inequality asks for a point $x$ of a set $K\subseteq\mathbb R^n$ at which a vector field $f$ points "into" $K$: $(x'-x)^T f(x)\ge 0$ for every $x'\in K$. It is the common form of the first-order optimality conditions of constrained optimization, of complementarity problems, and of equilibrium models in economics and traffic networks. In many of these models the feasible set itself depends on the decision: the admissible actions of one agent are restricted by the current state, as in the impulse-control problems of Bensoussan and Lions that motivated quasi-variational inequalities, where $K=K(x)$.

D. Chan and J. S. Pang, The generalized quasi-variational inequality problem (Math. Oper. Res. 7 (1982) 211–222), unify the quasi-variational inequality with the generalized (set-valued) variational inequality of Fang and Peterson (JOTA 1982). Their §§3–4 prove existence by fixed-point theorems for set-valued maps; §5 takes a different route and characterizes solutions as fixed points of a composite projection map. Theorem 5.3, the subject of this mission, gives conditions under which that map is a contraction, so that its fixed point exists, is unique, solves the problem, and is computed by plain fixed-point iteration from any starting point. It is the algorithmic result of the paper, and an early instance of the projection methods for strongly monotone quasi-variational inequalities studied since (e.g. Nesterov and Scrimali 2011).

Setting

Throughout, $\mathbb R^n$ carries the Euclidean inner product $x^T y$ and norm $\|x\|$.

Given point-to-set mappings $K$ and $f$ of $\mathbb R^n$ into itself, the generalized quasi-variational inequality problem $\mathrm{GQVI}(K,f)$ is to find vectors $x$ and $y$ with

$$x\in K(x),\qquad y\in f(x),\qquad (x'-x)^T y\ge 0\quad\text{for all }x'\in K(x).$$

When $f$ is point-to-point, $f(x)$ is read as the singleton $\{f(x)\}$.

For a set $S$ and a point $z$, the projection $P_S(z)$ is the point of $S$ nearest to $z$, $P_S(z)=\operatorname{sol}\min_{x\in S}\|x-z\|$; it exists and is unique when $S$ is nonempty, closed and convex.

Theorem 5.3 concerns the special structure in which the feasible set moves by translation: fix a nonempty closed convex set $\tilde K$ and a point-to-point mapping $m$, and put

$$K(x)=m(x)+\tilde K=\{m(x)+k : k\in\tilde K\}.$$

For a step length $\lambda>0$ and a point-to-point $f$, the projection map is

$$F_\lambda(x)=P_{K(x)}\bigl(x-\lambda f(x)\bigr).$$

The mappings $m$ and $f$ are assumed Lipschitz continuous with constants $\alpha$, $\beta$ ($\|m(x)-m(y)\|\le\alpha\|x-y\|$, $\|f(x)-f(y)\|\le\beta\|x-y\|$) and strongly monotone with constants $\gamma$, $\delta$ ($(x-y)^T(m(x)-m(y))\ge\gamma\|x-y\|^2$, $(x-y)^T(f(x)-f(y))\ge\delta\|x-y\|^2$).

Formalization targets

Goal: Theorem 5.3 (p. 221)

For each $\lambda>0$ with

$$\lambda^2\beta^2+2\lambda(\alpha\beta-\delta)-2(\gamma-\alpha)<0,$$

the map $F_\lambda$ is a contraction (Lipschitz with a constant $c<1$ independent of the points), it has a fixed point $\tilde x_\lambda$, the point $\tilde x_\lambda$ solves $\mathrm{GQVI}(K,f)$, and the iterates $x^{k+1}=F_\lambda(x^k)$ converge to $\tilde x_\lambda$ from every initial vector $x^0\in\mathbb R^n$. All four conclusions are stated together.

Milestones

1. Projection onto a translate (§5, proof of Theorem 5.3, first display, p. 221): $P_{K(x)}(y)=m(x)+P_{\tilde K}(y-m(x))$ for all $x,y$.
2. Lipschitz estimate (§5, proof of Theorem 5.3, last display, p. 221): for every $\lambda>0$,

$$\|F_\lambda(y^1)-F_\lambda(y^2)\|\le\Bigl[\alpha+\bigl(\lambda^2\beta^2+2\lambda(\alpha\beta-\delta)+(1+\alpha^2-2\gamma)\bigr)^{1/2}\Bigr]\|y^1-y^2\|.$$

3. Theorem 5.1 (p. 220): if every $K(x)$ is closed and convex, $(x^*,y^*)$ solves $\mathrm{GQVI}(K,f)$ if and only if $x^*=P_{K(x^*)}(x^*-y^*)$ and $y^*\in f(x^*)$.

Significance

The result. Theorem 5.3 turns an existence question into a computation: under Lipschitz and strong monotonicity assumptions, a quasi-variational inequality with translated feasible sets has exactly one solution reachable by projection iterations, each of which is a projection on the fixed set $\tilde K$ (a convex quadratic program when $\tilde K$ is polyhedral). The step-size window it gives is explicit in $\alpha,\beta,\gamma,\delta$, so it certifies a convergent method before any iteration is run. The closing remark of the paper (p. 222) reads each step as solving the GQVI under a zero-th order approximation of $K$, the viewpoint behind later splitting methods.

Formalizing it. The result is proved in the paper; the proof is short, but its constants and the equivalence of the two contraction conditions are easy to get wrong. A machine-checked version fixes the exact hypotheses (no sign conditions on the constants, Euclidean geometry), and produces reusable pieces: the translation identity for projections, nonexpansiveness of the Euclidean projection on a closed convex set, and the projection characterization of quasi-variational inequalities (Theorem 5.1).

Difficulty

Banach's fixed-point theorem does the last step; the work is the estimate. The naive bound, projection nonexpansiveness applied directly to $F_\lambda$, fails because the sets $K(y^1)$ and $K(y^2)$ differ: two projections on different sets are not controlled by the distance of the projected points alone. The translation identity separates the moving part $m(y^1)-m(y^2)$ from a projection on the one set $\tilde K$, at the cost of the additive term $\alpha$ in the constant. The remaining square must be expanded with the inner-product cross terms bounded by the monotonicity constants in the right directions, including the cross term between $f$ and $m$. Finally, the condition "bracket $<1$" is equivalent to the stated $\lambda$-condition only when $\alpha<1$, which must be derived from the hypotheses rather than assumed.

Formalization scope

• The space is EuclideanSpace ℝ (Fin n), with Mathlib's Euclidean norm and inner product; Fin n → ℝ (sup norm) would change every constant. No assumption $n\ge1$ is made; at $n=0$ all statements hold trivially.
• The constants $\alpha,\beta,\gamma,\delta$ are real numbers with no sign conditions, as in the paper; the Lipschitz and monotonicity hypotheses are the displayed inequalities for all $x,y$. For $n\ge1$ they force $\alpha,\beta\ge0$, $\gamma\le\alpha$, $\delta\le\beta$, and the $\lambda$-condition then forces $\alpha<1$ and $\delta>\alpha\beta$.
• $K(x)$ is the translate $\{m(x)+k : k\in\tilde K\}$ of a fixed set $\tilde K$, assumed nonempty, closed and convex. The projection is a nearest-point function proj that returns a junk value only when no nearest point exists; under the hypotheses of every statement using it, the nearest point exists and is unique, so proj is the paper's $P$. Theorem 5.1 is stated relationally (nearest-point predicate IsProj) to avoid junk values altogether.
• "Contraction" is Mathlib's ContractingWith c F with c : ℝ≥0: c < 1 and a Lipschitz bound with that single constant. A constant allowed to depend on the points, or a Lipschitz bound without c < 1, is not a contraction and would trivialize the goal; so would a projection whose junk value is reachable (e.g. with $\tilde K$ empty), which makes $F_\lambda$ unrelated to the paper's map. The fixed point must be linked to the GQVI and to the iteration from every starting point.
• The square root in the Lipschitz estimate is Real.sqrt; its radicand is nonnegative under the hypotheses when $n\ge1$.
• Useful infrastructure: Mathlib's ContractingWith.fixedPoint and ContractingWith.tendsto_iterate_fixedPoint (Banach), exists_norm_eq_iInf_of_complete_convex and norm_eq_iInf_iff_real_inner_le_zero (projection on convex sets), and the platform theorem VectorSpaceOpt.min_distance_convex_set. A general lemma that the Euclidean nearest-point map of a closed convex set is 1-Lipschitz is reusable well beyond this mission and is welcome as a separate contribution.

Selected references

• D. Chan and J. S. Pang, The generalized quasi-variational inequality problem, Mathematics of Operations Research 7(2) (1982) 211–222. https://doi.org/10.1287/moor.7.2.211
• S. C. Fang and E. L. Peterson, Generalized variational inequalities, Journal of Optimization Theory and Applications 38 (1982) 363–383. https://doi.org/10.1007/BF00935344
• Y. Nesterov and L. Scrimali, Solving strongly monotone variational and quasi-variational inequalities, Discrete and Continuous Dynamical Systems 31(4) (2011) 1383–1396. https://doi.org/10.3934/dcds.2011.31.1383

Motivation

Quadratic programming, the minimisation of a quadratic function subject to linear inequality constraints, is a basic subproblem of nonlinear optimisation (sequential quadratic programming, trust-region methods) and a model in its own right in portfolio selection, least-squares estimation and control. The classical solution methods are active-set methods: they keep a set of constraints treated as equalities, minimise over the resulting affine set, and then decide which constraint to drop or add (Gill, Murray and Wright, Practical Optimization, 1981; Fletcher, Practical Methods of Optimization, Vol. 2, 1981). Their finite-termination proofs need either a nondegeneracy assumption (linearly independent active constraints) or an anti-cycling rule, because under degeneracy the choice of the constraint to drop, made from Lagrange multiplier estimates, can cycle.

Calamai and Moré (Mathematical Programming 39, 1987) showed that the gradient projection method can take over the step that leaves a working set. Their Algorithm 6.1 alternates two kinds of step: an arbitrary non-increasing step that adds constraints to the working set until the equality-constrained subproblem is solved, and a single projected-gradient step once it is solved. Theorem 6.2 states that this algorithm terminates at a stationary point for every quadratic that is bounded below on the feasible polyhedron, with no nondegeneracy assumption and no anti-cycling rule. The same paper's Sections 2–4 (the subject of the first two missions of this series) supply the properties of the gradient projection step that the argument uses.

Timeline of the ingredients:

• 1964, 1966: Goldstein and Levitin–Polyak introduce the gradient projection method $x_{k+1} = P(x_k - \alpha_k\nabla f(x_k))$ for convex constraint sets.
• 1976: Bertsekas proves finite identification of the active constraints for bound constraints and the Armijo rule.
• 1981: Dunn uses the descent inequalities (2.4)–(2.5) in the analysis of the method.
• 1987: Calamai and Moré generalise the step rule to (2.1)–(2.2), prove convergence of projected gradients, identification of active constraints for general polyhedra, and finite termination of Algorithm 6.1.

Setting

Let $E$ be a finite-dimensional real inner product space (the paper's $\mathbb{R}^n$ with a general inner product). The feasible set is a polyhedron

$$\Omega = \{x \in E : \langle c_j, x\rangle \ge \delta_j,\ j = 1, \dots, m\},$$

with active set $A(x) = \{j : \langle c_j, x\rangle = \delta_j\}$. The objective is a quadratic function $f(x) = \tfrac12\langle x, Qx\rangle + \langle b, x\rangle + c_0$ with $Q$ self-adjoint but not necessarily positive semidefinite, so $f$ may be nonconvex. Its gradient $\nabla f$ is taken with respect to the inner product of $E$.

The projection into $\Omega$ is $P(x) = \operatorname{argmin}\{\|z - x\| : z \in \Omega\}$, and a point $x^* \in \Omega$ is stationary if $\langle\nabla f(x^*), x - x^*\rangle \ge 0$ for every $x \in \Omega$.

A gradient projection step from $x_k$ is $x_{k+1} = P(x_k - \alpha_k\nabla f(x_k))$ with $\alpha_k > 0$ satisfying the sufficient decrease condition (2.1) with a constant $\mu_1 \in (0,1)$, the condition (2.2) that $\alpha_k \ge \gamma_1$ or $\alpha_k \ge \gamma_2\bar\alpha_k > 0$ for some $\bar\alpha_k$ at which the decrease test (2.3) with constant $\mu_2 \in (0,1)$ fails, and the upper bound $\alpha_k \le \gamma_3$ (3.2).

A working set is a set $W \subseteq \{1, \dots, m\}$; problem (6.2) is $\min\{f(y) : \langle c_j, y\rangle = \delta_j,\ j \in W\}$, over an affine set that ignores the inequality constraints outside $W$.

Algorithm 6.1 produces iterates $x_k \in \Omega$ and working sets $W_k \subseteq A(x_k)$ from $x_0 \in \Omega$:

• (a) if $x_k$ is a global minimiser of (6.2) for $W_k$, then $x_{k+1}$ is a gradient projection step from $x_k$;
• (b) otherwise $x_{k+1} \in \Omega$, $f(x_{k+1}) \le f(x_k)$, $W_k \subseteq W_{k+1}$, and if $W_{k+1} = W_k$ then $x_{k+1}$ is a global minimiser of (6.2).

Formalization targets

Goal: Theorem 6.2

For every quadratic $f$ bounded below on $\Omega$, all constants $\gamma_1, \gamma_2 > 0$, $\mu_1, \mu_2 \in (0,1)$, $\gamma_3 \in \mathbb{R}$, and every run $(x_k, W_k, \alpha_k)_{k\ge0}$ of Algorithm 6.1,

$$\exists\, l \ge 0 :\quad \langle \nabla f(x_l), x - x_l\rangle \ge 0 \quad \text{for all } x \in \Omega.$$

The theorem makes no assumption on the boundedness of the iterates and none on the linear independence of the constraints.

Milestones

• Lemma 2.1(a): for nonempty closed convex $\Omega$, $z \in \Omega$ and any $x$, $\langle P(x) - x, z - P(x)\rangle \ge 0$.
• Eq. (2.5): for $x_k \in \Omega$, $\alpha_k > 0$ and $x_{k+1} = P(x_k - \alpha_k\nabla f(x_k))$,

$$\langle\nabla f(x_k), x_k - x_{k+1}\rangle \ge \frac{\|x_{k+1} - x_k\|^2}{\alpha_k}.$$

Significance

Theorem 6.2 separates the two roles an active-set method plays: solving equality-constrained subproblems, for which any method that does not increase $f$ may be used, and choosing the next working set, which the gradient projection step does. The consequence is a finitely terminating quadratic programming algorithm for nonconvex quadratics that needs neither nondegeneracy nor an anti-cycling rule, and a template for large-scale bound-constrained and linearly constrained solvers that combine projection steps with subspace minimisation.

The result is proved in the paper and is classical; no machine-checked version is known. A formalization produces a checked finite-termination theorem for an active-set method on degenerate problems, together with reusable pieces: the projection onto a polyhedron as a total function with its variational inequality, the descent estimate of a projected step, and a predicate describing active-set runs with working sets, which other active-set algorithms can reuse.

Difficulty

The obvious argument — each step decreases $f$ and there are finitely many working sets — fails on two counts. First, step (b) only guarantees $f(x_{k+1}) \le f(x_k)$, so $f$ values alone do not rule out infinitely many iterations; the nesting of working sets and the final clause of step (b) are what bound consecutive (b)-steps. Second, a gradient projection step taken at a solution of (6.2) can in principle leave $f$ unchanged, and nothing in the algorithm's rules says directly that it makes progress; that it does so at every non-stationary iterate is a property of the projection and of the step conditions (2.1)–(2.2), not of the algorithm. A further point is that the minimum of (6.2) is taken over an affine set, not over $\Omega$, and the link between the value at a step-(a) iterate and later iterates runs through the requirement $W_k \subseteq A(x_k)$.

Formalization scope

The space is a finite-dimensional real inner product space E; $\nabla f$ is Mathlib's gradient. Constraints are indexed by Fin m, working and active sets are Finset (Fin m), and a run is a predicate IsAlgorithm61Run on three sequences $x : \mathbb{N} \to E$, $W$, $\alpha$ indexed from $0$. The run does not stop by itself; "the algorithm terminates at a stationary iterate" is rendered as the existence of an index $l$ with $x_l$ stationary. The projection is argmin made total by a junk value $0$ that is unreachable when $\Omega$ is nonempty, closed and convex. The paper's "$\subset$" between working sets is inclusion. A quadratic function is $\tfrac12\langle x, Qx\rangle + \langle b, x\rangle + c_0$ with $Q$ symmetric and no definiteness assumption.

A trivializing formalization — a run predicate that forces $x_0$ to be stationary, or that no sequence satisfies — would make the goal empty; the step rules here are the paper's verbatim, and a run on $\Omega = [0,\infty) \subset \mathbb{R}$ with $f(x) = x$ starting at the non-stationary point $x_0 = 1$ satisfies the predicate. Replacing step (a) by "any step that strictly decreases $f$" would assume the central fact and is not acceptable.

A complete development needs: the variational inequality of the projection (Lemma 2.1(a)), the descent estimate (2.5), the characterisation of stationary points as fixed points of the projected step, and a finiteness argument over the finitely many subsets of Fin m. Proofs of the milestones, of these auxiliary facts, and of the goal are all welcome.

Selected references

• P. H. Calamai and J. J. Moré, Projected gradient methods for linearly constrained problems, Mathematical Programming 39 (1987) 93–116. https://doi.org/10.1007/BF02592073
• A. A. Goldstein, Convex programming in Hilbert space, Bulletin of the AMS 70 (1964) 709–710. https://doi.org/10.1090/S0002-9904-1964-11178-2
• E. S. Levitin and B. T. Polyak, Constrained minimization methods, USSR Computational Mathematics and Mathematical Physics 6 (1966) 1–50. https://doi.org/10.1016/0041-5553(66)90114-5
• D. P. Bertsekas, On the Goldstein–Levitin–Polyak gradient projection method, IEEE Transactions on Automatic Control 21 (1976) 174–184. https://doi.org/10.1109/TAC.1976.1101194
• J. C. Dunn, Global and asymptotic convergence rate estimates for a class of projected gradient processes, SIAM Journal on Control and Optimization 19 (1981) 368–400. https://doi.org/10.1137/0319022
• P. E. Gill, W. Murray and M. H. Wright, Practical Optimization, Academic Press, 1981. https://doi.org/10.1137/1.9781611975604

Motivation

Many computations in numerical linear algebra, statistics and computational physics need the trace of a matrix $A$ that is never formed explicitly: $A$ may be an inverse, a matrix function $f(B)$, or a product of large operators, and the only affordable access is a routine that returns $Av$ or $v^TAv$ for a given vector $v$. Monte Carlo trace estimators handle this setting: draw random vectors $z$ and average the quadratic forms $z^TAz$, each of which costs one matrix–vector product.

Avron and Toledo (J. ACM 2011) compare such estimators by the number of samples $M$ that guarantee relative error $\epsilon$ with probability $1-\delta$, and by the number of random bits each sample consumes. Hutchinson's estimator (Hutchinson 1990) and the Gaussian estimator need $\Omega(n)$ random bits per sample. Section 8 of the paper studies two estimators that sample only from the $n$ standard basis vectors and so need about $\log_2 n$ bits per sample, which allows the samples to be generated in advance. The plain version has a sample bound that depends on how uneven the diagonal of $A$ is; the mixed version first multiplies $A$ on both sides by a random orthogonal mixing matrix of the kind introduced by Ailon and Chazelle (2006) for the fast Johnson–Lindenstrauss transform and used by Avron, Maymounkov and Toledo (2010) in least-squares solvers. This mission formalizes the resulting sample bound, Theorem 8.4.

Setting

Let $n \ge 1$, let $A \in \mathbb{R}^{n\times n}$ be symmetric positive semi-definite, and let $e_1,\ldots,e_n$ be the standard basis of $\mathbb{R}^n$.

A random variable $T$ is an $(\epsilon,\delta)$-approximator of $\mathrm{trace}(A)$ if

$$\Pr\bigl(|T-\mathrm{trace}(A)| \le \epsilon\,\mathrm{trace}(A)\bigr) \ge 1-\delta$$

(Definition 4.1).

The unit vector estimator with $M$ samples is

$$U_M = \frac{n}{M}\sum_{i=1}^M z_i^TAz_i,$$

where $z_1,\ldots,z_M$ are independent uniform random samples from $\{e_1,\ldots,e_n\}$ (Definition 3.4). Each term $z_i^TAz_i$ is a diagonal entry of $A$ chosen uniformly at random. Its behaviour is governed by

$$r_D(A) = \frac{n\cdot\max_i A_{ii}}{\mathrm{trace}(A)},$$

which lies between $1$ and $n$.

A random mixing matrix is $\mathcal F = FD$, where the seed $F$ is a fixed orthogonal $n\times n$ matrix and $D$ is diagonal with i.i.d. Rademacher entries, $\Pr(D_{ii}=\pm1) = 1/2$ (Definition 3.5). The seed enters through

$$\eta = \max_{i,j}|F_{ij}|^2,$$

which satisfies $1/n \le \eta \le 1$; normalized DFT and Hadamard matrices attain $\eta = 1/n$, DCT and DHT matrices have $\eta = 2/n$ (p. 8:5).

The mixed unit vector estimator is

$$T_M = \frac{n}{M}\sum_{i=1}^M z_i^T\mathcal F A\mathcal F^T z_i,$$

with $z_1,\ldots,z_M$ as above, independent of $D$ (Definition 3.6). It is the unit vector estimator applied to $\mathcal FA\mathcal F^T$, whose trace equals $\mathrm{trace}(A)$.

Formalization targets

Goal: Theorem 8.4

For every orthogonal seed $F$, every symmetric positive semi-definite $A$, every $\epsilon > 0$, $\delta \in (0,1)$ and every $M \ge 1$,

$$M \ \ge\ 2n^2\eta^2\epsilon^{-2}\ln(4/\delta)\ln^2(4n^2/\delta) \quad\Longrightarrow\quad T_M \text{ is an } (\epsilon,\delta)\text{-approximator of } \mathrm{trace}(A).$$

Milestones

1. Lemma 8.1. For symmetric $A$, $\mathrm{E}(U_1) = \mathrm{trace}(A)$ and $\mathrm{Var}(U_1) = n\sum_{i}A_{ii}^2 - \mathrm{trace}^2(A)$.
2. Theorem 8.2. $U_M$ is an $(\epsilon,\delta)$-approximator of $\mathrm{trace}(A)$ whenever

$$M \ge \tfrac12\epsilon^{-2}\ln(2/\delta)\,r_D^2(A).$$

3. Lemma 8.3. For $U \in \mathbb{R}^{n\times m}$ with orthonormal columns and $\delta > 0$, with probability at least $1-\delta$,

$$|(\mathcal FU)_{ij}| \le \sqrt{2\eta\ln(2mn/\delta)} \quad\text{for all } i,j.$$

4. Proof of Theorem 8.4, p. 8:13. With probability at least $1-\delta/2$ over $D$, $0 \le (\mathcal FA\mathcal F^T)_{jj} \le 2\eta\ln(4n^2/\delta)\,\mathrm{trace}(A)$ for all $j$, and hence

$$r_D(\mathcal FA\mathcal F^T) \le 2n\eta\ln(4n^2/\delta).$$

Significance

Theorem 8.2 alone shows that the unit vector estimator can need order $n^2$ samples: when the trace is concentrated on one diagonal entry, $r_D(A) = n$. Theorem 8.4 removes the dependence on $A$ entirely. For a Fourier-type seed with $\eta = \Theta(1/n)$ the bound becomes $O(\epsilon^{-2}\ln(1/\delta)\ln^2(n/\delta))$ samples for every positive semi-definite $A$, while each sample still costs about $\log_2 n$ random bits, and the $n$ bits of $D$ are drawn once. Among the estimators of the paper this is the only one with both an $A$-independent sample bound and logarithmic randomness per sample (Table I, p. 8:5). Lemma 8.3 is a standalone statement about randomized orthogonal transforms that is used well beyond trace estimation, in the analysis of subsampled randomized Hadamard transforms, sketching-based least squares, and fast Johnson–Lindenstrauss embeddings.

All results of the mission are proved in the literature: Lemma 8.3 in the cited works, the rest in the paper. As far as the platform's catalogue shows, none has a machine-checked proof. The mission produces checked statements of the paper's Section 8 with their exact constants, a probability model for random sign matrices and uniform basis-vector sampling that other randomized linear-algebra missions can reuse, and, once proved, a checked instance of Hoeffding's inequality applied to a concrete estimator.

Difficulty

The obvious argument for Theorem 8.4 applies Theorem 8.2 to $\mathcal FA\mathcal F^T$. That matrix is random, so Theorem 8.2, which is a statement about a fixed matrix, cannot be applied directly: the proof must condition on $D$, use that the samples $z_i$ are independent of $D$, and combine a failure event over $D$ with a conditional failure event over the $z_i$, each with probability at most $\delta/2$. The second difficulty is Lemma 8.3: each entry $(\mathcal FU)_{ij} = \sum_k F_{ik}D_{kk}U_{kj}$ is a Rademacher sum whose coefficient vector has squared norm at most $\eta$, and the bound needs a sub-Gaussian tail for such sums together with a union bound over all $mn$ entries. Bounding the diagonal of $\mathcal FA\mathcal F^T$ through the diagonal of $A$ alone does not work: each mixed diagonal entry depends on all entries of $A$, including the off-diagonal ones.

Formalization scope

Everything is over $\mathbb{R}$. The paper allows complex unitary seeds; since the estimator uses the transpose $\mathcal F^T$, the mission takes $F$ real orthogonal ($F^TF = I$). Matrices are Matrix (Fin n) (Fin n) ℝ, "symmetric positive semi-definite" is Matrix.PosSemidef, and $n \ge 1$ is assumed throughout. The sample spaces are explicit product measures: indices $k_1,\ldots,k_M$ uniform on Fin n with $z_i = e_{k_i}$, the diagonal of $D$ with the $n$-fold Rademacher product law, and, for $T_M$, the product of the two, which makes $D$ and the $z_i$ independent as the paper assumes implicitly. Probabilities are Measure.real. $\eta$ and $\max_i A_{ii}$ are maxima over finite nonempty index sets; $r_D$ uses real division, whose value at $\mathrm{trace}(A) = 0$ is irrelevant because a positive semi-definite matrix with zero trace is $0$. The sample-count thresholds are exactly the paper's constants.

Deviations from the page, all recorded in the items' Formalization Notes: Definition 3.4 and Theorem 8.2 are stated for positive semi-definite rather than positive definite $A$ (the proof uses only $A_{ii} \ge 0$); Table I's entry $8\epsilon^{-2}\ln(4n^2/\delta)\ln(4/\delta)$ for the mixed estimator, which disagrees with Theorem 8.4, is not used; the proof of Theorem 8.4 prints the conditional failure probability as "$\le 1-\delta/2$" where $\delta/2$ is meant, and no statement copies it; Remark 8.5 ("for some small $C$") has no pinned constant and is not stated.

A formalization in which $D$ is an arbitrary orthogonal diagonal matrix, the $z_i$ are correlated with $D$, or the law of the estimator is assumed rather than constructed would make the goal either false or a restatement of its hypotheses; the product-measure model rules this out.

Needed infrastructure: Hoeffding's inequality for bounded i.i.d. sums (in Mathlib as sub-Gaussian moment generating function bounds), a sub-Gaussian tail for Rademacher linear combinations, conditioning on one factor of a product measure, and the spectral theorem for real symmetric matrices. The Rademacher sign model and the random-mixing-matrix entry bound are reusable beyond this mission. Proofs of any milestone, and alternative arguments for Lemma 8.3, are welcome.

Selected references

• H. Avron and S. Toledo, Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix, J. ACM 58(2), Article 8, 2011. https://doi.org/10.1145/1944345.1944349
• N. Ailon and B. Chazelle, Approximate nearest neighbors and the fast Johnson–Lindenstrauss transform, STOC 2006. https://doi.org/10.1145/1132516.1132597
• H. Avron, P. Maymounkov and S. Toledo, Blendenpik: Supercharging LAPACK's least-squares solver, SIAM J. Sci. Comput. 32(3), 2010. https://doi.org/10.1137/090767911
• M. F. Hutchinson, A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines, Comm. Statist. Simulation Comput. 19(2), 1990. https://doi.org/10.1080/03610919008812866
• W. Hoeffding, Probability inequalities for sums of bounded random variables, J. Amer. Statist. Assoc. 58, 1963. https://doi.org/10.1080/01621459.1963.10500830

Motivation

Line searches are the step-size rules inside most methods for smooth unconstrained minimization $\min_{x \in \mathbb{R}^n} f(x)$: steepest descent, conjugate gradient, quasi-Newton and limited-memory methods all choose a direction $d_k$ and then a step $\alpha_k$ along it. Classical Armijo and Wolfe rules are monotone: they require $f(x_{k+1}) < f(x_k)$. Grippo, Lampariello and Lucidi (SIAM J. Numer. Anal., 1986) observed that insisting on monotone decrease can slow a method down, and proposed comparing $f(x_{k+1})$ with the maximum of the last $M$ function values instead. That max-based rule discards good function values and depends strongly on $M$, and Dai showed that R-linearly convergent iterates can violate it for every fixed memory $M$.

Zhang and Hager (SIAM J. Optim., 2004) replaced the maximum by a weighted average of all previous function values. Their averaged nonmonotone line search is used in practical codes, for example with L-BFGS and in later nonmonotone spectral and conjugate-gradient methods. This mission formalizes the paper's first main result: global convergence to stationary points for nonconvex $f$.

Setting

Let $f : \mathbb{R}^n \to \mathbb{R}$ be continuously differentiable, with gradient $g_k = \nabla f(x_k)$ at the $k$-th iterate. The Nonmonotone Line Search Algorithm (NLSA) has parameters

$$0 \le \eta_{\min} \le \eta_{\max} \le 1, \qquad 0 < \delta < \sigma < 1 < \rho, \qquad \mu > 0.$$

It maintains weights $Q_k$ and reference values $C_k$:

$$Q_0 = 1,\ C_0 = f(x_0), \qquad Q_{k+1} = \eta_k Q_k + 1, \qquad C_{k+1} = \frac{\eta_k Q_k C_k + f(x_{k+1})}{Q_{k+1}}, \qquad (1.6)$$

with $\eta_k \in [\eta_{\min}, \eta_{\max}]$ chosen at each step. The iterates are $x_{k+1} = x_k + \alpha_k d_k$, where the step $\alpha_k > 0$ satisfies one of two rules, fixed for the whole run:

• the nonmonotone Wolfe conditions

$$f(x_k + \alpha_k d_k) \le C_k + \delta \alpha_k g_k^{\mathsf T} d_k \quad (1.4), \qquad \nabla f(x_k + \alpha_k d_k) d_k \ge \sigma g_k^{\mathsf T} d_k \quad (1.5);$$

• the nonmonotone Armijo conditions: $\alpha_k = \bar\alpha_k \rho^{h_k}$, where $\bar\alpha_k > 0$ is a trial step and $h_k$ is the largest integer such that (1.4) holds and $\alpha_k \le \mu$.

The choice $\eta_k = 0$ gives $C_k = f(x_k)$, the monotone line search; $\eta_k = 1$ gives $C_k = A_k = \frac{1}{k+1}\sum_{i \le k} f(x_i)$.

The direction assumption asks for constants $c_1, c_2 > 0$ with $g_k^{\mathsf T} d_k \le -c_1\|g_k\|^2$ (2.4) and $\|d_k\| \le c_2\|g_k\|$ (2.5) for all sufficiently large $k$. The level set is $\mathcal L = \{x : f(x) \le f(x_0)\}$, and $\bar{\mathcal L}$ is the set of points whose distance to $\mathcal L$ is at most $\mu d_{\max}$, where $d_{\max} = \sup_k \|d_k\|$.

Formalization targets

Goal: Theorem 2.2

Suppose $f$ is bounded from below, $g_k^{\mathsf T} d_k \le 0$ for every $k$, the direction assumption holds, and $\nabla f$ is Lipschitz continuous on $\mathcal L$ (Wolfe rule) or on $\bar{\mathcal L}$ (Armijo rule). Then

$$\liminf_{k \to \infty} \|\nabla f(x_k)\| = 0, \qquad (2.6)$$

and if $\eta_{\max} < 1$,

$$\lim_{k \to \infty} \nabla f(x_k) = 0, \qquad (2.7)$$

so every limit of a convergent subsequence of iterates is a stationary point. No convexity is assumed.

Milestones

• Lemma 1.1: $f_k \le C_k \le A_k$ along the run, and a Wolfe step and a largest Armijo exponent exist whenever $g_k^{\mathsf T} d_k < 0$ and $f$ is bounded below.
• Eq. (1.8): $Q_{j+1} = 1 + \sum_{i=0}^{j} \prod_{m=0}^{i} \eta_{j-m} \le j+2$.
• Lemma 2.1: the lower bounds (2.1) and (2.2) on accepted Wolfe and Armijo steps.
• Eqs. (2.8)–(2.9): $f_{k+1} \le C_k - \beta\|g_k\|^2$ with the explicit constant

$$\beta = \min\left\{\frac{\delta\mu c_1}{\rho}, \frac{2\delta(1-\delta)c_1^2}{L\rho c_2^2}, \frac{\delta(1-\sigma)c_1^2}{Lc_2^2}\right\}.$$

• Eq. (2.14): $\sum_k \|g_k\|^2 / Q_{k+1} < \infty$.
• Eq. (2.15): $Q_{k+1} \le 1/(1-\eta_{\max})$ when $\eta_{\max} < 1$.
• Corollary 2.3: the analogue of Theorem 2.2 when (2.5) is replaced by the growth condition $\|d_k\|^2 \le \tau_1 + \tau_2 k$ (2.16).

Significance

Theorem 2.2 is the convergence guarantee that makes the averaged reference value $C_k$ usable in practice: any direction method whose directions are uniformly gradient-related (for example L-BFGS with bounded Hessian approximations) inherits stationarity of its limit points when combined with this line search, for every choice of the weights $\eta_k$. The monotone Wolfe and Armijo results are the special case $\eta_k \equiv 0$. The same estimates, (2.8) and (2.15), are the input to the paper's second main result, R-linear convergence for strongly convex $f$ (Theorem 3.1, a separate mission in this series).

The result is proved in the paper. It has no machine-checked proof that this mission is aware of: the platform has monotone backtracking statements for convex problems, but no nonmonotone line search, no Wolfe conditions and no Zoutendijk-type global convergence theorem for nonconvex $f$. A complete development also provides reusable Lean statements of the Wolfe and Armijo conditions and of step-size lower bounds under local Lipschitz continuity of the gradient.

Difficulty

Each step is elementary, but the argument has several places where a naive formalization fails. The Lipschitz hypothesis is only local: on $\mathcal L$ for the Wolfe rule, on the $\mu d_{\max}$-neighbourhood $\bar{\mathcal L}$ for the Armijo rule. So the proof must first show that every iterate stays in $\mathcal L$ even though $f(x_k)$ is not monotone. This needs $f_k \le C_k$ and the monotonicity of $C_k$, and then that the Armijo rule's rejected trial point $x_k + \rho\alpha_k d_k$ lies in $\bar{\mathcal L}$. The Armijo lower bound uses the maximality of the integer exponent $h_k$ and a first-order Taylor bound along a segment. The passage from (2.14) to (2.6) and (2.7) uses the two growth bounds on $Q_{k+1}$. The first gives only $\liminf$, since $\sum \|g_k\|^2/(k+2) < \infty$ does not force $g_k \to 0$.

Formalization scope

The space is EuclideanSpace ℝ (Fin n) with the Euclidean norm, $f$ is ContDiff ℝ 1, and the paper's row vector $\nabla f(x)$ acting on $d$ is the inner product of Mathlib's gradient f x with $d$. $Q_k$, $C_k$ and $A_k$ are defined by recursion from the run. A run is an infinite sequence $(x_k, d_k, \alpha_k, \eta_k)$ satisfying the update, $\eta_k \in [\eta_{\min}, \eta_{\max}]$, and the chosen rule at every $k$. The stopping test is not modelled. The Armijo exponent ranges over $\mathbb{Z}$ (it may be negative since $\rho > 1$), and "largest" is IsGreatest. $d_{\max}$ and the distance to $\mathcal L$ are computed in $[0, \infty]$, so unbounded directions give $\bar{\mathcal L} = \mathbb{R}^n$.

The following repairs and readings of the printed statements are made:

1. Theorem 2.2 and Corollary 2.3 assume $\nabla f(x_k) d_k \le 0$ for every $k$. The printed theorem constrains $d_k$ only for large $k$, but its proof needs $f(x_{k+1}) \le C_k$ at every step, which is the hypothesis of Lemma 1.1. Without it an early ascent step could leave $\mathcal L$, where nothing is assumed. The direction assumption itself stays "for all sufficiently large $k$".
2. $\liminf_k \|\nabla f(x_k)\| = 0$ is stated as "for every $\varepsilon > 0$, $\|\nabla f(x_k)\| < \varepsilon$ for infinitely many $k$".
3. The final "Hence" sentence of Theorem 2.2 is stated under $\eta_{\max} < 1$, from which it is derived.
4. In Corollary 2.3, "positive constants $\tau_1, \tau_2$" is read as $\tau_1 > 0$, $\tau_2 \ge 0$, since the corollary itself treats $\tau_2 = 0$.
5. Lemma 2.1 is stated pointwise for one iteration. Its Armijo case makes explicit the fact $f(x_k) \le C_k$ that the paper's proof invokes.

The following formalizations would make the goal trivial or empty and are ruled out:

• a (2.6) written with Lean's real liminf, which is $0$ for a divergent sequence;
• a run class in which the step rule does not constrain $\alpha_k$ (Wolfe without (1.4), Armijo without maximality of $h_k$), or which no sequence satisfies. The constant run $f \equiv 0$, $d_k = 0$ satisfies both rules, so the class is nonempty;
• assuming $\nabla f$ globally Lipschitz or the directions bounded.

Welcome contributions: proofs of the milestones in the listed order, and general lemmas on Wolfe and Armijo steps under local gradient Lipschitz continuity, which are reusable for other line-search methods.

Selected references

• H. Zhang, W. W. Hager, A Nonmonotone Line Search Technique and Its Application to Unconstrained Optimization, SIAM J. Optim. 14(4):1043–1056, 2004. https://doi.org/10.1137/S1052623403428208
• L. Grippo, F. Lampariello, S. Lucidi, A Nonmonotone Line Search Technique for Newton's Method, SIAM J. Numer. Anal. 23(4):707–716, 1986. https://doi.org/10.1137/0723046
• Y.-H. Dai, On the Nonmonotone Line Search, J. Optim. Theory Appl. 112(2):315–330, 2002. https://doi.org/10.1023/A:1013653923062

Motivation

Many computations in numerical linear algebra, statistics and computational physics need the trace of a matrix $A$ that is never formed explicitly: $A$ may be an inverse, a matrix function $f(B)$, or a product of large operators, and the only affordable access is a routine that returns $Av$ for a given vector $v$. Examples include log-determinants and the generalized cross-validation criterion in statistics, counting eigenvalues in an interval, and charge densities in electronic-structure computations. Monte Carlo trace estimators handle this setting: draw random vectors $z$, and average the quadratic forms $z^TAz$, each of which costs one matrix–vector product.

Hutchinson (1990) introduced the estimator with Rademacher vectors and computed its variance. Avron and Toledo (J. ACM 2011) replaced variance statements with sample bounds: how many samples $M$ guarantee relative error $\epsilon$ with probability $1-\delta$. Section 6 of their paper proves one such bound for an entire class of estimators at once, the normalized Rayleigh-quotient trace estimators, which contains Hutchinson's estimator and the unit vector estimator. This mission formalizes that bound (Theorem 6.1).

Setting

Let $A \in \mathbb{R}^{n\times n}$ be symmetric positive semi-definite, with eigenvalues $0 \le \lambda_1 \le \cdots \le \lambda_n$ and rank $\mathrm{rank}(A)$. Write $\lambda_n$ for the largest eigenvalue and

$$\kappa_f(A) = \frac{\text{largest nonzero eigenvalue of }A}{\text{smallest nonzero eigenvalue of }A},$$

defined for $A \ne 0$; it is the condition number of $A$ on its range.

A normalized Rayleigh-quotient trace estimator of $A$ with $M$ samples is

$$R_M = \frac1M\sum_{i=1}^M z_i^TAz_i,$$

where $z_1,\ldots,z_M$ are independent random vectors in $\mathbb{R}^n$ with $z_i^Tz_i = n$ and $\mathrm{E}(z_i^TAz_i) = \mathrm{trace}(A)$ for each $i$ (Definition 3.2). The vectors need not be identically distributed. Hutchinson's vectors ($\pm1$ entries, i.i.d. uniform) and the vectors $\sqrt n\,e_k$ with $k$ uniform are instances.

A random variable $T$ is an $(\epsilon,\delta)$-approximator of $\mathrm{trace}(A)$ if

$$\Pr\bigl(|T-\mathrm{trace}(A)| \le \epsilon\,\mathrm{trace}(A)\bigr) \ge 1-\delta$$

(Definition 4.1).

Formalization targets

Goal: Theorem 6.1, in the form its proof establishes

For every nonzero symmetric positive semi-definite $A$, every $\epsilon>0$, $\delta\in(0,1)$, and every normalized Rayleigh-quotient estimator $R_M$ of $A$,

$$M \ \ge\ \frac{\ln(2/\delta)\cdot n^2\,\kappa_f^2(A)}{2\,\mathrm{rank}^2(A)\,\epsilon^2} \quad\Longrightarrow\quad R_M \text{ is an } (\epsilon,\delta)\text{-approximator of } \mathrm{trace}(A).$$

Milestones (the displayed steps of the proof, p. 8:10)

1. $\mathrm{trace}(A)\,\kappa_f(A) \ge \mathrm{rank}(A)\,\lambda_n$.
2. For every $z$ with $z^Tz = n$: $\ 0 \le z^TAz \le \lambda_n z^Tz = n\lambda_n \le \frac{n}{\mathrm{rank}(A)}\mathrm{trace}(A)\,\kappa_f(A)$.
3. For every $t>0$:

$$\Pr(|R_M-\mathrm{trace}(A)| \ge t) \le 2\exp\left(-\frac{2M^2\mathrm{rank}^2(A)t^2}{M n^2\mathrm{trace}^2(A)\kappa_f^2(A)}\right).$$

4. For every $\epsilon>0$:

$$\Pr(|R_M-\mathrm{trace}(A)| \ge \epsilon\,\mathrm{trace}(A)) \le 2\exp\left(-\frac{2M\mathrm{rank}^2(A)\epsilon^2}{n^2\kappa_f^2(A)}\right).$$

Significance

The result is distribution-free within the class: it needs only normalization and unbiasedness, so it covers Hutchinson's estimator, the unit vector estimator and any future normalized scheme with one argument. For well-conditioned matrices of full or nearly full rank the required number of samples is $O(\epsilon^{-2}\ln(1/\delta))$, independent of $n$. For ill-conditioned matrices the bound degrades with $\kappa_f^2(A)$, which is the reason the paper proves sharper estimator-specific bounds in Sections 7 and 8; Theorem 6.1 is the baseline those results are compared against (Table I, p. 8:5).

The theorem is proved in the paper. No machine-checked version of it, or of any sample bound for trace estimators, is known to exist. The formalization adds a precise statement of the class of estimators on a general probability space, a corrected threshold (see below), and a Lean development that connects Mathlib's spectral theorem for symmetric matrices with its Hoeffding inequality for independent bounded variables.

Difficulty

Each step is short on paper; the work is in the interfaces. The eigenvalue inequality of milestone 1 requires relating the number of nonzero eigenvalues (with multiplicity) to $\mathrm{rank}(A)$ and handling the maximum and minimum over the nonzero spectrum. Milestone 2 is the Rayleigh-quotient bound $z^TAz \le \lambda_n z^Tz$, which is a consequence of the spectral decomposition rather than a one-line identity. Milestone 3 applies Hoeffding's inequality to summands that are bounded only almost surely, are not identically distributed, and whose mean is fixed by hypothesis rather than computed; the two-sided bound must be assembled from two one-sided tails, and the event $|R_M - \mathrm{trace}(A)| \ge t$ must be rescaled to a statement about the sum $\sum_i z_i^TAz_i$. A naive attempt that fixes a particular distribution for the $z_i$ (Rademacher, say) proves a different, narrower theorem and does not settle the goal.

Formalization scope

• Matrices and spectrum. $A$ is Matrix (Fin n) (Fin n) ℝ with A.PosSemidef and A ≠ 0; eigenvalues are Mathlib's IsHermitian.eigenvalues. $\lambda_n$ is lambdaMax (the maximum eigenvalue) and $\kappa_f(A)$ is kappaF (maximum over minimum of the finite set of nonzero eigenvalues); both are Finset.max'/min'/sup' of nonempty finite sets. $\kappa_f(0)$ is a placeholder, and every statement assumes $A \ne 0$.
• Probability model. A general probability space $(\Omega, P)$ and random vectors z : Fin M → Ω → Fin n → ℝ satisfying IsNormalizedRayleighSample P A z: each $z_i$ measurable, the family mutually independent (iIndepFun), $z_i^Tz_i = n$ almost surely, and $\int z_i^TAz_i\,dP = \mathrm{trace}(A)$. The estimator is universally quantified over this class. Unbiasedness is required for the given $A$ only, as on the page. Probabilities are P.real of events; $M \ge 1$ is a natural number and $1/M$ is (M : ℝ)⁻¹.
• Correction of the printed statement. Theorem 6.1 and Table I print the threshold $\tfrac12\epsilon^{-2}n^{-2}\mathrm{rank}^2(A)\ln(2/\delta)\kappa_f^2(A)$. The last display of the proof gives $\ln(2/\delta)\,n^2\kappa_f^2(A)/(2\,\mathrm{rank}^2(A)\epsilon^2)$, with the exponents of $n$ and $\mathrm{rank}(A)$ swapped. The printed version is false: for $n=2$, $A=e_1e_1^T$, $z=\sqrt2\,e_k$ with $k$ uniform and $\epsilon=\delta=1/2$ it admits $M=1$, while $R_1\in\{0,2\}$. The goal states the proof's threshold.
• Indexing slip. The proof writes $0=\lambda_1=\cdots=\lambda_k$ with $k=n-\mathrm{rank}(A)+1$ and $\kappa_f(A)=\lambda_n/\lambda_k$, which would make $\lambda_k=0$; milestone 1 states the inequality with $\kappa_f$ as defined, not the indexing.
• Non-trivialization. The hypotheses of the class are satisfiable (by the unit vector and Hutchinson estimators), $A\ne0$ excludes the degenerate $\kappa_f(0)$, and all divisions in the statements have positive denominators, so no statement holds vacuously or through a junk value.
• Reusable parts. A proof of milestone 2 is a general Rayleigh-quotient bound for symmetric matrices; milestone 3 is a two-sided Hoeffding bound for independent, almost surely bounded, non-identically distributed summands, useful well beyond this mission. Contributions welcome: proofs of any milestone, and a sorry-free instance showing a concrete estimator (e.g. $\sqrt n\,e_k$) satisfies IsNormalizedRayleighSample.

Selected references

• H. Avron and S. Toledo, Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix, Journal of the ACM 58(2), Article 8, 2011. https://doi.org/10.1145/1944345.1944349
• M. F. Hutchinson, A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines, Communications in Statistics – Simulation and Computation 19(2), 433–450, 1990. https://doi.org/10.1080/03610919008812866
• W. Hoeffding, Probability inequalities for sums of bounded random variables, Journal of the American Statistical Association 58(301), 13–30, 1963. https://doi.org/10.1080/01621459.1963.10500830