Numerical Analysis

Motivation

Many problems in signal processing, statistics and machine learning take the form

$$\min_{x\in\mathcal H}\ \Psi(x)+\Phi(x),$$

where $\Phi$ is smooth and convex (a data-fit term) and $\Psi$ is convex but possibly nonsmooth or infinite-valued (an $\ell^1$ penalty, the indicator function of a constraint set). The forward-backward method alternates a gradient step on $\Phi$ with a proximal step on $\Psi$ and reduces the objective gap at rate $O(k^{-1})$ after $k$ iterations. Combining it with Nesterov's extrapolation scheme gives the accelerated forward-backward method, best known as FISTA, which improves the guaranteed rate to $O(k^{-2})$. FISTA and its variants are standard solvers for sparse recovery and image reconstruction.

Timeline:

• 1983. Nesterov introduces the extrapolation scheme for smooth convex minimization, with an $O(k^{-2})$ rate for function values (Nesterov 1983).
• 2009. Beck and Teboulle extend it to the composite problem above (FISTA), with the $O(k^{-2})$ rate (doi:10.1137/080716542).
• 2014. Su, Boyd and Candès read the scheme as a discretization of the ODE $\ddot x+\frac{\alpha}{t}\dot x+\nabla\Theta(x)=0$ (arXiv:1503.01243).
• 2014–2015. Chambolle and Dossal (doi:10.1007/s10957-015-0746-4) and, independently, Attouch, Chbani, Peypouquet and Redont (arXiv:1507.04782) prove weak convergence of the iterates for the variant with parameter $\alpha>3$. Before that, convergence of the iterates had been open for about two decades.
• 2015. May shows that for $\alpha>3$ the continuous-time gap is $o(t^{-2})$ (arXiv:1509.05598).
• 2016. Attouch and Peypouquet prove the discrete analogue: for $\alpha>3$ the function gap of the algorithm is $o(k^{-2})$ and the velocity is $o(k^{-1})$ (arXiv:1510.08740, SIAM J. Optim. 26(3):1824–1834).

Setting

Let $\mathcal H$ be a real Hilbert space with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$. Let:

1. $\Psi:\mathcal H\to\mathbb R\cup\{+\infty\}$ be proper (finite somewhere), lower semicontinuous and convex;
2. $\Phi:\mathcal H\to\mathbb R$ be convex and continuously differentiable, with gradient $\nabla\Phi$ Lipschitz continuous with constant $L$;
3. $\Theta=\Psi+\Phi$, and assume the set of minimizers $S=\operatorname{argmin}\Theta$ is nonempty.

For $s>0$, the proximal map $\operatorname{prox}_{s\Psi}(z)$ is the unique minimizer of $u\mapsto\Psi(u)+\frac1{2s}\|u-z\|^2$. Given $\alpha>0$, a step size $0<s<1/L$ and starting points $x_0,x_1$, algorithm (2) generates

$$y_k=x_k+\frac{k-1}{k+\alpha-1}(x_k-x_{k-1}),\qquad x_{k+1}=\operatorname{prox}_{s\Psi}\big(y_k-s\nabla\Phi(y_k)\big),\qquad k\ge1.$$

The common choice is $\alpha=3$; this mission concerns $\alpha>3$.

The proofs use the operator $G_s(y)=\frac1s\big(y-\operatorname{prox}_{s\Psi}(y-s\nabla\Phi(y))\big)$, the sequence $z_k=x_k+\frac{k-1}{\alpha-1}(x_k-x_{k-1})$, a minimizer $x^*$, and the quantity

$$\mathcal E(k)=\frac{2s}{\alpha-1}(k+\alpha-2)^2\big(\Theta(x_k)-\Theta(x^*)\big)+(\alpha-1)\|z_k-x^*\|^2 .$$

Write $\theta_k=\Theta(x_k)-\Theta(x^*)$ and $d_k=\frac1{2s}\|x_{k+1}-x_k\|^2$.

Formalization targets

Goal: Theorem 1 (p. 2)

For $\alpha>3$ and $0<s<1/L$,

$$\lim_{k\to\infty}k^2\big(\Theta(x_k)-\min\Theta\big)=0\qquad\text{and}\qquad\lim_{k\to\infty}k\,\|x_{k+1}-x_k\|=0.$$

The statement fixes no constants, so it is unaffected by later improvements of the explicit bounds below.

Milestones (in attack order)

1. (9), p. 2. If $sL\le1$, then for all $x,y$: $\Theta(y-sG_s(y))\le\Theta(x)+\langle G_s(y),y-x\rangle-\frac s2\|G_s(y)\|^2$.
2. (13), p. 3. For $\alpha\ge3$ and $k\ge1$: $\mathcal E(k+1)+2s\frac{\alpha-3}{\alpha-1}k\,\theta_k\le\mathcal E(k)$.
3. Fact 1, p. 3. $(\mathcal E(k))$ is nonincreasing and has a finite limit.
4. Fact 2, p. 3. $\theta_k\le\frac{(\alpha-1)\mathcal E(1)}{2s(k+\alpha-2)^2}$ and $\|z_k-x^*\|^2\le\frac{\mathcal E(1)}{\alpha-1}$ for $k\ge1$.
5. Fact 3, p. 3. For $\alpha>3$: $\sum_{k\ge1}k\,\theta_k\le\frac{(\alpha-1)\mathcal E(1)}{2s(\alpha-3)}$.
6. (14), p. 4. $\Theta(x_{k+1})+d_k\le\Theta(x_k)+\frac{(k-1)^2}{(k+\alpha-1)^2}d_{k-1}$ for $k\ge1$.
7. Fact 4, p. 4. For $\alpha>3$: $\sum_{k\ge1}k\,d_k\le\frac{\alpha(3\alpha-5)\mathcal E(1)}{4s(\alpha-1)(\alpha-3)}$.
8. Lemma 2, p. 4. For $\alpha>3$, $\lim_k\big[k^2d_k+(k+1)^2\theta_{k+1}\big]$ exists and is finite.

Significance

