Linear algebra

Motivation

The least-squares (LS) problem $\min_x \|Ax - b\|$ assumes that the data $A \in \mathbb{R}^{n\times m}$, $b \in \mathbb{R}^n$ are exact. In applications they rarely are: they come from measurements, from linearizations, or from models with neglected dynamics. A classical response is sensitivity analysis or regularization (Tikhonov), where a weight trades the size of the solution against the fit, and the choice of that weight is left to the user. El Ghaoui and Lebret (SIAM J. Matrix Anal. Appl. 18(4), 1997) take a deterministic view instead: the true data lie in a known ball around $(A, b)$, and the solution should minimize the residual it can be forced to have in the worst case over that ball. The paper shows that this robust least-squares (RLS) problem is solvable exactly, in the unstructured case by a second-order cone program (SOCP). The same worst-case idea, applied to regression, underlies the later equivalence between robustness and regularization (Xu, Caramanis and Mannor, 2009) and is a standard entry point to robust optimization (Ben-Tal, El Ghaoui and Nemirovski, Robust Optimization, 2009).

This mission formalizes the first main result of the paper, Theorem 3.1: the worst-case residual has a closed form, its minimizer is unique, and minimizing it is an SOCP.

Setting

Vectors carry the Euclidean norm $\|v\| = (\sum_i v_i^2)^{1/2}$. For a matrix $X$, $\|X\|_F = (\sum_{i,j} X_{ij}^2)^{1/2}$ is the Frobenius norm and $\|X\|$ the largest singular value, i.e. the smallest $c \ge 0$ with $\|Xv\| \le c\|v\|$ for all $v$.

Fix $A \in \mathbb{R}^{n\times m}$ and $b \in \mathbb{R}^n$. A perturbation is a pair $\Delta A \in \mathbb{R}^{n\times m}$, $\Delta b \in \mathbb{R}^n$, collected in the augmented matrix $\Delta = [\Delta A\ \Delta b] \in \mathbb{R}^{n\times(m+1)}$. For a bound $\rho \ge 0$ and $x \in \mathbb{R}^m$, the worst-case residual is (paper, eq. (1))

$$r(A,b,\rho,x) = \max_{\|[\Delta A\ \Delta b]\|_F \le \rho} \|(A+\Delta A)x - (b+\Delta b)\|,$$

and $x$ is an RLS solution if it minimizes $r(A,b,\rho,\cdot)$. The bound constrains the augmented matrix jointly, not $\Delta A$ and $\Delta b$ separately. The paper normalizes $\rho = 1$ and writes $r(A,b,x) = r(A,b,1,x)$. Finally, $[x;1] \in \mathbb{R}^{m+1}$ denotes $x$ stacked over $1$. In the Lean development these are RobustLS.Unstructured.eucNorm, frobNorm, specNorm, augment, stackOne, worstCaseResidual A b ρ x, its largest-singular-value variant worstCaseResidualSpec, and the SOCP constraint predicate SocpFeasible A b x λ τ.

Formalization targets

Goal: Theorem 3.1 (p. 1040)

For $n \ge 1$, every $A$, $b$:

$$r(A,b,x) = \|Ax-b\| + \sqrt{\|x\|^2+1} \quad \text{for all } x \in \mathbb{R}^m,$$

the problem $\min_{x \in \mathbb{R}^m} r(A,b,x)$ has exactly one solution $x_{\mathrm{RLS}}$, and it is the SOCP

$$\text{minimize } \lambda \quad\text{subject to}\quad \|Ax-b\| \le \lambda-\tau,\quad \|[x;1]\| \le \tau, \tag{15}$$

in the sense that $r(A,b,x)$ is the least $\lambda$ for which some $\tau$ makes $(x,\lambda,\tau)$ feasible.

Milestones

1. Eq. (16). Every perturbation with $\|[\Delta A\ \Delta b]\|_F \le 1$ has residual at most $\|Ax-b\| + \sqrt{\|x\|^2+1}$.
2. The worst-case perturbation. For a unit vector $u$ aligned with $Ax - b$ (arbitrary if $Ax = b$), the rank-one matrix $\Delta = u[x^T\ {-1}]/\sqrt{\|x\|^2+1}$ has $\|\Delta\|_F = \|\Delta\| = 1$ and attains the bound.
3. Spectral norm. The worst case over the larger ball $\|[\Delta A\ \Delta b]\| \le 1$ is the same value.
4. Strict convexity. $x \mapsto r(A,b,x)$ is strictly convex on $\mathbb{R}^m$.
5. The SOCP (15). For every $x$, $r(A,b,x)$ is the optimal $\lambda$ of (15) with $x$ fixed, and $x$ is an RLS solution exactly when it is the $x$-part of an optimal solution of (15).

Significance

The closed form replaces a maximization over a matrix ball of dimension $n(m+1)$ by two Euclidean norms. It shows that the RLS objective is the LS residual plus a penalty $\sqrt{\|x\|^2+1}$ that does not depend on $A$ or $b$, which is the starting point for the paper's Theorem 3.2 (the RLS solution is a Tikhonov-regularized LS solution with a data-dependent weight) and its analysis of continuity and conditioning. The SOCP formulation places the problem in the class solved by interior-point methods, at a cost the paper compares with one singular value decomposition of $A$. The spectral-norm statement says the worst case does not depend on which of the two standard matrix norms bounds the perturbation.

The result is proved in the paper; the proof is short. To the best of the planning survey (September 2026), no machine-checked proof exists, and Prove2Me has no statement about worst-case residuals or robust least squares. The mission produces a verified closed form that later missions of this series (Tikhonov form of the solution, structured and linear-fractional perturbations) and any formalization of robust regression can import.

Difficulty

The upper bound alone does not give the theorem: the statement is an equality, and the equality needs an explicit maximizer. The paper's printed maximizer is wrong by a sign: with $[x^T\ 1]$ in place of $[x^T\ {-1}]$ the perturbation does not attain the bound (for $A = 0$, $x = 0$, $b = e_1$ it gives residual $0$ instead of $2$), so a transcription of the printed proof fails. Two further points are silent in the paper. The operator norm of a rank-one matrix has to be computed from the definition of the largest singular value. Uniqueness of the minimizer needs existence first, which follows from growth of $r$ at infinity and is not stated. Working with the sSup definition of the worst case requires showing the set of residuals is bounded, which is milestone 1.

Formalization scope

