Numerical Analysis

Motivation

Line search methods for unconstrained minimization of a smooth function $f : \mathbb{R}^n \to \mathbb{R}$ choose a direction $d_k$ and a step $\alpha_k > 0$ and set $x_{k+1} = x_k + \alpha_k d_k$. Classical rules (Armijo, Wolfe) insist that every step decrease $f$. For quasi-Newton and conjugate gradient directions this monotonicity requirement often forces short steps, and nonmonotone line searches, which only ask for a decrease relative to some reference value built from past iterates, have been used since Grippo, Lampariello and Lucidi (1986) to let such methods take longer steps.

Zhang and Hager (SIAM J. Optim., 2004) replaced the maximum of recent function values used by Grippo et al. with a weighted average $C_k$ of all past function values. The paper proves two results: global convergence to stationary points (the companion mission) and, the subject of this mission, R-linear convergence of the function values when $f$ is strongly convex.

Timeline:

• 1986, Grippo, Lampariello, Lucidi: nonmonotone line search based on the maximum of the last $M$ function values; global convergence.
• 2002, Dai: R-linear convergence of the max-based scheme for strongly convex $f$.
• 2004, Zhang and Hager: the averaged reference value $C_k$; global convergence (Theorem 2.2) and R-linear convergence for strongly convex $f$ (Theorem 3.1).

Setting

Fix parameters $0 \le \eta_{\min} \le \eta_{\max} \le 1$, $0 < \delta < \sigma < 1 < \rho$ and $\mu > 0$. Write $g_k = \nabla f(x_k)$ and $\nabla f(x)d = \langle \nabla f(x), d\rangle$. The Nonmonotone Line Search Algorithm (NLSA) keeps weights $Q_k$ and reference values $C_k$:

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

with $\eta_k \in [\eta_{\min}, \eta_{\max}]$ chosen freely at each step. A step $\alpha_k$ is accepted either by the nonmonotone Wolfe conditions

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

or by the nonmonotone Armijo rule $\alpha_k = \bar\alpha_k \rho^{h_k}$, where $\bar\alpha_k > 0$ is a trial step and $h_k$ is the largest integer such that the first inequality holds and $\alpha_k \le \mu$. With $\eta_k = 0$ one recovers the monotone rules.

The direction assumption asks for constants $c_1, c_2 > 0$ with $g_k^{\mathsf T} d_k \le -c_1\|g_k\|^2$ and $\|d_k\| \le c_2\|g_k\|$. The function $f$ is strongly convex with constant $\gamma > 0$ if

$$f(x) \ge f(y) + \nabla f(y)(x - y) + \frac{1}{2\gamma}\|x - y\|^2\quad\text{for all } x, y.$$

Let $x^*$ be the minimizer, $\mathcal L = \{x : f(x) \le f(x_0)\}$, $d_{\max} = \sup_k\|d_k\|$, and $\bar{\mathcal L}$ the set of points within distance $\mu d_{\max}$ of $\mathcal L$.

Formalization targets

Goal: Theorem 3.1

Let $f$ be strongly convex with minimizer $x^*$, let $\nabla f$ be Lipschitz continuous on bounded sets, let $\eta_{\max} < 1$, let the directions satisfy the direction assumption at every iteration, and let $\alpha_k \le \mu$ for all $k$. Then there is $\theta \in (0,1)$ with

$$f(x_k) - f(x^*) \le \theta^k\big(f(x_0) - f(x^*)\big)\quad\text{for each } k.$$

The goal fixes no value of $\theta$: it asserts only the existence of a linear rate.

Milestones

In the paper's order of use:

1. Lemma 1.1: $f(x_k) \le C_k \le A_k$ when $g_k^{\mathsf T}d_k \le 0$ for each $k$.
2. $C_{k+1} \le C_k$, so all iterates lie in $\mathcal L$.
3. (3.4): $f(x) - f(x^*) \le \gamma\|\nabla f(x)\|^2$.
4. (2.15): $Q_{k+1} \le 1/(1 - \eta_{\max})$.
5. (3.6): $f(x_{k+1}) \le C_k - \beta\|g_k\|^2$, with

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

6. (3.7): $\|g_{k+1}\| \le b\|g_k\|$, $b = 1 + \mu c_2 L$.
7. (3.8): the explicit contraction

$$C_{k+1} - f(x^*) \le \theta\,(C_k - f(x^*)),\qquad \theta = 1 - \beta b_2(1-\eta_{\max}),\quad b_2 = \frac{1}{\beta + \gamma b^2}.$$

Here $L$ is a Lipschitz constant of $\nabla f$ on $\bar{\mathcal L}$. A further result on the same definitions is Theorem 3.2: if $f(x_k)$ converges R-linearly with ratio $\theta < \eta_{\min}$ inside a compact convex set on which $f$ is strongly convex, then the sufficient decrease condition with reference value $C_k$ holds for all large $k$.

Significance

Theorem 3.1 shows that averaging past function values costs nothing in the rate: on strongly convex functions the nonmonotone method keeps the linear rate of monotone descent, for any direction sequence satisfying the direction assumption (steepest descent, L-BFGS with bounded condition numbers, and so on). Theorem 3.2 is the converse side: for weights close enough to 1, the averaged test eventually accepts the steps of any R-linearly convergent iteration of this kind. The paper contrasts this with the max-based test of Grippo et al.

