Optimization

Motivation

Stochastic optimization methods often see only noisy gradients. For a nonconvex objective, a useful guarantee is that an output point has a small gradient: this is a first-order stationarity certificate, even when there is no claim that the point is globally optimal. A guarantee on the expected squared gradient leaves open how often one particular run succeeds. The two-phase randomized stochastic gradient method of Ghadimi and Lan (2013) addresses that question by generating several candidates and using further oracle calls to select one. The paper gives an explicit bound on the probability of failure and on the number of oracle calls needed for a prescribed accuracy and confidence.

The work sits within stochastic approximation, where the objective may be a population expectation or another function accessible through noisy observations. The 2013 paper develops first- and zeroth-order methods for smooth nonconvex problems. This mission concerns its two-phase first-order result, Theorem 2.4, rather than the Gaussian smoothing results in its later section. The source for every index here is the arXiv v1 preprint; the printed page and PDF page coincide.

Setting

The objective is a differentiable function $f:\mathbb R^n\to\mathbb R$ bounded below, with true optimal value $f^*=\inf_x f(x)$. Its gradient map $g=\nabla f$ is $L$-Lipschitz: $\|g(y)-g(x)\|\le L\|y-x\|$ for all $x,y$. The mission uses $L>0$. No convexity is assumed.

A stochastic first-order oracle returns a vector $G(x_k,\xi_k)$ when called at the current iterate $x_k$. Assumption A1 says that the oracle is unbiased at that iterate and has squared error of expectation at most $\sigma^2$. The noise may depend on earlier calls. In the formal model, the mean is conditional on the past, and the oracle and squared-error random variables are integrable. This makes explicit the probabilistic assumptions used when the paper takes expectations. The model takes $\sigma>0$ because the paper's displayed stepsize and post-optimization sample size otherwise degenerate under Lean's total division convention.

A randomized stochastic gradient run starts at $x_1$ and uses $x_{k+1}=x_k-\gamma_k G(x_k,\xi_k)$. It chooses an output index $R\in\{1,\ldots,N\}$ with the probability mass function in equation (2.3), independently of that run's oracle noise, and outputs $x_R$. For a parameter $\widetilde D>0$, the constant stepsize in equation (2.13) is $\gamma_k=\min\{1/L,\widetilde D/(\sigma\sqrt N)\}$. The paper writes $D_f=\sqrt{2(f(x_1)-f^*)/L}$ and

$$\mathcal B_N=\frac{LD_f^2}{N}+\left(\widetilde D+\frac{D_f^2}{\widetilde D}\right)\frac{\sigma}{\sqrt N}.$$

The two-phase method makes $S$ independent RSG runs, yielding $\bar x_1,\ldots,\bar x_S$. For each candidate it averages $T$ fresh oracle outputs into an estimated gradient $\widehat g_s$. It selects a candidate $\bar x^*$ whose estimated-gradient norm is minimal. The post-optimization samples for different candidates may be recycled; the paper explicitly permits this. The candidate itself is known when its post-optimization samples are taken, so the conditional oracle assumption is stated relative to a filtration that contains it at time zero.

Formalization targets

The first target is the one-run mean bound of Corollary 2.2, $L^{-1}\mathbb E\|g(x_R)\|^2\le\mathcal B_N$, followed by the one-run tail bound (2.19). Lemma 2.3(a) controls the squared norm of a sum of vector martingale differences. Equation (2.28) relates a selected candidate's true gradient to the best candidate and the estimation errors. Equation (2.29) factors the probability of simultaneous failure across independent runs and bounds it, while (2.30) bounds the other source of failure: every optimization run performing poorly, and at least one estimated gradient being inaccurate.

For arbitrary positive $S,N,T$, Theorem 2.4(a) states the stronger quantitative tail bound

$$\Pr\left\{\|g(\bar x^*)\|^2\ge2\left(4L\mathcal B_N+\frac{3\lambda\sigma^2}{T}\right)\right\}\le\frac{S+1}{\lambda}+2^{-S}\qquad(\lambda>0).$$

The goal, Theorem 2.4(b), specializes those parameters to an accuracy $\varepsilon>0$ and failure probability $0<\Lambda<1$. It sets

$$S=\left\lceil\log_2\frac2\Lambda\right\rceil,\qquad N=\left\lceil\max\left\{\frac{32L^2D_f^2}{\varepsilon},\left(\frac{32L(\widetilde D+D_f^2/\widetilde D)\sigma}{\varepsilon}\right)^2\right\}\right\rceil,\qquad T=\left\lceil\frac{24(S+1)\sigma^2}{\Lambda\varepsilon}\right\rceil.$$

An $(\varepsilon,\Lambda)$-solution is a point satisfying $\Pr\{\|g(\bar x^*)\|^2\le\varepsilon\}\ge1-\Lambda$. With those parameter values, the method computes such a point using at most $S(N+T)$ stochastic first-order oracle calls. The formulas and constants above are those of equations (2.24)–(2.27) of the preprint.

Significance

The result turns a mean stationarity bound into a confidence guarantee for one returned point, with an explicit accounting of both candidate generation and candidate assessment. It identifies how the number of independent runs, iterations per run, and gradient estimates scale with the requested accuracy and failure probability. Without the post-optimization phase, the paper's direct Markov estimate (2.19) gives a less favorable dependence on the failure probability for a single run.

The mathematical theorem is proved in the 2013 paper; the Lean declarations in this proposal are statements awaiting machine-checked proofs. A completed formalization would give reusable definitions for a dependent-noise gradient oracle, a random-output stochastic-gradient run, and a post-optimization selection procedure. The vector martingale estimate and finite-candidate inequality could also serve later stochastic-optimization developments. Contributions to the probability and measurability infrastructure, as well as the final rate proof, are within scope.

Difficulty

Selecting the candidate with the smallest estimated gradient need not select the one with the smallest true gradient. The estimates are random and are evaluated at candidates that are themselves random outputs of earlier oracle runs. A bound for a fixed point does not address this dependence. The argument needs a conditional oracle contract at each candidate and a separate guarantee that the optimization runs supply a sufficiently good candidate. The paper's intermediate display (2.31) applies a fixed-candidate tail estimate to the selected candidate even though the selection examines the same estimates; that step is not justified by the stated assumptions. The final bound can be phrased without treating (2.31) as a separate target.

Formalization scope

The space $\mathbb R^n$ is EuclideanSpace ℝ (Fin n). Iterations and oracle samples retain the paper's one-based indexing; candidate runs use Fin S, which reindexes the paper's $1,\ldots,S$. The objective's true infimum is carried by IsGLB (Set.range f) fstar, so $D_f$ cannot be weakened by substituting an arbitrary lower bound. Each candidate is generated by the RSG recursion and the mass function (2.3); the one-run performance bound is a theorem, not a hypothesis of the two-phase model. Independent candidate outputs are required for (2.29). The selected candidate may be any minimizer of the estimated-gradient norm, with any tie rule; no measurability of the selection is assumed.

All expectations are Bochner integrals, with integrability included in the oracle assumptions and in the one-run conclusion. Probabilities are ℝ≥0∞ values. The noise carrier is generalized from the paper's Borel subsets of Euclidean space to a measurable type. The general vector martingale lemma omits the paper's Polish-space and trivial-time-zero conditions; neither is needed for its squared-norm bound, and the latter would conflict with a candidate known at time zero. Its non-strict tail statement is restricted to a positive sum of variance bounds, because the printed version is false when every bound is zero. The goal reads the logarithm in (2.24) as base two, matching the paper's final $2^{-S}\le\Lambda/2$ comparison. Positive $\sigma$ and $\widetilde D$ avoid zero steps and an empty post-optimization average.

The formal goal states the probability guarantee. The $S(N+T)$ oracle-call bound follows directly from the algorithm's two loops and is recorded in the theorem's prose rather than in a separate operational semantics for calls. A complete proof may use Mathlib's conditional expectations, independent random variables, finite sums, and real logarithm and ceiling libraries. No assumption may simply postulate the one-run gradient bound or the final success probability.

Selected references