• Dimensions are Fin n, Fin m; $A$ is Matrix (Fin n) (Fin m) ℝ, $b$ and $x$ are functions Fin n → ℝ, Fin m → ℝ. The augmented matrix $[\Delta A\ \Delta b]$ is indexed by Fin m ⊕ Unit, and so is $[x;1]$.
• Vector norms are the Euclidean norm written as $\sqrt{\sum_i v_i^2}$ (eucNorm), never Mathlib's ‖·‖ on Fin n → ℝ, which is the sup norm. The Frobenius norm and the largest singular value are explicit definitions (frobNorm, specNorm); specNorm is the infimum of admissible operator constants.
• The maximum in (1) is sSup of the set of attained residuals. For $\rho \ge 0$ the set is nonempty and bounded, so this is the true maximum; milestones 1 and 2 state the bound and the attaining perturbation directly, so no statement relies on the value of sSup on an unbounded set.
• The paper's normalization $\rho = 1$ is kept; general $\rho > 0$ follows from the scaling $\phi(A,b,\rho) = \rho\,\phi(A/\rho,b/\rho,1)$ the paper records on p. 1039 and is not a target.
• The goal assumes $n \ge 1$. For $n = 0$ the only perturbation is the empty matrix, the worst case is $0$, and the closed form fails; the paper's setting ($Ax \simeq b$ with data $b \in \mathbb{R}^n$) has $n \ge 1$. Milestones 3–5 carry the same hypothesis.
• Milestone 2 states the corrected perturbation $[x^T\ {-1}]$; the printed $[x^T\ 1]$ is false.
• A trivializing formalization — an upper bound in place of the equality, a worst case over $\Delta A$ and $\Delta b$ bounded separately, or uniqueness among critical points only — is ruled out: the goal is the equality for the jointly bounded augmented matrix and ∃! of a global minimizer over all of $\mathbb{R}^m$.

Contributions welcome: lemmas on Frobenius and operator norms of rank-one matrices, the inequality $\|Mz\| \le \|M\|_F\|z\|$ in this explicit setting, and strict convexity of $x \mapsto \sqrt{\|x\|^2+1}$; these are reusable beyond the mission.

Selected references

• L. El Ghaoui and H. Lebret, Robust Solutions to Least-Squares Problems with Uncertain Data, SIAM Journal on Matrix Analysis and Applications 18(4):1035–1064, 1997. https://doi.org/10.1137/S0895479896298130
• A. Ben-Tal, L. El Ghaoui and A. Nemirovski, Robust Optimization, Princeton University Press, 2009. https://doi.org/10.1515/9781400831050
• H. Xu, C. Caramanis and S. Mannor, Robust Regression and Lasso, Journal of Machine Learning Research 10:1485–1510, 2009 (IEEE Trans. Inf. Theory 56(7), 2010). https://jmlr.org/papers/v10/xu09b.html
• M. S. Lobo, L. Vandenberghe, S. Boyd and H. Lebret, Applications of Second-Order Cone Programming, Linear Algebra and its Applications 284:193–228, 1998. https://doi.org/10.1016/S0024-3795(98)10032-0

Motivation

Many computations in scientific computing, statistics and machine learning need the trace of a matrix $A$ that is never formed explicitly: $A$ may be $f(B)$ for a large sparse $B$, an inverse $B^{-1}$, or a product of operators, and the only access to it is the ability to compute products $Az$ for chosen vectors $z$. Examples are log-determinant estimation in Gaussian process regression, counting triangles in graphs through $\mathrm{trace}(B^3)$, computing charge densities in electronic structure calculations, and generalized cross-validation in regularized regression. For such matrices the $n$ diagonal entries are not available, and computing them one at a time costs $n$ matrix–vector products.

Randomized trace estimators replace this by a small number $M$ of products: draw random vectors $z_1,\ldots,z_M$ from a fixed distribution with $\mathrm{E}(z_i^T A z_i) = \mathrm{trace}(A)$ and average the quadratic forms. Hutchinson (1989) introduced the estimator with Rademacher vectors and computed its variance; Silver and Röder (1997) used Gaussian vectors. Before Avron and Toledo (2011), the analyses of these estimators were variance computations, which do not say how many samples guarantee a given relative accuracy with a given probability. Avron and Toledo gave the first such sample bounds for several estimators, stated in terms of an $(\epsilon,\delta)$ guarantee; this mission formalizes their bound for the Gaussian estimator. Later work (Roosta-Khorasani and Ascher 2015; Cortinovis and Kressner 2022) sharpened these bounds and extended them to indefinite matrices.

Setting

Let $A \in \mathbb{R}^{n\times n}$ be symmetric positive semi-definite, and write $\tau = \mathrm{trace}(A)$. 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) is

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

Each term $z_i^TAz_i$ has expectation $\mathrm{trace}(A)$, so $G_M$ is unbiased. A randomized trace estimator $T$ is an $(\epsilon,\delta)$-approximator of $\mathrm{trace}(A)$ (Definition 4.1) if

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

that is, if its relative error is at most $\epsilon$ except on an event of probability at most $\delta$.

The analysis uses the eigenvalues $\lambda_1,\ldots,\lambda_n \ge 0$ of $A$, listed with multiplicity, and the polynomial $h(t) = \sum_{s=2}^{n}(-2)^s t^s \sum_{|S| = s}\prod_{i\in S}\lambda_i$, where $S$ ranges over subsets of $\{1,\ldots,n\}$; it satisfies $\prod_i(1-2\lambda_i t) = 1 - 2\tau t + h(t)$.

Formalization targets

Goal: Theorem 5.2 (corrected)

For every symmetric positive semi-definite $A$, every $0 < \epsilon \le 1/10$, every $0<\delta<1$ and every natural number $M$ with

$$M \ge 20\,\epsilon^{-2}\ln(2/\delta),$$

the estimator $G_M$ is an $(\epsilon,\delta)$-approximator of $\mathrm{trace}(A)$. The sample count depends on neither $n$ nor $A$.

Milestones

1. Lemma 5.1. For symmetric $A$, $\mathrm{E}(G_1) = \mathrm{trace}(A)$ and $\mathrm{Var}(G_1) = 2\|A\|_F^2$.
2. Eq. (1). For symmetric $A$ and every $t$ with $2\lambda_i t < 1$ for all $i$, the moment generating function of $Z = MG_M$ is

$$m_Z(t) = \prod_{i=1}^{n}(1-2\lambda_i t)^{-M/2} = (1 - 2\tau t + h(t))^{-M/2}.$$

3. Elementary symmetric sums (p. 8:8). For non-negative $x_1,\ldots,x_n$ and $1\le i\le n$, $\sum_{|S|=i}\prod_{j\in S}x_j \le (\sum_j x_j)^i$; hence $|h(t)| \le \sum_{j=2}^{n}(2\tau t)^j$ for $t\ge0$ when all $\lambda_i\ge0$.
4. Upper tail (pp. 8:8–8:9). If $\tau > 0$, $M\ge1$ and $0<\epsilon\le 0.1$, then

$$\Pr\bigl(G_M \ge \tau(1+\epsilon)\bigr) \le \exp(-M\epsilon^2/20).$$

5. Both tails (p. 8:9). If $\tau>0$, $0<\epsilon\le0.1$ and $M \ge 20\epsilon^{-2}\ln(2/\delta)$, then $\Pr(G_M \ge \tau(1+\epsilon)) \le \delta/2$ and $\Pr(G_M \le \tau(1-\epsilon)) \le \delta/2$.

Significance

The theorem gives a number of matrix–vector products, $O(\epsilon^{-2}\ln(1/\delta))$, that suffices for a relative-error guarantee on the trace of any positive semi-definite matrix, independent of its dimension and spectrum. It is the reference row of the paper's Table I, against which the Hutchinson, normalized Rayleigh-quotient and unit-vector estimators are compared, and it is the form in which trace estimation enters the analysis of randomized algorithms for log-determinants, spectral densities and matrix functions.

