Numerical Analysis

Motivation

Every numerical computation is performed in a finite arithmetic: a real number is stored with a fixed number of decimal digits, and each stored value carries an error. Before any method is analysed, a course in numerical methods has to say how large that error can be and how it grows when approximate values are combined. The answers are the two rounding bounds — $10^{1-t}$ for truncation and $\\tfrac12 10^{1-t}$ for symmetric rounding with $t$ digits — and the propagation formulas for sums, products and quotients.

This mission is the eighth in a series formalizing the lecture notes Métodos Numéricos by Sergio Roberto de Freitas (UFMS, 2000). It covers Chapter 2, Erros.

Setting

Every nonzero real $x$ can be written in normalized decimal form with $t$ retained digits,

$$x = f \\times 10^{e} + g \\times 10^{e-t}, \\qquad 0.1 \\le |f| < 1, \\quad 0 \\le g < 1,$$

where $e$ is the decimal exponent of $x$ and $f$ its normalized mantissa.

Truncated rounding keeps $\\tilde{x} = f\\times 10^{e}$, discarding the tail; the error is $e_x = x - \\tilde{x} = g\\times 10^{e-t}$.

Symmetric rounding keeps $\\tilde{x} = f\\times 10^{e}$ when $g < 0.5$ and $\\tilde{x} = f\\times10^{e} + 10^{e-t}$ when $g \\ge 0.5$; this is the familiar school rule of rounding the $(t+1)$-st digit up when it is at least $5$.

If $\\tilde{x}$ approximates $x$ and $\\tilde{y}$ approximates $y$, the absolute errors are $e_x = x - \\tilde{x}$ and $e_y = y - \\tilde{y}$, and the chapter computes the error of the four arithmetic operations in terms of them.

Target

The goal theorem is Proposição 2.4.3: with symmetric rounding to $t \\ge 1$ digits, the relative error of a stored positive number satisfies

$$\\frac{|x - \\tilde{x}|}{|\\tilde{x}|} \\;\\le\\; \\frac{1}{2}\\,10^{1-t}.$$

The milestones are the companion bound for truncated rounding, $10^{1-t}$ (Proposição 2.4.2), and the propagation identities of §2.5.1 for the sum and difference, for the product and for the quotient.

Significance

These bounds are the unit of account for every later error analysis: the machine epsilon of a $t$-digit decimal arithmetic is exactly the constant appearing in the goal, and the factor two between truncation and symmetric rounding is the reason the latter is the default. The propagation identities explain the two classical failure modes of floating-point computation — cancellation in a difference of nearly equal numbers, where the absolute errors survive while the result shrinks, and the amplification of a relative error by division by a small number.

Difficulty

The difficulty is entirely in the normalization. The bounds follow from two facts about the decimal decomposition — that the discarded part is at most one unit in the last retained place, and that the retained part is at least $10^{e-1}$ — and both have to be derived from the definition of the decimal exponent through the floor function and the logarithm, including the boundary case where $x$ is an exact power of ten and the case where symmetric rounding carries into a new leading digit. The propagation statements, by contrast, are algebraic identities.

Formalization scope

The decimal exponent is defined as $\\lfloor \\log_{10}|x| \\rfloor + 1$, and the two roundings are defined by scaling by $10^{t-e}$, applying the floor function (with an added $\\tfrac12$ in the symmetric case, so that halves round up) and scaling back. This is a total definition for every real number; the theorems are stated for $x > 0$ and $t \\ge 1$, which is the situation the source describes, and no claim is made about negative arguments or about $x = 0$.

Two points of deviation should be checked by the auditor. First, the source states both bounds with a strict inequality; for symmetric rounding equality is attained when the discarded tail is exactly one half unit in the last place, so the formal statement uses a non-strict inequality, while the truncation bound remains strict. Second, the propagation formulas for the product and the quotient are stated here as exact identities: the source's versions drop the second-order term $e_xe_y$ and the higher terms of a geometric series, which is a legitimate approximation but not an equality.

Selected references

• S. R. Freitas, Métodos Numéricos, Departamento de Computação e Estatística, Universidade Federal de Mato Grosso do Sul, 2000. Chapter 2, Erros, pp. 17–31. (Course notes supplied with this mission.)

Motivation

The conjugate-gradient method in Luenberger's Chapter 10 is one of the most enduring consequences of Hilbert-space geometry in numerical optimization. For a bounded self-adjoint coercive operator, it solves the quadratic first-order equation Q x = b using only operator applications, inner products, and a short recurrence. Luenberger develops the method from steepest descent and conjugate directions, then proves convergence in a general real Hilbert space rather than only for finite matrices. This mission formalizes that full setting. It also repairs a practical omission in the printed recursion: division formulas are undefined after exact convergence, so the formal algorithm explicitly stops and stutters once its search direction is zero.

Setting

Let H be a complete real inner-product space and Q : H →L[ℝ] H a bounded self-adjoint operator. Constants m and M satisfy 0 < m ≤ M and

$$m\lVert x\rVert^2 \le \langle x,Qx\rangle \le M\lVert x\rVert^2$$

for every x. The first inequality is coercivity; together with self-adjointness it supplies the positive Q-energy. For a right-hand side b and initial point x₀, the initial residual and direction are both b - Q x₀. A conjugate-gradient state records the current iterate, residual, and direction. If the direction is nonzero, the next state uses Luenberger's alpha and beta ratios. If the direction is zero, conjugateGradientStep returns the same state, so every natural-number iterate is total and all denominators occur only on the active branch.

Formalization targets

The root theorem VectorSpaceOpt.conjugate_gradient_converges states that there is a unique xStar satisfying Q xStar = b and that the iterate component of the guarded conjugate-gradient state tends to xStar in norm. Four milestones provide reusable structure. coercive_selfadjoint_bijective establishes existence and uniqueness for Q x = b from bounded self-adjoint coercivity. conjugate_directions_converge formalizes §10.6, Theorem 1: a complete sequence of nonzero pairwise Q-orthogonal directions produces residuals orthogonal to every earlier direction and iterates converging to the solution. cg_directions_conjugate_until_stop records the §10.8 invariants only before the explicit stopping time. cg_energy_contraction captures the uniform energy reduction factor derived from the bounds m and M.

The total algorithm is represented by conjugateGradientIterate, and its error functional is

$$E(x)=\langle x-x^*,Q(x-x^*)\rangle.$$

These definitions are proposed as mission-owned reusable objects in the shared VectorSpaceOpt namespace.

Significance

The mission gives a coordinate-free verification target for an algorithm usually presented through arrays and matrices. Its theorem applies directly to finite-dimensional symmetric positive-definite systems but also retains Luenberger's infinite-dimensional perspective. The guarded recursion is suitable for later executable specializations and makes exact termination a first-class semantic event. The coercivity and conjugate-directions milestones can be reused for Galerkin methods, preconditioned variants, and other Krylov algorithms, while the energy estimate provides a natural connection to condition-number convergence rates.

