Optimization

Why reduce scenarios

Multistage and two-stage stochastic programs are solved numerically on a discrete probability distribution: a finite set of scenarios $\omega_1,\dots,\omega_N$ with probabilities $p_1,\dots,p_N$. The size of the resulting optimization problem grows with $N$, and scenario sets produced by sampling or by historical data are often far too large to be solved directly. Scenario reduction replaces the original distribution by one supported on a small subset of the scenarios, chosen so that the optimal value and solutions of the stochastic program change as little as possible.

Stability theory for stochastic programs (Rachev and Römisch, 2002) shows that this change is controlled by a probability distance of Fortet–Mourier type, which for discrete measures is bounded by the value of a transportation problem. Dupačová, Gröwe-Kuska and Römisch (2003) turned this into a combinatorial problem and proposed greedy backward and forward heuristics. Heitsch and Römisch (2003) gave faster versions of both heuristics; the forward version, fast forward selection, is the subject of this mission. Implementations of these reduction heuristics are distributed with the GAMS modelling system (SCENRED) and are used in energy and finance applications of stochastic programming.

Setting

Let $E$ be a finite-dimensional real vector space with a norm $\|\cdot\|$, let $\omega_0\in E$, and let $h:[0,\infty)\to[0,\infty)$ be continuous and nondecreasing with $h(0)=0$. The cost between two points of $E$ is

$$c(\omega,\tilde\omega)=\max\bigl\{1,\,h(\|\omega-\omega_0\|),\,h(\|\tilde\omega-\omega_0\|)\bigr\}\,\|\omega-\tilde\omega\| .$$

It is nonnegative, symmetric, and zero on the diagonal.

The original distribution is $P=\sum_{i=1}^N p_i\delta_{\omega_i}$ with $p_i>0$ and $\sum_i p_i=1$. Deleting the scenarios in a set $J\subset\{1,\dots,N\}$ and assigning new weights $q_j\ge 0$, $\sum_{j\notin J}q_j=1$, to the kept ones gives $Q=\sum_{j\notin J}q_j\delta_{\omega_j}$. The distance $D(J;q)$ between $P$ and $Q$ is the optimal value of the transportation problem

$$D(J;q)=\min\Bigl\{\sum_{i=1}^N\sum_{j\notin J}c(\omega_i,\omega_j)\eta_{ij}:\ \eta_{ij}\ge 0,\ \sum_{i}\eta_{ij}=q_j,\ \sum_{j\notin J}\eta_{ij}=p_i\Bigr\}.$$

The reduction cost of deleting $J$ is

$$D_J=\sum_{i\in J}p_i\min_{j\notin J}c(\omega_i,\omega_j),$$

and the optimal reduction problem (8) minimizes $D_J$ over all $J$ with $\#J=N-n$, where $n$ is the number of scenarios to keep.

Forward selection builds the kept set greedily. With $J^{[0]}=\{1,\dots,N\}$ and $J^{[i]}=\{1,\dots,N\}\setminus\{u_1,\dots,u_i\}$, it chooses

$$u_i\in\arg\min_{u\in J^{[i-1]}}D_{J^{[i-1]}\setminus\{u\}},\qquad i=1,\dots,n. \tag{16}$$

Fast forward selection (Algorithm 2.4) computes the same choices through an updated cost matrix: $c^{[1]}_{ku}=c(\omega_k,\omega_u)$, $c^{[i]}_{ku}=\min\{c^{[i-1]}_{ku},c^{[i-1]}_{ku_{i-1}}\}$, $z^{[i]}_u=\sum_{k\in J^{[i-1]}\setminus\{u\}}p_kc^{[i]}_{ku}$, and $u_i\in\arg\min_{u\in J^{[i-1]}}z^{[i]}_u$.

Formalization targets

Goal: Theorem 2.5

For $1\le n\le N$ and every run $u_1,\dots,u_n$ of Algorithm 2.4, with any tie-breaking in the arg min,

$$u_i\ \text{satisfies (16)}\quad\text{and}\quad z^{[i]}_{u_i}=D_{J^{[i]}}\qquad(i=1,\dots,n).$$

Milestones

1. Theorem 2.1 (redistribution). For $J$ with at least one kept scenario, $D_J=\min_q D(J;q)$, and the minimum is attained at $\bar q_j=p_j+\sum_{i\in J,\,j(i)=j}p_i$ for every choice of nearest kept scenarios $j(i)$.
2. Eq. (10). $D_{\{1,\dots,N\}\setminus\{u\}}=\sum_{i=1}^Np_ic(\omega_i,\omega_u)$, so (8) with $\#J=N-1$ is problem (10).
3. Eq. (12). The sum $lb$ of the $N-n$ smallest single-deletion costs $p_l\min_{j\neq l}c(\omega_l,\omega_j)$, taken in the greedy order (11), is at most $D_J$ for every $J$ with $\#J=N-n$.
4. Optimality condition (p. 191). If each $l_i$ has a nearest other scenario outside $\{l_1,\dots,l_{N-n}\}\setminus\{l_i\}$, then $\{l_1,\dots,l_{N-n}\}$ solves (8).
5. Eq. (17), unrolled recursion. For any index sequence, $c^{[i]}_{ku}=\min_{j\notin J^{[i-1]}\setminus\{u\}}c(\omega_k,\omega_j)$ for $u\in J^{[i-1]}$.
6. Eq. (17), conclusion. For any index sequence, $z^{[i]}_u=D_{J^{[i-1]}\setminus\{u\}}$ for $u\in J^{[i-1]}$.

Significance

Theorem 2.5 certifies that the cheap update of Algorithm 2.4 (one pairwise minimum per matrix entry and step) produces exactly the greedy forward selection defined through the reduction costs, and that the running objective $z^{[i]}_{u_i}$ is the reduction cost of the scenarios deleted so far. Combined with Theorem 2.1, $z^{[i]}_{u_i}$ is the optimal transportation distance between $P$ and the best measure on the kept scenarios, which is the quantity practitioners monitor to decide how many scenarios to keep. The lower bound (12) and the optimality condition give a posteriori quality certificates for any reduced set.

All results of this mission are proved in the paper or in the works it cites (Dupačová et al., 2003); none is open. To the best of a platform search, none has been machine-checked. The mission provides a verified specification of a widely deployed algorithm, a formal link between a combinatorial set-covering objective and a finite transportation problem, and definitions (reduction cost, transportation plans with a partially free target marginal, greedy runs with arbitrary tie-breaking) reusable by the regular-tree missions of this series and by later scenario-tree construction papers.

Difficulty

The mathematics is elementary; the difficulty is bookkeeping. The recursion for $c^{[i]}$ refers to the previous step's column $u_{i-1}$, which itself was updated, so unrolling it to a minimum over $\{u,u_1,\dots,u_{i-1}\}$ is an induction on $i$ in which the index sets $J^{[i]}$, the 1-based step counter and the complement structure all move together. The natural first attempt, identifying $c^{[i]}_{ku}$ with the minimum over the complement of $J^{[i]}$, is off by one step: the correct set is the complement of $J^{[i-1]}\setminus\{u\}$, which contains $u$ itself. For Theorem 2.1 the lower bound requires using that $c(\omega_i,\omega_i)=0$ for kept scenarios and that every plan ships all of $p_i$ somewhere outside $J$; the attainment part requires constructing the plan explicitly from the choice $j(\cdot)$, including scenarios for which several kept scenarios are equally near.

