Optimization

Motivation

Decisions under uncertainty are often modelled by a finite set of scenarios: possible futures $s\in S$ with weights $p_s$, each with its own optimization problem. Solving every scenario separately is easy but useless in practice, because the decision taken today cannot depend on which scenario will be revealed tomorrow. Progressive hedging, introduced by Rockafellar and Wets in Scenarios and policy aggregation in optimization under uncertainty (IIASA Working Paper WP-87-119, 1987; Mathematics of Operations Research 16(1), 1991), restores this nonanticipativity by iterating between scenario-wise solves and an averaging step, guided by price systems that act as multipliers. It is the standard decomposition method for multistage stochastic programs and is implemented in widely used software such as PySP/mpi-sppy. This mission formalizes the paper's convergence theorem in the convex case, Theorem 5.1, together with the duality theory of §4 and the proof steps of §5 it rests on.

Setting

Let $S$ be a finite set of scenarios with weights $p_s>0$, $\sum_s p_s=1$. For each $s$ the scenario subproblem $(P_s)$ is to minimize $f_s(x)$ over $x\in C_s\subset\mathbb R^n$, where $C_s$ is nonempty and closed and $f_s$ is locally Lipschitz with bounded level sets $\{x\in C_s\mid f_s(x)\le\alpha\}$. The convex case is when every $f_s$ and $C_s$ is convex.

A decision vector is split by time, $x=(x_1,\dots,x_T)$. A policy is a map $X:S\to\mathbb R^n$; the policy space $\mathcal E$ carries the inner product $\langle X,Y\rangle=\sum_s p_s X(s)\cdot Y(s)$ and norm $\|X\|=\langle X,X\rangle^{1/2}$. At each time $t$ the scenarios are partitioned into bundles $A\in\mathcal A_t$ of scenarios that cannot yet be told apart. The aggregation $\hat X=JX$ replaces $X_t(s)$ by its $p$-weighted average over the bundle of $s$ at time $t$, and $K=I-J$. The implementable policies are $\mathcal N=\{X\mid X_t$ constant on each $A\in\mathcal A_t\}$, the price systems are $\mathcal M=\{W\mid JW=0\}=\mathcal N^\perp$, and the admissible policies are $\mathcal C=\{X\mid X(s)\in C_s\ \forall s\}$. With $F(X)=\sum_s p_s f_s(X(s))$ the problem is

$$(P)\qquad \text{minimize } F(X) \text{ subject to } X\in\mathcal C\cap\mathcal N .$$

Its Lagrangian is $L(X,W)=F(X)+\langle X,W\rangle$ on $\mathcal C\times\mathcal M$, the dual is $(D)$: maximize $G(W)=\inf_{X\in\mathcal C}L(X,W)$ over $W\in\mathcal M$ with $G(W)>-\infty$.

The progressive hedging algorithm fixes $r>0$, starts from any $X^0$ and any $W^0\in\mathcal M$, and in iteration $\nu$ sets $\hat X^\nu=JX^\nu$, lets $X^{\nu+1}(s)$ minimize

$$f_s(x)+x\cdot W^\nu(s)+\tfrac12 r\,|x-\hat X^\nu(s)|^2\quad\text{over } x\in C_s$$

for every scenario, and updates $W^{\nu+1}=W^\nu+rKX^{\nu+1}$.

Formalization targets

Goal: Theorem 5.1

In the convex case with exact minimization, the sequences $\{\hat X^\nu\}$ and $\{W^\nu\}$ are bounded if and only if (P) and (D) have optimal solutions. In that case there are optimal $X^*$ for (P) and $W^*$ for (D) with

$$\hat X^\nu\to X^*,\qquad W^\nu\to W^*,$$

and, with $\|(X,W)\|_r=(\|X\|^2+r^{-2}\|W\|^2)^{1/2}$,

$$\|(\hat X^{\nu+1},W^{\nu+1})-(X^*,W^*)\|_r\le\|(\hat X^{\nu},W^{\nu})-(X^*,W^*)\|_r$$

for all $\nu\ge0$, strictly unless $(\hat X^\nu,W^\nu)=(X^*,W^*)$, and the step lengths $\|(\hat X^{\nu+1},W^{\nu+1})-(\hat X^\nu,W^\nu)\|_r$ are nonincreasing from $\nu=1$ on.

Milestones

1. Propositions 3.1–3.3: existence for the scenario subproblems, the lower bound $\hat\alpha=\min_{\mathcal C}F=E\{\min(P_s)\}\le\min(P)$, well-posedness of every modified subproblem (unique solution in the convex case), and compact level sets of $F$ on $\mathcal C$.
2. Theorem 4.2 and Proposition 4.4 (in the convex case): the scenario-wise subgradient conditions $-W^*(s)\in\partial f_s(X^*(s))+N_{C_s}(X^*(s))$ are equivalent to a saddle point of $L$, and the analogous conditions are necessary and sufficient for the iterates.
3. Proposition 4.5 and Theorem 4.6: the perturbation function $\Phi(U)=\min\{F(X)\mid X\in\mathcal C,\ KX=U\}$ and the duality $-\infty<\min(P)=\sup(D)$, $\operatorname{argmax}(D)=-\partial\Phi(0)$.
4. Proposition 5.3: each iteration is the unique saddle point of a proximal saddle function.
5. From the proof of Theorem 5.1: firm nonexpansiveness (5.26) of the iteration map in $\|\cdot\|_r$, and the identification of its fixed points with primal–dual optimal pairs.

Significance

Theorem 5.1 is the guarantee behind progressive hedging in the convex case: solving only small scenario problems and averaging, the method finds a nonanticipative optimal policy and an optimal price system whenever these exist, with a distance to the solution that never increases. The dual sequence $W^\nu$ converges to a solution of (D), which gives the prices of information that the paper interprets economically. Boundedness of the iterates is an exact test for solvability.

The result is classical and fully proved on paper; its proof reduces the algorithm to Rockafellar's proximal point algorithm for a maximal monotone operator. No machine-checked version is known to exist. The mission produces a formal statement and, eventually, a proof of the full chain from scenario data to convergence, including the duality theorem of §4, which is of independent use for scenario-based stochastic programming. A related platform mission in this series treats the nonconvex case (limits of the algorithm satisfy first-order conditions).

Difficulty

The algorithm is not a proximal point method in the variable $X$: it acts on the pair $(\hat X,W)$, and the averaged policy $\hat X^\nu$ is in general not admissible while $X^\nu$ is in general not implementable. The step that fails in a direct attempt is showing that the map $(\hat X^\nu,W^\nu)\mapsto(\hat X^{\nu+1},W^{\nu+1})$ is firmly nonexpansive, which requires recognizing it as the resolvent of the monotone operator of a saddle function restricted to $\mathcal N\times\mathcal M$, in the rescaled norm $\|\cdot\|_r$ rather than in the original product norm. A second difficulty is duality: identifying the fixed points of the iteration with optimal pairs of (P) and (D) needs the convex duality theory of §4 for an extended-real-valued perturbation function without any constraint qualification.

Formalization scope

$\mathbb R^n$ is EuclideanSpace ℝ (Fin n); policies are functions S → EuclideanSpace ℝ (Fin n) with S a Fintype. Time periods are a monotone map assigning each coordinate to its period, and each $\mathcal A_t$ is a Setoid S; no refinement between periods is assumed. The standing assumptions of §3 are fields of the model structure, and the convex case is a separate hypothesis.