Unlike a matrix-only formalization, the Hilbert-space theorem cleanly separates the geometric reason for convergence from any storage representation. It therefore complements Mathlib's existing operator and orthogonality libraries and can serve as a specification against which finite implementations are later verified. It also preserves the book's unifying theme: optimization algorithms arise from the geometry of carefully chosen inner products rather than from coordinate manipulation alone.

Difficulty

The difficulty is medium to high. Algebraic invariants of the three-term recurrence involve several interacting orthogonality relations and require strict control of nonzero denominators. Infinite-dimensional convergence additionally uses density of the closed span of directions and comparison of the Q-energy with the ambient norm. The theorem must move between self-adjoint continuous linear maps, scalar inner products, filters on sequences, and function iteration. Exact termination creates a case split that informal accounts routinely ignore; the formal statement must show that the zero-direction branch is stable and already represents the solution.

Formalization scope

The proposal follows §10.6 and §10.8, pp. 291–296, and uses Chapter 10, Problem 10 on p. 309 for the coercive-invertibility dependency. All assumptions on Q, m, and M that §10.8 inherits from the preceding sections are repeated explicitly. The conjugate-directions milestone explicitly assumes every direction is nonzero and that the closed span of the directions is the whole Hilbert space. The conjugate-gradient invariants are asserted only for iterations before a zero direction occurs. Once it occurs, the state stutters by definition; the proposal never relies on Lean's totalized value for 0 / 0.

Luenberger's §10.7, Theorem 1 is not included as a literal milestone. As printed, its orthogonalization-of-moments statement omits self-adjointness of the auxiliary operator relative to the Q inner product and omits the linear-independence/nonbreakdown conditions needed to keep denominators nonzero. The mission instead isolates the Q-conjugacy invariant directly from §10.8. It does not claim finite-dimensional termination within dim H steps, floating-point stability, preconditioning, a sharp Chebyshev condition-number rate, or computability of equality tests on arbitrary Hilbert spaces.

Selected references

• David G. Luenberger, Optimization by Vector Space Methods, Wiley, 1969, Chapter 10, §10.6, Theorem 1, pp. 291–292; §10.8, Theorem 1, pp. 294–296; Problem 10, p. 309. Scan: https://sites.science.oregonstate.edu/~show/old/142_Luenberger.pdf
• Lean community, Mathlib documentation, continuously updated: https://leanprover-community.github.io/mathlib4_docs/ (real inner-product spaces, continuous linear maps, coercivity, closed spans, orthogonality, and filter convergence).

Motivation

Quadratic exterior penalties in Luenberger's §10.11 replace a constrained problem by a sequence of unconstrained minimizations. The method is simple enough to state in a few lines, yet Luenberger's convergence theorem is strikingly general: no convexity, differentiability, or convergence of the full minimizer sequence is required. If penalty weights increase to infinity and a subsequence of exact penalty minimizers converges, lower semicontinuity alone makes its limit feasible and optimal. This mission isolates that robust primal convergence result as a tractable companion to the more analytic conjugate-gradient and optimal-control missions. It offers a clean formalization target with direct relevance to nonlinear programming and approximation schemes.

Setting

Let X be a topological space, f : X → ℝ, and G : X → (Fin p → ℝ). Feasibility means G x i ≤ 0 for every component. Define the positive part componentwise and the squared violation by

$$G_i^+(x)=\max(0,G_i(x)), \qquad v(x)=\sum_i (G_i^+(x))^2.$$

For a positive weight K, the penalty objective is f x + K * v x. A sequence K n is positive, nondecreasing, and tends to +∞. The constrained problem is assumed to have a minimizer xStar, and for each n an exact global minimizer x n of the corresponding penalty objective is supplied. A limit point is represented explicitly by a strictly increasing index map phi for which x ∘ phi tends to x₀.

Formalization targets

The root VectorSpaceOpt.quadratic_penalty_cluster_point_converges formalizes §10.11, Theorem 1. Assuming lower semicontinuity of f and v, it concludes that every stated subsequential limit x₀ is feasible, has the same objective value as xStar, and globally minimizes f over the feasible set.

Three milestones split the exact source content into reusable statements. quadratic_penalty_basic_estimates is §10.11, Lemma 1: the attained penalty values are nondecreasing, are bounded above by f xStar, and the stronger weighted violation K n * v (x n) tends to zero. penalty_cluster_point_feasible combines convergence of violations with lower semicontinuity at a subsequential limit to recover all component inequalities. penalty_cluster_point_optimal combines lower semicontinuity of f, the uniform upper bound f (x n) ≤ f xStar, feasibility of the limit, and optimality of xStar to identify the limiting objective value and global constrained optimality.

Significance

The theorem captures the essential consistency guarantee behind one of the most widely used constraint-handling methods. Its assumptions separate optimization existence from convergence: minimizers of each auxiliary problem and at least one cluster point are assumed, while the theorem identifies what any such cluster point must be. The componentwise positive-part and violation definitions are reusable for augmented Lagrangians, exact penalties, barrier comparisons, and finite inequality systems. The basic-estimates lemma is particularly useful because it requires neither topology nor continuity and exposes a quantitative fact stronger than mere feasibility residual convergence.

Because the proof target is stated over an arbitrary topological space, the mission also clarifies which parts of penalty convergence are genuinely metric and which depend only on order, finite nonnegative sums, and lower semicontinuity. This abstraction is faithful to the source's vector-space viewpoint.

Difficulty

The mission has moderate difficulty and relatively low infrastructure risk. The main analytic interfaces are lower semicontinuity along a convergent subsequence and real filter convergence to both zero and infinity. The basic estimates require reasoning simultaneously about minimizers for changing objectives, monotonicity of the weights, and the asymptotic product K n * v (x n). The cluster-point theorem must extract componentwise feasibility from a finite sum of nonnegative squares without assuming continuity of G. Lean's IsMinOn does not itself assert membership in the feasible set, so feasibility of the known constrained minimizer is included separately rather than hidden in prose.

Formalization scope

The proposal covers the primal part of §10.11: Lemma 1 on p. 305 and Theorem 1 on p. 306. It makes “limit point” precise through a strictly monotone subsequence, avoiding any assumption that the full sequence converges. The weight sequence may have repeated values because the source only needs it to be nondecreasing, but every weight is positive and the sequence tends to atTop. Lower semicontinuity is required for f and the composite violation v, exactly as in the book; continuity or componentwise lower semicontinuity of G is not substituted. Existence of xStar and of every penalty minimizer is assumed rather than derived from compactness or coercivity.

The mission does not include §10.11, Lemma 2 or Theorem 2 on dual multipliers. Those results add convexity and continuity assumptions and naturally require careful treatment of an extended-real dual functional. It also does not address approximate minimizers, rates, boundedness of the sequence, existence of cluster points, equality constraints beyond their encoding as paired inequalities, or finite exactness. Keeping those extensions separate preserves the unusually weak hypotheses and clear conclusion of the cited primal theorem.

Selected references