Formalization scope

• Scenarios are ω : Fin N → E with E a finite-dimensional real normed space; the paper's closed set $\Omega\subset\mathbb R^s$ plays no role beyond containing the scenarios and is omitted. Scenarios need not be distinct.
• $h$ is a function ℝ → ℝ with the paper's assumptions imposed on $[0,\infty)$ (IsGrowthFunction); every theorem carries them, together with $p_i>0$ and $\sum_ip_i=1$.
• The functions $f_0$, $g$ and the stochastic program (1)–(2) that motivate $c$ appear in no statement.
• $D(J;q)$ is the paper's finite transportation problem (p. 188), not the Kantorovich functional on measures. Weights $q$ and plans $\eta$ are indexed by all of $\{1,\dots,N\}$ with entries at deleted indices fixed to $0$.
• $D_J$ requires a proof that the complement of $J$ is nonempty; minima are Finset.inf', never a real infimum with a default value.
• Algorithm 2.4 is a relation on sequences u : ℕ → Fin N with 1-based steps. $c^{[i]}$ is the printed recursion, extended to all indices; runs are any sequences satisfying the arg-min conditions, so every tie-breaking rule is covered.
• The paper's standing restriction $n<N$ is relaxed to $n\le N$ in Theorem 2.5; the statement remains true at $n=N$.
• A trivializing formalization is ruled out: defining $c^{[i]}$ or $z^{[i]}$ directly as the minimum over the selected set or as $D_{J^{[i-1]}\setminus\{u\}}$ would make Theorem 2.5 hold by definition, and proving it for one fixed tie-breaking rule would prove less than the paper; neither is done.
• Proofs of the milestones, alternative proofs of Theorem 2.1 via LP duality, and a verified executable implementation of Algorithm 2.4 are all welcome.

Selected references

• H. Heitsch, W. Römisch, Scenario Reduction Algorithms in Stochastic Programming, Computational Optimization and Applications 24 (2003), 187–206. https://doi.org/10.1023/A:1021805924152
• J. Dupačová, N. Gröwe-Kuska, W. Römisch, Scenario reduction in stochastic programming: an approach using probability metrics, Mathematical Programming 95 (2003), 493–511. https://doi.org/10.1007/s10107-002-0331-0
• S. T. Rachev, W. Römisch, Quantitative stability in stochastic programming: the method of probability metrics, Mathematics of Operations Research 27 (2002), 792–818. https://doi.org/10.1287/moor.27.4.792.304

Motivation

A monotone variational inequality asks for a point at which a monotone operator and a convex function are in equilibrium. It unifies convex minimization (where $F$ is a gradient), convex–concave saddle-point problems (where $F$ is the skew gradient of a Lagrangian), Nash equilibria of monotone games, and complementarity problems in economics and traffic assignment. In operations research, first-order methods for such problems are the workhorse behind large-scale saddle-point formulations of linear and conic programs, where only one operator evaluation and one projection or proximal step per iteration are affordable.

The classical method for Lipschitz monotone operators is Korpelevich's extragradient method (1976) and its proximal variant, Tseng's forward–backward–forward method (2000); both need two evaluations of $F$ per iteration. The reflected projected gradient method of Malitsky (SIAM J. Optim., 2015) uses one evaluation of $F$ but evaluates it at $2z^k-z^{k-1}$, a point that may lie outside the domain of $g$. Malitsky's Golden Ratio Algorithm (GRAAL), introduced in Golden Ratio Algorithms for Variational Inequalities (preprint 2018; published in Mathematical Programming, doi:10.1007/s10107-019-01416-w), uses one evaluation of $F$, always at a feasible point, and one proximal step per iteration. Its fixed-step version, Theorem 1 of that paper, is the subject of this mission; the explicit, adaptive-step version (Theorem 2) is a separate mission of this series.

Setting

Let $\mathcal E$ be a finite-dimensional real inner product space with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\| = \sqrt{\langle\cdot,\cdot\rangle}$. Let $g:\mathcal E\to(-\infty,+\infty]$ and write $\operatorname{dom} g = \{x : g(x)<+\infty\}$. Let $F:\operatorname{dom} g\to\mathcal E$. The variational inequality is

$$\text{find } z^*\in\mathcal E \quad\text{such that}\quad \langle F(z^*), z-z^*\rangle + g(z)-g(z^*)\ \ge\ 0\qquad \forall z\in\mathcal E. \tag{1}$$

The standing assumptions are:

• (C1) the solution set $S$ of (1) is nonempty;
• (C2) $g$ is proper (never $-\infty$, finite somewhere), convex, and lower semicontinuous;
• (C3) $F$ is monotone: $\langle F(u)-F(v),u-v\rangle\ge0$ for all $u,v\in\operatorname{dom} g$.

The proximal operator of $g$ is $\operatorname{prox}_g(z) = \operatorname{argmin}_x\{g(x)+\tfrac12\|x-z\|^2\}$. Let $\varphi = \frac{\sqrt5+1}{2}$ be the golden ratio, so that $\varphi^2 = 1+\varphi$. For a step $\lambda>0$ and arbitrary starting points $z^1,\bar z^0\in\mathcal E$, the Golden Ratio Algorithm generates, for $k\ge1$,

$$\bar z^k = \frac{(\varphi-1)z^k + \bar z^{k-1}}{\varphi},\qquad z^{k+1} = \operatorname{prox}_{\lambda g}\big(\bar z^k - \lambda F(z^k)\big). \tag{6}$$

The first line is a convex combination of the newest iterate and the previous average; the second is a forward–backward step taken from the average rather than from $z^k$.

Formalization targets

Goal: Theorem 1

If $F$ is $L$-Lipschitz on $\operatorname{dom} g$ ($L>0$), (C1)–(C3) hold, and $\lambda\in\big(0,\frac{\varphi}{2L}\big]$, then there is $z^*\in S$ with

$$z^k\to z^*\qquad\text{and}\qquad \bar z^k\to z^*\qquad(k\to\infty).$$

Both sequences converge, to one and the same solution. The goal is stated with the paper's exact step range; no rate is claimed, as the paper claims none.

Milestones

1. Eq. (4), the prox-inequality: for proper convex lsc $g$,

$$\bar x = \operatorname{prox}_g z \iff \langle\bar x - z, x-\bar x\rangle\ge g(\bar x)-g(x)\quad\forall x\in\mathcal E.$$

2. Eq. (12), an identity using only the averaging step of (6): for every point $z^*$,

$$\|z^{k+1}-z^*\|^2 = (1+\varphi)\|\bar z^{k+1}-z^*\|^2-\varphi\|\bar z^k-z^*\|^2+\tfrac1\varphi\|z^{k+1}-\bar z^k\|^2 .$$

3. Eq. (14), the energy inequality: for $z^*\in S$ and $k\ge2$,

$$(1+\varphi)\|\bar z^{k+1}-z^*\|^2+\tfrac\varphi2\|z^{k+1}-z^k\|^2\le(1+\varphi)\|\bar z^k-z^*\|^2+\tfrac\varphi2\|z^k-z^{k-1}\|^2-\varphi\|z^k-\bar z^k\|^2 .$$

4. Lemma 1 (Bauschke–Combettes, Theorem 5.5): a sequence that is Fejér monotone with respect to a nonempty set $C$ and whose cluster points all lie in $C$ converges to a point of $C$.