The result. The $O(k^{-2})$ rate of FISTA is often quoted as optimal for first-order methods. Theorem 1 shows that for $\alpha>3$ the worst-case rate along every single run is strictly better, $o(k^{-2})$, and that the steps $\|x_{k+1}-x_k\|$ decay faster than $1/k$. No better power is possible: by Attouch et al., Example 2.13 there is no $p>2$ with an $O(k^{-p})$ rate for every $\Phi$ and $\Psi$. The intermediate estimates (Facts 2–4) give explicit, quantitative bounds that are reused in the analysis of inexact and perturbed variants (Theorem 4 of the paper) and in mission II of this series, which proves weak convergence of the iterates.

Formalizing it. The result is proved on paper. As far as we know, no machine-checked proof of the $O(k^{-2})$ rate of FISTA in this Hilbert-space, extended-valued setting exists in Mathlib or on this platform, and neither does the $o(k^{-2})$ refinement. A complete development would provide reusable statements about proximal-gradient steps for functions valued in $\mathbb R\cup\{+\infty\}$. The page also has a factor slip in Fact 4 and Lemma 2 (see Formalization scope); checking it mechanically is one of the outputs.

Difficulty

The $O(k^{-2})$ bound follows from the monotonicity of $\mathcal E$ alone. That argument cannot give $o(k^{-2})$: it controls $k^2\theta_k$ only by the constant $\mathcal E(1)$, and summability of $k\theta_k$ (Fact 3) gives a decay of $k^2\theta_k$ only along a subsequence, not of the whole sequence. The missing ingredient is convergence of a weighted combination of function gaps and velocities, which is Lemma 2. The bookkeeping is delicate: every step carries explicit coefficients in $k$ and $\alpha$, and the inequalities must hold in the extended reals because $\Psi$ may be $+\infty$ at the starting point.

Formalization scope

• Space and functions. $\mathcal H$ is a real inner product space that is complete. $\Psi$ takes values in EReal and satisfies the published predicate IsProperClosedConvex (never $-\infty$, finite somewhere, lower semicontinuous, convex epigraph). $\Phi:\mathcal H\to\mathbb R$ is ConvexOn and ContDiff ℝ 1, and gradient Φ is LipschitzWith L for some L : ℝ≥0.
• Step size. $0<s<1/L$ is written 0 < s and s * L < 1, so that $L=0$ is allowed (s < 1 / L would be unsatisfiable when $L=0$). Display (9) uses the page's $s\le1/L$, written s * L ≤ 1.
• Proximal map. It is a map P with the published predicate IsProx s Ψ P, which determines $P=\operatorname{prox}_{s\Psi}$ for proper closed convex $\Psi$ and $s>0$.
• The run. It is required to follow (2) for every $k\ge1$, with $x_0,x_1$ arbitrary; the coefficient of $x_0$ vanishes at $k=1$.
• Extended reals. $\Theta$, $\mathcal E$ and the function-value statements live in EReal. No value is ever converted to $\mathbb R$ with toReal, which would turn $+\infty$ into $0$. "The limit exists" (Fact 1, Lemma 2) means convergence to a real number, because in EReal every monotone sequence converges. Infinite series of nonnegative terms are stated as bounds on every partial sum. $\min\Theta$ in the goal is ⨅ y, Θ y together with the hypothesis that a minimizer exists.
• Corrected slips. Fact 2 is stated with $\mathcal E(1)$ for $k\ge1$; the page writes $\mathcal E(0)$ for $k\ge0$, which needs an iterate $x_{-1}$. Fact 4 is stated for $\sum k\,d_k$; the page prints $\sum k\|x_{k+1}-x_k\|^2$ with the same constant, which is off by the factor $2s$ and false in general. Lemma 2 is stated for the bracket $k^2d_k+(k+1)^2\theta_{k+1}$ of its proof (16).
• Ruled out. The goal assumes nothing about $\mathcal E$, $z_k$, $\theta_k$ or $d_k$. A statement of the first limit for toReal values, or with an extended-real "limit exists", would be trivially weaker and is not what is asked.
• Welcome contributions. Proofs of any milestone; general facts about proximal maps of EReal-valued convex functions and about the descent property of $L$-smooth convex functions, which are reusable well beyond this mission.

Selected references

• H. Attouch, J. Peypouquet, The Rate of Convergence of Nesterov's Accelerated Forward-Backward Method is Actually Faster than $1/k^2$, SIAM J. Optim. 26(3):1824–1834, 2016. arXiv:1510.08740v4, doi:10.1137/15M1046095
• H. Attouch, Z. Chbani, J. Peypouquet, P. Redont, Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity, Math. Program. 168:123–175, 2018. arXiv:1507.04782
• A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imaging Sci. 2(1):183–202, 2009. doi:10.1137/080716542
• A. Chambolle, C. Dossal, On the convergence of the iterates of the "fast iterative shrinkage/thresholding algorithm", J. Optim. Theory Appl. 166:968–982, 2015. doi:10.1007/s10957-015-0746-4
• R. May, Asymptotic for a second order evolution equation with convex potential and vanishing damping term, 2015. arXiv:1509.05598
• Y. Nesterov, A method of solving a convex programming problem with convergence rate $O(1/k^2)$, Soviet Math. Dokl. 27:372–376, 1983. mathnet
• W. Su, S. Boyd, E. J. Candès, A differential equation for modeling Nesterov's accelerated gradient method: theory and insights, J. Mach. Learn. Res. 17(153):1–43, 2016. arXiv:1503.01243

Motivation

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

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

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

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

Setting

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

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

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

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

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

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

Formalization targets

Goal: Theorem 3.9, one stage squares the error

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

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

Milestones inside the proof

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

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

Companion result: Theorem 3.8

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

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

Significance

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

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

Difficulty

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

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

Formalization scope

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

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

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

Selected references

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

Motivation

The ellipsoid method is the algorithm that settled the polynomial-time solvability of linear programming (Khachiyan, 1979) and that underlies the equivalence of separation and optimization in combinatorial optimization (Grötschel, Lovász and Schrijver, 1981). Its origin is in nonsmooth convex optimization. In 1976 Yudin and Nemirovskii proposed a modified method of centered sections that localizes an optimum inside a sequence of ellipsoids, and in 1977 N. Z. Shor observed independently that the same scheme is a subgradient method with space dilation along the gradient, the family of methods he had developed since 1969. Section 3.8 of Shor's monograph Minimization Methods for Non-Differentiable Functions (Springer 1985) presents the method in this second form and proves its basic localization property.

