The OR Formalization Drive

Motivation

Quadratic programming, the minimization of a convex quadratic function subject to linear constraints, is the simplest nonlinear extension of linear programming. It arises in least-squares estimation with sign constraints, in portfolio selection, and as the subproblem solved at each iteration of Newton-type methods for general smooth convex programs. E. M. L. Beale's 1955 paper (DOI 10.1111/j.2517-6161.1955.tb00191.x) gave one of the first finite algorithms for it by extending Dantzig's simplex method: the method keeps the simplex tableau and adds free variables, linear functions of the original variables with no sign restriction, along which the quadratic stops decreasing.

Timeline:

• 1951: Dantzig publishes the simplex method for linear programming.
• 1952: Charnes introduces ε-perturbations to resolve degeneracy in the simplex method.
• 1955: Beale extends the simplex method to convex quadratic objectives and proves that the iteration terminates (§3 of the paper; the result formalized here).
• 1959: Beale's "On quadratic programming" (Naval Research Logistics Quarterly 6) develops the method further; Wolfe's simplex method for quadratic programming (Econometrica 27) appears the same year.

Setting

There are $n$ restricted variables $x_j \ge 0$ satisfying $m$ linearly independent linear equations, and a convex quadratic objective $C$. The iteration keeps $N = n - m$ nonbasic variables $z_1, \dots, z_N$, each either a restricted variable or a free variable, and writes every restricted variable as an affine function of them:

$$x_h = a_{h0} + \sum_{l=1}^{N} a_{hl} z_l. \qquad (2.3)$$

A restricted variable that is not nonbasic is basic. The associated solution sets every $z_l = 0$, so $x_h = a_{h0}$. The objective is written as

$$C = \sum_{k=0}^{N} \sum_{l=0}^{N} c_{kl} z_k z_l, \qquad z_0 = 1, \qquad (3.1)$$

with $(c_{kl})$ symmetric. Thus $c_{00}$ is the value of $C$ at the associated solution and $2c_{k0}$ is its linear coefficient in $z_k$. The number of nonbasic free variables is $s$.

One step chooses a nonbasic $z_p$ that can profitably be altered: a free one with $c_{p0} \ne 0$ if there is one, otherwise a restricted one with $c_{p0} < 0$. It orients $z_p$ so that it is to be increased, and increases it from $0$. It stops at the first of two events. Either a basic variable $x_q$ reaches $0$ (the ratio test (2.4)), and then $x_q$ becomes nonbasic in place of $z_p$. Or $C$ stops decreasing where the free variable $u_r = c_{p0} + \sum_l c_{pl} z_l$ vanishes (3.2), and then $u_r$ becomes nonbasic in place of $z_p$. The coefficients are then transformed by substituting for $z_p$ (eqs. (3.4)–(3.6)). $C$ is in standard form when it has no linear term in any free variable.

Formalization targets

Goal: the iteration terminates

From a tableau with symmetric $(c_{kl})$, positive semidefinite quadratic block $(c_{kl})_{k,l \ge 1}$ and consistent labels, there is no infinite run

$$T_0 \to T_1 \to T_2 \to \cdots$$

of steps along which every basic variable stays strictly positive in the associated solution. No bound on the number of steps is claimed, as in the paper.

Milestones

1. Eq. (3.7): the closed form of the transformed matrix, its symmetry, and the invariance $\sum c_{kl} z_k z_l = \sum c''_{kl} z'_k z'_l$.
2. Lemma 1: when a free variable enters, its row and column vanish off the diagonal, the index $0$ included.
3. Lemma 2: a slot whose row and column vanish off the diagonal keeps this property when another free variable enters.
4. The optimality criterion (p. 175): if no nonbasic variable can profitably be altered and $C$ is convex, then $c_{00}$ is the minimum over the feasible region.
5. In standard form, $c_{00} \le C(z)$ for every $z$ with the restricted nonbasic variables at $0$.
6. $C$ decreases at every step: $c'_{00} < c_{00}$.
7. A finite run never returns to a standard form with the same set of restricted nonbasic variables.
8. If $C$ is not in standard form and $s = s_0$, then within $s_0$ steps either standard form is reached or $s$ drops, and $s$ never exceeds $s_0$ on the way.

Significance

The theorem makes Beale's method an algorithm: a finite procedure that ends either at an optimal tableau (milestone 4) or with a ray along which $C$ decreases without bound. This finiteness is what later active-set methods for quadratic programming inherit.

The result was proved in 1955. What remains is to formalize it: a machine-checked account of a simplex-type method whose state includes variables that are created during the run and later discarded. Mathlib has no simplex-type algorithm for quadratic programming, and no machine-checked proof of this theorem is known. The pivot algebra (3.4)–(3.7) and the tableau model are reusable for other pivoting methods for quadratic programs.

Difficulty

The argument for linear programming does not carry over. There, the objective strictly decreases and a basis is a subset of a finite set of columns, so no basis repeats. Here each step may create a new free variable, and nothing bounds the number of distinct free variables that can occur. Tableaux are therefore not drawn from a finite set, and a strictly decreasing objective alone does not give termination. The paper states this itself: "there is no obvious limit to the number of free variables that may be involved". The difficulty is to bound the number of steps between returns to a well-behaved tableau, and this depends both on the rule that free variables are chosen first and on how the coefficient matrix evolves under repeated pivots.

Formalization scope

• Representation. The nonbasic variables occupy fixed slots Fin (N+1). Slot 0 is $z_0 = 1$; the nonbasic slot k : Fin N is index k.succ. A pivot stores the new nonbasic variable in the slot of the variable it replaces, so the paper's index $q$ in (3.4)–(3.7) is that slot. The tableau holds the labels (restricted $x_j$ or free), the rows of all $n$ restricted variables (a nonbasic one has the unit row), and $(c_{kl})$. Free variables carry no row, as in the paper.
• The pivot. pivotC is computed literally from (3.5) and then (3.6). Rows are transformed by the same substitution, as the paper states.
• The step. The step is a relation. It allows any profitable choice of $z_p$ subject to the free-first rule, and at a tie either outcome. No pricing rule is fixed, since the paper fixes none.
• Convexity. Convexity of $C$ is the symmetry of $(c_{kl})$ plus positive semidefiniteness of the block $(c_{kl})_{k,l \ge 1}$, assumed on the initial tableau.
• Added hypothesis. The one hypothesis not on the page is that every basic restricted variable is strictly positive in the associated solution of every tableau of the run. It replaces Charnes's ε-perturbations, by which the paper ensures "the $a_{h0}$ are always positive, and not zero". Positivity is required of basic variables only; nonbasic variables are $0$ in the associated solution.
• Out of scope. The link to the original equations (2.1) and phase 1 (artificial variables, the M-method) are not formalized: the iteration starts from a tableau already in the form (2.3).
• Ruling out a trivial goal. A step relation that never fires, or a positivity hypothesis that no tableau can meet after a step, would make the goal trivially true. A sorry-free check exhibits a convex instance with consistent labels, a step, and positive basic variables before and after it.

Contributions are welcome on every milestone. The algebraic milestones 1–3 are self-contained.