• Saeed Ghadimi and Guanghui Lan, Stochastic First- and Zeroth-Order Methods for Nonconvex Stochastic Programming, SIAM Journal on Optimization 23(4), 2013. arXiv:1309.5549v1.

Motivation

Stochastic programs and chance-constrained models in operations research contain expectations of the form

$$Z_0(x)=\mathbb E_0[h(x;\xi)],$$

where $\xi$ is a random input under an unknown distribution $P_0$, $x$ is a decision and $h$ is a known loss. In practice $P_0$ is seen only through i.i.d. data $\xi_1,\dots,\xi_n$. Distributionally robust optimization (DRO) replaces the unknown expectation by the worst case over a set of plausible distributions, and a recurring question is how large that set should be so that the robust value is a statistically valid bound, without being needlessly conservative.

The usual answer, for divergence-based sets on a finite support of size $k$, sets the radius from a $\chi^2_{k-1}$ quantile (Ben-Tal et al. 2013). This guarantee holds simultaneously for every function of the distribution, and so it is loose for a single expectation, and it degrades as $k$ grows. H. Lam (arXiv:1605.09349, Operations Research 2019) observed that if the ball is placed around the empirical distribution and measured with the Burg entropy, the right radius uses one degree of freedom, whatever the support. The DRO bounds then form an asymptotically exact confidence interval for $Z_0(x)$. The mechanism is the empirical likelihood method of Owen (Owen 1988; Owen 2001).

Timeline:

• 1988: Owen proves the empirical likelihood theorem for a mean (−2 log of the profile likelihood ratio tends to $\chi^2_1$), the nonparametric analogue of Wilks' theorem.
• 2013: Ben-Tal, den Hertog, De Waegenaere, Melenberg and Rennen calibrate $\varphi$-divergence balls on a finite support with $\chi^2_{k-1}$ quantiles.
• 2016: Lam (this paper) connects the empirical Burg-ball DRO to empirical likelihood, with pointwise (Theorem 2) and uniform (Theorem 4) coverage.
• 2016–2021: Duchi, Glynn and Namkoong (arXiv:1610.03425) develop the general $f$-divergence version and exact coverage of optimal values.

Setting

Data: $\xi_1,\xi_2,\dots$ are i.i.d. random elements of a measurable space $\Xi$ under $P_0$. A decision $x\in\mathbb R^m$ is fixed, and $h(x;\cdot):\Xi\to\mathbb R$ is measurable with $0<\operatorname{Var}_0 h(x;\xi)<\infty$.

The empirical Burg-entropy divergence ball of radius $\eta$ is the set of weight vectors on the $n$ data points

$$\mathcal U_n(\eta)=\Big\{w\in\mathbb R^n:\ -\frac1n\sum_{i=1}^n\log(nw_i)\le\eta,\ \sum_{i=1}^n w_i=1,\ w_i\ge0\Big\},$$

with $\log 0=-\infty$. It is the $\varphi$-divergence ball around the uniform weights $1/n$ for the Burg generator $\varphi(t)=-\log t+t-1$. Let $\chi^2_{1,1-\alpha}$ be the $(1-\alpha)$-quantile of the $\chi^2_1$ law, i.e. $P(Y^2\le\chi^2_{1,1-\alpha})=1-\alpha$ for $Y\sim N(0,1)$. The empirical DRO bounds are

$$\underline Z_n(x)=\min_{w\in\mathcal U_n(\chi^2_{1,1-\alpha}/(2n))}\sum_{i=1}^n h(x;\xi_i)w_i,\qquad \overline Z_n(x)=\max_{w\in\mathcal U_n(\chi^2_{1,1-\alpha}/(2n))}\sum_{i=1}^n h(x;\xi_i)w_i .$$

For a sample $z\in\mathbb R^n$ and a value $\mu$, the profile empirical likelihood statistic is

$$-2\log R(\mu)=\min\Big\{-2\sum_{i=1}^n\log(nw_i):\ \sum_i w_iz_i=\mu,\ \sum_i w_i=1,\ w_i\ge0\Big\},$$

equal to $+\infty$ when no weight vector is feasible.

Formalization targets

Goal: Theorem 2 (p. 11)

$$\lim_{n\to\infty}P\big(\underline Z_n(x)\le Z_0(x)\le\overline Z_n(x)\big)=1-\alpha .$$

The event is two-sided and joint. No rate is claimed.

Milestones

• Proposition 2 (p. 23). For every data vector and $\kappa\ge0$,

$$\min_{\mathcal U_n(\kappa/(2n))}\textstyle\sum_i w_iz_i\le\mu\le\max_{\mathcal U_n(\kappa/(2n))}\sum_i w_iz_i\iff -2\log R(\mu)\le\kappa .$$

• Theorem 1, the empirical likelihood theorem (p. 10). $-2\log R(Z_0(x))\Rightarrow\chi^2_1$.
• The pointwise steps of the paper's proof of Theorem 1 (the proof of Theorem 5 with a one-point index set): Lemma 3 (p. 38), $\max_{i\le n}|Y_i|=o(n^{1/2})$ a.s. for square-integrable i.i.d. $Y_i$; the closed form (49)–(50) (pp. 25–26) of the optimal weights; the multiplier bound (52) (p. 27); the expansion (68) (p. 30), $-2\log R=n s^{-1}\bar h^2+o_p(1)$.

Significance

Theorem 2 says that a robust optimization formulation with a data-driven Burg ball of radius $\chi^2_{1,1-\alpha}/(2n)$ yields an asymptotically exact $(1-\alpha)$ confidence interval for an expected loss, independently of the support size of $\xi$. Proposition 2 is the bridge: it identifies the DRO interval with the empirical likelihood confidence region, so any distributional result on $-2\log R$ transfers to the DRO bounds. This is the pointwise basis of the paper's later results on feasibility of chance and expectation constraints.

The result is proved in the paper, and Theorem 1 is classical. To our knowledge none of these statements has a machine-checked proof. Formalizing them requires the empirical likelihood theorem in Lean, which is reusable well beyond DRO (empirical likelihood tests and intervals for means). The Open platform theorem GenEmpLik.Coverage.theorem_3 (Duchi, Glynn and Namkoong) implies the goal when its decision set is the single point $\{x\}$ and its generator is $2\varphi$ with radius $\chi^2_{1,1-\alpha}$; it is a different, more general statement and is not proved. Likewise the Open GenEmpLik.Expansion.lemma_7 implies Lemma 3, and the Proved NumStochOpt.ListScheduling.lemma_8_1_i_max_over_sqrt_ae gives Lemma 3 after passing to $|Y_i|$.

Difficulty

Proposition 2 is a finite-dimensional convexity statement, but it needs attainment of the minimum defining $-2\log R$, which holds only because the objective blows up at the boundary of the simplex. The central difficulty is Theorem 1. The optimal weights are explicit only through a Lagrange multiplier $\lambda$ that solves a nonlinear equation; to expand $-2\log R$ one must first show $\lambda=O_p(n^{-1/2})$, which requires controlling the largest observation by Lemma 3 together with the law of large numbers for $s$, and then bound the cubic Taylor remainder of $\log(1+\lambda\tilde h_i)$ uniformly in $i$. A direct appeal to the central limit theorem for the sample mean does not suffice, since $-2\log R$ is an implicitly defined optimal value, not a smooth function of the mean.

Formalization scope