The inner product (2.2), norm (2.3) and $\|\cdot\|_r$ (5.2) are explicit weighted sums. Mathlib's norm on a function type is the sup norm; it is used only for topology (boundedness, convergence, local Lipschitz continuity, compactness), which coincides with that of (2.3) in finite dimension. Stating (5.3) or (5.26) with Mathlib's sup norm would change the inequalities and is ruled out. Optimal values $\min(P)$, $\sup(D)$, $G$, $\Phi$, $\hat\alpha$ and $\ell$ are EReal infima and suprema, so that infeasibility gives $+\infty$ and an unbounded dual objective gives $-\infty$; real-valued sInf would silently return $0$ there. A solution of (D) must have $G(W^*)>-\infty$. The algorithm is a predicate on sequences, with $X^0$ arbitrary and $W^0\in\mathcal M$; Proposition 3.2 shows the predicate is satisfiable.

Subgradients and normal cones are those of convex analysis; the normal cone is the published definition FirstOrderOpt.ConvexTheory.normalCone. Not formalized: the nonconvex (Clarke) parts of §4, the linear-quadratic case and its rate result, and inexact minimization.

A complete development needs: separable minimization over product sets, finite-dimensional compactness arguments, convex duality for a perturbation function on a subspace, and the convergence of the proximal point algorithm for maximal monotone operators (Rockafellar 1976, Theorem 1 and Proposition 1). The last two are reusable well beyond this mission; contributions of either as standalone results are welcome.

Selected references

• R. T. Rockafellar and R. J.-B. Wets, Scenarios and policy aggregation in optimization under uncertainty, IIASA Working Paper WP-87-119, 1987; Mathematics of Operations Research 16(1):119–147, 1991. https://pure.iiasa.ac.at/id/eprint/2933/ — https://doi.org/10.1287/moor.16.1.119
• R. T. Rockafellar, Monotone operators and the proximal point algorithm, SIAM Journal on Control and Optimization 14(5):877–898, 1976. https://doi.org/10.1137/0314056
• G. J. Minty, Monotone (nonlinear) operators in Hilbert space, Duke Mathematical Journal 29(3):341–346, 1962. https://doi.org/10.1215/S0012-7094-62-02933-2
• R. T. Rockafellar, Convex Analysis, Princeton University Press, 1970. https://doi.org/10.1515/9781400873173

Motivation

A constrained system can be written as a generalized equation: a smooth or linear expression must balance the normal cone of a feasible set. This form includes variational inequalities and the first-order conditions of many optimization problems. It is useful to know whether a small change in the right-hand side determines one nearby solution and how far that solution moves. In the linear finite-dimensional case, Robinson's Theorem 3.1 gives matrix conditions that guarantee something stronger: a unique solution for every right-hand side, with one Lipschitz bound valid over the entire space. The theorem also identifies an exact matrix criterion when the feasible set is a product with a nonnegative orthant.

The question is relevant whenever some coordinates of a system are free and others are constrained by inequalities. Eliminating the free coordinates produces a smaller matrix, the Schur complement. The resulting condition can be checked on the reduced block, and in the orthant case it becomes the familiar P-matrix criterion from linear complementarity. Robinson states both claims in a single theorem and relates them to the behavior of a generalized equation rather than to an isolated matrix property Robinson, pp. 51–52.

Setting

Let $r,s$ be positive integers and give $\mathbb R^{r+s}$ its Euclidean norm. Write $w=(w_1,w_2)$, with $w_1\in\mathbb R^r$ and $w_2\in\mathbb R^s$. Let $K\subseteq\mathbb R^s$ be nonempty, closed, and convex, and put $D=\mathbb R^r\times K$. For $w\in D$, the normal cone is

$$N_D(w)=\{v\in\mathbb R^{r+s}:\langle v,c-w\rangle\leq0\text{ for every }c\in D\}.$$

For $w\notin D$, $N_D(w)$ is empty. Thus the inclusion $y\in Aw+N_D(w)$ itself enforces feasibility. This is Robinson's indicator-subdifferential convention Robinson, pp. 43 and 51.

Partition the real square matrix $A$ into blocks $A_{11}$, $A_{12}$, $A_{21}$, and $A_{22}$, where $A_{11}$ is $r\times r$. When $A_{11}$ is nonsingular, its Schur complement in $A$ is

$$A/A_{11}=A_{22}-A_{21}A_{11}^{-1}A_{12}.$$

A real matrix $M$ is positive definite in the theorem's sense when $v^\top Mv>0$ for every nonzero $v$. Symmetry is not assumed. A P-matrix has strictly positive determinants for every nonempty principal submatrix. Write $\mathbb R^s_+$ for the vectors with all components nonnegative. These conventions are stated or used in Theorem 3.1 and its proof.

Formalization targets

General closed convex set

For $T_K(w)=Aw+N_{\mathbb R^r\times K}(w)$, the first target is the sufficient condition

$$\det A_{11}\ne0\quad\text{and}\quad v^\top(A/A_{11})v>0\;(v\ne0) \quad\Longrightarrow\quad T_K^{-1}:\mathbb R^{r+s}\to\mathbb R^{r+s}\text{ is Lipschitz}.$$

Here “is Lipschitz” includes that exactly one inverse value exists at every $y$. The assertion holds for every nonempty closed convex $K$; it is not restricted to polyhedral sets.

Nonnegative orthant

When $K=\mathbb R^s_+$, the second target is an equivalence:

$$T_{\mathbb R^s_+}^{-1}\text{ is everywhere defined, single-valued and Lipschitz} \quad\Longleftrightarrow\quad \det A_{11}\ne0\text{ and }A/A_{11}\text{ is a P-matrix}.$$

The goal theorem is the conjunction of these two statements, matching Theorem 3.1. Its milestones record the paper's local upper Lipschitz assertion for polyhedral multifunctions, the block reduction, the reduced inverse result, the complementarity characterization, and the necessity of nonsingularity.

Significance

The first part gives a global sensitivity guarantee for a generalized equation with an arbitrary nonempty closed convex lower-block constraint. For any two right-hand sides, their unique solutions differ by at most a fixed multiple of the distance between those right-hand sides. The second part tells when this guarantee holds for the orthant model: the P-matrix condition is necessary as well as sufficient. These are the precise conclusions of Robinson's theorem; the mission concerns their machine-checked formalization, not an unresolved mathematical conjecture.

A complete development would provide reusable definitions of the normal-cone generalized equation, an explicitly nonsymmetric notion of positive definiteness, P-matrices, and the global inverse property. It would also formalize the linear complementarity characterization quoted in the paper. The current draft has Lean-checked statements with proof placeholders. No proof of Theorem 3.1 or of its milestones is claimed here.

Difficulty

Nonsingularity of the full matrix $A$ alone does not decide whether the constrained inverse exists uniquely: the normal cone can change as $w_2$ reaches different faces of $K$. A pointwise uniqueness argument also does not by itself supply one Lipschitz constant on all of $\mathbb R^{r+s}$. For the orthant case, checking only the full determinant or leading principal minors would miss the criterion Robinson invokes; every nonempty principal minor matters. The difficulty is to connect a global statement about all right-hand sides with the reduced matrix condition while preserving the normal cone's behavior at boundary points Robinson, pp. 51–52.

Formalization scope

Lean represents $\mathbb R^{r+s}$ as a Euclidean space indexed by the disjoint union of $\mathrm{Fin}(r)$ and $\mathrm{Fin}(s)$. This gives the paper's Euclidean norm rather than the sup norm of a product type. The normal cone is empty off its set, and the inverse property explicitly requires existence for every right-hand side, uniqueness, and one global real Lipschitz bound. The Schur complement appears only with a nonsingular $A_{11}$ in the substantive claims; matrix inversion on a singular input would otherwise be totalized by Lean. Positive definiteness does not add a symmetry hypothesis. The goal retains $r,s>0$ as printed; the paper's later remark about zero-dimensional blocks is outside this theorem.

The polyhedral milestones use finite intersections of closed half-spaces. The reduced-operator milestone uses a nonempty closed convex $K$, matching the theorem's setting. The LCP milestone allows any finite dimension; at dimension zero, both its universal uniqueness statement and its empty-minor P-matrix condition hold. A weaker encoding in which the inverse is merely single-valued where it happens to exist would erase the global conclusion and is excluded here.