These results are proved in the paper. As far as the platform catalogue shows, no line search with the Wolfe conditions or a nonmonotone reference value has been formalized. The only related item is a monotone backtracking gradient descent rate, which is a different theorem. Machine-checked proofs would give reusable Lean statements of the nonmonotone Wolfe and Armijo rules and of the averaged reference value. They would also give the explicit constants $\beta$, $b$ and $\theta$ in checked form, and a checked record of the repair of the printed statement described below.

Difficulty

The obvious argument would show that $f(x_k) - f(x^*)$ contracts at each step. It fails: the method is nonmonotone, and $f(x_{k+1})$ may exceed $f(x_k)$. The quantity that contracts is $C_k - f(x^*)$, and only $C_k$ is controlled by the line search. The contraction must be derived by relating $\|g_k\|^2$ to $C_k - f(x^*)$ in two regimes, and the second regime needs a bound on $f(x_{k+1}) - f(x^*)$ from the gradient at the previous iterate. That bound requires the Lipschitz constant on a region containing every point the line search examines, which is why the region $\bar{\mathcal L}$ and the step bound $\mu$ enter. In Lean, the sufficient decrease (3.6) also rests on the step-length lower bounds of Lemma 2.1 for both rules, including the integer-exponent Armijo rule with its maximality condition.

Formalization scope

Conventions:

• The space is EuclideanSpace ℝ (Fin n), $f$ is ContDiff ℝ 1, and $\nabla f(x)d$ is ⟪gradient f x, d⟫_ℝ.
• A run is infinite, uses one rule throughout (Wolfe or Armijo), and has arbitrary directions subject to the stated hypotheses.
• $Q_k$ and $C_k$ are defined by recursion from the run.
• The Armijo exponent ranges over $\mathbb{Z}$ and "largest" is IsGreatest.
• $d_{\max}$ and distances are taken in $[0,\infty]$, so unbounded directions make $\bar{\mathcal L}$ the whole space.
• Strong convexity keeps the paper's constant $\gamma$, the inverse modulus.
• "Lipschitz on bounded sets" means that every bounded set admits a Lipschitz constant for $\nabla f$.

Repair of the printed statement: the paper's direction assumption holds only for all sufficiently large $k$, but Theorem 3.1 concludes (3.5) for each $k$, and its proof uses the assumption at every $k$. As printed the theorem is false: $d_0 = 0$ with the Wolfe step $\alpha_0 = 1$ gives $x_1 = x_0$, which violates (3.5) at $k = 1$ whenever $x_0 \ne x^*$. The goal therefore requires the direction assumption for every $k \ge 0$.

The step bound "there exists $\mu > 0$ with $\alpha_k \le \mu$" is stated with $\mu$ the algorithm's parameter. For the Armijo rule this holds by construction. The Wolfe rule does not use $\mu$, so nothing is lost.

Trivializing formalizations are ruled out:

• the printed hypotheses with "for each $k$" (false, as shown above);
• a $\theta$ allowed to equal 1, or an unspecified constant factor in front of $\theta^k$ (weaker than (3.5));
• a globally Lipschitz gradient, or a step bound different from the $\mu$ used to build $\bar{\mathcal L}$.

A complete development needs:

• the NLSA model and Lemma 1.1;
• the step-length lower bounds of Lemma 2.1 for both rules (via the descent lemma for Lipschitz gradients);
• the geometric bound on $Q_k$;
• boundedness of level sets of strongly convex functions.

The NLSA definitions are reusable by any later mission on nonmonotone line searches. Proofs of the milestones in any order are welcome, as are a proof of the counterexample to the printed statement and a proof that the platform theorem ConvexOptimization.strong_convexity_quadratic_lower_bound implies (3.4).

Selected references

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

Motivation

An initial value problem $y' = f(x,y)$, $y(x_0) = y_0$, is the standard way a law of motion, a growth model or a circuit equation is written down, and only a small minority of such problems can be integrated in closed form. Single-step methods advance the solution by small increments: Euler's method takes the tangent line at the current point, and the higher-order Runge-Kutta schemes refine the same idea. What justifies them is a pair of theorems: that the problem has a unique solution at all, and that the computed polygon approaches that solution as the step size decreases.

This mission is the seventh in a series formalizing the lecture notes Métodos Numéricos by Sergio Roberto de Freitas (UFMS, 2000). It covers Chapter 9, Métodos Numéricos para EDO's, up to the convergence of Euler's method.

Setting

Consider the initial value problem

$$y' = f(x,y), \\qquad y(x_0) = y_0,$$

with $f$ defined on a rectangle $D = [x_0-a, x_0+a]\\times[y_0-b, y_0+b]$.

Picard's method replaces the differential equation by the integral equation $y(x) = y_0 + \\int_{x_0}^{x} f(t, y(t))\\,dt$ and iterates it:

$$y_0(x) = y_0, \\qquad y_{k+1}(x) = y_0 + \\int_{x_0}^{x} f(t, y_k(t))\\,dt .$$

With $M = \\max_D |f|$, $N = \\max_D |\\partial f/\\partial y|$ and $h = \\min\\{a, b/M\\}$, the source states the bound $|\\varphi(x) - y_k(x)| \\le \\frac{M N^{k-1}}{k!}h^k$ on $[x_0-h, x_0+h]$, where $\\varphi$ is the exact solution.