Significance

The result. Theorem 1 shows that monotone variational inequalities with a Lipschitz operator can be solved with one operator evaluation and one proximal step per iteration, with $F$ evaluated only at points of $\operatorname{dom} g$, where it is defined. This matters when $F$ is expensive (a large matrix–vector product, a simulation) or undefined outside the feasible set (for instance an operator involving $\log x$ on the positive orthant). The analysis also explains the constant: the averaging weight $\varphi$ is the largest $c$ with $1/c\ge c-1$, and the step bound $\varphi/(2L)$ follows from it. The fixed-step analysis is the template for the explicit, adaptive-step EGRAAL of the same paper (Theorem 2), which needs only local Lipschitz continuity of $F$.

The formalization. The theorem has a published proof, and no machine-checked version of it or of GRAAL is known. Mathlib contains the golden ratio, Lipschitz conditions, lower semicontinuity and cluster points, but no proximal operator of an extended-valued function, no prox-inequality and no Fejér-monotonicity convergence lemma. This mission produces those pieces and a complete convergence proof for a first-order VI method, which are reusable for projected gradient, forward–backward, extragradient and reflected-gradient analyses.

Difficulty

The naive approach, to show that $\|z^k-z^*\|$ decreases, fails: GRAAL is not Fejér monotone in $z^k$, because the forward step is taken from the average $\bar z^k$ and uses $F(z^k)$ rather than $F$ at the new point. The quantity that decreases is an energy mixing $\|\bar z^k-z^*\|^2$ with the successive difference $\|z^k-z^{k-1}\|^2$, and both the averaging identity and the Lipschitz estimate must produce matching coefficients for the cross terms to cancel. The energy inequality alone gives only boundedness and vanishing successive differences; convergence of the whole sequence, and the fact that the limit solves (1) when $g$ is merely lower semicontinuous and extended-valued, is a separate step. On the formal side, $g$ takes the value $+\infty$, so the prox-inequality and the variational inequality must be handled in extended arithmetic without letting $\infty-\infty$ decide anything.

Formalization scope

• $\mathcal E$ is a type E with [NormedAddCommGroup E] [InnerProductSpace ℝ E] [FiniteDimensional ℝ E].
• $g$ is E → EReal. (C2) is IsProperConvexLSC g: never $\bot$, somewhere finite, convex epigraph $\{(x,t)\in E\times\mathbb R: g(x)\le t\}$, and LowerSemicontinuous g on all of E. $\operatorname{dom} g$ is effDom g = {x | g x ≠ ⊤}.
• $F$ is a total function E → E; monotonicity and the Lipschitz bound $\|F(u)-F(v)\|\le L\|u-v\|$ are required on effDom g only. The step range is 0 < λ, λ ≤ φ / (2 * L) with 0 < L and φ = Real.goldenRatio.
• $S$ is solutionSet g F: points of effDom g satisfying (1) for every $z\in E$, evaluated in EReal.
• The proximal step is the argmin predicate IsProxPoint (fun x => λ * g x) w z⁺, not a choice function, so no junk value is involved. A run of (6) is IsGRAALRun g F λ z zbar on sequences ℕ → E indexed as in the paper: $z^1$ and $\bar z^0$ are free and the entry $z^0$ is unused.
• The conclusion is ∃ zs ∈ solutionSet g F, Tendsto z atTop (𝓝 zs) ∧ Tendsto zbar atTop (𝓝 zs).

The hypotheses of the goal are jointly satisfiable, so the theorem is not vacuous: for $g\equiv0$ and $F\equiv0$ every point is a solution and constant sequences form a run of (6); a formalization under which IsGRAALRun has no instances, or in which $S$ may be empty, is ruled out. Two hypotheses are added to printed statements and flagged in their notes: $C\neq\emptyset$ in Lemma 1, which is false without it, and $z^1\in\operatorname{dom} g$ in Eq. (14), needed at $k=2$ because the paper's $F$ is only defined on $\operatorname{dom} g$.

Welcome contributions: existence and uniqueness of the proximal point of a proper convex lsc function in finite dimensions; the prox-inequality; Fejér-monotonicity lemmas; the energy inequality; and the final convergence argument. The prox and Fejér infrastructure is independent of the golden ratio and is shared with the second mission of this series.

Selected references

• Y. Malitsky, Golden Ratio Algorithms for Variational Inequalities, preprint, Optimization Online 6598, 2018. https://optimization-online.org/wp-content/uploads/2018/05/6598.pdf ; published in Mathematical Programming. https://doi.org/10.1007/s10107-019-01416-w
• H. H. Bauschke, P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer, 2011 (2nd ed. 2017). https://doi.org/10.1007/978-3-319-48311-5
• G. M. Korpelevich, The extragradient method for finding saddle points and other problems, Ekonomika i Matematicheskie Metody 12 (1976) 747–756.
• P. Tseng, A modified forward–backward splitting method for maximal monotone mappings, SIAM J. Control Optim. 38 (2000) 431–446. https://doi.org/10.1137/S0363012998338806
• Y. Malitsky, Projected reflected gradient methods for monotone variational inequalities, SIAM J. Optim. 25 (2015) 502–520. https://doi.org/10.1137/14097238X

Motivation

Many estimation problems in statistics, signal processing and bioinformatics minimize a composite function $f = g + h$: a smooth convex loss $g$ plus a convex but nonsmooth penalty or constraint $h$, such as the lasso's $\ell_1$ norm or the indicator of a convex set. Proximal Newton-type methods handle such problems by minimizing, at each iterate $x_k$, a model $\hat f_k = \hat g_k + h$ in which $g$ is replaced by its second-order Taylor expansion. Widely used solvers of this kind (glmnet, newGLMNET, QUIC) never solve these model subproblems exactly; they stop an inner iterative solver early by some heuristic. Lee, Sun and Saunders (arXiv:1206.1623v13, 2014) proposed an adaptive stopping rule for the inner solver and proved that it preserves fast local convergence. This mission formalizes that local convergence theory (§3.4 of the paper).

Timeline:

• 1982: Dembo, Eisenstat and Steihaug introduce inexact Newton methods for smooth equations and prove local linear and superlinear convergence under a relative-residual condition with forcing terms $\eta_k$ (doi:10.1137/0719025).
• 1996: Eisenstat and Walker propose self-adjusting forcing terms that avoid oversolving (doi:10.1137/0917003).
• 2012–2014: Lee, Sun and Saunders transfer the relative-residual condition to composite functions, replacing gradients by composite gradient steps, and prove Theorems 3.10 and 3.11.
• 2016: Byrd, Nocedal and Oztoprak analyze inexact proximal Newton methods for $\ell_1$-regularized problems under an additional sufficient-descent condition on the subproblem (doi:10.1007/s10107-015-0941-y).

Setting

Work in $\mathbb R^n$ with the Euclidean inner product. The smooth part $g:\mathbb R^n\to\mathbb R$ is twice continuously differentiable and strongly convex with constant $m>0$: $g(y)\ge g(x)+\nabla g(x)^T(y-x)+\frac m2\|x-y\|^2$ for all $x,y$. Its gradient $\nabla g$ is Lipschitz with constant $L_1$, its Hessian $\nabla^2 g$ is Lipschitz with constant $L_2$, and $\nabla^2 g(x)\preceq MI$ for a constant $M>0$. The nonsmooth part $h$ is proper, closed and convex, and may take the value $+\infty$; it is given by its domain $D$ and its values on $D$. The problem is $\min_x f(x)=g(x)+h(x)$, and $x^\star$ denotes its (unique) optimal solution.