• The sample is a sequence ξ : ℕ → Ω → Ξ, measurable, mutually independent and identically distributed, indexed from $0$; the size-$n$ sample is fun i : Fin n => h x (ξ i ω).
• The moment assumption is MemLp (h x ∘ ξ 0) 2 together with 0 < variance. Positive variance cannot be dropped: for a constant loss the coverage is $1$.
• The quantile is pinned by gaussianReal 0 1 {y | y^2 ≤ q} = ofReal (1 - α) with $0<\alpha<1$; it is not a free parameter, and it is not the normal quantile.
• The upper bound is the published GenEmpLik.Expansion.robustMean burg (q/2) (its ball has radius $\rho/n$, so $\rho=q/2$); the lower bound is its infimum twin robustLower. The ball uses the published probUncertaintySet and phiDiv.
• The Burg generator is $+\infty$ at $t\le0$, and the statistic elStat minimizes over strictly positive weights with the empty infimum equal to $+\infty$. Lean's junk value $\log 0=0$ therefore never enters: a formalization with $w_i\ge0$ and real logarithms would admit boundary weights with spurious finite divergence and is ruled out. For $n\ge1$ the value sets of the DRO bounds are nonempty and bounded, so the real infimum and supremum are not junk values.
• The decision set $\Theta$ is dropped (the statements are pointwise in $x$). Theorem 1 is stated as convergence of distribution functions at every $t$, which is equivalent to $\Rightarrow\chi^2_1$ since the limit is continuous. Proposition 2 is stated for general $\kappa\ge0$ and every data vector, as displayed on p. 11. Milestones (49)–(50), (52), (68) are the one-point ($\Theta=\{x\}$) specialisations of the proof of Theorem 5.
• Probabilities of events are outer probabilities; no measurability of events is assumed, as the paper ignores these issues (p. 14).
• Not posed: the headline Theorem 4 (uniform calibration over $\Theta$ by the quantile of a $\chi^2$-process supremum), Theorem 5 and Lemma 2, which need Donsker classes and weak convergence in $\ell^\infty(\Theta)$, for which Mathlib has no substrate.

Contributions welcome: proofs of Proposition 2 and Lemma 3 (from the Proved reference), the KKT analysis (49)–(52), and a Lean empirical likelihood theorem.

Selected references

• H. Lam, Recovering Best Statistical Guarantees via the Empirical Divergence-Based Distributionally Robust Optimization, arXiv:1605.09349v1, 2016; Operations Research 67(4), 2019. https://arxiv.org/abs/1605.09349v1
• A. B. Owen, Empirical likelihood ratio confidence intervals for a single functional, Biometrika 75(2), 1988. https://doi.org/10.1093/biomet/75.2.237
• A. B. Owen, Empirical Likelihood, Chapman & Hall/CRC, 2001. https://doi.org/10.1201/9781420036152
• A. Ben-Tal, D. den Hertog, A. De Waegenaere, B. Melenberg, G. Rennen, Robust solutions of optimization problems affected by uncertain probabilities, Management Science 59(2), 2013. https://doi.org/10.1287/mnsc.1120.1641
• J. C. Duchi, P. W. Glynn, H. Namkoong, Statistics of robust optimization: a generalized empirical likelihood approach, Mathematics of Operations Research 46(3), 2021. https://arxiv.org/abs/1610.03425

Motivation

Least-squares regression, $\min_X\|AX-B\|_F$ for a tall matrix $A\in\mathbb R^{n\times d}$ with $n\gg d$, is among the most frequently solved problems in numerical computation. Direct methods (QR or SVD) cost $O(nd^2)$ time, which is prohibitive when $n$ is in the millions. Sketching replaces the problem by a much smaller one: choose a matrix $S\in\mathbb R^{t\times n}$ with $t\ll n$, form $SA$ and $SB$, and solve $\min_X\|S(AX-B)\|_F$ instead. The question is when the solution of the small problem is nearly optimal for the original one.

Clarkson and Woodruff, Low Rank Approximation and Regression in Input Sparsity Time (STOC 2013; JACM 2017), showed that a very sparse random sketch (one non-zero per column, the CountSketch matrix) can be applied in time proportional to the number of non-zero entries of $A$ and still yields $(1+\varepsilon)$-approximate regression and low-rank approximation. §7 of that paper isolates the deterministic part of the argument: two conditions on a fixed sketch $S$ that guarantee a near-optimal solution, whatever distribution produced $S$.

Timeline. Sarlós (FOCS 2006) introduced sketch-and-solve regression with fast Johnson–Lindenstrauss transforms. Drineas, Mahoney, Muthukrishnan and Sarlós (Numer. Math. 2011) stated the structural conditions for single-response regression (their Lemma 1); Theorem 36 of Clarkson–Woodruff generalizes it to multiple responses. Clarkson and Woodruff then supplied the sparse sketches meeting the conditions in input-sparsity time.

Setting