• David G. Luenberger, Optimization by Vector Space Methods, Wiley, 1969, Chapter 10, §10.11, Lemma 1 and Theorem 1, pp. 305–306. Scan: https://sites.science.oregonstate.edu/~show/old/142_Luenberger.pdf
• Lean community, Mathlib documentation, continuously updated: https://leanprover-community.github.io/mathlib4_docs/ (lower semicontinuity, finite sums, Fin-indexed vectors, subsequences, global minima on sets, and filter convergence).

Motivation

A tall $N\times n$ matrix with independent standard normal entries is, with high probability, a near-isometry: all its singular values lie in $[\sqrt N-\sqrt n,\ \sqrt N+\sqrt n]$ up to Gaussian fluctuations. This is Gordon's theorem together with Gaussian concentration (Vershynin 2012, Theorem 5.32 and Corollary 5.35), the non-asymptotic counterpart of the Bai–Yin law and the model result for the whole theory of non-asymptotic random matrices. In randomized numerical linear algebra the same two ingredients give the expected spectral norm of $S\Gamma T$ (Halko–Martinsson–Tropp 2011, Proposition 10.1 (10.2)), its second moment (Tropp–Webber 2023, Lemma B.1), and the lower tail of $\sigma_{\min}$ of a Gaussian sketch, which controls the inverse moments $\mathbb E\,\sigma_{\min}^{-2p}$ that appear in the analysis of sketching algorithms with rank truncation. None of these results is available in Mathlib; this mission proves them.

Setting

A standard Gaussian matrix $A\in\mathbb R^{N\times n}$ has independent $\mathcal N(0,1)$ entries; its law is $\gamma_{N,n}$. Its extreme singular values are

$$s_{\min}(A)=\inf_{\|x\|_2=1}\|Ax\|_2,\qquad s_{\max}(A)=\sup_{\|x\|_2=1}\|Ax\|_2=\|A\|,$$

the infimum and supremum over the unit sphere of $\mathbb R^n$ ($n\ge1$); $\|A\|$ is the spectral norm. For a wide matrix $G\in\mathbb R^{k\times t}$, $k\le t$, the usual $\sigma_{\min}(G)$ is $s_{\min}(G^{\mathsf T})$. The Frobenius norm is $\|A\|_F=(\sum_{i,j}A_{ij}^2)^{1/2}$; a function $h$ on matrices is $L$-Lipschitz if $|h(X)-h(Y)|\le L\|X-Y\|_F$. A chi-square variable with $d$ degrees of freedom is $\chi_d^2=\sum_{i=1}^d g_i^2$ for independent standard normals $g_i$. The letter $e$ denotes Euler's number.

Formalization targets

Goal — deviation of the extreme singular values (Vershynin, Corollary 5.35)

For every $t\ge0$, with probability at least $1-2e^{-t^2/2}$,

$$\sqrt N-\sqrt n-t\ \le\ s_{\min}(A)\qquad\text{and}\qquad s_{\max}(A)\ \le\ \sqrt N+\sqrt n+t .$$

This is the goal because it is the statement users of Gaussian matrices quote: a quantitative near-isometry with explicit sub-Gaussian deviations. It follows from the expectation bounds (Gordon) and concentration, both of which are milestones.

Milestones

• Gaussian concentration for Lipschitz functions of a Gaussian matrix: $\mathbb P\{h(G)\ge\mathbb E h(G)+Lt\}\le e^{-t^2/2}$ (HMT Proposition 10.3; Vershynin Proposition 5.34).
• Gordon–Chevet bound. $\mathbb E\|S\Gamma T\|\le\|S\|\|T\|_F+\|S\|_F\|T\|$ (HMT Proposition 10.1 (10.2), A.2) and its second-moment form $(\mathbb E\|S\Gamma T\|^2)^{1/2}\le\|S\|\|T\|_F+\|S\|_F\|T\|$ (Tropp–Webber Lemma B.1).
• Gordon's theorem. $\sqrt N-\sqrt n\le\mathbb E s_{\min}(A)\le\mathbb E s_{\max}(A)\le\sqrt N+\sqrt n$ (Vershynin Theorem 5.32).
• One-sided tail of $s_{\min}$. $\mathbb P\{s_{\min}(A)\le\sqrt N-\sqrt n-u\}\le e^{-u^2/2}$ for $u\ge0$ (Davidson–Szarek Theorem II.13; one-sided form of Vershynin Corollary 5.35).
• Chi-square lower tail. $\mathbb P\{\chi_d^2\le u\}\le(eu/d)^{d/2}$ for $0\le u\le d$ (Dasgupta–Gupta Lemma 2.2).
• Small-ball bound for $s_{\min}$. $\mathbb P\{s_{\min}(A)\le s\}\le n\bigl(\frac{e\,ns^2}{N-n+1}\bigr)^{(N-n+1)/2}$ for $s\ge0$, $ns^2\le N-n+1$, from the invertibility-via-distance lemma of Rudelson–Vershynin (Lemma 3.5), rotational invariance and the chi-square lower tail.

Significance

The results. Gordon's theorem is sharp: $\mathbb E s_{\min}$ and $\mathbb E s_{\max}$ are pinned to within $O(1)$ of the Bai–Yin limits for every $N,n$. Gaussian concentration upgrades it to tail bounds with the optimal sub-Gaussian rate, and the Chevet-type bound is the same mechanism for $S\Gamma T$, where it yields the expected spectral error of the randomized SVD. The small-ball bound covers the other regime, $s\to0$, with polynomial decay of the exact order $N-n+1$; together the two tails make every negative moment $\mathbb E s_{\min}(A)^{-2p}$, $2p<N-n+1$, finite and explicitly bounded, which is what the early-truncated analysis needs.

Formalizing them. Mathlib has the standard Gaussian measure on finite-dimensional inner product spaces but neither Gaussian concentration for Lipschitz functions, nor Slepian's or Gordon's comparison inequalities, nor any result on singular values of random matrices. These are foundational pieces of high-dimensional probability whose formalization would be reused far beyond this series.

Difficulty

The upper bound in Gordon's theorem is proved by Slepian's comparison lemma for the Gaussian process $(u,v)\mapsto\langle Au,v\rangle$ on $S^{n-1}\times S^{N-1}$, and the lower bound by Gordon's minimax comparison inequality for $\max_v\min_u\langle Au,v\rangle$; both require a theory of Gaussian processes indexed by continuous parameter sets, or at least a discretization argument with careful limits. Gaussian concentration (Tsirelson–Ibragimov–Sudakov) is usually proved through the Gaussian logarithmic Sobolev inequality or the Ornstein–Uhlenbeck semigroup, with a smoothing step to pass from smooth to Lipschitz functions. Both are long developments with no Mathlib support. The elementary alternatives do not reach the stated constants: an $\varepsilon$-net argument bounds $\mathbb E\|A\|$ by $C(\sqrt N+\sqrt n)$ with $C>1$, and a union bound over rows for the small-ball estimate is exactly how the stated bound is obtained, but it needs the fact that the distance from a Gaussian vector to a fixed $(n-1)$-dimensional subspace is a $\chi_{N-n+1}$ variable, which again rests on rotational invariance.

Formalization scope