Euler's method with step $h$ produces $y_{i+1} = y_i + h f(x_i, y_i)$, $x_i = x_0 + ih$. Its local error is $O(h^2)$ per step and its global error $O(h)$.

Target

The goal theorem is the global error statement: for an initial value problem whose exact solution $\\varphi$ is twice continuously differentiable on $[x_0, X]$ and whose right-hand side is Lipschitz in $y$ there, there is a constant $C$ such that for every step size $h \\in (0,1]$ and every grid point $x_0 + ih \\le X$,

$$|y_i - \\varphi(x_0+ih)| \\le C h .$$

This is the precise form of the source's assertion $E = O(h)$.

The milestones are the results the chapter states on the way: local existence and uniqueness for the initial value problem under a Lipschitz condition in $y$ (Teorema 9.4.1), the Picard error bound, the convergence of the Picard iterates to the solution, and the one-step Taylor identity behind Euler's method.

Significance

The global error bound is what makes Euler's method a method rather than a heuristic: it converts an arbitrary choice of step size into a guaranteed accuracy, and it exhibits the characteristic loss of one order between the local error $O(h^2)$ and the global error $O(h)$ that recurs for every single-step scheme. The existence and uniqueness theorem is the prerequisite for all of it — without uniqueness there is no well-defined object for the numerical solution to approximate — and the Picard iteration is both the constructive proof of that theorem and a method in its own right.

Difficulty

The global error proof has to handle the interaction of two error sources: the local truncation error at each step and the amplification of previous errors by the Lipschitz constant, which gives the discrete Gronwall recursion $e_{i+1} \\le (1+hL)e_i + Ch^2$. Carrying that recursion to the closed-form bound, and doing it uniformly in $h$, is the core of the mission. The Picard bound requires the iterates to remain inside the rectangle so that the hypotheses on $f$ continue to apply — that is what the choice $h = \\min\\{a, b/M\\}$ is for — and an induction on $k$ with the factorial denominator.

Formalization scope

The unknown is a real function of a real variable and $f$ is a function of two real variables. A solution is a function $\\varphi$ satisfying $\\varphi(x_0)=y_0$ and having derivative $f(x,\\varphi(x))$ at every point of the relevant interval. The Lipschitz condition in the second variable is stated explicitly, $|f(x,u)-f(x,v)| \\le L|u-v|$, in place of the source's hypothesis that $\\partial f/\\partial y$ is bounded; for the global-error goal it is assumed for all real $u,v$ and all $x$ in the interval, and in the Picard milestones only inside the rectangle. Uniqueness is stated as agreement of any two solutions on the interval, not as uniqueness of a function on all of $\\mathbb{R}$, since a solution is unconstrained outside the interval. Euler's iterates and the Picard iterates are explicit recursive definitions; the Picard integral is the interval integral, which is total, so no integrability side condition appears in the definition. The constant $C$ in the goal is asserted to exist and to be positive, uniformly over step sizes in $(0,1]$; the restriction $h \\le 1$ merely normalizes the range of steps considered.

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 9, Métodos Numéricos para EDO's, pp. 183–213. (Course notes supplied with this mission.)

Motivation

Link-based ranking of web pages rests on a small number of linear-algebraic facts. Kleinberg's HITS algorithm (Hyperlink-Induced Topic Search), introduced in Kleinberg 1999, assigns each page of a hyperlinked collection two scores: an authority score, high for pages that many good hubs point to, and a hub score, high for pages that point to many good authorities. The scores are computed by a simple alternating iteration, and the paper's justification for using them is a convergence theorem: the iteration settles, and its limits are the principal eigenvectors of two symmetric matrices built from the link structure. HITS, together with PageRank (Brin–Page 1998), founded the study of link analysis, and the eigenvector characterization is what later work on the stability and uniqueness of such rankings starts from (for example Ng–Zheng–Jordan 2001 and Farahat et al. 2006).

This mission formalizes that convergence theorem, §3 of the paper, together with the steps of its proof.

Setting

A collection of $n$ hyperlinked pages is a directed graph $G=(V,E)$ with $V=\{p_1,\dots,p_n\}$; an edge $(p,q)\in E$ means that page $p$ links to page $q$. Its adjacency matrix $A$ is the $n\times n$ matrix with $A_{ij}=1$ if $(p_i,p_j)\in E$ and $A_{ij}=0$ otherwise. Self-links are allowed.

Each page $p$ carries an authority weight $x^{\langle p\rangle}$ and a hub weight $y^{\langle p\rangle}$, collected into vectors $x,y\in\mathbb R^n$. Two operations update them:

$$\mathcal I:\ x^{\langle p\rangle}\leftarrow\sum_{q:(q,p)\in E}y^{\langle q\rangle},\qquad \mathcal O:\ y^{\langle p\rangle}\leftarrow\sum_{q:(p,q)\in E}x^{\langle q\rangle}.$$

Normalizing a vector rescales it so that its squares sum to $1$.