Timeline:

• 1965. A. Yu. Levin proposes the method of centered sections (cuts through the center of gravity of a polyhedron); each cut removes at least a fixed fraction of the volume, but computing centers of gravity is impractical for $n > 3$.
• 1969–1972. Shor introduces subgradient methods with space dilation along the gradient (SDG methods).
• 1976. Yudin and Nemirovskii replace the polyhedron by a minimal ellipsoid containing a half-ellipsoid, obtaining a geometric volume decrease depending only on the dimension ([Yudin–Nemirovskii 1976]).
• 1977. Shor shows that the same method is an SDG algorithm with coefficient $\beta = \sqrt{(n-1)/(n+1)}$ ([Shor 1977]).
• 1979. Khachiyan applies the method to linear inequalities with integer data, obtaining the first polynomial-time algorithm for linear programming ([Khachiyan 1979]).

Setting

Let $E_n$ be $n$-dimensional Euclidean space with inner product $(x, y)$, and $n > 1$. For a unit vector $\xi$ and a number $\alpha$, the operator of space dilation along $\xi$ with coefficient $\alpha$ is $R_\alpha(\xi) = I + (\alpha - 1)\xi\xi^T$: it multiplies the component of a vector along $\xi$ by $\alpha$ and leaves the orthogonal component unchanged.

Let $g : E_n \to E_n$ be a vector field, not necessarily continuous. The problem is to find a point $x^*$ with

$$(g(x), x - x^*) \ge 0 \quad \text{for all } x \in E_n,$$

where it is known that such an $x^*$ exists in the closed ball $S(x_0, R)$ of radius $R > 0$ about a given point $x_0$. Put $\beta = \sqrt{(n-1)/(n+1)}$ and $r = n/\sqrt{n^2-1}$. The algorithm (3.57)–(3.60) starts from $x_0$, $B_0 = I_n$, $h_0 = R/(n+1)$, and at iteration $k+1$ stops if $g(x_k) = 0$, and otherwise sets

$$\xi_k = \frac{B_k^T g(x_k)}{\|B_k^T g(x_k)\|}, \quad x_{k+1} = x_k - h_k B_k \xi_k, \quad B_{k+1} = B_k R_\beta(\xi_k), \quad h_{k+1} = r h_k .$$

With $A_k = B_k^{-1}$, the localizing ellipsoid is $\Phi_k = \{x : \|A_k(x - x_k)\| \le (n+1)h_k\}$, and the dimension-dependent ratio is

$$q_n = \sqrt{\frac{n-1}{n+1}}\left(\frac{n}{\sqrt{n^2-1}}\right)^n < 1 .$$

Three problems produce such a field: minimizing a convex $f$ on a ball (the field (3.62), a subgradient inside the ball and the outward radial direction outside); the convex program $\min f_0$ s.t. $f_i \le 0$ (the field (3.65), a subgradient of the objective at feasible points and of a most violated constraint otherwise); and a convex–concave saddle point problem (the field $\{g_f^x, -g_f^y\}$).

Formalization targets

Goal: Theorem 3.14 (p. 86)

For every $k$,

$$\|A_k(x_k - x^*)\| \le h_k (n+1), \tag{3.61}$$

that is, $x^* \in \Phi_k$. The statement holds for every field $g$ satisfying the monotonicity condition at $x^*$; nothing about continuity or convexity is assumed.

Milestones

1. Eq. (3.4): $\|R_\alpha(\xi)x\| = \sqrt{\|x\|^2 + (\alpha^2-1)(x,\xi)^2}$ for unit $\xi$.
2. Volume of $\Phi_k$ (p. 87): $(n+1)h_k = R r^k$ and $v(\Phi_k) = v_0 R^n r^{nk}/\det A_k$, $v_0$ the volume of the unit ball.
3. Volume ratio (p. 87–88): $v(\Phi_{k+1}) = q_n\, v(\Phi_k)$ with $q_n < 1$, the volumes being positive and finite.
4. Eq. (3.62): the ball field satisfies $(g(x), x - x^*) \ge 0$.
5. Eq. (3.65): the convex-programming field satisfies $(g(x), x - x^*) \ge 0$.
6. Saddle point field (p. 90): $(g(z), z - z^*) \ge 0$.

Significance

The result. Theorem 3.14 with the volume identity says that after $k$ steps the solution is confined to an ellipsoid of volume $q_n^k$ times that of the initial ball, for any field of the above kind. Milestones 4–6 turn this into localization guarantees for constrained convex minimization, general convex programming and convex–concave saddle points, with a rate that depends only on the dimension. The same localization underlies the complexity bounds of the ellipsoid method for linear programming and the polynomial equivalence of separation and optimization.

Formalizing it. The results are classical and proved in the book. No machine-checked proof of the space-dilation form of the method is known to exist. The platform already contains a proved version of the Bertsimas–Tsitsiklis form (LinearOptimization.ellipsoid_update_halfspace_subset, LinearOptimization.ellipsoid_update_volume_lt: a half-ellipsoid $E(z, D) \cap \{a^Tx \ge a^Tz\}$ is covered by an updated ellipsoid whose volume is smaller by a factor below $e^{-1/(2(n+1))}$), and Khachiyan's feasibility algorithm (SmaleNinth.khachiyan_ellipsoid_decides). Those statements are about a center/shape-matrix update and give a volume inequality; this mission is about the iterates of Shor's matrix recursion $B_{k+1} = B_k R_\beta(\xi_k)$ and the exact ratio $q_n$. Relating the two parametrizations ($D_k = (n+1)^2 h_k^2 B_k B_k^T$) is a welcome side result.

Difficulty

