Numerical Analysis

Motivation

Large scale unconstrained optimization often makes an exact Newton step expensive. At each iterate, Newton's method asks for a solution of a linear system involving the Hessian, and the cost of solving that system can dominate the cost of evaluating the objective. Dembo and Steihaug's 1983 truncated Newton method uses conjugate gradients to stop the linear solve early, either when the residual is small enough or when the current direction has insufficient positive curvature. The resulting direction is paired with a line search. The question for this mission is whether those early exits still give a globally reliable optimization method under the paper's assumptions. The answer, in Theorem 2.1, is that every resulting run has gradient norm tending to zero, and a positive definite Hessian at one limit point makes the whole sequence converge to that point.

Setting

The objective is a twice continuously differentiable function $f:\mathbb R^n\to\mathbb R$. Write $g(x)=\nabla f(x)$ for its gradient and $H(x)=Dg(x)$ for its Hessian. For every initial point $x_0$, the paper assumes that the sublevel set $L(x_0)=\{x:f(x)\le f(x_0)\}$ is bounded. The formalization uses the Euclidean norm and its inner product, so the operator norm of $H$ is also Euclidean.

At a major iterate $x_k$, the minor conjugate-gradient iteration starts at $p_0=0$ with residual and direction $r_0=d_0=-g(x_k)$. It maintains a scalar $\delta_i$, initialized to $\|r_0\|^2$, alongside the usual conjugate-gradient state. Before taking its next conjugate-gradient step, it checks whether $\langle d_i,H(x_k)d_i\rangle\le\varepsilon\delta_i$. If so, it returns $d_0=-g(x_k)$ when $i=0$ and the current approximation $p_i$ otherwise. If this first test fails, it takes a conjugate-gradient step and checks whether $\|r_{i+1}\|\le\eta_k\|g(x_k)\|$. If so, it returns $p_{i+1}$. A TNCG direction is the direction returned by the first exit that occurs. The forcing value $\eta_k$ sets the residual tolerance.

The major iteration chooses a positive step length $\lambda_k$. With $p_k$ the minor iteration's output and $0<\alpha<1/2$, $\alpha<\beta<1$, the line search requires both sufficient decrease and a curvature condition:

$$f(x_k+\lambda_k p_k)\le f(x_k)+\alpha\lambda_k\langle g(x_k),p_k\rangle, \qquad \langle g(x_k+\lambda_kp_k),p_k\rangle\ge\beta\langle g(x_k),p_k\rangle.$$

The update is $x_{k+1}=x_k+\lambda_kp_k$. A run is a sequence satisfying this rule from its specified starting point. The paper also states an alternative function-value curvature condition; this mission defines it because the local unit-step theorem discusses both choices, although the main run uses the displayed gradient condition.

Formalization targets

Global convergence

For every starting point there is a TNCG run, and every run satisfies

$$\|g(x_k)\|\longrightarrow0.$$

If $x^*$ is a limit point of that run and $H(x^*)$ is positive definite, then

$$x_k\longrightarrow x^*.$$

The existence clause is the formal meaning of the paper's assertion that the iterates are well defined. The six milestones follow the paper's appendix: the conjugate-gradient relations of Theorem A.1; the uniform descent and direction bounds of Lemma A.2; the line-search existence and directional-derivative limit used in the proof of Theorem A.3; Theorem A.3 itself; and eventual unit-step admissibility in Theorem A.4. The companion Theorem 2.2 addresses the eventual residual exit and unit-step behavior when a run converges to a point with positive definite Hessian. These results are all stated in the original paper.

Significance

The target shows that the computational shortcut in the inner linear solve preserves a global first-order guarantee: no TNCG run can maintain a gradient norm bounded away from zero. It also identifies when the sequence has one limit, rather than merely having stationary accumulation points. The local companion result describes the regime where the truncated method can take full steps and the residual stopping test governs the minor iteration.

The mathematical results are proved in the 1983 paper, while this mission's Lean statements are open proof targets. A complete formalization would connect the finite-dimensional conjugate-gradient recurrence, the line-search argument on bounded sublevel sets, and the convergence statement in one machine-checked development. The conjugate-gradient identities and line-search facts could be reused in other optimization algorithms; the exact first-exit model and the global conclusion are specific to this method.

Difficulty

The minor iteration is not an ordinary positive definite conjugate-gradient solve. The Hessian may be indefinite away from a solution, and the iteration must stop at insufficient curvature before using a later residual state. A convergence argument that assumes positive definiteness at every iterate would therefore miss the case the algorithm was designed to handle. The line search must work for the returned direction, including the first iteration's steepest-descent fallback. The limit-point conclusion adds a separate issue: a vanishing gradient norm alone does not force an entire sequence to converge to a chosen limit point.

Formalization scope

Lean represents $\mathbb R^n$ as EuclideanSpace ℝ (Fin n), uses its Euclidean norm and inner product, and defines $H(x)$ as the derivative of the gradient. The objective is $C^2$ and every sublevel set is bounded, as specified at the start of the paper. The TNCG recurrence is total, but an admissible direction must come from its first actual exit, with the curvature test preceding the residual test. At $i=0$, the curvature exit returns $-g$, not the zero initial approximation. This prevents a false stationary run. The residual test is written as a product inequality, equivalent to the paper's quotient when the residual test is reached.

Step lengths are positive. A run is extended through stationary points by zero directions, giving an infinite sequence on which a limit can be stated. Theorem 2.1 imposes no extra restriction on the forcing sequence and no upper bound on step lengths. Theorem A.1 corrects the printed Krylov-space index from $H^{k-1}g$ to $H^kg$. Theorem A.4 adds $g(x^*)=0$, needed for the conclusion and true in the paper's applications. Theorem 2.2 requires nonnegative forcing values and guards its residual-exit conclusion at stationary iterates. These corrections are recorded item by item in the moderation notes. Contributions to the conjugate-gradient identities, bounded-level-set line search, and isolated-limit-point argument are all within scope.

Selected references

• R. S. Dembo and T. Steihaug, Truncated-Newton algorithms for large-scale unconstrained optimization, Mathematical Programming 26 (1983), 190–212. DOI: 10.1007/BF02592055.

Motivation

Truncated-Newton methods minimize a smooth function $f:\mathbb{R}^n\to\mathbb{R}$ by taking, at each iterate $x_k$, a search direction obtained by solving the Newton equation $H(x_k)p=-g(x_k)$ only approximately, with the linear conjugate gradient (CG) method. Here $g$ is the gradient and $H$ the Hessian of $f$. The approach was analysed by Dembo and Steihaug (Math. Programming 26, 1983), building on the inexact Newton framework of Dembo, Eisenstat and Steihaug (SIAM J. Numer. Anal. 19, 1982), and Newton–CG with curvature tests remains a standard large-scale method (Nocedal and Wright, Numerical Optimization, Ch. 7).

Away from a minimizer the Hessian need not be positive definite, and CG on an indefinite system can produce directions along which $f$ curves downward. The Dembo–Steihaug method stops CG when it meets such a direction. A practical question follows: started at or near a stationary point that is not a local minimizer, will the method ever notice the nonpositive curvature? Theorem 2.4 of the paper answers it for exact CG: it does, except in a degenerate case described by the eigenspaces of $H$.

Setting

Fix a symmetric linear operator $H$ on $\mathbb{R}^n$, with the Euclidean inner product $u^{\mathsf T}v$ and norm $\|\cdot\|$, and a vector $g\in\mathbb{R}^n$. In the paper these are $H(x_k)$ and $g(x_k)$ at the current iterate.

The TNCG minor iteration (p. 194) runs CG on $Hp=-g$ from $p_0=0$. Its state at iteration $i$ is $(p_i,r_i,d_i,\delta_i)$:

$$p_0=0,\quad r_0=-g,\quad d_0=r_0,\quad \delta_0=r_0^{\mathsf T}r_0,$$

and, writing $q_i=Hd_i$,