The procedure Iterate$(G,k)$ starts from $x_0=y_0=z=(1,\dots,1)$. In round $i=1,\dots,k$ it applies $\mathcal I$ to $y_{i-1}$, obtaining $x_i'$; applies $\mathcal O$ to $x_i'$, obtaining $y_i'$; and normalizes $x_i'$ and $y_i'$ to get $x_i$ and $y_i$. It returns $(x_k,y_k)$.

For a real symmetric matrix $M$, list its eigenvalues with multiplicity in order of decreasing absolute value, $\lambda_1(M),\lambda_2(M),\dots,\lambda_n(M)$. The paper's standing Assumption (†) is

$$|\lambda_1(M)|>|\lambda_2(M)|,$$

so that the eigenspace of $\lambda_1(M)$ is a line. A principal eigenvector $\omega_1(M)$ is a unit vector spanning that line; it is determined up to sign.

Formalization targets

Goal: Theorem 3.2 with the limits of Theorem 3.1

Assume (†) for $A^{\top}A$ and for $AA^{\top}$. Then there are $x^*,y^*\in\mathbb R^n$ with

$$x_k\to x^*,\qquad y_k\to y^*,$$ $$A^{\top}A\,x^*=\lambda_1(A^{\top}A)\,x^*,\quad AA^{\top}y^*=\lambda_1(AA^{\top})\,y^*,\quad \|x^*\|_2=\|y^*\|_2=1.$$

The authority scores converge to the principal eigenvector of the co-citation matrix $A^{\top}A$, and the hub scores to that of the bibliographic coupling matrix $AA^{\top}$.

Milestones: the proof of Theorem 3.1

1. The operations are linear: $\mathcal I(y)=A^{\top}y$ and $\mathcal O(x)=Ax$.
2. Closed forms: for $k\ge1$, $x_k$ is the unit vector along $(A^{\top}A)^{k-1}A^{\top}z$ and $y_k$ the unit vector along $(AA^{\top})^kz$.
3. Normalized power iteration: for symmetric $M$ with positive dominant simple eigenvalue and $v$ not orthogonal to $\omega_1(M)$, $M^kv/\|M^kv\|$ converges to a principal eigenvector.
4. A symmetric matrix with non-negative entries satisfying (†) has a principal eigenvector with non-negative entries.
5. $z$ is not orthogonal to $\omega_1(AA^{\top})$.
6. If $\lambda_1(A^{\top}A)\ne0$, then $A^{\top}z$ is not orthogonal to $\omega_1(A^{\top}A)$.
7. Theorem 3.1: $(x_k)$ and $(y_k)$ converge.

A companion item states Proposition 6.1 (§6): $AA^{\top}$ and $A^{\top}A$ have the same characteristic polynomial, and $A$ maps each eigenvector of $A^{\top}A$ with nonzero eigenvalue to an eigenvector of $AA^{\top}$ with the same eigenvalue.

Significance

The theorem turns a heuristic, mutually reinforcing update into a well-defined object: the HITS scores of a graph are the principal eigenvectors of $A^{\top}A$ and $AA^{\top}$, independent of the number of rounds run, and computable by any eigenvector algorithm. It also locates where the definition is fragile: when (†) fails, the limit depends on the starting vector, which is the starting point of the later stability analyses cited above.

The result has been proved since 1999; its proof in the paper is a short appeal to the power method and to a Perron-type non-negativity fact. What this mission adds is a machine-checked version of the whole chain, from the combinatorial update rules on a directed graph to the spectral characterization, with the hypotheses made explicit. Two generic power-method theorems are already proved on the platform (Trefethen–Bau, Lecture 27) and are included as references; to our knowledge no formalization of HITS, of Assumption (†), or of the non-negativity of principal eigenvectors of non-negative symmetric matrices exists in Mathlib or on the platform.

Difficulty

The obvious argument — "the power method converges" — is only half of the proof. The power method converges from a start vector that has a nonzero component along the dominant eigenvector, and the start vectors here are fixed: $z$ for the hubs and $A^{\top}z$ for the authorities. Showing that they are not orthogonal to the principal eigenvectors needs a Perron-type fact (the principal eigenvector of a non-negative symmetric matrix can be chosen non-negative), which in turn needs that the dominant eigenvalue is the spectral radius itself and not its negative. For the authorities the start vector is $A^{\top}z$, not $z$, and the argument must also use $\lambda_1(A^{\top}A)\ne0$: when $A^{\top}A$ has a zero dominant eigenvalue the iterates are $0$. A second difficulty is the passage from the rescaled power iterates $\lambda_1^{-k}M^kv$ of the referenced platform theorem to the normalized iterates of the algorithm, and from a statement on an abstract inner product space to one about matrices on $\mathbb R^n$.

Formalization scope