Matrices are sampled from the product measure $\gamma_{N,n}$ on arrays. $s_{\min}$ is sMin, the real infimum over the unit sphere of $\mathbb R^n$; $n\ge1$ is assumed so the sphere is nonempty, and $N\ge n$ is not assumed for Gordon's theorem (for $N<n$ one has $s_{\min}=0$ and the lower bound is trivial, as in the source). $s_{\max}$ is Mathlib's $\ell_2$ operator norm. Expectation bounds are stated together with integrability of their integrands, so the Lean convention that non-integrable functions integrate to $0$ cannot trivialize them. Probabilities are measures in $[0,\infty]$ compared with $\operatorname{ofReal}$ of the real bound; the goal's lower bound $1-2e^{-t^2/2}$ uses truncated subtraction and is vacuous for $t\le\sqrt{2\log2}$, exactly as in the source. Gaussian concentration is stated with $L>0$ and $t\ge0$, which the sources leave implicit (for $L=0$ or $t<0$ the event has probability one). Results that the sketching literature states for wide $k\times t$ matrices are stated here for tall $N\times n$ matrices via $\sigma_{\min}(G)=s_{\min}(G^{\mathsf T})$, with $(N,n)=(t,k)$.

Welcome contributions: Gaussian concentration for Lipschitz functions on $\mathbb R^d$ (any proof), Slepian's and Gordon's inequalities, the $1$-Lipschitz property of $s_{\min}$ and $s_{\max}$, the chi-square Chernoff bounds, and the law of the distance from a Gaussian vector to a subspace.

Selected references

• N. Halko, P.-G. Martinsson, J. A. Tropp, Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions, SIAM Review 53(2), 217–288, 2011. https://arxiv.org/abs/0909.4061
• J. A. Tropp, R. J. Webber, Randomized algorithms for low-rank matrix approximation: design, analysis, and applications, 2023. https://arxiv.org/abs/2306.12418
• R. Vershynin, Introduction to the non-asymptotic analysis of random matrices, in Compressed Sensing: Theory and Applications, Cambridge University Press, 2012. https://arxiv.org/abs/1011.3027
• L. Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika 12, 134–139, 1918. https://doi.org/10.1093/biomet/12.1-2.134
• S. Dasgupta, A. Gupta, An elementary proof of a theorem of Johnson and Lindenstrauss, Random Structures & Algorithms 22, 60–65, 2003. https://doi.org/10.1002/rsa.10073
• M. Rudelson, R. Vershynin, The Littlewood–Offord problem and invertibility of random matrices, Advances in Mathematics 218, 600–633, 2008. https://arxiv.org/abs/math/0703503
• Y. Chen, J. Dongarra, Condition numbers of Gaussian random matrices, SIAM J. Matrix Anal. Appl. 27(3), 603–620, 2005. https://doi.org/10.1137/040616413
• A. Edelman, Eigenvalues and condition numbers of random matrices, SIAM J. Matrix Anal. Appl. 9(4), 543–560, 1988. https://doi.org/10.1137/0609045
• Y. Gordon, Some inequalities for Gaussian processes and applications, Israel J. Math. 50, 265–289, 1985. https://doi.org/10.1007/BF02759761
• K. R. Davidson, S. J. Szarek, Local operator theory, random matrices and Banach spaces, Handbook of the Geometry of Banach Spaces, Vol. 1, 317–366, 2001. https://doi.org/10.1016/S1874-5849(01)80010-3
• M. Ledoux, The Concentration of Measure Phenomenon, AMS Mathematical Surveys and Monographs 89, 2001. https://doi.org/10.1090/surv/089

Shared definitions: the items of this mission import Definitions.Def_GaussianMatrix_basic, a draft of Mission I; submit that proposal first so the import resolves at compile time.

Series: I. Rotation Invariance and Moment Identities · II. Inverse Moments of the Pseudoinverse · III. Extreme Singular Values and Concentration. All three missions share the definition bundle GaussianMatrix_basic, published with Mission I, and the namespace GaussianMatrix.

Motivation

Randomized algorithms for low-rank approximation — the randomized SVD of Halko, Martinsson and Tropp (HMT 2011), generalized Nyström, randomized subspace iteration, one-pass sketching (Tropp–Webber 2023) — multiply the input matrix by a random test matrix $\Omega$ with independent standard normal entries. Every error bound for these methods rests on a small number of facts about such standard Gaussian matrices: they are invariant under rotations, disjoint blocks of a rotated copy are independent, and the second and fourth moments of a linearly transformed Gaussian matrix are given by explicit formulas. These facts are classical, but none of them is available in Mathlib, so a Lean development of randomized numerical linear algebra currently has to carry every one of them as an explicit hypothesis of its main theorems. This mission removes those hypotheses.

Setting

For $p,m\in\mathbb N$, a standard Gaussian matrix of size $p\times m$ is a random matrix $\Gamma\in\mathbb R^{p\times m}$ whose $pm$ entries are independent $\mathcal N(0,1)$ random variables. Its law, written $\gamma_{p,m}$, is the product measure $\bigotimes_{i<p}\bigotimes_{j<m}\mathcal N(0,1)$ on real $p\times m$ arrays.

For a real matrix $A$, $\|A\|_F^2=\sum_{i,j}A_{ij}^2$ is the squared Frobenius norm; $\operatorname{rank}A$ is the dimension of its column space. A matrix $V_1\in\mathbb R^{n\times k}$ has orthonormal columns when $V_1^{\mathsf T}V_1=I_k$; a square such matrix is orthogonal. For $t\ge1$, a matrix $S\in\mathbb R^{n\times t}$ with independent $\mathcal N(0,1/t)$ entries is written $S=t^{-1/2}\Omega$ with $\Omega\sim\gamma_{n,t}$, so that $SS^{\mathsf T}=\tfrac1t\Omega\Omega^{\mathsf T}$ is an unbiased sketch of the identity.

Formalization targets

Goal — Gaussian approximate matrix multiplication

For fixed $M\in\mathbb R^{a\times n}$, $N\in\mathbb R^{b\times n}$ and $S\in\mathbb R^{n\times t}$ with independent $\mathcal N(0,1/t)$ entries,

$$\mathbb E\,\bigl\|MSS^{\mathsf T}N^{\mathsf T}-MN^{\mathsf T}\bigr\|_F^2=\frac{\|M\|_F^2\|N\|_F^2+\|MN^{\mathsf T}\|_F^2}{t}\le\frac2t\|M\|_F^2\|N\|_F^2 .$$

This is the exact mean-squared error of the sketched product $MSS^{\mathsf T}N^{\mathsf T}$ as an estimator of $MN^{\mathsf T}$, a direct consequence of Isserlis' pairing formula for Gaussian fourth moments. It is the goal because it is the one fourth-moment computation the series needs and because every other identity in the mission is a step towards it.

Milestones

