Operations Research

Motivation

A variational inequality asks for a point $x$ of a set $K\subseteq\mathbb R^n$ at which a vector field $f$ points "into" $K$: $(x'-x)^T f(x)\ge 0$ for every $x'\in K$. It is the common form of the first-order optimality conditions of constrained optimization, of complementarity problems, and of equilibrium models in economics and traffic networks. In many of these models the feasible set itself depends on the decision: the admissible actions of one agent are restricted by the current state, as in the impulse-control problems of Bensoussan and Lions that motivated quasi-variational inequalities, where $K=K(x)$.

D. Chan and J. S. Pang, The generalized quasi-variational inequality problem (Math. Oper. Res. 7 (1982) 211–222), unify the quasi-variational inequality with the generalized (set-valued) variational inequality of Fang and Peterson (JOTA 1982). Their §§3–4 prove existence by fixed-point theorems for set-valued maps; §5 takes a different route and characterizes solutions as fixed points of a composite projection map. Theorem 5.3, the subject of this mission, gives conditions under which that map is a contraction, so that its fixed point exists, is unique, solves the problem, and is computed by plain fixed-point iteration from any starting point. It is the algorithmic result of the paper, and an early instance of the projection methods for strongly monotone quasi-variational inequalities studied since (e.g. Nesterov and Scrimali 2011).

Setting

Throughout, $\mathbb R^n$ carries the Euclidean inner product $x^T y$ and norm $\|x\|$.

Given point-to-set mappings $K$ and $f$ of $\mathbb R^n$ into itself, the generalized quasi-variational inequality problem $\mathrm{GQVI}(K,f)$ is to find vectors $x$ and $y$ with

$$x\in K(x),\qquad y\in f(x),\qquad (x'-x)^T y\ge 0\quad\text{for all }x'\in K(x).$$

When $f$ is point-to-point, $f(x)$ is read as the singleton $\{f(x)\}$.

For a set $S$ and a point $z$, the projection $P_S(z)$ is the point of $S$ nearest to $z$, $P_S(z)=\operatorname{sol}\min_{x\in S}\|x-z\|$; it exists and is unique when $S$ is nonempty, closed and convex.

Theorem 5.3 concerns the special structure in which the feasible set moves by translation: fix a nonempty closed convex set $\tilde K$ and a point-to-point mapping $m$, and put

$$K(x)=m(x)+\tilde K=\{m(x)+k : k\in\tilde K\}.$$

For a step length $\lambda>0$ and a point-to-point $f$, the projection map is

$$F_\lambda(x)=P_{K(x)}\bigl(x-\lambda f(x)\bigr).$$

The mappings $m$ and $f$ are assumed Lipschitz continuous with constants $\alpha$, $\beta$ ($\|m(x)-m(y)\|\le\alpha\|x-y\|$, $\|f(x)-f(y)\|\le\beta\|x-y\|$) and strongly monotone with constants $\gamma$, $\delta$ ($(x-y)^T(m(x)-m(y))\ge\gamma\|x-y\|^2$, $(x-y)^T(f(x)-f(y))\ge\delta\|x-y\|^2$).

Formalization targets

Goal: Theorem 5.3 (p. 221)

For each $\lambda>0$ with

$$\lambda^2\beta^2+2\lambda(\alpha\beta-\delta)-2(\gamma-\alpha)<0,$$

the map $F_\lambda$ is a contraction (Lipschitz with a constant $c<1$ independent of the points), it has a fixed point $\tilde x_\lambda$, the point $\tilde x_\lambda$ solves $\mathrm{GQVI}(K,f)$, and the iterates $x^{k+1}=F_\lambda(x^k)$ converge to $\tilde x_\lambda$ from every initial vector $x^0\in\mathbb R^n$. All four conclusions are stated together.

Milestones

1. Projection onto a translate (§5, proof of Theorem 5.3, first display, p. 221): $P_{K(x)}(y)=m(x)+P_{\tilde K}(y-m(x))$ for all $x,y$.
2. Lipschitz estimate (§5, proof of Theorem 5.3, last display, p. 221): for every $\lambda>0$,

$$\|F_\lambda(y^1)-F_\lambda(y^2)\|\le\Bigl[\alpha+\bigl(\lambda^2\beta^2+2\lambda(\alpha\beta-\delta)+(1+\alpha^2-2\gamma)\bigr)^{1/2}\Bigr]\|y^1-y^2\|.$$

3. Theorem 5.1 (p. 220): if every $K(x)$ is closed and convex, $(x^*,y^*)$ solves $\mathrm{GQVI}(K,f)$ if and only if $x^*=P_{K(x^*)}(x^*-y^*)$ and $y^*\in f(x^*)$.

Significance

The result. Theorem 5.3 turns an existence question into a computation: under Lipschitz and strong monotonicity assumptions, a quasi-variational inequality with translated feasible sets has exactly one solution reachable by projection iterations, each of which is a projection on the fixed set $\tilde K$ (a convex quadratic program when $\tilde K$ is polyhedral). The step-size window it gives is explicit in $\alpha,\beta,\gamma,\delta$, so it certifies a convergent method before any iteration is run. The closing remark of the paper (p. 222) reads each step as solving the GQVI under a zero-th order approximation of $K$, the viewpoint behind later splitting methods.

Formalizing it. The result is proved in the paper; the proof is short, but its constants and the equivalence of the two contraction conditions are easy to get wrong. A machine-checked version fixes the exact hypotheses (no sign conditions on the constants, Euclidean geometry), and produces reusable pieces: the translation identity for projections, nonexpansiveness of the Euclidean projection on a closed convex set, and the projection characterization of quasi-variational inequalities (Theorem 5.1).

Difficulty

Banach's fixed-point theorem does the last step; the work is the estimate. The naive bound, projection nonexpansiveness applied directly to $F_\lambda$, fails because the sets $K(y^1)$ and $K(y^2)$ differ: two projections on different sets are not controlled by the distance of the projected points alone. The translation identity separates the moving part $m(y^1)-m(y^2)$ from a projection on the one set $\tilde K$, at the cost of the additive term $\alpha$ in the constant. The remaining square must be expanded with the inner-product cross terms bounded by the monotonicity constants in the right directions, including the cross term between $f$ and $m$. Finally, the condition "bracket $<1$" is equivalent to the stated $\lambda$-condition only when $\alpha<1$, which must be derived from the hypotheses rather than assumed.

Formalization scope

• The space is EuclideanSpace ℝ (Fin n), with Mathlib's Euclidean norm and inner product; Fin n → ℝ (sup norm) would change every constant. No assumption $n\ge1$ is made; at $n=0$ all statements hold trivially.
• The constants $\alpha,\beta,\gamma,\delta$ are real numbers with no sign conditions, as in the paper; the Lipschitz and monotonicity hypotheses are the displayed inequalities for all $x,y$. For $n\ge1$ they force $\alpha,\beta\ge0$, $\gamma\le\alpha$, $\delta\le\beta$, and the $\lambda$-condition then forces $\alpha<1$ and $\delta>\alpha\beta$.
• $K(x)$ is the translate $\{m(x)+k : k\in\tilde K\}$ of a fixed set $\tilde K$, assumed nonempty, closed and convex. The projection is a nearest-point function proj that returns a junk value only when no nearest point exists; under the hypotheses of every statement using it, the nearest point exists and is unique, so proj is the paper's $P$. Theorem 5.1 is stated relationally (nearest-point predicate IsProj) to avoid junk values altogether.
• "Contraction" is Mathlib's ContractingWith c F with c : ℝ≥0: c < 1 and a Lipschitz bound with that single constant. A constant allowed to depend on the points, or a Lipschitz bound without c < 1, is not a contraction and would trivialize the goal; so would a projection whose junk value is reachable (e.g. with $\tilde K$ empty), which makes $F_\lambda$ unrelated to the paper's map. The fixed point must be linked to the GQVI and to the iteration from every starting point.
• The square root in the Lipschitz estimate is Real.sqrt; its radicand is nonnegative under the hypotheses when $n\ge1$.
• Useful infrastructure: Mathlib's ContractingWith.fixedPoint and ContractingWith.tendsto_iterate_fixedPoint (Banach), exists_norm_eq_iInf_of_complete_convex and norm_eq_iInf_iff_real_inner_le_zero (projection on convex sets), and the platform theorem VectorSpaceOpt.min_distance_convex_set. A general lemma that the Euclidean nearest-point map of a closed convex set is 1-Lipschitz is reusable well beyond this mission and is welcome as a separate contribution.

Selected references

• D. Chan and J. S. Pang, The generalized quasi-variational inequality problem, Mathematics of Operations Research 7(2) (1982) 211–222. https://doi.org/10.1287/moor.7.2.211
• S. C. Fang and E. L. Peterson, Generalized variational inequalities, Journal of Optimization Theory and Applications 38 (1982) 363–383. https://doi.org/10.1007/BF00935344
• Y. Nesterov and L. Scrimali, Solving strongly monotone variational and quasi-variational inequalities, Discrete and Continuous Dynamical Systems 31(4) (2011) 1383–1396. https://doi.org/10.3934/dcds.2011.31.1383

Motivation

Quadratic programming, the minimisation of a quadratic function subject to linear inequality constraints, is a basic subproblem of nonlinear optimisation (sequential quadratic programming, trust-region methods) and a model in its own right in portfolio selection, least-squares estimation and control. The classical solution methods are active-set methods: they keep a set of constraints treated as equalities, minimise over the resulting affine set, and then decide which constraint to drop or add (Gill, Murray and Wright, Practical Optimization, 1981; Fletcher, Practical Methods of Optimization, Vol. 2, 1981). Their finite-termination proofs need either a nondegeneracy assumption (linearly independent active constraints) or an anti-cycling rule, because under degeneracy the choice of the constraint to drop, made from Lagrange multiplier estimates, can cycle.

Calamai and Moré (Mathematical Programming 39, 1987) showed that the gradient projection method can take over the step that leaves a working set. Their Algorithm 6.1 alternates two kinds of step: an arbitrary non-increasing step that adds constraints to the working set until the equality-constrained subproblem is solved, and a single projected-gradient step once it is solved. Theorem 6.2 states that this algorithm terminates at a stationary point for every quadratic that is bounded below on the feasible polyhedron, with no nondegeneracy assumption and no anti-cycling rule. The same paper's Sections 2–4 (the subject of the first two missions of this series) supply the properties of the gradient projection step that the argument uses.

Timeline of the ingredients:

• 1964, 1966: Goldstein and Levitin–Polyak introduce the gradient projection method $x_{k+1} = P(x_k - \alpha_k\nabla f(x_k))$ for convex constraint sets.
• 1976: Bertsekas proves finite identification of the active constraints for bound constraints and the Armijo rule.
• 1981: Dunn uses the descent inequalities (2.4)–(2.5) in the analysis of the method.
• 1987: Calamai and Moré generalise the step rule to (2.1)–(2.2), prove convergence of projected gradients, identification of active constraints for general polyhedra, and finite termination of Algorithm 6.1.

Setting

Let $E$ be a finite-dimensional real inner product space (the paper's $\mathbb{R}^n$ with a general inner product). The feasible set is a polyhedron

$$\Omega = \{x \in E : \langle c_j, x\rangle \ge \delta_j,\ j = 1, \dots, m\},$$

with active set $A(x) = \{j : \langle c_j, x\rangle = \delta_j\}$. The objective is a quadratic function $f(x) = \tfrac12\langle x, Qx\rangle + \langle b, x\rangle + c_0$ with $Q$ self-adjoint but not necessarily positive semidefinite, so $f$ may be nonconvex. Its gradient $\nabla f$ is taken with respect to the inner product of $E$.

The projection into $\Omega$ is $P(x) = \operatorname{argmin}\{\|z - x\| : z \in \Omega\}$, and a point $x^* \in \Omega$ is stationary if $\langle\nabla f(x^*), x - x^*\rangle \ge 0$ for every $x \in \Omega$.

A gradient projection step from $x_k$ is $x_{k+1} = P(x_k - \alpha_k\nabla f(x_k))$ with $\alpha_k > 0$ satisfying the sufficient decrease condition (2.1) with a constant $\mu_1 \in (0,1)$, the condition (2.2) that $\alpha_k \ge \gamma_1$ or $\alpha_k \ge \gamma_2\bar\alpha_k > 0$ for some $\bar\alpha_k$ at which the decrease test (2.3) with constant $\mu_2 \in (0,1)$ fails, and the upper bound $\alpha_k \le \gamma_3$ (3.2).

A working set is a set $W \subseteq \{1, \dots, m\}$; problem (6.2) is $\min\{f(y) : \langle c_j, y\rangle = \delta_j,\ j \in W\}$, over an affine set that ignores the inequality constraints outside $W$.

Algorithm 6.1 produces iterates $x_k \in \Omega$ and working sets $W_k \subseteq A(x_k)$ from $x_0 \in \Omega$:

• (a) if $x_k$ is a global minimiser of (6.2) for $W_k$, then $x_{k+1}$ is a gradient projection step from $x_k$;
• (b) otherwise $x_{k+1} \in \Omega$, $f(x_{k+1}) \le f(x_k)$, $W_k \subseteq W_{k+1}$, and if $W_{k+1} = W_k$ then $x_{k+1}$ is a global minimiser of (6.2).

Formalization targets

Goal: Theorem 6.2

For every quadratic $f$ bounded below on $\Omega$, all constants $\gamma_1, \gamma_2 > 0$, $\mu_1, \mu_2 \in (0,1)$, $\gamma_3 \in \mathbb{R}$, and every run $(x_k, W_k, \alpha_k)_{k\ge0}$ of Algorithm 6.1,

$$\exists\, l \ge 0 :\quad \langle \nabla f(x_l), x - x_l\rangle \ge 0 \quad \text{for all } x \in \Omega.$$

The theorem makes no assumption on the boundedness of the iterates and none on the linear independence of the constraints.

Milestones

• Lemma 2.1(a): for nonempty closed convex $\Omega$, $z \in \Omega$ and any $x$, $\langle P(x) - x, z - P(x)\rangle \ge 0$.
• Eq. (2.5): for $x_k \in \Omega$, $\alpha_k > 0$ and $x_{k+1} = P(x_k - \alpha_k\nabla f(x_k))$,

$$\langle\nabla f(x_k), x_k - x_{k+1}\rangle \ge \frac{\|x_{k+1} - x_k\|^2}{\alpha_k}.$$

Significance

Theorem 6.2 separates the two roles an active-set method plays: solving equality-constrained subproblems, for which any method that does not increase $f$ may be used, and choosing the next working set, which the gradient projection step does. The consequence is a finitely terminating quadratic programming algorithm for nonconvex quadratics that needs neither nondegeneracy nor an anti-cycling rule, and a template for large-scale bound-constrained and linearly constrained solvers that combine projection steps with subspace minimisation.

The result is proved in the paper and is classical; no machine-checked version is known. A formalization produces a checked finite-termination theorem for an active-set method on degenerate problems, together with reusable pieces: the projection onto a polyhedron as a total function with its variational inequality, the descent estimate of a projected step, and a predicate describing active-set runs with working sets, which other active-set algorithms can reuse.

Difficulty

The obvious argument — each step decreases $f$ and there are finitely many working sets — fails on two counts. First, step (b) only guarantees $f(x_{k+1}) \le f(x_k)$, so $f$ values alone do not rule out infinitely many iterations; the nesting of working sets and the final clause of step (b) are what bound consecutive (b)-steps. Second, a gradient projection step taken at a solution of (6.2) can in principle leave $f$ unchanged, and nothing in the algorithm's rules says directly that it makes progress; that it does so at every non-stationary iterate is a property of the projection and of the step conditions (2.1)–(2.2), not of the algorithm. A further point is that the minimum of (6.2) is taken over an affine set, not over $\Omega$, and the link between the value at a step-(a) iterate and later iterates runs through the requirement $W_k \subseteq A(x_k)$.

Formalization scope

The space is a finite-dimensional real inner product space E; $\nabla f$ is Mathlib's gradient. Constraints are indexed by Fin m, working and active sets are Finset (Fin m), and a run is a predicate IsAlgorithm61Run on three sequences $x : \mathbb{N} \to E$, $W$, $\alpha$ indexed from $0$. The run does not stop by itself; "the algorithm terminates at a stationary iterate" is rendered as the existence of an index $l$ with $x_l$ stationary. The projection is argmin made total by a junk value $0$ that is unreachable when $\Omega$ is nonempty, closed and convex. The paper's "$\subset$" between working sets is inclusion. A quadratic function is $\tfrac12\langle x, Qx\rangle + \langle b, x\rangle + c_0$ with $Q$ symmetric and no definiteness assumption.

A trivializing formalization — a run predicate that forces $x_0$ to be stationary, or that no sequence satisfies — would make the goal empty; the step rules here are the paper's verbatim, and a run on $\Omega = [0,\infty) \subset \mathbb{R}$ with $f(x) = x$ starting at the non-stationary point $x_0 = 1$ satisfies the predicate. Replacing step (a) by "any step that strictly decreases $f$" would assume the central fact and is not acceptable.

A complete development needs: the variational inequality of the projection (Lemma 2.1(a)), the descent estimate (2.5), the characterisation of stationary points as fixed points of the projected step, and a finiteness argument over the finitely many subsets of Fin m. Proofs of the milestones, of these auxiliary facts, and of the goal are all welcome.

Selected references

• P. H. Calamai and J. J. Moré, Projected gradient methods for linearly constrained problems, Mathematical Programming 39 (1987) 93–116. https://doi.org/10.1007/BF02592073
• A. A. Goldstein, Convex programming in Hilbert space, Bulletin of the AMS 70 (1964) 709–710. https://doi.org/10.1090/S0002-9904-1964-11178-2
• E. S. Levitin and B. T. Polyak, Constrained minimization methods, USSR Computational Mathematics and Mathematical Physics 6 (1966) 1–50. https://doi.org/10.1016/0041-5553(66)90114-5
• D. P. Bertsekas, On the Goldstein–Levitin–Polyak gradient projection method, IEEE Transactions on Automatic Control 21 (1976) 174–184. https://doi.org/10.1109/TAC.1976.1101194
• J. C. Dunn, Global and asymptotic convergence rate estimates for a class of projected gradient processes, SIAM Journal on Control and Optimization 19 (1981) 368–400. https://doi.org/10.1137/0319022
• P. E. Gill, W. Murray and M. H. Wright, Practical Optimization, Academic Press, 1981. https://doi.org/10.1137/1.9781611975604

Motivation

Many combinatorial optimization problems on graphs are linear programs over a polytope whose inequality description is unknown. The stable set polytope is the standard example: maximizing a linear function over it is the maximum weight stable set problem, which is NP-hard, and no complete inequality description is known for general graphs. A productive line of work, begun in V. Chvátal's 1975 paper On certain polytopes associated with graphs (J. Combin. Theory Ser. B 18 (1975) 138–154), asks instead how such descriptions behave under graph operations: if descriptions are known for small graphs, can one write one down for a graph built from them?

Section 5 of that paper answers this for substitution, the operation that replaces a vertex of one graph by a whole second graph. Substitution contains three familiar constructions as special cases: duplicating a vertex, forming the join of two graphs, and forming the lexicographic product (composition). Duplication is one of the two ingredients of Lovász's proof of the perfect graph theorem (Lovász 1972); substitution in general is the operation under which perfection is preserved, and graphs built from simple pieces by substitution are a recurring source of classes with tractable stable set polytopes.

Setting

All graphs are finite, undirected and loopless. A stable set of a graph $G=(V,E)$ is a set of vertices no two of which are adjacent. Write $S(G)\subseteq\mathbb R^V$ for the set of incidence vectors of stable sets (the zero–one vectors $x$ with $\{u:x_u=1\}$ stable), and

$$P(G)=\operatorname{conv}S(G)$$

for the stable set polytope. A finite system of linear inequalities in the variables $(x_u:u\in V)$ is a defining linear system of $P(G)$ when its set of solutions is exactly $P(G)$.

Let $G_1=(V_1,E_1)$ and $G_2=(V_2,E_2)$ be graphs with $V_1\cap V_2=\emptyset$, and let $v\in V_1$. The graph $G$ obtained from $G_1$ by substituting $G_2$ for $v$ has vertex set $(V_1-\{v\})\cup V_2$. Its edges are the edges of $G_1-v$, the edges of $G_2$, and every edge joining a vertex of $G_2$ to a neighbour of $v$ in $G_1$. In Lean the vertex type is the disjoint sum {u : V₁ // u ≠ v} ⊕ V₂ and the graph is substitute G₁ v G₂.

Formalization targets

Goal: Theorem 5.1

For $k\in\{1,2\}$ let

$$-x_u\le 0\ (u\in V_k),\qquad \sum_{u\in V_k}a_{iu}x_u\le b_i\ (i\in J_k)$$

be a defining linear system of $P(G_k)$, with $J_1,J_2$ finite index sets and real coefficients, and put $a^+_{iv}=\max\{a_{iv},0\}$ for $i\in J_1$. Then

$$-x_u\le 0\ \ (u\in V_2\cup(V_1-\{v\})),\qquad a^+_{iv}\sum_{u\in V_2}a_{ju}x_u+b_j\sum_{u\in V_1-\{v\}}a_{iu}x_u\le b_ib_j\ \ (i\in J_1,\ j\in J_2)\tag{5.1}$$

is a defining linear system of $P(G)$. The statement fixes no particular system for $G_1$ or $G_2$: any defining systems of the two pieces produce one for $G$, with $|J_1|\cdot|J_2|$ rows besides nonnegativity.

Milestones

1. Validity of (5.1) (§5, p. 145): every $x\in S(G)$ satisfies (5.1), hence so does every point of $P(G)$.
2. Proposition 2.1 (pp. 139–140): for a finite nonempty set $S$ of solutions of a system with nonnegativity rows $-x_u\le 0$, the solution set equals $\operatorname{conv}S$ if and only if for every integer vector $c$ the value $\max\{cx:x\in S\}$ equals the minimum of the associated dual linear program, the minimum being attained.
3. Decomposition of the optimum (§5, pp. 145–146): for an integer vector $c$ on $V_2\cup W$, $W=V_1-\{v\}$, with $d_u=\max\{c_u,0\}$,

$$\max\{cx:x\in S(G)\}=\max\{m_0,\ m_1+m_2\},$$

where $m_0$ and $m_1$ are the maxima of $\sum_{u\in W}d_ux_u$ over $x\in S(G_1)$ with $x_v=0$ and $x_v=1$ respectively, and $m_2$ is the maximum of $\sum_{u\in V_2}d_ux_u$ over $S(G_2)$.

Significance

Theorem 5.1 gives an explicit construction: from a polyhedral description of $P(G_1)$ and $P(G_2)$ it writes one of $P(G)$, row by row, with no loss. Specialized to $G_1=K_2$ it gives Corollary 5.2 of the paper, a defining linear system for the join $G_1+G_2$; applied repeatedly it gives defining systems for lexicographic products, and applied with $G_2=\overline{K_2}$ it describes the effect of duplicating a vertex. Applied to clique systems, whose coefficients are 0 and 1, the rows of (5.1) are again clique inequalities of $G$, so the class of graphs whose stable set polytope is described by nonnegativity and clique inequalities is closed under substitution.

The result has been proved in print since 1975. As far as a search of the Prove2Me catalogue shows, none of it, including Proposition 2.1 and the substitution operation itself, has a machine-checked statement or proof. The mission asks for a formal proof of the theorem and of the two combinatorial and polyhedral steps it rests on. The definitions of $S(G)$, $P(G)$ and graph substitution, and the LP characterization of Proposition 2.1, are reusable by every other mission on stable set polytopes and on polyhedral descriptions of 0–1 sets.

Difficulty

That every point of $P(G)$ satisfies (5.1) is a short case check on stable sets of $G$. The difficulty is the reverse inclusion: that no point outside $P(G)$ satisfies (5.1). The first idea, taking a point that satisfies (5.1) and splitting it directly into a point of $P(G_1)$ and a point of $P(G_2)$, fails: (5.1) couples the two input systems through products of their coefficients and right-hand sides, and a fractional solution of (5.1) carries no evident decomposition into the two pieces. Nothing is assumed about the signs of the input coefficients, so the rows of (5.1) can mix positive and negative terms, and the positive part $a^+_{iv}$ in place of $a_{iv}$ is what keeps the system valid when $a_{iv}<0$.

The polyhedral step behind Proposition 2.1, relating a convex hull of finitely many points to an inequality system through linear programming duality, is not available in Mathlib in this form and has to be built.

Formalization scope

• Graphs are SimpleGraph on a Fintype with decidable equality; the substituted graph lives on {u : V₁ // u ≠ v} ⊕ V₂, which builds in $V_1\cap V_2=\emptyset$.
• $S(G)$ is the set of real incidence vectors of finite stable sets (IsIndepSet); $P(G)$ is convexHull ℝ (S G), never the solution set of an inequality system.
• A linear system is a finite index type $J$ with real a : J → V → ℝ, b : J → ℝ. The nonnegativity rows $-x_u\le0$ are kept as a separate conjunct ∀ u, 0 ≤ x u everywhere; Proposition 2.1 is false without them. "Defining linear system" is set equality of the solution set with $P(G)$.
• No sign conditions on the $a_{iu}$ or $b_i$ are assumed; the paper assumes none.
• Implicit hypothesis made explicit: $V_2\ne\emptyset$ ([Nonempty V₂]) in Theorem 5.1. The paper's graphs have nonempty vertex sets and its proof picks a vertex of $G_2$; with $V_2=\emptyset$, $J_2=\emptyset$ and $V_1\ne\{v\}$, (5.1) is just $x\ge0$ and the theorem fails. The validity milestone does not need it.
• In Proposition 2.1 the set $S$ is assumed nonempty, which the paper's $\max\{cx:x\in S\}$ presupposes. "max = min" is stated as a lower bound for every feasible dual vector plus a feasible dual vector attaining the maximum.
• In the decomposition milestone each maximum is a real sSup over a finite set that always contains the zero vector or the incidence vector of $\{v\}$, so no junk value of sSup can occur.
• A trivializing formalization is excluded: $P(G)$ is the convex hull of stable-set vectors rather than a set defined through the same inequalities, and the goal is the full set equality, not the validity inclusion alone.

Contributions welcome: a proof of Proposition 2.1 (the reusable core), the combinatorial decomposition, the validity case check, and the assembly of the goal.

Selected references

• V. Chvátal, On certain polytopes associated with graphs, J. Combin. Theory Ser. B 18 (1975) 138–154. https://doi.org/10.1016/0095-8956(75)90041-6
• L. Lovász, Normal hypergraphs and the perfect graph conjecture, Discrete Math. 2 (1972) 253–267. https://doi.org/10.1016/0012-365X(72)90006-4
• J. Edmonds, Maximum matching and a polyhedron with 0,1-vertices, J. Res. Nat. Bur. Standards 69B (1965) 125–130. https://doi.org/10.6028/jres.069B.013
• F. Harary, Graph Theory, Addison-Wesley, 1969.

Motivation

Minimizing a smooth function subject to linear inequality constraints is the core subproblem of much of nonlinear optimization: bound-constrained problems, quadratic programs, and the subproblems of sequential quadratic programming and augmented Lagrangian methods all have this form. Methods for these problems are usually built from two parts, one that decides which constraints hold with equality at the solution and one that solves the resulting equality-constrained problem quickly. The first part only pays off if the decision stabilizes after finitely many iterations; otherwise the fast local method never gets to run.

Calamai and Moré (Math. Programming 39, 1987) proved that this stabilization is a property of the limit point, not of the algorithm. Any feasible sequence that converges and whose projected gradients tend to zero identifies the active constraints of a nondegenerate limit in finitely many steps. This is the result that later active-set and gradient-projection methods for bound-constrained and linearly constrained problems invoke to justify switching to a fast local phase.

Timeline.

• 1976: Bertsekas proves finite identification of the active set for the gradient projection method with an Armijo step on bound constraints, at a local minimizer satisfying strict complementarity and second-order sufficiency.
• 1984: Gafni and Bertsekas (SIAM J. Control Optim. 22) prove a similar result for two-metric projection methods, under an assumption that excludes the choice of the gradient as search direction.
• 1987: Calamai and Moré remove the second-order condition, allow a general polyhedral feasible set and a general inner product, and make the result independent of the method generating the sequence (Theorem 4.1); they extend it to binding sets defined by multiplier estimates (Theorem 4.2).

Setting

Let $E$ be a finite-dimensional real inner product space (the paper's $\mathbb{R}^n$ with a general inner product) and let $f : E \to \mathbb{R}$ be continuously differentiable on the feasible set, with gradient $\nabla f$ taken with respect to the inner product of $E$.

The feasible set is a polyhedral set

$$\Omega = \{x \in E : \langle c_j, x\rangle \ge \delta_j,\ j = 1, \dots, m\}$$

for constraint normals $c_j \in E$ and scalars $\delta_j$. The active set at $x$ is $A(x) = \{j : \langle c_j, x\rangle = \delta_j\}$.

A direction $v$ is feasible at $x \in \Omega$ if $x + \tau v \in \Omega$ for all sufficiently small $\tau > 0$. The tangent cone $T(x)$ is the closure of the set of feasible directions. The projected gradient is the point of $T(x)$ closest to $-\nabla f(x)$:

$$\nabla_\Omega f(x) = \operatorname{argmin}\{\|v + \nabla f(x)\| : v \in T(x)\}.$$

A point $x^* \in \Omega$ is stationary if $\langle \nabla f(x^*), x - x^*\rangle \ge 0$ for all $x \in \Omega$. It is a Kuhn–Tucker point if $\nabla f(x^*) = \sum_{j \in A(x^*)} \lambda^*_j c_j$ with $\lambda^*_j \ge 0$. It is nondegenerate if the active normals $\{c_j : j \in A(x^*)\}$ are linearly independent and the multipliers satisfy $\lambda^*_j > 0$ for every $j \in A(x^*)$.

A Lagrange multiplier estimate is a map $x \mapsto \lambda(x) \in \mathbb{R}^m$. It defines the binding set $B(x) = \{j \in A(x) : \lambda_j(x) \ge 0\}$. The estimate is consistent if $\lambda_j(x_k) \to \lambda_j(x^*)$ whenever $x_k \to x^*$, the point $x^*$ is a nondegenerate Kuhn–Tucker point, and $A(x_k) = A(x^*)$ for every $k$.

Formalization targets

Goal: Theorem 4.1 (finite identification of the active set)

Let $\{x_k\}$ be an arbitrary sequence in $\Omega$ converging to $x^*$. If $\|\nabla_\Omega f(x_k)\| \to 0$ and $x^*$ is nondegenerate, then

$$A(x_k) = A(x^*) \quad \text{for all sufficiently large } k.$$

The sequence need not come from any particular algorithm. The goal asserts eventual equality of the index sets, not inclusion.

Milestones

• Lemma 3.1. At $x \in \Omega$: $-\langle\nabla f(x), \nabla_\Omega f(x)\rangle = \|\nabla_\Omega f(x)\|^2$; $\min\{\langle \nabla f(x), v\rangle : v \in T(x), \|v\| \le 1\} = -\|\nabla_\Omega f(x)\|$; and $x$ is stationary if and only if $\nabla_\Omega f(x) = 0$.
• Lemma 3.3. The map $x \mapsto \|\nabla_\Omega f(x)\|$ is lower semicontinuous on $\Omega$.
• Tangent cone of a polyhedron (p. 105). For $x \in \Omega$, $T(x) = \{v : \langle c_j, v\rangle \ge 0,\ j \in A(x)\}$.
• Eq. (4.3). For polyhedral $\Omega$, a point $x^* \in \Omega$ is stationary if and only if it is a Kuhn–Tucker point.
• Theorem 4.2. Assume the binding sets come from a consistent estimate whose value at $x^*$ is the Kuhn–Tucker multiplier vector, and assume the hypotheses of Theorem 4.1. Then $B(x_k) = B(x^*)$ for all sufficiently large $k$.

Significance

The result. Theorem 4.1 separates identification from convergence. Any method that keeps its iterates feasible and drives the projected gradient to zero inherits finite identification, whatever its step-size rule or search direction. After identification the constrained problem is locally an unconstrained problem on the affine subspace $\{x : \langle c_j, x\rangle = \delta_j,\ j \in A(x^*)\}$, so Newton-type or conjugate-gradient methods can take over. Theorem 4.2 carries the same conclusion to methods that drop constraints according to the signs of multiplier estimates. The companion missions of this series use the result: the gradient projection method drives the projected gradients to zero (mission I), and a gradient projection algorithm for quadratic programs terminates finitely (mission III).

Formalizing it. The theorems are proved in the paper. No machine-checked version of the projected gradient, of tangent cones of polyhedra with their active-set description, or of finite active-set identification is known to exist. Formalization adds a reusable account of tangent cones and polar cones of polyhedral sets and of the equivalence between stationarity and the Kuhn–Tucker conditions for linear constraints, together with a method-independent identification theorem stated at the level of generality of the paper.

Difficulty

Two different limits are involved. Convergence $x_k \to x^*$ is enough to show that no inactive constraint of $x^*$ is active at $x_k$ for large $k$. The hard direction is the converse: a constraint active at $x^*$ might be inactive at infinitely many $x_k$, approached from the interior. Convergence of the points alone cannot rule this out. The projected gradient is also not continuous, because the tangent cone changes when a new constraint becomes active. So the hypothesis $\|\nabla_\Omega f(x_k)\| \to 0$ cannot be passed to the limit naively. Both nondegeneracy conditions matter: without linear independence, or with a zero multiplier, the statement fails.

Formalization scope

The space is a real inner product space E with [FiniteDimensional ℝ E], and $\nabla f$ is Mathlib's gradient. "Continuously differentiable on $\Omega$" means DifferentiableAt ℝ f x for every $x \in \Omega$ together with ContinuousOn (gradient f) Ω. The constraints are indexed by Fin m. $\Omega$ is polyhedron c δ, and $A(x)$ is activeSet c δ x : Finset (Fin m).

The tangent cone is defined as the closure of the feasible directions, not by the polyhedral formula, which is a milestone. The projected gradient is the nearest point of $T(x)$ to $-\nabla f(x)$, chosen by a choice function that returns $0$ only when no nearest point exists. That never happens at a point of a polyhedral set.

Nondegeneracy is bundled as IsNondegenerate c δ f x*: $x^* \in \Omega$, the family $(c_j)_{j \in A(x^*)}$ is linearly independent, and positive multipliers represent $\nabla f(x^*)$. "For all sufficiently large $k$" is ∀ᶠ k in Filter.atTop.

In Theorem 4.2 the paper leaves one condition implicit: the estimate at $x^*$ must be the Kuhn–Tucker multiplier vector, $\nabla f(x^*) = \sum_{j \in A(x^*)} \lambda_j(x^*) c_j$. Without it the statement is false, so it is an explicit hypothesis. Consistency is required only along feasible sequences and only in the coordinates $j \in A(x^*)$.

A formalization that assumes $A(x_k) \subseteq A(x^*)$, assumes the active sets are eventually constant, weakens nondegeneracy to nonnegative multipliers, or concludes only inclusion is not the paper's theorem and does not satisfy this mission.

Needed infrastructure: tangent cones of convex sets, the Moreau decomposition into a closed convex cone and its polar, Farkas' lemma in a general inner product space, and orthogonal projections onto subspaces spanned by linearly independent vectors. The polyhedral tangent-cone and Kuhn–Tucker results are reusable beyond this mission. Proofs of any milestone are welcome, as are auxiliary lemmas on polyhedral cones.

Selected references

• P. H. Calamai and J. J. Moré, Projected gradient methods for linearly constrained problems, Mathematical Programming 39 (1987) 93–116. https://doi.org/10.1007/BF02592073
• D. P. Bertsekas, On the Goldstein–Levitin–Polyak gradient projection method, IEEE Transactions on Automatic Control 21 (1976) 174–184. https://doi.org/10.1109/TAC.1976.1101194
• E. M. Gafni and D. P. Bertsekas, Two-metric projection methods for constrained optimization, SIAM Journal on Control and Optimization 22 (1984) 936–964. https://doi.org/10.1137/0322061
• E. H. Zarantonello, Projections on convex sets in Hilbert space and spectral theory, in: Contributions to Nonlinear Functional Analysis, Academic Press, 1971, 237–424.

Motivation

Minimizing a smooth function over a closed convex set $\Omega\subseteq\mathbb R^n$ on which projection is cheap (a box, a ball, a simplex) is a routine subproblem in large-scale optimization. Box-constrained minimization is the inner solver of augmented Lagrangian methods, and bound-constrained least squares, image restoration and density estimation all have this form. The classical gradient projection method of Goldstein and of Levitin and Polyak needs only gradients and projections, but with constant or monotone Armijo step lengths it is slow.

Spectral projected gradient (SPG) methods, introduced by Birgin, Martínez and Raydan (paper), combine three ingredients. The first is the projection. The second is the Barzilai–Borwein (spectral) step length $\alpha_{k+1}=\langle s_k,s_k\rangle/\langle s_k,y_k\rangle$, an inverse Rayleigh quotient of the average Hessian along the last step. The third is the nonmonotone line search of Grippo, Lampariello and Lucidi, which compares a trial value with the worst of the last $M$ objective values instead of the current one. The paper defines two variants. This mission concerns SPG1, which backtracks along the projection arc $\lambda\mapsto P(x_k-\lambda g(x_k))$, as in Bertsekas's analysis of the Armijo rule for gradient projection. The companion mission concerns SPG2, which backtracks along a fixed feasible direction.

Timeline:

• 1964–1966: Goldstein; Levitin and Polyak introduce gradient projection.
• 1976: Bertsekas analyses the Armijo rule along the projection arc (IEEE TAC).
• 1986: Grippo, Lampariello and Lucidi introduce the nonmonotone line search for unconstrained problems.
• 1988: Barzilai and Borwein propose the two-point step size. Raydan (1997) combines it with nonmonotone search in the unconstrained case.
• 2000: Birgin, Martínez and Raydan define SPG1 and SPG2 for convex constraints (SIAM J. Optim. 10(4)).
• 2003: the same authors publish the convergence analysis that the proof of Theorem 2.2 adapts, in the inexact setting (IMA J. Numer. Anal. 23).

Setting

Let $\Omega\subseteq\mathbb R^n$ be nonempty, closed and convex, with the Euclidean inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$. Let $f$ have continuous partial derivatives on an open set $U\supseteq\Omega$, and write $g(x)=\nabla f(x)$. The orthogonal projection $P(z)$ is the unique point of $\Omega$ nearest to $z$. The scaled projected gradient is $g_t(x)=P(x-t\,g(x))-x$ for $x\in\Omega$ and $t>0$. A point $\bar x$ is a constrained stationary point if $\langle g(\bar x),x-\bar x\rangle\ge0$ for all $x\in\Omega$.

The parameters are an integer $M\ge1$, reals $0<\alpha_{\min}<\alpha_{\max}$, a sufficient-decrease constant $\gamma\in(0,1)$ and safeguards $0<\sigma_1<\sigma_2<1$. Algorithm SPG1 (Algorithm 2.1) starts from $x_0\in\Omega$ and $\alpha_0\in[\alpha_{\min},\alpha_{\max}]$. At iteration $k=0,1,\dots$ it does the following.

1. Stop test. If $\|P(x_k-g(x_k))-x_k\|=0$, stop: $x_k$ is stationary.
2. Backtracking along the projection arc. Set $\lambda=\alpha_k$. While the trial point $x_+=P(x_k-\lambda g(x_k))$ fails

$$f(x_+)\le\max_{0\le j\le\min\{k,M-1\}}f(x_{k-j})+\gamma\langle x_+-x_k,g(x_k)\rangle,\qquad(1)$$

replace $\lambda$ by any $\lambda_{\rm new}\in[\sigma_1\lambda,\sigma_2\lambda]$. When (1) holds, set $\lambda_k=\lambda$ and $x_{k+1}=x_+$. 3. Spectral step. With $s_k=x_{k+1}-x_k$, $y_k=g(x_{k+1})-g(x_k)$ and $b_k=\langle s_k,y_k\rangle$, set $\alpha_{k+1}=\alpha_{\max}$ if $b_k\le0$, and otherwise $\alpha_{k+1}=\min\{\alpha_{\max},\max\{\alpha_{\min},\langle s_k,s_k\rangle/b_k\}\}$.

The first trial of each backtracking is the spectral step $\alpha_k$, not $1$. The sufficient-decrease term in (1) is $\gamma\langle x_+-x_k,g(x_k)\rangle=\gamma\langle g(x_k),g_\lambda(x_k)\rangle$, with no factor $\lambda$.

In Lean these objects are written as follows:

• the projection is a function P with the predicate IsProjOnto Ω P;
• $g_t$ is scaledProjGrad P f t;
• stationarity is IsConstrainedStationary Ω f;
• the maximum in (1) is nonmonotoneRef f x M k;
• test (1) is SPG1Test;
• an infinite run is IsSPG1Run Ω f P M αmin αmax γ σ₁ σ₂ x α.

Formalization targets

Goal: Theorem 2.2, accumulation points are stationary

For every infinite run $(x_k,\alpha_k)$ of SPG1 and every accumulation point $\bar x$ of $(x_k)$,

$$\langle g(\bar x),x-\bar x\rangle\ge0\qquad\text{for all }x\in\Omega.$$

The statement fixes no parameter values, and it assumes neither convexity of $f$ nor a bounded level set.

Milestones

• Lemma 2.1 (ii). For $\bar x\in\Omega$ and $t\in(0,\alpha_{\max}]$, $g_t(\bar x)=0$ if and only if $\bar x$ is a constrained stationary point.
• Lemma 2.1 (i). For $x\in\Omega$ and $t\in(0,\alpha_{\max}]$,

$$\langle g(x),g_t(x)\rangle\le-\tfrac1t\|g_t(x)\|_2^2\le-\tfrac1{\alpha_{\max}}\|g_t(x)\|_2^2.$$

• Lemma 2.2 (i). For $x\in\Omega$ and $z\in\mathbb R^n$, the map $s\mapsto\|P(x+sz)-x\|/s$ is nonincreasing on $s>0$.
• Lemma 2.2 (ii). For every $x\in\Omega$ there is $s_x>0$ such that $f(P(x-tg(x)))-f(x)\le\gamma\langle g(x),g_t(x)\rangle$ for all $t\in[0,s_x]$.
• Theorem 2.2, first clause (SPG1 is well defined). At a point where Step 1 does not stop, every admissible backtracking sequence starting at $\alpha\in[\alpha_{\min},\alpha_{\max}]$ reaches a trial point satisfying (1). The statement is for an arbitrary reference value $R\ge f(x)$, which covers the maximum in (1).

Significance

Theorem 2.2 is the global convergence guarantee for SPG1. It holds without monotone decrease of $f$ and with no restriction on the spectral step beyond the safeguards. Lemma 2.2 carries Bertsekas's curvilinear Armijo analysis, stated for monotone gradient projection, over to the nonmonotone spectral setting. The projection-arc search is the natural one when $\Omega$ is a box or a polyhedron: there the arc is piecewise linear and each trial point is feasible by construction.

Status: the theorem is proved in the literature. This paper's proof reads "Use Lemma 2.2 with the proof technique of [7]", and Lemma 2.2 is quoted from Bertsekas's Nonlinear Programming (Lemma 2.3.1 and Theorem 2.3.3 (a)). No Lean formalization of this theorem, of the Armijo analysis along the projection arc, or of the monotonicity of $\|P(x+sz)-x\|/s$ is known. The mission produces a formal proof and a reusable Lean interface for projection-based first-order methods on convex sets.

Difficulty

For monotone descent methods, the usual argument shows that $f(x_k)$ decreases, so the total decrease is finite and the per-iteration decrease tends to zero. That argument fails here, because $f(x_k)$ need not decrease. Only the window maximum $\max_{0\le j\le\min\{k,M-1\}}f(x_{k-j})$ is nonincreasing, and a small decrease of this maximum does not by itself give a small decrease at the iterates that approach a given accumulation point $\bar x$.

Along the projection arc there is a second obstacle. The decrease predicted by (1) is $\gamma\langle g(x_k),g_{\lambda_k}(x_k)\rangle$, and $g_\lambda(x_k)$ depends nonlinearly on $\lambda$: for $\lambda<\alpha_k$ the trial point is not a rescaling of the first one. So small accepted steps do not translate into small multiples of a fixed direction, as they do for SPG2. The step lengths $\lambda_k$ may also tend to zero, $f$ is $C^1$ only on a neighbourhood of $\Omega$, and no Lipschitz constant for $g$ is available.

Formalization scope

• Space and data. The space is EuclideanSpace ℝ (Fin n) with inner ℝ and the 2-norm. $f$ is a total function EuclideanSpace ℝ (Fin n) → ℝ with ContDiffOn ℝ 1 f U on an open U ⊇ Ω, and $g$ is Mathlib's gradient f. Every trial point is a projection, so the algorithm evaluates $f$ and $g$ only at points of $\Omega$.

• Iteration and trials. Iterations are indexed from $0$. The backtracking choice (2) is universally quantified. At each iteration, a run carries a finite trial list with $\lambda^{(0)}=\alpha_k$ and $\lambda^{(i+1)}\in[\sigma_1\lambda^{(i)},\sigma_2\lambda^{(i)}]$; test (1) fails at every trial but the last and holds at the last.

• Step size. $\alpha_{k+1}$ is given by Step 3 exactly.

• Accumulation point. An accumulation point is MapClusterPt x̄ atTop x.

• Lemma 2.2 (i). The paper names the domain $[0,\infty)$ but defines $h$ only for $s>0$, so the milestone is stated on $(0,\infty)$.

• Excluded simplifications. None of the following is SPG1:

  • a run predicate that accepts any positive step;
  • a run predicate that starts backtracking at $1$;
  • a run predicate that uses SPG2's test $\gamma\lambda\langle d_k,g(x_k)\rangle$;
  • a run predicate that lets $\alpha_{k+1}$ range freely over $[\alpha_{\min},\alpha_{\max}]$.

  Nor is a goal that states $g_t(\bar x)=0$ instead of the variational inequality, or one that adds convexity of $f$, a Lipschitz gradient or a bounded level set.

• Non-vacuity. The hypotheses of the goal are satisfiable. Take $f(x)=\|x\|^2$, $\Omega=\mathbb R^n$, $M=1$, $\alpha_{\min}=1/8$, $\alpha_{\max}=1/4$, $\gamma=1/2$, $\sigma_1=1/10$, $\sigma_2=9/10$ and $v\ne0$. Then the iterates $x_k=2^{-k}v$ with $\alpha_k=1/4$ form an infinite run with accumulation point $0$.

• Infrastructure. A complete development needs:

  • the variational characterization of the projection (Mathlib has it in the iInf form, norm_eq_iInf_iff_real_inner_le_zero) and the nonexpansiveness of the projection;
  • a first-order expansion of a $C^1$ function along curves in $\Omega$;
  • the bookkeeping of the nonmonotone reference value.

  The projection lemmas, including Lemma 2.2 (i), are reusable for any gradient projection method and are welcome as separate contributions.

Selected references

• E. G. Birgin, J. M. Martínez, M. Raydan, Nonmonotone spectral projected gradient methods on convex sets, SIAM J. Optim. 10(4) (2000) 1196–1211; authors' updated version, July 2004. https://doi.org/10.1137/S1052623497330963, https://www.ime.unicamp.br/~martinez/bmr.pdf
• E. G. Birgin, J. M. Martínez, M. Raydan, Inexact spectral projected gradient methods on convex sets, IMA J. Numer. Anal. 23 (2003) 539–559. https://doi.org/10.1093/imanum/23.4.539
• D. P. Bertsekas, Nonlinear Programming, Athena Scientific, 1995, Section 2.3.
• D. P. Bertsekas, On the Goldstein–Levitin–Polyak gradient projection method, IEEE Trans. Automat. Control 21 (1976) 174–184. https://doi.org/10.1109/TAC.1976.1101194
• J. Barzilai, J. M. Borwein, Two-point step size gradient methods, IMA J. Numer. Anal. 8 (1988) 141–148. https://doi.org/10.1093/imanum/8.1.141
• L. Grippo, F. Lampariello, S. Lucidi, A nonmonotone line search technique for Newton's method, SIAM J. Numer. Anal. 23 (1986) 707–716. https://doi.org/10.1137/0723046
• M. Raydan, The Barzilai and Borwein gradient method for the large scale unconstrained minimization problem, SIAM J. Optim. 7 (1997) 26–33. https://doi.org/10.1137/S1052623494266365