Selected references

• E. M. L. Beale, On Minimizing a Convex Function Subject to Linear Inequalities, Journal of the Royal Statistical Society, Series B 17(2):173–184, 1955. https://doi.org/10.1111/j.2517-6161.1955.tb00191.x
• A. Charnes, Optimality and Degeneracy in Linear Programming, Econometrica 20(2):160–170, 1952. https://doi.org/10.2307/1907845
• G. B. Dantzig, Maximization of a Linear Function of Variables Subject to Linear Inequalities, in T. C. Koopmans (ed.), Activity Analysis of Production and Allocation, Wiley, 1951, pp. 339–347.
• E. M. L. Beale, On Quadratic Programming, Naval Research Logistics Quarterly 6(3):227–243, 1959. https://doi.org/10.1002/nav.3800060305
• P. Wolfe, The Simplex Method for Quadratic Programming, Econometrica 27(3):382–398, 1959. https://doi.org/10.2307/1909468

Motivation

A make-to-stock manufacturer who also sells on a spot market must decide, at every moment, whether to keep producing and whether to accept or reject incoming orders at the prevailing price. Caldentey and Wein (Revenue Management of a Make-to-Stock Queue, Operations Research 54(5), 2006) study this problem in heavy traffic. The limit is a two-dimensional singular control problem for a diffusion: the inventory level and the logarithm of the price move jointly as a correlated Brownian motion, and the controls push the inventory only when it reaches one of two free boundaries. The optimal boundaries are characterized by an elliptic free-boundary problem that the authors could not solve in closed form.

The paper's way forward is an approximation: change the direction of reflection on the boundary so that the stationary distribution of the controlled process becomes an explicit exponential. Proposition 2 states that exponential form, and it turns the free-boundary problem into a calculus-of-variations problem for the two boundary curves. Explicit stationary densities of reflected diffusions in two dimensions are rare; the classical condition for an exponential stationary density of a reflected Brownian motion, and the characterization of the stationary law by a basic adjoint relation, are due to Harrison and Williams, Multidimensional reflected Brownian motions having exponential stationary distributions, Annals of Probability 15, 1987, the reference the paper cites. This mission formalizes the analytic core of Proposition 2: the exponential density satisfies that relation for the reflection field the proposition singles out.

Setting

Points of the plane are $(x,y)$, with $x$ the inventory level and $y$ the logarithm of the price. The limiting process $(\mathcal X,\mathcal Y)$ has drift $(\theta,0)$ and covariance matrix

$$\Sigma=\begin{pmatrix}\sigma^2&\sigma\delta\varrho\\ \sigma\delta\varrho&\delta^2\end{pmatrix},\qquad \sigma>0,\ \delta>0,\ -1<\varrho<1,$$

so its generator is

$$\Gamma=\theta\frac{\partial}{\partial x}+\frac{\sigma^2}{2}\frac{\partial^2}{\partial x^2}+\sigma\delta\varrho\frac{\partial^2}{\partial x\,\partial y}+\frac{\delta^2}{2}\frac{\partial^2}{\partial y^2}.$$

Two curves bound the region where the process lives: the rejection boundary $x=\eta(y)$ (below it, orders are rejected) and the idleness boundary $x=\xi(y)$ (above it, production stops). For $y_{\min}<y_{\max}$ the region is

$$\Omega=\{(x,y):\ y_{\min}<y<y_{\max},\ \eta(y)<x<\xi(y)\},$$

and its boundary splits into four pieces: $x=\eta(y)$, $x=\xi(y)$, $y=y_{\min}$, $y=y_{\max}$. Write $\vec n$ for the inward unit normal on $\partial\Omega$ and $dl$ for arc length. A reflection field $\vec v$ on $\partial\Omega$ gives the direction in which the process is pushed back into $\Omega$. The basic adjoint relation (BAR) of the paper, equation (43), is

$$\int_\Omega \Gamma f\,\pi_\Omega\,ds+\frac12\int_{\partial\Omega}\vec v\cdot\nabla f\,\pi_\Omega\,dl=0\quad\text{for all test functions } f,$$

and the paper cites Harrison and Williams for the fact that the stationary distribution $\pi_\Omega$ of the reflected process satisfies it. Proposition 2 introduces the eigen-decomposition $\Sigma=V'EV$ ($V$ a rotation whose rows are eigenvectors, $E$ diagonal), the whitening map $T=E^{-1/2}V$ and $\Omega^*=T(\Omega)$, and assumes that $T\vec v$ is normal to $\partial\Omega^*$. The exponents are

$$m_x=\frac{2\theta}{\sigma^2(1-\varrho^2)},\qquad m_y=\frac{-2\varrho\theta}{\sigma\delta(1-\varrho^2)}.\tag{47}$$

Formalization targets

Goal: the exponential density satisfies the BAR under conormal reflection

For $\eta,\xi$ continuously differentiable with $\eta<\xi$ on $[y_{\min},y_{\max}]$, $\pi(x,y)=e^{m_xx+m_yy}$, and every $C^2$ function $f$ on $\mathbb R^2$,

$$\int_\Omega \Gamma f\;\pi\,ds+\frac12\int_{\partial\Omega}(\Sigma\vec n)\cdot\nabla f\;\pi\,dl=0 .$$

The boundary integral is written out on the four pieces, with $\vec n\,dl$ equal to $(1,-\eta'(y))\,dy$, $(-1,\xi'(y))\,dy$, $(0,1)\,dx$ and $(0,-1)\,dx$ respectively. The normalizing constant is left out because the relation is linear in $\pi$.

Milestones

1. The interior equation: $\Gamma^*\pi=-\theta\pi_x+\frac{\sigma^2}{2}\pi_{xx}+\sigma\delta\varrho\,\pi_{xy}+\frac{\delta^2}{2}\pi_{yy}=0$ everywhere.
2. The zero-flux identity: $\frac12\Sigma\nabla\pi=(\theta,0)\,\pi$ everywhere.
3. The meaning of the hypothesis: with $T=E^{-1/2}V$, $(Tv)\cdot(Tw)=v\cdot\Sigma^{-1}w$, and for $n\neq0$, $Tv$ is orthogonal to $T$ of every vector orthogonal to $n$ exactly when $v$ is a multiple of $\Sigma n$.
4. The normalizing constant: $\pi$ is integrable on $\Omega$ and a unique $K_\Omega>0$ makes $K_\Omega\pi$ integrate to one.

Significance

For the operations model, Proposition 2 is what makes the problem computable. Once the stationary density is explicit, the long-run average cost of any pair of boundary curves is an explicit integral, and optimizing over $(\eta,\xi)$ becomes a variational problem with Euler–Lagrange equations; the paper's proposed policy and its numerical comparisons all rest on it.

For formalization, the mission produces a machine-checked version of a statement whose proof the paper does not contain (it is in an online companion) and whose hypothesis is stated only in words. The formal statements fix exactly which reflection field makes the claim true, which the prose leaves ambiguous. None of the statements has, to our knowledge, a machine-checked proof anywhere; the result itself is classical in spirit (an integration by parts on a planar region), but no divergence theorem on a region between two graphs with an anisotropic operator is currently available as a ready-made statement.