The proximal mapping of $h$ is $\operatorname{prox}_h(v)=\arg\min_y h(y)+\frac12\|y-v\|^2$. The composite gradient step with step length $t>0$ is

$$G_{tf}(x)=\tfrac1t\big(x-\operatorname{prox}_{th}(x-t\nabla g(x))\big),$$

with $G_f=G_{1f}$; it vanishes exactly at minimizers of $f$ and plays the role of the gradient. The step $G_{f/M}$ is the unit step on $f/M=g/M+h/M$. The model at $x_k$ is $\hat f_k=\hat g_k+h$ with $\hat g_k(y)=g(x_k)+\nabla g(x_k)^T(y-x_k)+\frac12(y-x_k)^T\nabla^2 g(x_k)(y-x_k)$.

The inexact proximal Newton method with unit step lengths produces $x_{k+1}=x_k+\Delta x_k$, where the direction $\Delta x_k$ is any point satisfying the adaptive stopping condition

$$\|G_{\hat f_k/M}(x_k+\Delta x_k)\|\le\eta_k\,\|G_{f/M}(x_k)\|\qquad(2.24)$$

for a forcing term $\eta_k\ge0$. The Eisenstat–Walker choice is

$$\eta_k=\min\Big\{\frac m2,\ \frac{\|G_{\hat f_{k-1}/M}(x_k)-G_{f/M}(x_k)\|}{\|G_{f/M}(x_{k-1})\|}\Big\}.\qquad(2.25)$$

Formalization targets

Goal: Theorem 3.10 (p. 17)

1. There are $\bar\eta\in(0,m/2)$, $\delta>0$ and $r\in[0,1)$ such that, whenever $0\le\eta_k\le\bar\eta$ for all $k$ and $\|x_0-x^\star\|<\delta$,

$$\|x_{k+1}-x^\star\|\le r\,\|x_k-x^\star\|\quad\text{for all }k.$$

2. For every forcing sequence with $\eta_k\ge0$, $\eta_k\to0$, there is $\delta>0$ such that every run with $\|x_0-x^\star\|<\delta$ converges to $x^\star$ q-superlinearly: for every $\varepsilon>0$, eventually $\|x_{k+1}-x^\star\|\le\varepsilon\|x_k-x^\star\|$.

Both parts are asserted together. The statement fixes no constant beyond the existence of $\bar\eta$, $\delta$ and $r$.

Milestones

• §2.1, property 3: $G_f(x)=0$ if and only if $x$ minimizes $f$.
• Lemma 2.2: $\|G_f(x)\|\le(L_1+1)\|x-x^\star\|$.
• Lemma 3.8: $\|G_f(x)-G_{\hat f_k}(x)\|\le\frac{L_2}2\|x-x_k\|^2$.
• Lemma 3.9: $(x-y)^T(G_{tf}(x)-G_{tf}(y))\ge\frac m2\|x-y\|^2$ for $0<t\le1/L_1$.
• Theorem 3.11: with the forcing terms (2.25), the method converges q-superlinearly from every start sufficiently close to $x^\star$.

Significance

Theorem 3.10 justifies stopping the inner solver of a proximal Newton method at a relative accuracy that is set by the current optimality measure $\|G_{f/M}(x_k)\|$: a constant small forcing term keeps linear convergence, and forcing terms that decay to zero recover superlinear convergence, with no sufficient-descent condition on the subproblem and for a generic nonsmooth $h$. Theorem 3.11 shows that the self-adjusting choice (2.25) achieves the superlinear regime automatically. Together they are the composite analogue of the inexact Newton theory used in most large-scale smooth solvers.

The results are proved in the paper. No machine-checked proof of any of them is known, and the platform contains no proximal mapping, composite gradient step or inexact Newton condition. A formalization adds a checked proximal-operator toolkit (existence and nonexpansiveness of $\operatorname{prox}$, the optimality characterization of $G_f$, strong monotonicity of $G_{tf}$) and a precise form of the theorem: the paper's proofs mix two scalings of the composite step and cite a lemma where another is meant, so the formal proof settles which constants are valid.

Difficulty

The obvious argument compares the inexact step with the exact proximal Newton step and treats the gap as a perturbation. For composite functions this fails: the exact step is defined by a nonsmooth inclusion, and the stopping condition bounds a residual of the model's composite gradient step, not the distance to the model's minimizer. The link between the two is strong monotonicity of the composite gradient step (Lemma 3.9), which requires controlling the proximal mapping of a general closed convex $h$ jointly with the curvature of $g$; for $h=0$ it is immediate, and for general $h$ it is the central step. A second difficulty is that the threshold on $\eta_k$ is not scale invariant: a threshold below $m/2$ chosen arbitrarily does not give convergence, so the admissible $\bar\eta$ has to come out of the analysis.

Formalization scope

The space is EuclideanSpace ℝ (Fin n). The nonsmooth part is a pair $(D,h)$: $D$ nonempty and convex, $h$ convex on $D$, and the extended function ($h$ on $D$, $+\infty$ off $D$) lower semicontinuous; indicator functions of closed convex sets are included. The proximal mapping is a total function chosen among the minimizers over $D$, which exist uniquely under these hypotheses. $G_{f/M}$, $G_{\hat f_k/M}$ and $G_{tf}$ are functions of the split $(g,D,h)$ and a scalar, never of $f$ alone. The Hessian is the derivative of the gradient map, measured in operator norm. Sequences are indexed from $k=0$; a run requires $x_0\in D$ and $x_k+\Delta x_k\in D$. Rates are stated without quotients.

"$x_0$ sufficiently close to $x^\star$" is an existential radius chosen before the run; assuming $x_k\to x^\star$, letting the radius depend on the run, or reading part 1 as "for every $\bar\eta<m/2$" (which is false: $g(x)=2x^2$, $h=0$, $\eta_k\equiv\frac32$ diverges) are ruled out. The forcing sequence of Theorem 3.10 is fixed in advance; that of Theorem 3.11 depends on the iterates through (2.25), with a free first term $\eta_0\in[0,m/2]$.

A complete development needs existence, uniqueness and firm nonexpansiveness of the proximal mapping of an extended-valued closed convex function, the subgradient characterization of $\operatorname{prox}$, and a second-order Taylor bound for $C^2$ functions with Lipschitz Hessian. These are reusable well beyond this mission. Contributions of any of them, of the milestones, or of alternative proofs of Lemma 3.9 are welcome.

Selected references

• J. D. Lee, Y. Sun, M. A. Saunders, Proximal Newton-type methods for minimizing composite functions, arXiv:1206.1623v13, 2014; SIAM J. Optim. 24(3), 2014. https://arxiv.org/abs/1206.1623
• R. S. Dembo, S. C. Eisenstat, T. Steihaug, Inexact Newton methods, SIAM J. Numer. Anal. 19(2), 1982. https://doi.org/10.1137/0719025
• S. C. Eisenstat, H. F. Walker, Choosing the forcing terms in an inexact Newton method, SIAM J. Sci. Comput. 17(1), 1996. https://doi.org/10.1137/0917003
• R. H. Byrd, J. Nocedal, F. Oztoprak, An inexact successive quadratic approximation method for L-1 regularized optimization, Math. Program. 157, 2016. https://doi.org/10.1007/s10107-015-0941-y