Pages are Fin n; the link relation is a decidable relation E : Fin n → Fin n → Prop; adjMatrix E is the real 0/1 matrix (Mathlib's SimpleGraph.adjMatrix is undirected and loopless and is not used). The operations $\mathcal I$, $\mathcal O$ are defined as neighbour sums, not as matrix products, so that milestone 1 has content. normalize v divides by the Euclidean norm and sends $0$ to $0$; no statement is true only because of that convention. hitsIter E k is $(x_k,y_k)$, with the $\mathcal O$ step applied to the unnormalized $x_k'$ as on the page. Convergence is in the product topology of $\mathbb R^n$.

Assumption (†) is encoded as Dagger M: an eigenvalue $\lambda_1$ with one-dimensional eigenspace, $|\mu|<|\lambda_1|$ for every other eigenvalue $\mu$, and $\lambda_1\ne0$. The last clause is the paper's own reading ("as dictated by Assumption (†)", p. 11); without it, the one-page graph with no link satisfies (†) vacuously and Theorem 3.2 fails. Theorems 3.1 and 3.2 assume (†) for both $A^{\top}A$ and $AA^{\top}$, as the paper does. "The principal eigenvector" is a predicate, IsPrincipalEigenvector, true of both unit vectors in the principal eigenspace. Added hypotheses, each disclosed in its item: $\lambda_1>0$ in the power-iteration milestone (the claim fails for $\lambda_1<0$; it costs nothing for the positive semidefinite $A^{\top}A$, $AA^{\top}$); the closed form is stated for $k\ge1$. The non-negativity of $x^*,y^*$ is true but not part of Theorem 3.2 and is not asserted.

A formalization in which $\mathcal I$ and $\mathcal O$ are matrix products by definition, or in which Theorem 3.1 alone is the goal (a sequence converging to $0$ would satisfy it), is ruled out: the goal pairs convergence with membership of the limit in the principal eigenspace as a unit vector.

Useful infrastructure, reusable beyond this mission: normalized power iteration for real symmetric matrices, and the non-negativity of the principal eigenvector of a non-negative symmetric matrix under a simple dominant eigenvalue. Contributions of either, and of the bridge between Matrix.toLin'-eigenspaces and Matrix.IsHermitian.eigenvalues, are welcome.

Selected references

• J. M. Kleinberg, Authoritative sources in a hyperlinked environment, Journal of the ACM 46(5):604–632, 1999. https://doi.org/10.1145/324133.324140
• S. Brin and L. Page, The anatomy of a large-scale hypertextual Web search engine, Computer Networks and ISDN Systems 30:107–117, 1998. https://doi.org/10.1016/S0169-7552(98)00110-X
• A. Y. Ng, A. X. Zheng and M. I. Jordan, Stable algorithms for link analysis, Proc. SIGIR 2001, 258–266. https://doi.org/10.1145/383952.384003
• A. Farahat, T. LoFaro, J. C. Miller, G. Rae and L. A. Ward, Authority rankings from HITS, PageRank, and SALSA: existence, uniqueness, and effect of initialization, SIAM J. Sci. Comput. 27(4):1181–1201, 2006. https://doi.org/10.1137/S1064827502412875
• L. N. Trefethen and D. Bau III, Numerical Linear Algebra, SIAM, 1997, Lecture 27. https://doi.org/10.1137/1.9780898719574

Motivation

The subgradient method for a nonsmooth convex function converges, but slowly: when the level sets of the objective are elongated, the subgradient is nearly orthogonal to the direction towards the minimum, and the method zigzags. For smooth functions the remedy is a change of metric (Newton and quasi-Newton methods); for nonsmooth functions no Hessian exists to supply one. N. Z. Shor's answer was to learn a metric from the subgradients themselves: after each step, stretch the space along the latest (transformed) subgradient, so that components of future subgradients parallel to it are damped. These subgradient methods with space dilation along the gradient (SDG methods) are the ancestors of Shor's r-algorithm and of the ellipsoid method, which the book (p. 49) describes as a special case of the same family and which Khachiyan later used to show that linear programming is solvable in polynomial time.

Section 3.4 of Shor's monograph (Springer 1985, translated by K. C. Kiwiel and A. Ruszczyński, doi:10.1007/978-3-642-82118-9) proves that, under a two-sided condition on the objective, a suitable SDG method decreases function values at the speed of a geometric progression whose ratio is invariant under nonsingular linear changes of variables. This mission formalizes that chain of results.

Setting

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

$$R_\alpha(\xi)\,x = x + (\alpha - 1)(x,\xi)\,\xi ,$$

which multiplies the component of $x$ along $\xi$ by $\alpha$ and leaves the orthogonal complement fixed.

Let $f : E_n \to \mathbb{R}$ and let $g : E_n \to E_n$ be a generalized gradient: a subgradient of $f$ when $f$ is convex, an almost-gradient when $f$ is almost differentiable. The SDG method starts from $x_0$ and a nonsingular operator $B_0 = A_0^{-1}$. At step $k = 0, 1, \dots$: if $g(x_k) = 0$ it stops; otherwise it forms the transformed gradient $\tilde g_k = B_k^* g(x_k)$, the direction $\xi_{k+1} = \tilde g_k/\|\tilde g_k\|$, and

$$x_{k+1} = x_k - h_{k+1} B_k \xi_{k+1}, \qquad B_{k+1} = B_k R_{1/\alpha_{k+1}}(\xi_{k+1}), \qquad A_{k+1} = R_{\alpha_{k+1}}(\xi_{k+1}) A_k ,$$

with a stepsize $h_{k+1}$ and a dilation coefficient $\alpha_{k+1}$. So $A_k$ is the accumulated space transformation, $B_k = A_k^{-1}$, and each step is a subgradient step for $\varphi_k(y) = f(B_k y)$ in the variables $y = A_k x$.

The quantitative results assume, for a point $x^*$ and the ball $S_d = \{x : \|x - x^*\| \le d\}$, the two-sided condition

$$N\,[f(x) - f(x^*)] \le (g(x),\, x - x^*) \le M\,[f(x) - f(x^*)], \qquad x \in S_d, \quad M > N > 0. \qquad (3.18)$$

For a convex function the lower inequality holds with $N = 1$; the upper one bounds how far $f$ is from a positively homogeneous function around $x^*$.

Formalization targets

Goal: Theorem 3.4

Under (3.18), with $B_0 = I$, $x_0 \in S_d$, stepsizes $h_{k+1} = \frac{2MN}{M+N}\frac{f(x_k)-f(x^*)}{\|\tilde g_k\|}$, a constant coefficient $1 < \alpha \le \frac{M+N}{M-N}$, and $G$ a bound for $\|g\|$ on $S_d$: there are $c > 0$ and indices $k_1 < k_2 < \cdots$ with

$$f(x_{k_p}) - f(x^*) \le c\,\alpha^{-k_p/n},$$

and for every $k \ge 1$

$$\min_{0 \le i \le k-1}\,[f(x_i) - f(x^*)] \le \frac{G\sqrt{k(\alpha^2-1)}\,d}{N\sqrt{\alpha^{2k/n}-1}} .$$

Milestones

1. Eq. (3.4): $\|R_\alpha(\xi)x\| = \sqrt{\|x\|^2 + (\alpha^2-1)(x,\xi)^2}$.
2. Theorem 3.1: if $\|g(x_k)\| \le d$ and $1+\delta \le \alpha_k \le \alpha^*$, then $\|\tilde g_{k_p}\| < c\,(\prod_{j\le k_p}\alpha_j)^{-1/n}$ along a subsequence.
3. Theorem 3.2: for constant $\alpha > 1$ and $B_0 = I$, $\min_{0\le r\le k-1}\|\tilde g_r\| \le d\sqrt{k(\alpha^2-1)}/\sqrt{\alpha^{2k/n}-1}$.
4. Theorem 3.3: under (3.18) and the rules above, $\|A_k(x_k - x^*)\| \le d$ for all $k$.

Significance

The result shows that one fixed rule, depending only on $M$, $N$ and $n$, yields linear convergence of function values for every objective satisfying (3.18), at a ratio $\alpha^{-1/n}$ that does not deteriorate when the problem is badly scaled: the method, and hence its rate, is invariant under nonsingular linear changes of variables. This is the property the ellipsoid method inherits (Section 3.8 of the book), and the same space-dilation machinery drives the r-algorithm (Section 3.7), still used for large nonsmooth problems such as Lagrangian duals of integer programs.

The theorems are proved in the book; none of them, and no space-dilation method, has a machine-checked proof to our knowledge, and the platform has no statement about variable-metric subgradient methods. The mission provides a reusable formal model of the SDG iteration, the eigenvalue-growth arguments behind Theorems 3.1–3.2, and the one-step invariant of Theorem 3.3. The formalization also corrects two points of the printed text (see Formalization scope).

Difficulty

Theorem 3.3 is a one-step computation, but Theorems 3.1 and 3.2 are not: they relate the size of the transformed gradients to the growth of the singular values of $A_k$, whose determinant is $\prod_j \alpha_j$. A bound on $\|\tilde g_k\|$ at a single step says nothing, since the dilations can concentrate in few directions; the argument has to control the largest singular value of $A_k$ over many steps against the geometric-mean lower bound $(\det A_k)^{1/n}$. The dimension $n$ enters the rate exactly through this comparison. Theorem 3.4 then needs the invariant of Theorem 3.3 to keep every iterate inside $S_d$, where (3.18) and the bound $G$ are available.

Formalization scope

• $E_n$ is EuclideanSpace ℝ (Fin n) with $n \ge 1$ in the rate statements. Operators are continuous linear maps; $B_0$ is a continuous linear equivalence and $A_0$ its inverse. The state $(x_k, B_k, A_k)$ is produced by a defined recursion sdg, not assumed; $\tilde g_k$ is gTilde.
• The stepsize rule receives the index, the current point and $\tilde g_k$; the dilation coefficient at step $k$ is α (k+1). When $g(x_k) = 0$ the state is repeated, which encodes the book's stop; no division by zero is used.
• The book states Theorems 3.1–3.4 for almost differentiable $f$ with $g$ an almost-gradient; the proofs use only the bounds on $\|g\|$ and (3.18), so the Lean statements quantify over every map $g$ with those properties. $G$ is any bound for $\|g\|$ on $S_d$ in place of the maximum.
• Theorems 3.2–3.4 take $B_0 = I$, as their proofs do; Theorem 3.1 allows any nonsingular $B_0$.
• Two corrections to the printed statements, both following the proofs: Theorem 3.4's record bound carries the factor $1/N$ that the proof derives, and the record minima in Theorems 3.2 and 3.4 range over the indices $0, \dots, k-1$ that the proof controls, rather than $1, \dots, k$.
• Constants $c$ and subsequences are existential and chosen after the data of the run, before the index $p$. A statement placing $c$ after $p$, or dropping (3.18) on $S_d$, would be trivially true or false and is ruled out.

Useful infrastructure, reusable for the r-algorithm and ellipsoid chapters: identities for $R_\alpha(\xi)$, determinants and singular values of products of rank-one dilations, and the invariance of the SDG iteration under a linear change of variables. Proofs of any milestone, and alternative arguments for Theorem 3.1, are welcome.

Selected references

• N. Z. Shor, Minimization Methods for Non-Differentiable Functions, Springer Series in Computational Mathematics 3, Springer, 1985, §§3.2–3.4, pp. 49–62 (translated by K. C. Kiwiel and A. Ruszczyński). https://doi.org/10.1007/978-3-642-82118-9

Motivation

Many problems in optimization and equilibrium modelling reduce to a system of equations $F(x) = 0$ with $F : \mathbb{R}^n \to \mathbb{R}^n$ that is continuous and locally Lipschitz but not differentiable. Complementarity problems rewritten through the min or Fischer–Burmeister functions, Karush–Kuhn–Tucker systems of nonlinear programs, and the gradients of augmented Lagrangians all have this form. Newton's method $x^{k+1} = x^k - F'(x^k)^{-1}F(x^k)$ cannot be applied verbatim, because $F'(x^k)$ need not exist.

Qi and Sun (Math. Programming 58 (1993) 353–367) replaced the Jacobian by an arbitrary element of Clarke's generalized Jacobian and proved that the resulting method converges under a condition they called semismoothness, extending Mifflin's notion (SIAM J. Control Optim. 15 (1977)) from functionals to maps. The paper has two convergence results. Theorem 3.2 is local: near a semismooth root with nonsingular generalized Jacobian the method converges superlinearly. Theorem 3.3, the target of this mission, is global in the Newton–Kantorovich sense: explicit constants on a ball $S$ around the starting point guarantee that the iteration never leaves $S$, that $F$ has exactly one zero in $S$, and that the iterates converge to it with a computable error bound. The authors describe it as "an extension of the classical Newton-Kantorovich theorem" (p. 361), which is stated for smooth maps in Ortega and Rheinboldt's monograph.

Timeline.

• 1948: Kantorovich proves semilocal convergence of Newton's method for smooth operators in Banach spaces.
• 1975–1983: Clarke introduces the generalized gradient and generalized Jacobian of a locally Lipschitz map (Optimization and Nonsmooth Analysis, Wiley 1983).
• 1977: Mifflin defines semismooth functionals.
• 1990: Pang proves convergence of a B-derivative Newton method under a strong Fréchet derivative at the solution.
• 1993: Qi and Sun (this paper) prove local and global convergence of the generalized-Jacobian Newton method for semismooth maps.

Setting

Work in $\mathbb{R}^n$ with the Euclidean norm $\|\cdot\|$; for a linear map $V$, $\|V\|$ is the induced operator norm. Let $F : \mathbb{R}^n \to \mathbb{R}^n$ be locally Lipschitz. Write $D_F$ for the set of points where $F$ is differentiable and $JF(y)$ for the derivative at $y \in D_F$.

• The generalized Jacobian of $F$ at $x$ is

$$\partial F(x) = \mathrm{co}\Big\{\lim_{i\to\infty} JF(x_i) : x_i \to x,\ x_i \in D_F\Big\}.$$

• The directional derivative is $F'(x;h) = \lim_{t\downarrow 0}\,(F(x+th)-F(x))/t$.
• $F$ is semismooth at $x$ if it is Lipschitz near $x$ and, for every $h$, $V h'$ has a limit as $V \in \partial F(x+th')$, $h' \to h$, $t \downarrow 0$. Semismoothness implies that $F'(x;h)$ exists.
• A run of the nonsmooth Newton method (3.2) from $x^0$ is a pair of sequences $(x^k)$, $(V_k)$ with $V_k \in \partial F(x^k)$ and $V_k(x^{k+1} - x^k) = -F(x^k)$ for all $k$; any element of $\partial F(x^k)$ may be chosen.