The obvious approach, tracking the ellipsoid through the center and shape matrix and invoking a minimum-volume covering argument, is not what the algorithm computes: here the iterate is updated through the factor $B_k$ and the stepsize $h_k$ is fixed in advance, independent of the field, so the induction must be carried out in the transformed coordinates $z_k = A_k(x_k - x^*)$ in which the ellipsoid is a ball. The difficulty is that the monotonicity condition gives only the sign of one inner product, $(z_k, \xi_k) \ge 0$, while the norm of $z_{k+1}$ depends on both $(z_k,\xi_k)$ and $\|z_k\|$; the constants $\beta$ and $r$ are exactly those for which the resulting quadratic estimate closes. For the volume identity, the main technical step is computing $\det R_\beta(\xi) = \beta$ and the Lebesgue measure of a linear image of a ball in EuclideanSpace.

Formalization scope

• $E_n$ is EuclideanSpace ℝ (Fin n); matrices act through Matrix.toEuclideanLin; $B_k^*$ is the transpose. $R_\alpha(\xi)$ is the matrix $I + (\alpha-1)\xi\xi^T$ (the book's property 10); for unit $\xi$ this is the operator of the book's definition.
• The algorithm is the definition ellipsoidMethod g R x₀ : ℕ → EllState n, with state $(x_k, B_k, h_k)$, $B_0 = I$, $h_0 = R/(n+1)$. If $g(x_k) = 0$ the state is repeated from then on (the book stops); the normalization in (3.57) is only performed when $g(x_k) \ne 0$. $A_k$ is the matrix inverse of $B_k$, which is nonsingular.
• Theorem 3.14 is stated for $n > 1$, $R > 0$, $x^*$ with $\|x_0 - x^*\| \le R$ and $(g(x), x - x^*) \ge 0$ for all $x$, for every $k$. The book's additional assumption that $g(x) \ne 0$ for $x \ne x^*$ is not used by its proof and is omitted. The iterates are those of the recursion; a statement about an arbitrary ellipsoid containing $x^*$, or about an arbitrary invertible matrix in place of $B_k$, would not be this theorem and is ruled out by the definitions.
• Volumes are Lebesgue measure in ℝ≥0∞. The ellipsoid is given for an arbitrary center (the book writes $x^*$ in one place and $x_k$ in another; the volume is the same). The volume formula requires the first $k$ iterations to have been performed ($g(x_j) \ne 0$, $j < k$); the ratio requires iteration $k+1$ to be performed.
• The printed chain on p. 87 has misprints (exponents $2$ and $1$ on $n/\sqrt{n^2-1}$ where $n$ is meant, and $(n-1)/(n+1)$ for $\sqrt{(n-1)/(n+1)}$); the statement follows the value of $q_n$ given on p. 88. The estimate for $x \notin S(x_0,R)$ before (3.62) has a sign misprint; only the conclusion is stated.
• Needed infrastructure: determinant of a rank-one perturbation of the identity ($\det(I + c\,\xi\xi^T) = 1 + c\|\xi\|^2$, available in Mathlib as the matrix determinant lemma), the measure of a linear image (MeasureTheory.Measure.addHaar_image_linearMap), and elementary real inequalities for $q_n < 1$. The space-dilation lemmas are reusable in the other Shor missions on SDG methods and the $r$-algorithm.

Selected references

• N. Z. Shor, Minimization Methods for Non-Differentiable Functions, Springer Series in Computational Mathematics 3, Springer, 1985, §3.8. https://doi.org/10.1007/978-3-642-82118-9
• D. B. Yudin and A. S. Nemirovskii, Informational complexity and efficient methods for the solution of convex extremal problems, Ekonomika i Matematicheskie Metody 12 (1976), 357–369 (English translation: Matekon 13 (1977), 25–45).
• N. Z. Shor, Cut-off method with space extension in convex programming problems, Cybernetics 13 (1977), 94–96.
• L. G. Khachiyan, A polynomial algorithm in linear programming, Soviet Mathematics Doklady 20 (1979), 191–194.
• R. G. Bland, D. Goldfarb and M. J. Todd, The ellipsoid method: a survey, Operations Research 29 (1981), 1039–1091. https://doi.org/10.1287/opre.29.6.1039
• M. Grötschel, L. Lovász and A. Schrijver, Geometric Algorithms and Combinatorial Optimization, Springer, 1988. https://doi.org/10.1007/978-3-642-97881-4

Motivation

Linear programming is usually described as "solvable in polynomial time", but that statement is about rational inputs given exactly. In Smale's list of problems for the 21st century (Smale 1998), Problem 9 asks for a polynomial-time algorithm over the reals deciding the feasibility of $Ax \ge y$, and Smale explicitly calls for "models which process approximate inputs and which permit round-off computations". Real data such as $\sqrt 2$, entries of a discrete cosine transform, or even $1/3$ in floating point can only be accessed approximately.

Bastounis, Hansen and Vlačić pose the extended Smale's 9th problem: in a model where the algorithm can only query approximations of the input to any requested accuracy, can one compute minimisers of linear programming, basis pursuit and Lasso to $K$ correct digits? Their Main Theorem I (Theorem 3.4) shows that the answer depends on $K$ in a sharp way: for a suitable class of well-conditioned, bounded inputs, no algorithm at all (not only no efficient one) produces $K$ correct digits, while $K-1$ digits are computable (but not in bounded time) and $K-2$ digits are computable in polynomial time.

This mission targets the first, impossibility, half of Theorem 3.4(i) for linear programming.

Setting

Linear program. For $A \in \mathbb R^{m\times N}$, $y\in\mathbb R^m$ and $c = \mathbf 1_N=(1,\dots,1)$, the solution set is

$$\Xi(y,A) = \operatorname*{argmin}_{x\in\mathbb R^N}\ \langle x, c\rangle \quad\text{subject to}\quad Ax = y,\ x\ge 0 .$$

It is a subset of $M_N = \mathbb R^N$ with the $\ell^p$ norm, $p\in[1,\infty]$. An input is a pair $\iota = (y,A)$, and the evaluations of $\iota$ are its coordinates $y_i$ and entries $A_{ij}$.

Extended model ($\Delta_1$-information). Let $D_n = \{k2^{-n} : k\in\mathbb Z\}$. An oracle representation of $\iota$ is a family $\tilde\iota = (\tilde\iota_{j,n})$, indexed by evaluations $j$ and accuracies $n = 1,2,\dots$, with $\tilde\iota_{j,n}\in D_n + iD_n$ and $|\tilde\iota_{j,n} - f_j(\iota)|\le 2^{-n}$. An algorithm must succeed on every oracle representation of every input.

General algorithm. To make impossibility results independent of the machine model, the paper uses general algorithms (Definition 9.3): a map $\Gamma$ from inputs to $M\cup\{\mathrm{NH}\}$ ($\mathrm{NH}$ = no output) together with a nonempty set $\Lambda_\Gamma(\iota)$ of evaluations read on $\iota$. This set is finite whenever $\Gamma$ halts. The output is determined by the values read, and any input that agrees on those values reads the same set. Turing machines and BSS machines with an oracle are special cases; general algorithms can even solve the halting problem.

Error and breakdown epsilon. The error is $\operatorname{dist}(\Gamma(\iota),\Xi(\iota)) = \inf_{\xi\in\Xi(\iota)} d(\Gamma(\iota),\xi)$, with distance $\infty$ from $\mathrm{NH}$. The strong breakdown epsilon $\varepsilon_B^s$ is the supremum of all $\varepsilon\ge 0$ such that every general algorithm has error $>\varepsilon$ on some input (Definition 9.17).

Formalization targets

Goal: Theorem 3.4(i), deterministic part, for LP

For every integer $K\ge1$, all dimensions $4\le m<N$ and every $p\in[1,\infty]$ there is a nonempty class $\Omega_{m,N}$ of inputs $(y,A)$ with nonempty solution sets, $\|y\|_\infty\le 2$ and $\|A\|_{\max}=1$, such that

$$\neg\,\exists\,\Gamma\ \text{general algorithm on oracle representations}:\quad \forall\,\tilde\iota,\ \ \operatorname{dist}_{\ell^p}\big(\Gamma(\tilde\iota),\,\Xi(\iota)\big)\le 10^{-K}.$$

Milestones

1. Lemma 11.1: the explicit solution sets of the LP inputs $(y_1e_1, A(\alpha,\beta,m,N))$.
2. Proposition 10.5 (ii), deterministic part: two input sequences that converge in evaluation to a common input and whose solutions stay $\kappa$ apart force $\varepsilon_B^s\ge\kappa/2$ for a suitable choice of $\Delta_1$-information.
3. §9.6, (i) ⇒ (ii): a lower bound on $\varepsilon_B^s$ for one specific $\Delta_1$-information transfers to the problem with all oracle representations.
4. Proposition 9.32 (i) (deterministic consequence via Proposition 10.1): $\varepsilon_B^s>10^{-K}$ for LP on a suitable $\Omega_{m,N}$.

Significance

The theorem shows that for LP with inexact input, being non-computable in Turing's sense does not rule out a finer complexity theory. The paper builds a "$K$ / $K-1$ / $K-2$ digits" classification on this. It also explains why established solvers can return wrong answers with a success flag on small, well-conditioned LPs (§4 of the paper), and it bears on computer-assisted proofs that rely on inexact LP, such as the Flyspeck proof of the Kepler conjecture.

The result is proved on paper. As far as the proposer knows, it has not been machine-checked. This mission formalizes the deterministic impossibility part for LP and puts in place reusable infrastructure: general algorithms, breakdown epsilons and $\Delta_1$-information. That infrastructure is the base for later missions on the randomised parts of Theorem 3.4(i)–(ii), the weak breakdown epsilon (iii), the exit-flag theorem (Theorem 5.1), and basis pursuit and Lasso.

Difficulty

The obvious objection is that LP is in P for rational inputs, so some rounding scheme ought to work. It fails because an algorithm must halt after reading finitely many approximations. Two inputs that agree to that accuracy but have minimisers far apart then receive the same output. Setting this up needs a notion of algorithm strong enough to cover every computational model, a precise $\Delta_1$-information model in which the adversary controls the approximations, and explicit LP geometry in which an arbitrarily small perturbation of $A$ moves the minimiser by a fixed amount.

Formalization scope

• Inputs are $(y,A)\in(\mathrm{Fin}\,m\to\mathbb R)\times\mathrm{Matrix}(\mathrm{Fin}\,m)(\mathrm{Fin}\,N)\,\mathbb R$. Evaluations are complex-valued, as in Definition 9.2. Outputs lie in PiLp p (Fin N → ℝ).
• A general algorithm is a structure with an output run : Ω → Option M (none = NH) and a read set queried, satisfying the axioms (i)–(iii) of Definition 9.3.
• Errors take values in $[0,\infty]$ (ℝ≥0∞), and the error of NH is $\infty$. The infimum over an empty solution set is $\infty$. The goal also requires nonempty solution sets, so no junk value enters.
• Oracle accuracies are indexed by $n\in\{1,2,\dots\}$ (ℕ+). An oracle input is stored as a pair (input, oracle family), and algorithms can read only the oracle family.
• Out of scope: randomised algorithms, the positive statements (iii)–(iv), runtime, and the condition-number bounds $\mathrm{Cond}(AA^*)\le3.2$, $C_{FP}\le4$, $\mathrm{Cond}(\Xi)\le179$.

Selected references

• A. Bastounis, A. C. Hansen, V. Vlačić, The extended Smale's 9th problem — On computational barriers and paradoxes in estimation, regularisation, computer-assisted proofs, and learning, preprint (2021).
• S. Smale, Mathematical problems for the next century, Math. Intelligencer 20 (1998). https://doi.org/10.1007/BF03025291

Motivation

Many computations need the trace of a matrix $A$ that is never formed explicitly and can only be applied to vectors: the trace of a matrix function such as $\mathrm{trace}(A^{-1})$ or $\log\det A = \mathrm{trace}(\log A)$ in statistics and lattice QCD, the Frobenius norm $\|B\|_F^2 = \mathrm{trace}(B^TB)$ of an operator, or the number of triangles of a graph. The standard tool is Monte-Carlo estimation, introduced by M. F. Hutchinson (Hutchinson 1989): average $M$ quadratic forms $z_i^TAz_i$ over random sign vectors $z_i$. Each sample costs one matrix–vector product, uses one random bit per entry, and needs only additions and subtractions.

Before Avron and Toledo (2011), only the variance of such estimators had been analysed. A small variance does not say how many samples guarantee a given relative error with a given probability. Avron and Toledo gave the first bounds of this kind for several estimators. This mission covers the bound for Hutchinson's estimator, their Theorem 7.1.

Setting

A Rademacher random variable takes the values $+1$ and $-1$, each with probability $1/2$. Let $A \in \mathbb{R}^{n\times n}$ be symmetric positive semi-definite. Draw $M \ge 1$ independent random vectors $z_1, \ldots, z_M \in \mathbb{R}^n$ whose $Mn$ entries are independent Rademacher variables. Hutchinson's trace estimator is

$$H_M = \frac{1}{M}\sum_{i=1}^{M} z_i^TAz_i .$$

A single sample $z^TAz$ is an unbiased estimator of $\mathrm{trace}(A)$ (Lemma 2.1 of the paper, due to Hutchinson). For symmetric $A$ its variance is $2\bigl(\|A\|_F^2 - \sum_i A_{ii}^2\bigr)$, twice the squared Frobenius mass of $A$ off the diagonal.

Given $\epsilon > 0$ and $\delta \in (0,1)$, a random estimator $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 .$$

The rank $\mathrm{rank}(A)$ is the number of nonzero eigenvalues $\lambda_j$ of $A$, counted with multiplicity.

Formalization targets

Goal: Theorem 7.1, sample bound for $H_M$

For every symmetric positive semi-definite $A$, every $0 < \epsilon \le 1/2$ and every $0 < \delta < 1$,

$$M \ge 6\epsilon^{-2}\ln\!\left(\frac{2\,\mathrm{rank}(A)}{\delta}\right) \quad\Longrightarrow\quad H_M \text{ is an } (\epsilon,\delta)\text{-approximator of } \mathrm{trace}(A).$$

The bound depends on $A$ only through its rank. It does not depend on the dimension $n$, on the condition number, or on how the trace is spread over the diagonal.

Milestones

1. Lemma 7.2 (Achlioptas 2001, Lemma 5). For a unit vector $\alpha \in \mathbb{R}^n$ and $S = \frac1M\sum_{i=1}^M(\alpha^Tz_i)^2$, for every $\epsilon > 0$,

$$\Pr(|S - 1| \ge \epsilon) \le 2\exp\!\left(-\frac{M}{2}\left(\frac{\epsilon^2}{2} - \frac{\epsilon^3}{3}\right)\right).$$

2. Per-direction bound (proof of Theorem 7.1, p. 8:11). Let $r \ge 1$ and $0 < \epsilon \le 1/2$. If $M \ge 6\epsilon^{-2}\ln(2r/\delta)$, then $\Pr(|S - 1| \ge \epsilon) \le \delta/r$.
3. From directions to the trace (proof of Theorem 7.1, p. 8:11). Write $A = U\Lambda U^T$ and $y_i = U^Tz_i$. If $|\frac1M\sum_i y_{ij}^2 - 1| \le \epsilon$ for every $j$ with $\lambda_j \ne 0$, then $|H_M - \mathrm{trace}(A)| \le \epsilon\,\mathrm{trace}(A)$. This step is deterministic.
4. Lemma 2.1 (Hutchinson). $\mathrm{E}(z^TAz) = \mathrm{trace}(A)$, and for symmetric $A$, $\mathrm{Var}(z^TAz) = 2(\|A\|_F^2 - \sum_i A_{ii}^2)$.

Significance

Theorem 7.1 gives a practitioner an explicit number of matrix–vector products after which Hutchinson's method is guaranteed to reach relative accuracy $\epsilon$ with confidence $1-\delta$. No bound was available before, although the method had been in wide use for two decades. The bound exceeds the one the same paper proves for Gaussian test vectors ($20\epsilon^{-2}\ln(2/\delta)$) by a $\ln(\mathrm{rank}(A))$ factor. The authors conjecture that this factor is not needed. Later work removed it: Roosta-Khorasani and Ascher (2015) proved a rank-free bound for the Rademacher case, and Cortinovis and Kressner (2022) extended sample bounds to indefinite matrices. Sample bounds of this type underlie variance-reduced estimators such as Hutch++ (Meyer et al. 2021).

The theorem is proved in the literature. To the platform's knowledge it has no machine-checked proof. Formalizing it produces a checked Rademacher concentration inequality for averages of squared linear forms (Achlioptas' lemma, which is also a core lemma of database-friendly Johnson–Lindenstrauss projections), and a checked spectral reduction from a quadratic-form estimator to its eigen-directions. Neither is currently in Mathlib.