All matrices are real. For a vector $v$ let $\|v\|_2^2=\sum_iv_i^2$; for a matrix $M$ let $\|M\|_F^2=\sum_{i,j}M_{ij}^2$. Fix $A\in\mathbb R^{n\times d}$ of rank at most $r$, a right-hand side $B\in\mathbb R^{n\times d'}$ and a sketch $S\in\mathbb R^{t\times n}$. Let $C(A)=\{Ax\}$ be the column space and $U\in\mathbb R^{n\times k}$ a matrix with orthonormal columns ($U^\top U=I$) spanning $C(A)$.

• $Y^*$ is a least-squares minimizer: it minimizes $\|AY-B\|_F^2$ (problem (8)).
• $\tilde Y$ is a sketched minimizer: it minimizes $\|S(AY-B)\|_F^2$ (problem (7)).
• $S$ is a subspace embedding for $A$ with parameter $\eta$ if $\bigl|\|SAx\|_2^2-\|Ax\|_2^2\bigr|\le\eta\|Ax\|_2^2$ for every $x$.
• $S$ satisfies the approximate matrix multiplication event (Lemma 32) for a pair $(U,W)$ with error parameter $e$ if $\|U^\top S^\top SW-U^\top W\|_F^2\le e^2\|U\|_F^2\|W\|_F^2$.

Formalization targets

Goal: Theorem 36 (p. 24)

If $S$ is a subspace embedding for $A$ with $0\le\eta\le1/2$ and satisfies the event of Lemma 32 for $(U,B-AY^*)$ with error parameter $\sqrt{\varepsilon/r}$, then every sketched minimizer satisfies

$$\|A\tilde Y-B\|_F^2\le(1+4\varepsilon)\|AY^*-B\|_F^2 .$$

This is the explicit form the paper's proof establishes; the printed $(1+\epsilon)$ bound on norms follows after rescaling $\varepsilon$.

Milestones

• Fact 34 (Pythagorean theorem): $C^\top D=0\Rightarrow\|C+D\|_F^2=\|C\|_F^2+\|D\|_F^2$.
• Fact 35 (normal equations): $C^-D$ is a minimizer; every minimizer $X^*$ satisfies $C^\top(CX^*-D)=0$ and $\|CX-D\|_F^2=\|C(X-X^*)\|_F^2+\|CX^*-D\|_F^2$.
• Lemma 37: for $A$ with orthonormal columns, $\|A(\tilde Y-Y^*)\|_F\le2\sqrt\varepsilon\|B-AY^*\|_F$.
• Theorem 39 (affine embeddings, an independent companion result of §7.5): under an embedding with parameter $\varepsilon\le1/2$ and the event of Lemma 32 with parameter $\varepsilon/\sqrt r$, $\|S(AX-B)\|_F^2-\|S\tilde B\|_F^2$ equals $\|AX-B\|_F^2-\|\tilde B\|_F^2$ up to $\pm2\varepsilon\|AX-B\|_F^2$ for all $X$, and $\|S(AX-B)\|_F^2=(1\pm3\varepsilon)\|AX-B\|_F^2$ if $S$ also preserves $\|\tilde B\|_F^2$ to $1\pm\varepsilon$.

Significance

Theorem 36 separates the algorithm design for sketched regression from its correctness proof: any random construction (CountSketch, subsampled randomized Hadamard transforms, leverage-score sampling, Gaussian matrices) yields a $(1+O(\varepsilon))$-approximate regression solver once it is shown to satisfy the two conditions with constant probability. Clarkson and Woodruff use it for their input-sparsity-time regression algorithm (Theorem 38) and, through Theorem 39 and affine embeddings, for low-rank approximation. The same conditions recur throughout the randomized numerical linear algebra literature (see Woodruff's survey Sketching as a Tool for Numerical Linear Algebra, 2014).

The results are proved in the paper (§7.3, §7.5 and Appendix A, p. 37). To our knowledge none has a machine-checked proof. Formalizing them produces a reusable, distribution-free layer of sketching theory: the probabilistic statements about specific sketches can then be proved separately and composed with it.

Difficulty

The obvious argument applies the subspace embedding to the residual $A\tilde Y-B$, but that vector does not lie in $C(A)$, so an embedding of $C(A)$ alone says nothing about it; embedding $C(A)$ augmented with the columns of $B$ instead costs a sketch dimension that grows with $d'$. The hypotheses of Theorem 36 control the optimal residual $B-AY^*$ only through the approximate-product condition, which is far weaker than norm preservation, and the difficulty is to transfer that control to the sketched solution (Lemma 37) with the constants $2\sqrt\varepsilon$ and $1+4\varepsilon$. A second obstacle is the passage from a general $A$ of rank at most $r$ to an orthonormal basis $U$ of $C(A)$: it changes $\|A\|_F$, which is why the approximate-product condition is stated for $U$.

Formalization scope

Matrices are Matrix (Fin m) (Fin n) ℝ; squared norms are explicit sums sqNorm and frobSq; the column space is the range of Matrix.mulVecLin; rank is Matrix.rank. Everything is deterministic: $S$ is an arbitrary fixed matrix and no probability space appears. The committed conventions are:

1. The subspace embedding is in squared form. In Theorem 36 and Lemma 37 its parameter $\eta$ plays the paper's $\epsilon_0^2$ (so $\epsilon_0\le1/\sqrt2$ becomes $\eta\le1/2$), which is how the proof of Lemma 37 (p. 37) uses it; in Theorem 39 the parameter is $\varepsilon$ itself, as that proof uses it.
2. "$S$ satisfies Lemma 32" is the event inside the probability, for the pair $(U,B-AY^*)$ (resp. $(U,\tilde B)$) that the proofs use, with $\le$ in place of the printed $<$ (a weaker hypothesis).
3. $\tilde Y$, $Y^*$, $X^*$ are arbitrary minimizers; "the solution" need not be unique. The Moore–Penrose inverse is specified by the four Penrose equations.
4. $r\ge1$ guards $\varepsilon/r$ and $\varepsilon/\sqrt r$; "$n>d$" (§7) concerns running times and is not assumed.

A trivializing formalization is ruled out: $\tilde Y$ is never defined as $Y^*$, the approximate-product hypothesis is an inequality and not the exact identity $U^\top S^\top S(B-AY^*)=0$, $\eta$ is non-negative, and the goal is stated for general $A$ of rank at most $r$, not only for orthonormal $A$. The hypotheses are jointly satisfiable (for instance $S=I$).

A complete development needs the Frobenius inner product and its Pythagorean identity, the normal equations for matrix least squares, the existence of the Moore–Penrose inverse, and the spectral-norm bound for a symmetric matrix from its quadratic form; all are reusable beyond this mission. Proofs of any milestone, and alternative arguments for Theorem 36, are welcome.

Selected references

• K. L. Clarkson and D. P. Woodruff, Low Rank Approximation and Regression in Input Sparsity Time, arXiv:1207.6365v4, 2013; J. ACM 64(6), 2017. https://arxiv.org/abs/1207.6365
• P. Drineas, M. W. Mahoney, S. Muthukrishnan and T. Sarlós, Faster Least Squares Approximation, Numerische Mathematik 117(2), 2011. https://arxiv.org/abs/0710.1435
• T. Sarlós, Improved Approximation Algorithms for Large Matrices via Random Projections, FOCS 2006. https://doi.org/10.1109/FOCS.2006.37
• D. P. Woodruff, Sketching as a Tool for Numerical Linear Algebra, Foundations and Trends in Theoretical Computer Science 10(1–2), 2014. https://arxiv.org/abs/1411.4357

Motivation

Large optimization problems often split into blocks owned by different agents (subsystems of a power grid, vehicles in a formation, scenarios of a stochastic program) that are coupled only through a few linear constraints. Distributed optimization methods solve such problems by letting each agent work on its own block and exchanging a small amount of information, typically a vector of prices for the coupling constraints. The classical method of this kind is dual decomposition (Everett 1963; Bertsekas and Tsitsiklis, Parallel and Distributed Computation, 1989): fix the prices, let every agent solve its own subproblem, then update the prices by a gradient or Newton step on the concave dual function. Dual decomposition is reliable for strictly convex problems but has no convergence guarantee for nonconvex ones.

Houska, Frasch and Diehl (SIAM J. Optim. 26 (2016)) proposed ALADIN (augmented Lagrangian based alternating direction inexact Newton), which adds a proximal term to the agents' subproblems and replaces the price update by a coupled equality-constrained quadratic program. ALADIN is used in distributed model predictive control and in power-system optimization. Remark 4.2 and Appendix A of the paper explain where it sits relative to the classical method: if the proximal term is switched off ($\rho = 0$), the price update of ALADIN is exactly a regularized Newton step of dual decomposition. Lemma 5 is the precise statement of that equivalence, and it is the goal of this mission.

Setting

The problem (1.1) has $N$ blocks $x_i\in\mathbb R^n$:

$$\min_x\ \sum_{i=1}^N f_i(x_i)\quad\text{s.t.}\quad \sum_{i=1}^N A_i x_i = b,\qquad h_i(x_i)\le 0,$$

with $f_i:\mathbb R^n\to\mathbb R$, $h_i:\mathbb R^n\to\mathbb R^{n_h}$ twice continuously differentiable, $A_i\in\mathbb R^{m\times n}$ and $b\in\mathbb R^m$. One iteration of Algorithm 2 (ALADIN) starts from a primal point $x$ and prices $\lambda\in\mathbb R^m$ and does the following.

1. Each agent solves the decoupled problem $\min_{y_i} f_i(y_i) + \lambda^\top A_i y_i + \frac\rho2\|y_i-x_i\|^2_{\Sigma_i}$ subject to $h_i(y_i)\le 0$, returning a minimizer $y_i$ and a multiplier $\kappa_i\ge 0$. With $\rho = 0$ this is the dual-decomposition subproblem (A.2).
2. Each agent forms the active-constraint Jacobian $C^*_i$ (rows $\nabla(h_i)_j(y_i)^\top$ for active $j$, zero rows otherwise), chooses an approximation $C_i$ and a symmetric Hessian approximation $H_i$, and computes the modified gradient $g_i = \nabla f_i(y_i) + (C^*_i - C_i)^\top\kappa_i$.
3. For a parameter $\mu>0$, a coordinator solves the coupled QP (3.3)

$$\min_{\Delta y, s}\ \sum_i\Big\{\tfrac12\Delta y_i^\top H_i\Delta y_i + g_i^\top\Delta y_i\Big\} + \lambda^\top s + \tfrac\mu2\|s\|_2^2\ \ \text{s.t.}\ \ \sum_i A_i(y_i+\Delta y_i) = b + s\ \mid\ \lambda_{\mathrm{QP}},\quad C_i\Delta y_i = 0,$$

and the new prices are $\lambda^+ = \lambda + \alpha_3(\lambda_{\mathrm{QP}}-\lambda)$.

On the dual-decomposition side, $\nabla V(\lambda) = \sum_i A_i y_i - b$ (A.3) is the gradient of the dual function at $\lambda$, and the (inexact) dual Newton step (A.4) with a symmetric negative semidefinite scaling matrix $M$ and step size $\alpha\in(0,1]$ is

$$\lambda^+_{\mathrm{DD}} = \lambda - \alpha\Big(M - \frac1\mu I\Big)^{-1}\nabla V(\lambda),$$

where $1/\mu$ plays the role of a Levenberg–Marquardt regularization parameter.

Formalization targets

Goal: Lemma 5 (p. 1123)

Assume every $f_i$ strictly convex, every $h_i$ convex, $C_i$ of full row rank, $H_i$ positive definite, $\rho = 0$ and $\alpha_3 = \alpha$. Then for every $\mu>0$, with

$$M = -\sum_{i=1}^N A_i\Big[H_i^{-1} - H_i^{-1}C_i^\top\big[C_iH_i^{-1}C_i^\top\big]^{-1}C_iH_i^{-1}\Big]A_i^\top,$$

the matrix $M - \frac1\mu I$ is invertible and

$$\lambda^+ = \lambda^+_{\mathrm{DD}} .$$

Milestones (the steps of the proof, pp. 1123–1124)

1. Slack elimination: $\min_s \frac\mu2\|s\|_2^2 - (\lambda_{\mathrm{QP}}-\lambda)^\top s = -\frac1{2\mu}\|\lambda_{\mathrm{QP}}-\lambda\|_2^2$, attained only at $s = (\lambda_{\mathrm{QP}}-\lambda)/\mu$.
2. Stationarity: with $\rho = 0$, $0 = \nabla f_i(y_i) + A_i^\top\lambda + (C^*_i)^\top\kappa_i = g_i + A_i^\top\lambda + C_i^\top\kappa_i$.
3. Explicit minimization over $\Delta y$: for $H\succ 0$ and $C$ of full row rank, $-Pw$ with $P = H^{-1} - H^{-1}C^\top(CH^{-1}C^\top)^{-1}CH^{-1}$ minimizes $\frac12\Delta^\top H\Delta + w^\top\Delta$ on $\{C\Delta = 0\}$, with value $-\frac12 w^\top Pw\le 0$.
4. Argmax display: for $M$ symmetric negative semidefinite, $\lambda - (M - \frac1\mu I)^{-1}g$ is the unique maximizer of $\frac12(\nu-\lambda)^\top M(\nu-\lambda) + \nu^\top g - \frac1{2\mu}\|\nu-\lambda\|_2^2$.

Significance

The result. Lemma 5 places ALADIN on a known map. With $\rho=0$ it inherits the convergence theory of dual decomposition for strictly convex problems for small step sizes, independently of how crude the approximations $H_i$, $C_i$ are, as long as $H_i\succ 0$. It also reinterprets the augmented-Lagrangian parameter $\mu$ as the inverse of a dual Levenberg–Marquardt parameter, and it isolates the proximal term $\frac\rho2\|y_i - x_i\|^2_{\Sigma_i}$ as the only ingredient that distinguishes ALADIN from dual decomposition and makes the nonconvex analysis of the paper possible.

Formalizing it. The lemma is proved in the paper; to our knowledge no machine-checked proof exists. A formal proof certifies an identity whose printed form contains a typo (see below), and produces reusable linear-algebra facts: the closed-form solution of an equality-constrained strictly convex quadratic program via the reduced inverse $P$, the semidefiniteness of $P$, and the closed-form maximizer of a regularized concave quadratic.

Difficulty

The proof is a chain of exact eliminations, and the work is in carrying each one out without losing a sign or a hypothesis. The step most easily done wrong is the minimization over $\Delta y$: the constraint $C_i\Delta y_i = 0$ couples the block coordinates, the minimizer is expressed through the inverse of $C_iH_i^{-1}C_i^\top$ (invertible only because $C_i$ has full row rank). A second trap is the sign convention for $\lambda_{\mathrm{QP}}$: it must agree with footnote 4 ($\lambda + \mu s - \lambda_{\mathrm{QP}} = 0$), or the conclusion flips. The tempting shortcut of treating the max–min dual of (3.3) abstractly runs into the fact that strong duality and attainment are not free; working from the KKT system of (3.3) avoids that.

Formalization scope

Vectors are Fin n → ℝ; blocks are indexed by Fin N; matrices are Mathlib Matrix over Fin index types; the gradient is the vector of partial derivatives computed with fderiv. The step-1 output $(y_i,\kappa_i)$ is a local minimizer of the decoupled problem together with KKT multipliers (the paper assumes LICQ for the lower-level constraints). The QP (3.3) and $\lambda_{\mathrm{QP}}$ enter through the KKT system of (3.3), not through a closed-form expression: defining $\lambda_{\mathrm{QP}}$ or $\Delta y$ by the formula of the conclusion would make the goal true by definition and is ruled out. $\nabla V(\lambda)$ is the formula (A.3); that it is the gradient of the dual function is a cited fact the lemma does not use. The invertibility of $M - \frac1\mu I$ is part of the conclusion, because Mathlib's matrix inverse returns $0$ on a singular matrix.

Correction of (A.6). As printed, (A.6) has the bracket $[C_iH_iC_i^\top]^\dagger$. This is a typo for $[C_iH_i^{-1}C_i^\top]^\dagger$, as in (A.5) on the same page and as the explicit minimization over $\Delta y$ forces; with the printed bracket the identity fails already for $N = 1$, $H_1 = \mathrm{diag}(1,2)$, $C_1 = (1,1)$, $A_1 = I$. The mission uses the corrected matrix and writes the pseudo-inverse as an inverse, which it is under the hypotheses.

The hypotheses of strict convexity of $f_i$, convexity of $h_i$ and $\alpha\in(0,1]$ are kept as on the page although the identity does not use them. The related platform draft GoldfarbIdnani.DualQP.properties_2_6_to_2_9 (not proved) states properties of the same reduced inverse; it is not imported. Proofs of the four milestones and of the goal are welcome; the milestones are self-contained matrix and calculus facts and can be attacked independently.

Selected references

• B. Houska, J. Frasch, M. Diehl, An augmented Lagrangian based algorithm for distributed nonconvex optimization, SIAM J. Optim. 26(2) (2016), 1101–1127. https://doi.org/10.1137/140975991
• D. P. Bertsekas, J. N. Tsitsiklis, Parallel and Distributed Computation: Numerical Methods, Prentice-Hall, 1989.
• H. Everett, Generalized Lagrange multiplier method for solving problems of optimum allocation of resources, Oper. Res. 11(3) (1963), 399–417. https://doi.org/10.1287/opre.11.3.399
• D. P. Bertsekas, Nonlinear Programming, 2nd ed., Athena Scientific, 1999.

Motivation

Quasi-Newton methods minimize a smooth function $f:\mathbb R^n\to\mathbb R$ using gradients only, while building up an approximation to the inverse of the Hessian matrix from the changes of gradient they observe. They are the default algorithms of unconstrained nonlinear optimization: BFGS and its limited-memory variant L-BFGS sit inside most numerical optimization libraries and most large-scale fitting codes in statistics and machine learning.

The family starts with William C. Davidon's Argonne report ANL-5990 of 1959, Variable Metric Method for Minimization, rejected by a journal in 1957 and published only in 1991 as the first article of the SIAM Journal on Optimization, with a preface by the author (DOI 10.1137/0801001). The report's body describes a flowchart algorithm; its Appendix describes "a simplified method embodying some of the ideas" of that algorithm, in one and a half pages. That simplified method is what is now called the symmetric rank-one (SR1) update, and the Appendix contains its quadratic-termination property: on a quadratic, after at most $n$ steps the trial matrix equals the inverse Hessian, and the next step lands at the minimum.

Timeline.

• 1952: Hestenes and Stiefel's conjugate gradient method minimizes a strictly convex quadratic in at most $n$ steps (J. Res. NBS 49); Davidon cites it as [2].
• 1959: Davidon's ANL-5990, including the Appendix method treated here.
• 1963: Fletcher and Powell simplify the body's method into the DFP update and prove its quadratic termination with exact line searches (Comput. J. 6).
• 1967: Broyden's survey of quasi-Newton methods discusses the symmetric rank-one update (Math. Comp. 21); textbooks later analyse it under the name SR1.
• 1980: Nocedal's limited-storage BFGS, whose termination on quadratics with exact line searches is a separate Prove2Me mission (Math. Comp. 35).

Setting

Vectors are elements of $\mathbb R^n$ and matrices are real $n\times n$. Fix a symmetric positive definite matrix $G$ and a point $\xi$, and let

$$f(x)=\tfrac12\,(x-\xi)^{\mathsf T}G\,(x-\xi)+c ,$$

a quadratic with constant Hessian $G$ and minimum point $\xi$. Its gradient is $\nabla(x)=G(x-\xi)$, display (6) of the Appendix.

The method keeps a point $x_k$ and a symmetric trial matrix $H_k$, which plays the role of $G^{-1}$. Write $\nabla_k=\nabla(x_k)$. One iteration is

$$x_{k+1}=x_k-H_k\nabla_k \quad\text{(2)},\qquad H_{k+1}=H_k+a_k\,(H_k\nabla_{k+1})(H_k\nabla_{k+1})^{\mathsf T}\quad\text{(3)},$$

a full step with no line search followed by a symmetric rank-one update. With the two scalars

$$N_k=\nabla_{k+1}^{\mathsf T}H_k\nabla_{k+1},\qquad M_k=\nabla_{k+1}^{\mathsf T}H_k\nabla_k ,$$

footnote 2 of the Appendix fixes the coefficient on a quadratic to $a_k=(M_k-N_k)^{-1}$, which makes the paper's quality measure $\Delta$ of (4) vanish. Write $S_k=x_{k+1}-x_k$ for the step and $D_k=\nabla_{k+1}-\nabla_k$ for the change of gradient. The paper writes $x$, $H$, $\nabla$ for the current quantities and $x^{+}$, $H^{+}$, $\nabla^{+}$ for the updated ones, and uses the letter $N$ both for the number of variables and for the scalar $N_k$; here the number of variables is $n$.

Formalization targets

Goal: termination in $n$ steps

If $H_0$ is symmetric and $M_k\ne N_k$ for every $k<n$, then

$$H_n=G^{-1}\qquad\text{and}\qquad x_{n+1}=\xi .$$

This is the sentence after (8) on p. 17: "After no more than $N$ steps (for which $\Delta=0$), $H$ will equal $G^{-1}$ and the following step will be to the exact minimum."

Milestones

1. Display (6): the gradient of $f$ is $G(x-\xi)$.
2. Display (8): with $a=(M-N)^{-1}$ and $M\ne N$, $H^{+}GS=H^{+}D=S$.
3. Display (7): if $HGu=u$ then $H^{+}Gu=u$, for any coefficient $a$.
4. The sentence "so that $S$ becomes another such eigenvector", read along the run: if $M_i\ne N_i$ for $i<k$, then $H_kGS_j=S_j$ for every $j<k$.

Further statements of the Appendix

• Condition 1 (p. 16): if $H$ is positive definite and $R^{-1}\le\det H^{+}/\det H\le R$ with $R>1$, then $H^{+}$ is positive definite.
• The constraint remark (p. 17): for any function and any coefficients, if $H_0$ is symmetric and $H_0b=0$, then every step is perpendicular to $b$ and $b\cdot x_k$ is conserved.
• Table (5) (p. 16): for general (non-quadratic) $f$, the coefficient $a$ that minimizes $\Delta$ subject to condition 1, and the minimum value, on each of five ranges of $M$.

Significance

The termination theorem says that on a quadratic the rank-one update recovers the exact inverse Hessian from $n$ gradient differences, without line searches and without positive definiteness of the trial matrix, and then takes the Newton step. This is the property that justifies the name "variable metric": the metric $H_k$ converges to the one that makes the problem trivial. The secant equation (8) and the persistence of earlier secant equations (7) are the template for the analysis of every later quasi-Newton update, and the SR1 update remains in use in trust-region methods because it does not force positive definiteness.

The result is classical and its proof is short; what this mission adds is a machine-checked version stated in the paper's own form. In particular the hypothesis is the paper's: only that each of the first $n$ coefficients is defined. Textbook statements of SR1 termination often assume in addition that the steps are linearly independent; the paper does not, and the formal goal does not either. No Lean formalization of SR1, DFP or BFGS termination was found in Mathlib or on Prove2Me. On Prove2Me, the related open targets are L-BFGS termination (Nocedal 1980) and conjugate-gradient termination (Hestenes–Stiefel 1952); the definitions here are independent of theirs.

Difficulty

Each update changes $H$ by a rank-one term, and nothing in a single step forces $H_n$ to be $G^{-1}$: the coefficient $a_k$ depends on the current iterate and the update could, a priori, destroy what earlier steps achieved, as general rank-one updates do. The obvious count — $n$ steps, $n$ conditions — does not by itself give $n$ independent conditions on $H_n$; the steps produced by the method are not chosen to be independent, and the only hypothesis available is that the denominators $M_k-N_k$ are nonzero. The difficulty is to get the full conclusion from that hypothesis alone, without assuming independence, positive definiteness of $H_0$, or a line search.

Formalization scope

Vectors are Fin n → ℝ with 0-based indices, matrices Matrix (Fin n) (Fin n) ℝ; $u^{\mathsf T}Hv$ is u ⬝ᵥ (H *ᵥ v) and $(Hg)(Hg)^{\mathsf T}$ is vecMulVec (H *ᵥ g) (H *ᵥ g). The definition file DavidonVM.Termination.Method provides quad, grad, Nval, Mval, the run predicate IsRunWith for an arbitrary gradient field and coefficients, and IsQuadRun for the quadratic with $a_k=(M_k-N_k)^{-1}$; DavidonVM.Termination.Delta provides the update and the quantity $\Delta$ of (4).

Committed conventions, each a reading of something the paper leaves implicit:

• The function is an explicit quadratic with constant symmetric positive definite $G$ ("in the neighborhood of a minimum"); the run's gradient field is $G(x-\xi)$, and milestone 1 ties it to the quadratic. The goal uses positive definiteness; milestones use only symmetry where that suffices; the secant milestone needs nothing on $G$.
• Only $H_0$ is assumed symmetric (the paper, p. 4: "It is to be symmetric"); positive definiteness of $H_0$ is not assumed.
• "Steps for which $\Delta=0$" is read through footnote 2 as $M_k\ne N_k$. Lean's real inverse gives $0^{-1}=0$, so without this hypothesis the update would silently leave $H$ unchanged.
• The stopping rule "provided that $N$ is greater than some preassigned $\varepsilon$" is not modelled; the run continues for every $k$, and the goal looks at steps $0,\dots,n$ only.
• $G^{-1}$ is Mathlib's matrix inverse, the true inverse since $G$ is positive definite.

A statement of the goal with the extra hypothesis "$S_0,\dots,S_{n-1}$ are linearly independent", or with $H_k G S_j=S_j$ assumed, would be a weaker theorem than the paper's and is ruled out. The case $n=0$ is trivial; the hypotheses are jointly satisfiable for $n=1$ (checked locally).

Welcome contributions: proofs of the four milestones and the goal; general facts about rank-one updates (determinant and positive definiteness, Sherman–Morrison form of the inverse) that the further statements need and that are reusable for DFP, BFGS and trust-region analyses.

Selected references

• W. C. Davidon, Variable Metric Method for Minimization, SIAM J. Optim. 1(1) (1991) 1–17 (Argonne report ANL-5990, 1959). https://doi.org/10.1137/0801001
• M. R. Hestenes and E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Nat. Bur. Standards 49 (1952) 409–436. https://doi.org/10.6028/jres.049.044
• R. Fletcher and M. J. D. Powell, A rapidly convergent descent method for minimization, Comput. J. 6 (1963) 163–168. https://doi.org/10.1093/comjnl/6.2.163
• C. G. Broyden, Quasi-Newton methods and their application to function minimisation, Math. Comp. 21 (1967) 368–381. https://doi.org/10.1090/S0025-5718-1967-0224273-2
• J. Nocedal, Updating quasi-Newton matrices with limited storage, Math. Comp. 35 (1980) 773–782. https://doi.org/10.1090/S0025-5718-1980-0572855-7

Why lower bounds matter for progressive hedging

Progressive hedging separates a stochastic optimization problem by scenario, solves those smaller problems, and coordinates their first-stage decisions through weighted aggregation and price updates. It is useful when the full extensive formulation is too large to solve directly. For stochastic mixed-integer programs, its iterates can provide feasible decisions, but a feasible cost alone does not reveal how far that cost is from the global optimum. A lower bound gives the missing comparison: the difference between a feasible cost and a valid lower bound is an optimality gap. Branch-and-bound methods also use lower bounds to discard regions whose best possible cost cannot improve the current solution. These are the reasons the authors seek a bound available during a progressive-hedging run, rather than only after convergence. Gade et al., author manuscript, §§1–3.

The paper establishes lower bounds for two-stage stochastic mixed-integer programs, extends them to scenario bundles, and studies their behavior computationally. The mission concerns the mathematical statement that the price vector produced at each completed iteration is eligible for a decomposed lower bound. The paper was published in Mathematical Programming in 2016; the theorem and page labels here follow the authors' SAND2013-9195J manuscript. Gade et al., author manuscript.

Finite-scenario two-stage setting

Let $\Xi$ be a finite set of scenarios and let $p_\xi>0$ be their probabilities, with $\sum_{\xi\in\Xi}p_\xi=1$. A first-stage decision $x\in\mathbb R^{n_1}$ is made before the scenario is known. Its first $p_1$ coordinates are nonnegative integers; the others are unrestricted reals. It incurs cost $c^\top x$ and must satisfy $Ax\ge b$. After scenario $\xi$ is observed, a second-stage decision $y(\xi)\in\mathbb R^{n_2}$ incurs cost $g(\xi)^\top y(\xi)$. Its first $p_2$ coordinates are nonnegative integers, and it must satisfy $Wy(\xi)\ge r(\xi)-T(\xi)x$. The recourse matrix $W$ is common to all scenarios. The scenario feasible set $X(\xi)$ contains pairs $(x,y)$ satisfying these restrictions for one $\xi$. Gade et al., §2, pp. 3–5.

In the extensive scenario formulation, each scenario has a copy $x(\xi)$ of the first-stage vector, and a common vector $\hat x$ enforces implementability through $p_\xi x(\xi)-p_\xi\hat x=0$. Its optimum is $z^*$. Algorithm 1 starts with zero prices $w^0(\xi)=0$, solves unpenalized scenario problems, and at iteration $\nu$ computes $\hat x^\nu=\sum_\xi p_\xi x^\nu(\xi)$ and $w^\nu(\xi)=w^{\nu-1}(\xi)+\rho(x^\nu(\xi)-\hat x^\nu)$. Later scenario solves add both $w^\nu(\xi)^\top x$ and $(\rho/2)\|x-\hat x^\nu\|_2^2$ to their objectives. Here $\rho$ is the algorithm's scalar parameter. Gade et al., Algorithm 1, p. 6.

Formalization targets

Scenario price bound

For any scenario prices $w(\xi)$, define the decomposed dual value by

$$D_\xi(w(\xi))=\inf_{(x,y)\in X(\xi)}\left(c^\top x+g(\xi)^\top y+w(\xi)^\top x\right),\qquad D(w)=\sum_{\xi\in\Xi}p_\xi D_\xi(w(\xi)).$$

Proposition 1 states that prices with $\sum_\xi p_\xi w(\xi)=0$ satisfy $D(w)\le z^*$. The price invariant in §3.1 states that the update maintains this weighted zero sum. The mission's goal combines these source claims at every completed iteration $N\ge1$ of Algorithm 1:

$$\sum_{\xi\in\Xi}p_\xi w^N(\xi)=0,\qquad D(w^N)\le z^*.$$

The formal goal is the lower-bound inequality; the zero-sum statement and Proposition 1 are its two milestones. Gade et al., Proposition 1, p. 6, and §3.1, p. 7.

Bundle extension

For a partition $B$ of the scenarios, a bundle $\beta$ has probability $P_\beta=\sum_{\xi\in\beta}p_\xi$. A bundle subproblem uses one first-stage vector for its scenarios and second-stage weights $p_\xi/P_\beta$. Proposition 3 gives the corresponding inequality $D_B(w)\le z^*$ when $\sum_{\beta\in B}P_\beta w(\beta)=0$. The proposal also records the bundle price invariant and the resulting inequality at each completed bundled iteration as companion statements. Gade et al., §3.3, pp. 9–10.

What the results provide

The scenario bound turns each price vector into an independent set of optimization problems whose weighted values certify a lower bound on the original mixed-integer optimum. The bound does not require the scenario first-stage decisions to agree at that iteration. The bundle result gives the same kind of certificate when several scenarios are solved together. These are the paper's known results, not new open mathematical conjectures. Gade et al., §§3.1 and 3.3.

Formalizing them requires a reusable interface for finite-support stochastic mixed-integer models, their scenario feasible sets, implementability equations, price systems, and finite prefixes of both algorithms. It also separates the general price-feasibility inequality from the algorithmic invariant. The mission currently proposes statements with open Lean proofs; it does not claim a machine-checked proof of these paper results.

Main mathematical difficulty

An iteration of progressive hedging solves scenario problems with a quadratic proximal term, while the lower-bound calculation uses scenario problems without that term. Their objective values cannot simply be identified. The relevant connection is the price vector: to apply Proposition 1 at an arbitrary iteration, its weighted zero-sum condition must follow from the actual aggregation and price-update rules. The bundle version adds normalization by $P_\beta$, so every bundle must have positive mass and the bundle weights must cover the full scenario distribution. Gade et al., §§3.1 and 3.3.

Formalization scope

Lean represents vectors as Fin n → ℝ, linear constraints componentwise, and the extensive formulation with the weighted implementability equation printed in (15). The dimensions satisfy $p_1\le n_1$ and $p_2\le n_2$, as required by the paper's spaces $\mathbb Z_+^{p_i}\times\mathbb R^{n_i-p_i}$. The finite scenario weights are positive and sum to one. The first integer coordinates are nonnegative; every remaining real coordinate is free. The paper's assumptions that the SMIP has an attained finite optimum and each $X(\xi)$ is nonempty are explicit in the goal and Propositions 1 and 3. The scenario and bundle subproblem values, and $z^*$, live in the extended reals: an unbounded-below subproblem retains $-\infty$ rather than acquiring an arbitrary real infimum. Gade et al., §§2–3.

An algorithm run is a finite prefix through its last aggregation and price update, with exact minimizers recorded as conditions on its sequences; a stopping test truncates the prefix. The proximal term is the Euclidean sum of squared coordinate differences. No positivity condition on $\rho$ is added to the lower-bound statements, since the page does not require one for these claims. Bundles form a partition, but equal bundle sizes are not required for the inequalities. Algorithm 2's printed w^ν(i) at initialization and Σ_i P_β x^ν(β) at aggregation are read with bundle indices. A bundle recourse vector is represented as a function on all scenarios, with only its coordinates inside the bundle constrained or used.

The exact equations of the scenario model and runs prevent a trivial encoding that assumes dual feasibility or the desired bound as part of a run. Contributions to the model, the price invariants, the extended-real weak-duality arguments, and the bundle extension are all within scope. The paper's convexified limit result (Proposition 2), its cited Theorem 1, and multi-stage Algorithm 3 are outside this mission.

Selected references

• D. Gade, G. Hackebeil, S. M. Ryan, J.-P. Watson, R. J-B Wets, and D. L. Woodruff, Obtaining Lower Bounds from the Progressive Hedging Algorithm for Stochastic Mixed-Integer Programs, author manuscript SAND2013-9195J; published in Mathematical Programming, 2016. OSTI manuscript; DOI.

Motivation

Benders decomposition solves an optimization problem in two groups of variables by fixing one group, solving the remaining subproblem, and turning what the subproblem's solution proves into a constraint, the Benders cut, on the fixed variables. In the classical method of Benders (1962) the subproblem is a linear program and the cut comes from its LP dual. Many problems that decompose naturally have subproblems that are not linear programs: scheduling, satisfiability, 0-1 and constraint programs. Hooker and Ottosson (Math. Program. 96, 2003) replaced the LP dual by the inference dual, the problem of inferring the strongest bound on the objective from the constraints, and obtained a Benders scheme in which any sound inference method produces cuts. This logic-based Benders decomposition has since become a standard technique for planning and scheduling problems that combine an assignment master problem with combinatorial subproblems.

The correctness of the generic scheme rests on two short theorems of the paper. Theorem 1 says that if every cut is valid, then each of the algorithm's three ways of stopping reports the right answer. Theorem 2 says that the algorithm stops when the fixed variables range over a finite set. This mission formalizes both theorems, together with the facts used in the proof of Theorem 1.

Setting

An optimization problem (1) has a domain $D$, a feasible set $S$ and a real-valued objective $f$; its value is $\inf_{x\in S} f(x)$. Optimal values lie in $[-\infty,+\infty]$: an infeasible minimization has value $+\infty$, an unbounded one $-\infty$, and the reverse holds for maximization. For propositions $P,Q$ about $x$, $P$ implies $Q$ with respect to $D$, written $P\xrightarrow{D}Q$, if $Q(x)$ holds at every $x\in D$ where $P(x)$ holds. The inference dual (2) of (1) is

$$\max\ \beta\quad\text{s.t.}\quad x\in S\xrightarrow{D} f(x)\ge\beta .$$

In Benders decomposition the variables are pairs $(x,y)$ with $x\in D_x$, $y\in D_y$, and problem (6) is $\min f(x,y)$ over $(x,y)\in S$. Fixing $y$ at a trial value $\bar y$ gives the subproblem (7), $\min f(x,\bar y)$ over $(x,\bar y)\in S$, with inference dual (8): $\max\beta$ s.t. $(x,\bar y)\in S\xrightarrow{D_x} f(x,\bar y)\ge\beta$. A bounding function $\beta_{\bar y}:D_y\to[-\infty,+\infty]$ defines the cut $z\ge\beta_{\bar y}(y)$. The cut satisfies

(B1) if every feasible $(x,y)$ of (6) satisfies $f(x,y)\ge\beta_{\bar y}(y)$, and

(B2) if $\beta_{\bar y}(\bar y)=\beta$, the value obtained in the subproblem dual.

The cuts generated so far form the master problem (9): $\min z$ s.t. $z\ge\beta_{y^k}(y)$ for every cut $k$, $y\in D_y$.

The generic Benders algorithm (Figure 1 of the paper) starts with $\bar z=-\infty$ and some $\bar y\in D_y$. While the subproblem dual at $\bar y$ has a feasible $\beta>\bar z$, it formulates $\beta_{\bar y}$ with $\beta_{\bar y}(\bar y)=\beta$ and adds its cut. If the master is then infeasible, it stops and reports (6) infeasible. Otherwise it takes an optimal solution $(\bar z,\bar y)$ of the master and repeats. When the loop exits, the optimal value of (6) is reported as $\bar z$.

In Lean, LogicBenders.Generic.Setting defines all of these objects: optVal, subVal, IsDualFeasible, ValidCut, masterLHS, IsMasterOptimal and MasterInfeasible, and a Run with the predicates IsRun, StopsAtWhile and StopsAtMaster.

Formalization targets

Goal: Theorem 1 (p. 9)

Assume (B1) for every cut the run has added. Then three statements hold.

1. If the run stops at the While test with a finite master value $\bar z$, then (6) has an optimal solution $(\bar x,\bar y)$ with $f(\bar x,\bar y)=\bar z$.
2. If it stops with an infeasible master, then $S=\varnothing$.
3. If it stops because the subproblem dual has no real feasible $\beta$, then (6) is unbounded.

$$\text{(B1)}\ \Longrightarrow\ \bigl[\text{finite stop}\Rightarrow \exists\,\bar x:\ f(\bar x,\bar y)=\bar z=\min_{S} f\bigr]\wedge\bigl[\text{master infeasible}\Rightarrow S=\varnothing\bigr]\wedge\bigl[\text{dual infeasible}\Rightarrow \inf_S f=-\infty\bigr].$$

Milestones (statements used in the proof)

• §3, p. 6: strong inference duality, $\inf_{x\in S}f(x)=\sup\{\beta: x\in S\xrightarrow{D}f(x)\ge\beta\}$.
• The subproblem value at any $\bar y$ bounds the value of (6) from above.
• Under (B1), the value $\bar z$ of a master optimum bounds the value of (6) from below.
• At a terminating master optimum, $\beta^*=\bar z$.
• Under (B1), an infeasible master implies that (6) is infeasible.
• An infeasible subproblem dual makes the subproblem, and hence (6), unbounded.

Companion: Theorem 2 (p. 10)

If (B1) and (B2) hold, $D_y$ is finite, and the subproblem dual is solved to optimality, then no run of the algorithm continues forever.

Further items

Lemma 3 (p. 20) gives a criterion for a 0-1 inequality $ax\ge\alpha$ to imply the branching clause (29). The §4 statement (p. 7) characterizes when a feasible linear system implies $cx\ge\beta$.

Significance

Theorem 1 is what allows any sound inference method to supply Benders cuts. Correctness needs only (B1), so the method that derives the cuts (resolution, constraint propagation, a branch-and-bound proof, LP duality) never has to be re-verified at the level of the decomposition. The paper's specializations to satisfiability, 0-1 programming and machine scheduling, and much of the later literature on logic-based Benders, rely on this theorem when they check validity of their cuts and nothing more. Theorem 2 is the matching termination guarantee, and the example on p. 10 shows that the finiteness of $D_y$ cannot be dropped.

Both theorems are proved in the paper, and to the extent known no machine-checked version exists. A formalization makes the paper's implicit conventions precise: values in $[-\infty,+\infty]$, what "the subproblem dual is infeasible" means when $-\infty$ is always a feasible bound, and exactly which cuts the master contains at each iteration. It also exposes two gaps in the printed proofs. Theorem 1's first clause uses the attainment of the subproblem optimum without stating it. Theorem 2's printed proof claims the subproblem value increases at every iteration, which is false, although the theorem itself is true. The definitions form a small reusable layer for later formalizations of specific logic-based Benders schemes.

Difficulty

The mathematics is elementary, and the difficulty lies in stating it faithfully. Several readings of the statements are wrong. A "run" whose $(\bar z,\bar y)$ are not master optima for exactly the cuts generated so far makes Theorem 1 false. Optimal values taken as real infima turn the infeasible case into junk. Defining "dual infeasible" over the extended reals makes clause 3 vacuous. For Theorem 2, the obvious argument is the printed one, that $\beta$ (or $\bar z$) increases strictly. It fails: with three values of $y$ and cuts equal to $-100$ away from their own $\bar y$, the subproblem value can decrease between iterations while $\bar z$ stays fixed. The printed argument therefore cannot be transcribed, and the theorem needs a proof the page does not give.

Formalization scope

• Domains and values. Domains are Lean types: X stands for $D_x$, Y for $D_y$ and D for $D$. The feasible set is S : Set (X × Y) and the objective is f : X → Y → ℝ. Optimal values are EReal infima and suprema, with $+\infty$ for infeasible and $-\infty$ for unbounded problems.
• Dual feasibility and the master problem. Dual feasibility ranges over EReal, so $\beta=-\infty$ is always feasible. "The subproblem dual is infeasible" is therefore stated as "no real $\beta$ is feasible". A master optimum $(z,\bar y)$ means $z=\max_{j\le k}\beta^{(j)}(\bar y)<+\infty$ and $z$ is at most that maximum at every $y$; the value $z=-\infty$ means the master is unbounded.
• Iterations. Iterations are 0-based: iteration $k$ works at ybar k, has bound zbar k, and adds cut k. The structure IsRun encodes Figure 1 exactly: zbar 0 = ⊥, every iteration finds a dual-feasible beta k > zbar k with cut k (ybar k) = beta k, and the next $(\bar z,\bar y)$ is a master optimum for exactly the cuts 0, …, k. An unconstrained "run", or a dual infeasibility that $-\infty$ satisfies vacuously, would trivialize the goal, and both are ruled out by these definitions.
• Added hypotheses. Theorem 1's first clause assumes that the subproblem optimum at the terminal $\bar y$ is attained whenever it is finite. The printed proof uses this fact ("some optimal solution $\bar x$"), and the clause is false without it. Lemma 3 adds $J_0\cap J_1=\varnothing$, which branching guarantees. No other hypotheses are added.
• Contributions welcome. Proofs of all items, especially a correct proof of Theorem 2.

Selected references

• J. N. Hooker and G. Ottosson, Logic-based Benders decomposition, Mathematical Programming 96(1):33–60, 2003. https://doi.org/10.1007/s10107-003-0375-9 (statements cited from the authors' revised manuscript of November 2000).
• J. F. Benders, Partitioning procedures for solving mixed-variables programming problems, Numerische Mathematik 4:238–252, 1962. https://doi.org/10.1007/BF01386316