Operations Research

Motivation

Many problems in optimization and equilibrium modelling reduce to a system of equations $F(x) = 0$ whose map $F : \mathbb R^n \to \mathbb R^n$ is Lipschitz but not differentiable: reformulations of nonlinear complementarity problems through the componentwise minimum or the Fischer–Burmeister function, Karush–Kuhn–Tucker systems of constrained programs, and gradients of augmented Lagrangians all have kinks. Newton's method, $x^{k+1} = x^k - F'(x^k)^{-1}F(x^k)$, is the standard fast local solver for smooth systems, but it needs a derivative at every iterate.

Qi and Sun (Math. Programming 58, 1993) replaced the Jacobian by an arbitrary element of Clarke's generalized Jacobian and showed that the resulting method converges locally superlinearly under a regularity condition they called semismoothness, extending Mifflin's notion for functionals (Mifflin, SIAM J. Control Optim. 15, 1977) to vector-valued maps. This theorem is the foundation of the family of semismooth Newton methods used in complementarity, variational inequalities and PDE-constrained optimization.

Timeline. Robinson (1988) and Pang (Math. OR 15, 1990) studied Newton methods built on B-derivatives, with convergence proved under a strong Fréchet derivative at the solution; Kummer (1988) gave an abstract framework for Newton methods for nonsmooth equations; Qi and Sun (1993) proved local superlinear convergence for the generalized-Jacobian iteration under semismoothness and nonsingularity of $\partial F(x^*)$, with order $1+p$ under $p$-order semismoothness.

Setting