Difficulty

The estimator is not an average of independent copies of a bounded variable with a small range: a single sample $z^TAz$ can deviate from $\mathrm{trace}(A)$ by an amount comparable to $\|A\|_F\sqrt n$. Chebyshev's inequality with the variance of Lemma 2.1 only gives a polynomial dependence on $1/\delta$. Hoeffding's inequality applied to $z^TAz$ directly gives a dependence on $n$ and on the size of the entries of $A$. Unlike the Gaussian case, the estimator cannot be written as a weighted sum of independent chi-squared variables, because rotating a Rademacher vector does not give another Rademacher vector. The central difficulty is Lemma 7.2: a tail bound for $(\alpha^Tz)^2$ that holds uniformly over every unit direction $\alpha$, including directions in which $\alpha^Tz$ is far from Gaussian (for $\alpha = e_1$ it is a constant).

Formalization scope

Matrices are Matrix (Fin n) (Fin n) ℝ, and positive semi-definiteness is Matrix.PosSemidef. The Rademacher law is $\tfrac12(\delta_1 + \delta_{-1})$ on ℝ. The sample space is Fin M → Fin n → ℝ with the product of $Mn$ copies of this law (hutchinsonSampleMeasure), so the law of the estimator is constructed, not assumed. $H_M(\omega) = (M:\mathbb{R})^{-1}\sum_i \omega_i\cdot(A\omega_i)$. Probabilities are Measure.real, and the $(\epsilon,\delta)$-approximator is Definition 4.1 verbatim. $\mathrm{rank}(A)$ is Matrix.rank. In Lemma 7.2 the i.i.d. copies $Q_i$ are the functions $(\alpha\cdot z_i)^2$ on this product space.