Motivation

The maximum-flow problem asks how much of a commodity can be sent from a source to a sink through a network whose edges carry capacities. It is a basic model of operations research. Transportation, scheduling, bipartite matching and image segmentation reduce to it, and it is the inner step of many combinatorial algorithms.

Goldberg and Tarjan introduced the push–relabel (preflow) method in A New Approach to the Maximum-Flow Problem (J. ACM 35(4), 1988). Ford–Fulkerson-type algorithms augment along whole source–sink paths. The push–relabel method instead moves excess flow across single edges, guided by integer distance labels on the vertices. Whatever order its local operations are applied in, it is correct and performs $O(n^2 m)$ of them (§3 of the paper). Section 4 shows that one particular order, processing the active vertices first-in, first-out, cuts the dominant term, the number of nonsaturating pushes, to $O(n^3)$. The method and its FIFO and highest-label variants are the standard practical maximum-flow codes.

Timeline:

• 1956: Ford and Fulkerson, augmenting paths and max-flow min-cut.
• 1970–72: Dinic, and Edmonds and Karp, give polynomial augmenting-path bounds.
• 1974: Karzanov introduces preflows and obtains $O(n^3)$.
• 1982: Shiloach and Vishkin give a parallel $O(n^2 \log n)$ preflow algorithm with a first-in, first-out flavour.
• 1988: Goldberg and Tarjan, the generic push–relabel method, the FIFO bound of this mission, and $O(nm \log(n^2/m))$ with dynamic trees.

Setting

A flow network has a finite vertex set $V$ with $n = |V|$, a capacity $c(v,w) \ge 0$ on every ordered pair, a source $s$ and a sink $t \neq s$. The edges are the pairs with $c(v,w) > 0$, and there are no loops. A preflow is a function $f$ on vertex pairs with $f(v,w) \le c(v,w)$ and $f(v,w) = -f(w,v)$. Its excess $e(v) = \sum_u f(u,v)$ must be nonnegative at every $v \neq s$. The residual capacity is $r_f(v,w) = c(v,w) - f(v,w)$. A labeling $d$ assigns each vertex a value in $\mathbb{N} \cup \{\infty\}$. A vertex $v \notin \{s,t\}$ is active if $d(v) < \infty$ and $e(v) > 0$.

The two basic operations (Fig. 1 of the paper) are:

• push$(v,w)$, applicable when $v$ is active, $r_f(v,w) > 0$ and $d(v) = d(w)+1$. It sends $\delta = \min(e(v), r_f(v,w))$ from $v$ to $w$. The push is saturating if $r_f(v,w) = 0$ afterwards and nonsaturating otherwise.
• relabel$(v)$, applicable when $v$ is active and $d(v) \le d(w)$ for every residual edge $(v,w)$. It sets $d(v) \leftarrow \min\{d(w)+1 : r_f(v,w) > 0\}$.

The algorithm starts by saturating every edge leaving $s$, with $d(s) = n$ and $d(v) = 0$ for $v \neq s$.

In the first-in, first-out algorithm (§4), each vertex $v$ scans a fixed list $L(v)$ of its neighbours through a current edge. The push/relabel$(v)$ operation pushes through the current edge if possible. Otherwise it advances the current edge, or, at the end of the list, returns to the first edge and relabels $v$. Active vertices wait in a queue $Q$, initially $\{v \in V - \{s,t\} : c(s,v) > 0\}$. The discharge operation removes the front vertex $v$ and repeats push/relabel$(v)$ until $e(v) = 0$ or $d(v)$ increases. Every vertex that becomes active meanwhile is appended to $Q$, and $v$ is appended too if it is still active. Passes over the queue are defined inductively. Pass 1 consists of the discharges of the initially queued vertices. Pass $i+1$ consists of the discharges of vertices added during pass $i$.

Formalization targets

Goal: Corollary 4.4 (p. 931)

For every network, every edge-list order, every initial queue order, and every run of the FIFO algorithm,

$$\#\{\text{nonsaturating pushes}\} \le 4n^3 .$$

The constant is the printed one.

Milestones

• Lemma 4.1 (p. 929): the push/relabel operation relabels only when relabeling is applicable.
• Lemma 3.5 (p. 926): from any vertex with positive excess, the source is reachable in the residual graph.
• Lemma 3.7 (p. 927): at any time, $d(v) \le 2n-1$ for every vertex.
• Lemma 3.8 (p. 927): at most $2n-1$ relabelings per vertex and at most $(2n-1)(n-2) < 2n^2$ in total.
• Lemma 4.3 (p. 930): at most $4n^2$ passes over the queue.

Significance

Corollary 4.4 is the combinatorial core of Theorem 4.5, which states that the FIFO algorithm runs in $O(n^3)$ time. Theorem 4.2 shows that the remaining work of the implementation is $O(nm)$ plus constant time per nonsaturating push. The bound of Corollary 4.4 is therefore what separates the $O(n^3)$ FIFO method from the $O(n^2 m)$ bound of the generic method, which matters on dense networks. The same pass-counting argument is reused for the parallel algorithm of §6 and underlies later analyses of highest-label and wave variants.

The results are proved in the paper. Formalizing them adds an analysis of a push–relabel algorithm, which the platform does not yet have. Its existing network-flow material states max-flow min-cut and Ford–Fulkerson termination in an arc-based model with nonnegative flows (the Introduction to Linear Optimization missions). The mission builds a precise operational model of the FIFO implementation, with edge lists, current edges and a queue carrying pass numbers, and states an explicit operation count for it. A companion mission in this series treats the generic algorithm's correctness and its $(2n-1)(n-2) + 2nm + 4n^2m$ operation bound.

Difficulty

The obvious argument is the potential-function count of §3, over the sum of the labels of active vertices. It yields only $4n^2 m$ and does not use the queue discipline at all. The $4n^3$ bound has to charge nonsaturating pushes to passes over the queue, and then bound the number of passes by the total growth of the labels. Neither step is visible in the generic algorithm, because both depend on the order in which vertices are processed.

Making this rigorous requires invariants of the implementation that the paper uses silently:

• a vertex is in the queue exactly when it is active, and at most once;
• pass numbers are nondecreasing along the queue;
• current edges only move forward between relabelings.

Lemma 4.1 in particular depends on the current-edge scan: an edge passed over earlier stays inadmissible until $v$ is relabeled.

Formalization scope

The Lean development works in namespace GoldbergTarjan.FIFO. Vertices form a type V with [Fintype V] [DecidableEq V], and $n$ is Fintype.card V. Capacities are c : V → V → ℝ with c ≥ 0 and c v v = 0. Flows are antisymmetric real functions on all ordered pairs, not nonnegative arc flows. Excess is computed from the flow. Labels are in ℕ∞, and the empty minimum in relabel is ⊤.

The state of the algorithm consists of the flow, the labels, the current-edge index cur v into the edge list L v, and the queue Q : List (V × ℕ), each entry tagged with its pass number. Push/relabel (Fig. 3) is a total function, and a discharge (Fig. 4) is a relation carrying the number of push/relabel operations it performs. A run consists of the states S 0, …, S K with S 0 the initial state and consecutive states related by one discharge. The printed variant of Fig. 4, which stops as soon as $v$ is relabeled, is the one formalized. Counts are natural numbers over all push/relabel operations of all discharges. The number of passes is the largest pass tag of a discharged entry.