$$\alpha_i=\frac{r_i^{\mathsf T}r_i}{d_i^{\mathsf T}q_i},\ \ p_{i+1}=p_i+\alpha_id_i,\ \ r_{i+1}=r_i-\alpha_iq_i,\ \ \beta_i=\frac{r_{i+1}^{\mathsf T}r_{i+1}}{r_i^{\mathsf T}r_i},\ \ d_{i+1}=r_{i+1}+\beta_id_i,\ \ \delta_{i+1}=r_{i+1}^{\mathsf T}r_{i+1}+\beta_i^2\delta_i .$$

At iteration $i$ the method first applies the curvature test (2.1), $d_i^{\mathsf T}Hd_i\le\varepsilon\,\delta_i$, and stops if it holds; otherwise it computes $r_{i+1}$ and applies the truncation test (2.2), $\|r_{i+1}\|/\|g\|\le\eta$, and stops if it holds. With $\varepsilon=\eta=0$, (2.1) reads $d_i^{\mathsf T}Hd_i\le0$ and (2.2) reads $r_{i+1}=0$. The method exits through (2.1) at $i$ when no test fired at any iteration $j<i$ and (2.1) holds at $i$.

For $\mu\in\mathbb{R}$ the eigenspace is $E(\mu,H)=\{v: Hv=\mu v\}$ (A.21). The span of vectors $v_1,\dots,v_m$ is written $[v_1,\dots,v_m]$; the Krylov space of $g$ is $[g,Hg,\dots,H^{m}g]$.

Formalization targets

Goal: Theorem 2.4 (p. 197)

Let $\eta_k=\varepsilon=0$. If $H$ is not positive definite, then either the minor iteration exits through (2.1) at some iteration $i$, so that

$$d_i^{\mathsf T}Hd_i\le 0 \quad\text{at an iteration } i \text{ the method reaches},$$

or

$$g^{\mathsf T}v=0\quad\text{for every } v \text{ with } Hv=\mu v,\ \mu\le0 .$$

Milestones

1. Theorem A.1, (A.2) and (A.5), at $\varepsilon=0$ (p. 206). If $d_i^{\mathsf T}Hd_i>0$ for $i=0,\dots,k$, then $d_i^{\mathsf T}Hd_j=0$ for $i\ne j\le k$, and $[d_0,\dots,d_k]=[g,Hg,\dots,H^kg]$.
2. Theorem A.5 (p. 210). If $g$ has a nonzero projection on each eigenspace of $H$, $k$ is the number of distinct eigenvalues, and vectors $d_0,\dots,d_{k-1}$ are $H$-conjugate with $(d_i,Hd_i)>0$ and span $[g,\dots,H^{k-1}g]$, then $H$ is positive definite.
3. The first sentence of the proof of Theorem 2.4 (p. 211). If $g$ is not orthogonal to any eigenspace of $H$, and exact CG passes the curvature test up to iteration $i$ and ends with $r_{i+1}=0$, then $H$ is positive definite.

Significance

The theorem is the theoretical justification, given in the paper on p. 197, for using the direction of negative curvature $d_i$ produced in case (ii) of the minor iteration: with exact arithmetic and no truncation, CG either produces a direction of nonpositive curvature or the gradient lies entirely in the span of eigenvectors with positive eigenvalues. The paper draws the consequence that if $x_k\to x^*$ and the method terminates through (2.2), then, outside this rare case, $x^*$ is a strong local minimizer. Theorem A.5 is a self-contained statement about conjugate families in Krylov spaces, usable in any analysis of CG on indefinite systems.

The result is proved in the paper (Appendix, pp. 206–211), with the CG relations of Theorem A.1 cited from Hestenes and Stiefel. No machine-checked proof of these statements is known to this mission. Formalizing them produces a Lean development of CG on a symmetric, possibly indefinite operator, with its conjugacy and Krylov-span relations, which classical treatments state only for positive definite systems.

Difficulty

The CG relations of Theorem A.1 are usually proved for positive definite $H$, where every denominator is positive. Here $H$ is indefinite, and the relations hold only as long as the curvature along the directions has stayed positive; the induction must carry that hypothesis and must also show the residuals do not vanish early. The passage from "CG ended with $r=0$" to a statement about all eigenspaces is the second obstacle: the CG run sees only the Krylov space of $g$, so one must relate that space to the eigenspaces that $g$ meets, using symmetry of $H$ and a count of distinct eigenvalues. A direct argument via "the directions span $\mathbb{R}^n$" fails, because when $g$ misses an eigenspace the Krylov space is a proper subspace.

Formalization scope

$\mathbb{R}^n$ is EuclideanSpace ℝ (Fin n); $H$ is a continuous linear map that is self-adjoint (IsSelfAdjoint H), as the paper's $H$ is a Hessian. Positive definiteness is $v^{\mathsf T}Hv>0$ for all $v\ne0$. The CG recursion is total and carries no stopping test; division by zero returns $0$ in Lean, so states after a would-be stop are junk, and the statements read the iteration only through first-exit predicates. The test (2.2) is encoded in the product form $\|r_{i+1}\|\le\eta\|g\|$.

Theorem 2.4 is stated for an arbitrary symmetric $H$ and vector $g$ rather than at an iterate of a run of the full algorithm. It uses nothing else about $f$ or the run, and every such pair occurs as the Hessian and gradient of an admissible objective. The paper's (A.5) prints $H^{k-1}g$ as the last Krylov vector; the correct $H^kg$ is used.

A trivializing formalization exists and is ruled out: the first alternative is not "$d_i^{\mathsf T}Hd_i\le0$ for some $i$" over the unstopped recursion, which holds automatically once a residual vanishes (the recursion then continues with $d=0$); it requires that iteration $i$ is reached with no earlier exit. The goal does not assume conjugacy, the Krylov span identity, or the hypotheses of Theorem A.5.

A complete development needs the spectral theorem for self-adjoint operators on finite-dimensional spaces (in Mathlib), finiteness of the set of eigenvalues, and the CG invariants. The CG lemmas are reusable for any Newton–CG or trust-region CG analysis. Proofs of the milestones, and of auxiliary CG invariants such as (A.3)–(A.4), are welcome.

Selected references

• R. S. Dembo and T. Steihaug, Truncated-Newton algorithms for large-scale unconstrained optimization, Mathematical Programming 26 (1983) 190–212. https://doi.org/10.1007/BF02592055
• R. S. Dembo, S. C. Eisenstat and T. Steihaug, Inexact Newton methods, SIAM Journal on Numerical Analysis 19 (1982) 400–408. https://doi.org/10.1137/0719025
• M. R. Hestenes and E. Stiefel, Methods of conjugate gradients for solving linear systems, Journal of Research of the National Bureau of Standards 49 (1952) 409–436. https://doi.org/10.6028/jres.049.044
• T. Steihaug, The conjugate gradient method and trust regions in large scale optimization, SIAM Journal on Numerical Analysis 20 (1983) 626–637. https://doi.org/10.1137/0720042
• J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed., Springer, 2006. https://doi.org/10.1007/978-0-387-40065-5

Motivation

Minimizing a smooth function subject to simple bounds on the variables, $l\le x\le u$, is one of the most common problems in nonlinear optimization: it arises on its own (nonnegativity constraints, box-constrained least squares, parameter estimation with physical ranges) and as the subproblem of augmented Lagrangian and penalty methods. The classical approaches are active-set methods, which follow the boundary of the box and must guess which bounds are tight, and projected-gradient methods, which project trial points onto the box.

Coleman and Li proposed an alternative that never touches the boundary: an interior-reflective Newton method. Every iterate is strictly feasible; a new affine scaling $D(x)$ measures how far each variable is from the bound its gradient points toward; and the line search follows a reflective path, which bounces off the faces of the box instead of stopping at them or being projected back. The method is the basis of the trust-region-reflective algorithm in MATLAB's Optimization Toolbox (fmincon, lsqnonlin), so its convergence theory underlies a widely used solver.

Timeline:

• 1992–1994: Coleman and Li introduce the interior-reflective approach and analyse its convergence (Math. Programming 67, 1994), proving first-order global convergence (Theorem 8) and local quadratic convergence (Theorem 13).
• 1996: Coleman and Li extend the approach to trust-region subproblems (SIAM J. Optim. 6, 1996), the version implemented in MATLAB.