Corrections of printed statements, each recorded in the item's Formalization Note:

• Theorem 7.1 is printed without a range on $\epsilon$. Its proof needs $\frac{M}{2}(\frac{\epsilon^2}{2} - \frac{\epsilon^3}{3}) \ge \frac{M\epsilon^2}{6}$, which holds exactly when $\epsilon \le 1/2$. Without a range the printed statement is false: take $A = \frac1n\mathbf 1\mathbf 1^T$, $n = 10^4$, $\epsilon = 10$ and $\delta = 10^{-4}$. The condition admits $M = 1$, while $\Pr(H_1 > 11) \approx 9\cdot10^{-4} > \delta$. The goal and milestone 2 are therefore stated for $0 < \epsilon \le 1/2$.
• Lemma 2.1 is printed for an arbitrary $n\times n$ matrix. The variance formula fails for non-symmetric $A$: for $A = \begin{pmatrix}0&1\\0&0\end{pmatrix}$ the variance is $1$, not $2$. Symmetry is assumed for the variance part only.
• Proof of Theorem 7.1. The proof writes $\Lambda = UAU^T$ together with $y_i = U^Tz_i$. These fit together only for $A = U\Lambda U^T$, which is the convention of milestone 3.

For $A = 0$ the threshold involves $\ln 0$. Lean's Real.log 0 = 0 turns the condition into $M \ge 0$, which agrees with the paper's reading $\ln 0 = -\infty$. The conclusion then holds because $H_M = \mathrm{trace}(A) = 0$, so no extra hypothesis is added. A formalization that took the law of the samples as a hypothesis could make that hypothesis unsatisfiable and the theorem vacuous; the constructed product space rules this out.