In Lean these are NonsmoothNewton.Global.clarkeJac, dirDeriv, SemismoothAt and IsNewtonRun, over EuclideanSpace ℝ (Fin n).

Fix $x^0$, $r \ge 0$, the closed ball $S = \{x : \|x - x^0\| \le r\}$, and constants $\beta, \gamma, \delta$ with $\alpha = \beta(\gamma+\delta)$.

Formalization targets

Goal: Theorem 3.3 (global convergence)

Assume $F$ is semismooth at every point of $S$ and, for all $x, y \in S$ and $V \in \partial F(x)$: $V$ is nonsingular,

$$\|V^{-1}\| \le \beta,\qquad \|V(y-x) - F'(x;y-x)\| \le \gamma\|y-x\|,\qquad \|F(y)-F(x)-F'(x;y-x)\| \le \delta\|y-x\|,$$

with $\alpha < 1$ and $\beta\|F(x^0)\| \le r(1-\alpha)$. Then every run of (3.2) from $x^0$ stays in $S$, every $V_k$ is nonsingular, $F$ has a unique zero $x^*$ in $S$, $x^k \to x^*$, and

$$\|x^k - x^*\| \le \frac{\alpha}{1-\alpha}\,\|x^k - x^{k-1}\|,\qquad k = 1, 2, \dots \tag{3.4}$$