All constants are explicit, exactly as printed:

• $2n-1$ (Lemmas 3.7, 3.8);
• $(2n-1)(n-2)$ and $2n^2$ (Lemma 3.8);
• $4n^2$ (Lemma 4.3);
• $4n^3$ (Corollary 4.4).

No asymptotic notation is used, and no $m \ge n-1$ assumption is made.

A model without current edges, where relabeling happens whenever no push applies, would make Lemma 4.1 vacuous and change the algorithm. Pass numbers that are not propagated by the "added during pass $i$" rule would make the pass count arbitrary. Both are ruled out by the definitions. A sorry-free check, outside the proposal, exhibits a three-vertex network with two legal discharges, two passes and no nonsaturating push, so the run hypotheses are satisfiable.

Reusable beyond this mission are the network, preflow, push and relabel definitions and Lemma 3.5, which is about an arbitrary preflow. Contributions welcome: invariants of FIFO runs (preflow, valid labeling, queue = active set, cur within bounds), proofs of the milestones, and the reduction of Corollary 4.4 to Lemma 4.3.

Selected references

• A. V. Goldberg, R. E. Tarjan, A New Approach to the Maximum-Flow Problem, Journal of the ACM 35(4):921–940, 1988. https://doi.org/10.1145/48014.61051
• A. V. Karzanov, Determining the maximal flow in a network by the method of preflows, Soviet Math. Doklady 15:434–437, 1974.
• Y. Shiloach, U. Vishkin, An O(n² log n) parallel max-flow algorithm, Journal of Algorithms 3(2):128–146, 1982. https://doi.org/10.1016/0196-6774(82)90013-X
• L. R. Ford, D. R. Fulkerson, Maximal flow through a network, Canadian Journal of Mathematics 8:399–404, 1956. https://doi.org/10.4153/CJM-1956-045-5
• J. Edmonds, R. M. Karp, Theoretical improvements in algorithmic efficiency for network flow problems, Journal of the ACM 19(2):248–264, 1972. https://doi.org/10.1145/321694.321699

Motivation

Many estimation problems in statistics, machine learning and signal processing minimize a composite function, the sum of a smooth loss and a convex but nonsmooth regularizer or constraint: the lasso and $\ell_1$-regularized logistic regression, the graphical lasso for sparse inverse covariance estimation, and constrained least squares, where the nonsmooth part is the indicator function of a convex set. First-order proximal gradient methods (ISTA, FISTA, SpaRSA) are the standard tools, and their convergence is at best linear. Practical solvers such as glmnet, newGLMNET and QUIC instead minimize a local quadratic model of the smooth part plus the nonsmooth part at every iteration, and in practice they need far fewer iterations.

Lee, Sun and Saunders (arXiv:1206.1623, SIAM J. Optim. 2014) put these methods into one framework, proximal Newton-type methods, and proved that they inherit the local convergence rates of Newton and quasi-Newton methods for smooth problems. This mission formalizes the exact-subproblem half of their analysis, ending with q-superlinear convergence of proximal quasi-Newton methods whose Hessian approximations satisfy the Dennis–Moré criterion.

Timeline. Dennis and Moré (1974) characterized superlinear convergence of quasi-Newton methods for smooth equations and minimization by what is now called the Dennis–Moré condition. Tseng and Yun (2009) analyzed coordinate gradient descent for composite problems with a scaled quadratic model. Byrd, Nocedal and Oztoprak (2013) studied inexact proximal Newton methods for $\ell_1$-regularized problems. Lee, Sun and Saunders (2012–2014) proved quadratic and superlinear local convergence for a generic closed convex $h$.

Setting

The problem is

$$\min_{x\in\mathbb R^n} f(x) := g(x) + h(x). \qquad (1.1)$$

The smooth part $g:\mathbb R^n\to\mathbb R$ is twice continuously differentiable and strongly convex with constant $m>0$, meaning $g(y)\ge g(x)+\nabla g(x)^T(y-x)+\tfrac m2\|x-y\|^2$ for all $x,y$ (Definition 3.2). Its gradient $\nabla g$ and Hessian $\nabla^2 g$ are Lipschitz continuous with constants $L_1$ and $L_2$. The nonsmooth part $h$ is a proper closed convex function that may take the value $+\infty$. Its effective domain $D=\operatorname{dom} h$ is nonempty and convex, and $x^\star$ denotes the optimal solution of (1.1), which is unique by strong convexity.

At an iterate $x_k$ the method chooses a symmetric positive definite matrix $H_k$ and computes the search direction $\Delta x_k$, the minimizer of the model subproblem

$$\Delta x_k=\arg\min_d\ \nabla g(x_k)^Td+\tfrac12 d^TH_kd+h(x_k+d). \qquad (2.9)$$

The predicted decrease is $\lambda_k=\nabla g(x_k)^T\Delta x_k+h(x_k+\Delta x_k)-h(x_k)$. A step length $t$ satisfies the sufficient descent condition (2.19) if $f(x_k+t\Delta x_k)\le f(x_k)+\alpha t\lambda_k$ for a fixed $\alpha\in(0,\tfrac12)$. A backtracking line search with factor $\beta\in(0,1)$ takes $t_k=\beta^{j}$ for the least $j\ge0$ that passes, so the unit step is tried first. The update is $x_{k+1}=x_k+t_k\Delta x_k$ (Algorithm 1). With $H_k=\nabla^2 g(x_k)$ this is the proximal Newton method. With any other choice of $H_k$ it is a proximal quasi-Newton method. The sequence $\{H_k\}$ satisfies the Dennis–Moré criterion if

$$\frac{\|(H_k-\nabla^2 g(x^\star))(x_{k+1}-x_k)\|}{\|x_{k+1}-x_k\|}\to0. \qquad (3.2)$$

Formalization targets

Goal: Theorem 3.7

If $mI\preceq H_k\preceq MI$ for all $k$, with $0<m\le M$, and $\{H_k\}$ satisfies (3.2), then every run of Algorithm 1 from any $x_0\in D$ satisfies

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

The goal fixes no rate constant. It asserts only the shape of the convergence.

Milestones

In the order the proof uses them:

1. Proposition 2.4: $\lambda\le-\Delta x^TH\Delta x$ and $f(x+t\Delta x)\le f(x)+t\lambda+O(t^2)$.
2. Proposition 2.5: $x$ is optimal if and only if $\Delta x=0$ at $x$.
3. Lemma 2.6: every $t\le\min\{1,(2m/L_1)(1-\alpha)\}$ satisfies (2.19).
4. Theorem 3.1 (global convergence), restated under the assumptions of §3.3: $x_k\to x^\star$.
5. Lemma 3.3: the proximal Newton method eventually accepts the unit step.
6. Theorem 3.4: the proximal Newton method converges q-quadratically, with eventually

$$\|x_{k+1}-x^\star\|\le\frac{L_2}{2m}\|x_k-x^\star\|^2 .$$

7. Lemma 3.5 / A.1: under (3.2) the unit step is eventually accepted.
8. Proposition 3.6: $\|\Delta x_1-\Delta x_2\|\le\sqrt{(1+\bar\theta)/m}\,\|(H_2-H_1)\Delta x_1\|^{1/2}\|\Delta x_1\|^{1/2}$, with $\bar\theta$ depending only on the eigenvalue bounds.