A complete development needs:

• the product Rademacher measure and moment generating functions of Rademacher sums;
• a Chernoff bound for averages of i.i.d. bounded variables on a product space;
• the spectral theorem for real symmetric matrices, with rank equal to the number of nonzero eigenvalues;
• a finite union bound.

The Rademacher concentration results are reusable beyond this mission. Proofs of any milestone are welcome, as are alternative proofs of Lemma 7.2.

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 18(3), 1059–1076, 1989. https://doi.org/10.1080/03610918908812806
• D. Achlioptas, Database-friendly random projections, Proceedings of PODS 2001, 274–281. https://doi.org/10.1145/375551.375608
• F. Roosta-Khorasani and U. Ascher, Improved bounds on sample size for implicit matrix trace estimators, Foundations of Computational Mathematics 15, 1187–1212, 2015. https://doi.org/10.1007/s10208-014-9220-1
• A. Cortinovis and D. Kressner, On randomized trace estimates for indefinite matrices with an application to determinants, Foundations of Computational Mathematics 22, 875–903, 2022. https://doi.org/10.1007/s10208-021-09525-9
• R. A. Meyer, C. Musco, C. Musco and D. P. Woodruff, Hutch++: Optimal stochastic trace estimation, SOSA 2021, 142–155. https://doi.org/10.1137/1.9781611976496.16

Motivation

Many matrices in scientific computing are available only implicitly: one can multiply a vector by $A$, for instance by solving a linear system or applying a matrix function, but the entries of $A$ are never formed. Quantities such as the trace must then be estimated from a small number of matrix–vector products. Randomized trace estimators do exactly this: draw random vectors $z$, compute the quadratic forms $z^TAz$, and average.

Avron and Toledo (J. ACM 58(2), 2011) gave the first sample bounds of the form "$M$ samples suffice for relative error $\epsilon$ with probability $1-\delta$" for the standard estimators. Among their results is a case in which the estimate is not merely approximate but exact: when $A$ is a projection matrix, its trace equals its rank, an integer, and rounding a Gaussian trace estimate recovers that integer with high probability. Computing the rank of a projection arises, for example, when charge densities are computed in electronic structure calculations without diagonalization (Bekas, Kokiopoulou and Saad 2007).

Setting

Let $n \ge 1$ and $A \in \mathbb{R}^{n\times n}$. A projection matrix here is an orthogonal projection: $A$ is symmetric and $A^2 = A$. Equivalently, $A$ is symmetric and every eigenvalue of $A$ is $0$ or $1$; its rank $\mathrm{rank}(A)$ is the number of eigenvalues equal to $1$.

Fix a number of samples $M \ge 1$. Let $z_1,\ldots,z_M \in \mathbb{R}^n$ be random vectors whose $Mn$ entries are independent standard normal random variables. The Gaussian trace estimator (Definition 3.1 of the paper) is

$$G_M = \frac{1}{M}\sum_{i=1}^{M} z_i^T A z_i .$$

Its expectation is $\mathrm{trace}(A)$. For $x \in \mathbb{R}$, $\mathrm{round}(x)$ denotes the nearest integer to $x$. For $k \ge 1$, a random variable $X$ has the $\chi^2$ distribution with $k$ degrees of freedom, $X \sim \chi^2(k)$, if it has the law of $g_1^2 + \cdots + g_k^2$ for independent standard normal $g_1, \ldots, g_k$.

Formalization targets

Goal: Lemma 5.3 (p. 8:9)

For a projection matrix $A$, a failure probability $\delta > 0$, and every integer $M \ge 1$ with $M \ge 24\,\mathrm{rank}(A)\ln(2/\delta)$,

$$\Pr\bigl(\mathrm{round}(G_M) \ne \mathrm{rank}(A)\bigr) \le \delta .$$

The number of samples depends on the rank and on $\delta$ only; there is no accuracy parameter, and none on $n$.

Milestones, in the order the proof uses them