This mission is the first of two drawn from the 1994 paper; the second treats the local quadratic rate.

Setting

The problem (1.1) is

$$\min_{x\in\mathbb R^n} f(x)\quad\text{subject to}\quad l\le x\le u,$$

with $l\in(\mathbb R\cup\{-\infty\})^n$, $u\in(\mathbb R\cup\{+\infty\})^n$, $l\le u$. The feasible box is $\mathcal F=\{x:l\le x\le u\}$ and its strict interior is $\operatorname{int}(\mathcal F)=\{x:l<x<u\}$. Write $g(x)=\nabla f(x)$ and $H(x)=\nabla^2 f(x)$, and $g_k=g(x_k)$, $H_k=H(x_k)$ along a sequence.

The standing assumption (p. 191): for the starting point $x_1$, the level set $\mathcal L=\{x\in\mathcal F: f(x)\le f(x_1)\}$ is compact, and $f$ is twice continuously differentiable on an open set $D\supseteq\mathcal F$.

The scaling vector $v(x)$ (Definition 2) has components $v_i=x_i-u_i$ if $g_i<0$ and $u_i<\infty$; $v_i=x_i-l_i$ if $g_i\ge0$ and $l_i>-\infty$; $v_i=-1$ if $g_i<0$ and $u_i=\infty$; $v_i=1$ if $g_i\ge0$ and $l_i=-\infty$. The scaling matrix is $D(x)=\operatorname{diag}(|v(x)|^{1/2})$, and $D_k=D(x_k)$. The system $D(x)^2g(x)=0$ is equivalent to the first-order optimality conditions of (1.1).

The reflective transformation $R$ (Fig. 6) folds $\mathbb R^n$ onto $\mathcal F$ coordinate by coordinate: a coordinate with one finite bound is reflected about it, and a coordinate with two finite bounds is folded periodically (a triangle wave between $l_i$ and $u_i$). The reflective path from $x$ in the direction $s$ is $p(\alpha)=R(x+\alpha s)-x$: it moves along $s$ and reverses the corresponding component of the direction whenever it meets a face. A breakpoint is a stepsize at which the path lies on the boundary.

The interior-reflective algorithm (Fig. 9), with $0<\sigma_l<\sigma_u<1$ and $\rho>0$, starts at $x_1\in\operatorname{int}(\mathcal F)$; at each step it chooses a descent direction $s_k$ and a stepsize $\alpha_k>0$ that is not a breakpoint and satisfies

$$f(x_{k+1})<f(x_k)+\sigma_l\big(\alpha_kg_k^{T}s_k+\tfrac12\alpha_k^2\min(s_k^{T}H_ks_k,0)\big)\qquad(3.4)$$

and either

$$f(x_{k+1})>f(x_k)+\sigma_u\big(\alpha_kg_k^{T}s_k+\tfrac12\alpha_k^2\min(s_k^{T}H_ks_k,0)\big)\qquad(3.5)$$

or $\alpha_k>\rho$; then $x_{k+1}=x_k+p_k(\alpha_k)$.

The directions $\{s_k\}$ are constraint-compatible (Definition 3) if $\{D_k^{-2}s_k\}$ is bounded, and consistent (Definition 4) if $s_k^{T}g_k\to0$ implies $D_kg_k\to0$.

Formalization targets

Goal: Theorem 8

If $\{x_k\}$ is generated by the algorithm of Fig. 9 and $\{s_k\}$ is both consistent and constraint-compatible, then

$$D_k^2g_k\to0\qquad\text{and}\qquad\alpha_k^2\min(s_k^{T}H_ks_k,0)\to0 .$$

No constants are involved; the statement is the paper's.

Milestones

1. The iterates stay in $\mathcal L$ and $f(x_k)$ strictly decreases (proof of Theorem 8, p. 207).
2. Theorem 3: constraint compatibility bounds the breakpoints $BR_k(j)=|v_{kj}|/|s_{kj}|$ away from zero.
3. Lemma 7, in its one-sided form: if $\alpha_k\to0$, then $f(x_{k+1})-f(x_k)\le\alpha_kg_k^{T}s_k+C\alpha_k^2$ eventually.
4. $\alpha_kg_k^{T}s_k\to0$ and $\alpha_k^2\min(s_k^{T}H_ks_k,0)\to0$ (proof of Theorem 8, p. 208).
5. $g_k^{T}s_k\to0$ for constraint-compatible directions (proof of Theorem 8, p. 208).

Companion results: Theorem 2 (the acceptable stepsize interval exists, so the algorithm is well defined), and Theorems 5 (1)–(2) and 6 (1)–(2) (the scaled steepest-descent directions $-D_k^2g_k$ and $-D_k^2\operatorname{sgn}(g_k)$ are constraint-compatible and consistent).

Significance

Theorem 8 is the global first-order convergence theorem of the interior-reflective approach. With Theorems 5 and 6 it shows that the simplest instance, scaled steepest descent along the reflective path, reaches first-order points from any strictly feasible start; the paper's later second-order and quadratic results for the trust-region variant are built on it. Constraint compatibility and consistency isolate exactly what a search direction must satisfy, so the theorem applies to any direction rule that can be checked against these two conditions.

The result has been proved since 1994 but, to our knowledge, not machine-checked. A formalization requires a precise model of the reflective path, which the paper describes by an algorithm (Fig. 4) and a coordinate formula (Fig. 6), and a careful treatment of points where the scaling degenerates. The formal statements also make explicit a point the printed text leaves loose: Lemma 7 is printed as an equality but only an inequality holds, and only the inequality is used.

Difficulty

The argument has the shape of a standard sufficient-decrease proof, but the step that is routine for straight-line searches fails here. The path is not a line: as $\alpha$ grows, components of the direction are reversed at breakpoints, so $f(x_k+p_k(\alpha))-f(x_k)$ is not $\alpha g_k^{T}s_k+O(\alpha^2)$ in general, since a reflected component changes the linear term by an amount of order $\alpha$. Lemma 7 must show that, for constraint-compatible directions, only components whose reflection decreases the linear term can be reflected at small steps, and the Taylor remainders must be controlled uniformly over a run, on a neighbourhood of the compact level set inside the open set where $f$ is $C^2$. The scaling also degenerates at the boundary, where $v$ vanishes; the iterates must be kept strictly interior for $D_k^{-2}$ to be meaningful.

Formalization scope

Points are EuclideanSpace ℝ (Fin n); the bounds are EReal vectors with $l_i\neq+\infty$, $u_i\neq-\infty$ and $l\le u$. The box and the level set are the published LewisTorczon.BoundPS.box and levelSet. The gradient is Mathlib's gradient, and $s^{T}H(x)s$ is the inner product of the derivative of the gradient map applied to $s$ with $s$. $D(x)$ is never formed as a matrix: $D^2g$, $Dg$, $D^{-2}s$ and $D^2\operatorname{sgn}(g)$ are written componentwise. The reflective path is $R(x+\alpha s)-x$ with $R$ the transformation of Fig. 6, which the paper describes as the same method as Fig. 4 in another notation (p. 199); "not a breakpoint" is the requirement that the new iterate be interior. Descent direction is footnote 3's: strict decrease along the path for all small $\alpha>0$. Sequences start at index 0, which is the paper's $x_1$. The standing compactness and smoothness assumption is carried as explicit hypotheses.

The goal does not assume $g_k^{T}s_k\to0$, $\alpha_k\to0$, or anything about breakpoints; a version of Theorem 8 that did, or one whose algorithm projected onto the box instead of reflecting, or whose run predicate no sequence can satisfy, would not be the paper's theorem. Lemma 7 keeps the paper's assumption $\alpha_k\to0$ and records the one-sided upper bound established by its proof. Theorem 3 uses the breakpoint distance in the direction of motion and applies its bound when that distance equals $|v_{kj}|/|s_{kj}|$.