Significance

The result. Theorem 3.7 is the composite counterpart of the Dennis–Moré theorem. It says that the rate of a proximal quasi-Newton method is governed by how well $H_k$ approximates the Hessian of the smooth part along the steps actually taken, whatever the nonsmooth part is. It covers proximal BFGS-type methods for $\ell_1$-regularized and constrained problems, and it explains why solvers built on these methods reach high accuracy in few iterations. Theorem 3.4 gives the corresponding quadratic rate when the exact Hessian is used.

Formalizing it. The results are proved in the paper. None of them has a machine-checked proof: the platform currently has Newton's method only for smooth objectives (Boyd–Vandenberghe's quadratic phase in the mission Convex Optimization V: Newton's Method), and nothing on proximal or composite Newton-type methods. The formalization also settles two defects of the printed text. Theorem 3.1 is false as printed, since it lacks an upper bound on $H_k$: with $g(x)=x^2/2$, $h=0$ and $H_k=2^{k+1}$ the iterates stall at about $0.289\,x_0$. It is therefore stated here under the assumptions of §3.3. Proposition 3.6 uses an undefined constant $m_1$ (read as $m$), and the first-order inequalities in its printed proof contain a typo. The explicit constant $L_2/(2m)$ in Theorem 3.4 is the one the paper's proof derives.

Difficulty

The difficulty is the nonsmooth part. For smooth $g$ the Newton step solves a linear system, and the classical analysis works with that closed form. Here $\Delta x_k$ is defined only as the minimizer of a nonsmooth subproblem, and every estimate on it has to come from the optimality of that minimizer, i.e. from the firm nonexpansiveness of scaled proximal maps in a norm that changes with $H_k$. A natural first idea is to apply the smooth Dennis–Moré argument to $\nabla f$. It fails because $f$ is not differentiable and may be $+\infty$ outside $D$. The superlinear rate also depends on the line search eventually accepting the unit step. That acceptance comes only from a third-order Taylor bound combined with (2.15) and the Dennis–Moré residual, and a line search that may return any admissible step does not give it.

Formalization scope

• Space. The space is EuclideanSpace ℝ (Fin n), and matrices are continuous linear operators. $\nabla g$ is Mathlib's gradient, and $\nabla^2 g$ is fderiv ℝ (gradient g). "Positive definite" includes symmetry, and $mI\preceq H\preceq MI$ is stated through quadratic forms of a symmetric $H$.
• The nonsmooth part. $h$ is encoded by its domain $D$ and its values on $D$: $D$ is nonempty and convex, $h$ is convex on $D$, and the $+\infty$-extension of $h$ is lower semicontinuous. The objective $f$ is extended-valued. Only comparisons are made in EReal, never arithmetic. Replacing $h$ by a real-valued function on all of $\mathbb R^n$ would exclude indicator functions and is not the paper's setting.
• Algorithm. The search direction is a predicate: it minimizes (2.9) over $\{d : x+d\in D\}$. Backtracking is the least-$j$ rule with factor $\beta\in(0,1)$, the convention of Boyd and Vandenberghe, whom the paper cites for its line search. Runs are infinite and indexed from $k=0$, with no stopping test.
• Rates. $o(\cdot)$ and (3.2) are stated without quotients: for every $\varepsilon>0$ the inequality holds eventually.
• Ruling out trivial versions. A line search allowed to return any step satisfying (2.19) would make Theorems 3.4 and 3.7 false. Taking $x^\star$ to be an arbitrary point instead of the minimizer, or letting $\bar\theta$ in Proposition 3.6 depend on the data, would empty the statements. None of these readings is used.
• Infrastructure. A complete development needs: existence and uniqueness of minimizers of strongly convex, lower semicontinuous extended functions; first-order optimality for the subproblem; firm nonexpansiveness of scaled proximal maps in the $H$-norm; and second- and third-order Taylor bounds from Lipschitz derivatives. These pieces are reusable across proximal methods. Contributions of any milestone, of these supporting lemmas, or of the goal directly are welcome.

Selected references

• J. D. Lee, Y. Sun, M. A. Saunders, Proximal Newton-type methods for minimizing composite functions, arXiv:1206.1623v13 (2014); SIAM J. Optim. 24(3), 2014. https://arxiv.org/abs/1206.1623
• J. E. Dennis, J. J. Moré, A characterization of superlinear convergence and its application to quasi-Newton methods, Math. Comp. 28 (1974), 549–560. https://doi.org/10.1090/S0025-5718-1974-0343581-1
• P. Tseng, S. Yun, A coordinate gradient descent method for nonsmooth separable minimization, Math. Program. 117 (2009), 387–423. https://doi.org/10.1007/s10107-007-0170-0
• R. H. Byrd, J. Nocedal, F. Oztoprak, An inexact successive quadratic approximation method for convex L-1 regularized optimization, Math. Program. 157 (2016), 375–396; arXiv:1309.3529. https://arxiv.org/abs/1309.3529
• S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004. https://web.stanford.edu/~boyd/cvxbook/

Motivation

The maximum-flow problem asks how much of a commodity can be sent from a source to a sink through a network whose edges have capacities. It is a basic model in operations research (transportation, scheduling, bipartite matching) and a standard subroutine in combinatorial optimization.

Classical algorithms, from Ford and Fulkerson (1956) through Edmonds–Karp and Dinic (1970–1972) and Karzanov (1974), increase a feasible flow along augmenting paths or blocking flows. Goldberg and Tarjan, A New Approach to the Maximum-Flow Problem (J. ACM 35(4), 1988, doi:10.1145/48014.61051), replaced this global view by a local one: the push-relabel method maintains a preflow, which may violate conservation at intermediate vertices, and moves excess along edges toward vertices with smaller distance labels. The generic method, with the basic operations applied in any order, is the starting point of the FIFO, highest-label and dynamic-tree implementations analysed later in the same paper, and push-relabel codes remain among the fastest practical maximum-flow solvers.

This mission formalizes §2–§3 of the paper: the generic algorithm is correct, and it stops after a number of basic operations bounded by an explicit polynomial in the numbers of vertices and edges, whatever order of operations is chosen.

Setting

A flow network has a finite vertex set $V$ with $n = |V|$, a source $s$ and a sink $t \ne s$, and a capacity $c(v,w) \ge 0$ for every ordered pair of vertices, positive exactly on the edges $E = \{(v,w) : c(v,w) > 0\}$; $m = |E|$, and there are no loops, $c(v,v) = 0$.

Flows are real functions on all vertex pairs. A function $f$ satisfies the capacity constraint if $f(v,w) \le c(v,w)$ and antisymmetry if $f(v,w) = -f(w,v)$ for all pairs. The excess of $v$ is $e(v) = \sum_{u} f(u,v)$. A flow also has $e(v) = 0$ for $v \notin \{s,t\}$; a preflow only $e(v) \ge 0$ for $v \ne s$. The value of a flow is $|f| = \sum_v f(v,t)$, and a maximum flow is a flow of maximum value.