• Rotational invariance. $U\Gamma V\sim\gamma_{p,m}$ for orthogonal $U\in\mathbb R^{p\times p}$, $V\in\mathbb R^{m\times m}$ (HMT §10.1).
• Blocks of a rotated Gaussian matrix. For $V_1\in\mathbb R^{n\times k}$ with orthonormal columns, $V_1^{\mathsf T}\Omega\sim\gamma_{k,t}$; for orthonormal $V_1,V_2$ with $V_1^{\mathsf T}V_2=0$, the blocks $V_1^{\mathsf T}\Omega$ and $V_2^{\mathsf T}\Omega$ are independent (HMT, proof of Theorem 10.5).
• Second moment. $\mathbb E\|S\Gamma T\|_F^2=\|S\|_F^2\|T\|_F^2$ (HMT Proposition 10.1 (10.1); Tropp–Webber Lemma B.1).
• Fourth moment. $\mathbb E\|S\Gamma T\|_F^4=(\|S\|_F^2\|T\|_F^2)^2+2\|S^{\mathsf T}S\|_F^2\|TT^{\mathsf T}\|_F^2$, the Isserlis formula for the Gaussian quadratic form $\operatorname{tr}(\Gamma^{\mathsf T}S^{\mathsf T}S\,\Gamma\,TT^{\mathsf T})$; its special case $\mathbb E\|\Omega\|_F^4=(kt)^2+2kt$ is the chi-square second moment (HMT Proposition A.8).
• Full rank. $\operatorname{rank}\Gamma=\min\{p,m\}$ almost surely (HMT, proof of Theorem 10.5 and Proposition A.5).

Significance

The results. Rotational invariance and block independence are the mechanism behind every conditioning argument in randomized linear algebra: writing the right singular vectors of the input as $[V_1\ V_2]$, the sketch splits into independent standard Gaussian blocks $\Omega_1=V_1^{\mathsf T}\Omega$ and $\Omega_2=V_2^{\mathsf T}\Omega$, and one may condition on $\Omega_1$ while treating $\Omega_2$ as fresh. The second-moment identity then turns the deterministic range-finder bound $\|\Sigma_2\Omega_2\Omega_1^{\dagger}\|_F$ into $\|\Sigma_2\|_F^2\,\mathbb E\|\Omega_1^{\dagger}\|_F^2$, and the approximate-multiplication identity gives the one-pass error terms. Almost-sure full rank is what makes $\Omega_1^{\dagger}$ a genuine pseudoinverse.

Formalizing them. Mathlib has the one-dimensional Gaussian measure, product measures, independence of coordinate projections, and (recently) the standard Gaussian measure on a finite-dimensional inner product space together with its invariance under linear isometries. None of the matrix-level statements above is in Mathlib. Once proved here they are reusable by every later mission in the series and by any formalization of randomized low-rank approximation, whose expectation-level inputs (e.g. $\mathbb E\|\Sigma_2\Omega_2\Omega_1^{\dagger}\|_F^2=\frac{k}{t-k-1}\|\Sigma_2\|_F^2$) follow from these milestones together with the inverse moments of Mission II.

Difficulty

Rotational invariance is not a computation with densities on $\mathbb R^{pm}$ that Mathlib supports out of the box: the clean route is through characteristic functions or through the basis-independence of the standard Gaussian on an inner product space, transported along the linear isometry $\Gamma\mapsto U\Gamma V$ of $\mathbb R^{p\times m}$ with the Frobenius inner product. Independence of the blocks requires identifying $[V_1\ V_2]^{\mathsf T}\Omega$ with a product measure and then projecting onto disjoint coordinate sets, with a completion of $[V_1\ V_2]$ to an orthogonal matrix when $k+r<n$. The moment identities are finite computations, but the fourth-moment identity needs the Gaussian moments $\mathbb E g^4=3$, $\mathbb E g^2=1$ organized through Isserlis' pairing formula over a four-fold sum of indices; the obvious entrywise expansion produces $O((pm)^4)$ terms that must be collapsed symbolically. Almost-sure full rank is the statement that the zero set of a nonzero polynomial ($\det$ of a maximal square submatrix) is Lebesgue-null, which must be established for the product Gaussian measure.

Formalization scope

The sample space is the function type $\mathbb R^{p}\to\mathbb R^{m}\to\mathbb R$ with the product $\sigma$-algebra; a sample is viewed as a matrix through the identity map $\mathrm{of}$. Expectations are real Bochner integrals: for an identity with a nonzero right-hand side this forces integrability, and the integrands here are polynomials in Gaussian entries. Laws are stated as equalities of pushforward measures; independence is Mathlib's IndepFun. Dimensions may be $0$; all statements remain true in that case. The hypothesis $V_1^{\mathsf T}V_1=I_k$ silently forces $k\le n$ (vacuous otherwise), which matches the source.