Let $F : \mathbb R^n \to \mathbb R^m$ be locally Lipschitz. By Rademacher's theorem $F$ is differentiable on a set $D_F$ of full measure; write $JF(y)$ for the Jacobian at $y \in D_F$. The generalized Jacobian 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 convex hull of all limits of Jacobians along sequences of differentiability points converging to $x$. The one-sided 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$, the limit of $Vh'$ over $V \in \partial F(x+th')$, $h' \to h$, $t \downarrow 0$ exists. For $0 < p \le 1$, $F$ is $p$-order semismooth at $x$ if in addition $Vh - F'(x;h) = O(\|h\|^{1+p})$ for $V \in \partial F(x+h)$, $h \to 0$.

For $m = n$, the nonsmooth Newton method is

$$x^{k+1} = x^k - V_k^{-1}F(x^k), \qquad V_k \in \partial F(x^k), \tag{3.2}$$

where any element of $\partial F(x^k)$ may be chosen at each step. A run 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$. A root $x^*$ ($F(x^*) = 0$) is regular when every $V \in \partial F(x^*)$ is nonsingular.

Formalization targets

Goal: Theorem 3.2, local superlinear convergence

Let $F$ be locally Lipschitz, $F(x^*) = 0$, $F$ semismooth at $x^*$, and every $V \in \partial F(x^*)$ nonsingular. Then there is $\delta > 0$ such that every $V \in \partial F(y)$ with $\|y - x^*\| < \delta$ is nonsingular, a Newton step from such a $y$ stays within $\delta$ of $x^*$, and every run with $\|x^0 - x^*\| < \delta$ satisfies

$$x^k \to x^*, \qquad \|x^{k+1} - x^*\| = o(\|x^k - x^*\|).$$

The goal asserts only the shape of the convergence (superlinear) and fixes no constants.

Stronger: Theorem 3.2, order $1 + p$

If moreover $F$ is $p$-order semismooth at $x^*$, $0 < p \le 1$, there are $\delta > 0$ and $C$ with

$$\|x^{k+1} - x^*\| \le C\|x^k - x^*\|^{1+p}$$

for every run started within $\delta$ of $x^*$.

Milestones

The milestones follow the paper's route: Proposition 2.1 (the limit in the definition of semismoothness is the directional derivative), Lemma 2.2 (Lipschitz continuity of $F'(x;\cdot)$ and its realisation by an element of $\partial F(x)$), Theorem 2.3 (semismoothness is equivalent to $Vh - F'(x;h) = o(\|h\|)$ and to the corresponding condition at differentiability points), the Remark's expansion (2.17), Proposition 3.1 (uniform invertibility near a regular point), the order-$(1+p)$ sentence of Theorem 3.2, and Corollary 2.5 (strong Fréchet differentiability implies semismoothness).

Significance

The theorem gives a locally superlinearly convergent method for Lipschitz equations with no smoothness beyond semismoothness at the root. Convex, smooth and subsmooth functions are semismooth, as are sums and scalar products of semismooth functions (the paper, citing Mifflin), and later work showed that the complementarity and KKT reformulations on which semismooth Newton solvers are built are semismooth as well; the order-$(1+p)$ variant gives local quadratic convergence for strongly semismooth maps. Mission II of this series treats the paper's global convergence theorem on a ball, and Mission III the semismoothness of augmented Lagrangian gradients, which supplies the application.

The results are proved in the paper. No machine-checked version of the generalized Jacobian, of semismoothness or of the nonsmooth Newton method is known to exist in Mathlib or on this platform; the platform's formalized Newton results concern one-dimensional $C^2$ functions (MetodosNumericos.newton_local_convergence) and smooth convex minimization. A complete development would provide the first formal library for Clarke's generalized Jacobian and semismooth maps.

Difficulty

The classical Newton proof compares $F(x^k)$ with its linearization $JF(x^*)(x^k - x^*)$ and uses continuity of the Jacobian at $x^*$. Here neither is available: $F$ need not be differentiable at $x^*$ or at any iterate, the element $V_k$ is chosen arbitrarily from a set, and $V_k$ need not be close to any fixed linear map. The comparison has to go through the directional derivative $F'(x^*; \cdot)$, which is only positively homogeneous, not linear. The analytic content therefore sits in Section 2: showing that semismoothness, defined through a limit over a set-valued map, controls $Vh - F'(x;h)$ uniformly in the direction, and that $F(x+h) - F(x) - F'(x;h)$ is small. Both rest on Clarke's mean-value inclusion and on compactness and upper semicontinuity of $\partial F$, none of which is in Mathlib. The superlinear rate also requires a uniform bound on $\|V^{-1}\|$ in a whole neighbourhood, not just at $x^*$.

Formalization scope

Everything lives in the namespace NonsmoothNewton.Local. Section 2 results are stated for maps between finite-dimensional real normed spaces $E \to G$ (the paper's $\mathbb R^n \to \mathbb R^m$ is the Euclidean instance); Section 3 results use EuclideanSpace ℝ (Fin n). Conventions fixed by the Lean statements:

• $JF$ is fderiv; the generalized Jacobian is the convex hull (no closure) of limits of fderiv along sequences $x_i \to x$ of differentiability points.
• $F'(x;h)$ is the one-sided limit over $t \downarrow 0$, never the two-sided lineDeriv; its value is a limUnder, used only where existence is a hypothesis or a consequence.
• Nonsingular means IsUnit in the ring of continuous linear endomorphisms; $\|V^{-1}\| \le C$ is a two-sided inverse of operator norm at most $C$.
• A run of (3.2) is encoded by the linear equation $V_k(x^{k+1} - x^k) = -F(x^k)$ with $V_k \in \partial F(x^k)$; all choices of $V_k$ are quantified, and $\delta$ is chosen before the run.
• Pinned asymptotics. The goal's rate is the proof's display (3.3), stated as IsLittleO along atTop; the printed Theorem 3.2 states only well-definedness and convergence. "Order $1+p$" is pinned as $\|x^{k+1}-x^*\| \le C\|x^k-x^*\|^{1+p}$ with $\delta$ and $C$ uniform over runs. Every $o(\|h\|)$ in (2.8), (2.9) and (2.17) is its $\varepsilon$–$\delta$ form with a non-strict inequality $\le \varepsilon\|h\|$, and every $O(\|h\|^{1+p})$ is an explicit constant and radius.
• The standing assumptions "$F$ locally Lipschitzian" of Sections 2 and 3 are hypotheses of every statement.
• The strong Fréchet derivative of Corollary 2.5 is Mathlib's HasStrictFDerivAt, which corrects the misprint $F(x)$ for $F(z)$ in the paper's display (2.16).

A trivializing formalization is ruled out: the update is not written with a junk inverse (which would make a singular step "well defined"), the generalized Jacobian is the paper's nonempty set rather than one that could be empty, and the theorem quantifies over every run rather than asserting that some run converges.

A complete development needs Clarke's mean-value inclusion (2.2), compactness and upper semicontinuity of $\partial F$ for locally Lipschitz maps (via Rademacher's theorem, available in Mathlib), and perturbation bounds for inverses of linear maps. The generalized-Jacobian and semismoothness layer is reusable beyond this mission, in particular for Missions II and III of this series. Contributions of proofs of any milestone, and of general lemmas about $\partial F$, 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, 1983 (SIAM reprint 1990). https://doi.org/10.1137/1.9781611971309
• R. Mifflin, Semismooth and semiconvex functions in constrained optimization, SIAM Journal on Control and Optimization 15 (1977) 959–972. https://doi.org/10.1137/0315061
• 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
• J. M. Ortega, W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, 1970 (SIAM reprint 2000). https://doi.org/10.1137/1.9780898719468

Motivation

Connecting a set of terminals by a network of direct links of least total length is one of the oldest problems of combinatorial optimization. R. C. Prim's 1957 paper in the Bell System Technical Journal (DOI) was motivated by the rate structure for Bell System leased-line services, in which the charge for connecting a set of terminals depends on the length of a shortest network connecting them. The paper states two local construction principles, P1 and P2, and shows that any sequence of their applications produces a shortest network, first for points in the plane and then for arbitrary connected labelled graphs with arbitrary real edge lengths. The paper's §V specialization of the principles, growing a single fragment, is what is now called Prim's algorithm, and its §IV statement is the form of the minimum spanning tree theorem used throughout network design, clustering and approximation algorithms.

Timeline. O. Borůvka (1926) solved the problem for an electrical network in Moravia; V. Jarník (1930) gave the single-fragment procedure; J. B. Kruskal (1956, Proc. AMS 7, 48–50) proved that adding globally shortest links avoiding cycles yields a shortest spanning tree; Prim (1957) gave the more permissive principles P1 and P2, which contain both the Jarník procedure and Kruskal's rule as special orders of application; E. W. Dijkstra (1959) rediscovered the single-fragment procedure.

Setting

Let $V$ be a finite set of $N$ terminals and $G$ a simple graph on $V$, the labelled graph whose edges are the possible links. Each edge $e$ carries a real length $w(e)$; lengths may be negative, zero, or tie. For a finite set $F$ of links, $H(F)$ denotes the graph on $V$ whose edges are the links of $F$.

• A spanning subtree of $G$ is a set $F$ of edges of $G$ such that $H(F)$ is a tree on $V$. Its length is $\ell_w(F) = \sum_{e \in F} w(e)$.
• A shortest spanning subtree (SSS) is a spanning subtree of least length among all spanning subtrees of $G$. Prim's dictionary is "shortest connection network (SCN) ↔ shortest spanning subtree (SSS)". $L(G,w)$ denotes that least length.
• Given the links $F$ made so far, the connected components of $H(F)$ are the isolated terminals (one terminal) and isolated fragments (two or more terminals).
• Principle 1: any isolated terminal $t$ can be connected to a nearest neighbor, a $G$-neighbor $n$ with $w(\{t,n\}) \le w(\{t,m\})$ for all $G$-neighbors $m$ of $t$.
• Principle 2: any isolated fragment $C$ can be connected to a nearest neighbor $n \notin C$ by a shortest available link $\{u,n\}$, $u \in C$; equivalently $\{u,n\}$ is a shortest edge of $G$ with one end in $C$ and the other outside.
• A construction is a sequence of links $e_0, e_1, \dots$, each an application of P1 or P2 with respect to the links before it. It is complete when it has $N-1$ links.

Only edges of $G$ are possible links; in Prim's distance table a missing edge has length $\infty$.

Formalization targets

Goal (§IV, p. 1396)

For every finite connected graph $G$ and every $w$,

$$\Bigl(\exists\ \text{a complete construction}\Bigr) \ \wedge\ \Bigl(\forall\ \text{complete constructions } e_0,\dots,e_{N-2}:\ \{e_0,\dots,e_{N-2}\} \text{ is a SSS of } G\Bigr).$$

This is the sentence "P1 and P2 will provide a SSS for any connected labelled graph with any set of real edge lengths." It fixes nothing about the order of applications, the component chosen, or the tie-breaking.

Milestones

1. Counting (§II, p. 1392): after any construction with $k$ links, $H(F)$ is acyclic with $N-k$ components; a complete construction is a spanning subtree; a construction with fewer than $N-1$ links can be extended.
2. Necessary Condition 1 (p. 1392): every terminal of a SSS is linked in it to at least one nearest neighbor.
3. Necessary Condition 2 (p. 1392): every fragment $S$ of a SSS, $\emptyset \ne S \ne V$, is linked in it to a nearest neighbor by a shortest available link.
4. Distinct lengths (§III, p. 1393): if the edge lengths are pairwise distinct, every link of every construction belongs to every SSS.
5. Continuity (§III, p. 1394): $w \mapsto L(G,w)$ is continuous.

Significance

The goal is the correctness theorem of a whole family of greedy minimum spanning tree procedures at once: Jarník–Prim (one growing fragment), Kruskal (globally shortest link first) and Borůvka-style interleavings all produce sequences of P1/P2 applications. Because lengths are arbitrary reals, it also covers maximum spanning trees by a sign change (p. 1397) and graphs that are not complete.

The result is classical and fully proved in the literature. What this mission adds is a machine-checked statement in exactly Prim's generality. Mathlib has spanning trees of connected graphs (SimpleGraph.Connected.exists_isTree_le) and the edge count of trees, but no minimum spanning tree theory. Existing Prove2Me items on minimum spanning trees are either restricted to complete graphs with distance matrices or state a cut property in existence form at a single vertex; none states Prim's principles or his necessary conditions.

Difficulty

The obvious argument, "each link P1 or P2 adds belongs to the shortest network", uses a unique shortest network, and that fails with ties: when two links tie, a P1/P2 link need not lie in a given SSS. Prim's own treatment of ties (§III) is an informal perturbation argument; the formal statement must hold for every tie-breaking choice made during a construction, not only for a generic perturbed instance. Negative lengths remove the easy reading "shortest connected spanning subgraph": the minimum must range over trees only. The statements also involve the component structure of $H(F)$ as it changes during a construction, and tree paths in an arbitrary, not necessarily complete, graph.

Formalization scope

Namespace ShortestConnection.Principles, Mathlib SimpleGraph. Conventions:

• $V$ is a Fintype with decidable equality; $G$ is a SimpleGraph V (at most one link per pair, no loops, which is Prim's setting). Lengths are w : Sym2 V → ℝ; only values on edges of $G$ matter.
• Link sets are Finset (Sym2 V); linkGraph F is SimpleGraph.fromEdgeSet F. A spanning subtree requires ↑F ⊆ G.edgeSet and (linkGraph F).IsTree.
• An isolated fragment is a whole connected component of linkGraph F; the P2 condition is a single inequality against every $G$-edge leaving it, which is equivalent to "nearest neighbor and shortest link" in Prim's sense.
• A construction is a List (Sym2 V) checked entrywise against l.take i; complete means length Fintype.card V - 1 (natural subtraction, used only for nonempty $V$).
• $L$ is sInf of the lengths of spanning subtrees; continuity is in the product topology.

Implicit hypotheses made explicit: $G$ connected (hence $V \ne \emptyset$) wherever an SSS or a complete construction is involved; at least two terminals for Necessary Condition 1; $S$ nonempty and $S \ne V$ for Necessary Condition 2; pairwise distinct edge lengths only in milestone 4, as in the paper's temporary assumption.

The goal's existence clause rules out a vacuous formalization in which no complete construction exists; the step predicates are defined from lengths and components only, never through shortest spanning subtrees, and they are not restricted to one growing fragment or to the globally shortest link.

Needed infrastructure: tree exchange (adding an edge to a spanning tree creates one cycle; removing any other cycle edge yields a spanning tree), component counts under edge addition, and minima of finitely many continuous functions. The exchange and counting lemmas are reusable for any matroid-greedy or spanning-tree mission. Contributions of intermediate lemmas, and proofs of the milestones in any order, are welcome.

Selected references

• R. C. Prim, Shortest Connection Networks And Some Generalizations, Bell System Technical Journal 36 (1957), 1389–1401. https://doi.org/10.1002/j.1538-7305.1957.tb01515.x
• J. B. Kruskal, On the shortest spanning subtree of a graph and the traveling salesman problem, Proceedings of the AMS 7 (1956), 48–50. https://doi.org/10.1090/S0002-9939-1956-0078686-7
• V. Jarník, O jistém problému minimálním, Práce Moravské Přírodovědecké Společnosti 6 (1930), 57–63.
• O. Borůvka, O jistém problému minimálním, Práce Moravské Přírodovědecké Společnosti 3 (1926), 37–58.
• R. L. Graham, P. Hell, On the history of the minimum spanning tree problem, Annals of the History of Computing 7 (1985), 43–57. https://doi.org/10.1109/MAHC.1985.10011

Motivation

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

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

Timeline:

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

Setting

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

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

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

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

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

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

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

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

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

Formalization targets

Goal: quadratic termination of SQN

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

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

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

Milestones

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

Significance

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

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

Difficulty

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

Formalization scope

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

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

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

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

Selected references

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

Motivation

Newton's method is the standard second-order algorithm for unconstrained minimization, but without safeguards it has no global guarantee: far from a minimizer the Newton step can increase the objective, and at a point where the Hessian is indefinite the step can head for a saddle point or a maximum. The usual repairs (line search, trust regions, Levenberg–Marquardt damping) come with convergence proofs, but for nonconvex objectives those proofs typically give no rate at all, or only the rate of the gradient method.

Nesterov and Polyak (Math. Program. 108 (2006) 177–205) proposed to regularize the second-order Taylor model of the objective with a cubic term and to take as the next iterate a global minimizer of the regularized model. They showed that the resulting method has a global worst-case rate of convergence to points satisfying the second-order necessary conditions, for every objective with a Lipschitz continuous Hessian and without any convexity. That rate, $O(k^{-2/3})$ for the gradient norm, is better than the $O(k^{-1/2})$ of the gradient method. It became the reference point for the complexity theory of nonconvex second-order optimization: adaptive variants (Cartis, Gould and Toint, Math. Program. 127 (2011) 245–295) and lower bounds showing that $O(\epsilon^{-3/2})$ iterations are optimal among second-order methods (Carmon, Duchi, Hinder and Sidford, Math. Program. 184 (2020) 71–120) are stated against it.

This mission formalizes the general convergence result of that paper, Theorem 1 of Section 3, together with the properties of the cubic step from Section 2 on which it rests.

Setting

Let $F \subseteq \mathbb{R}^n$ be a closed convex set with nonempty interior, and let $f$ be twice differentiable on $F$ with gradient $f'(x)$ and Hessian $f''(x)$. A starting point $x_0 \in \operatorname{int} F$ is fixed, and $F$ is assumed to contain the level set $\mathcal{L}(f(x_0)) = \{x \in \mathbb{R}^n : f(x) \le f(x_0)\}$ in its interior. Assumption 1: the Hessian is Lipschitz continuous on $F$ in the spectral norm, $\|f''(x) - f''(y)\| \le L\|x - y\|$ for all $x, y \in F$, with $L > 0$.

For a parameter $M > 0$ the cubic model of $f$ at $x$ is

$$m_{M,x}(y) = \langle f'(x), y - x\rangle + \tfrac12 \langle f''(x)(y - x), y - x\rangle + \tfrac{M}{6}\|y - x\|^3 .$$

The cubic-regularized Newton step $T_M(x)$ is any global minimizer of $m_{M,x}$ over $\mathbb{R}^n$; it exists because the model is continuous and coercive. Write $r_M(x) = \|x - T_M(x)\|$ and $\bar f_M(x) = f(x) + \min_y m_{M,x}(y)$.

The cubic regularization of Newton method (3.3) fixes $L_0 \in (0, L]$, starts at $x_0$ and, for $k \ge 0$, chooses $M_k \in [L_0, 2L]$ such that $f(T_{M_k}(x_k)) \le \bar f_{M_k}(x_k)$, then sets $x_{k+1} = T_{M_k}(x_k)$. The choice $M_k = L$ always passes the test.

Write $\lambda_n(A)$ for the smallest eigenvalue of a symmetric matrix $A$. The measure of local optimality is

$$\mu_M(x) = \max\Big\{ \sqrt{\tfrac{2}{L + M}\,\|f'(x)\|},\ -\tfrac{2}{2L + M}\,\lambda_n(f''(x)) \Big\}.$$

It is nonnegative and vanishes exactly when $f'(x) = 0$ and $f''(x) \succeq 0$.

Formalization targets

Goal: Theorem 1, inequality (3.4)

If $f(x) \ge f^*$ for all $x \in F$, then every run of method (3.3) satisfies, for every $k \ge 1$,

$$\min_{1 \le i \le k} \mu_L(x_i) \le \frac{8}{3}\cdot\left(\frac{3\,(f(x_0) - f^*)}{2k\cdot L_0}\right)^{1/3}.$$

The constant $8/3$ and the exponent $1/3$ are the paper's; the statement holds for every admissible choice of the parameters $M_k$ and of the global minimizers $x_{k+1}$.

Milestones, in attack order

1. Lemma 1 (2.2): $\|f'(y) - f'(x) - f''(x)(y - x)\| \le \tfrac12 L\|y - x\|^2$ on $F$.
2. Eq. (2.5): $f'(x) + f''(x)(T - x) + \tfrac12 M\|T - x\|(T - x) = 0$ for $T = T_M(x)$.
3. Proposition 1 (2.7): $f''(x) + \tfrac12 M r_M(x) I \succeq 0$.
4. Lemma 2 (2.8): $\langle f'(x), x - T_M(x)\rangle \ge 0$ when $f(x) \le f(x_0)$.
5. Lemma 4 (2.11): $f(x) - \bar f_M(x) \ge \tfrac{M}{12} r_M(x)^3$.
6. Lemma 4 (2.12): for $M \ge L$, $T_M(x) \in F$ and $f(T_M(x)) \le \bar f_M(x)$.
7. Lemma 3 (2.9): $\|f'(T_M(x))\| \le \tfrac12(L + M) r_M(x)^2$ when $T_M(x) \in F$.
8. Lemma 5: $\mu_M(T_M(x)) \le r_M(x)$.
9. Theorem 1, first claim: $\sum_{i \ge 0} r_{M_i}(x_i)^3 \le \tfrac{12}{L_0}(f(x_0) - f^*)$.
10. Theorem 1, second claim: $\lim_{i\to\infty} \mu_L(x_i) = 0$.

Significance

Inequality (3.4) is a global, dimension-free complexity bound for reaching approximate second-order stationarity. It controls both the gradient norm, $\min_{1\le i\le k}\|f'(x_i)\| = O(k^{-2/3})$, and the most negative curvature, $\max\{0, -\lambda_n(f''(x_i))\} = O(k^{-1/3})$, along the best iterate, from a single scalar potential $f(x_0) - f^*$. The second claim of Theorem 1 gives the asymptotic counterpart: every limit point satisfies the second-order necessary conditions. Section 4 of the paper derives its faster rates for star-convex and gradient-dominated functions from the same Section 2 lemmas.

The result is proved on paper and widely cited; to our knowledge no machine-checked proof of it or of the Section 2 lemmas exists. A formalization adds a checked statement of the method with its exact constants, and reusable facts about global minimizers of cubic models (Proposition 1 in particular) that the companion missions on star-convex, gradient-dominated and locally quadratic convergence also rely on.

Difficulty

Most steps are short inequalities, but two are not. Proposition 1 is a statement about a global minimizer of a nonconvex function: the first- and second-order conditions of a local minimizer give only $f''(x) + \tfrac12 M r I + \tfrac{M}{2r}(T - x)(T - x)^\top \succeq 0$, which is weaker. The natural first attempt, "take the second-order optimality condition of the model at $T$", therefore fails. The paper proves it in Section 5.1 through a one-dimensional dual characterization of the minimizer.

The second is Lemma 2's second claim, used for (2.12): showing that $T_M(x)$ stays in $F$ requires a boundary argument along the segment from $x$ to $T_M(x)$, since the Taylor bounds are only available inside $F$. The remaining work is calculus in $\mathbb{R}^n$: the integral form of Taylor's theorem for the gradient under a Lipschitz Hessian, and eigenvalue perturbation for the second entry of $\mu$.

Formalization scope

The space is EuclideanSpace ℝ (Fin n) for arbitrary n : ℕ. The gradient and Hessian are maps g and H with HasGradientAt f (g x) x and HasFDerivAt g (H x) x at every x ∈ F. At boundary points of $F$ this asks for two-sided derivatives, a mild strengthening of "twice differentiable on $F$". The Lipschitz condition uses the operator norm, which is the spectral norm. $T_M(x)$ is represented by the predicate IsCubicStep (global minimizer of cubicModel), and every lemma is stated for every such minimizer. The run predicate IsCubicNewtonRun is 0-based. It writes $\bar f_{M_k}(x_k)$ as $f(x_k)$ plus the model value at $x_{k+1}$, which is the minimum because $x_{k+1}$ attains it. $\lambda_n$ is lamMin, the Rayleigh-quotient infimum over the unit sphere, which equals the smallest eigenvalue for the (symmetric) Hessian. The lower bound $f^*$ is required on $F$ only. The minimum over $1 \le i \le k$ is written as the existence of an index attaining the bound.

A stationary point of the cubic model is not an admissible step, and the run must keep the test $M_k \in [L_0, 2L]$ and the acceptance test. Replacing the step by any point with $f(x_{k+1}) \le f(x_k)$ makes the goal false, and dropping the square root in $\mu_M$ makes Lemma 5 false. The statements rule out all three. Lemma 5 carries the hypothesis $T_M(x) \in F$, which its printed proof uses and which holds at every iterate.

A complete development needs the Taylor bounds (2.2)–(2.3) for vector-valued derivatives on convex sets, and first- and second-order optimality for the cubic model. It also needs a proof of Proposition 1 (Section 5.1 or any other correct argument) and eigenvalue perturbation via Rayleigh quotients. The cubic-model lemmas and Proposition 1 are reusable across the whole series. Proofs of any milestone, alternative proofs of Proposition 1, and general Mathlib-level lemmas about Rayleigh quotients are welcome.

Selected references

• Yu. Nesterov and B. T. Polyak, Cubic regularization of Newton method and its global performance, Mathematical Programming, Ser. A 108 (2006) 177–205. https://doi.org/10.1007/s10107-006-0706-8
• C. Cartis, N. I. M. Gould and Ph. L. Toint, Adaptive cubic regularisation methods for unconstrained optimization. Part I: motivation, convergence and numerical results, Mathematical Programming 127 (2011) 245–295. https://doi.org/10.1007/s10107-009-0286-5
• Y. Carmon, J. C. Duchi, O. Hinder and A. Sidford, Lower bounds for finding stationary points I, Mathematical Programming 184 (2020) 71–120. https://doi.org/10.1007/s10107-019-01406-y
• Yu. Nesterov, Introductory Lectures on Convex Optimization: A Basic Course, Kluwer, 2004. https://doi.org/10.1007/978-1-4419-8853-9

Motivation

The Critical-Path Method (CPM) was introduced by J. E. Kelley, Jr. (Remington Rand) and M. R. Walker (du Pont) in Critical-Path Planning and Scheduling (Proc. Eastern Joint Computer Conference, 1959, pp. 160–173, doi:10.1145/1460299.1460318). Together with PERT, developed at the same time for the Polaris programme, it became the standard way to plan and schedule large projects in construction, maintenance and engineering, and it is taught in every introductory operations research course.

The paper reduces project scheduling to arithmetic on a directed acyclic graph: the earliest and latest times of the project's events are computed by two recursions, and the jobs whose timing has no slack, the critical jobs, are singled out by an equation. Its central structural claim is that critical jobs, when they exist, form a path from the start of the project to its end. The paper states this without proof ("a detailed development being reserved for a separate paper", p. 161). This mission formalizes that claim and the facts about the two recursions on which it rests.

Setting

A project network has $n+1$ events labelled $0,1,\dots,n$ with $n \ge 1$: event $0$ is the origin and event $n$ the terminus. A job is an arrow from an event $i$ to an event $j$, written job $(i,j)$; the jobs form a finite set $P$ of ordered pairs of events. Two standing assumptions of the paper (pp. 161–162) are part of the model:

1. every job has $i < j$ (events are labelled so that the head of an arrow has the larger label);
2. origin precedes and terminus follows every event: for every event $k$ there are chains of jobs from $0$ to $k$ and from $k$ to $n$.

Each job has a real duration $y_{ij}$. The earliest event times $t^{(0)}$ are given by display (1) of the paper,

$$t_0^{(0)} = 0,\qquad t_j^{(0)} = \max\,[\,y_{ij} + t_i^{(0)} \mid i<j,\ (i,j)\in P\,],\quad 1\le j\le n,$$

and, for a project completion time $\lambda \ge t_n^{(0)}$, the latest event times $t^{(1)}$ by display (2),

$$t_n^{(1)} = \lambda,\qquad t_i^{(1)} = \min\,[\,t_j^{(1)} - y_{ij} \mid i<j,\ (i,j)\in P\,],\quad 0\le i\le n-1.$$

The maximum time available for job $(i,j)$ is $t_j^{(1)} - t_i^{(0)}$. The job is critical if this equals its duration, $t_j^{(1)} - t_i^{(0)} = y_{ij}$, and a floater if it exceeds it. A critical path is a contiguous path of critical jobs from origin to terminus: events $0 = v_0, v_1, \dots, v_k = n$ with every $(v_{r-1}, v_r)$ a critical job of $P$.

In the Lean development these are ProjectNetwork n (with field P), earliest N y, latest N y λ, maxTimeAvailable, IsCritical, IsFloater and IsCriticalPath, in the namespace CriticalPath.Events.

Formalization targets

Goal: critical jobs force $\lambda = t_n^{(0)}$ and a critical path (p. 163)

For every project network, durations $y$ and completion time $\lambda \ge t_n^{(0)}$,

$$\bigl(\exists (i,j)\in P,\ t_j^{(1)} - t_i^{(0)} = y_{ij}\bigr) \;\Longrightarrow\; \lambda = t_n^{(0)} \ \wedge\ \exists\ \text{a critical path}.$$

This is the paper's "A project will contain critical jobs only when $\lambda = t_n^{(0)}$. If a project does contain critical jobs, then it also contains at least one contiguous path of critical jobs through the project diagram from origin to terminus." Only the "only when" direction is asserted, as on the page.

Milestones

1. Display (1), pp. 162–163. $t^{(0)}$ is the least vector $t$ with $t_0 = 0$ and $y_{ij} \le t_j - t_i$ for every job.
2. Display (2), p. 163. For $\lambda \ge t_n^{(0)}$, $t_n^{(1)} = \lambda$ and $t^{(1)}$ is the greatest vector $t$ with $t_n \le \lambda$ and $y_{ij} \le t_j - t_i$ for every job.
3. Critical or floater, p. 163. For $\lambda \ge t_n^{(0)}$, $t_i^{(0)} \le t_i^{(1)}$ for every event, and every job is critical or a floater: $t_j^{(1)} - t_i^{(0)} \ge y_{ij}$.
4. Delay of a critical job, p. 163. Lengthening a critical job by $\delta \ge 0$ raises $t_n^{(0)}$ by exactly $\delta$.

Significance

The result. The theorem is what makes the method's name meaningful: it says that the jobs without slack are not scattered but line up along an origin–terminus path, and that such jobs exist only when the project is scheduled at its earliest possible completion time. Project managers use this to decide which jobs to watch, which to expedite, and which may slip; the delay statement (milestone 4) is the quantitative form of that advice. The characterisations of (1) and (2) as least and greatest feasible schedules are the bridge between CPM and linear programming: they identify $t^{(0)}$ and $t^{(1)}$ with extreme solutions of the system of difference constraints $y_{ij} \le t_j - t_i$, which the paper's own §3 uses to build the project cost curve.

Formalizing it. The results are classical and folklore, but the paper proves none of them, and textbook treatments usually define the critical path as a longest path, which makes the goal a tautology. This mission states the claims with the paper's own definitions: criticality by the float equation, event times by the recursions. To the best of current knowledge no machine-checked version of these statements for activity-on-arrow networks exists; the platform has a related activity-on-node development (Brucker and Knust, Complex Scheduling) in which the critical path is defined as a longest path.

Difficulty

The recursions (1) and (2) are local: each event looks only at its immediate predecessors or successors. The goal is global: from one critical job it asserts a statement about the whole completion time and a whole origin–terminus path. The float equation $t_j^{(1)} - t_i^{(0)} = y_{ij}$ mixes a quantity computed forward from the origin with one computed backward from the terminus, and neither recursion alone says anything about the other. The naive reading "a critical job lies on a longest path" is not available as a definition: it is, in substance, what has to be established from the recursions. The formal overhead is the well-founded recursion on the labels, in both directions, and the bookkeeping of lists of events forming a path.

Formalization scope

Events are Fin (n + 1), origin 0, terminus Fin.last n, with 1 ≤ n. Jobs are a Finset of ordered pairs, so there is at most one job per ordered pair. The standing assumptions (labels increase along jobs; origin precedes and terminus follows every event, via Relation.ReflTransGen) are fields of the structure ProjectNetwork and are never dropped. Durations and times are real numbers; durations are a function Fin (n+1) → Fin (n+1) → ℝ read only on jobs of P, with no sign condition, as in the paper's deterministic case.

The event times are defined by the recursions (1) and (2) themselves, by well-founded recursion on the label with Finset.sup'/Finset.inf' over the predecessor/successor set; these sets are nonempty by the standing assumptions, so no fallback value exists. The latest times are defined for every real $\lambda$; the paper's assumption $\lambda \ge t_n^{(0)}$ is a hypothesis of every theorem that uses them.

Disclosed readings: "earliest time occurance" (milestone 1) and "latest time … relative to a fixed project completion time" (milestone 2) are read as least and greatest vectors satisfying the job constraints $y_{ij} \le t_j - t_i$ (the paper's constraint (8), p. 165); milestone 3 is the fact implicit in the dichotomy "critical or floater"; "comparable delay" (milestone 4) is read as an exact delay of $\delta$ in $t_n^{(0)}$ for $\delta \ge 0$.

A trivializing formalization is ruled out: defining a critical job or path through longest paths, or taking $t^{(0)}$ and $t^{(1)}$ as arbitrary functions satisfying (1) and (2), would make the goal a restatement of its definitions; here criticality is the float equation and the times are computed by the recursions. Dropping the reachability assumptions would make (2) ill-defined at events without successors.

Contributions welcome: proofs of the milestones, general lemmas on longest paths in finite labelled DAGs and on difference constraints $y_{ij} \le t_j - t_i$, which are reusable for the companion mission on the project cost curve.

Selected references

• J. E. Kelley, Jr. and M. R. Walker, Critical-Path Planning and Scheduling, Papers presented at the December 1–3, 1959, Eastern Joint IRE-AIEE-ACM Computer Conference, pp. 160–173, 1959. doi:10.1145/1460299.1460318
• J. E. Kelley, Jr., Critical-Path Planning and Scheduling: Mathematical Basis, Operations Research 9(3), pp. 296–320, 1961. doi:10.1287/opre.9.3.296
• P. Brucker and S. Knust, Complex Scheduling, 2nd ed., Springer, 2012. doi:10.1007/978-3-642-23929-8

Motivation

Many operational systems move between a finite number of states at random times: a machine alternates between working and repair, a queue between occupancy levels, an inventory between stock positions. When the time spent in a state is not exponential and not a fixed period, neither discrete-time Markov decision processes nor continuous-time Markov chains describe the system faithfully. William S. Jewell's 1963 paper Markov-Renewal Programming. I extends Howard's Markov decision processes to Markov-renewal processes (also called semi-Markov processes), in which the time between transitions is a random variable whose law depends on the current state, the next state, and the decision taken. The resulting model, now called a semi-Markov decision process, is standard in maintenance, queueing control and reliability.

Timeline:

• 1954: Lévy, Smith and Takács independently introduce Markov-renewal and semi-Markov processes; Pyke later surveys them.
• 1960: Howard, Dynamic Programming and Markov Processes, introduces policy iteration for finite discrete-time Markov decision processes.
• 1962: Blackwell, Discrete Dynamic Programming, shows that for the discounted discrete-time problem a stationary policy is optimal among all policies.
• 1963: Jewell formulates Markov-renewal programming, with a continuous discount factor α, and carries Howard's algorithm and Blackwell's stationarity result over to it. Part II of the paper treats the undiscounted (infinite-return) models.

Setting

A Markov-renewal program has a finite set of states $S$ (the paper's $i = 1, \dots, N$) and a finite, nonempty set of alternatives (the paper's $z = 1, \dots, Z$), each available in every state. For each alternative $z$ and states $i, j$ it specifies:

• a transition probability $p^z_{ij} \ge 0$, with $\sum_j p^z_{ij} = 1$;
• a sojourn-time distribution $F^z_{ij}$, the law of the time $\tau$ between entering $i$ and moving to $j$, with $\tau \ge 0$ and $F^z_{ij}(0) = 0$;
• for each continuous discount factor $\alpha > 0$, a real number $\rho^z_{ij}(\alpha)$, the expected discounted return earned during that transition.

The Laplace–Stieltjes transform $\tilde f^z_{ij}(s) = \int_0^\infty e^{-st}\,dF^z_{ij}(t)$ is the expected discount $\mathbb E[e^{-s\tau}]$ over one interval. The average one-step return is $\rho^z_i(\alpha) = \sum_j p^z_{ij}\rho^z_{ij}(\alpha)$, and for a vector of returns $v$ the test quantity is

$$\rho^z_i(\alpha) + \sum_j p^z_{ij}\,\tilde f^z_{ij}(\alpha)\,v_j .$$

A stationary policy is a map $d : S \to A$; a nonstationary policy is a sequence $\pi = (\pi_0, \pi_1, \dots)$ of such maps, $\pi_k$ being used at the $k$-th transition. The $n$-step return $V^\pi_i(n)$ of a policy, with boundary rewards $V_i(0,\alpha)$, is the test quantity of $\pi_0(i)$ applied to the $(n-1)$-step return of the shifted policy; the optimal $n$-step return $V_i(n,\alpha)$ of equation (6) replaces $\pi_0(i)$ by a maximum over $z$. The value-determination equations (15) of a stationary policy $d$ are $v_i = \rho^{d(i)}_i(\alpha) + \sum_j p^{d(i)}_{ij}\tilde f^{d(i)}_{ij}(\alpha) v_j$.

The algorithm of Fig. 1 alternates two steps: solve (15) for the current policy, then in every state pick an alternative maximizing the test quantity, keeping the old alternative if it still attains the maximum. It stops when two successive policies are identical.

Formalization targets

Goal (p. 947)

For every $\alpha > 0$: (15) has a unique solution for every stationary policy, and every run $(d_k, v_k)$ of Fig. 1 reaches $d_{K+1} = d_K$ with $K < Z^N$, where

$$\lim_{n\to\infty} V^{d_K}_i(n) = (v_K)_i,\qquad \lim_{n\to\infty} V^\pi_i(n) \le (v_K)_i \ \ \forall \pi,\qquad \lim_{n\to\infty} V_i(n,\alpha) = (v_K)_i ,$$

for every state $i$ and all boundary rewards, every limit existing. This is the paper's "the algorithm of Fig. 1 produces an optimal, stationary policy that is as good as any optimal, nonstationary policy".

Milestones

1. p. 945: $0 \le p^z_{ij}\tilde f^z_{ij}(s) < 1$ for $s > 0$.
2. Claim (a): (15) has exactly one solution for each stationary policy.
3. Eq. (15): the $n$-step return of a stationary policy converges to a solution of (15).
4. p. 946: $I - \tilde q(\alpha)$ is invertible and $([I-\tilde q(\alpha)]^{-1})_{ii} \ge 1$.
5. Claim (b): a change of policy raises the return of some state and lowers none.
6. Claim (c): a policy reproduced by the improvement step is optimal among stationary policies.
7. Claim (d): a run of Fig. 1 terminates within $Z^N$ cycles.
8. Eq. (14): $V_i(n,\alpha)$ converges, independently of the boundary rewards, to a solution of $v_i = \max_z\{\rho^z_i(\alpha) + \sum_j p^z_{ij}\tilde f^z_{ij}(\alpha) v_j\}$.
9. p. 946: some stationary policy's return dominates the limiting return of every nonstationary policy.

Significance

The result says that the infinite-step discounted Markov-renewal program is solved exactly, in finitely many cycles, by a finite-dimensional algorithm, and that the answer is a stationary policy. As the paper notes, this matters operationally because a nonstationary policy is hard to follow. The sojourn distributions enter only through the numbers $\tilde f^z_{ij}(\alpha)$, so the same algorithm serves any sojourn-time law. For fixed $\alpha$ the model is a discounted Markov decision process whose discount factor depends on the transition, which contains Howard's and Blackwell's constant-discount problem as the case of intervals of fixed length (p. 943).

The results are classical and proved in the literature. The paper itself refers the proofs to Howard and Blackwell. Machine-checked versions exist on this platform for finite stochastic shortest path and constant-discount problems (Bertsekas, Dynamic Programming and Optimal Control, Prop. 7.2.2 and 7.3.1, mission Dynamic Programming and Optimal Control VII). No formal treatment of Markov-renewal programs, of transition-dependent discounting, or of the retention rule of Fig. 1 is known to exist. The mission produces a verified policy-iteration theorem for semi-Markov decision processes, with an explicit termination bound and the comparison against nonstationary policies.

Difficulty

The discount over one transition, $\tilde f^z_{ij}(\alpha)$, varies with $i$, $j$ and $z$, so the problem is not a constant-$\gamma$ contraction of textbook form; the relevant bound is that every row of $\tilde q(\alpha)$ sums to less than one, which rests on $F^z_{ij}(0) = 0$. Entries in $[0, 1)$ alone, the paper's stated justification of Claim (a), do not make $I - \tilde q(\alpha)$ invertible: the $2 \times 2$ matrix with every entry $1/2$ has entries in $[0,1)$, yet $I$ minus it is singular.

Finite termination is not automatic either. If the improvement step may switch between tied maximizers, the iterates can cycle forever between two policies with equal returns; the retention rule of Fig. 1 excludes this, and Claim (b) must deliver a strict increase in some state, with no decrease anywhere, to rule out revisiting a policy. Comparing with nonstationary policies requires controlling returns of arbitrary policy sequences, whose limits must be shown to exist, not assumed.

Formalization scope

• States and alternatives are finite types; alternatives are nonempty. Every alternative is available in every state.
• $F^z_{ij}$ is a probability measure on $\mathbb R$ with no mass on $(-\infty, 0]$. The transform is integrated over $(0, \infty)$, which carries all the mass.
• The one-transition returns $\rho^z_{ij}(\alpha)$ are arbitrary real numbers, a generalization of the paper's Stieltjes integral (4), which is not formalized. The reward functions $R^z_{ij}(t\mid\tau)$ do not appear.
• Returns of policies are defined by the one-step recursion (the policy form of (6)); the Markov-renewal process is not built as a stochastic process.
• Policies are the paper's: deterministic and Markov, nonstationary ones indexed by the number of transitions made. Randomized and history-dependent policies are not in the comparison class.
• The following informal words are read as follows. "Solve the set of simultaneous equations": (15) has exactly one solution. "Strictly increases the expected return of at least one state": no state's return decreases and one strictly increases, under the hypothesis that the policy changed. "No other policy can lead to higher expected returns" in Claim (c): no stationary policy. "Terminates in a finite number of cycles": two successive policies coincide at some cycle $K < Z^N$. "If there is no improvement in the test quantity, retain the same alternative": the old alternative is kept whenever it attains the maximum. "Optimal" and "as good as any nonstationary policy": the limiting return of the returned policy dominates that of every policy from every state, for all boundary rewards. "$\lim_{n\to\infty} V_i(n,\alpha)$": the limit is proved to exist. "$\max_z$": a maximum over the finite nonempty set of alternatives.
• Eq. (14) is printed with $v_i(\alpha)$ inside the sum over $j$; the formalization uses $v_j(\alpha)$, as (6), (15) and Fig. 1 do.
• Every statement fixes one $\alpha > 0$. The undiscounted models (16)–(19), the finite-time and mixed-horizon models (9)–(13), and the infinite-time case of the stationarity result are out of scope.
• The return of a policy is never defined through a matrix inverse, whose Mathlib value for a singular matrix is $0$; (15) is a predicate, and the goal asserts unique solvability, so a vacuous reading through junk inverses or assumed limits is excluded.

Reusable infrastructure: bounds for substochastic matrices with row sums below one (invertibility, Neumann series, nonnegative inverse), convergence of iterated Bellman operators with transition-dependent discount, and the policy-iteration termination argument with a tie-breaking rule. Proofs of any milestone and alternative arguments are welcome.

Selected references

• W. S. Jewell, Markov-Renewal Programming. I: Formulation, Finite Return Models, Operations Research 11(6), 938–948, 1963. https://doi.org/10.1287/opre.11.6.938
• W. S. Jewell, Markov-Renewal Programming. II: Infinite Return Models, Example, Operations Research 11(6), 949–971, 1963. https://doi.org/10.1287/opre.11.6.949
• R. A. Howard, Dynamic Programming and Markov Processes, MIT Press, 1960.
• D. Blackwell, Discrete Dynamic Programming, Annals of Mathematical Statistics 33(2), 719–726, 1962. https://doi.org/10.1214/aoms/1177704593
• R. Pyke, Markov Renewal Processes: Definitions and Preliminary Properties, Annals of Mathematical Statistics 32(4), 1231–1242, 1961. https://doi.org/10.1214/aoms/1177704863
• D. P. Bertsekas, Dynamic Programming and Optimal Control, Vol. I, 3rd ed., Athena Scientific, 2005, Section 7.2–7.3.

Motivation

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

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

Timeline:

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

Setting

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

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

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

Formalization targets

Goal: Theorem 6:3

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

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

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

Milestones

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

Significance

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

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

Difficulty

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

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

Formalization scope

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

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

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

Selected references

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