A complete development needs Taylor's theorem with uniform remainder along piecewise linear paths in a box, and basic facts about the reflection $R$ (it maps into $\mathcal F$, is the identity on $\mathcal F$, and is continuous). These are reusable for any reflective or projected line search. Proofs of the milestones in any order are welcome, as are independent proofs of the companion results.

Selected references

• T. F. Coleman and Y. Li, On the convergence of interior-reflective Newton methods for nonlinear minimization subject to bounds, Mathematical Programming 67 (1994), 189–224. https://doi.org/10.1007/BF01582221
• T. F. Coleman and Y. Li, An interior trust region approach for nonlinear minimization subject to bounds, SIAM Journal on Optimization 6 (1996), 418–445. https://doi.org/10.1137/0806023
• R. M. Lewis and V. Torczon, Pattern search algorithms for bound constrained minimization, SIAM Journal on Optimization 9 (1999), 1082–1099 (source of the published box and level-set definitions reused here). https://doi.org/10.1137/S1052623496300507

Motivation

Smooth optimization with nonnegative variables appears when the variables are Lagrange multipliers for inequality constraints or when an objective incorporates an augmented Lagrangian or an exact penalty. Bertsekas studies how to retain a Newton-like convergence rate in this setting while using an iteration that projects a scaled gradient step onto the nonnegative orthant. His 1982 paper also discusses large optimal-control examples, where repeatedly solving a quadratic subproblem may be costly. These are the applications motivating the method, rather than assumptions of the theorem. Bertsekas, 1982.

The paper contrasts a partly diagonal scaling matrix with two other approaches: a fully diagonal projected gradient step, which generally has a linear rate, and a constrained Newton step defined by a quadratic program. The main theoretical question is whether the simpler projected iteration can still identify binding coordinates and converge superlinearly near the solution. The mission covers the nonnegative-orthant method in §2 and its stated convergence results. The paper's later extension to general linear constraints and its computational examples concern different objects. Bertsekas, 1982.

Setting

Fix a dimension $n\ge1$ and a continuously differentiable function $f:\mathbb R^n\to\mathbb R$. Problem (1) minimizes $f(x)$ over the nonnegative orthant $\mathbb R_+^n=\{x:x^i\ge0\text{ for every }i\}$. The positive part $[z]^+$ replaces each negative coordinate of $z$ by zero. A feasible $x$ is critical when every partial derivative $\partial_i f(x)$ is nonnegative and $\partial_i f(x)=0$ wherever $x^i>0$. This is the paper's componentwise first-order condition.

For a feasible iterate $x_k$, the binding set is $B(x_k)=\{i:x_k^i=0\}$. The paper selects a larger working set $I_k^+$ using a fixed $\varepsilon>0$ and the projected-gradient residual $w_k=\|x_k-[x_k-M\nabla f(x_k)]^+\|$, where $M$ is fixed, diagonal and positive definite. More precisely, $I_k^+$ contains the indices with $0\le x_k^i\le\min(\varepsilon,w_k)$ and $\partial_i f(x_k)>0$. The matrix $D_k$ is symmetric positive definite and has zero off-diagonal entries in the rows indexed by $I_k^+$. The projected arc and update are

$$p_k=D_k\nabla f(x_k),\qquad x_k(a)=[x_k-a p_k]^+,\qquad x_{k+1}=x_k(\beta^{m_k}).$$

Here $0<\beta<1$, and $m_k$ is the first nonnegative integer passing the paper's two-sum Armijo test (37), with $0<\sigma<1/2$. One sum uses the scaled gradient outside $I_k^+$; the other uses the actual projected displacement on $I_k^+$. These details define the algorithm whose rate is at issue. Bertsekas, 1982, pp. 228–229.

Formalization targets

The first targets establish the fixed-point and descent properties of the projected arc, positivity of the Armijo right-hand side, well-defined step selection, criticality of limit points, and local attraction with finite binding-set identification. Proposition 3 identifies both the working set and the actual binding set:

$$I_k^+=B(x_k)=B(x^*)\quad\text{for all sufficiently late }k.$$

The goal is Proposition 4. Let $f$ be convex and $C^2$, let $x^*$ be the unique minimizer on $\mathbb R_+^n$ satisfying Assumption (C), and suppose the Hessian quadratic form has uniform positive lower and finite upper bounds on the initial, unrestricted sublevel set. The scaling matrix is $D_k=H_k^{-1}$, where $H_k$ retains Hessian entries except for off-diagonal entries touching $I_k^+$. Then

$$x_k\longrightarrow x^*,\qquad \forall c>0\ \exists K\ \forall k\ge K:\ \|x_{k+1}-x^*\|\le c\|x_k-x^*\|.$$

If the Hessian is Lipschitz in a neighborhood of $x^*$, the same proposition states an at-least-quadratic error bound: $\|x_{k+1}-x^*\|\le C\|x_k-x^*\|^2$ eventually for some $C>0$. The separate final milestone records the paper's claim that the initial unit trial is eventually accepted. Bertsekas, 1982, pp. 233–236.

Significance

Proposition 4 says that this specific orthant-projected Newton iteration converges to the unique constrained optimum with a superlinear rate under its stated smoothness and curvature conditions. Proposition 3 supplies a distinct finite-identification statement: the working and binding coordinates eventually match those at the optimum. These facts specify the behavior of an algorithm that does not define its direction by solving a quadratic program at every step. The paper proves the preceding propositions and states Proposition 4 with its proof left to the reader, referring to its earlier discussion and standard unconstrained Newton results. Bertsekas, 1982, pp. 234–236.

Formalizing the claims would provide machine-checked statements and, once solved, proofs for the exact projected arc, two-part line search, matrix selection, identification result and rate conclusion. The local mission items are open proof obligations; their compilation checks the definitions and theorem types, not the mathematical claims. The geometry definitions can also support other orthant-constrained algorithms, while the enlarged working set and Armijo rule belong to this paper's method.

Difficulty

A positive definite matrix by itself does not make projection along $[x-aD\nabla f(x)]^+$ a descent move. An off-diagonal coupling can push coordinates against the boundary in a way that defeats the usual unconstrained descent calculation; the paper gives such a situation before Proposition 1. The partly diagonal condition is therefore substantive. A second obstacle is that the exact active set $I^+(x)$ can jump at a boundary point: iterates approaching that point from the interior need not have the same indexed rows. The enlarged $I_k^+$ and its dependence on the residual are central to the finite-identification and rate claims. Bertsekas, 1982, pp. 225–229.

Formalization scope

Lean represents $\mathbb R^n$ as EuclideanSpace ℝ (Fin n), with the Euclidean norm and zero-based coordinates. The dimension is positive. gradient and the derivative of gradient represent first and second derivatives; the paper's $M$ is diag μ with every $\mu^i>0$. The run predicate includes a feasible initial point and the first acceptable integer $m_k$, and the matrix choices are indexed by iteration. Propositions 1–3 use the paper's explicit admissibility condition. Proposition 4 constructs $D_k$ from $H_k$ and does not assume its invertibility, admissibility, convergence or eventual active-set equality.

Assumption (C) includes local $C^2$ smoothness, curvature bounds on directions zero at the binding coordinates, and strict complementarity. A local minimum is required to be feasible as well as locally minimal on the orthant. Limit points use subsequential convergence. Superlinearity uses a uniform eventual error inequality, which also covers an iterate that reaches the solution exactly; a ratio with a zero denominator would distort this case. The paper's Proposition 4 display mistakenly binds a direction to the level set while leaving the Hessian's point free. The formalization states the intended reading: every point in the unrestricted initial level set and every direction satisfy the Hessian bounds. The displayed level set has no $x\ge0$ restriction, so neither does the Lean hypothesis.

A definition that accepts any Armijo exponent or a goal that assumes the eventual identification or convergence conclusion would erase the paper's claim. Contributions can build the matrix and projected-arc lemmas, the convergence and identification proofs, and the final rate proof. The local inverse-Hessian observation before Proposition 4 is discussed in the notes but has no separate milestone until its full local setting can be captured without weakening it.

Selected references