The definition bundle GaussianMatrix_basic supplies the law $\gamma_{p,m}$, $\|\cdot\|_F^2$, $\|\cdot\|_F$, the spectral norm (Mathlib's $\ell_2$ operator norm), the smallest singular value and the pseudoinverse $G^{\dagger}=G^{\mathsf T}(GG^{\mathsf T})^{-1}$; it is shared by the whole series. Welcome contributions: the Gaussian integration-by-parts and Isserlis lemmas for finite index sets, the identification of $\gamma_{p,m}$ with Mathlib's stdGaussian on the Frobenius inner product space, and the "polynomial zero sets are null" lemma for product measures.

Selected references

• N. Halko, P.-G. Martinsson, J. A. Tropp, Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions, SIAM Review 53(2), 217–288, 2011. https://arxiv.org/abs/0909.4061
• J. A. Tropp, R. J. Webber, Randomized algorithms for low-rank matrix approximation: design, analysis, and applications, 2023. https://arxiv.org/abs/2306.12418
• R. Vershynin, Introduction to the non-asymptotic analysis of random matrices, in Compressed Sensing: Theory and Applications, Cambridge University Press, 2012. https://arxiv.org/abs/1011.3027
• L. Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika 12, 134–139, 1918. https://doi.org/10.1093/biomet/12.1-2.134
• S. Dasgupta, A. Gupta, An elementary proof of a theorem of Johnson and Lindenstrauss, Random Structures & Algorithms 22, 60–65, 2003. https://doi.org/10.1002/rsa.10073
• M. Rudelson, R. Vershynin, The Littlewood–Offord problem and invertibility of random matrices, Advances in Mathematics 218, 600–633, 2008. https://arxiv.org/abs/math/0703503
• Y. Chen, J. Dongarra, Condition numbers of Gaussian random matrices, SIAM J. Matrix Anal. Appl. 27(3), 603–620, 2005. https://doi.org/10.1137/040616413
• A. Edelman, Eigenvalues and condition numbers of random matrices, SIAM J. Matrix Anal. Appl. 9(4), 543–560, 1988. https://doi.org/10.1137/0609045
• Y. Gordon, Some inequalities for Gaussian processes and applications, Israel J. Math. 50, 265–289, 1985. https://doi.org/10.1007/BF02759761
• K. R. Davidson, S. J. Szarek, Local operator theory, random matrices and Banach spaces, Handbook of the Geometry of Banach Spaces, Vol. 1, 317–366, 2001. https://doi.org/10.1016/S1874-5849(01)80010-3
• M. Ledoux, The Concentration of Measure Phenomenon, AMS Mathematical Surveys and Monographs 89, 2001. https://doi.org/10.1090/surv/089

Series: I. Rotation Invariance and Moment Identities · II. Inverse Moments of the Pseudoinverse · III. Extreme Singular Values and Concentration. All three missions share the definition bundle GaussianMatrix_basic, published with Mission I, and the namespace GaussianMatrix.

Problem

Solve a consistent system $Ax=b$, with $A\in\mathbb{C}^{m\times n}$ of full rank and $m\ge n\ge 1$. Kaczmarz's method (1937) projects the current iterate onto the solution hyperplane of one equation at a time. The cyclic version converges, but its rate depends on the order of the rows and has no clean bound in terms of a condition number.

Strohmer and Vershynin (2009) pick row $i$ at random with probability $\|a_i\|_2^2/\|A\|_F^2$ and prove

$$\mathbb{E}\,\|x_k-x\|_2^2 \le \bigl(1-\kappa(A)^{-2}\bigr)^k\,\|x_0-x\|_2^2,\qquad \kappa(A)=\frac{\|A\|_F}{\sigma_{\min}(A)}.$$

The rate does not depend on the number of equations $m$.

Setting

• Rows. $a_i\in\mathbb{C}^n$ is the conjugate of row $i$, so equation $i$ reads $\langle a_i,x\rangle=b_i$.
• Condition number. $\sigma_{\min}(A)=\inf_{\|z\|_2=1}\|Az\|_2$ and $\kappa(A)=\|A\|_F/\sigma_{\min}(A)$ (Demmel). It satisfies $\sqrt n\le\kappa(A)\le\sqrt n\,\|A\|_2/\sigma_{\min}(A)$.
• Algorithm 1. From any $x_0$, draw row $r$ independently with probability $p_r=\|a_r\|_2^2/\|A\|_F^2$ and set

$$x_{k+1}=x_k+\frac{b_r-\langle a_r,x_k\rangle}{\|a_r\|_2^2}\,a_r .$$

• Expectation. $\mathbb{E}\|x_k-x\|_2^2$ is a finite sum over the $m^k$ possible row sequences.

Targets

• Theorem 2 (goal): the bound above, for every $x_0$ and $k$.
• Theorem 3: some $x_0\ne x$ has $\mathbb{E}\|x_k-x\|_2^2\ge(1-2k/\kappa(A)^2)\|x_0-x\|_2^2$ for all $k\ge1$, so $\kappa(A)^{-2}$ is sharp up to a constant. Caveat: the source's proof only reaches the unsquared bound $\mathbb{E}\|x_k-x\|_2\ge\dots$ and then cites Jensen, which goes the wrong way. Its estimates give the squared bound with $4$ in place of $2$. Both forms are milestones.
• Sharpness (§3.2): Theorem 2 is an equality when $\kappa(A)=\sqrt n$.
• Iteration count (§2.1): $k\ge 2\log\varepsilon/\log(1-\kappa(A)^{-2})$ steps give $\mathbb{E}\|x_k-x\|_2^2\le\varepsilon^2\|x_0-x\|_2^2$.

Proof idea

A single projection can barely reduce the error, when the error is almost orthogonal to the chosen row. On average it always does. If $Z=a_j/\|a_j\|_2$ is drawn with probability $p_j$, then

$$\mathbb{E}\,|\langle z,Z\rangle|^2=\frac{\|Az\|_2^2}{\|A\|_F^2}\ge\kappa(A)^{-2}\|z\|_2^2 .$$

Pythagoras for one projection then gives $\mathbb{E}\|x_{k+1}-x\|_2^2\le(1-\kappa(A)^{-2})\,\|x_k-x\|_2^2$; induct on $k$.

Formalization

• Vectors live in EuclideanSpace ℂ (Fin n) and matrices in Matrix (Fin m) (Fin n) ℂ. Mathlib's inner product is conjugate-linear in its first argument, which is why $a_i$ is a conjugated row.
• Full rank is stated as injectivity of $z\mapsto Az$.
• The expectation expErrSq is an explicit finite sum, so no measure theory is needed. The tower identity is a milestone.
• A zero row makes the step the identity (Lean's division by zero) and has probability $0$, so it is harmless.
• Mathlib has no Kaczmarz iteration, scaled condition number or $\sigma_{\min}$ lower bound; the mission builds them.

History

• 1937: Kaczmarz introduces the cyclic method and proves convergence, with no rate.
• 1970: Gordon, Bender and Herman rediscover it as ART for tomography.
• 2009: Strohmer and Vershynin give the first rate in terms of a condition number.
• 2010 onward: extensions to noisy systems, block methods, SGD and sketch-and-project (Needell; Needell–Tropp; Needell–Srebro–Ward; Gower–Richtárik).

References

• S. Kaczmarz, Angenäherte Auflösung von Systemen linearer Gleichungen, Bull. Int. Acad. Polon. Sci. Lett. A 35 (1937), 355–357.
• T. Strohmer and R. Vershynin, A randomized Kaczmarz algorithm with exponential convergence, J. Fourier Anal. Appl. 15 (2009), 262–278. arXiv:math/0702226
• J. Demmel, The probability that a numerical analysis problem is difficult, Math. Comp. 50 (1988), 449–480. DOI
• D. Needell, Randomized Kaczmarz solver for noisy linear systems, BIT Numer. Math. 50 (2010), 395–403. arXiv:0902.0958
• D. Needell, N. Srebro and R. Ward, Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm, Math. Program. 155 (2016), 549–573. arXiv:1310.5715
• R. M. Gower and P. Richtárik, Randomized iterative methods for linear systems, SIAM J. Matrix Anal. Appl. 36 (2015), 1660–1690. arXiv:1506.03296

Motivation

A standard way to prove that a real number $\alpha$ is irrational is to produce integer linear forms $b+a\alpha$ that are nonzero but arbitrarily small: if $\alpha=p/q$ were rational, then $bq+ap$ would be a nonzero integer of absolute value below $1$ once the form is smaller than $1/q$. This is the shape of every hypergeometric construction of linear forms in odd zeta values — Rivoal's proof that infinitely many $\zeta(2n+1)$ are irrational and Zudilin's proof that at least one of $\zeta(5),\zeta(7),\zeta(9),\zeta(11)$ is irrational both produce such forms and read off irrationality (of at least one member of a finite set) from a determinant condition.

The same reduction is useful in the other direction: it isolates exactly what a construction has to supply — small forms — from the arithmetic that consumes them. This mission formalizes that abstract layer: the criteria that turn small integer forms into irrationality, together with the positivity and quadrature lemmas used to certify that a form is nonzero.

The material is distilled from a research note on $\zeta(9)$ (Xu, 2026). That note does not prove the irrationality of $\zeta(9)$, and nothing in this mission depends on whether it can be: every statement below is a statement about real numbers, integer linear forms, real polynomials, and finite sums, with $\zeta(9)$ and every other specific constant removed.

Setting

All objects live over $\mathbb{R}$.

• An integer linear form in $x$ is a number $b+a\,x$ with $a,b\in\mathbb{Z}$; the pair $(b,a)$ is its coefficient vector. Two forms are independent when their coefficient vectors have nonzero cross determinant, $b_1a_2\neq b_2a_1$.
• Irrational x is Mathlib's predicate: $x\notin\mathbb{Q}$ as a real number.
• The moment-matching hypothesis for a linear functional $L$ on real polynomials, a five-point node vector $y$ and a weight vector $w$, is $L(X^m)=\sum_{j}w_j\,y_j^{\,m}$ for every $m\le 4$. A functional satisfying it is exact on a polynomial $p$ when $L\,p=\sum_j w_j\,p(y_j)$.
• A positive weight vector has $w_j>0$ for all $j$; a node vector is injective when $y$ is injective on Fin 5.
• Polynomial.taylor u₀ p is the Taylor expansion of $p$ about $u_0$; its coefficients are nonnegative when $($taylor u₀ p$).\mathrm{coeff}\ i\ge 0$ for every $i$.
• Matrix.mulVec M v is the usual matrix–vector product over Fin 5; $\sum'$ denotes tsum over a Summable family.

Formalization targets

Goal — the one-form criterion

$$ \bigl(\forall \varepsilon>0,\ \exists, b,a\in\mathbb{Z}:\ b+ax\neq 0\ \wedge\ |b+ax|<\varepsilon\bigr)\ \Longrightarrow\ \text{$x$ irrational.}

The goal is the weakest non-vacuous statement in the family: it assumes one form at a time and no rate. ### Stronger — the two-form criterion

\bigl(\forall \varepsilon>0,\ \exists, b_1a_1b_2a_2\in\mathbb{Z}:\ b_1a_2\neq b_2a_1\ \wedge\ |b_1+a_1x|<\varepsilon\ \wedge\ |b_2+a_2x|<\varepsilon\bigr)\ \Longrightarrow\ \text{$x$ irrational.} $$

Supporting targets

1. Moment-matching quadrature — matching the five moments $m\le 4$ implies exactness on every polynomial of degree at most $4$.
2. Weighted average is interior — with positive weights summing to $1$, a non-constant five-tuple has its weighted average strictly between its minimum and maximum.
3. Mediant is interior — the ratio $\sum w_ia_i\,/\,\sum w_ib_i$ with $w,b>0$ lies strictly between the extreme values of $a_j/b_j$.
4. Positive matrices — an entrywise positive $5\times5$ matrix sends every nonzero nonnegative vector to a strictly positive vector.
5. Taylor-sign kernel sum — nonnegative Taylor coefficients at a lower bound of a sequence, positive summable weights, and one positive sample force a strictly positive weighted sum.
6. Five-sample nonvanishing — under moment matching with positive weights and injective nodes, a nonzero polynomial of degree $\le 4$ whose five sampled values share a sign has $L\,p\neq 0$.

Targets 1–6 correspond to the mission's milestones; the goal and the two-form criterion close the mission.

Significance

The results. The two criteria are the exact statements that a linear-form construction has to feed, and they are what turns "small forms exist" into irrationality without any analytic input. The supporting lemmas are the standard certificates used to show a form is nonzero — which is the other half of the argument, and the half that finite checks can actually settle.

Formalizing them. All eight statements are elementary and already have informal proofs; each also has a locally compiled Lean proof (lake env lean, exit 0, no sorry) against Lean 4.33.1 and Mathlib revision 0df444a3, held by the mission captain and published in the companion repository. What this mission adds is platform verification plus reusable infrastructure: the moment-matching quadrature lemma, the weighted-average and mediant inequalities, and the positivity lemmas are stated in a form that transfers to any setting where five-point data is certified by moments. Alternative proofs, generalizations to $n$-point quadrature, and sharper variants are welcome contributions.

Difficulty

The integrality step, not the estimate. In the one-form criterion the obvious move — take $\varepsilon=1/|q|$ — leaves the real inequality $|b+ax|<1/|q|$, which says nothing until the form is rewritten as $(bq+ap)/q$ with $bq+ap\in\mathbb{Z}$; only then does $|\,\cdot\,|<1$ force vanishing and contradict nonzeroness. Writing that rewrite in Lean means carrying the cast from $\mathbb{Z}$ through field_simp and back through exact_mod_cast, which is where naive attempts break.

Moment matching needs a degree bound, not interpolation. The quadrature lemma is not "five values determine a degree-$4$ polynomial": the hypothesis is about the functional $L$ on the five monomials, and the proof must expand an arbitrary $p$ in the monomial basis (as_sum_range_C_mul_X_pow' with natDegree < 5) and commute two finite sums.

Sign conditions are load-bearing. In target 6, the shared-sign hypothesis is what turns a vanishing weighted sum into vanishing samples; the root-counting step then needs injective nodes and positive weights. Dropping either silently makes the statement false, and both are easy to forget.

Formalization scope

Everything is over $\mathbb{R}$; no complex numbers appear. The quadrature statements are fixed at five nodes (Fin 5) and degree $\le 4$, as in the source note; the functional $L$ is a Polynomial ℝ →ₗ[ℝ] ℝ, not a measure. natDegree (not degree) is the degree notion. The infinite sum in target 5 is tsum with an explicit Summable hypothesis. Matrices are Matrix (Fin 5) (Fin 5) ℝ with mulVec; irrationality is Mathlib's Irrational.

Ruled out: a quadrature statement in which the weights are unconstrained by positivity but the conclusion is strengthened to a lower bound — target 1 assumes only moment matching, and any strengthening must add hypotheses rather than reinterpret the existing ones. A "criterion" whose hypothesis is vacuous for every real $x$ is likewise out of scope: both criteria are satisfiable hypotheses, not vacuous ones.

Infrastructure needed: the polynomial expansion and evaluation lemmas (as_sum_range, eval_eq_sum_range'), Finset sum rearrangement, Matrix.mulVec, Summable.tsum_lt_tsum_of_nonneg, and irrational_iff_ne_rational. The quadrature lemma, the mediant inequality, and the positivity lemmas are reusable beyond this mission.

Selected references

• Y. Xu, Research Notes on ζ(9): Constructions, Computations, and Open Problems, v0.1, Zenodo, 2026. https://doi.org/10.5281/zenodo.22951155
• W. Zudilin, Arithmetic of linear forms involving odd zeta values, J. Théor. Nombres Bordeaux 16:1 (2004), 251–291. https://arxiv.org/abs/math/0206176
• T. Rivoal, La fonction zêta de Riemann prend une infinité de valeurs irrationnelles aux entiers impairs, C. R. Acad. Sci. Paris Sér. I Math. 331 (2000), 267–270. https://arxiv.org/abs/math/0008051

Motivation

Finding the eigenvalues of a large matrix or linear operator by computing its characteristic polynomial is numerically unworkable: the roots of a degree-$n$ polynomial are exponentially sensitive to small coefficient perturbations, and no closed-form root formula exists once $n\ge5$. The power method avoids the polynomial entirely. Introduced in essentially its modern form by Müntz (1913) and von Mises and Pollaczek-Geiringer (1929), and analyzed rigorously alongside its shifted and inverse variants throughout the mid-20th century (Wilkinson, The Algebraic Eigenvalue Problem, 1965), it remains, in the guise of one power iteration per step, the engine inside PageRank, spectral clustering, and the Lanczos/Arnoldi methods used to find eigenpairs of matrices too large to diagonalize directly.

Setting

Let $E$ be a finite-dimensional inner product space over $\mathbb{k}\in\{\mathbb{R},\mathbb{C}\}$, with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$, and let $T:E\to E$ be a self-adjoint (symmetric) linear operator: $\langle Tx,y\rangle=\langle x,Ty\rangle$ for all $x,y\in E$. The spectral theorem for finite-dimensional self-adjoint operators gives an orthonormal basis $e_0,\dots,e_{n-1}$ of $E$ ($n=\dim E$) consisting of eigenvectors of $T$, with real eigenvalues $\lambda_0,\dots,\lambda_{n-1}$ satisfying $Te_i=\lambda_ie_i$.

Call $\lambda_{i_0}$ dominant if $|\lambda_j|<|\lambda_{i_0}|$ for every $j\neq i_0$ — it is then the unique eigenvalue of largest magnitude. Given a starting vector $x_0\in E$ with coordinates $x_0=\sum_ic_ie_i$ in the eigenbasis, the power iterates are $T^kx_0$ for $k=0,1,2,\dots$, and the Rayleigh quotient of $T$ at a nonzero vector $x$ is

$$R_T(x)=\frac{\operatorname{Re}\langle x,Tx\rangle}{\|x\|^2},$$

which recovers $\lambda_i$ exactly when $x$ is the eigenvector $e_i$.

Formalization targets

Iterate expansion

$$T^kx_0=\sum_i\bigl(c_i\lambda_i^{\,k}\bigr)\,e_i$$

Rewriting the $k$-th power iterate in the eigenbasis: applying $T$ $k$ times raises each coordinate's eigenvalue factor to the $k$-th power, since $T$ acts diagonally on the eigenbasis. This is the algebraic core the rest of the argument rescales and takes limits of.

Rescaled convergence

$$\lambda_{i_0}^{-k}\,T^kx_0\;\longrightarrow\;c_{i_0}\,e_{i_0}\quad(k\to\infty)$$

Given a dominant eigenvalue $\lambda_{i_0}\neq0$ and $c_{i_0}\neq0$, dividing the expansion above by $\lambda_{i_0}^k$ leaves the $i_0$-th term fixed at $c_{i_0}e_{i_0}$ while every other term is multiplied by $(\lambda_j/\lambda_{i_0})^k\to0$, since $|\lambda_j/\lambda_{i_0}|<1$ for $j\neq i_0$. This is the precise sense in which the power iterates "align" with the dominant eigenvector.

Goal — Rayleigh quotient convergence

$$R_T\bigl(T^kx_0\bigr)\;\longrightarrow\;\lambda_{i_0}\quad(k\to\infty)$$

The practical output of the power method: the Rayleigh quotient of the (unrescaled) iterates converges to the dominant eigenvalue itself, giving a numerically computable estimator that needs no knowledge of $\lambda_{i_0}$ in advance. This is the weakest statement that captures "the power method converges to the dominant eigenvalue" without hard-coding a convergence rate, so it is the mission's goal.

Significance

The power method is the template every practical large-scale eigenvalue algorithm departs from: shifted inverse iteration, Rayleigh quotient iteration (with locally cubic convergence), the QR algorithm, and Krylov subspace methods (Lanczos, Arnoldi) all begin from the same diagonal-power argument formalized here, then add a trick — a shift, a change of subspace, an orthogonalization step — to accelerate or extend it. The result itself is classical and completely settled mathematically; there is no open question in the convergence theory of the basic power method under the dominant-eigenvalue hypothesis used here. What this mission contributes is a machine-checked version of that classical argument built directly on Mathlib's existing finite-dimensional spectral theorem (LinearMap.IsSymmetric.eigenvalues/eigenvectorBasis) — as of this writing, Mathlib's InnerProductSpace/Spectrum.lean and Rayleigh.lean files contain the spectral decomposition itself, and a Rayleigh quotient for ContinuousLinearMap, but not this convergence statement.

Difficulty

The obvious first attempt is to bound $\|T^kx_0-\lambda_{i_0}^kc_{i_0}e_{i_0}\|$ by a naive sum of norms and take limits termwise; this works for the rescaled sequence (Milestone 2) but does not by itself give the Rayleigh-quotient limit, because $R_T$ is invariant only under nonzero scalar rescaling, not under limits taken carelessly — one has to first establish that the limit vector $c_{i_0}e_{i_0}$ is nonzero (using $c_{i_0}\neq0$), then invoke continuity of $R_T$ away from $0$ to transport the Tendsto from the rescaled sequence to $R_T(T^kx_0)=R_T(\lambda_{i_0}^{-k}T^kx_0)$. Getting the degenerate case $n=1$ right is the other trap: with only one eigenvalue, the dominance hypothesis is vacuous, and if that eigenvalue is allowed to be $0$ the rescaling $\lambda_{i_0}^{-k}$ divides by zero and the rescaled-convergence statement becomes false — the formalization must therefore assume $\lambda_{i_0}\neq0$ explicitly rather than deriving it from dominance alone.

Formalization scope

The mission works with a general RCLike 𝕜 field (real or complex $E$), a LinearMap.IsSymmetric operator on a FiniteDimensional inner product space, and Mathlib's own eigenvalues/ eigenvectorBasis (which already fixes the eigenbasis and a specific, decreasing-by-value ordering of eigenvalues — the formalization does not re-derive the spectral theorem). Dominance is stated by magnitude (|\lambda_j| < |\lambda_{i_0}|), not by position in Mathlib's ordering, since the dominant eigenvalue need not be the largest by value (it could be the most negative). The starting vector $x_0$ is arbitrary subject to $c_{i_0}\neq0$; no normalization ($\|x_0\|=1$) is imposed, since the Rayleigh quotient and the rescaled limit are both scale-invariant/ scale-equivariant. A trivializing formalization is ruled out explicitly: without both $\lambda_{i_0}\neq0$ and $c_{i_0}\neq0$, the $n=1$, $T=0$ counterexample above makes the rescaled-convergence statement false, so these are load-bearing hypotheses, not decoration. Contributions on the two milestones (the algebraic iterate expansion, and the rescaled-limit argument) are especially welcome, since they are reusable building blocks for any future mission on shifted/inverse power iteration or Rayleigh quotient iteration.

Selected references

• R. von Mises and H. Pollaczek-Geiringer, Praktische Verfahren der Gleichungsauflösung, ZAMM, 1929.
• J. H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University Press, 1965.
• L. N. Trefethen and D. Bau III, Numerical Linear Algebra, SIAM, 1997 (Lecture 27: the power method).