The result is proved in the paper; to our knowledge it has not been formalized in any proof assistant. The mission produces a machine-checked version with the constant $20$ and the range of $\epsilon$ made explicit, and with the misprints of the printed argument resolved (see Formalization scope). It also produces reusable pieces: the moment generating function of a Gaussian quadratic form, and the bound on elementary symmetric sums by powers of the power sum. Sharper constants, the removal of the restriction $\epsilon\le 0.1$, or a direct formalization of the lower tail through a $\chi^2$ tail bound are welcome as further theorems.

Difficulty

Unbiasedness and the variance formula do not give the result: Chebyshev's inequality with $\mathrm{Var}(G_M) = 2\|A\|_F^2/M \le 2\tau^2/M$ yields $M \ge 2\epsilon^{-2}\delta^{-1}$, with a polynomial rather than logarithmic dependence on $1/\delta$. A logarithmic bound needs exponential moments of $G_M$, and the exponential moment of $z^TAz$ is finite only for $t$ below $1/(2\lambda_{\max})$; the argument has to choose $t$ inside that range uniformly in the spectrum, using only $\lambda_{\max}\le\tau$. The second obstacle is distributional: $z^TAz$ is not a sum of independent terms in the coordinates of $z$, and reducing it to a weighted sum of independent $\chi^2$ variables requires the rotation invariance of the standard Gaussian vector. The paper proves only the upper tail with an explicit constant and states that the lower tail follows "using the same technique"; that step has to be supplied.

Formalization scope

Matrices are Matrix (Fin n) (Fin n) ℝ; "symmetric positive semi-definite" is A.PosSemidef, "symmetric" is A.IsHermitian, and eigenvalues are Matrix.IsHermitian.eigenvalues. The sample space is Fin M → Fin n → ℝ with the product measure of $Mn$ copies of gaussianReal 0 1, so the law of the samples is constructed, not assumed; $G_M(\omega) = (M:\mathbb{R})^{-1}\sum_i \omega_i\cdot(A\omega_i)$. Probabilities are Measure.real; the moment generating function is Mathlib's mgf; the variance is Mathlib's variance, and Lemma 5.1 asserts square integrability so that neither the integral nor the variance takes its default value. The sample count is a natural number $M\ge1$; $\epsilon$ and $\delta$ are real.

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

• Theorem 5.2 is printed without a range for $\epsilon$ and is false without one (for rank-one $A$, $\epsilon=100$, $\delta=e^{-400}$ the threshold allows $M=1$, while $\Pr(\chi^2_1>101)\approx e^{-50}>\delta$). The proof gives its key bound "for $\epsilon\le0.1$"; the goal is stated for $0<\epsilon\le 1/10$.
• Eq. (1) is printed for $|\lambda_i t|\le\frac12$, which admits $1-2\lambda_it = 0$, where the moment generating function is infinite. It is stated for $2\lambda_i t<1$. The page's sum over subsets of "the set $\Lambda$ of eigenvalues" is taken over index sets, so repeated eigenvalues count with multiplicity.
• The last paragraph of the proof prints $\Pr(G_M\le\tau(1+\epsilon))\le\delta/2$ for the upper tail and $\Pr(|G_M-\tau|\le\tau(1+\epsilon))\le\delta$ for the conclusion; the intended statements are $\Pr(G_M\ge\tau(1+\epsilon))\le\delta/2$ and $\Pr(|G_M-\tau|>\epsilon\tau)\le\delta$. Milestone 5 states the upper and lower tail bounds.
• Lemma 5.1 is followed by the remark that it "also applies when $A$ is non-symmetric"; this is false for the variance and is not formalized. Definition 3.1 says "positive-definite"; the estimator is defined for every matrix and each theorem carries its own hypothesis.
• The one-sided tail milestones assume $\mathrm{trace}(A)>0$; for $A=0$ their events are certain, while the goal holds trivially.

A formalization that makes the goal trivial is ruled out: the estimator's law is the explicit product Gaussian measure rather than a hypothesis, $M$ ranges over all natural numbers above the threshold, and the approximator predicate is evaluated on the genuine event $|G_M - \mathrm{trace}(A)|\le\epsilon\,\mathrm{trace}(A)$.