Milestones (steps of the proof, p. 360)

1. The first step: $\|x^1 - x^0\| \le \beta\|F(x^0)\| \le r(1-\alpha)$, so $x^1 \in S$.
2. One-step contraction: for consecutive Newton points with $x^{k-1}, x^k \in S$, $\|x^{k+1} - x^k\| \le \alpha\|x^k - x^{k-1}\|$.
3. All iterates remain in $S$, with $\|x^{k+1} - x^k\| \le r\alpha^k(1-\alpha)$.
4. A run that stays in $S$ and converges has uniformly bounded $\|V_k\|$, and its limit is a zero of $F$.
5. $F$ has at most one zero in $S$.

Significance

The theorem certifies, from data checkable on a single ball, that a solution exists, where it is, that it is unique there, and how far the current iterate is from it; (3.4) is an a-posteriori stopping criterion. It holds for every choice of $V_k \in \partial F(x^k)$, which is what implementations need, since they compute one element of $\partial F$ and not the whole set. The result underlies the global and semilocal analysis of semismooth Newton methods for complementarity problems, variational inequalities and nonsmooth KKT systems, a line that continued through the 1990s and 2000s.

The theorem is proved on paper; no machine-checked version is known. The platform has no Newton–Kantorovich theorem, smooth or nonsmooth, and no Clarke generalized Jacobian. A complete formalization would supply both, and the definitions layer here (generalized Jacobian, one-sided directional derivative, semismoothness, Newton runs) is shared with the two companion missions on Qi and Sun's local convergence theorem and on semismoothness of augmented Lagrangian gradients.