• Dimitri P. Bertsekas, Projected Newton Methods for Optimization Problems with Simple Constraints, SIAM Journal on Control and Optimization 20(2), 221–246, 1982. DOI: 10.1137/0320018.

Motivation

Optimization over sets of matrices of fixed rank comes up in low-rank matrix completion, model reduction, and the low-rank approximation of solutions of large matrix equations. A standard approach treats the constraint set as a Riemannian manifold and runs gradient or Newton-type methods on it (Absil, Mahony & Sepulchre 2008). Each iteration takes a step in the tangent space and then has to return to the manifold. A map that does this to first order is called a retraction, and its cost is often what decides whether a manifold method is practical.

Absil and Malick (2012) study retractions defined by projection: step to $X+Z$ in the ambient space, then take the nearest point of the manifold. Their Proposition 3.2 shows that for any $C^k$ submanifold ($k\ge2$) this projective retraction is a retraction. Section 3.2 makes it explicit for the manifold of fixed-rank matrices: near any matrix of rank $r$, the nearest matrix of rank exactly $r$ is the truncated singular value decomposition (Proposition 3.3). Section 4.4 also gives a closed form for a second retraction on the same manifold, the orthographic retraction (Proposition 4.11).

Setting

Fix natural numbers $n$, $m$ and $r\ge1$. The space $\mathbb R^{n\times m}$ of real $n\times m$ matrices carries the Frobenius norm $\|X\|^2=\sum_{i,j}X_{ij}^2=\operatorname{trace}(X^\top X)$ (3.6). The fixed-rank set is

$$\mathcal R_r=\{X\in\mathbb R^{n\times m}:\ \operatorname{rank}(X)=r\},$$

a smooth submanifold of $\mathbb R^{n\times m}$. It is not closed: its closure is the set of matrices of rank at most $r$.

A singular value decomposition (3.5) of $X$ is a factorization $X=U\Sigma V^\top$ in which $U=[u_1,\dots,u_n]\in\mathbb R^{n\times n}$ and $V=[v_1,\dots,v_m]\in\mathbb R^{m\times m}$ are orthogonal and $\Sigma\in\mathbb R^{n\times m}$ is zero off its diagonal. The diagonal of $\Sigma$ holds the singular values of $X$ in nonincreasing order,

$$\sigma_1(X)\ge\sigma_2(X)\ge\cdots\ge\sigma_{\min\{n,m\}}(X)\ge0,$$

and $\sigma_i(X)=0$ for $i>\min\{n,m\}$. A matrix has rank $r$ exactly when $\sigma_r(X)>0=\sigma_{r+1}(X)$. The truncated SVD of rank $r$ is $\hat X=\sum_{i=1}^r\sigma_i(X)u_iv_i^\top$ (3.7).

For a set $Q$ and a point $X$, the projection $P_Q(X)$ is the set of nearest points: the $Y\in Q$ with $\|X-Y\|\le\|X-W\|$ for all $W\in Q$. For a set that is not closed, $P_Q(X)$ may be empty or contain several points.

Formalization targets

Goal: Proposition 3.3

Let $\bar X\in\mathcal R_r$. For every $X$ with $\|X-\bar X\|<\sigma_r(\bar X)/2$ and every singular value decomposition $X=U\Sigma V^\top$,

$$P_{\mathcal R_r}(X)=\Bigl\{\sum_{i=1}^r\sigma_i(X)\,u_iv_i^\top\Bigr\}.$$

The projection exists, is unique, and is the truncated SVD. The radius $\sigma_r(\bar X)/2$ and the strict inequality are those of the paper.

Milestones

The paper's proof goes through four claims, which are the milestones in attack order:

1. Weyl's bound (§3.2, proof of Proposition 3.3, citing Horn–Johnson 7.3.8): $|\sigma_i(\bar X)-\sigma_i(X)|\le\|X-\bar X\|$ for every $i\ge1$.
2. Eckart–Young (3.7): for every singular value decomposition of $X$, $\hat X$ is a nearest matrix to $X$ of rank at most $r$.
3. The gap (3.8): if $\operatorname{rank}\bar X=r$ and $\|X-\bar X\|<\sigma_r(\bar X)/2$, then $\sigma_{r+1}(X)<\sigma_r(\bar X)/2<\sigma_r(X)$.
4. Uniqueness under a gap (§3.2, proof of Proposition 3.3, citing Helmke–Moore): if $\sigma_r(X)>\sigma_{r+1}(X)$, then $\hat X$ is the only nearest matrix of rank at most $r$.

Further result: Proposition 4.11

Write $X=U\left[\begin{smallmatrix}\Sigma_0&0\\0&0\end{smallmatrix}\right]V^\top$ (4.7), with $\Sigma_0$ the positive diagonal of nonzero singular values, and a tangent vector $Z=U\left[\begin{smallmatrix}A&C\\B&0\end{smallmatrix}\right]V^\top$ (4.8). If $\Sigma_0+A$ is invertible, the orthographic retraction is

$$R(X,Z)=U\begin{bmatrix}\Sigma_0+A&C\\B&B(\Sigma_0+A)^{-1}C\end{bmatrix}V^\top .$$

Here $R(X,Z)$ is the nearest point to $X+Z$ of $\mathcal R_r\cap(X+Z+N_{\mathcal R_r}(X))$.

Significance

The result. Proposition 3.3 reduces the projective retraction on $\mathcal R_r$ to one truncated singular value decomposition. With Proposition 3.2 this gives an explicit, computable retraction for Riemannian optimization on fixed-rank matrices. It is used, for instance, in low-rank matrix completion algorithms. The point is local: $\mathcal R_r$ is not closed, so far from $\mathcal R_r$ the nearest matrix of rank exactly $r$ need not exist, and the proposition gives an explicit radius on which it does exist and is unique. Proposition 4.11 gives a second retraction that needs only products of matrices and one $r\times r$ inverse.

Formalizing it. All results here are proved, on paper. To the best of a search of Mathlib and the Prove2Me catalog, none is machine-checked. Mathlib has singular values of linear maps between finite-dimensional inner product spaces (LinearMap.singularValues), but no singular value decomposition in matrix form, no Weyl perturbation inequality for singular values, and no Eckart–Young theorem. A complete development would add these, and the Eckart–Young theorem with its uniqueness case is a standard result of numerical linear algebra in its own right.

Difficulty

The proof in the paper is short only because it cites three facts, and each of them is a real piece of matrix analysis.

• Weyl's bound needs the variational (min–max) description of singular values, which Mathlib does not have for singular values.
• Eckart–Young requires comparing $\|X-Y\|$ with the singular values of $X$ for every $Y$ of rank at most $r$, not only for those diagonal in the same bases. The obvious approach, writing $Y$ in the singular bases of $X$, fails because $Y$ need not be diagonal there.
• Uniqueness under the gap requires showing that every minimizer is diagonal in some singular bases of $X$ and then using the gap to fix its support. Without the gap uniqueness fails: for $X=I_2$ and $r=1$ every $uu^\top$ with $\|u\|=1$ is a nearest point.

A further subtlety: Proposition 3.3 must hold for any singular value decomposition of $X$. Singular vectors are not unique, so the proof has to show that the truncation does not depend on the choice under the gap (3.8).

Formalization scope

Matrices are Matrix (Fin n) (Fin m) ℝ with Mathlib's Frobenius norm (open scoped Matrix.Norms.Frobenius). The projection is the published platform predicate RandomGradFree.Nonsmooth.IsMetricProjection, and $P_{\mathcal R_r}(X)$ is the set {Y | IsMetricProjection (rankSet r) X Y}. The set equality with a singleton states existence, uniqueness and the formula together.

Conventions committed to:

• sv X i is $\sigma_i(X)$, 1-based as on the page. It is Mathlib's 0-based singularValues of Matrix.toEuclideanLin X at i - 1, so the index $0$ is meaningless and every statement uses indices $\ge1$.
• IsSVD X U S V is the predicate of (3.5). The singular value decomposition is a hypothesis of the theorems, never a choice made inside them, so the theorems hold for every singular value decomposition.
• truncSVD r U S V is $\sum_{i\le r}\Sigma_{ii}u_iv_i^\top$, written with the diagonal of S. That these entries are the singular values $\sigma_i(X)$ is a fact to be proved, not part of the definition.
• The hypothesis $r\ge1$ is added, because $\sigma_r$ needs it (and $\mathcal R_0=\{0\}$). It is the only hypothesis of Proposition 3.3 not printed on the page.
• In Proposition 4.11, matrices use the block index types Fin r ⊕ Fin p and Fin r ⊕ Fin q. The tangent vector and the normal space are taken in the form the page gives, and invertibility of $\Sigma_0+A$ is an explicit hypothesis. In the paper's proof it comes from "$Z$ in a neighborhood of the origin".

A formalization in which the radius is replaced by a smaller one, the projection is taken onto the matrices of rank at most $r$, or the singular value decomposition is chosen inside the statement is a different theorem and is ruled out.

Wanted contributions, all reusable beyond this mission:

• existence of a singular value decomposition in matrix form, and the identification of its diagonal with singularValues;
• Weyl's inequality for singular values;
• the Eckart–Young theorem and its uniqueness case;
• the Schur-complement rank formula for $2\times2$ block matrices, used for Proposition 4.11.

Selected references

• P.-A. Absil and J. Malick, Projection-like retractions on matrix manifolds, SIAM J. Optim. 22(1):135–158, 2012. https://doi.org/10.1137/100802529 (authors' version HAL hal-00651608v2: https://hal.science/hal-00651608v2)
• P.-A. Absil, R. Mahony and R. Sepulchre, Optimization Algorithms on Matrix Manifolds, Princeton University Press, 2008. https://press.princeton.edu/absil
• C. Eckart and G. Young, The approximation of one matrix by another of lower rank, Psychometrika 1(3):211–218, 1936. https://doi.org/10.1007/BF02288367
• R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press, 1985 (§7.3, §7.4). https://doi.org/10.1017/CBO9780511810817
• U. Helmke and J. B. Moore, Optimization and Dynamical Systems, Springer, 1994 (Ch. 5). https://doi.org/10.1007/978-1-4471-3467-1

Motivation

The r-algorithm of N. Z. Shor is a variable-metric method for minimizing functions that are not differentiable everywhere, such as maxima of smooth functions, penalty functions and Lagrangian dual functions of integer and decomposition problems. It combines a subgradient step with a space dilation along the difference of two successive almost-gradients, and the book reports it competitive with the best variable-metric and conjugate-gradient methods on "gully"-shaped test problems (Shor 1985, §3.6, pp. 68–77). Its convergence theory is limited. Section 3.7 of Shor's book gives a convergence proof only for an idealized version with exact directional minimization, on a class of piecewise smooth functions, and the book itself calls a weakening of the assumptions and rate estimates desirable (p. 85). This mission formalizes that section.

Timeline. Shor and Zhurbenko introduced the r-algorithm in 1971 (Kibernetika, no. 3, 51–59). Shor analysed the convergence of its exact-line-search version in 1975 (Kibernetika, no. 4, 48–53). Section 3.7 of the 1979 Russian monograph, translated in 1985, gives that analysis for piecewise smooth functions. The book states on p. 85 that weaker assumptions and rate estimates would be desirable.

Setting

Let $E_n$ be $n$-dimensional Euclidean space with inner product $(x, y)$, $n \ge 1$. For a unit vector $\xi$ and $\alpha > 0$ the operator of space dilation is $R_\alpha(\xi)x = x + (\alpha - 1)(x, \xi)\xi$.

Widths. For a compact convex set $W$ and a unit vector $\eta$, the width in direction $\eta$ is $d_\eta(W) = \max_{x \in W}(\eta, x) - \min_{x \in W}(\eta, x)$. The width is $d(W) = \min_{\|\eta\| = 1} d_\eta(W)$ and the diameter is $D(W) = \max_{\|\eta\| = 1} d_\eta(W)$. With $e_\rho(W) = \inf_{z \in W}|(\rho, z)|$ put $K_\rho(W) = d_\rho(W)/e_\rho(W)$ (or $+\infty$ when $e_\rho(W) = 0$) and $p(W) = \inf_{\|\rho\| = 1} K_\rho(W)$.

The class $K$. $E_n$ is partitioned into closed sets $\bar D_1, \dots, \bar D_m$, each the closure of its interior, whose interiors are disjoint and homeomorphic to an open ball or open halfspace; $f_i$ is continuously differentiable on an open set containing $\bar D_i$; and $f = f_i$ on $\bar D_i$. At a point $x$, the set of almost-gradients $G_f(x)$ is the set of gradients $\nabla f_i(x)$ of the pieces with $x \in \bar D_i$. For $\delta, \varepsilon > 0$, $\bar P_{\delta,\varepsilon}(x)$ is the closed convex hull of all $G_f(y)$, $\|y - x\| < \delta$, together with the $\varepsilon$-balls around the elements of $G_f(x)$.

The $r(\alpha)$-algorithm. Fix $\alpha > 1$, $\beta = 1/\alpha$. Start from $x_0$, $\tilde g_0 = 0$ and a nonsingular $B_0$. At iteration $k + 1$ choose $g_f(x_k) \in G_f(x_k)$ with $(B_k^* g_f(x_k), \tilde g_k) \le 0$; set $g_k^* = B_k^* g_f(x_k)$, $r_k = g_k^* - \tilde g_k$, $\xi_{k+1} = r_k/\|r_k\|$, $B_{k+1} = B_k R_\beta(\xi_{k+1})$, $\tilde g_{k+1} = R_\beta(\xi_{k+1}) g_k^*$, and $x_{k+1} = x_k - h_{k+1}B_{k+1}\tilde g_{k+1}$ with $h_{k+1} \ge 0$ such that $f$ does not increase along the step and some almost-gradient at $x_{k+1}$ makes a non-acute angle with $\tilde g_{k+1}$ in the transformed metric. The standing assumption is $f(x) \to +\infty$ as $\|x\| \to \infty$ (3.50).

Formalization targets

Goal: convergence to an isolated local minimum (Theorem 3.13)

Assume (3.50) and $\|x_{k+1} - x_k\| \to 0$ (3.52). If $x^*$ is an isolated local minimum, the component of $\{x : f(x^*) \le f(x) \le f(x_0)\}$ containing $x_0$ also contains $x^*$, and no other point $z$ of that component has linearly dependent $G_f(z)$, then

$$\lim_{k \to \infty} x_k = x^* .$$

Milestones

• Lemma 3.2 (p. 80): for $B = SO$ with minimum eigenvalue $\lambda(B)$ of $S$, $\lambda(B)d(W) \le d(BW) \le \lambda(B)D(W)$.
• Lemma 3.3 (p. 80): for $z_1, z_2 \in W$, $0 < \beta \le 1$ and $\gamma = \|z_1 - z_2\|/d(W) \ge 1$,

$$d\Big(R_\beta\big(\tfrac{z_1 - z_2}{\|z_1 - z_2\|}\big)W\Big) \ge \frac{d(W)}{\sqrt{1 + (1-\beta^2)/(\beta^2\gamma^2)}} .$$

• Theorem 3.11 (p. 82): for every $\sqrt[n]{\beta} < v < 1$, $\varepsilon, \delta > 0$ and $r \ge 1$ there is $\bar k > r$ with

$$p\big(\bar P_{\delta,\varepsilon}(x_{\bar k})\big) \ge \sqrt{\frac{v^2\sqrt[n]{\alpha^2} - 1}{\alpha^2 - 1}} .$$

• Theorem 3.12 (p. 84): the level set $\{f = f_\infty\}$, $f_\infty = \lim_k f(x_k)$, contains a point $x^*$ with $G_f(x^*)$ linearly dependent.

Significance

Theorem 3.13 is the only convergence theorem the book gives for the r-algorithm. Theorems 3.11 and 3.12 hold for any function of the class $K$, convex or not, and say that iterates of the exact version cannot stall at a point where the local almost-gradients are linearly independent: linear dependence of $G_f(x)$ is a generalized stationarity condition that includes $0 \in \operatorname{conv} G_f(x)$. Lemmas 3.2 and 3.3 are statements about widths of convex bodies under linear maps and single dilations, and apply to any analysis of space-dilation methods, including the ellipsoid method.