1. Law of $MG_M$. $MG_M \sim \chi^2(M\,\mathrm{rank}(A))$.
2. $\chi^2$ tail bound (cited by the paper from Li, Hastie and Church 2007). For $X \sim \chi^2(k)$, $k\ge1$ and $0 < \epsilon \le \tfrac12$,

$$\Pr(|X - k| \ge \epsilon k) \le 2\exp(-k\epsilon^2/6).$$

3. Tail of $G_M$. For $\mathrm{rank}(A) \ge 1$ and $0 < \epsilon \le \tfrac12$,

$$\Pr(|G_M - \mathrm{rank}(A)| \ge \mathrm{rank}(A)\epsilon) \le 2\exp(-M\,\mathrm{rank}(A)\epsilon^2/6).$$

4. Eq. (2). If moreover $M \ge 6\,\mathrm{rank}(A)^{-1}\epsilon^{-2}\ln(2/\delta)$, then $\Pr(|G_M - \mathrm{rank}(A)| \ge \mathrm{rank}(A)\epsilon) \le \delta$.
5. Trace of a projection. $\mathrm{trace}(A) = \mathrm{rank}(A)$.

Significance

The lemma turns a randomized estimator into an exact algorithm with a controlled failure probability: the rank of an implicitly given projection is obtained from $O(\mathrm{rank}(A)\log(1/\delta))$ matrix–vector products, independently of the dimension $n$. It is also an instance where the paper's general relative-error bound for the Gaussian estimator (Theorem 5.2, whose sample count grows like $\epsilon^{-2}$) is improved by exploiting the spectrum of $A$: an absolute error below $\tfrac12$ is a relative error $\epsilon = 1/(2\,\mathrm{rank}(A))$, for which the general bound would require a number of samples quadratic in the rank, whereas Lemma 5.3 needs only a linear number.

The result is proved in the paper. The mission produces a machine-checked version of it, together with two pieces of reusable substrate: the exact $\chi^2$ law of a Gaussian quadratic form in an orthogonal projection, and a two-sided $\chi^2$ tail bound with explicit constant $1/6$ on the range $0<\epsilon\le 1/2$. To the best of available knowledge, none of these statements is formalized in Mathlib; a platform mission on the Johnson–Lindenstrauss lemma states a $\chi^2$ concentration bound with a different exponent, $(\epsilon^2-\epsilon^3)/4$, on the open range $0<\epsilon<1/2$, which does not cover the value $\epsilon = 1/2$ needed here when $\mathrm{rank}(A)=1$.

Difficulty

The deterministic part is short; the probabilistic part is not. The step "$y = Uz$ has independent standard normal entries because $U$ is orthogonal" is the rotation invariance of the standard Gaussian measure on $\mathbb{R}^n$, and it must be combined with the independence of the $M$ samples to identify the law of a sum of $M\,\mathrm{rank}(A)$ squares; this is a statement about product measures and pushforwards, not about moments. The $\chi^2$ tail bound is quoted by the paper without proof. Its constant $1/6$ is not the constant of the usual textbook $\chi^2$ estimates, and the bound is false outside a restricted range of $\epsilon$ (see below), so a generic sub-exponential concentration inequality with unspecified constants does not deliver it. Finally, the rounding step requires matching Mathlib's round with the event $|G_M - \mathrm{rank}(A)| < \tfrac12$.

Formalization scope

All declarations live in the namespace TraceEstimation.ProjectionRank. Matrices are Matrix (Fin n) (Fin n) ℝ. The sample space of $G_M$ is Fin M → Fin n → ℝ with the product measure Measure.pi (fun _ => Measure.pi (fun _ => gaussianReal 0 1)), so the law of the estimator is constructed, not assumed. Probabilities are Measure.real of events. "Projection matrix" is A.IsHermitian ∧ A * A = A (orthogonal projection); a non-symmetric idempotent also has eigenvalues $0$ and $1$, but the paper's proof diagonalizes $A$ by a unitary matrix, which requires symmetry. round is Mathlib's round : ℝ → ℤ; the rank is Matrix.rank. The $\chi^2$ distribution is not defined: its role is played by the pushforward of a product of standard normals under $g \mapsto \sum_l g_l^2$.

Corrections of the printed text, both in the $\chi^2$ tail bound (milestone 2):

• The paper prints $\Pr(|X - k| \le \epsilon k) \le 2\exp(-k\epsilon^2/6)$. The inner $\le$ is a misprint for $\ge$; the next display applies the bound with $\ge$.
• The paper gives no range for $\epsilon$. The bound fails for $\epsilon = 1$ and large $k$, since $\Pr(X \ge 2k)$ decays like $e^{-k(1-\ln 2)/2}$, slower than $e^{-k/6}$. The mission states it for $0<\epsilon\le 1/2$, and milestones 3 and 4 inherit that range; the proof uses only $\epsilon = 1/(2\,\mathrm{rank}(A)) \le 1/2$.

The goal keeps the paper's hypotheses: $\delta > 0$ with no upper bound, and $\mathrm{rank}(A) = 0$ allowed (both are true cases). A statement about the event $|G_M - \mathrm{rank}(A)| < \tfrac12$, or Eq. (2) alone, is not the goal; the goal is about $\mathrm{round}(G_M)$. The sample measure is a probability measure, so the bound $\le \delta$ is not vacuous.

Contributions welcome beyond the milestones: rotation invariance of the standard Gaussian vector under orthogonal matrices, the moment generating function of $\chi^2(k)$, and the spectral fact $\mathrm{trace}(A) = \mathrm{rank}(A)$ for idempotents, each reusable well outside this mission.

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
• P. Li, T. Hastie and K. Church, Nonlinear estimators and tail bounds for dimension reduction in $l_1$ using Cauchy random projections, in Learning Theory (COLT 2007), Lecture Notes in Computer Science 4539, Springer, 514–529, 2007 (the version the paper cites); journal version in Journal of Machine Learning Research 8, 2497–2532, 2007, https://jmlr.org/papers/v8/li07b.html
• C. Bekas, E. Kokiopoulou and Y. Saad, An estimator for the diagonal of a matrix, Applied Numerical Mathematics 57(11–12), 1214–1229, 2007. https://doi.org/10.1016/j.apnum.2007.01.003
• 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