A complete development needs rotation invariance of the standard Gaussian on $\mathbb{R}^n$ (Mathlib's stdGaussian_map), the moment generating function of a squared standard normal, independence of products of Gaussian vectors, and a Chernoff bound from the moment generating function. The Gaussian quadratic-form results (milestones 1 and 2) are reusable for the other Gaussian estimators of the paper, such as the rank estimator of Lemma 5.3. Proofs of any milestone, alternative proofs of the lower tail, and helper lemmas on $\chi^2$ moment generating functions are welcome.

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
• R. N. Silver and H. Röder, Calculation of densities of states and spectral functions by Chebyshev recursion and maximum entropy, Physical Review E 56(4), 4822–4829, 1997. https://doi.org/10.1103/PhysRevE.56.4822
• 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

Motivation

In high-dimensional linear regression one observes $y = X\beta^* + w \in \mathbb R^n$ with a design matrix $X\in\mathbb R^{n\times M}$ whose number of columns $M$ may far exceed the sample size $n$. The two standard estimators of a sparse $\beta^*$, the Lasso (Tibshirani, 1996) and the Dantzig selector (Candès and Tao, 2007), both come with error bounds of order $s\log M/n$ for an $s$-sparse $\beta^*$, but only under a condition on $X$: since $X$ has a non-trivial kernel when $M>n$, some form of restricted invertibility is unavoidable.

Bickel, Ritov and Tsybakov (arXiv:0801.1095, Ann. Statist. 2009) introduced the restricted eigenvalue (RE) condition, which asks for invertibility of $X$ only on a cone of approximately sparse vectors. It is weaker than the conditions used before it and has since become the default assumption in the sparse-estimation literature. Section 4 of the paper relates RE to the earlier conditions:

• 2005–2007: Candès and Tao (arXiv:math/0506081) analyse the Dantzig selector under a uniform uncertainty principle involving restricted eigenvalues and restricted correlations of $X$; the condition $\phi_{\min}(2s)>\theta_{s,2s}$ is Assumption 1 below with $c_0=1$.
• 2006–2009: Meinshausen and Yu (arXiv:math/0605584) analyse the Lasso under a lower bound on sparse eigenvalues of order $s\log n$.
• 2006: Donoho, Elad and Temlyakov (doi:10.1109/TIT.2005.860430) use mutual coherence for sparse recovery; 2007: Bunea, Tsybakov and Wegkamp (doi:10.1214/07-EJS008) use coherence-type conditions for the Lasso.
• 2009: Bickel, Ritov and Tsybakov show (Lemma 4.1 and Section 4) that each of these conditions implies RE.

This mission formalizes those implications.

Setting

Fix integers $n\ge1$ and $M\ge2$ and a matrix $X\in\mathbb R^{n\times M}$ with columns $x_1,\dots,x_M$. The Gram matrix is $\Psi_n = X^TX/n$. For $\delta\in\mathbb R^M$ and $J\subseteq\{1,\dots,M\}$, $\delta_J$ is the vector equal to $\delta$ on $J$ and $0$ off $J$; $|\cdot|_1$, $|\cdot|_2$ are the $\ell_1$ and Euclidean norms; $\mathcal M(\delta)$ is the number of non-zero coordinates of $\delta$; $J_0^c$ is the complement of $J_0$.

The cone condition for $J_0$ and $c_0>0$ is

$$|\delta_{J_0^c}|_1\le c_0\,|\delta_{J_0}|_1. \tag{4.1}$$

Assumption RE$(s,c_0)$ holds with constant $\kappa>0$ if $|X\delta|_2\ge\kappa\sqrt n\,|\delta_{J_0}|_2$ for every $J_0$ with $|J_0|\le s$ and every $\delta\ne0$ satisfying (4.1). For $m\ge s$, let $J_1$ be a set of $m$ indices outside $J_0$ carrying the $m$ largest $|\delta_j|$, and $J_{01}=J_0\cup J_1$; Assumption RE$(s,m,c_0)$ replaces $|\delta_{J_0}|_2$ by $|\delta_{J_{01}}|_2$.

The restricted eigenvalues are $\phi_{\min}(u)$ and $\phi_{\max}(u)$, the minimum and maximum of $x^T\Psi_nx/|x|_2^2$ over $x$ with $1\le\mathcal M(x)\le u$. The restricted correlations $\theta_{m_1,m_2}$ are the maximum of $c_1^TX_{I_1}^TX_{I_2}c_2/(n|c_1|_2|c_2|_2)$ over disjoint index sets $I_1,I_2$ with $|I_i|\le m_i$ and non-zero $c_i\in\mathbb R^{I_i}$. Two constants are attached to them:

$$\kappa_1(s,c_0)=\sqrt{\phi_{\min}(2s)}\Big(1-\frac{c_0\theta_{s,2s}}{\phi_{\min}(2s)}\Big),\qquad \kappa_2(s,m,c_0)=\sqrt{\phi_{\min}(s+m)}\Big(1-c_0\sqrt{\tfrac{s\,\phi_{\max}(m)}{m\,\phi_{\min}(s+m)}}\Big).$$

$P_{01}$ is the orthogonal projector in $\mathbb R^n$ onto the span of the columns $x_j$, $j\in J_{01}$.

Formalization targets

Goal: Lemma 4.1 (ii)

For integers $1\le s\le M/2$, $m\ge s$, $s+m\le M$ and $c_0>0$, if Assumption 2 $m\,\phi_{\min}(s+m)>c_0^2\,s\,\phi_{\max}(m)$ holds, then $\kappa_2(s,m,c_0)>0$, RE$(s,c_0)$ and RE$(s,m,c_0)$ hold with constant $\kappa_2(s,m,c_0)$, and for every $J_0$ with $|J_0|\le s$ and every $\delta$ satisfying (4.1)

$$\frac1{\sqrt n}|P_{01}X\delta|_2\ \ge\ \kappa_2(s,m,c_0)\,|\delta_{J_{01}}|_2 .$$

Assumption 2 involves no correlations, only extreme eigenvalues of small principal submatrices of $\Psi_n$.

Lemma 4.1 (i)

For $1\le s\le M/2$ and $c_0>0$, Assumption 1 $\phi_{\min}(2s)>c_0\theta_{s,2s}$ implies the same conclusions with $m=s$ and constant $\kappa_1(s,c_0)$.

Coherence-type conditions (Section 4)

For $1\le s\le M$ and $c_0>0$, each of

$$\phi_{\min}(s)>2c_0\theta_{s,1}\sqrt s,\qquad \phi_{\min}(s)>2c_0\theta_{1,1}s,\qquad \operatorname{diag}\Psi_n=1\ \text{and}\ \theta_{1,1}<\frac1{(1+2c_0)s}$$

(Assumptions 3, 4, 5) implies RE$(s,c_0)$, with the constants $\kappa^2=\phi_{\min}(s)-2c_0\theta_{s,1}\sqrt s$, $\phi_{\min}(s)-2c_0\theta_{1,1}s$ and $1-(1+2c_0)\theta_{1,1}s$ respectively.

The milestones are the steps of the proof in Appendix A — the projection inequality (A.1), the block bound (A.2), the shelling bound (A.3), the Candès–Tao correlation bound used for part (i) — followed by part (i) and the three coherence-type implications.

Significance

RE$(s,c_0)$ with $c_0=3$ and $c_0=1$ is the hypothesis of the paper's prediction and $\ell_1$ bounds for the Lasso and the Dantzig selector (Theorems 5.1, 6.1, 7.1, 7.2), and RE$(s,m,c_0)$ is the hypothesis of its $\ell_p$ bounds. Assumptions 1–5 are stated through quantities that are standard in compressed sensing and random matrix theory, so known bounds for $\phi_{\min}$, $\phi_{\max}$ and $\theta$ of random designs transfer, through this mission's theorems, to every result stated under RE. Lemma 4.1 also shows that RE is weaker than the Candès–Tao condition used for the Dantzig selector.

The results are proved in the paper; parts of Lemma 4.1's proof (the correlation bound for part (i)) are cited from Candès and Tao without proof. None of these results is formalized: the platform has pairwise-incoherence and restricted-nullspace statements from Wainwright's textbook (a different conclusion and normalization) and restricted isometry definitions, but neither restricted eigenvalues $\phi_{\min}(u),\phi_{\max}(u)$, restricted correlations $\theta_{m_1,m_2}$, nor the RE condition in this form.

Difficulty

The naive attempt to bound $|X\delta|_2$ from below splits $\delta=\delta_{J_0}+\delta_{J_0^c}$ and applies an eigenvalue bound to each part. This fails: $\delta_{J_0^c}$ can have up to $M-s$ non-zero coordinates, and no condition on $s$- or $2s$-sparse submatrices controls $|X\delta_{J_0^c}|_2$ directly. The cone condition bounds only the $\ell_1$ norm of $\delta_{J_0^c}$, while eigenvalue conditions speak about $\ell_2$ norms of sparse vectors; bridging the two with the right constant $\sqrt{s/m}$, and keeping track of how the leading block $J_{01}$ interacts with the rest through the projector $P_{01}$, is where the work lies. For part (i), the interaction between disjoint sparse blocks has to be controlled by $\theta_{s,2s}$ rather than by $\phi_{\max}$.

Formalization scope

• Representation. $X$ is Matrix (Fin n) (Fin M) ℝ; vectors are Fin M → ℝ and Fin n → ℝ; $|X\delta|_2=(\sum_i (X\delta)_i^2)^{1/2}$. The projector $P_{01}$ is Mathlib's orthogonal projection on EuclideanSpace ℝ (Fin n) onto the span of the columns indexed by $J_{01}$.
• RE through a witness. RE X s c0 κ asserts the RE inequality with constant $\kappa$ for all admissible $J_0$ and $\delta$. The paper's $\kappa(s,c_0)$ is the largest such $\kappa$ (the minimum is attained), so "RE holds with $\kappa(s,c_0)\ge\kappa_2$" is exactly "$\kappa_2>0$ is a witness". This avoids a real infimum over an empty set when $J_0=\emptyset$.
• Ties. Every admissible choice of $J_1$ (the $m$ largest $|\delta_j|$ outside $J_0$) is quantified over.
• Restricted eigenvalues and correlations are sInf/sSup over nonempty bounded sets (a basis vector for $\phi$; two disjoint singletons for $\theta$, since $M\ge2$), so they equal the paper's attained min/max. $u$, $s$, $m$ are natural numbers; $s\le M/2$ is written $2s\le M$.
• Corrections of the printed statement. (1) Lemma 4.1 says the RE assumptions "hold with $\kappa(s,c_0)=\kappa(s,m,c_0)=\kappa_2(s,m,c_0)$" (and likewise with $\kappa_1$); the proof gives only the lower bound, and the lower bound is what is stated. (2) The paper calls $P_{01}$ "the projector in $\mathbb R^M$"; it acts on $\mathbb R^n$. (3) The Section 4 claims "Assumption 3/4/5 implies RE$(s,c_0)$" are stated with the explicit constant produced by the displayed argument, a labelled strengthening. (4) The Candès–Tao bound is stated with the hypotheses the proof uses: the blocks are disjoint, of sizes at most $s$ and $2s$, and $\phi_{\min}(2s)>0$.
• Ruling out trivializations. RE quantifies over all $J_0$ with $|J_0|\le s$ and all non-zero $\delta$ in the cone, and bounds the full $|X\delta|_2$, not $|X\delta_{J_0}|_2$; no hypothesis restricts $X$ beyond the stated assumptions. The hypotheses are satisfiable: for $n=M=4$, $X=2I$ (so $\Psi_n=I$), $s=1$, $m=2$, $c_0=1$, Assumption 2 reads $2>1$.
• Infrastructure. A sparse-vector library (restriction, support, sorting coordinates into blocks) and facts about orthogonal projections onto column spans are needed; both are reusable for the other missions of this series and for compressed-sensing results. Proofs of any milestone, and alternative arguments, are welcome.

Selected references

• P. J. Bickel, Y. Ritov, A. B. Tsybakov, Simultaneous analysis of Lasso and Dantzig selector, Ann. Statist. 37(4), 1705–1732, 2009. arXiv:0801.1095v3. https://arxiv.org/abs/0801.1095
• E. Candès, T. Tao, The Dantzig selector: statistical estimation when p is much larger than n, Ann. Statist. 35(6), 2313–2351, 2007. https://arxiv.org/abs/math/0506081
• N. Meinshausen, B. Yu, Lasso-type recovery of sparse representations for high-dimensional data, Ann. Statist. 37(1), 246–270, 2009. https://arxiv.org/abs/math/0605584
• F. Bunea, A. B. Tsybakov, M. H. Wegkamp, Sparsity oracle inequalities for the Lasso, Electron. J. Statist. 1, 169–194, 2007. https://doi.org/10.1214/07-EJS008
• D. L. Donoho, M. Elad, V. N. Temlyakov, Stable recovery of sparse overcomplete representations in the presence of noise, IEEE Trans. Inform. Theory 52(1), 6–18, 2006. https://doi.org/10.1109/TIT.2005.860430

Motivation

Quasi-Newton methods minimize a smooth function $f$ on $\mathbb{R}^n$ by moving along $d_k = -H_k g_k$, where $g_k$ is the gradient and $H_k$ is an approximation of the inverse Hessian built from observed gradient differences. The BFGS update is the most widely used way of building $H_k$, but it stores a dense $n \times n$ matrix, which is prohibitive for large $n$.

Nocedal's 1980 paper (Math. Comp. 35, 773–782) proposed keeping only the last $m$ correction pairs and rebuilding the matrix from a simple initial matrix $H_0$ at every step. The resulting method, called SQN in the paper, is now known as L-BFGS, and it is the default large-scale unconstrained optimizer in many numerical libraries and in machine learning. The paper's main theoretical claim is that this truncation does not destroy the finite termination of BFGS on quadratics.

Timeline:

• 1970: Broyden, Fletcher, Goldfarb and Shanno introduce the BFGS update (references [1] and [5] of the paper).
• 1977: Nazareth relates BFGS to conjugate gradients (Argonne Tech. Memo 282, reference [7]); his form of preconditioned conjugate gradients is the one the paper uses.
• 1977–1978: Shanno studies the memoryless BFGS update, the case $m = 1$ (reference [11]; journal version Math. Oper. Res. 3, 1978).
• 1980: Nocedal defines the special BFGS matrices and the SQN method and states that on quadratics with exact line searches SQN is identical to preconditioned conjugate gradients, hence has quadratic termination.
• 1989: Liu and Nocedal (Math. Programming 45) study the method, now called L-BFGS, for large-scale problems.
• 1998: Kolda, O'Leary and Nazareth (SIAM J. Optim. 8) treat limited-memory and update-skipping BFGS variants with exact line searches on quadratics.

Setting

Let $A$ be a symmetric positive definite $n \times n$ matrix and $b \in \mathbb{R}^n$, and let $f(x) = \tfrac12 x^T A x + b^T x$, a strictly convex quadratic with gradient $g(x) = Ax + b$ and unique minimizer $x^\ast = -A^{-1} b$.

Exact line search. Along a direction $d \neq 0$ from $x$, the step $\alpha = -g(x)^T d / d^T A d$ minimizes $f(x + \alpha d)$.

BFGS update. For a pair $(s, y)$ with $\rho = 1/y^T s$ and $v = I - \rho y s^T$, the BFGS update of $H$ is

$$\bar H = v^T H v + \rho s s^T .$$

Special BFGS matrices. Fix $H_0$ symmetric positive definite and a number $m \ge 1$ of stored corrections. Given pairs $(s_j, y_j)$, the special matrix $H_K$ is $H_0$ updated by the pairs $j = K - \min(K, m), \dots, K-1$, oldest first (the paper's (4)–(5)). Only the $m$ most recent pairs enter, and the matrix is rebuilt from $H_0$.

SQN. Starting from $x_0$, with $g_i = g(x_i)$:

$$d_i = -H_i g_i, \qquad x_{i+1} = x_i + \alpha_i d_i, \qquad s_i = x_{i+1} - x_i,\quad y_i = g_{i+1} - g_i,$$

with $\alpha_i$ the exact step and $H_{i+1}$ the special matrix built from the last $\min(i+1, m)$ pairs.

PCG with fixed preconditioner $H_0$. $d_0 = -H_0 g_0$, $x_{i+1} = x_i + \alpha_i d_i$, $d_{i+1} = -H_0 g_{i+1} + \beta_{i+1} d_i$ with $\beta_{i+1} = y_i^T H_0 g_{i+1} / y_i^T d_i$.

Formalization targets

Goal: quadratic termination of SQN

For every $n$, every symmetric positive definite $A$ and $H_0$, every $b$, $x_0$ and every $m \ge 1$,

$$\exists\, k \le n : \quad A x_k + b = 0 ,$$

where $x_k$ are the SQN iterates. The statement fixes no constant beyond the dimension bound $n$.

Milestones

1. Property (a): the special matrices are positive definite whenever $y_i^T s_i > 0$ for all $i$.
2. Eq. (7): along conjugate steps, $v_i y_i = 0$ and $v_i y_j = y_j$ for $i > j$.
3. Eq. (6): along conjugate steps, $H_k y_j = s_j$ for the $m$ most recent $j$ (when $k > m$).
4. Eq. (10): the special matrix equals $m$ sum-form BFGS corrections applied to $H_0$.
5. Eq. (15): the PCG directions satisfy $d_i^T y_j = 0$ for $i \neq j$.
6. Eq. (16): $g_i^T H_0 g_j = 0$ for $i \neq j$ and $g_i^T d_j = 0$ for $j < i$.
7. The PCG with fixed preconditioner $H_0$ reaches the minimizer in at most $n$ steps.
8. SQN and this PCG produce identical iterates and directions at every step.

Significance

The result shows that storing only $m$ correction pairs costs nothing on quadratics: for any $m \ge 1$, SQN terminates in at most $n$ steps, like full BFGS and conjugate gradients. It explains why L-BFGS with small $m$ is competitive, and it is the model case for later analyses of limited-memory methods (their linear convergence on uniformly convex functions, and their relation to Krylov methods). Property (b) is the reason one expects efficiency to grow with $m$: the matrix satisfies the secant equation on the $m$ most recent directions.

The claims are classical and generally accepted, but the paper argues them in a few lines ("it is straightforward to show"), deferring the PCG facts (15)–(16) to a reference. No machine-checked proof of the termination of BFGS, L-BFGS or preconditioned conjugate gradients is known to this mission. A formalization would provide a verified model of L-BFGS on quadratics and a reusable development of conjugate-direction methods with a preconditioner.

Difficulty

The obvious route, "SQN is BFGS and BFGS terminates", fails: SQN discards old corrections, so the classical BFGS argument (hereditary secant conditions on all past directions) does not apply once more than $m$ steps have been taken. The paper asserts the identity of SQN with preconditioned conjugate gradients in one sentence ("using a similar argument as for the SCG"), and the PCG relations it relies on are quoted from a technical report. The other difficulty is bookkeeping: the window of stored pairs shifts, the matrix is a nested product, and the runs must remain meaningful after the minimizer is reached.

Formalization scope

Vectors are Fin n → ℝ, matrices Matrix (Fin n) (Fin n) ℝ, $x^T y$ is dotProduct, and $s y^T$ is Matrix.vecMulVec. Symmetric positive definiteness is Matrix.PosDef. Indices are 0-based, as in the paper. The exact line search is the closed-form step $-g^T d / d^T A d$. The iterations have no stopping rule: once the gradient vanishes the direction and step are zero and the iterate stays at the minimizer (Lean's $0/0 = 0$). Past that point the zero pair stored by SQN leaves the BFGS step unchanged. The hypotheses are exactly the paper's: $A \succ 0$, $H_0 \succ 0$, $m \ge 1$ and exact line searches. $H_0$ need not be diagonal.

Two misprints are corrected and flagged in the items: the denominator of $\beta$ in (13) is $y_{i-1}^T d_{i-1}$ (as in (12) and p. 778), and the second relation of (16) is stated for $j < i$ (as used on p. 778), since it fails for $i < j$.

Ruled out: SQN is defined through its own matrices (4)–(5), rebuilt from $H_0$ and the last $m$ pairs. It is not defined through the PCG recurrence, not by one BFGS update of the previous matrix, and not with a stop rule that returns $-A^{-1}b$. The standing assumption $y_k^T s_k > 0$ is not a hypothesis of any statement about a run (it fails after termination and would make the goal vacuous). With $m = 0$ SQN is steepest descent and the goal is false, so $m \ge 1$ is required.

Needed infrastructure: algebra of rank-one updates and of Matrix.PosDef under congruence, conjugate-direction lemmas for quadratics, and the fact that $n+1$ mutually $H_0$-orthogonal vectors in $\mathbb{R}^n$ include a zero vector. The PCG results (milestones 5–7) are reusable beyond this mission. Proofs of any milestone, or of the goal directly, are welcome.

Selected references

• J. Nocedal, Updating Quasi-Newton Matrices with Limited Storage, Mathematics of Computation 35(151), 1980, 773–782. https://doi.org/10.1090/s0025-5718-1980-0572855-7
• D. F. Shanno, Conjugate gradient methods with inexact searches, Mathematics of Operations Research 3(3), 1978, 244–256. https://doi.org/10.1287/moor.3.3.244
• L. Nazareth, A Relationship Between the BFGS and Conjugate Gradient Algorithms, ANL-AMD Tech. Memo 282 (rev.), Argonne National Laboratory, 1977 (reference [7] of Nocedal 1980; no online copy located).
• T. G. Kolda, D. P. O'Leary, L. Nazareth, BFGS with update skipping and varying memory, SIAM Journal on Optimization 8(4), 1998, 1060–1083. https://doi.org/10.1137/S1052623496306450
• D. C. Liu, J. Nocedal, On the limited memory BFGS method for large scale optimization, Mathematical Programming 45, 1989, 503–528. https://doi.org/10.1007/BF01589116

Motivation

The conjugate gradient method (cg-method) of Hestenes and Stiefel is the standard iterative solver for linear systems $Ax=k$ with a symmetric positive definite matrix $A$. It is used for the large sparse systems of finite-element and finite-difference discretizations, as the inner solver of Newton-type and interior-point methods in optimization, and as the prototype of the Krylov subspace methods. Its original 1952 paper (Hestenes and Stiefel, J. Res. NBS 49(6), 1952) already presented it as two things at once: a direct method that reaches the exact solution in at most $n$ steps, and a method of successive approximations whose intermediate estimates are useful in their own right.

The second view needs a guarantee that the intermediate estimates actually approach the solution. The method is built to decrease the $A$-weighted error $f(x)=(h-x,A(h-x))$, and the residual $|k-Ax_i|$ need not decrease (Section 18 of the paper, p. 432, notes that it can increase at every step). Theorem 6:3 of the paper supplies the guarantee in the plain Euclidean length: the distance $|h-x_i|$ from the estimate to the solution decreases strictly at every step, by an exactly computable amount. This mission formalizes that theorem together with the relations from Sections 5 and 6 of the paper on which its proof rests.

Timeline:

• 1952: Hestenes and Stiefel introduce the method and prove, in one paper, finite termination (Theorems 4:2 and 5:2), the monotone decrease of the error function $f$ (Theorem 6:1), and the monotone decrease of the Euclidean error (Theorem 6:3). The later literature on cg as an iterative method for large sparse systems takes these properties as its starting point.

Setting

Let $A$ be a real $n\times n$ matrix that is symmetric and positive definite, let $k\in\mathbb{R}^n$, and let $h$ be the solution of $Ah=k$. Write $(x,y)=x_1y_1+\cdots+x_ny_n$ and $|x|^2=(x,x)$. From an arbitrary starting point $x_0$, the cg-method (5:1) computes estimates $x_i$, residuals $r_i$ and direction vectors $p_i$ by

$$p_0=r_0=k-Ax_0,\quad a_i=\frac{|r_i|^2}{(p_i,Ap_i)},\quad x_{i+1}=x_i+a_ip_i,\quad r_{i+1}=r_i-a_iAp_i,\quad b_i=\frac{|r_{i+1}|^2}{|r_i|^2},\quad p_{i+1}=r_{i+1}+b_ip_i .$$

The error vector of $x_i$ is $y_i=h-x_i$. The error function (4:5) is $f(x)=(h-x,A(h-x))$, which is nonnegative and vanishes only at $x=h$. The Rayleigh quotient (4:12) of a vector $z\neq 0$ is $\mu(z)=(z,Az)/|z|^2$. The Lean development names these cgIter A k x₀ i (with fields .x, .r, .p), cgAlpha for $a_i$, errorFun A h x and rayleigh A z.

Formalization targets

Goal: Theorem 6:3

For every step that the method performs, that is, every $i$ with $r_i\neq 0$,

$$|y_i|^2-|y_{i+1}|^2=\frac{f(x_{i+1})+f(x_i)}{\mu(p_i)}\qquad\text{and}\qquad |y_{i+1}|<|y_i| .$$

The paper writes the step from $x_{i-1}$ to $x_i$; the Lean statement shifts the index by one. The goal holds for every dimension $n$, every symmetric positive definite $A$, every $k$ and every $x_0$.

Milestones

1. Theorems 4:2 and 5:2: some $m\le n$ has $x_m=h$.
2. Theorem 5:3, (5:6a): $(p_i,p_j)=|r_j|^2|p_i|^2/|r_i|^2$ for $i\le j$.
3. Theorem 6:1, (6:1): $f(x_i)-f(x_{i+1})=a_i|r_i|^2=\mu(p_i)|x_i-x_{i+1}|^2$.
4. Theorem 6:1, (6:2): $f(x_i)-f(x_j)=\sum_{l=i}^{j-1}a_l|r_l|^2$ for $i<j$.
5. Section 6, (6:6): $(y_{i+1},x_{i+1}-x_i)=f(x_{i+1})/\mu(p_i)$.

Significance

The theorem is what makes an early stop of the cg-method safe in the norm a user usually cares about. Every intermediate estimate is closer to the solution, in Euclidean distance, than the previous one, and the identity (6:5) states by how much. It also separates the cg-method from methods that minimize the residual: the $A$-norm error, the Euclidean error and the residual behave differently, and only the first two are monotone along cg.

The results are proved in the 1952 paper. They have not been formalized: the Prove2Me library has no statement of the conjugate gradient recursion (5:1), and Mathlib has none either. What this mission adds is a machine-checked version of the paper's Section 6 argument for the recursion exactly as printed, including the case analysis at termination that the paper leaves implicit, and a reusable Lean definition of the cg iteration with its basic identities.

Difficulty

The obvious argument does not reach the conclusion. The method decreases $f(x)=(y,Ay)$ at every step, but a decrease in this $A$-weighted norm does not imply a decrease in the Euclidean norm: for a single step along an arbitrary direction, even the best step for $f$ can lengthen the Euclidean error. So the theorem cannot be proved one step at a time from the local minimization property. It depends on how the current direction relates to all the later directions of the same run, and those relations in turn rest on the mutual orthogonality of the residuals and the conjugacy of the directions, which are established by an induction over the whole run.

A second difficulty is bookkeeping at the end of the run. The recursion divides by $|r_i|^2$ and by $(p_i,Ap_i)$, which vanish after termination. Every milestone has to hold, or be guarded, past that point, and the goal needs the hypothesis $r_i\neq 0$ exactly because the strict inequality fails once $x_i=h$.

Formalization scope

Vectors are Fin n → ℝ, the scalar product is dotProduct (⬝ᵥ), $Ax$ is Matrix.mulVec (*ᵥ), and the standing assumption is A.PosDef, which in Mathlib includes symmetry. The solution $h$ is a variable with the hypothesis A *ᵥ h = k. Indices are 0-based. The cg recursion is the definition cgIter, which computes (5:1b)–(5:1f) literally and in order; it has no stopping rule, and Lean's convention $t/0=0$ makes it stay at $h$ with $r_i=p_i=0$ once $r_m=0$. Lengths appear squared, as $(y,y)$. The milestones are stated for every index without a termination guard, because both sides of each identity vanish after termination; only the goal carries $r_i\neq 0$.

Two formalizations would trivialize the goal and are ruled out. The goal does not assume termination ($x_m=h$) or any bound on $i$: it quantifies over every cg run and every step that takes place. And it is about the Euclidean length $|h-x_i|$, not the error function $f$ (that is Theorem 6:1, a different and weaker statement) and not the residual.

A complete development needs Theorem 5:1 (orthogonality of residuals, conjugacy of directions) for the literal recursion, the identities (5:2) and (5:3c), and the positivity of $(p,Ap)$ for $p\neq 0$. These are reusable for any further work on the cg-method, including the sister mission on finite termination. Proofs of any milestone, and alternative proofs of Theorem 6:3 through the Krylov-subspace characterization, are welcome.

Selected references

• M. R. Hestenes and E. Stiefel, Methods of Conjugate Gradients for Solving Linear Systems, J. Res. Natl. Bur. Stand. 49(6), 409–436, 1952. https://doi.org/10.6028/jres.049.044 (publisher's scan: https://nvlpubs.nist.gov/nistpubs/jres/049/jresv49n6p409_A1b.pdf)

Motivation

Solving a linear system $Ax = k$ with a large symmetric positive definite matrix $A$ is a basic task of scientific computing: it arises from discretized elliptic equations, least-squares problems and the Newton steps of optimization methods. In 1952 Magnus Hestenes and Eduard Stiefel published the conjugate gradient method (cg-method) (J. Res. Natl. Bur. Stand. 49 (1952) 409–436). The method uses $A$ only through matrix–vector products and stores a few vectors. The paper's abstract states its central property in one sentence: "The solution is given in $n$ steps."

The paper obtains this property from a more general scheme, the method of conjugate directions (cd-method), which also contains Gaussian elimination as a special case. Its argument has two parts. Every cd-method with nonzero directions reaches the solution within $n$ steps (Theorem 4:2). The cg-method is a cd-method (Theorem 5:2), which follows from the orthogonality and conjugacy relations of Theorem 5:1. This mission formalizes that chain.

Setting

Vectors are real $n$-tuples, with scalar product $(x, y) = x_1y_1 + \cdots + x_ny_n$ and squared length $|x|^2 = (x, x)$. The matrix $A$ is real, $n \times n$, symmetric and positive definite, which is the paper's standing assumption (p. 410). The solution $h$ satisfies $Ah = k$. The residual of an estimate $x$ is $r = k - Ax$. Two vectors $x, y$ are conjugate when $(x, Ay) = 0$.

The cg-method (eq. (3:1), p. 411) starts from an arbitrary estimate $x_0$ and sets $p_0 = r_0 = k - Ax_0$. Given $x_i$, $r_i$, $p_i$, it computes

$$a_i = \frac{|r_i|^2}{(p_i, Ap_i)},\quad x_{i+1} = x_i + a_i p_i,\quad r_{i+1} = r_i - a_i Ap_i,\quad b_i = \frac{|r_{i+1}|^2}{|r_i|^2},\quad p_{i+1} = r_{i+1} + b_i p_i.$$

In Lean the iterates are cgIter A k x₀ i, a structure with fields x, r, p. The step length is cgA.

The cd-method (Section 4, p. 412) chooses an arbitrary first direction $p_0$ and then sets $x_{i+1} = x_i + a_i p_i$ with $a_i = (p_i, r_i)/(p_i, Ap_i)$ and $r_i = k - Ax_i$. Each new direction $p_{i+1}$ may be any vector conjugate to $p_0, \dots, p_i$. Because the directions are free, a cd-run is a property of sequences: IsCDRun A k x r p.

Formalization targets

Goal: finite termination of the cg-method

For every $n$, every symmetric positive definite $A$, every $k$ and $h$ with $Ah = k$, and every initial estimate $x_0$,

$$\exists\, m \le n:\quad x_m = h,$$

where $x_m$ is the $m$-th cg iterate. This is the statement of the abstract and of Section 3 (p. 410): "one will reach an estimate $x_m$ ($m \le n$) at which $r_m = 0$. This estimate is the desired solution $h$."

Milestones

1. Theorem 4:1 (p. 412). For every cd-run, the directions are mutually conjugate (4:3a). The residual $r_i$ is orthogonal to $p_0, \dots, p_{i-1}$ (4:3b). The products $(p_i, r_j)$ are the same for all $j \le i$ (4:3c). Hence $a_i = (p_i, r_0)/(p_i, Ap_i)$ (4:4).
2. Theorem 4:2 (p. 412). Every cd-run whose directions $p_0, \dots, p_{n-1}$ are nonzero has $x_m = h$ for some $m \le n$.
3. Theorem 5:1 (p. 414), in four items. For the cg-method:
   • (5:3a) $(r_i, r_j) = 0$ for $i \ne j$;
   • (5:3b) $(p_i, Ap_j) = 0$ for $i \ne j$;
   • (5:3c) $(p_i, r_j) = 0$ for $i < j$ and $(p_i, r_j) = |r_i|^2$ for $i \ge j$;
   • (5:3d) $(r_i, Ap_i) = (p_i, Ap_i)$, and $(r_i, Ap_j) = 0$ for $i \ne j, j + 1$.
4. Theorem 5:5, eq. (5:10) (p. 416):

$$a_i = \frac{|r_i|^2}{(p_i, Ap_i)} = \frac{(p_i, r_i)}{(p_i, Ap_i)} = \frac{(p_i, r_0)}{(p_i, Ap_i)}.$$

5. Theorem 5:2, first sentence (p. 415). The cg-method is a cd-method: its iterates satisfy IsCDRun.

Significance

Finite termination is the property that distinguishes the conjugate gradient method from stationary iterations such as Jacobi or Gauss–Seidel. It explains why the method can be used as a direct solver in exact arithmetic. It is the starting point for the later theory of Krylov subspace methods. The relations of Theorem 5:1 are the checks the paper recommends for monitoring a computation (eq. (3:3)). They are also the input to the paper's further results on the monotone decrease of the error (Section 6) and on the connection with orthogonal polynomials (Sections 14–18).

The result has been proved since 1952 and appears in every numerical linear algebra textbook. As far as a search of the Prove2Me catalog shows (September 2026), it has no machine-checked proof there. The only conjugate gradient material on the platform is a Hilbert-space convergence result for a different recurrence. This mission adds a faithful formal version of the original recursion (3:1) and of the cd-method. It also adds the complete termination argument, organized as in the paper. The definitions and the Theorem 5:1 relations are reusable by any later formalization of Krylov methods, including the second mission of this series on the decrease of the error $|h - x_i|$.

Difficulty

The obvious argument says that the residuals are mutually orthogonal, so at most $n$ of them are nonzero. That argument is only as good as the orthogonality, and Theorem 5:1 must be established by a simultaneous induction over four families of relations. The recursion defines $r_{i+1}$ by an update, not as $k - Ax_{i+1}$, so even $r_i = k - Ax_i$ needs a proof. Orthogonality of $r_{i+1}$ to the earlier residuals needs the conjugacy of the earlier directions, and conjugacy of $p_{i+1}$ needs the orthogonality of the earlier residuals. Neither family can be proved first.

A second obstacle is the passage from orthogonality to termination. A cd-run with a zero direction stalls, so Theorem 4:2 needs the directions $p_0, \dots, p_{n-1}$ to be nonzero. The cg directions become zero exactly when the solution is reached, so applying Theorem 4:2 to cg needs a case split at the first vanishing residual.

Formalization scope

Vectors are Fin n → ℝ, the scalar product is dotProduct (⬝ᵥ), and $Ax$ is Matrix.mulVec (*ᵥ). The hypothesis on $A$ is Mathlib's Matrix.PosDef, which includes symmetry. The solution enters only through the hypothesis A *ᵥ h = k. Indices are 0-based, as in the paper.

The cg iteration is total and has no stopping test. Once $r_m = 0$, Lean's convention $x/0 = 0$ gives $p_m = 0$ and $a_m = 0$. From then on the iteration stays at $x_m$, with zero residuals and directions. For this reason the relations of Theorem 5:1 are stated for all indices without guards: after termination they hold trivially.

Two formulations would make the goal trivial, and neither is used. One is a stopping test or step that refers to $h$ or to $A^{-1}k$. The other replaces (3:1b), (3:1d) or (3:1e) by the equivalent formulas (3:2a), (3:2b) or $r_{i+1} = k - Ax_{i+1}$, which would move Theorem 5:5 and part of Theorem 5:2 into the definition. The iteration is (3:1) literally.

The cd-method's implicit hypothesis, that the directions $p_0, \dots, p_{n-1}$ are nonzero, is an explicit binder of Theorem 4:2. Without it the statement fails (take $p_0 = 0$). The hypothesis is satisfiable: the cg run with nonzero residuals is one example.

A complete development needs the standard facts that mutually conjugate nonzero vectors are linearly independent and that $n$ independent vectors span $\mathbb{R}^n$. It also needs the induction behind Theorem 5:1. Contributions welcome beyond the milestones include the converse half of Theorem 5:2, the relation (5:2) expressing $p_k$ through $r_0, \dots, r_k$, and Theorem 4:5 (the cd-method computes $A^{-1}$).

Selected references

• M. R. Hestenes and E. Stiefel, Methods of Conjugate Gradients for Solving Linear Systems, Journal of Research of the National Bureau of Standards 49(6), 409–436, 1952. https://doi.org/10.6028/jres.049.044
• L. Fox, H. D. Huskey and J. H. Wilkinson, Notes on the solution of algebraic linear simultaneous equations, Quarterly Journal of Mechanics and Applied Mathematics 1(1), 149–173, 1948 (the cd-method from a different point of view; cited in the paper's Section 4 footnote). https://doi.org/10.1093/qjmam/1.1.149
• G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013, §11.3 (textbook account of the method).