Solvers can contribute the block-normal-cone equivalence, the global inverse statement for the reduced operator, the polyhedral Lipschitz result, the P-matrix characterization, and the final theorem. The normal-cone and matrix definitions can also serve later generalized-equation missions that use the same conventions.

Selected references

• S. M. Robinson, Strongly Regular Generalized Equations, Mathematics of Operations Research 5(1), 43–62, 1980. DOI 10.1287/moor.5.1.43.

Motivation

Many problems of optimization and equilibrium can be written as one inclusion. Some examples are the Karush–Kuhn–Tucker conditions of a nonlinear program, nonlinear complementarity problems, variational inequalities and the equilibrium conditions of economic and traffic models. In each of them a point $x$ has to satisfy $0 \in f(x) + \partial\psi_C(x)$. In applications the data $f$ depend on parameters: estimated coefficients, prices, or the iterate of an algorithm. One then wants to know whether a solution persists under a perturbation of the parameter, whether it stays locally unique, and how fast it moves. For a smooth equation $f(p, x) = 0$ the classical implicit-function theorem answers all three questions. Inequality constraints make the solution set nonsmooth, and the classical theorem no longer applies.

S. M. Robinson's paper Strongly regular generalized equations (Math. Oper. Res. 5 (1980) 43–62) introduced the condition of strong regularity and proved an implicit-function theorem for generalized equations under it. The condition and the theorem became the starting point of the stability theory of variational inequalities and nonlinear programs. Later work characterises strong regularity for nonlinear programs (Robinson's own §4, through the strong second-order sufficient condition), and develops it in the theory of metric regularity and of single-valued Lipschitz localisations (Dontchev–Rockafellar, Implicit Functions and Solution Mappings, Springer 2009/2014). Its sensitivity results underlie the convergence analysis of Newton and sequential quadratic programming methods for these problems.

This mission formalizes §1–§2 of the paper: the definition, the implicit-function theorem (Theorem 2.1) with its proof, and the three results derived from it.

Setting

Let $X$ be a real normed linear space and $X'$ its topological dual, the continuous linear functionals on $X$. Let $C \subseteq X$ be closed and convex. The normal-cone operator of $C$ is