All four results and the goal are proved in the book, the goal as a corollary without a written proof. None of them is formalized. The formalization adds a machine-checked definition of the class $K$ and of the algorithm as a relation on sequences covering every admissible choice of almost-gradients and stepsizes, a verified proof of the book's argument, and a check of the corollary step, which the book leaves to the reader.

Difficulty

The obvious argument for descent methods is that the function value drops by a fixed amount at every step unless the gradient is small. That argument fails here. The r-algorithm can make null steps, with $x_{k+1} = x_k$ while $B_k$ and $\tilde g_k$ change, and its steps are measured in a metric that degenerates as the dilations accumulate ($\det B_k = \beta^k \det B_0$). A small step does not certify near-stationarity, and a large accumulated dilation does not certify progress.

The difficulty is to relate the geometry of the transformed almost-gradient sets $B_k^*\bar P_{\delta,\varepsilon}(x_k)$ to the decay forced on $B_k$. This needs width estimates for convex bodies under linear maps and single dilations, which Mathlib does not have. Theorem 3.12 also needs a uniform lower bound, near a compact level set, on the Gram determinants of the piece gradients, and the book's proof of it is only sketched. Theorem 3.13 is stated as a corollary with no proof at all. Deriving it requires showing that the iterates cannot leave the prescribed component of $\{f(x^*) \le f \le f(x_0)\}$, and that their accumulation points reduce to $x^*$.

Formalization scope

• $E_n$ is EuclideanSpace ℝ (Fin n) with $n \ge 1$; operators are continuous linear maps; $B_k^*$ is the adjoint. $K_\rho$ and $p$ are valued in $[0, +\infty]$, so $K_\rho(W) = +\infty$ is represented exactly. Widths are sInf/sSup of support values, equal to the book's minima and maxima on nonempty compact sets.
• A function of class $K$ is given with its representation (pieces, open domains, smooth $f_i$), and $G_f(x)$ is the set of gradients of the incident pieces, as on p. 79. Linear dependence of $G_f(x)$ is the negation of linear independence of that set.
• A run of the algorithm is a predicate on sequences $x_k, \tilde g_k, B_k, g_f(x_k), h_{k+1}$; every theorem holds for every run. Step (3) divides by $\|r_k\|$, so $r_k \ne 0$ is part of the run: a run that reaches a point where no admissible choice gives $r_k \ne 0$ has no continuation, as in the book. The stepsize $h_{k+1} = \operatorname{argmin}$ is one admissible choice, not the definition.
• The standing assumption (3.50) is a hypothesis of Theorem 3.11 as well as of 3.12 and 3.13, because the book imposes it for the rest of the section on p. 79.
• In Lemma 3.3, $\beta > 0$ is added (the bound divides by $\beta^2$). "Body" means nonempty interior. An isolated local minimum is a local minimum with a neighbourhood containing no other local minimum.
• Runs exist, so the hypotheses are not vacuous. For $f(t) = |t_1| + 2|t_2|$ (four quadrant pieces), started at the minimizer $x_0 = 0$ with null steps, one can choose $g_f(x_{k+1}) = -g_f(x_k)$ at every step, and $0 \notin G_f$ keeps $r_k \neq 0$. A formalization under which no infinite run exists would make every theorem vacuous. The run predicate also forbids the degenerate reading $h_{k+1} < 0$, under which the monotonicity condition (a) would hold vacuously on an empty segment.
• Needed infrastructure: widths of compact convex sets under linear maps (via the singular values of $B$), the effect of a rank-one dilation on widths, determinant bounds for products of dilations, and Gram-determinant continuity. The first three apply to any space-dilation method and are welcome as independent lemmas.

Selected references

• N. Z. Shor, Minimization Methods for Non-Differentiable Functions, Springer Series in Computational Mathematics 3, Springer 1985, §3.7, pp. 77–85. doi:10.1007/978-3-642-82118-9
• N. Z. Shor, N. G. Zhurbenko, A minimization method using the operation of space dilation in the direction of the difference of two successive gradients, Kibernetika (Kiev), no. 3, 51–59, 1971 (listed in the bibliography of Shor 1985; no online copy known).
• N. Z. Shor, The analysis of convergence of a gradient type method with space dilation in the direction of the difference of two successive gradients, Kibernetika (Kiev), no. 4, 48–53, 1975 (listed in the bibliography of Shor 1985; no online copy known).

Motivation

Stochastic differential equations whose coefficients grow faster than linearly (for instance stochastic volatility models such as the 3/2-model) cannot be simulated with the classical explicit Euler–Maruyama scheme: Hutzenthaler, Jentzen and Kloeden showed that its moments diverge when the drift or diffusion grows superlinearly (Proc. R. Soc. A, 2011). Implicit schemes avoid the divergence but require solving a nonlinear equation at every step. Tamed Euler schemes keep the scheme explicit and damp the coefficients by a factor depending on the step size.

• 2012: Hutzenthaler, Jentzen and Kloeden introduce a tamed Euler scheme for SDEs with superlinearly growing drift (Ann. Appl. Probab. 22, 2012).
• 2013: Sabanis extends the taming analysis for superlinearly growing drift (Electron. Commun. Probab. 18, no. 10, 2013, MR3070913).
• 2015: Hutzenthaler and Jentzen survey explicit schemes for non-globally Lipschitz coefficients (Mem. Amer. Math. Soc. 236, no. 1112, 2015); for the 3/2-model their results give $\mathcal L^p$-convergence without rate only for $p<1/2$, as Sabanis (2016, p. 2) notes.
• 2016: Sabanis (Ann. Appl. Probab. 26(4), 2016, arXiv:1308.1796v4) treats schemes with varying coefficients $b_n,\sigma_n$, covering superlinearly growing diffusion coefficients, and proves $\mathcal L^p$ convergence (Theorem 1), the $\mathcal L^p$ rate 1/2 (Theorem 2) and the uniform $\mathcal L^q$ rate 1/2 (Theorem 3). A motivating example is the $d$-dimensional analogue of the 3/2-model of stochastic volatility, $dX=\lambda X(\mu-|X|)\,dt+\xi|X|^{3/2}\,dW$, used for pricing VIX options.

Theorem 3 is the target here; Theorems 1 and 2 are the targets of companion missions.

Setting

Fix a filtered probability space $(\Omega,\{\mathcal F_t\}_{t\ge0},\mathcal F,P)$ satisfying the usual conditions, a $d_1$-dimensional Wiener martingale $W$, a horizon $T>0$, and exponents $p_0,p_1\ge2$. For $x\in\mathbb R^d$, $|x|$ is the Euclidean norm; for a $d\times d_1$ matrix $A$, $|A|$ is the Hilbert–Schmidt norm; $xy$ is the scalar product. The coefficients are Borel functions $b:[0,\infty)\times\mathbb R^d\to\mathbb R^d$ and $\sigma:[0,\infty)\times\mathbb R^d\to\mathbb R^{d\times d_1}$, and the SDE is

$$dX(t)=b(t,X(t))\,dt+\sigma(t,X(t))\,dW(t),\qquad t\in[0,T],\qquad(2.1)$$

with an $\mathcal F_0$-measurable initial value $X(0)$. With $\kappa_n(t)=\lfloor nt\rfloor/n$, the scheme (2.2) is

$$dX_n(t)=b_n(t,X_n(\kappa_n(t)))\,dt+\sigma_n(t,X_n(\kappa_n(t)))\,dW(t),\qquad X_n(0)=X(0),$$

and in this mission $b_n,\sigma_n$ are the tamed coefficients of Model 2 with $\alpha=1/2$:

$$b_n(t,x)=\frac{b(t,x)}{1+n^{-1/2}|x|^l},\qquad \sigma_n(t,x)=\frac{\sigma(t,x)}{1+n^{-1/2}|x|^l}.$$