Difficulty

The interior equation and the zero-flux identity are finite computations with the exponential. The difficulty is the goal: it is an integration-by-parts identity on a curved planar region with an anisotropic second-order operator. The obvious first step, "apply Green's identity", presupposes a divergence theorem on a region bounded by two graphs $x=\eta(y)$, $x=\xi(y)$ and two horizontal segments, with the boundary integral written in the parametrization of each piece and the orientation of every normal tracked. Mathlib has the divergence theorem on rectangular boxes, not on such regions, and the moving limits $\eta(y)$, $\xi(y)$ are exactly where the terms in $\eta'$ and $\xi'$ of the boundary integral come from.

The second trap is the reflection field. The page describes the modification as substituting the inward unit normal $\vec n$ for $\vec v$; with $\vec v=\vec n$ the identity is false as soon as $\Sigma$ is not a multiple of the identity (on a random instance the residual is of order one). Only the conormal field $\Sigma\vec n$, which is what the hypothesis of Proposition 2 selects, gives a true statement.

Formalization scope

Everything lives in the namespace MakeToStockRM.ExpDensity. The plane is ℝ × ℝ with the inventory first; partial derivatives are Fréchet derivatives applied to (1, 0) and (0, 1), and the mixed partial is $\partial_x(\partial_y f)$. Parameters satisfy $\sigma>0$, $\delta>0$, $|\varrho|<1$; $\theta$ is any real number, and $\theta=0$ (then $\pi\equiv1$) is allowed.

This is the analytic, pinned-down content of Proposition 2. The identification "the BAR characterizes the stationary law of the reflected diffusion" (Harrison–Williams 1987) is out of scope: Mathlib has no reflected Brownian motion. Relative to the page, the formalization commits to the following:

• The reflection field is $\vec v=\Sigma\vec n$ with $\vec n$ the inward unit normal and $dl$ arc length. The hypothesis "$T\vec v$ is normal to $\partial\Omega^*$" fixes only the direction of $\vec v$ (milestone 3); the length $\Sigma\vec n$ is the one for which the BAR holds. The page's phrase "substituting the inward unit normal $\vec n$ for $\vec v$" is inconsistent with the proposition's own hypothesis and is not followed.
• The boundary curves are $C^1$ on $\mathbb R$ with $\eta<\xi$ on $[y_{\min},y_{\max}]$, and $y_{\min}<y_{\max}$, so $\Omega$ is a nonempty bounded region; the paper assumes this implicitly.
• Test functions are all $C^2$ functions on $\mathbb R^2$, which are bounded with bounded derivatives on the closure of $\Omega$ (the paper's "twice continuous and bounded").
• The constant $K_\Omega$ is dropped from the goal and treated in milestone 4.

The goal quantifies over every $C^2$ test function; restricting to functions supported inside $\Omega$ would delete the boundary term and reduce the goal to milestone 1, and that trivialization is ruled out. The second half of Proposition 2 ("(45)–(46) is equivalent to (48)–(49)"), Proposition 1, the heavy-traffic limit, the HJB equation and the proposed policy are not formalized: their normalizations or proofs are only in the online companion.

A complete development needs a divergence theorem on regions between two $C^1$ graphs, which is reusable for any planar PDE statement on such regions. Contributions of that lemma, and of the four milestones, are welcome.

Selected references

• R. Caldentey, L. M. Wein, Revenue Management of a Make-to-Stock Queue, Operations Research 54(5):859–875, 2006. https://doi.org/10.1287/opre.1060.0289
• J. M. Harrison, R. J. Williams, Multidimensional reflected Brownian motions having exponential stationary distributions, Annals of Probability 15(1):115–137, 1987. https://doi.org/10.1214/aop/1176992259
• F. John, Partial Differential Equations, 4th ed., Springer, 1982. https://doi.org/10.1007/978-1-4684-9333-7

Motivation

In 1955 E. M. L. Beale showed how Dantzig's simplex method, which was built for linear objectives, can be carried over to certain nonlinear convex objectives that are minimized subject to linear inequalities (Beale 1955). Section 4 of that paper treats one such objective: the sum of the $t$ largest of a set of $g$ linear forms. Beale's motivation comes from the theory of games: "if the enemy has to choose $t$ out of a set of $g$ possible actions, and $L_f$ represents his average gain through using the $f$th", then the defender wants to minimize the sum of the $t$ largest $L_f$.

The same objective can be written as a linear program. One introduces a bound $u$ and requires every sum of $t$ forms to be at most $u$. That formulation has $\binom{g}{t}$ constraints, which is unwieldy once $t>1$ and $g$ is large. Beale's alternative works with the nonlinear objective directly, and he needs a test that tells him when the current basic solution is already optimal. This mission formalizes that test, Theorem 1 of the paper.

The objective reappears in later work under other names: the sum of the $k$ largest components of a vector, the "top-$k$ sum", and $k$ times the conditional value-at-risk of an empirical distribution. Beale's paper is an early source for its optimality conditions.

Setting

There are real variables $z_l$, indexed by $l$ in a finite set (possibly empty), and $u_1,\dots,u_s$. Two linear forms in these variables are given,

$$A=A_0+\sum_l A_l z_l+\sum_{f=1}^{s}\varphi_f u_f,\qquad L_0=c_{00}+\sum_l c_{0l} z_l+\sum_{f=1}^{s}\theta_f u_f,$$

together with $s$ further forms

$$L_f=L_0-u_f\qquad(f=1,\dots,s).$$

For an integer $\tau\ge0$ the objective is

$$C=A+\bigl(\text{sum of the }\tau\text{ largest of }L_0,L_1,\dots,L_s\bigr).$$

The sum of the $\tau$ largest of $s+1$ numbers is the largest total of any $\tau$ of them. Ties do not make it ambiguous.

The feasible region is fixed by a set $F$ of indices. The variables $z_l$ with $l\in F$ and all the $u_f$ are free, and every other $z_l$ is restricted to $z_l\ge0$. At the origin $z=0$, $u=0$ all $s+1$ forms are equal to $c_{00}$, so the origin is where $C$ fails to be differentiable. In Beale's algorithm the origin is the current basic solution: the $u_f$ measure how far the "borderline" forms sit from a chosen critical form, and $A$ collects the forms that are certainly among the largest.

Write $a_l=A_l+\tau c_{0l}$ and $w_f=\varphi_f+\tau\theta_f$.

Formalization targets

Goal: Theorem 1 (a), p. 179

For $\tau\le s$, $C$ is minimized over the feasible region when all the $z_l$ and $u_f$ vanish if and only if

$$\begin{aligned} &a_l\ge0\ \text{for all } l, \qquad a_l=0\ \text{for all } l\in F,\\ &0\le w_f\le1\ \text{for all } f,\qquad \tau-1\le\sum_{f=1}^{s}w_f\le\tau . \end{aligned}\tag{4.5}$$

"Minimized" means a global minimum: $C(0,0)\le C(z,u)$ at every feasible point.

Milestones

1. Convexity (p. 179). $C$ is a convex function of $(z,u)$ for $\tau\le s+1$.
2. Descent rules (second half of Theorem 1 (a), p. 179). When a condition of (4.5) fails, a stated move of one variable, or of all $u_f$ together, lowers $C$ below $C(0,0)$ for every small enough step. There are six moves: $z_l\uparrow$ if $a_l<0$; $z_l\downarrow$ if $a_l>0$ and $l\in F$; $u_f\uparrow$ if $w_f<0$; $u_f\downarrow$ if $w_f>1$; all $u_f\uparrow$ if $\sum w_f<\tau-1$; all $u_f\downarrow$ if $\sum w_f>\tau$.
3. The rearrangement identity (proof of Theorem 1 (a), p. 180). If $1\le\tau\le s$, $u'_1\le\dots\le u'_s$ and $u'_\tau\le0$, then

$$C=A_0+\tau c_{00}+\sum_l a_l z'_l+\sum_{f=1}^{\tau}(w_f-1)(u'_f-u'_\tau)+\sum_{f=\tau+1}^{s}w_f(u'_f-u'_\tau)+\Bigl\{\sum_{f=1}^{s}w_f-\tau\Bigr\}u'_\tau .$$

4. Theorem 1 (b) (p. 180). For $\tau=s+1$, the origin is a minimum if and only if (4.5) holds and $w_f=1$ for every $f$. Otherwise some value of $u_f$ with the sign opposite to $w_f-1$ lowers $C$.

Significance

Theorem 1 is the optimality test of Beale's simplex method for the sum-of-largest objective. The algorithm on pp. 178–179 changes nonbasic variables one at a time. When no single change is profitable it applies Theorem 1: either (4.5) holds and the current solution is optimal, or one of the six descent rules names the variable to change next. The test is exact even though the objective is not differentiable at the current point. It is a closed-form description of the subdifferential of a top-$\tau$ sum at a point where all the forms tie. The theorem is also the base case of the multi-group generalization that Beale mentions on p. 181.

The paper proves Theorem 1 by hand. To our knowledge neither the theorem nor the rearrangement identity behind it has been formalized in any proof assistant. The mission produces:

• a checked statement and proof of the test, including the degenerate cases $\tau=0$ and $s=0$, which the paper does not discuss separately;
• the boundary case $\tau=s+1$;
• a reusable Lean definition of the sum of the $\tau$ largest entries of a finite real family, with its convexity.

Difficulty

Necessity, the "only if" direction, is the part the paper calls obvious: each descent rule changes $C$ linearly for small steps. Two features still have to be handled explicitly. The step must be small only in rule-dependent ways, and the ordering of the forms changes along the moves of rules 4 and 6.

Sufficiency is where the work lies. The naive argument, "the directional derivative in every coordinate direction is non-negative, so the origin is a minimum", fails because $C$ is not differentiable at the origin. Nonnegative derivatives along the coordinate axes do not control mixed directions in which several $u_f$ move by different amounts, which reorders the forms. Which $\tau$ forms are the largest then depends on the point, and the paper settles the configurations in which $L_0$ is among the $\tau$ largest by an informal appeal to the "essential symmetry" between $L_0$ and the other forms. A formal proof cannot leave that appeal informal: the forms are parametrised relative to $L_0$ (each $L_f$ is $L_0-u_f$), so the symmetry is a change of variables that has to be written down and shown to preserve (4.5).

Formalization scope

• Data. The variables are z : Fin r → ℝ (any r, including $0$) and u : Fin s → ℝ. The paper's $u_f$ for $f=1,\dots,s$ is Lean's u f for $f=0,\dots,s-1$. The coefficients $(A_0,A_l,\varphi_f,c_{00},c_{0l},\theta_f)$ form a structure Forms r s.
• Forms. The family $L_0,\dots,L_s$ is Fin (s+1) → ℝ, with index $0$ for $L_0$ and index f.succ for $L_0-u_f$. The free set $F$ is a Finset (Fin r), and $\tau$ is a natural number cast to $\mathbb R$ wherever it multiplies a coefficient.
• Sum of the largest. sumLargest τ v is the maximum over $\tau$-element subsets $S$ of $\sum_{i\in S}v_i$ (Finset.sup' over powersetCard). It is the junk $0$ for $\tau$ larger than the number of entries, a case no statement uses.
• Minimality. "Minimized when all variables vanish" is the global statement $C(0,0)\le C(z,u)$ for all $(z,u)$ with $z_l\ge0$ for $l\notin F$. It is not a local minimum, and the sign constraints on restricted $z_l$ are kept: they are why the first condition of (4.5) is an inequality.
• Descent. "$C$ can be decreased by moving $x$ from zero" is a strict decrease for all step sizes in some interval $(0,\varepsilon)$, with every other variable at zero.
• No trivialization. The goal is an equivalence with no hypothesis beyond $\tau\le s$. Neither direction can be satisfied vacuously, and the cases $\tau=0$ and $s=0$ are included, as on the page.
• Added hypotheses. The rearrangement milestone assumes $\tau\ge1$, because the paper's $u'_\tau$ does not exist at $\tau=0$. Its second line uses $c_{0l}$ where the page misprints $c_l$.

Needed infrastructure:

• basic lemmas on sumLargest: its value at a constant family, at a family sorted by a monotone shift, and under adding a common constant;
• the change of variables behind the paper's symmetry between $L_0$ and the other forms.

These lemmas are reusable for any top-$k$-sum or empirical-CVaR objective. Contributions are welcome at any level: lemmas about sumLargest, any of the milestones, or an alternative sufficiency proof through convexity and one-sided directional derivatives.

Not in scope: the pivoting rules (4.2)–(4.4), the degeneracy discussion on pp. 180–181, and the multi-group generalization, which the paper says is "cumbersome to state" and does not state.

Selected references

• E. M. L. Beale, On Minimizing a Convex Function Subject to Linear Inequalities, Journal of the Royal Statistical Society, Series B 17(2), 173–184, 1955. https://doi.org/10.1111/j.2517-6161.1955.tb00191.x
• G. B. Dantzig, A. Orden and P. Wolfe, The generalized simplex method for minimizing a linear form under linear inequality restraints, Pacific Journal of Mathematics 5(2), 183–195, 1955. https://doi.org/10.2140/pjm.1955.5.183
• R. T. Rockafellar and S. Uryasev, Optimization of conditional value-at-risk, Journal of Risk 2(3), 21–41, 2000. https://doi.org/10.21314/JOR.2000.038

Motivation

A linear program is solved with known data, but in planning problems the data are often only known in distribution when the main decision is taken: demands, yields and requirements are revealed later, and a corrective action is taken after they are. E. M. L. Beale's 1955 paper On Minimizing a Convex Function Subject to Linear Inequalities formulates this situation in its §5, "Linear Programming with Random Coefficients", as what is now called a two-stage stochastic linear program with recourse. Beale's motivating example is the transportation problem of Hitchcock (1941) with random requirements at the destinations, where every unit of shortage or excess incurs a loss. The same model was put forward in the same year by Dantzig, Linear Programming under Uncertainty (Management Science, 1955), as the paper's note added in proof acknowledges.

Timeline:

• 1955. Beale (§5, Theorems 2 and 3) and Dantzig independently introduce two-stage linear programs with random data; Beale proves that the expected cost is convex in the first-stage decision, and that the cost is convex in the random data for fixed decision.
• 1967. Walkup and Wets, Stochastic Programs with Recourse, study the domain of the expected recourse function and its properties under fixed recourse.
• 1974. Wets, Stochastic Programs with Fixed Recourse: The Equivalent Deterministic Program, gives the systematic treatment of convexity, finiteness and polyhedrality of the expected recourse function, now textbook material (Birge and Louveaux, Introduction to Stochastic Programming, Ch. 3).

Setting

Constants $c\in\mathbb R^n$, $f\in\mathbb R^p$ and an $m\times p$ matrix $D=(d_{ik})$ are given. The data $A=(\alpha_{ij})$, an $m\times n$ matrix, and $\beta\in\mathbb R^m$ are random variables on a probability space $(\Omega,P)$: their distribution is known when the first-stage decision $x\in\mathbb R^n$, $x\ge0$, is chosen, and their values are known when the second-stage decision $y\in\mathbb R^p$, $y\ge0$, is chosen. The cost is

$$C=c'x+f'y,\qquad Ax+Dy=\beta. \qquad(5.3),(5.4)$$

For a right-hand side $b\in\mathbb R^m$ the second-stage value is

$$Q(b)=\min\{f'y : y\ge0,\ Dy=b\},$$

and for fixed data the cost of a first-stage decision is $C(x)=c'x+Q(\beta-Ax)$. The expected cost is

$$E(C)(x)=\int_\Omega \bigl(c'x+Q(\beta(\omega)-A(\omega)x)\bigr)\,dP(\omega).$$

The problem is to choose $x\ge0$ minimising $E(C)$. In Lean the value is secondStageValue D f b, the cost is cost c f D A β x, and the expected cost is expectedCost P c f D A β x, all in the namespace BealeConvexMin.RandomLP.

Formalization targets

Goal: Theorem 2 (p. 182)

Assume that for every $x\ge0$ the second-stage minimum is attained for almost every outcome and that $\omega\mapsto C(x,\omega)$ is integrable. Then

$$E(C)(\lambda_1x_1+\lambda_2x_2)\le\lambda_1E(C)(x_1)+\lambda_2E(C)(x_2)\qquad(x_1,x_2\ge0,\ \lambda_1,\lambda_2\ge0,\ \lambda_1+\lambda_2=1),$$

that is, $E(C)$ is convex on the non-negative orthant. The statement fixes no distribution class: it is claimed for any known distribution of $(A,\beta)$.

Milestones

1. Pointwise convexity (last display of the proof of Theorem 2, p. 182): for fixed data $(A,\beta)$, with the minimum attained at every $x\ge0$,

$$C(\lambda_1x_1+\lambda_2x_2)\le\lambda_1C(x_1)+\lambda_2C(x_2).$$

2. Theorem 3 (p. 182): for fixed $x$, the cost $(A,\beta)\mapsto c'x+Q(\beta-Ax)$ is jointly convex on every convex set of data on which the second-stage minimum is attained.
3. Eqs. (5.5)–(5.6) (p. 182): for a finitely supported distribution, $A=A_r$ and $\beta=\beta_r$ with probability $p_r$, the value $E(C)(x)$ is the minimum of $c'x+\sum_r p_r f'y_r$ over non-negative $y_r$ with $A_rx+Dy_r=\beta_r$ for all $r$; minimising $E(C)$ is then a linear program.

Significance

The result. Theorem 2 is the basic structural fact of two-stage stochastic linear programming: the first-stage problem is a convex program in $x$, whatever the distribution of the data. It is what makes local optimality global for the first-stage problem, what justifies cutting-plane and decomposition methods that approximate $E(C)$ from below by supporting hyperplanes, and what makes sample-average approximations convex programs. Theorem 3, joint convexity in the data, gives through Jensen's inequality the comparison between the stochastic problem and its mean-value problem that Beale draws on p. 182. The discrete reformulation (5.5)–(5.6) is the deterministic-equivalent linear program used for finitely many scenarios.

Formalizing it. The theorems are proved in the paper, and their content is classical. The mission produces machine-checked statements of the model with its implicit hypotheses made explicit (attainment of the second stage, integrability of the cost), and proofs of the three results in Lean. The platform already has related statements in other models (finite scenario sets with extended-real recourse, and a complete-recourse, finite-second-moment version); none has Beale's hypotheses, and none states convexity of $c'x+E\,Q$ for an arbitrary distribution.

Difficulty

The mathematics is short; the difficulty is in the encoding. The second-stage value is a minimum that may fail to exist: the second stage may be infeasible for some $x$ and some outcomes, or unbounded below. A real-valued infimum then takes an arbitrary default value, and convexity would become a statement about that default. Similarly, the mean value only exists when the cost is integrable. A faithful statement has to carry attainment and integrability exactly where the paper tacitly assumes them, on the domain $x\ge0$ the paper uses, and no stronger condition (such as complete recourse or moment bounds) that the paper does not make. In the discrete reformulation, the minimum over the whole family $(y_r)_r$ has to be matched with the probability-weighted sum of per-scenario minima.

Formalization scope

• Vectors are Fin n → ℝ, matrices Matrix (Fin m) (Fin n) ℝ, inner products dotProduct, and $y\ge0$ is the componentwise order. The random data are functions A : Ω → Matrix (Fin m) (Fin n) ℝ and β : Ω → Fin m → ℝ on a measurable space with a probability measure P; no measurability of the data is assumed beyond integrability of the cost.
• The second-stage value is the real infimum of $f'y$ over the feasible set. It equals $0$ on an infeasible or unbounded-below second stage, so each theorem assumes attainment of the minimum where it is evaluated (the paper's "value of $y$ that minimizes $C$"). The goal assumes attainment for almost every outcome at every $x\ge0$.
• $E(C)$ is the Bochner integral, which is $0$ for a non-integrable integrand, so the goal assumes integrability of $C(x,\cdot)$ at every $x\ge0$ (the paper's "mean value $E(C)$").
• Convexity is claimed on $\{x : x\ge0\}$, the paper's domain, not on all of $\mathbb R^n$. Theorem 3 is stated for fixed non-negative $x$ (the model's first-stage domain) and on every convex set of data on which the minimum is attained, since the paper names no domain.
• A formalization in which the value is an unconstrained infimum without attainment, or the expectation is taken without integrability, is trivially convex on the region where the default values apply and does not state Beale's theorem; such variants are ruled out.
• Reusable beyond this mission: basic facts on the optimal value of a parametric linear program in its right-hand side and cost data, and convexity of integrals of pointwise-convex integrands. Proofs of the milestones and of the goal, and alternative formulations in extended reals, are welcome.

Selected references

• E. M. L. Beale, On Minimizing a Convex Function Subject to Linear Inequalities, Journal of the Royal Statistical Society, Series B 17(2):173–184, 1955. https://doi.org/10.1111/j.2517-6161.1955.tb00191.x
• G. B. Dantzig, Linear Programming under Uncertainty, Management Science 1(3–4):197–206, 1955. https://doi.org/10.1287/mnsc.1.3-4.197
• F. L. Hitchcock, The Distribution of a Product from Several Sources to Numerous Localities, Journal of Mathematics and Physics 20:224–230, 1941. https://doi.org/10.1002/sapm1941201224
• D. W. Walkup and R. J.-B. Wets, Stochastic Programs with Recourse, SIAM Journal on Applied Mathematics 15(5):1299–1314, 1967. https://doi.org/10.1137/0115113
• R. J.-B. Wets, Stochastic Programs with Fixed Recourse: The Equivalent Deterministic Program, SIAM Review 16(3):309–339, 1974. https://doi.org/10.1137/1016053
• J. R. Birge and F. Louveaux, Introduction to Stochastic Programming, 2nd ed., Springer, 2011. https://doi.org/10.1007/978-1-4614-0237-4

Motivation

Conic quadratic programs (second-order cone programs) minimize a linear objective subject to linear constraints and constraints of the form $\|A_\ell x-b_\ell\|_2\le c_\ell^Tx-d_\ell$. They model robust linear programs with ellipsoidal uncertainty, truss topology design, contact problems with Coulomb friction, and convex quadratically constrained quadratic programs. In theory they are no harder than linear programs of the same size; in practice, linear programming software handles far larger instances than conic quadratic solvers did at the time of writing (Ben-Tal & Nemirovski 2001, pp. 193–195). This raises a question about geometry rather than algorithms: can a second-order cone be replaced by a polyhedral cone of moderate size without losing much accuracy?

The obvious answer — circumscribe the cone by a polyhedral cone with many facets — fails: the number of facets must grow exponentially in the dimension, even for constant accuracy. Ben-Tal and Nemirovski showed that auxiliary variables change the picture completely: a projection of a polyhedral cone can approximate the Lorentz cone with size only $O(k\ln(1/\varepsilon))$. The construction is now standard; it underlies, for instance, the lifted linear-programming branch-and-bound algorithm for mixed-integer conic quadratic programs of Vielma, Ahmed & Nemhauser 2008.

Setting

For $y\in\mathbb R^k$ let $\|y\|_2=\sqrt{y_1^2+\dots+y_k^2}$. The $(k+1)$-dimensional Lorentz cone is

$$L^k=\{(y,t)\in\mathbb R^k\times\mathbb R\mid t\ge\|y\|_2\}.$$

Fix $\varepsilon>0$. A polyhedral $\varepsilon$-approximation of $L^k$ is a linear map

$$\Pi(y,t,u):\mathbb R^k\times\mathbb R\times\mathbb R^{p}\to\mathbb R^{q}$$

such that

1. if $(y,t)\in L^k$, there is $u\in\mathbb R^p$ with $\Pi(y,t,u)\ge0$ (componentwise);
2. if $\Pi(y,t,u)\ge0$ for some $u$, then $\|y\|_2\le(1+\varepsilon)t$.

Equivalently, the polyhedral cone $\{(y,t,u)\mid\Pi(y,t,u)\ge0\}$ projects onto a cone lying between $L^k$ and its $(1+\varepsilon)$-extension. The size of the approximation is $p+q$: the number of auxiliary variables plus the number of linear inequalities (an equation counts as two).

The construction in the paper uses a tower of variables: for $k=2^\theta$, the coordinates $y_1,\dots,y_k$ form generation $0$, each consecutive pair of generation $\ell-1$ has a successor in generation $\ell$ ($y_i^\ell$ has parents $y_{2i-1}^{\ell-1},y_{2i}^{\ell-1}$), and the single variable of generation $\theta$ is $t$. It also uses an explicit linear system (8) in variables $\xi^j,\eta^j$, $j=0,\dots,\nu$, with trigonometric coefficients $\cos(\pi/2^{j+1})$, $\sin(\pi/2^{j+1})$, $\tan(\pi/2^{\nu+1})$, whose accuracy is $\delta(\nu)=1/\cos(\pi/2^{\nu+1})-1$.

Formalization targets

Goal: Theorem 1.1

There is an absolute constant $C$ such that for every positive integer $k$ and every $\varepsilon\in(0,1]$, $L^k$ admits a polyhedral $\varepsilon$-approximation with

$$p_k+q_k\le C\,k\ln\frac{2}{\varepsilon}.$$

The constant is not fixed; the goal asserts only the order of growth, which is what the paper claims.

Milestones

1. §2, Eq. (5). For $k=2^\theta$, $\theta\ge1$: $(y,t)$ extends to a tower solving $\sqrt{[y_{2i-1}^{\ell-1}]^2+[y_{2i}^{\ell-1}]^2}\le y_i^\ell$ for all $i,\ell$ if and only if $\|y\|_2\le t$.
2. §2, Eqs. (6)–(7). Placing polyhedral $\varepsilon_\ell$-approximations of $L^2$ on every level of the tower yields a polyhedral approximation of $L^k$ with $1+\varepsilon=\prod_{\ell=1}^\theta(1+\varepsilon_\ell)$.
3. Proposition 2.1 (i), (ii) and Eq. (9). System (8) is a polyhedral $\delta(\nu)$-approximation of $L^2$, and $\delta(\nu)=O(4^{-\nu})$.
4. Proof of Theorem 1.1, system (10), property 3. System (8) with parameter $\nu_\ell$ on level $\ell$ of the tower approximates $L^{2^\theta}$ with quality $\beta=\prod_{\ell=1}^\theta 1/\cos(\pi/2^{\nu_\ell+1})-1$.
5. Proof of Theorem 1.1, choice of $\nu_\ell$. With $\nu_\ell=\lfloor c\,\ell\ln(2/\varepsilon)\rfloor$: $\beta\le\varepsilon$ and $\sum_\ell 2^{\theta-\ell}\nu_\ell\le C\,2^\theta\ln(2/\varepsilon)$.

Significance

The theorem shows that conic quadratic constraints are, up to a factor logarithmic in the accuracy, no more expensive to express as linear constraints than they are in their native form. Consequences listed in the paper include approximating convex quadratically constrained quadratic programs, robust counterparts of linear programs with ellipsoidal uncertainty, and problems with low-dimensional cones (Coulomb friction, $k\le3$; truss design, $k\le2$) by linear programs of comparable size. Together with the matching lower bound of §3 of the same paper (a separate mission of this series), it pins down the size of the best polyhedral approximation up to constants. The recursive halving of dimensions through the tower of 3-dimensional cones is a reusable device for other rotation-invariant cones.

The result is proved in the paper; as far as is known it has not been machine-checked. This mission produces a formal proof of the construction, including the trigonometric estimate $\delta(\nu)=O(4^{-\nu})$ and the explicit linear encoding with its size count. Explicit values of the absolute constants are welcome as additional results.

Difficulty

The planar estimate is the core. Part (ii) of Proposition 2.1 must hold for every solution of the inequality system (8), not only for the solution one would write down for a given point of $L^2$; an argument that tracks only the intended solution proves part (i) and nothing about part (ii). The accuracy must also come out as $1/\cos(\pi/2^{\nu+1})-1$, geometric in $\nu$; a bound that decays only polynomially in $\nu$ would give size $\mathrm{poly}(1/\varepsilon)$ instead of $\ln(1/\varepsilon)$. The naive idea of approximating $L^k$ directly by tangent hyperplanes is ruled out by the exponential facet count mentioned above; the auxiliary variables are indispensable. The second difficulty is bookkeeping: packaging $k-1$ copies of system (8) on a tower of depth $\theta=\log_2 k$ into a single linear map, counting its variables and inequalities exactly, handling $k$ that is not a power of two, and summing the accuracies so that the total size is $O(k\ln(2/\varepsilon))$ rather than $O(k\ln k\ln(1/\varepsilon))$.

Formalization scope

• Vectors of $\mathbb R^k$ are Fin k → ℝ. The norm $\|y\|_2$ is written out as eucNorm y = Real.sqrt (∑ i, y i ^ 2); the norm Mathlib attaches to Fin k → ℝ is the sup norm, under which the cone would be polyhedral and the theorem trivial.
• A polyhedral approximation is an $\mathbb R$-linear map (Fin k → ℝ) × ℝ × (Fin p → ℝ) →ₗ[ℝ] (Fin q → ℝ) and $\ge0$ is the componentwise order. Linearity is essential: with an arbitrary map, $\Pi(y,t)=t-\|y\|_2$ would be an exact approximation with $p=0$, $q=1$. Affine maps are not allowed either; the paper's approximations are homogeneous.
• The paper's absolute constants $O(1)$ are existential constants quantified before $k$, $\varepsilon$ and $\theta$. The goal requires $k\ge1$ and $\varepsilon\in(0,1]$, as in the paper; $\ln$ is Real.log.
• System (8) and system (10) are stated as propositions with the absolute values written out; their parameters $\nu$, $\nu_\ell$ are required to be positive integers, as in the paper (at $\nu=0$ the coefficient $\tan(\pi/2)$ would be evaluated as $0$ by Lean).
• Tower variables are indexed Y ℓ i with 0-based i, so the parents of Y ℓ i are Y (ℓ-1) (2i) and Y (ℓ-1) (2i+1); the milestones on (6)–(7) and (10) are stated on solution sets rather than on an explicit linear map. The size counts of (10) (properties 1–2) are not separate milestones; the arithmetic milestone on $\nu_\ell$ records the bound on $\sum_\ell 2^{\theta-\ell}\nu_\ell$ to which they reduce.
• $\delta(\nu)=O(1/4^\nu)$ is stated as $\exists C>0,\ \forall\nu\ge1,\ \delta(\nu)\le C/4^\nu$.

A complete development needs: elementary trigonometry of $\pi/2^{j}$ (available in Mathlib), rotations in the plane, finite products and sums over $\{1,\dots,\theta\}$, and a way to assemble many small linear systems into one linear map with an exact count of rows and columns. The last piece, and the tower of variables with the reduction from arbitrary $k$ to a power of two, are reusable for other lifted polyhedral approximations. Contributions of any milestone, of explicit linear encodings of (8) and (10), and of the extension from $k=2^\theta$ to all $k$ are welcome.

Selected references

• A. Ben-Tal and A. Nemirovski, On Polyhedral Approximations of the Second-Order Cone, Mathematics of Operations Research 26(2):193–205, 2001. https://doi.org/10.1287/moor.26.2.193.10561
• J. P. Vielma, S. Ahmed and G. L. Nemhauser, A lifted linear programming branch-and-bound algorithm for mixed-integer conic quadratic programs, INFORMS Journal on Computing 20(3):438–450, 2008. https://doi.org/10.1287/ijoc.1070.0256
• A. Ben-Tal and A. Nemirovski, Robust convex optimization, Mathematics of Operations Research 23(4):769–805, 1998. https://doi.org/10.1287/moor.23.4.769
• A. Ben-Tal and A. Nemirovski, Lectures on Modern Convex Optimization, SIAM, 2001. https://doi.org/10.1137/1.9780898718829

Motivation

Many problems in convex optimization, variational inequalities and signal processing reduce to a monotone inclusion: find a point $x$ at which the sum of several monotone operators contains $0$. When the sum has two terms, the classical operator-splitting methods (Douglas–Rachford, forward–backward, forward–backward–forward) solve it by iterating a fixed-point map that uses each operator separately, through its resolvent or through a forward (explicit) step. Problems with three terms, for instance a smooth loss plus two nonsmooth regularizers or constraints, are common in practice, and before 2015 no fixed-point map was known that handled three operators one at a time without a product-space reformulation.

Davis and Yin (Set-Valued Var. Anal. 25 (2017) 829–858; preprint arXiv:1504.01032) introduced such a map, now called Davis–Yin three-operator splitting. It contains Douglas–Rachford splitting ($C = 0$) and forward–backward splitting ($B = 0$) as special cases, and it has become a standard building block of first-order methods for composite optimization. This mission formalizes Section 2 of the paper: the fixed-point encoding, the averagedness of the map, and the weak and strong convergence of the resulting iteration.

Setting

Let $H$ be a real Hilbert space. A set-valued operator $A : H \to 2^H$ is monotone if $\langle x - y, u - v\rangle \ge 0$ for all $u \in Ax$, $v \in Ay$, and maximal monotone if its graph is not properly contained in the graph of another monotone operator. Its domain is $\operatorname{dom}(A) = \{x : Ax \ne \emptyset\}$ and the zero set of an operator $M$ is $\operatorname{zer}(M) = \{x : 0 \in Mx\}$. A single-valued $C : H \to H$ is $\beta$-cocoercive ($\beta > 0$) if $\beta\|Cx - Cy\|^2 \le \langle Cx - Cy, x - y\rangle$ for all $x, y$.

Problem (1.1) is: given maximal monotone $A, B$ and $\beta$-cocoercive $C$, find

$$x \in H \quad\text{with}\quad 0 \in Ax + Bx + Cx .$$

For $\gamma > 0$ the resolvent $J_{\gamma A} = (I + \gamma A)^{-1}$ is the map with $x \in J_{\gamma A}x + \gamma A(J_{\gamma A}x)$. The Davis–Yin operator (Eq. (1.2)) is

$$T := J_{\gamma A} \circ (2J_{\gamma B} - I - \gamma C \circ J_{\gamma B}) + I - J_{\gamma B}.$$

Algorithm 1 starts from $z^0 \in H$ and, for relaxation parameters $\lambda_k > 0$, iterates

$$x_B^k = J_{\gamma B}(z^k),\qquad x_A^k = J_{\gamma A}(2x_B^k - z^k - \gamma Cx_B^k),\qquad z^{k+1} = z^k + \lambda_k(x_A^k - x_B^k),$$

so that $z^{k+1} = (1 - \lambda_k)z^k + \lambda_k Tz^k$. A sequence converges weakly, $u_k \rightharpoonup u$, if $\langle u_k, y\rangle \to \langle u, y\rangle$ for every $y \in H$.

Formalization targets

Goal: Theorem 2.1 (Main convergence theorem)

Fix $\varepsilon \in (0,1)$, $\gamma \in (0, 2\beta\varepsilon)$, $\alpha = 1/(2-\varepsilon)$ and $\lambda_k \in (0, 1/\alpha)$ with $\sum_k \tau_k = \infty$, where $\tau_k = \lambda_k(1-\lambda_k) + \lambda_k(1-\alpha)/\alpha$, and $\inf_k \lambda_k > 0$. If $\operatorname{Fix} T \ne \emptyset$, there is $z^* \in \operatorname{Fix} T$ with $z^k \rightharpoonup z^*$ and

$$Cx_B^k \to Cx^* \ \ (\forall x^* \in \operatorname{zer}(A+B+C)),\qquad x_B^k \rightharpoonup J_{\gamma B}(z^*) \in \operatorname{zer}(A+B+C),\qquad x_A^k \rightharpoonup J_{\gamma B}(z^*),$$

and if $A$ or $B$ is uniformly monotone on every nonempty bounded subset of its domain, or $C$ is demiregular at every zero of $A + B + C$, then $x_B^k$ and $x_A^k$ converge strongly to a common point of $\operatorname{zer}(A + B + C)$.

Milestones

In the order the proof uses them: Lemma 2.1 (the identities for one application of $T$), Lemma 2.2 ($\operatorname{zer}(A+B+C) = J_{\gamma B}(\operatorname{Fix} T)$), Lemma 2.3 (inequality (2.1)), Proposition 2.1 ($T$ is $2\beta/(4\beta-\gamma)$-averaged, inequality (2.2)), Remark 2.1 (the strengthened inequality (2.4)), Corollary 2.1 Parts 1–3 (Fejér monotonicity, vanishing residual, weak convergence of $z^k$), Corollary 2.1 Part 4 (the residual rates $\|Tz^k - z^k\|^2 \le \|z^0 - z^*\|^2/(\underline\tau(k+1))$ and $o(1/(k+1))$), and Eqs. (2.6)–(2.7) (the per-step descent inequality and its summed form).

Significance

Theorem 2.1 is the basic convergence guarantee for three-operator splitting: it certifies that the computable sequences $x_B^k$, $x_A^k$, not only the auxiliary sequence $z^k$, approach a solution of (1.1). In infinite dimensions this is the delicate part: for Douglas–Rachford splitting ($C = 0$) weak convergence of the shadow sequence $J_{\gamma B}(z^k)$ was only established by Svaiter in 2011. The result underlies the convergence of the many algorithms obtained from it by specialization (Douglas–Rachford, forward–backward, and the three-block methods of Section 4 of the paper), and the averagedness coefficient of Proposition 2.1 reduces, for $B = 0$, to the best known one for forward–backward splitting.

All statements of this mission are proved in the paper, partly by appeal to Bauschke and Combettes' monograph (Krasnosel'skiĭ–Mann convergence, the demiclosedness of maximal monotone graphs). None of them has a machine-checked proof: Mathlib has no maximal monotone operators, resolvents, averaged maps or Krasnosel'skiĭ–Mann theorem. The mission therefore produces both a formal proof of the Davis–Yin theorem and a first body of monotone-operator theory in Lean.

Difficulty

The fixed-point part is standard once $T$ is known to be averaged: Krasnosel'skiĭ–Mann theory and Opial's argument give $z^k \rightharpoonup z^*$. The obstacle is transferring this to $x_B^k = J_{\gamma B}(z^k)$. Resolvents are nonexpansive but not weakly continuous, so $z^k \rightharpoonup z^*$ does not imply $J_{\gamma B}(z^k) \rightharpoonup J_{\gamma B}(z^*)$; the naive argument fails at exactly this step. Identifying the weak cluster points of $x_B^k$ requires a closedness property of sums of maximal monotone operators under mixed weak and strong convergence, fed by the strong convergence of $Cx_B^k$, which in turn needs the extra term of (2.4) that (2.2) discards. Strong convergence in Part 2 needs yet another argument for each of the three alternative hypotheses.

Formalization scope

• $H$ is an arbitrary real Hilbert space (NormedAddCommGroup, InnerProductSpace ℝ, CompleteSpace); a finite-dimensional space would identify weak and strong convergence and change the theorems.
• Operators $A, B$ are H → Set H; $C$ is single-valued H → H. The resolvents are not constructed: $J_A, J_B$ are maps satisfying the resolvent inclusion $\gamma^{-1}(x - Jx) \in A(Jx)$, which for maximal monotone operators determines them uniquely and exists by Minty's theorem.
• Weak convergence is $\langle u_k, y\rangle \to \langle u, y\rangle$ for every $y$; strong convergence is norm convergence. Iterates are indexed from $0$.
• The printed hypothesis $\alpha = 1/(2-\varepsilon) < 2\beta/(4\beta-\gamma)$ of Corollary 2.1 and Theorem 2.1 contradicts $\gamma < 2\beta\varepsilon$ (it is a typo for $>$) and is not assumed. The printed $\tau_k = (1-\lambda_k/\alpha)\lambda_k/\alpha$ is replaced by the $\tau_k$ of the proof (p. 836), a weaker hypothesis.
• Uniform monotonicity uses a nondecreasing $\varphi : [0,\infty) \to [0,+\infty]$ with $\varphi(0) = 0$ that vanishes only at $0$, as the proof requires; with $\varphi \equiv 0$ allowed, Part 2(a) would be false.
• The O-constant of Corollary 2.1 Part 4 is explicit, $\|z^0 - z^*\|^2/\underline\tau$, and the little-$o$ is stated as $(k+1)\|Tz^k - z^k\|^2 \to 0$. Eq. (2.7) is stated with a uniform lower bound $\underline\lambda \le \lambda_i$ in place of the printed $\lambda_k$, with summability part of the conclusion.
• A formalization with $T$ an arbitrary averaged map, with resolvents replaced by arbitrary nonexpansive maps, or with the contradictory comparison of $\alpha$ kept as a hypothesis would make the theorem vacuous or different; all three are ruled out.

A complete development needs the basic theory of monotone operators (monotonicity of resolvents' graphs, firm nonexpansiveness of resolvents, weak-to-strong closedness of maximal monotone graphs), Krasnosel'skiĭ–Mann iteration with Opial's lemma, and weak sequential compactness of bounded sets in Hilbert space. All of this is reusable far beyond this mission, and contributions of any of these pieces as separate theorems are welcome.

Selected references

• D. Davis and W. Yin, A Three-Operator Splitting Scheme and its Optimization Applications, Set-Valued and Variational Analysis 25 (2017) 829–858. https://doi.org/10.1007/s11228-017-0421-z
• H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer, 2011. https://doi.org/10.1007/978-1-4419-9467-7
• B. F. Svaiter, On weak convergence of the Douglas–Rachford method, SIAM J. Control Optim. 49 (2011) 280–287. https://doi.org/10.1137/100788100
• D. Davis and W. Yin, Convergence rate analysis of several splitting schemes, in: Splitting Methods in Communication, Imaging, Science, and Engineering, Springer, 2016. https://arxiv.org/abs/1406.4834