$$\partial\psi_C(x) = \begin{cases} \{y \in X' : y(c - x) \le 0 \text{ for all } c \in C\}, & x \in C, \\ \emptyset, & x \notin C. \end{cases}$$

A generalized equation is the inclusion $0 \in f(x) + \partial\psi_C(x)$, for a function $f$ from an open set $\Omega \subseteq X$ to $X'$. For $C = X$ it is the equation $f(x) = 0$. For $C$ the nonnegative orthant of $\mathbb R^n$ it is the nonlinear complementarity problem. For a general closed convex $C$ it is the variational inequality over $C$.

Let $x_0 \in \Omega$ solve the generalized equation, and let $f$ be Fréchet differentiable at $x_0$. Linearising $f$ gives the multifunction

$$T x = f(x_0) + f'(x_0)(x - x_0) + \partial\psi_C(x).$$

The generalized equation is strongly regular at $x_0$ with associated Lipschitz constant $\lambda$ if there are neighbourhoods $U$ of $0$ in $X'$ and $V$ of $x_0$ with the following property: for every $y \in U$ the inclusion $y \in Tx$ has exactly one solution $s(y)$ in $V$, and $\|s(y_1) - s(y_2)\| \le \lambda \|y_1 - y_2\|$ on $U$. The condition depends only on $f(x_0)$, $f'(x_0)$ and $C$.

In the parametric problem the function is $f(p, x)$, with $p$ in a topological space $P$ and base point $p_0$. Its partial derivative in $x$ is $f'(p, x)$.

Formalization targets

Goal: Theorem 2.1 (p. 45)

Assume that $f'$ exists on $P \times \Omega$, that $f$ and $f'$ are continuous at $(p_0, x_0)$, and that $x_0$ solves $0 \in f(p_0, x) + \partial\psi_C(x)$, strongly regularly with constant $\lambda$. Then for every $\varepsilon > 0$ there are neighbourhoods $N_\varepsilon$ of $p_0$ and $W_\varepsilon$ of $x_0$, and a function $x : N_\varepsilon \to W_\varepsilon$, with two properties. First, $x(p)$ is the unique solution in $W_\varepsilon$ of $0 \in f(p, x) + \partial\psi_C(x)$. Second,

$$\|x(p) - x(q)\| \le (\lambda + \varepsilon)\,\|f(p, x(q)) - f(q, x(q))\| \qquad (p, q \in N_\varepsilon).$$

Milestones: the steps of the proof (p. 46)

Write $L$ for the linearisation at $(p_0, x_0)$, $r(p, x) = f(p_0, x_0) + f'(p_0, x_0)(x - x_0) - f(p, x)$ for the residual, and $\Phi_p(x) = V \cap L^{-1}[r(p, x)]$. The milestones are four claims of the proof:

1. On the closed ball $V_\varepsilon$, the fixed points of $\Phi_p$ are exactly the solutions of the perturbed inclusion.
2. $\|\Phi_p(x_1) - \Phi_p(x_2)\| \le \lambda\delta \|x_1 - x_2\|$ on $V_\varepsilon$.
3. $\Phi_p$ maps $V_\varepsilon$ into itself.
4. The contraction principle on $V_\varepsilon$, with the bound $\|x(p) - x\| \le (1 - \lambda\delta)^{-1} \|\Phi_p(x) - x\|$ (2.5).

Companion results

• Corollary 2.2 (pp. 46–47). If $P$ lies in a normed space and $\|f(p, x) - f(q, x)\| \le \nu \|p - q\|$, then $x(\cdot)$ is Lipschitzian with modulus $\nu(\lambda + \varepsilon)$.
• Theorem 2.3 (p. 47). Let $\Phi_p(x_0)$ be the unique local solution of the linear generalized equation $0 \in f(p, x_0) + f'(p_0, x_0)(x - x_0) + \partial\psi_C(x)$. Then $\|x(p) - \Phi_p(x_0)\| \le \alpha_\varepsilon(p)\|p - p_0\|$, with $\alpha_\varepsilon(p) \to 0$.
• Theorem 2.4 (p. 48). This is a perturbation lemma for linear generalized equations $0 \in Ax + a + \partial\psi_C(x)$. Near $(A_0, a_0)$ the localised inverse stays single-valued, and it is Lipschitzian with modulus $\lambda(1 - \lambda\|A - A_0\|)^{-1}$.
• The case $C = X$ (remark on p. 45). There strong regularity is equivalent to $f'(x_0)$ having a continuous linear inverse.
• Strong regularity is the weakest possible condition (remark on p. 47). For the canonical perturbation $f(p,x) = f(x_0) + f'(x_0)(x - x_0) - p$, $p \in X'$, locally unique solvability with a Lipschitzian solution map $x(\cdot)$ near $p_0 = 0$ is exactly strong regularity at $x_0$.

Significance

Theorem 2.1 reduces the stability of a nonsmooth problem to a property of one linearised problem at the base point. Once strong regularity is checked at $x_0$, the solution exists, is locally unique and moves in a Lipschitz manner for every nearby parameter. Corollary 2.2 and Theorem 2.3 make this quantitative. Theorem 2.3 shows that a linear complementarity problem or a quadratic program approximates the nonlinear solution to first order. Theorem 2.4 is the analogue of Banach's perturbation lemma, which keeps an operator invertible under small perturbations. The later sections of the paper give checkable criteria for strong regularity: a Schur-complement test, a reduced-form test over polyhedral sets, and second-order conditions for nonlinear programs. Those criteria are the subject of the companion missions II–IV.

The theorems have been proved since 1980 and have been restated in textbooks. No machine-checked version of strong regularity or of this implicit-function theorem is known. Mathlib has the classical implicit- and inverse-function theorems for $C^1$ maps between Banach spaces, and the contraction mapping principle. It has nothing on normal-cone operators, generalized equations or Lipschitz localisations. The formalization adds this layer. It also shows that Theorem 2.1 needs no completeness of $X$.

Difficulty

The obvious approach is to apply the classical implicit-function theorem. It fails because $x \mapsto f(x) + \partial\psi_C(x)$ is set-valued and is not differentiable in any useful sense. The inverse is localised only to a neighbourhood $V$ and is single-valued only there. So every step must keep the iterates inside sets on which strong regularity says something. The neighbourhoods have to be chosen in the right order: $\delta$ from $\varepsilon$, then $U$ and $V$, then the ball radius $\rho$, then the parameter neighbourhood. The estimate (2.4) has to come out with the constant $\lambda + \varepsilon$ and no larger. The page applies the contraction principle in a space it calls only normed. A faithful proof must also obtain the fixed point without completeness of $X$, for example through the completeness of $X'$.

Formalization scope

• $X'$ is StrongDual ℝ X, the normal cone is a Set (StrongDual ℝ X) that is empty off $C$, and an inclusion $0 \in f(x) + \partial\psi_C(x)$ is written $-f(x) \in \partial\psi_C(x)$.
• Strong regularity is a predicate on $(C, f(x_0), f'(x_0), x_0, \lambda)$ with real $\lambda$. It requires existence, uniqueness in $V$ and the Lipschitz bound.
• $f$ is a total function $P \times X \to X'$ whose values off $P \times \Omega$ carry no hypothesis. The conclusion therefore asserts $W_\varepsilon \subseteq \Omega$.
• Differentiability is HasFDerivAt at every point of $\Omega$, for every $p$. Continuity of $f'$ at $(p_0, x_0)$ is joint and in operator norm.
• Uniqueness is uniqueness within $W_\varepsilon$, never global. Global uniqueness is false in general. Dropping uniqueness, or quantifying $\varepsilon$ after the neighbourhoods, would trivialize the goal, and the statements are written to exclude both.
• Theorem 2.1, Corollary 2.2 and Theorem 2.3 do not assume completeness, as printed. Theorem 2.4 assumes a Banach space, as printed.
• The contraction-principle milestone assumes completeness, under which the principle as invoked holds.
• Corollary 2.2 asks for its Lipschitz hypothesis on $N_\varepsilon \times V_\varepsilon$, sets that Theorem 2.1 produces. It is therefore stated on fixed neighbourhoods given in advance.
• Theorem 2.4 also asserts $\lambda\|A - A_0\| < 1$, which its proof secures.

A complete development needs the following:

• the mean-value inequality for Fréchet-differentiable maps on convex sets;
• the contraction principle on a closed ball;
• a fixed-point argument that obtains the fixed point through the complete space $X'$.

The definitions of normal cone and strong regularity are reusable for variational inequalities in general. Contributions are welcome at any level: proofs of the milestones, the goal from the milestones, or the companion results from the goal.

Selected references

• S. M. Robinson, Strongly regular generalized equations, Mathematics of Operations Research 5(1):43–62, 1980. https://doi.org/10.1287/moor.5.1.43
• A. L. Dontchev and R. T. Rockafellar, Implicit Functions and Solution Mappings, 2nd ed., Springer, 2014. https://doi.org/10.1007/978-1-4939-1037-3

Motivation

Maximum likelihood estimates, least-squares fits and sample-average approximations of stochastic programs solve $0 = \bar f^\nu(x)$, where $\bar f^\nu$ averages a random integrand over $\nu$ observations. Their classical asymptotic theory rests on the implicit function theorem and needs a smooth, unconstrained problem.

When the estimate is constrained to a set, for example a nonnegativity constraint, a simplex or a polyhedron, the first-order conditions become a generalized equation

$$0 \in f(z, x) + N(x),$$

where $N$ is a multifunction such as the normal cone of the constraint set. The same form describes optimality conditions of stochastic programs and variational inequalities. Aitchison and Silvey (1958) treated equality-constrained maximum likelihood. Huber (1967) allowed nonsmooth estimating functions but required an open parameter domain. Dupačová and Wets (1988) and Shapiro (1989) derived limit laws for solutions of stochastic programs under smoothness of the expected gradient. King and Rockafellar (1993) gave a general theory that needs neither smoothness of the expected map nor single-valuedness of $N$: the limit law of the normalized error is the image of a Gaussian under a contingent derivative, a positively homogeneous and generally nonlinear map, so the limit is generally not normal.

Setting

Let $Z$ be a separable Banach space with norm $\|\cdot\|$, and let $|\cdot|$ be the Euclidean norm on $\mathbb R^n$ and $\mathbb R^m$. A multifunction $G : Z \rightrightarrows \mathbb R^n$ assigns a set $G(z) \subseteq \mathbb R^n$ to each $z$. Its graph is $\operatorname{gph} G$ and its inverse is $G^{-1}(x) = \{z \mid x \in G(z)\}$.

For sets $A_t$ indexed by $t \downarrow 0$, the upper limit $\limsup A_t$ consists of the points $x$ with $x = \lim x_k$, $x_k \in A_{t_k}$ for some $t_k \downarrow 0$, and the lower limit of the points reachable along every such sequence. The contingent derivative of $G$ at $(z, x) \in \operatorname{gph} G$ is the multifunction $DG(z|x)$ with

$$\operatorname{gph} DG(z|x) = \limsup_{t \downarrow 0} t^{-1}\big[\operatorname{gph} G - (z,x)\big].$$

$G$ is proto-differentiable when this upper limit equals the lower limit, and semi-differentiable when $t^{-1}[G(z + t w') - x] \to DG(z|x)(w)$ as $t \downarrow 0$ and $w' \to w$. A single-valued $g$ is B-differentiable at $z$ when $t^{-1}[g(z + tw') - g(z)] \to Dg(z)(w)$ in the same sense.

The deterministic problem is $0 \in f(z, x) + N(x)$ with $f : Z \times \mathbb R^n \to \mathbb R^m$, data $z$ and solution map $J(z)$. At a reference pair $(z^*, x^*)$ set $F = f(z^*, \cdot) + N$. The analytical assumptions M.1–M.4 are as follows. $f$ is jointly continuous and B-differentiable in each variable, with the $z$-derivative $D_z f(z^*, x^*)$ strong (uniform in $x$ near $x^*$). $N$ is closed and proto-differentiable. $F$ is subinvertible: $0 \in F(x^*)$, and a closed-graph, convex-valued selection of $F^{-1}$ near $0$ passes through $x^*$. The contingent derivative $DF^{-1}(0|x^*)$ is at most single-valued.

The statistical problem has i.i.d. random elements $s_1, s_2, \dots$ of a measurable space $S$ and an integrand $f : U \times S \to \mathbb R^m$ on a compact neighborhood $U$ of $x^*$. It satisfies the probabilistic assumptions P.1–P.4: continuity in $x$, measurability in $s$, a finite second moment at one point, and a Lipschitz bound $|f(x_1,s) - f(x_2,s)| \le a(s)|x_1 - x_2|$ with $E a(s_1)^2 < \infty$. The M-estimate $x^\nu$ is a measurable solution of

$$0 \in \bar f^\nu(x) + N(x), \qquad \bar f^\nu(x) = \frac1\nu\sum_{i=1}^\nu f(x, s_i),$$

and the true equation is $0 \in Ef(x) + N(x)$ with $F = Ef + N$.

Formalization targets

Goal: Theorem 2.7 (asymptotic distribution of M-estimates)

Under P.1–P.4 on a compact neighborhood $U$ of $x^*$, B-differentiability of $Ef$ at $x^*$, and M.2–M.4 for $F = Ef + N$, every sequence of measurable solutions $x^\nu$ of (2.5) with $x^\nu \to x^*$ almost surely satisfies

$$\sqrt\nu\,[x^\nu - x^*] \xrightarrow{\ \mathcal D\ } DF^{-1}(0|x^*)(-w^*), \qquad w^* \sim \mathcal N\big(0, \operatorname{cov} f(x^*, s_1)\big).$$

The goal fixes the limit law completely: the map is the contingent derivative of $F^{-1}$, and the Gaussian has the covariance of the integrand at $x^*$.

Milestones

• Theorem 2.4 gives bounds in probability, $P\{|x^\nu - x^*| > \delta\} \le P\{\alpha\lambda\|z^\nu - z^*\| > \delta\}$. Its proof uses the upper-Lipschitz property $U \cap F^{-1}(y) \subseteq x^* + \lambda|y|B$ (a display of the proof).
• Theorem 2.6 is the abstract limit theorem: if $\tau_\nu^{-1}[z^\nu - z^*] \to_{\mathcal D} w$, then $\tau_\nu^{-1}[x^\nu - x^*] \to_{\mathcal D} DF^{-1}(0|x^*)(-D_z f(z^*,x^*)(w))$. Its proof uses two displays: semi-differentiability of the localized solution map, with $DJ(z^*|x^*)(w) = DF^{-1}(0|x^*)(-D_z f(z^*,x^*)(w))$, and a Lipschitz bound $|x - x^*| \le \lambda\|z - z^*\|$ on $U \cap J(z)$.
• Proposition A1, Corollary A2 and Theorem A3 concern the space $C_m(U)$ under P.1–P.4. The integrand and the empirical means are random elements of $C_m(U)$, and $\sqrt\nu(\bar f^\nu - Ef)$ converges in distribution to a Gaussian element of $C_m(U)$.

The three displays (Theorem 2.4's upper-Lipschitz inclusion, and Theorem 2.6's semi-differentiability and Lipschitz bound) are statements the paper cites from King and Rockafellar, Sensitivity analysis for nonsmooth generalized equations ([12]: Proposition 2.1, Theorem 4.1, Remark 4.3). They are cited results, not this paper's own, and are milestones because the proofs of Theorems 2.4 and 2.6 rest on them.

Significance

Theorem 2.7 gives the limit law of constrained and nonsmooth M-estimates in a form that can be computed. When $N$ is the normal cone of a polyhedron, $DF^{-1}(0|x^*)$ is piecewise linear, and the limit is the solution of a random linear complementarity or quadratic problem driven by a Gaussian vector. This underlies the asymptotic theory of sample-average approximation in stochastic programming, where the paper applies it to stochastic programs (§3) and to piecewise linear-quadratic tracking problems (§4). Theorem 2.6 separates the deterministic sensitivity analysis from the probability, so any data sequence with a known limit law yields a limit law for the solutions.

The result is proved in the paper modulo the cited theorems of [12] and [11], but none of it is machine-checked. Mathlib has the real-valued i.i.d. central limit theorem, Gaussian measures on Banach spaces and convergence in distribution. It has no multivariate or Banach-space central limit theorem, no contingent derivatives and no set-valued implicit function theorem. Formalizing the mission produces these, along with a checked version of the cited sensitivity results.

Difficulty

The classical argument linearizes $F$ at $x^*$, inverts the Jacobian and applies the delta method. Here $F$ is set-valued and its derivative is only positively homogeneous. There is no Jacobian to invert, and the solution map need not be differentiable or even single-valued away from $x^*$. The replacement for the implicit function theorem is the semi-differentiability of the localized solution map under M.1–M.4. Proving it means controlling both the upper and the lower set limits of difference quotients of solution sets, and subinvertibility is what supplies existence of nearby solutions.

The probabilistic side cannot work coordinate by coordinate either. The estimate solves an equation in the whole function $\bar f^\nu$, so convergence of $\bar f^\nu$ at finitely many points is not enough. The central limit theorem must hold in the sup norm on $C_m(U)$, which requires tightness of the empirical process, and only then can the deterministic sensitivity result be composed with it.

Formalization scope

Points live in EuclideanSpace ℝ (Fin n), and $Z$ is a real normed space with [CompleteSpace Z] [SeparableSpace Z] where the paper says "separable Banach". Set limits are Kuratowski limits along filters ($t \downarrow 0$ is 𝓝[>] 0, and $(t,w') \to (0^+,w)$ is the product filter). The contingent derivative is defined by (2.2) alone. M.4's printed sum formula equals it under M.1, and this is not assumed. "B-differentiable" is read as the limit (2.4). Products carry Lean's max norm. $s_1$ is s 0, empirical means sum over Finset.range ν, and (2.5) is required for $\nu \ge 1$. $F = Ef + N$ is empty off $U$. Convergence in distribution is Mathlib's TendstoInDistribution, with the limit on its own probability space. The law of $w^*$ is fixed through linear functionals: $\langle\ell, w^*\rangle \sim \mathcal N(0, \operatorname{Var}\langle\ell, f(x^*,s_1)\rangle)$. In Appendix A1–A3 the integrand is S → C(↥U, Rn m) with the Borel σ-algebra, and "Gaussian" is IsGaussian.

Explicit choices relative to the printed text:

1. The paper states that an almost surely convergent sequence of solutions "converges to the point $x^*$" (Theorems 2.6 and 2.7). This is false when the true equation has a second solution: $f(z,x) = x^2 - x - z$, $N \equiv \{0\}$, $z^* = x^* = 0$ satisfies M.1–M.4 with $J(0) = \{0, 1\}$. The formalization assumes $x^\nu \to x^*$ almost surely instead.
2. $0 \in F(x^*)$ (presupposed by $DF^{-1}(0|x^*)$) is explicit in Theorem 2.4 and the upper-Lipschitz display.
3. The threshold for "all sufficiently small $\delta$" in Theorem 2.4 is chosen with $U$ and $\lambda$, before the random elements.
4. Proposition A1 and Corollary A2 carry P.1–P.4, as stated or inherited on the Appendix page, though their measurability conclusions use only P.1.
5. In the limit theorems, the single-valued map $DF^{-1}(0|x^*)$ is a function $L$ whose values lie in the contingent derivative at every point.

The conclusion of the goal names its limit: the image of the stated Gaussian under a map $L$ whose values lie in the contingent derivative of $F^{-1}$. A statement asserting only that $\sqrt\nu(x^\nu - x^*)$ converges in distribution to some limit, or replacing $DF^{-1}(0|x^*)$ by a linear map, is a different and weaker theorem. A sorry-free check confirms that the goal's hypotheses can all be met (a degenerate instance with $f(x,s) = x$, $N \equiv \{0\}$).

A complete development needs Kuratowski set convergence and contingent derivatives (reusable across set-valued analysis), a central limit theorem in $C(K)$ for Lipschitz-indexed processes (reusable for empirical-process theory), and a continuous-mapping argument for random closed sets. Contributions to any of these layers are welcome, as are proofs of the cited [12] statements.

Selected references

• A. J. King and R. T. Rockafellar, Asymptotic theory for solutions in statistical estimation and stochastic programming, Mathematics of Operations Research 18(1) (1993). https://doi.org/10.1287/moor.18.1.148
• A. J. King and R. T. Rockafellar, Sensitivity analysis for nonsmooth generalized equations, Mathematical Programming 55 (1992) 193–212. https://doi.org/10.1007/BF01581199
• A. J. King, Generalized delta theorems for multivalued mappings and measurable selections, Mathematics of Operations Research 14(4) (1989) 720–736. https://doi.org/10.1287/moor.14.4.720
• J. Dupačová and R. J.-B. Wets, Asymptotic behavior of statistical estimators and of optimal solutions of stochastic optimization problems, Annals of Statistics 16(4) (1988) 1517–1549. https://doi.org/10.1214/aos/1176351052
• A. Shapiro, Asymptotic properties of statistical estimators in stochastic programming, Annals of Statistics 17(2) (1989) 841–858. https://doi.org/10.1214/aos/1176347146
• P. J. Huber, The behavior of maximum likelihood estimates under nonstandard conditions, Proc. Fifth Berkeley Symp. Math. Statist. Probab. 1 (1967) 221–233. https://projecteuclid.org/euclid.bsmsp/1200512988
• A. Araujo and E. Giné, The Central Limit Theorem for Real and Banach Valued Random Variables, Wiley, 1980.

Motivation

Many problems in statistics, engineering and operations research ask for the input level $x$ at which an unknown response function is largest, when the response can only be measured with noise: the dose that maximizes a yield, the setting that maximizes throughput in a simulation, the parameter that maximizes an expected reward. Robbins and Monro (1951) had shown how to find, by sequential noisy measurements, the root of an unknown increasing function. Kiefer and Wolfowitz (1952) adapted the idea to maximization: at each step, measure the response at two points $z_n\pm c_n$ placed symmetrically about the current estimate, and move the estimate in the direction of the observed difference. Their scheme is the origin of finite-difference stochastic approximation, the gradient-free ancestor of simultaneous-perturbation and zeroth-order methods in simulation optimization and machine learning.

Timeline. 1951: Robbins and Monro, root finding, convergence in mean square. 1952: Wolfowitz, convergence in probability for Robbins–Monro under weaker conditions; Kiefer and Wolfowitz, the present paper, convergence in probability of the maximization scheme. 1954: Blum proved almost-sure convergence of both schemes under related conditions (doi:10.1214/aoms/1177728794). Later work (Spall 1992, among many) extended the finite-difference idea to many dimensions.

Setting

For each level $x\in\mathbb R$, an observation taken at $x$ is a random number with distribution $H(\cdot\mid x)$; the family depends measurably on $x$. The regression function is the mean observation,

$$M(x)=\int_{-\infty}^{\infty}y\,dH(y\mid x),$$

and the observations have uniformly bounded variance: $\int (y-M(x))^2\,dH(y\mid x)\le S$ for all $x$. The function $M$ is unimodal about an unknown point $\theta$: strictly increasing for $x<\theta$ and strictly decreasing for $x>\theta$.

Two sequences of positive numbers $a_n$ (step sizes) and $c_n$ (probe widths) satisfy

$$c_n\to0,\qquad \sum a_n=\infty,\qquad \sum a_nc_n<\infty,\qquad \sum a_n^2c_n^{-2}<\infty,$$

for example $a_n=n^{-1}$, $c_n=n^{-1/3}$. Starting from an arbitrary number $z_1$, the Kiefer–Wolfowitz process is

$$z_{n+1}=z_n+a_n\,\frac{y_{2n}-y_{2n-1}}{c_n},$$

where, given everything observed so far, $y_{2n-1}$ and $y_{2n}$ are independent with laws $H(\cdot\mid z_n-c_n)$ and $H(\cdot\mid z_n+c_n)$.

The paper imposes three regularity conditions on $M$:

1. Condition 1 (Lipschitz near $\theta$): for some $\beta,B>0$, $|x'-\theta|+|x''-\theta|<\beta$ and $x'\neq x''$ imply $|M(x')-M(x'')|<B|x'-x''|$.
2. Condition 2 (bounded oscillation): for some $\rho,R>0$, $|x'-x''|<\rho$ implies $|M(x')-M(x'')|<R$.
3. Condition 3 (no flat regions away from $\theta$): for every $\delta>0$ there is $\pi(\delta)>0$ such that $|z-\theta|>\delta$ implies $|M(z+\varepsilon)-M(z-\varepsilon)|/\varepsilon>\pi(\delta)$ for all $0<\varepsilon<\delta/2$.

In the Lean development these objects are regFun H ($M$), SecondMomentBound H S, Unimodal, Cond1, Cond2, Cond3, StepSizes a c and the structure IsKWProcess, all in the namespace KieferWolfowitz.Convergence.

Formalization targets

Goal: convergence in probability

$$\text{for every }\eta>0,\qquad P\{|z_n-\theta|>\eta\}\longrightarrow 0\quad(n\to\infty).$$

This is the theorem of the paper (stated on p. 462, proved in §3 in the form (3.22)). It fixes no rate and no constants, so it is the stable form of the result.

Milestones, in the order of the paper's proof

1. (3.6) The one-step identity for $b_n=E(z_n-\theta)^2$: $b_{n+1}=b_n+2\frac{a_n}{c_n}E\,U_n(z_n)+\frac{a_n^2}{c_n^2}E(y_{2n}-y_{2n-1})^2$, with $U_n(z)=(z-\theta)(M(z+c_n)-M(z-c_n))$.
2. (3.8) $0\le U_n^+(z)<2Bc_n^2$ whenever $c_n<\frac12\beta$.
3. (3.11)–(3.12) $E\{(y_{2n}-y_{2n-1})^2\mid z_n\}\le 2S+R^2$ whenever $c_n<\rho/2$.
4. (3.15) $\sum (a_n/c_n)N_n$ converges, where $N_n=E\,U_n^-(z_n)$.
5. (3.18) $\liminf_n E\{K_n|z_n-\theta|\}=0$, with $K_n=|(M(z_n+c_n)-M(z_n-c_n))/c_n|$.
6. After (3.19) Along any $n_1<n_2<\cdots$ with $E\{K_{n_j}|z_{n_j}-\theta|\}\to0$, $z_{n_j}\to\theta$ in probability.
7. (3.28) If $N_0$ satisfies (3.25)–(3.27) for $s>0$, then $E\{(z_n-\theta)^2\mid z_{N_0}\}<(z_{N_0}-\theta)^2+s$ for all $n>N_0$.

Significance

The result shows that the maximizer of a function observable only through noise can be located consistently using function evaluations alone, without gradients and without a parametric model of $M$. The conditions are qualitative: no differentiability, no specific noise distribution, and the noise may change shape with $x$. This is the template for gradient-free stochastic optimization: simulation optimization with finite-difference gradient estimates, simultaneous-perturbation methods, and zeroth-order methods in learning are all analysed by refining the same bias–variance balance between the probe width $c_n$ and the step $a_n$.

The theorem has been proved since 1952 and is not open. To our knowledge no machine-checked proof of it exists. A formalization produces a checked proof of a founding result of stochastic approximation and, along the way, the probabilistic infrastructure such proofs share: processes defined through conditional laws given a filtration, recursions for second moments, conditional Chebyshev arguments, and the passage from subsequence convergence to full convergence by a restart bound. Alternative proofs (via supermartingale convergence, which would give almost-sure convergence under stronger hypotheses) are welcome as separate developments but are not the goal.

Difficulty

The obvious approach is to show that $E(z_n-\theta)^2$ decreases. It does not: the drift term $2\frac{a_n}{c_n}E\,U_n(z_n)$ is negative only on average and only once the probes straddle $\theta$ correctly, while the noise term $\frac{a_n^2}{c_n^2}e_n$ is always positive and the finite-difference estimate of the slope has a bias of order $c_n$. A one-step contraction is therefore unavailable, and summability arguments alone give smallness only along a subsequence; turning that into convergence of the whole sequence is where the difficulty lies, and it is also where Condition 3, which forbids flat stretches of $M$ away from $\theta$, cannot be dispensed with. The argument needs integrability of quantities the paper treats as finite without comment, and the conditional statements require the conditional law of the observation pair, not just its unconditional moments.

Formalization scope

• Observations. $H$ is a Markov kernel Kernel ℝ ℝ. Every $H(\cdot\mid x)$ has $\int y^2\,dH<\infty$, so $M$ and the variance in (2.2) are genuine integrals.
• The process. On a probability space with a filtration $(\mathcal F_n)$, $z_n$ is $\mathcal F_n$-measurable, the pair $(y_{2n-1},y_{2n})$ is $\mathcal F_{n+1}$-measurable, and $P(y_{2n-1}\in A, y_{2n}\in B\mid\mathcal F_n)=H(A\mid z_n-c_n)H(B\mid z_n+c_n)$ a.s. for all Borel $A,B$. This is the reading of "independent chance variables with respective distributions $H(y\mid z_n-c_n)$ and $H(y\mid z_n+c_n)$". The starting point $z_1$ is a deterministic number.
• Indices are 0-based: Lean's z 0 is $z_1$, a n, c n are $a_{n+1}$, $c_{n+1}$.
• Condition 1 as printed has a strict inequality that fails at $x'=x''$, which makes it unsatisfiable; the formalization adds $x'\neq x''$. Condition 3's infimum over $\frac12\delta>\varepsilon>0$ is stated pointwise, which is equivalent since $\pi(\delta)$ is existential, with $\pi(\delta)$ uniform in $z$.
• Explicit readings of loose phrases. "For $n$ sufficiently large" in (3.10)–(3.12) is the hypothesis $c_n<\rho/2$; conditioning "on $z_n$" or "on $z_{N_0}=z$" is conditioning on $\mathcal F_n$ or $\mathcal F_{N_0}$; "converges stochastically" is Mathlib's TendstoInMeasure; "$\liminf=0$" for a nonnegative sequence is "below every $\varepsilon>0$ infinitely often".
• Integrability of every expectation in a milestone is either a hypothesis or part of the conclusion, so no statement holds by Lean's convention that the integral of a non-integrable function is $0$.
• Ruled out. A formalization in which the observations are arbitrary random variables not tied to $H$, in which $\sum a_n=\infty$ is dropped, or in which $\pi(\delta)$ may depend on $z$, would make the goal false or a different theorem; the hypotheses here are jointly satisfiable (checked for $H(\cdot\mid x)=\delta_{-|x|}$, $\theta=0$, $a_n=n^{-1}$, $c_n=n^{-1/3}$).
• Not formalized. The remark on p. 463 that truncating $z_n\pm c_n$ to an interval $[C_1,C_2]$ preserves the conclusion (asserted without proof), the motivational discussion (a)–(c), and §4's further problems.
• Distinct from the Kiefer–Wolfowitz equivalence theorem of optimal experimental design (1960), which is already on the platform under BanditAlgorithm.kiefer_wolfowitz_equivalence; the two results share authors and nothing else.

Useful infrastructure, reusable beyond this mission: conditional expectations of functions of a pair with a given conditional law, tail-sum bounds, and the passage from a restart bound plus subsequence convergence in probability to full convergence. Contributions of any milestone, of intermediate lemmas such as (3.7), (3.9), (3.13) and (3.29), and of the final assembly are welcome.

Selected references

• J. Kiefer and J. Wolfowitz, Stochastic estimation of the maximum of a regression function, Ann. Math. Statist. 23 (1952), 462–466. https://doi.org/10.1214/aoms/1177729392
• H. Robbins and S. Monro, A stochastic approximation method, Ann. Math. Statist. 22 (1951), 400–407. https://doi.org/10.1214/aoms/1177729586
• J. Wolfowitz, On the stochastic approximation method of Robbins and Monro, Ann. Math. Statist. 23 (1952), 457–461. https://doi.org/10.1214/aoms/1177729391
• J. R. Blum, Approximation methods which converge with probability one, Ann. Math. Statist. 25 (1954), 382–386. https://doi.org/10.1214/aoms/1177728794
• J. C. Spall, Multivariate stochastic approximation using a simultaneous perturbation gradient approximation, IEEE Trans. Automat. Control 37 (1992), 332–341. https://doi.org/10.1109/9.119632

Motivation

Newton's method for a nonlinear system $F(x) = 0$, with $F:\mathbb R^n\to\mathbb R^n$, solves the linear system $F'(x_k)s_k = -F(x_k)$ at every step. For large systems that solve is itself iterative (a Krylov method such as GMRES), and it is stopped early. The resulting inexact Newton methods, introduced by Dembo, Eisenstat and Steihaug (SIAM J. Numer. Anal. 19 (1982)), accept any step with $\|F(x_k)+F'(x_k)s_k\|\le\eta_k\|F(x_k)\|$, where the forcing term $\eta_k\in[0,1)$ controls how accurately the linear system is solved. Their theory is local: it applies once the iterates are near a solution with invertible derivative.

Practical solvers (Newton–Krylov codes in large-scale simulation, nonlinear solver libraries such as PETSc's SNES and SUNDIALS' KINSOL) combine such inexact steps with a globalization, most often backtracking along the step. Eisenstat and Walker (SIAM J. Optim. 4 (1994)) gave the global convergence theory for this combination: what can be said about the iterates from an arbitrary starting point, with no assumption that a solution exists or that $F'$ is invertible anywhere.

Timeline. Dembo, Eisenstat and Steihaug (1982) proved local convergence of inexact Newton methods. Dembo and Steihaug (Math. Program. 26 (1983)) studied truncated Newton methods for unconstrained minimization. Brown and Saad (1990) studied globalized Newton–Krylov methods with line searches and model trust regions, under the inner-product norm. Eisenstat and Walker (1994) gave the general framework treated here, for an arbitrary norm. Their 1996 paper (SIAM J. Sci. Comput. 17) proposed the forcing-term choices that are now standard.

Setting

Let $E$ be $\mathbb R^n$ with an arbitrary norm $\|\cdot\|$, and let $F:E\to E$ be continuously differentiable with derivative $F'(x)$. A point $x_*$ is a limit point of $(x_k)$ if every ball $N_\delta(x_*)=\{y:\|y-x_*\|<\delta\}$ contains $x_k$ for infinitely many $k$.

Algorithm GIN (global inexact Newton method). Fix $t\in(0,1)$. At each $k$, find a level $\eta_k\in[0,1)$ and a step $s_k$ with

$$\|F(x_k)+F'(x_k)s_k\|\le\eta_k\|F(x_k)\| \quad (2.1),\qquad \|F(x_k+s_k)\|\le[1-t(1-\eta_k)]\,\|F(x_k)\| \quad (2.2),$$

and set $x_{k+1}=x_k+s_k$. Condition (2.1) says that $s_k$ reduces the norm of the local linear model by the factor $\eta_k$. Condition (2.2) asks that $\|F\|$ itself decrease by a fixed fraction $t$ of that predicted reduction.

Algorithm MR (minimum reduction method). Fix $\eta_{\max}\in[0,1)$ and $0<\theta_{\min}<\theta_{\max}<1$. At step $k$, choose $\bar\eta_k\in[0,\eta_{\max}]$ and a curve $\sigma_k$ with $\|F(x_k)+F'(x_k)\sigma_k(\eta)\|\le\eta\|F(x_k)\|$ for $\bar\eta_k\le\eta\le1$ (5.1). Start at $\eta_k=\bar\eta_k$. While (2.2) fails for $s_k=\sigma_k(\eta_k)$, replace $\eta_k$ by $1-\theta(1-\eta_k)$ for some $\theta\in[\theta_{\min},\theta_{\max}]$. Then set $x_{k+1}=x_k+\sigma_k(\eta_k)$.

Algorithm INB (inexact Newton backtracking). Choose $\bar\eta_k\in[0,\eta_{\max}]$ and an inexact Newton step $\bar s_k$ at level $\bar\eta_k$. While (2.2) fails, shorten the step, $s_k\leftarrow\theta s_k$, and raise the level, $\eta_k\leftarrow1-\theta(1-\eta_k)$. INB is MR with the backtracking curve $\sigma_k(\eta)=\frac{1-\eta}{1-\bar\eta_k}\bar s_k$ (6.1).

An algorithm does not break down if it produces an infinite sequence of iterates, in particular if every while-loop exits.

Formalization targets

Goal: Theorem 6.1 (global convergence of Algorithm INB)

If Algorithm INB does not break down and $x_*$ is a limit point of $(x_k)$ at which $F'(x_*)$ is invertible, then

$$F(x_*)=0,\qquad x_k\to x_*,\qquad s_k=\bar s_k\ \text{and}\ \eta_k=\bar\eta_k\ \text{for all sufficiently large }k.$$

The statement assumes no solution, bounded level set, Lipschitz derivative or particular norm. The last clause says that backtracking eventually stops, so the local rate is governed by the forcing terms $\bar\eta_k$.

Milestones

In attack order:

1. Lemmas 1.1 and 1.2. Continuity of $y\mapsto F'(y)^{-1}$ at an invertible point, and a uniform linearization error $\|F(z)-F(y)-F'(y)(z-y)\|\le\varepsilon\|z-y\|$ near $x$.
2. Theorem 3.3. If $F(x_k)\to0$, the steps satisfy (2.1) with a fixed $\eta$, and $\|F(x_k)\|$ is nonincreasing, then an invertible limit point is a zero and the limit.
3. Theorem 3.4. A GIN run with $\sum_k(1-\eta_k)=\infty$ has $F(x_k)\to0$, plus the conclusion of Theorem 3.3 at invertible limit points.
4. Theorem 3.5. A GIN run converges to a limit point near which $\|s_k\|\le\Gamma(1-\eta_k)\|F(x_k)\|$ (3.2).
5. Lemma 5.1. The while-loop terminates, with $1-\eta_k\ge\min\{1-\bar\eta_k,\theta_{\min}\delta/(\Gamma\|F(x_k)\|)\}$.
6. MR runs are GIN runs (§5).
7. Theorem 5.2 for MR. Under $\|\sigma_k(\eta)\|\le\Gamma(1-\eta)\|F(x_k)\|$ near a limit point $x_*$ (5.6): $F(x_*)=0$, $x_k\to x_*$, and $\eta_k=\bar\eta_k$ eventually.
8. INB runs are MR runs with the curve (6.1), and $s_k=\sigma_k(\eta_k)$ throughout the loop (§6).

Further items, which are not milestones: Corollary 6.2 (exact Newton with backtracking takes full Newton steps eventually), Theorem 5.2 for Algorithm TL, Lemma 3.1 (existence of acceptable GIN steps), and Proposition 2.1 (the Goldstein–Armijo alpha condition implies (2.2) in the Euclidean norm).

Significance

The result. Theorem 6.1 is the global convergence guarantee for the inexact Newton backtracking method that Newton–Krylov solvers implement. It separates three outcomes: the iterates diverge, they accumulate only at points where $F'$ is singular, or they converge to a solution with invertible derivative and eventually take the unmodified inexact Newton steps. In the third case the local theory of Dembo, Eisenstat and Steihaug applies from some iteration on, so the forcing terms alone set the convergence rate. Theorem 5.2 is the template: §§6–8 of the paper derive the convergence of backtracking, equality-curve and dogleg-type methods from it by verifying (5.6).

Formalizing it. The results are proved on paper and none of them has a machine-checked proof that we know of. The mission produces a library of algorithm-run predicates with explicit while-loops for inexact Newton methods. It also checks the paper's reduction chain (INB is a run of MR, MR is a run of GIN) and the global convergence theorems in an arbitrary finite-dimensional norm.

Difficulty

The obvious argument fails at both ends. Sufficient decrease (2.2) alone gives monotonicity of $\|F(x_k)\|$, but not convergence of the iterates. On the page, $F(x)=x^2-1$ admits sequences satisfying (2.2) with both $\pm1$ as limit points. Convergence of $\|F(x_k)\|$ to $0$ needs $\sum(1-\eta_k)=\infty$, and nothing in the algorithm states this. In Theorem 5.2 it has to be derived from the exit level of the while-loop, which depends on a neighbourhood of $x_*$ where the linearization error is uniformly controlled. A solver must combine a limit-point argument (only infinitely many iterates are near $x_*$, not all of them) with the loop's worst-case backtracking factor $\theta_{\min}$. The invertibility of $F'(x_*)$ enters only through the bound (5.6) for the backtracking curve, which must be established uniformly in $k$.

Formalization scope

• Space and norm. $E$ is a finite-dimensional real normed space ([NormedAddCommGroup E] [NormedSpace ℝ E] [FiniteDimensional ℝ E]). This is exactly "$\mathbb R^n$ with an arbitrary norm". Fin n → ℝ (sup norm) and EuclideanSpace would each fix one norm. Only Proposition 2.1 assumes an inner product space, as the page does.
• Derivative. $F'$ is fderiv ℝ F, with ContDiff ℝ 1 F (the paper's standing assumption). Invertibility is ContinuousLinearMap.IsInvertible, and $F'(x)^{-1}$ is ContinuousLinearMap.inverse.
• Runs. An algorithm that does not break down is a predicate on infinite sequences (IsGINRun, IsMRRun, IsINBRun, plus IsTLRun, IsENBRun). The steps and levels the algorithm says to "find" or "choose" are data constrained only by the stated conditions. A while-loop is recorded by its number of passes $m_k$ and factors $\theta_{k,j}\in[\theta_{\min},\theta_{\max}]$. Every trial before the last fails the loop's test and the last one passes it. Termination is part of the run.
• Limit point is MapClusterPt xstar atTop x. "For all sufficiently large $k$" is ∀ᶠ k in atTop. The paper's "whenever $x_k$ is sufficiently near $x_*$ [and $k$ is sufficiently large]" is ∃ Γ, ∃ δ > 0, [∃ K,] ∀ k [≥ K], x k ∈ ball xstar δ → …. $\Gamma$ precedes $k$ ("independent of $k$"), and the "$k$ large" clause appears only in (3.2), where the page has it.
• Hypotheses as printed. Theorems 3.5 and 5.2 do not assume $F'(x_*)$ invertible. Theorem 5.2 does not assume $\sum(1-\eta_k)=\infty$. Lemma 5.1 is stated for one iteration of the loop shared by MR and TL, for every choice of factors.
• Trivializing readings ruled out. A run predicate without the "rejected" clause would allow needless backtracking, and one without the "accepted" clause would drop (2.2). Both clauses are present. A sorry-free check (F = id on $\mathbb R$) confirms that the INB, GIN and ENB run predicates and the goal's hypotheses are satisfiable, so the goal is not vacuous.
• Infrastructure. Continuity of operator inversion and uniform differentiability on neighbourhoods are in Mathlib. The run predicates and trial-level recursions are reusable for any line-search or backtracking analysis. Each milestone is stated independently; contributions in any order are welcome.

Selected references

• S. C. Eisenstat and H. F. Walker, Globally Convergent Inexact Newton Methods, SIAM J. Optim. 4(2) (1994) 393–422. https://doi.org/10.1137/0804022
• R. S. Dembo, S. C. Eisenstat and T. Steihaug, Inexact Newton Methods, SIAM J. Numer. Anal. 19(2) (1982) 400–408. https://doi.org/10.1137/0719025
• R. S. Dembo and T. Steihaug, Truncated-Newton algorithms for large-scale unconstrained optimization, Math. Program. 26 (1983) 190–212. https://doi.org/10.1007/BF02592055
• P. N. Brown and Y. Saad, Hybrid Krylov Methods for Nonlinear Systems of Equations, SIAM J. Sci. Stat. Comput. 11(3) (1990) 450–481. https://doi.org/10.1137/0911026
• S. C. Eisenstat and H. F. Walker, Choosing the Forcing Terms in an Inexact Newton Method, SIAM J. Sci. Comput. 17(1) (1996) 16–32. https://doi.org/10.1137/0917003
• J. E. Dennis and R. B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, SIAM Classics in Applied Mathematics 16 (1996). https://doi.org/10.1137/1.9781611971200