The conditions used are: A-2, local boundedness of $b$ on balls; A-4, the coercivity bound $2xb(t,x)+(p_0-1)|\sigma(t,x)|^2\le K(1+|x|^2)$; A-5, $\mathbb E|X(0)|^{p_0}<\infty$; and A-6, the global monotonicity

$$2(x-y)(b(t,x)-b(t,y))+(p_1-1)|\sigma(t,x)-\sigma(t,y)|^2\le L|x-y|^2$$

together with the polynomial Lipschitz bound $|b(t,x)-b(t,y)|\le L(1+|x|^l+|y|^l)|x-y|$, with $l,L>0$. The $\mathfrak p$-condition asks $l\le\frac{p_0-2}{4}$ and a moment exponent $p$ with $0<p<p_1$ and $p\le\frac{p_0}{2l+1}$.

Formalization targets

Goal: Theorem 3 (p. 6)

Under A-2, A-4–A-6 and the $\mathfrak p$-condition, for every $0<q<p$ there is $C$ independent of $n$ with

$$\mathbb E\Big[\sup_{0\le t\le T}|X(t)-X_n(t)|^q\Big]\le Cn^{-q/2},\qquad n\ge1.\qquad(2.14)$$

The supremum is inside the expectation, which distinguishes it from Theorem 2.

Milestones

1. Lemma 2 (p. 8): the moments of $X$ and of $X_n$ up to order $p_0$ are bounded on $[0,T]$ uniformly in $n$ (for general coefficients $b_n,\sigma_n$ under A-1–A-5, B-2, B-3).
2. Lemma 5 (p. 19): the Gyöngy–Krylov maximal inequality. If nonnegative continuous adapted $f,g$ satisfy $\mathbb E[f_\tau\mathbb 1_{\{g_0\le c\}}]\le\mathbb E[g_\tau\mathbb 1_{\{g_0\le c\}}]$ for all $c>0$ and stopping times $\tau\le T$, then

$$\mathbb E\Big[\sup_{t\le\tau}f_t^\gamma\Big]\le\frac{2-\gamma}{1-\gamma}\,\mathbb E\Big[\sup_{t\le\tau}g_t^\gamma\Big],\qquad\gamma\in(0,1).$$

3. Lemma 4 (p. 16): $\sup_{0\le t\le T}\mathbb E|X_n(t)-X_n(\kappa_n(t))|^p\le Cn^{-p/2}$.
4. Lemma 3 (p. 15): the taming errors $\mathbb E\int_0^T|b-b_n|^p(s,X_n(\kappa_n(s)))\,ds$ and the analogue for $\sigma$ are $\le Cn^{-\alpha p}$.

Significance

Theorem 3 gives the optimal strong rate 1/2 for an explicit scheme, uniformly over the time interval, for SDEs whose diffusion coefficient may grow superlinearly. Uniform-in-time error bounds are what pathwise functionals need (running maxima, barrier options, hitting times), and via Borel–Cantelli they give the almost-sure rate of Corollary 1 in the paper. Lemma 5 is a general-purpose tool: a domination inequality of Lenglart type, used throughout stochastic analysis to pass from bounds at stopping times to bounds on running suprema.

The results are proved in the paper; to our knowledge none of them has been formalized. Formalizing them would produce the first machine-checked convergence rate of a numerical scheme for SDEs, and a Lean proof of a Lenglart-type inequality, which Mathlib does not currently contain.

Difficulty

The error $X-X_n$ is controlled by the monotonicity condition A-6, which is one-sided: it bounds $2(x-y)(b(t,x)-b(t,y))$ from above but gives no Lipschitz bound on $b$ or $\sigma$ with a constant independent of $|x|,|y|$. The first idea, estimating $\mathbb E\sup_t|X-X_n|^p$ directly by applying the Burkholder–Davis–Gundy inequality to the martingale part of the error equation, therefore does not close: it produces terms of the same order as the quantity being estimated, weighted by polynomial factors in $|X|$ and $|X_n|$, and the argument that gives Theorem 2 controls $\mathbb E|X(t)-X_n(t)|^p$ only for fixed $t$. Passing from fixed times to the supremum is what forces the loss from $p$ to $q<p$ in Theorem 3, and it is where Lemma 5 enters. On the formal side, the Itô integral against Brownian motion, Itô's formula for functions of multidimensional Itô processes and the Burkholder–Davis–Gundy inequality are not in Mathlib.

Formalization scope

Processes are indexed by $t\in[0,\infty)$ (ℝ≥0) with values in EuclideanSpace ℝ (Fin d); diffusion values live in EuclideanSpace ℝ (Fin d × Fin d₁), whose norm is the Hilbert–Schmidt norm. The stochastic integral is the Ethier–Kurtz platform relation EthierKurtz.HasBrownianItoIntegral, applied coordinatewise, with integrands cut off after $T$. Solutions are hypotheses: the theorems apply to any solution $X$ of (2.1) and any family $(X_n)_{n\ge1}$ of solutions of (2.2) on one probability space with one $W$ and one $X(0)$; existence is not claimed. Processes are required to be adapted to the filtration augmented by null sets (weaker than a complete filtration), the filtration is right-continuous, and nothing is required after $T$. Expectations and suprema are computed in $[0,\infty]$, so no statement holds through an integral defaulting to $0$. Constants $C$ depend on everything except $n$ (and, in Theorem 3, may depend on $q$); only $n\ge1$ is used, and $q>0$, $p>0$ follow the paper's $\mathcal L^p$ convention.

Two hypotheses are added to the printed statements, both because the printed proofs need them:

• $p_1>2$ in Theorem 3. The proof applies Itô's formula for $p\ge2$; an admissible $p<2$ is handled through $p'=2$, which satisfies the $\mathfrak p$-condition exactly when $p_1>2$.
• B-3 in Lemma 4. The proof uses moment bounds uniform in $n$ (Lemma 2), which need B-3. Model 2 satisfies B-3.

Lemma 2 drops the dependence clause $C=C(p,T,K,\mathbb E|X(0)|^p)$ and asserts only finiteness. Lemma 5 assumes neither that $g$ is nondecreasing nor anything beyond the printed hypotheses.

A trivializing formalization is ruled out: the supremum sits inside a Lebesgue integral in $[0,\infty]$, so the goal cannot hold because of a junk value, and the hypotheses are satisfiable (zero coefficients, constant solutions), which a sorry-free check confirms.

Needed infrastructure: Itô's formula for multidimensional Itô processes, the Burkholder–Davis–Gundy inequality (for Lemma 4), Gronwall's inequality in integral form, and stopping-time localisation for continuous adapted processes. Contributions to any of these are welcome, and all are reusable far beyond this mission. Lemma 5 is independent of the SDE part and can be attacked first.

Selected references

• S. Sabanis, Euler approximations with varying coefficients: the case of superlinearly growing diffusion coefficients, Ann. Appl. Probab. 26(4), 2083–2105, 2016. https://doi.org/10.1214/15-AAP1140 — arXiv:1308.1796v4, https://arxiv.org/abs/1308.1796v4
• M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients, Proc. R. Soc. A 467, 1563–1576, 2011. https://doi.org/10.1098/rspa.2010.0348
• M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients, Ann. Appl. Probab. 22(4), 1611–1641, 2012. https://doi.org/10.1214/11-AAP803
• M. Hutzenthaler, A. Jentzen, Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients, Mem. Amer. Math. Soc. 236, no. 1112, 2015 (reference [5] of the paper).
• S. Sabanis, A note on tamed Euler approximations, Electron. Commun. Probab. 18, no. 10, 2013. MR3070913
• I. Gyöngy, N. V. Krylov, On the rate of convergence of splitting-up approximations for SPDEs, in Stochastic Inequalities and Applications, Progr. Probab. 56, 301–321, Birkhäuser, 2003 (cited by the paper for Lemma 5). MR2073438
• N. V. Krylov, Introduction to the Theory of Diffusion Processes, Transl. Math. Monographs 142, AMS, 1995 (cited by the paper for Lemma 5). MR1311478