The residual capacity is $r_f(v,w) = c(v,w) - f(v,w)$; pairs with $r_f(v,w) > 0$ are the edges of the residual graph $G_f$. A valid labeling is $d : V \to \mathbb{N} \cup \{\infty\}$ with $d(s) = n$, $d(t) = 0$ and $d(v) \le d(w) + 1$ on every residual edge. A vertex $v$ is active if $v \notin \{s,t\}$, $d(v) < \infty$ and $e(v) > 0$.

The two basic operations (Fig. 1 of the paper) are:

• Push$(v,w)$, applicable when $v$ is active, $r_f(v,w) > 0$ and $d(v) = d(w)+1$: send $\delta = \min(e(v), r_f(v,w))$, i.e. $f(v,w) \mathrel{+}= \delta$, $f(w,v) \mathrel{-}= \delta$. It is saturating if $r_f(v,w) = 0$ afterwards and nonsaturating otherwise.
• Relabel$(v)$, applicable when $v$ is active and $d(v) \le d(w)$ for every residual edge $(v,w)$: set $d(v) \leftarrow \min\{d(w)+1 : (v,w) \in E_f\}$ ($\infty$ if there is none).

The generic algorithm (Fig. 2) starts from the preflow that saturates every edge leaving $s$ and is zero elsewhere, with the simple labeling $d(s) = n$, $d(v) = 0$ otherwise, and applies applicable basic operations in any order while one exists. An execution with $K$ basic operations is a sequence of states $(f_0,d_0),\dots,(f_K,d_K)$ from the initial state, each obtained from the previous one by one applicable operation.

Formalization targets

Goal: Theorems 3.11 and 3.4

Assume the paper's standing assumption $m \ge n-1$. For every execution with $K$ basic operations,

$$K \le (2n-1)(n-2) + 2nm + 4n^2 m,$$

and if no basic operation applies in the final state, then $f_K$ is a maximum flow. The paper states the bound as $O(n^2m)$ and proves it as "immediate from Lemmas 3.8, 3.9, and 3.10"; the goal states the sum of those three printed bounds. Since every execution is this short, no order of operations runs forever.

Milestones

In the order the proof uses them: Lemma 2.1 (at an active vertex a push or a relabel applies); Lemma 3.1 (the labeling stays valid); Theorem 3.2 (Ford–Fulkerson: a flow is maximum iff $t$ is unreachable from $s$ in $G_f$); Lemma 3.3 (under a valid labeling $t$ is unreachable from $s$); Lemma 3.5 (from any vertex with positive excess, $s$ is reachable); Lemma 3.6 (labels never decrease; a relabeling increases the label); Lemma 3.7 ($d(v) \le 2n-1$ throughout); Theorem 3.4 (termination with finite labels gives a maximum flow); Lemma 3.8 ($\le 2n-1$ relabelings per vertex, $\le (2n-1)(n-2) < 2n^2$ in total); Lemma 3.9 ($\le 2nm$ saturating pushes); Lemma 3.10 ($\le 4n^2m$ nonsaturating pushes, under $m \ge n-1$). A further, non-milestone item states the unnumbered invariant that every $f_k$ is a preflow.

Significance

The generic bound shows that push-relabel terminates in a polynomial number of steps without any rule for choosing the next operation; the specific orderings of §4–§5 of the paper (first-in first-out, $O(n^3)$; dynamic trees, $O(nm\log(n^2/m))$) refine only the count of nonsaturating pushes, and reuse Lemmas 3.1–3.9 unchanged. The correctness argument, a valid labeling excludes augmenting paths, is the template for the push-relabel minimum-cost flow and assignment algorithms that followed.

These results are proved in the paper and are textbook material. Their machine-checked counterparts are, as far as is known here, not on the Prove2Me platform: the platform's network-flow statements (from Introduction to Linear Optimization, e.g. LinearOptimization.max_flow_min_cut) use a different model, with arc-indexed nonnegative flows and extended-real capacities, and contain nothing about preflows, labels or operation counts. This mission produces a formal account of the antisymmetric-flow model, of Ford–Fulkerson in that model, and of the amortized counting arguments, with the constants the paper prints.

Difficulty

The correctness half is short once the invariants are in place; the difficulty is in the counting. The label bound (Lemma 3.7) is a statement about the whole execution, and it depends on a structural fact about preflows (Lemma 3.5) whose truth rests on antisymmetry and on the nonnegativity of excesses. The obvious first idea for the push counts, bounding pushes per edge or per vertex locally, fails for nonsaturating pushes: flow pushed across a pair can be pushed back later, and nothing local limits how often this happens, so Lemma 3.10 holds only as an amortized statement over the entire execution and depends on both earlier counts. Saturating pushes on a pair can also recur, in both directions, and Lemma 3.9 has to control the interaction between the two directions.

Formally, all of this is reasoning about arbitrary interleavings of operations, with labels in $\mathbb{N} \cup \{\infty\}$ and real-valued flows.

Formalization scope

• Vertices form a finite type with decidable equality; $n$ is its cardinality, $s \ne t$, so $n \ge 2$ and the natural-number subtractions $2n-1$ and $n-2$ are exact. Capacities are a real function on all pairs, nonnegative, zero on the diagonal; $E$ is its support and $m$ its cardinality.
• Flows and preflows are antisymmetric real functions on all pairs (not nonnegative arc flows); the excess is computed from $f$, never stored. A maximum flow is a flow whose value is at least that of every flow.
• Labels live in ℕ∞, with $\infty + 1 = \infty$; the relabel value is an infimum, which is $\infty$ on the empty set.
• An execution is a sequence of states σ : ℕ → State V with a length $K$, starting at the Fig. 2 state with the simple labeling (the paper's own assumption for its proofs), each step an applicable push or relabel. "Terminates" means that no basic operation applies, the loop guard of Fig. 2. The three counts are cardinalities of the sets of step indices of each kind.
• Explicit constants: $2n-1$ per-vertex relabelings, $(2n-1)(n-2) < 2n^2$ total relabelings, $2nm$ saturating pushes, $4n^2m$ nonsaturating pushes, label bound $2n-1$, and the total $(2n-1)(n-2)+2nm+4n^2m$. The standing assumption $m \ge n-1$ appears only on Lemma 3.10 and the goal.
• A trivializing formalization is ruled out: the step relation fixes the pushed amount $\delta = \min(e(v), r_f(v,w))$ and the new label exactly as in Fig. 1, termination is the loop guard rather than "the result is a flow", and a sorry-free check exhibits a concrete network $s \to a \to t$ with a two-step execution (relabel $a$, then push $(a,t)$), so the run hypotheses are satisfiable.

Welcome contributions: proofs of the invariants (preflow, valid labeling, label monotonicity), of Ford–Fulkerson for antisymmetric flows (reusable beyond this mission), and of the counting lemmas. The FIFO bound of §4 is the subject of a companion mission.

Selected references

• A. V. Goldberg, R. E. Tarjan, A New Approach to the Maximum-Flow Problem, Journal of the ACM 35(4):921–940, 1988. doi:10.1145/48014.61051
• L. R. Ford, D. R. Fulkerson, Flows in Networks, Princeton University Press, 1962.
• J. Edmonds, R. M. Karp, Theoretical improvements in algorithmic efficiency for network flow problems, Journal of the ACM 19(2):248–264, 1972. doi:10.1145/321694.321699
• R. K. Ahuja, T. L. Magnanti, J. B. Orlin, Network Flows: Theory, Algorithms, and Applications, Prentice Hall, 1993.