Difficulty

The Newton map is set-valued: $x^{k+1}$ depends on the choice of $V_k$, so the Banach fixed-point theorem for a single contraction does not apply directly, and the argument must hold for every sequence of choices. The contraction estimate needs both consecutive steps to start inside $S$, so containment in $S$ and the geometric decay of the steps must be established together. Identifying the limit as a zero requires a uniform bound on $\|V_k\|$, which is not a hypothesis: it has to come from local Lipschitz continuity through the structure of the generalized Jacobian. Uniqueness uses an element $V^* \in \partial F(x^*)$, whose existence rests on Rademacher's theorem. Finally, the directional derivative in the hypotheses is only meaningful because semismoothness makes it exist.

Formalization scope

The space is EuclideanSpace ℝ (Fin n), so the norm is Euclidean, as in the paper (p. 356), and $\|V^{-1}\|$ is the operator norm. Local Lipschitzness is global (LocallyLipschitz F), the standing assumption of Section 3. Conventions:

• $\partial F(x)$ is the convex hull of limits of fderiv along sequences in $D_F$; no closure is taken (the limit set is compact for locally Lipschitz $F$).
• "$V$ nonsingular, $\|V^{-1}\| \le \beta$" is the existence of a two-sided inverse $W$ with $\|W\| \le \beta$; Ring.inverse is not used.
• $F'(x;h)$ is dirDeriv, a limUnder along $t \to 0^+$, used only at points of $S$, where semismoothness makes the limit exist.
• The run is a relation, not a function; every choice of $V_k$ is covered, and nonsingularity of each $V_k$ is a conclusion.
• The radius condition $r \ge 0$ is an explicit hypothesis. Without it, $r < 0$ and $\beta < 0$ would satisfy all other hypotheses vacuously while $x^0 \notin S$.
• The paper's third inequality sits under "for any $V \in \partial F(x)$"; it is stated without $V$, which is equivalent because $\partial F(x) \ne \emptyset$.
• (3.4) is stated at index $k+1$ for $k \ge 0$, avoiding natural-number subtraction.
• The paper's statements contain no $o(\cdot)$ or $O(\cdot)$; all constants are explicit and fixed before the quantifiers over points of $S$.

The goal is the full four-part conclusion: containment, existence, uniqueness, and convergence with (3.4). A formalization that proves only that the iterates converge to some zero, or that treats $\partial F(x)$ as possibly empty so that the hypotheses become vacuous, is not the theorem. The hypotheses are satisfiable by genuinely nonsmooth maps, e.g. $F(x) = x + \tfrac{1}{10}|x| - c$ on $\mathbb{R}$ with a ball containing the kink.

A complete development needs: nonemptiness and local boundedness of the generalized Jacobian (Rademacher's theorem, available in Mathlib as LipschitzWith.ae_differentiableAt); geometric-series and Cauchy-sequence arguments in a complete space. The generalized-Jacobian facts are reusable well beyond this mission. Contributions of any of the milestones, or of these general facts as separate lemmas, are welcome.

Selected references

• L. Qi, J. Sun, A nonsmooth version of Newton's method, Mathematical Programming 58 (1993) 353–367. https://doi.org/10.1007/BF01581275
• F. H. Clarke, Optimization and Nonsmooth Analysis, Wiley, New York, 1983. https://doi.org/10.1137/1.9781611971309
• R. Mifflin, Semismooth and semiconvex functions in constrained optimization, SIAM J. Control Optim. 15 (1977) 959–972. https://doi.org/10.1137/0315061
• J. M. Ortega, W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, 1970. https://doi.org/10.1137/1.9780898719468
• J.-S. Pang, Newton's method for B-differentiable equations, Mathematics of Operations Research 15 (1990) 311–341. https://doi.org/10.1287/moor.15.2.311