Numerical Analysis

Completed

Optimization·Captain: mikedeng1

Proximal Alternating Linearized Minimization for Nonconvex and Nonsmooth Problems II: The Proximal Map of the Nonnegative Sparsity Constraint Is Hard Thresholding of the Positive PartResearch Paper

Motivation

Nonnegative matrix factorization (NMF) approximates a data matrix $A\in\mathbb R^{m\times n}$ by a product $XY$ of two nonnegative factors. Adding a cap on the number of nonzero entries of each factor gives sparse NMF, used in signal separation, text mining and image analysis because sparse nonnegative factors are easier to interpret. Bolte, Sabach and Teboulle (Math. Program. 146 (2014)) introduced the algorithm PALM (proximal alternating linearized minimization) for problems of the form $f(x)+g(y)+H(x,y)$ with nonconvex, nonsmooth $f$ and $g$ , and proved that its bounded runs converge to critical points when the objective has the Kurdyka–Łojasiewicz property. Sparse NMF is the paper's main application (§4).

PALM is only implementable when the proximal map of each nonsmooth term can be computed. For sparse NMF the nonsmooth term is the indicator of the nonconvex set of nonnegative matrices with at most $s$ nonzero entries. Proposition 4.1 of the paper computes that proximal map in closed form. This mission formalizes Proposition 4.1 and the steps of its proof. The convergence theory of PALM is a separate mission of the same series.

Setting

Matrices are real $m\times n$ arrays $X=(X_{ij})$ . The squared Frobenius norm is $\|X\|_F^2=\sum_{i,j}X_{ij}^2=\langle X,X\rangle$ . The sparsity count $\|X\|_0$ is the number of nonzero entries of $X$ , and $X\ge0$ means $X_{ij}\ge0$ for every $(i,j)$ .

For a set $C$ , the indicator $\delta_C$ is $0$ on $C$ and $+\infty$ off $C$ . For a natural number $s$ , set

f:=\delta_{X\ge0}+\delta_{\|X\|_0\le s},

the function that is $0$ on nonnegative matrices with at most $s$ nonzero entries and $+\infty$ elsewhere.

For $\sigma:\mathbb R^{m\times n}\to(-\infty,+\infty]$ and $t>0$ , the proximal map (2.2) of the paper is the set

\operatorname{prox}^\sigma_t(U)=\operatorname{argmin}\Big\{\sigma(X)+\frac t2\|X-U\|_F^2\Big\}.

Throughout, $\operatorname{argmin}\{\varphi(X):X\in C\}$ is the set of all points of $C$ where $\varphi$ attains its smallest value on $C$ .

The hard-thresholding operator of Definition 4.1 is

T_s(U)=\operatorname{argmin}_V\big\{\|U-V\|_F^2:\|V\|_0\le s\big\},

which is multi-valued when the $s$ largest entries of $U$ in absolute value are not uniquely determined. The projection onto the nonnegative orthant is $P_+(U)=\max\{0,U\}$ , componentwise.

In the proof, $U$ is fixed and $\mathcal I^+=\{(i,j):U_{ij}\ge0\}$ , $\mathcal I^-=\{(i,j):U_{ij}<0\}$ , with partial squared norms $\|X\|_\pm^2=\sum_{(i,j)\in\mathcal I^\pm}X_{ij}^2$ .

Formalization targets

Goal: Proposition 4.1 (Proximal map formula), p. 28

For every $U\in\mathbb R^{m\times n}$ and every $s\in\mathbb N$ ,

\operatorname{prox}^f_1(U)=\operatorname{argmin}\Big\{\tfrac12\|X-U\|_F^2 : X\ge0,\ \|X\|_0\le s\Big\}=T_s\big(P_+(U)\big).

Both equalities are equalities of sets, and both are part of the goal.

Milestones (in the order the proof uses them)

Relations (i)–(iii) in the proof (pp. 28–29): $\|X\|_F^2=\|X\|_+^2+\|X\|_-^2$ ; $\|X-U\|_+^2+\|X\|_-^2=\|X-P_+(U)\|_F^2$ ; and $\|X\|_-^2=0$ iff $X$ vanishes on $\mathcal I^-$ .
The equality of the minimizer sets (4.1) and (4.2) (p. 29):

\operatorname{argmin}\Big\{\|X-U\|_+^2+\|X\|_-^2-2\!\!\sum_{(i,j)\in\mathcal I^-}\!\!X_{ij}U_{ij}:X\ge0,\|X\|_0\le s\Big\}=\operatorname{argmin}\big\{\|X-U\|_+^2:X|_{\mathcal I^-}=0,\ X\ge0,\ \|X\|_0\le s\big\}.

The constraint $X\ge0$ in (4.2) can be removed without changing its set of solutions (p. 29).

Significance

The result. Proposition 4.1 says that the proximal step for the nonnegative sparsity constraint is computed by clipping the negative entries of $U$ to zero and then hard-thresholding: keeping $s$ entries of largest magnitude and zeroing the rest. The paper remarks that this costs $O(mn)$ operations. It is what turns PALM into the explicit algorithm PALM-Sparse NMF, and the same formula applies to any method that needs a projection onto nonnegative sparse matrices (projected gradient methods, iterative hard thresholding with sign constraints).

Formalizing it. The result is proved in the paper; to our knowledge no machine-checked proof exists. The platform has a projection-onto-sparse-vectors statement without nonnegativity and stated as a single inclusion, which is a different theorem. The mission produces a set-level statement that handles ties, the degenerate cases $s=0$ and $s\ge mn$ , and the interaction of the sign and the sparsity constraints, all of which the informal proof passes over quickly.

Difficulty

The obvious argument is: "project onto the nonnegative orthant, then onto the sparse matrices". Composing two projections onto nonconvex sets is not in general the projection onto their intersection, so this argument is not valid as it stands; the equality holds here because of the specific coordinate structure, and the proof must show it. The second obstacle is that everything is set-valued: when several entries of $P_+(U)$ tie for the $s$ -th largest magnitude, $T_s(P_+(U))$ has several elements, and the equality must account for every one of them, in both directions. The step "(4.1) = (4.2)" and the removal of $X\ge0$ are dismissed in the paper as "a simple contradiction argument" and "arguing in a similar way"; each needs an explicit modification of a feasible point that preserves sparsity and does not increase the objective.

Formalization scope

Matrices are Matrix (Fin m) (Fin n) ℝ, with 0-based indices (harmless relabelling of $\{1,\dots,m\}\times\{1,\dots,n\}$ ). $\|M\|_F^2$ is the Frobenius inner product $\langle M,M\rangle$ of the published definition CaiCandesShen.ProximalLimit.Basic. The indicators and $f$ take values in EReal, with $+\infty=\top$ ; membership in $\operatorname{prox}^f_1(U)$ is the inequality $f(X)+\frac12\|X-U\|_F^2\le f(W)+\frac12\|W-U\|_F^2$ in EReal for every $W$ . The proximal map carries the weight $t/2$ of (2.2), not $1/(2t)$ . $P_+$ is defined by the componentwise formula. No restriction is placed on $s$ .

Every operator in the statement is a set of minimizers. A formalization that selects one element (for instance "keep the $s$ largest entries with a fixed tie-break") proves only one inclusion of a weaker statement and does not count.

The development needs only finite sums and the argmin-set definition shared by all items; no further library is required. Lemmas about hard thresholding on matrices (the description of $T_s(U)$ by $s$ largest entries, nonnegativity of $T_s$ of a nonnegative matrix) are reusable well beyond this mission, and contributions of such lemmas are welcome.

Theorem numbers and page numbers refer to the author version of the paper (36 pages, printed page = PDF page), not to the Springer typesetting.

Selected references

J. Bolte, S. Sabach, M. Teboulle, Proximal alternating linearized minimization for nonconvex and nonsmooth problems, Mathematical Programming 146 (2014), 459–494. https://doi.org/10.1007/s10107-013-0701-9
D. D. Lee, H. S. Seung, Learning the parts of objects by non-negative matrix factorization, Nature 401 (1999), 788–791. https://doi.org/10.1038/44565
J.-F. Cai, E. J. Candès, Z. Shen, A singular value thresholding algorithm for matrix completion, SIAM J. Optim. 20 (2010), 1956–1982. https://doi.org/10.1137/080738970 (source of the Frobenius inner product definition reused here)

6 thms2 active usersReviewed

Completed

Optimization·Captain: mikedeng1

Sparsity Constrained Nonlinear Optimization: Optimality Conditions and Algorithms II: Iterative Hard Thresholding for Sparse Least Squares Converges to an L-Stationary PointResearch Paper

Motivation

Many estimation problems ask for a vector with few nonzero entries that fits data well. The best-known case is compressed sensing: recover a sparse $x\in\mathbb R^n$ from $m\ll n$ noisy linear measurements $b\approx Ax$ . Fitting such data under an explicit sparsity budget leads to the sparse least-squares problem

\min\ \|Ax-b\|^2\quad\text{s.t.}\quad \|x\|_0\le s,

a nonconvex problem that is NP-hard in general. A standard heuristic for it is iterative hard thresholding (IHT): take a gradient step on the least-squares term and keep the $s$ entries of largest magnitude. Blumensath and Davies analysed it first under the normalization $\|A\|_2<1$ (Blumensath, Davies 2008) and then with an adaptive step size (Blumensath, Davies 2010).

Beck and Eldar (SIAM J. Optim. 2013, arXiv:1203.4580) place IHT inside a general theory of necessary optimality conditions for minimizing a smooth function under a sparsity constraint. In that theory IHT is a fixed-point method for a condition they call $L$ -stationarity, and for least squares with a mild condition on $A$ they prove that the iterates converge with a fixed step size $1/L$ tied to the Lipschitz constant of the gradient. That convergence theorem is the goal of this mission.

Setting

Fix integers $0<s<n$ . For $x\in\mathbb R^n$ let $\|x\|_0$ be the number of nonzero entries of $x$ , let $I_1(x)=\{i: x_i\neq0\}$ be its support and $I_0(x)=\{i:x_i=0\}$ the complement, and let

C_s=\{x\in\mathbb R^n:\ \|x\|_0\le s\}.

Write $M_s(x)$ for the $s$ -th largest of $|x_1|,\dots,|x_n|$ , counted with multiplicity; $M_s(x)=0$ exactly when $\|x\|_0<s$ .

The problem is

(\mathrm P)\qquad \min\ f(x)\quad\text{s.t.}\quad x\in C_s,

with $f:\mathbb R^n\to\mathbb R$ continuously differentiable and bounded below (Assumption 1). Where stated, $\nabla f$ is Lipschitz with constant $L(f)$ (Assumption 2).

The projection onto $C_s$ is the set of nearest points $P_{C_s}(y)=\operatorname{argmin}_{x\in C_s}\|x-y\|^2$ . Since $C_s$ is not convex this set can have more than one element: $P_{C_2}((2,1,1)^T)=\{(2,1,0)^T,(2,0,1)^T\}$ . For $L>0$ , a point $x^*\in C_s$ is $L$ -stationary if

[\mathrm{NC}_L]\qquad x^*\in P_{C_s}\Big(x^*-\tfrac1L\nabla f(x^*)\Big).

A point $x^*\in C_s$ is a basic feasible (BF) vector if $\nabla f(x^*)=0$ when $\|x^*\|_0<s$ , and $\nabla_if(x^*)=0$ for $i\in I_1(x^*)$ when $\|x^*\|_0=s$ .

The IHT method with stepsize $1/L$ produces any sequence with

x^0\in C_s,\qquad x^{k+1}\in P_{C_s}\Big(x^k-\tfrac1L\nabla f(x^k)\Big)\quad(k=0,1,2,\dots).

A matrix $A\in\mathbb R^{m\times n}$ is $s$ -regular if every $s$ of its columns are linearly independent. The least-squares objective is $f_{LI}(x)=\|Ax-b\|^2$ , whose gradient is Lipschitz with $L(f_{LI})=2\lambda_{\max}(A^TA)$ .

Formalization targets

Goal: Theorem 3.2

For $f=f_{LI}$ with $A$ $s$ -regular, $L>L(f)$ and any IHT run $(x^k)$ , $\exists\,x^*:\quad x^*\ \text{is $L$-stationary}\quad\text{and}\quad x^k\to x^* .$

Milestones, in the order the proof uses them

Lemma 2.3 (descent lemma): $f(x)\le f(y)+\langle\nabla f(y),x-y\rangle+\frac L2\|x-y\|^2$ for $L\ge L(f)$ . Already proved on the platform; included as a reference.
Lemma 2.4: for $x\in C_s$ and $y\in P_{C_s}(x-\frac1L\nabla f(x))$ with $L>L(f)$ ,

f(x)-f(y)\ge\frac{L-L(f)}2\|x-y\|^2.

Lemma 3.1: along an IHT run with $L>L(f)$ the same decrease holds for consecutive iterates, $f(x^k)$ is nonincreasing, $\|x^k-x^{k+1}\|\to0$ , and every move strictly decreases $f$ .
Corollary 3.1: $f(x^k)$ converges.
Lemma 2.2: $[\mathrm{NC}_L]$ holds iff $\|x^*\|_0\le s$ , $\nabla_if(x^*)=0$ on $I_1(x^*)$ and $|\nabla_if(x^*)|\le L\,M_s(x^*)$ on $I_0(x^*)$ .
Theorem 3.1: every accumulation point of an IHT run with $L>L(f)$ is $L$ -stationary (general $f$ ).
Corollary 2.1: an $L$ -stationary point is BF.
Lemma 2.1: for $f_{LI}$ with $A$ $s$ -regular, the set of BF points is finite.

Significance

Theorem 3.1 alone gives only subsequential information, and it says nothing when the iterates have no accumulation point. Theorem 3.2 turns it into convergence of the whole sequence for sparse least squares, with an explicit admissible step size $1/L$ for any $L>2\lambda_{\max}(A^TA)$ . It removes the normalization $\|A\|_2<1$ that the original analysis of IHT required, and it identifies the limit as a point satisfying a necessary optimality condition that is strictly stronger than basic feasibility.

These results are proved on paper. To our knowledge none of them has a machine-checked proof. A completed mission would give a formal convergence theorem for a hard-thresholding algorithm that is used in practice, and a reusable layer for sparsity-constrained optimization: the set-valued projection onto $C_s$ , its coordinate characterization (Lemma 2.2), and the order statistic $M_s$ . The sibling missions of this series (coordinate-wise minima and sparse-simplex methods) use the same objects.

Difficulty

The projection onto $C_s$ is set-valued and discontinuous: a tiny perturbation of $y$ can change which $s$ entries are kept. So the usual argument for projected gradient methods on convex sets, which passes to the limit in a continuous fixed-point map, does not apply to Theorem 3.1. The limit has to be taken coordinate by coordinate, separately on the support of the accumulation point and off it.

For Theorem 3.2, sufficient decrease and $\|x^k-x^{k+1}\|\to0$ do not by themselves give convergence of a sequence. Two further facts are needed. The iterates must stay bounded, which uses $s$ -regularity on every support of size at most $s$ . The candidate limits must be isolated, which uses the finiteness of the BF points (Lemma 2.1). Neither holds for a general smooth $f$ , which is why the goal is specific to least squares.

Formalization scope

Vectors live in EuclideanSpace ℝ (Fin n). Indices are $0,\dots,n-1$ instead of $1,\dots,n$ , and $\|\cdot\|$ is the Euclidean norm. The gradient is Mathlib's gradient f with $\nabla_if(x)$ = gradient f x i, and $f$ is assumed ContDiff ℝ 1. The standing hypotheses $0<s<n$ appear in every statement. Assumption 1 (∃ γ, ∀ x, γ ≤ f x) is a hypothesis wherever $f$ is general, also where the proof does not use it, because the paper makes it throughout. It is omitted for $f_{LI}\ge0$ . Assumption 2 takes the form of a constant Lf : ℝ≥0 with LipschitzWith Lf (gradient f) together with Lf < L. The statements hold for every such constant, and for $f_{LI}$ this is equivalent to $L>2\lambda_{\max}(A^TA)$ . $M_s(x)$ is the entry at position $s-1$ of the absolute values sorted in nonincreasing order. $Ax$ is Matrix.toEuclideanLin A x. An accumulation point is MapClusterPt xstar atTop x, and convergence is Tendsto x atTop (𝓝 xstar).

$P_{C_s}(y)$ is formalized as the full set of minimizers. An IHT run is any sequence of selections from it, and every theorem quantifies over all runs. Formalizing the projection as one chosen top- $s$ vector would change $[\mathrm{NC}_L]$ and Theorem 3.1, and is ruled out. The theorems assume the strict inequality $L>L(f)$ that the paper's results use, not the algorithm box's input condition $L\ge L(f)$ . A sanity check confirms that the hypotheses of Theorem 3.1 are jointly satisfiable, so the run predicate is not vacuous.

Theorem numbers and page references are those of the arXiv version (arXiv:1203.4580v1), which differs in typesetting from the SIAM publication.

The objects sit in two definition files: the shared SparseNLO.CWOpt.Setting ( $\|x\|_0$ , $I_1$ , $I_0$ , $C_s$ , $M_s$ , BF vectors) and SparseNLO.IHT.Setting (the projection, $L$ -stationarity, $s$ -regularity, $f_{LI}$ , IHT runs). A complete development needs the following:

the closed-form description of $P_{C_s}$ (keep $s$ largest entries);
continuity of $M_s$ ;
boundedness of the least-squares sublevel sets restricted to supports of size $\le s$ ;
a standard argument that a sequence with vanishing steps and finitely many isolated cluster candidates converges.

Proofs of any milestone are welcome, as are reusable lemmas on top- $s$ projections and order statistics.

Selected references

A. Beck, Y. C. Eldar, Sparsity Constrained Nonlinear Optimization: Optimality Conditions and Algorithms, SIAM J. Optim. 23(3), 2013. https://doi.org/10.1137/120869778 (arXiv:1203.4580v1, https://arxiv.org/abs/1203.4580)
T. Blumensath, M. E. Davies, Iterative Thresholding for Sparse Approximations, J. Fourier Anal. Appl. 14, 2008. https://doi.org/10.1007/s00041-008-9035-z
T. Blumensath, M. E. Davies, Normalized Iterative Hard Thresholding: Guaranteed Stability and Performance, IEEE J. Sel. Topics Signal Process. 4(2), 2010. https://doi.org/10.1109/JSTSP.2010.2042411
D. P. Bertsekas, Nonlinear Programming, 2nd ed., Athena Scientific, 1999 (descent lemma).

11 thms2 active usersReviewed

Completed

Optimization·Captain: mikedeng1

Proximal Alternating Linearized Minimization for Nonconvex and Nonsmooth Problems I: Every Bounded PALM Sequence of a KL Function Has Finite Length and Converges to a Critical PointResearch Paper

Motivation

Many problems in data analysis and signal processing are posed as minimizing a sum of three terms over two blocks of variables: a nonsmooth function of the first block, a nonsmooth function of the second block, and a smooth coupling term. Nonnegative and sparse matrix factorization, blind deconvolution and dictionary learning all have this form, with constraints or sparsity penalties as the nonsmooth parts. The functions involved are usually nonconvex, so convexity-based convergence theory does not apply, and the classical alternating minimization (Gauss–Seidel) scheme requires solving an exact nonconvex subproblem in each block.

Bolte, Sabach and Teboulle (Math. Program. 146 (2014)) introduced Proximal Alternating Linearized Minimization (PALM), which replaces each exact block minimization by one proximal-gradient step, and proved that every bounded sequence it generates converges to a critical point whenever the objective has the Kurdyka–Łojasiewicz property. The result became a standard template for the convergence analysis of first-order methods on nonconvex, nonsmooth problems.

Timeline.

1963 and 1993: Łojasiewicz proves the gradient inequality for real-analytic functions; Kurdyka extends it to definable functions.
2007: Bolte, Daniilidis, Lewis extend it to nonsmooth subanalytic functions.
2009: Attouch and Bolte prove convergence of the proximal point method under the Łojasiewicz inequality.
2010: Attouch, Bolte, Redont and Soubeyran prove it for the exact proximal alternating scheme.
2014: the present paper treats PALM under the general KL property, with only linearized (proximal-gradient) block steps.

Setting

Let $f:\mathbb R^n\to(-\infty,+\infty]$ and $g:\mathbb R^m\to(-\infty,+\infty]$ , let $H:\mathbb R^n\times\mathbb R^m\to\mathbb R$ , and set

\Psi(x,y)=f(x)+g(y)+H(x,y).

The product $\mathbb R^n\times\mathbb R^m$ carries the Euclidean norm $\|(x,y)\|^2=\|x\|^2+\|y\|^2$ . Assumption A: $f$ and $g$ are proper (never $-\infty$ , finite somewhere) and lower semicontinuous, and $H$ is $C^1$ . Assumption B:

$\inf\Psi$ , $\inf f$ , $\inf g$ are finite;
for each $y$ , $x\mapsto\nabla_xH(x,y)$ is Lipschitz with modulus $L_1(y)$ ; for each $x$ , $y\mapsto\nabla_yH(x,y)$ is Lipschitz with modulus $L_2(x)$ ;
along the generated sequence, $\lambda_1^-\le L_1(y^k)\le\lambda_1^+$ and $\lambda_2^-\le L_2(x^k)\le\lambda_2^+$ for positive constants $\lambda_i^\pm$ ;
$\nabla H$ is Lipschitz on bounded sets.

The proximal map of $\sigma$ with parameter $t>0$ is the set

\operatorname{prox}^\sigma_t(x)=\operatorname*{argmin}_u\Big\{\sigma(u)+\tfrac t2\|u-x\|^2\Big\}.

Fix $\gamma_1,\gamma_2>1$ . A sequence $z^k=(x^k,y^k)$ is generated by PALM if, from an arbitrary $(x^0,y^0)$ , with $c_k=\gamma_1L_1(y^k)$ and $d_k=\gamma_2L_2(x^{k+1})$ ,

x^{k+1}\in\operatorname{prox}^f_{c_k}\Big(x^k-\tfrac1{c_k}\nabla_xH(x^k,y^k)\Big),\qquad y^{k+1}\in\operatorname{prox}^g_{d_k}\Big(y^k-\tfrac1{d_k}\nabla_yH(x^{k+1},y^k)\Big).

The limiting subdifferential $\partial\sigma(x)$ consists of the limits $v$ of Fréchet subgradients $v^k$ at points $x^k\to x$ with $\sigma(x^k)\to\sigma(x)$ ; $z$ is critical if $0\in\partial\Psi(z)$ . Let $\Phi_\eta$ be the class of continuous concave functions $\varphi:[0,\eta)\to\mathbb R_+$ with $\varphi(0)=0$ , $C^1$ on $(0,\eta)$ and $\varphi'>0$ there. A proper lsc $\sigma$ has the KL property at $\bar u$ if there are $\eta>0$ , a neighbourhood $U$ of $\bar u$ and $\varphi\in\Phi_\eta$ with

\varphi'\big(\sigma(u)-\sigma(\bar u)\big)\operatorname{dist}\big(0,\partial\sigma(u)\big)\ge1\quad\text{whenever }u\in U,\ \sigma(\bar u)<\sigma(u)<\sigma(\bar u)+\eta,

and $\sigma$ is a KL function if it has this property at every point where $\partial\sigma$ is nonempty.

Formalization targets

Goal: Theorem 3.1 (p. 21)

If $\Psi$ is a KL function, Assumptions A and B hold and $(z^k)$ is a bounded sequence generated by PALM, then

\sum_{k=1}^\infty\|z^{k+1}-z^k\|<\infty\qquad\text{and}\qquad z^k\to z^*\ \text{with}\ 0\in\partial\Psi(z^*).

Milestones, in proof order

Lemma 3.1 (descent lemma) and Proposition 2.1 ( $\partial\Psi=(\nabla_xH+\partial f,\nabla_yH+\partial g)$ ): already on the platform, reused.
Lemma 3.2: one proximal-gradient step with $t>L_h$ gives $h(u^+)+\sigma(u^+)\le h(u)+\sigma(u)-\tfrac12(t-L_h)\|u^+-u\|^2$ .
Lemma 3.3: $\tfrac{\rho_1}2\|z^{k+1}-z^k\|^2\le\Psi(z^k)-\Psi(z^{k+1})$ with $\rho_1=\min\{(\gamma_1-1)\lambda_1^-,(\gamma_2-1)\lambda_2^-\}$ , and $\sum_k\|z^{k+1}-z^k\|^2<\infty$ .
Lemma 3.4: explicit vectors $(A^k_x,A^k_y)\in\partial\Psi(z^k)$ with $\|A^k_x\|+\|A^k_y\|\le(2M+3\rho_2)\|z^k-z^{k-1}\|$ , $\rho_2=\max\{\gamma_1\lambda_1^+,\gamma_2\lambda_2^+\}$ .
Lemma 3.5: the limit point set $\omega(z^0)$ is nonempty, compact, connected and contained in $\operatorname{crit}\Psi$ , attracts the sequence, and $\Psi$ is finite and constant on it.
Lemma 3.6: on a compact set where $\sigma$ is constant, the KL property holds with a single $(\varepsilon,\eta,\varphi)$ .

Significance

The result. Theorem 3.1 gives convergence of the whole sequence, not only subsequential convergence, for a method whose steps are explicit: a gradient step on $H$ followed by a proximal step on $f$ or $g$ . The KL hypothesis holds for semi-algebraic, subanalytic and o-minimal data, which covers $\ell_0$ penalties, rank constraints and polynomial couplings. The paper applies the theorem to sparse nonnegative matrix factorization (Theorem 4.1), and the same scheme (sufficient decrease, relative error, KL) underlies many later convergence proofs for nonconvex splitting and alternating methods.

Formalizing it. The result is proved on paper; no machine-checked proof is known to exist. A formalization produces reusable components: the nonconvex proximal-gradient decrease lemma, the structure of limit point sets of sequences with vanishing steps, and the uniformized KL lemma, which is independent of PALM and is used in essentially every KL-based convergence proof.

Difficulty

Lemmas 3.2–3.5 are estimates and compactness arguments; the core difficulty is in passing from them to convergence. Square-summability of the steps (Lemma 3.3) does not give summability of the steps themselves, and a bounded sequence with vanishing steps may wander along a continuum of limit points. The KL inequality supplies the missing control, but only near each critical point separately, with point-dependent $\eta$ and $\varphi$ ; the iterates are known only to approach the compact set $\omega(z^0)$ , not a single point. Combining Lemma 3.6 with the decrease and relative-error estimates requires handling the concave function $\varphi$ along the sequence and an extended-valued $\Psi$ whose subdifferential is a limiting construction. Lemma 3.5(i) also needs $\Psi(z^{k_q})\to\Psi(z^*)$ along subsequences, which lower semicontinuity alone does not give.

Formalization scope

$\mathbb R^n$ , $\mathbb R^m$ are EuclideanSpace ℝ (Fin n), EuclideanSpace ℝ (Fin m), and the product is the $L^2$ product WithLp 2 (… × …) of the published setting ProxAltMin.Conv.Setting, not Lean's default product with the sup norm. Values in $(-\infty,+\infty]$ are EReal. The prox map is a set with weight $t/2$ , exactly as in (2.2). Iterates are selections from it, and the starting point is arbitrary, so $\Psi(z^0)=+\infty$ is allowed; inequalities between values of $\Psi$ are stated with the difference moved across. $L_1,L_2$ are given modulus functions and $\gamma_1,\gamma_2$ fixed constants. The limiting subdifferential and the KL notions are the published definitions NonconvexSplitting.Shared.LimitingSubdiff and NonconvexSplitting.ADMMKL.KLProperty, with $\eta$ real and the KL inequality read for every $v\in\partial\sigma(u)$ (equivalent to the $\operatorname{dist}$ form). Assumption A together with B(iv) is AssumptionH; B(i)–(iii) and the PALM run are defined in PALM.Conv.Setting. The satisfiability of all hypotheses of the goal has been checked on a concrete instance.

The goal is not weakened: it asserts summability of $\|z^{k+1}-z^k\|$ and convergence of the full sequence to one critical point, not that every cluster point is critical (that is Lemma 3.5(i)). The constants $\rho_1$ and $2M+3\rho_2$ are explicit as in the paper, with $M$ chosen once for all $k$ .

Theorem numbers and pages are those of the author version (source PDF), whose page numbers coincide with its printed pages; they may differ from the Springer typesetting.

Welcome contributions: proofs of any milestone; a general-space version of Lemma 3.2; a proof of the published Proposition 2.1 (ProxAltMin.Conv.proposition_3, still open), which Lemma 3.4 relies on.

Selected references

J. Bolte, S. Sabach, M. Teboulle, Proximal alternating linearized minimization for nonconvex and nonsmooth problems, Math. Program. 146 (2014) 459–494. https://doi.org/10.1007/s10107-013-0701-9
H. Attouch, J. Bolte, P. Redont, A. Soubeyran, Proximal alternating minimization and projection methods for nonconvex problems: an approach based on the Kurdyka–Łojasiewicz inequality, Math. Oper. Res. 35 (2010) 438–457. https://doi.org/10.1287/moor.1100.0449
H. Attouch, J. Bolte, On the convergence of the proximal algorithm for nonsmooth functions involving analytic features, Math. Program. 116 (2009) 5–16. https://doi.org/10.1007/s10107-007-0133-5
J. Bolte, A. Daniilidis, A. Lewis, The Łojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems, SIAM J. Optim. 17 (2007) 1205–1223. https://doi.org/10.1137/050644641
R. T. Rockafellar, R. J.-B. Wets, Variational Analysis, Springer, 1998. https://doi.org/10.1007/978-3-642-02431-3

17 thms2 active usersReviewed

Completed

Optimization·Captain: mikedeng1

Sparsity Constrained Nonlinear Optimization: Optimality Conditions and Algorithms III: Accumulation Points of the Greedy Sparse-Simplex Method Are Coordinate-Wise MinimaResearch Paper

Motivation

Many estimation problems in signal processing and statistics look for a vector with few nonzero entries that fits data well. Compressed sensing (Donoho 2006) recovers a sparse $x$ from linear measurements by minimizing $\|Ax-b\|^2$ over sparse vectors. Sub-wavelength optical imaging and phase retrieval (Shechtman, Eldar, Szameit, Segev 2011) lead to the same task with a quartic, nonconvex objective, because the measurements are quadratic in the unknown image. Beck and Eldar (SIAM J. Optim. 2013) treat the common abstraction: minimize a general smooth function under a cardinality constraint.

For linear least squares, greedy methods such as matching pursuit (Mallat, Zhang 1993) and orthogonal matching pursuit (Tropp 2004) add one coordinate at a time, and iterative hard thresholding (Blumensath, Davies 2008) is a projected gradient method. For a general nonconvex objective the question is which optimality condition an algorithm can actually guarantee at its limit points. Beck and Eldar introduce the greedy sparse-simplex method, a coordinate-descent scheme that never leaves the feasible set, and prove that its limit points satisfy coordinate-wise optimality, the strongest of the necessary conditions studied in their paper. This mission formalizes that convergence theorem.

Setting

Let $n$ and $s$ be integers with $0<s<n$ , and let $f:\mathbb R^n\to\mathbb R$ be continuously differentiable and bounded below: there is $\gamma$ with $f(x)\ge\gamma$ for all $x$ (Assumption 1, a standing assumption of the paper). The problem is

(\mathrm P)\qquad \min\ f(x)\quad\text{s.t.}\quad \|x\|_0\le s,

where the $\ell_0$ norm $\|x\|_0$ counts the nonzero components of $x$ . Write $I_1(x)=\{i: x_i\ne0\}$ for the support of $x$ , $I_0(x)=\{i: x_i=0\}$ for its complement, $C_s=\{x:\|x\|_0\le s\}$ for the feasible set, and $e_i$ for the $i$ -th unit vector.

A feasible $x^*$ is a coordinate-wise (CW) minimum of (P) (Definition 2.4) if

Case I: $\|x^*\|_0<s$ and $f(x^*)\le f(x^*+te_i)$ for every $i$ and every $t\in\mathbb R$ ; or
Case II: $\|x^*\|_0=s$ and $f(x^*)\le f(x^*-x^*_ie_i+te_j)$ for every $i\in I_1(x^*)$ , every $j$ and every $t\in\mathbb R$ .

Below the sparsity level no single coordinate can be changed profitably; at the sparsity level no support coordinate can be profitably swapped for any coordinate carrying any value.

The greedy sparse-simplex method starts from $x^0\in C_s$ . If $\|x^k\|_0<s$ , it minimizes $f$ over all points $x^k+te_i$ (all $i$ , all $t$ ); if $\|x^k\|_0=s$ , it minimizes $f$ over all points $x^k-x^k_ie_i+te_j$ with $i\in I_1(x^k)$ , any $j$ and any $t$ . If the best candidate strictly improves on $f(x^k)$ it becomes $x^{k+1}$ ; otherwise the method stops. A run is any sequence $(x^k)$ produced this way, for any choice of minimizers.

Formalization targets

Goal: Theorem 3.3

For every run $(x^k)$ of the greedy sparse-simplex method and every accumulation point $x^*$ of $(x^k)$ ,

x^*\ \text{is a CW-minimum of (P)}.

No Lipschitz assumption on $\nabla f$ is made.

Milestones

Lemma 3.2. $f(x^{k+1})\le f(x^k)$ for every $k$ , with equality if and only if $x^{k+1}=x^k$ and $x^k$ is a CW-minimum.
Case $\|x^*\|_0=s$ (§4.2.1). $f(x^*)\le f(x^*-x^*_ie_i+te_j)$ for all $i\in I_1(x^*)$ , all $j$ , all $t$ .
(4.2). If $\|x^*\|_0<s$ and $i\in I_1(x^*)$ , then $f(x^*)\le f(x^*+te_i)$ for all $t$ .
(4.5) and the display after it. If $\|x^*\|_0<s$ and $i\in I_0(x^*)$ , then $f(x^*)\le f(x^*+te_i)$ for all $t$ .

Milestones 2–4 together with closedness of $C_s$ give the goal.

Significance

The theorem certifies what a cheap, feasible-point method achieves on a nonconvex, discontinuously constrained problem: its limit points cannot be improved by changing a single coordinate or by swapping a single support coordinate. Combined with Theorem 2.4 of the same paper (every CW-minimum is an $L_2(f)$ -stationary point), it gives Corollary 3.2: when $\nabla f$ is Lipschitz, the limit points are $L_2(f)$ -stationary, where $L_2(f)$ is a local Lipschitz constant that can be much smaller than the global one used by iterative hard thresholding, and the method needs no knowledge of either constant. For quartic objectives arising in phase retrieval each inner step reduces to minimizing a scalar quartic, so the guarantee applies to a practical algorithm.

The result is proved in the paper (§4.2.1). As far as a search of the Prove2Me library shows, neither coordinate-wise optimality for sparsity constrained problems nor any convergence theorem for sparse-simplex methods has a machine-checked proof. The mission produces a formal model of the method that is faithful to its selection rules and stopping rule, and a formal proof of the convergence theorem. The companion missions of the series formalize Theorem 2.4 and the partial sparse-simplex method.

Difficulty

The obvious argument passes to the limit in the inequality "the next iterate is at least as good as every candidate move from the current iterate" along a subsequence converging to $x^*$ . This works directly when the candidate moves from $x^{k_n}$ converge to the candidate moves from $x^*$ . It fails in the case $\|x^*\|_0<s$ while the iterates near $x^*$ sit at the sparsity level $s$ : there the method cannot simply add a coordinate $i\in I_0(x^*)$ , only swap one out, and the candidate set from $x^{k_n}$ is not the candidate set from $x^*$ . The discontinuity of $\|\cdot\|_0$ means that the regime of the iterates (below or at the sparsity level) need not be that of the limit.

Formalization scope

Vectors live in EuclideanSpace ℝ (Fin n); indices are $0,\dots,n-1$ in Lean and $1,\dots,n$ on the page, and $te_i$ is EuclideanSpace.single i t. The function $f$ is ContDiff ℝ 1, Assumption 1 is the hypothesis ∃ γ, ∀ x, γ ≤ f x, and 0 < s, s < n are hypotheses, as in (P). No Lipschitz hypothesis appears anywhere.

A run is a predicate IsGreedyRun f s x on a sequence x : ℕ → EuclideanSpace ℝ (Fin n): $x^0\in C_s$ , and each step is either a move to a point that strictly decreases $f$ and attains the minimum of $f$ over all candidate points, or a STOP, in which $f(x^k)$ is at most every candidate value and the sequence stays put, $x^{k+1}=x^k$ (the convention behind Lemma 3.2). The argmin selections are therefore arbitrary, and every theorem holds for every run. An accumulation point is MapClusterPt xstar atTop x. Membership $x^*\in C_s$ is part of the definition of a CW-minimum and has to be proved.

The predicate is not vacuous: a sorry-free check shows that on $n=2$ , $s=1$ , $f(x)=(x_1-1)^2+x_2^2$ , the sequence $x^0=0$ , $x^k=e_1$ for $k\ge1$ is a run with one move followed by a STOP. A formalization in which the STOP branch is missing, or in which accumulation points are replaced by limits, would prove a weaker or empty statement and is ruled out by this encoding.

The needed infrastructure is small: closedness of $C_s$ , the behaviour of supports along convergent sequences (supports eventually contain the support of the limit), and continuity arguments for limits of inequalities. A lemma that the support of nearby points contains the support of the limit is reusable for any $\ell_0$ -constrained problem. Proofs of the milestones as stated, and alternative proofs of the goal, are welcome.

Theorem and page numbers refer to the arXiv version (arXiv:1203.4580v1), whose printed and PDF pages coincide; the published version is SIAM J. Optim. 23(3), 2013.

Selected references

A. Beck, Y. C. Eldar, Sparsity Constrained Nonlinear Optimization: Optimality Conditions and Algorithms, SIAM Journal on Optimization 23(3), 2013. https://doi.org/10.1137/120869778 — arXiv:1203.4580v1 (2012), https://arxiv.org/abs/1203.4580
T. Blumensath, M. E. Davies, Iterative Thresholding for Sparse Approximations, Journal of Fourier Analysis and Applications 14(5), 2008. https://doi.org/10.1007/s00041-008-9035-z
D. L. Donoho, Compressed Sensing, IEEE Transactions on Information Theory 52(4), 2006. https://doi.org/10.1109/TIT.2006.871582
S. Mallat, Z. Zhang, Matching Pursuits with Time-Frequency Dictionaries, IEEE Transactions on Signal Processing 41(12), 1993. https://doi.org/10.1109/78.258082
Y. Shechtman, Y. C. Eldar, A. Szameit, M. Segev, Sparsity Based Sub-Wavelength Imaging with Partially Incoherent Light via Quadratic Compressed Sensing, Optics Express 19(16), 2011. https://doi.org/10.1364/OE.19.014807
J. A. Tropp, Greed Is Good: Algorithmic Results for Sparse Approximation, IEEE Transactions on Information Theory 50(10), 2004. https://doi.org/10.1109/TIT.2004.834793

7 thms2 active usersReviewed

Convex Optimization Operations Research Optimization·Captain: mikedeng1

Efficiency of Coordinate Descent Methods on Huge-Scale Optimization Problems 4: Accelerated Coordinate Descent ACDM Has Expected Error at Most (n/(k + 1))²·(2‖x₀ − x‖₁² + (f(x₀) − f)/n²)Research Paper

Motivation

Coordinate descent methods update one block of variables at a time. For problems with millions or billions of variables, where even one full gradient is too expensive to compute or store, a step that touches a single block can be orders of magnitude cheaper than a full gradient step. Nesterov's Efficiency of coordinate descent methods on huge-scale optimization problems (CORE Discussion Paper 2010/2; journal version SIAM J. Optim. 22 (2012) 341–362, doi:10.1137/100802001) gave the first global complexity bounds for randomized coordinate descent on smooth convex functions: choosing the block at random makes the expected progress per step comparable to that of a full gradient step divided by the number of blocks.

Section 5 of the paper asks whether the multistep acceleration of the full gradient method (Nesterov 1983) carries over to the randomized coordinate setting, and answers it with the accelerated coordinate descent method ACDM. This mission formalizes that answer: Theorem 6, the $O(n^2/k^2)$ expected rate of ACDM, together with the lemmas of its proof.

Timeline. Nesterov (1983) introduced the accelerated full-gradient method with rate $O(1/k^2)$ . Nesterov (CORE DP 2010/2; SIAM J. Optim. 2012) proved the $O(n/k)$ rate of random block coordinate descent and the $O(n^2/k^2)$ rate of its accelerated variant ACDM. Later work (Lee and Sidford 2013; Fercoq and Richtárik 2015, arXiv:1312.5799; Allen-Zhu, Qu, Richtárik and Yuan 2016, arXiv:1512.09103) made the accelerated scheme efficient per iteration and sharpened its constants.

Setting

The space is $\mathbb R^N=\mathbb R^{n_1}\times\cdots\times\mathbb R^{n_n}$ , split into $n\ge1$ blocks. A point $x$ has blocks $x^{(i)}\in\mathbb R^{n_i}$ , and $U_ih$ is the point whose $i$ -th block is $h$ and whose other blocks vanish. Each block carries a Euclidean norm $\|h\|_{(i)}^2=\langle B_ih,h\rangle$ with $B_i\succ0$ .

The objective $f:\mathbb R^N\to\mathbb R$ is convex and differentiable, and has a coordinate-wise Lipschitz gradient: the partial gradient $f'_i(x)=U_i^T\nabla f(x)$ satisfies $\|f'_i(x+U_ih)-f'_i(x)\|^*_{(i)}\le L_i\|h\|_{(i)}$ with constants $L_i>0$ . The method weighs blocks by these constants through the norm

\|x\|_1^2=\sum_{i=1}^nL_i\|x^{(i)}\|_{(i)}^2 .

The function $f$ is strongly convex in this norm with parameter $\sigma\ge0$ if $f(y)\ge f(x)+\langle\nabla f(x),y-x\rangle+\frac\sigma2\|y-x\|_1^2$ for all $x,y$ ; $\sigma=0$ is plain convexity. A point $x^*$ minimizes $f$ , and $f^*=f(x^*)$ .

The optimal coordinate step is $T_i(x)=x-\frac1{L_i}U_if'_i(x)^\#$ , where $s^\#$ maximizes $\langle s,x\rangle-\frac12\|x\|^2$ . ACDM $(x_0)$ keeps two sequences $x_k,v_k$ with $v_0=x_0$ and deterministic scalars $a_0=1/n$ , $b_0=2$ . At step $k$ it computes $\gamma_k\ge1/n$ from $\gamma_k^2-\gamma_k/n=(1-\gamma_k\sigma/n)a_k^2/b_k^2$ , sets $\alpha_k=\frac{n-\gamma_k\sigma}{\gamma_k(n^2-\sigma)}$ and $\beta_k=1-\gamma_k\sigma/n$ , forms $y_k=\alpha_kv_k+(1-\alpha_k)x_k$ , draws a block $i_k$ uniformly at random, and updates

x_{k+1}=T_{i_k}(y_k),\qquad v_{k+1}=\beta_kv_k+(1-\beta_k)y_k-\frac{\gamma_k}{L_{i_k}}U_{i_k}f'_{i_k}(y_k)^\#,

$b_{k+1}=b_k/\sqrt{\beta_k}$ , $a_{k+1}=\gamma_kb_{k+1}$ . The quantity $\phi_k=\mathbb E\,f(x_k)$ is the expectation over the first $k$ draws.

Formalization targets

Goal: Theorem 6, (5.3)

For every $k\ge0$ , with $D=2\|x_0-x^*\|_1^2+\frac1{n^2}(f(x_0)-f^*)$ ,

\phi_k-f^*\le\Big(\frac n{k+1}\Big)^2D,

and, when $\sigma>0$ ,

\phi_k-f^*\le\sigma D\Big[\Big(1+\tfrac{\sqrt\sigma}{2n}\Big)^{k+1}-\Big(1-\tfrac{\sqrt\sigma}{2n}\Big)^{k+1}\Big]^{-2}\le\Big(\frac n{k+1}\Big)^2D .

Milestones

(5.1), well-definedness: $\gamma_k\ge1/n$ is the unique such root, $0<\alpha_k,\beta_k\le1$ .
(5.2): $\gamma_k^2-\gamma_k/n=\beta_ka_k^2/b_k^2=(\beta_k\gamma_k/n)(1-\alpha_k)/\alpha_k$ .
The potential inequality: $2a_{k+1}^2(\phi_{k+1}-f^*)+b_{k+1}^2\mathbb E r_{k+1}^2\le2a_k^2(\phi_k-f^*)+b_k^2\mathbb E r_k^2$ , $r_k=\|v_k-x^*\|_1$ .
(5.4) $b_{k+1}\ge b_k+\frac\sigma{2n}a_k$ and (5.5) $a_{k+1}\ge a_k+\frac1{2n}b_k$ .
$a_k\ge\frac1{\sqrt\sigma}[Q_1^{k+1}-Q_2^{k+1}]$ and $b_k\ge Q_1^{k+1}+Q_2^{k+1}$ , $Q_{1,2}=1\pm\frac{\sqrt\sigma}{2n}$ .
$(1+t)^k-(1-t)^k\ge2kt$ for $t\ge0$ , hence $Q_1^{k+1}-Q_2^{k+1}\ge\frac{k+1}n\sqrt\sigma$ .

Significance

The result. Theorem 6 shows that randomized coordinate descent can be accelerated: the expected error falls as $n^2/k^2$ instead of the $n/k$ of the plain method, so after $k=n\cdot m$ iterations (about $m$ full-gradient equivalents) the error is $O(1/m^2)$ , the rate of the accelerated full-gradient method. For $\sigma>0$ the same scheme converges linearly with a ratio governed by $\sqrt\sigma/n$ rather than $\sigma/n$ . This is the starting point of the accelerated randomized coordinate methods (APPROX, NU-ACDM, accelerated proximal coordinate gradient) now used in large-scale machine learning.

Formalizing it. The theorem is proved in the paper; to our knowledge neither it nor any accelerated coordinate method has a machine-checked proof. The Lean development pins down a scheme whose correctness rests on a delicate parameter identity and on several conventions left implicit in the text (the range of $\sigma$ , the reading of (5.3) at $\sigma=0$ ).

Difficulty

The plain coordinate descent analysis tracks $f(x_k)$ alone and fails to accelerate: one coordinate step gains only $\frac1n$ of a full gradient step, and a momentum term built naively from coordinate steps does not keep the expected progress. The accelerated proof needs a potential mixing the function gap with a squared distance $\|v_k-x^*\|_1^2$ , weighted by the deterministic sequences $a_k^2$ and $b_k^2$ ; the weights must make the cross terms in $\|y_k-x^*\|_1^2$ and $f(y_k)$ cancel exactly after taking the expectation in the drawn block, which is what the coupled choice of $\gamma_k,\alpha_k,\beta_k$ achieves. The second difficulty is quantitative: the growth of $a_k$ follows from coupled recursions for $a_k$ and $b_k$ , not from a closed form.

Formalization scope

Space. $\mathbb R^N$ is Blocks E, the dependent product of finite-dimensional real inner-product spaces E i over Fin n (indices $0,\dots,n-1$ ); $B_i$ is absorbed into the inner product of E i. The dual norm is the operator norm of E i →L[ℝ] ℝ.
Randomness. Draws are explicit sequences Fin k → Fin n; $\phi_k$ and $\mathbb E r_k^2$ are finite averages over all $n^k$ sequences with weight $n^{-k}$ . No measure theory is needed.
$s^\#$ . Theorems hold for every selection satisfying (1.8), not a fixed choice.
Explicit hypotheses. The paper's standing assumptions become hypotheses: $n\ge1$ ; $L_i>0$ ; $f$ differentiable with (2.2); a minimizer $x^*$ exists. The convexity parameter $\sigma$ is any valid one (not necessarily the largest) with $0\le\sigma<n^2$ . The paper states $\sigma\ge0$ only; $\sigma<n^2$ is needed because $\alpha_k$ divides by $n^2-\sigma$ , and by footnote 2 ( $\sigma\le1$ ) it excludes only $n=1$ , $\sigma=1$ .
The case $\sigma=0$ of (5.3). The paper's middle expression is $0\cdot[0]^{-2}$ at $\sigma=0$ , read as a limit. Lean would evaluate it to $0$ , which would turn the statement into the false claim $\phi_k\le f^*$ ; the chain is therefore stated for $\sigma>0$ , and the outer bound $(n/(k+1))^2D$ separately for all $\sigma\ge0$ . Likewise the bound on $a_k$ involving $1/\sqrt\sigma$ is stated for $\sigma>0$ .
The coefficient $\gamma_k$ is the closed-form larger root of the quadratic in step 1; a milestone certifies it is the unique root $\ge1/n$ .
No trivialization. A formalization in which $\sigma=0$ is allowed in the middle bound, or in which $\phi_k$ is computed with a draw distribution other than uniform, or in which the coefficients depend on the draws, is not this theorem.
Corrections of the printed text. The update of $v_{k+1}$ prints $f'_i(y_k)^\#$ ; the index is $i_k$ . The middle term $b_k^2r_k^2$ of the potential inequality, written after the expectation in $\xi_{k-1}$ , is $b_k^2\mathbb E r_k^2$ .
Reusable parts. The block encoding (Blocks, partialGrad, CoordLipschitz, coordStep, wnorm, expect) is shared verbatim with the other missions of this series. The coefficient lemmas (5.2), (5.4), (5.5) and the growth bounds are pure real analysis and welcome as standalone contributions; so is a Lean proof of the block descent inequality (2.3)–(2.4) for Euclidean blocks, which the potential inequality needs.

Selected references

Yu. Nesterov, Efficiency of coordinate descent methods on huge-scale optimization problems, CORE Discussion Paper 2010/2, Université catholique de Louvain, 2010. Journal version: SIAM J. Optim. 22(2) (2012) 341–362. https://doi.org/10.1137/100802001
Yu. Nesterov, A method for solving the convex programming problem with convergence rate O(1/k²), Soviet Math. Dokl. 27 (1983) 372–376.
Y. T. Lee, A. Sidford, Efficient accelerated coordinate descent methods and faster algorithms for solving linear systems, FOCS 2013. https://arxiv.org/abs/1305.1922
O. Fercoq, P. Richtárik, Accelerated, parallel and proximal coordinate descent, SIAM J. Optim. 25(4) (2015) 1997–2023. https://arxiv.org/abs/1312.5799
Z. Allen-Zhu, Z. Qu, P. Richtárik, Y. Yuan, Even faster accelerated coordinate descent using non-uniform sampling, ICML 2016. https://arxiv.org/abs/1512.09103
S. Bubeck, Convex optimization: algorithms and complexity, Found. Trends Mach. Learn. 8 (2015), §3.7 (Nesterov's accelerated gradient descent, formalized on Prove2Me as ConvexOptAlg.NesterovSmooth.theorem_3_19). https://arxiv.org/abs/1405.4980

11 thms2 active usersReviewed

Completed

Optimization·Captain: mikedeng1

Sparsity Constrained Nonlinear Optimization: Optimality Conditions and Algorithms I: Every Coordinate-Wise Minimum Is an L2(f)-Stationary PointResearch Paper

Motivation

Many problems in signal processing, statistics and machine learning ask for a minimizer of a smooth loss among vectors with few nonzero entries: compressed sensing with a nonlinear measurement model, sparse regression with a non-quadratic loss, sparse phase retrieval. The common abstraction is the sparsity constrained problem

(\mathrm P)\qquad \min\ f(x)\quad\text{s.t.}\quad \|x\|_0\le s,

where $f$ is continuously differentiable but not necessarily convex. The constraint set is a finite union of linear subspaces, so (P) is nonconvex and in general NP-hard, and the classical optimality theory (KKT conditions, convex duality) gives no usable characterization of its solutions.

Beck and Eldar (SIAM J. Optim. 2013, arXiv:1203.4580) replaced that missing theory with a hierarchy of necessary optimality conditions of increasing strength, and matched each condition with an algorithm whose limit points satisfy it. This mission formalizes the comparison between the two strongest conditions of that hierarchy: coordinate-wise minimality and $L_2(f)$ -stationarity.

Setting

Vectors live in $\mathbb R^n$ with the Euclidean norm. For $x\in\mathbb R^n$ :

the $\ell_0$ norm $\|x\|_0$ is the number of nonzero components of $x$ ;
the support is $I_1(x)=\{i: x_i\ne0\}$ , and $I_0(x)=\{i: x_i=0\}$ is its complement;
$C_s=\{x:\|x\|_0\le s\}$ is the feasible set of (P), for an integer $0<s<n$ ;
$M_s(x)$ is the $s$ -th largest of $|x_1|,\dots,|x_n|$ , counted with multiplicity, so that $M_s(x)=0$ exactly when $\|x\|_0<s$ .

The objective $f:\mathbb R^n\to\mathbb R$ is continuously differentiable and bounded below (Assumption 1). Assumption 2 asks that $\nabla f$ be Lipschitz with some constant $L(f)$ .

The local Lipschitz condition (2.15) with constant $L_2$ asks that for every pair $i\ne j$ , every $x$ and every $d$ supported in $\{i,j\}$ ,

\big\|\nabla_{i,j}f(x)-\nabla_{i,j}f(x+d)\big\|\le L_2\|d\|,

where $\nabla_{i,j}f(x)=(\nabla_if(x),\nabla_jf(x))$ . The paper's local Lipschitz constant $L_2(f)=\max_{i\ne j}L_{i,j}(f)$ satisfies it, and $L_2(f)\le L(f)$ , often by a large factor.

A point $x^*\in C_s$ is a basic feasible (BF) vector (Definition 2.1) if $\nabla f(x^*)=0$ when $\|x^*\|_0<s$ , and $\nabla_if(x^*)=0$ on $I_1(x^*)$ when $\|x^*\|_0=s$ .

A point $x^*\in C_s$ is a coordinate-wise (CW) minimum (Definition 2.4) if either $\|x^*\|_0<s$ and $f(x^*)\le f(x^*+te_i)$ for every $i$ and every $t\in\mathbb R$ (Case I), or $\|x^*\|_0=s$ and

f(x^*)\le f(x^*-x^*_ie_i+te_j)\qquad\text{for every } i\in I_1(x^*),\ j\in\{1,\dots,n\},\ t\in\mathbb R

(Case II): no swap of a support coordinate for any coordinate, with any value, decreases $f$ .

Formalization targets

Goal: Theorem 2.4

Every CW-minimum $x^*$ of (P) satisfies

|\nabla_i f(x^*)|\ \begin{cases}\le L_2(f)\,M_s(x^*), & i\in I_0(x^*),\\ =0, & i\in I_1(x^*),\end{cases}

that is, $x^*$ is an $L_2(f)$ -stationary point. The constant $L_2(f)$ multiplies $M_s(x^*)$ exactly.

Milestones

Lemma 2.5. Every CW-minimum is a BF vector.
Lemma 2.6 (Local Descent Lemma). Under (2.15), for every $d$ with at most two nonzero components,

f(x+d)\le f(x)+\nabla f(x)^Td+\frac{L_2(f)}{2}\|d\|^2 .

(2.18). At a CW-minimum with $\|x^*\|_0=s$ , for $i\in I_0(x^*)$ , an index $m$ with $|x^*_m|=M_s(x^*)$ and $\sigma=\operatorname{sgn}(x^*_m\nabla_if(x^*))$ ,

f(x^*-x^*_me_m-\sigma x^*_me_i)\le f(x^*)-\sigma x^*_m\nabla_if(x^*)+L_2(f)(x^*_m)^2 .

Significance

Optimality conditions for (P) matter because no efficient algorithm finds a global minimizer; what algorithms can guarantee is a limit point satisfying some necessary condition, and the strength of that condition is the guarantee. The paper's hierarchy runs: optimal $\Rightarrow$ CW-minimum $\Rightarrow$ $L_2(f)$ -stationary $\Rightarrow$ $L$ -stationary for $L\ge L(f)$ $\Rightarrow$ BF vector. Theorem 2.4 is the step that ties the coordinate-wise notion, which needs no Lipschitz constant to state, to the stationarity notion behind iterative hard thresholding. Combined with "every optimal solution is a CW-minimum", it shows that every optimal solution of (P) is $L_2(f)$ -stationary, with the local constant instead of $L(f)$ . It is the reason the greedy and partial sparse-simplex methods of the same paper, whose limit points are CW-minima or $L_2(f)$ -stationary, are at least as strong as iterative hard thresholding.

The result is proved on paper. To our knowledge none of these notions (BF vectors, CW-minima, the local Lipschitz constant, the local descent lemma) has a machine-checked formalization. The mission produces a formal statement of the hierarchy's central step together with reusable definitions of the sparse-optimization objects, which the companion missions on the same paper (iterative hard thresholding, the greedy and partial sparse-simplex methods) state their results against.

Difficulty

The local descent lemma has no proof in the paper ("it is not difficult to see"). The global descent lemma integrates the gradient along the segment $[x,x+d]$ and applies the global Lipschitz bound; the local version must instead observe that along a direction supported on two coordinates only the two corresponding gradient components enter, and that the restricted Lipschitz bound (2.15) controls exactly those. Formally this requires the integral form of the mean value theorem for $C^1$ functions on a Euclidean space, together with the bookkeeping that $x+\tau d$ stays in the class of perturbations allowed by (2.15).

In Theorem 2.4 the natural first attempt, perturbing a single zero coordinate $i$ and comparing with $f(x^*)$ , fails when $\|x^*\|_0=s$ : the perturbed point has $s+1$ nonzeros and is infeasible, and Case II of the definition only compares $x^*$ with points obtained by swapping a support coordinate out. The bound on $|\nabla_if(x^*)|$ must therefore be extracted from swaps alone, and the degenerate case $\nabla_if(x^*)=0$ breaks the printed identity $\|x^*_me_m+\sigma x^*_me_i\|^2=2(x^*_m)^2$ in the paper's chain (2.18).

Formalization scope

All statements are in the namespace SparseNLO.CWOpt and work on EuclideanSpace ℝ (Fin n); indices run over Fin n, so $0,\dots,n-1$ instead of $1,\dots,n$ . The gradient is Mathlib's gradient, $te_i$ is EuclideanSpace.single i t, and $f$ is ContDiff ℝ 1. Assumption 1 is carried as a hypothesis on every statement, as the paper's standing assumption; Assumption 2 is carried wherever the page states it, with a free constant $L_f$ .

$M_s(x)$ is the entry at position $s-1$ of the list of $|x_i|$ sorted in nonincreasing order; the hypotheses $0<s<n$ of (P) make it meaningful. The local Lipschitz constant is not computed as a maximum: $L_2$ is any real constant satisfying (2.15), with "at most two nonzero components" read as "supported in some $\{i,j\}$ , $i\neq j$ " (the reading of Example 2.1). Every statement therefore holds for every valid $L_2$ , in particular for $L_2(f)$ . The minima in (2.11)–(2.12) are written as " $f(x^*)\le$ every value", which is equivalent. The goal states the displayed condition (2.16); its identification with Definition 2.3's $L$ -stationarity goes through Lemma 2.2, which belongs to a companion mission.

A trivializing formalization is ruled out: the CW-minimum predicate requires feasibility in $C_s$ and is satisfiable (for $n=2$ , $s=1$ , $f(x)=(x_1-1)^2+x_2^2$ , the point $(1,0)$ is a CW-minimum), $M_s$ is the $s$ -th largest absolute value, and $L_2$ is constrained by (2.15) rather than left free.

Contributions welcome: proofs of Lemma 2.6 (reusable for any coordinate-descent analysis on two coordinates), Lemma 2.5, (2.18) and the goal; general lemmas on $M_s$ (it is attained by some coordinate, and it vanishes exactly when $\|x\|_0<s$ ).

Theorem numbers and pages follow source.pdf, the arXiv version arXiv:1203.4580v1, not the SIAM typesetting.

Selected references

A. Beck, Y. C. Eldar, Sparsity Constrained Nonlinear Optimization: Optimality Conditions and Algorithms, SIAM Journal on Optimization 23(3), 2013. https://doi.org/10.1137/120869778 — arXiv:1203.4580v1 (2012), https://arxiv.org/abs/1203.4580
T. Blumensath, M. E. Davies, Iterative Hard Thresholding for Compressed Sensing, Applied and Computational Harmonic Analysis 27(3), 2009. https://doi.org/10.1016/j.acha.2009.04.002

5 thms2 active usersReviewed

Convex Optimization Operations Research Optimization·Captain: mikedeng1

Efficiency of Coordinate Descent Methods on Huge-Scale Optimization Problems 1: Random Block Coordinate Descent RCDM(α, x₀) Has Expected Error at Most 2·S_α·R²_{1−α}(x₀)/(k + 4)Research Paper

Motivation

Many large optimization problems arising in machine learning, statistics and network analysis have so many variables that computing a single full gradient is already expensive, while a single partial derivative, or a gradient with respect to a small block of variables, is cheap. Coordinate descent methods exploit this: at each iteration they update one block of variables only. They are among the oldest methods of numerical optimization, but for a long time their worst-case efficiency was not understood, because deterministic rules for choosing the next coordinate (cyclic order, greedy choice) are hard to analyse globally.

Nesterov (2010/2012) showed that choosing the coordinate at random, with probabilities tied to the coordinate-wise Lipschitz constants of the gradient, yields clean global complexity bounds. This paper started the modern theory of randomized coordinate descent, which was then extended to composite objectives, parallel and accelerated variants, and became a standard tool for huge-scale problems. This mission formalizes the first of its results: the expected sublinear rate of the basic method RCDM $(\alpha,x_0)$ (Theorem 1 of the CORE Discussion Paper 2010/2).

Setting

The variable $x\in\mathbb R^N$ is split into $n\ge1$ blocks, $\mathbb R^N=\mathbb R^{n_1}\times\cdots\times\mathbb R^{n_n}$ . Write $x^{(i)}\in\mathbb R^{n_i}$ for the $i$ -th block and $U_ih$ for the point whose $i$ -th block is $h$ and whose other blocks vanish. Each block space carries a norm $\|\cdot\|_{(i)}$ with dual norm $\|s\|^*_{(i)}=\max_{\|h\|_{(i)}=1}\langle s,h\rangle$ .

The objective $f:\mathbb R^N\to\mathbb R$ is convex and differentiable, and its set $X_*$ of minimizers is nonempty and bounded; $f^*$ is its optimal value. The partial gradient $f'_i(x)=U_i^T\nabla f(x)$ is the $i$ -th block of the gradient. The gradient is coordinate-wise Lipschitz with constants $L_i>0$ :

\|f'_i(x+U_ih_i)-f'_i(x)\|^*_{(i)}\le L_i\|h_i\|_{(i)}\qquad(2.2)

for all $x$ , $i$ and $h_i\in\mathbb R^{n_i}$ .

For a linear functional $s$ , $s^\#$ is any maximizer of $\langle s,x\rangle-\frac12\|x\|^2$ . The optimal coordinate step is $T_i(x)=x-\frac1{L_i}U_if'_i(x)^\#$ . For $\alpha\in\mathbb R$ let $S_\alpha=\sum_{i=1}^nL_i^\alpha$ and $p_\alpha^{(i)}=L_i^\alpha/S_\alpha$ . The method RCDM $(\alpha,x_0)$ starts at $x_0$ and, for $k\ge0$ , draws $i_k$ independently with $\Pr(i_k=i)=p^{(i)}_\alpha$ and sets $x_{k+1}=T_{i_k}(x_k)$ . Its expected objective value is $\varphi_k=E_{\xi_{k-1}}f(x_k)$ , the expectation over the draws $\xi_{k-1}=(i_0,\dots,i_{k-1})$ .

The weighted norms are $\|x\|_\beta=\big[\sum_iL_i^\beta\|x^{(i)}\|^2_{(i)}\big]^{1/2}$ and $\|g\|^*_\beta=\big[\sum_iL_i^{-\beta}(\|g^{(i)}\|^*_{(i)})^2\big]^{1/2}$ , and the size of the initial level set is

R_\beta(x_0)=\max_x\Big\{\max_{x_*\in X_*}\|x-x_*\|_\beta:\ f(x)\le f(x_0)\Big\}.

Formalization targets

Goal: Theorem 1

For every $k\ge0$ ,

\varphi_k-f^*\le\frac2{k+4}\cdot\Big[\sum_{j=1}^nL_j^\alpha\Big]\cdot R^2_{1-\alpha}(x_0).

The statement holds for every real $\alpha$ , every choice of the vectors $s^\#$ , and every block decomposition with arbitrary block norms. With $\alpha=0$ (uniform sampling) it reads $\varphi_k-f^*\le\frac{2n}{k+4}R_1^2(x_0)$ .

Milestones

$\|s^\#\|=\|s\|_*$ (§1, after (1.8)).
The block descent inequality (2.3).
The guaranteed decrease of an optimal coordinate step, $f(x)-f(T_i(x))\ge\frac1{2L_i}(\|f'_i(x)\|^*_{(i)})^2$ (2.4).
Lemma 2: the weighted Lipschitz bound $\|\nabla f(x)-\nabla f(y)\|^*_{1-\alpha}\le S_\alpha\|x-y\|_{1-\alpha}$ (2.9) and the quadratic upper bound (2.10).
The expected one-step decrease (2.13).
The recursion $\varphi_k-\varphi_{k+1}\ge\frac1C(\varphi_k-f^*)^2$ with $C=2S_\alpha R^2_{1-\alpha}(x_0)$ .
Lemma 1, a block-diagonal Loewner bound for positive semidefinite matrices, which the paper states in §1 but does not use later.

Significance

Theorem 1 is the first global efficiency estimate for a coordinate descent method on general smooth convex functions. Its constant is governed by $S_\alpha$ , an average of the coordinate Lipschitz constants, rather than by the global Lipschitz constant of the gradient; for $\alpha=1$ and Euclidean blocks this makes the method competitive with the full-gradient method even when partial derivatives are not cheap, and much faster when they are. The later results of the paper (linear convergence under strong convexity, high-probability bounds, the constrained and accelerated variants, adaptive Lipschitz estimates) reuse the objects and the one-step estimates of this mission.

The result is proved in the paper; nothing in it is open. As far as we know it has not been machine-checked. The platform contains a related item, ConvexOptAlg.CoordDescent.theorem_6_7 (Bubeck's monograph, Theorem 6.7, open), which treats scalar Euclidean coordinates, $\alpha\ge0$ , and the weaker factor $2/(t-1)$ ; together with its proved one-step lemmas (thm_6_7_coord_step, thm_6_7_expected_decrease) it covers the special case $n_i=1$ of milestones 3 and 5. This mission states the result in the paper's generality: arbitrary blocks and block norms, non-Euclidean dual norms through $s^\#$ , every real $\alpha$ , and the constant $2/(k+4)$ .

Difficulty

The bound is a statement about an expectation, not about individual runs: the recursion on $\varphi_k$ needs Jensen's inequality, $E[(f(x_k)-f^*)^2]\ge(\varphi_k-f^*)^2$ , and the fact that every run stays in the initial level set, so that $R_{1-\alpha}(x_0)$ controls $f(x_k)-f^*$ along every path. Lemma 2 is the step that links the coordinate-wise condition (2.2) to the weighted full gradient, and it is easy to misjudge: it is false without convexity, since $f(x)=x_1x_2$ on $\mathbb R^2$ satisfies (2.2) with arbitrarily small constants. The non-Euclidean block norms add a layer of convex analysis (the dual norm, the vector $s^\#$ and its identities) on top of the probabilistic bookkeeping.

Formalization scope

$\mathbb R^N$ is the dependent product Blocks E $=\prod_{i:\,\texttt{Fin }n}E_i$ of finite-dimensional real normed spaces; indices are 0-based. The norm Lean puts on the product is never used; statements use only the block norms and the weighted norms (2.7).
The dual norm is the operator norm of E i →L[ℝ] ℝ; the partial gradient is the derivative of $f$ composed with the inclusion of block $i$ .
$s^\#$ is an arbitrary selection satisfying (1.8) (IsSharpSelection); all results quantify over every selection.
Random draws are explicit index sequences, so $\varphi_k$ is a finite sum over $\{0,\dots,n-1\}^k$ weighted by products of the probabilities (2.5). No measure theory is involved.
$R_{1-\alpha}(x_0)$ is never computed as a supremum: the statements assume an upper bound $R$ (every point of the level set is within $R$ of every minimizer), which is equivalent when the maximum is finite.
Explicit hypotheses the paper leaves implicit: $n\ge1$ , $L_i>0$ , and convexity of $f$ in Lemma 2 (the standing assumption of §2, used in its proof). $\alpha$ is any real number; the remark " $S_\alpha$ … with $\alpha\ge0$ " on p. 6 is notation and is not imposed.
The recursion of milestone 6 is stated multiplied by $C$ , so that no statement divides by a quantity that may vanish.
A trivializing formalization is ruled out: replacing $\varphi_k$ by $f(x_k)$ along a single path, taking $R_{1-\alpha}(x_0)$ as a supremum that defaults to $0$ on unbounded level sets, or fixing one particular $s^\#$ would each change the theorem, and none is used. A sorry-free check confirms that all hypotheses of the goal hold for $n=1$ , $f(x)=x^2/2$ , $x_0=1$ , $R=1$ .

Contributions welcome: proofs of the one-step estimates (2.3), (2.4), (2.13), which need a block mean-value argument and the convex analysis of $s^\#$ ; Lemma 2, which needs the convex conjugate argument of its proof; and the final summation argument. The block encoding and the identities for $s^\#$ are reusable for every other mission of this paper.

Selected references

Yu. Nesterov, Efficiency of coordinate descent methods on huge-scale optimization problems, CORE Discussion Paper 2010/2, Université catholique de Louvain, 2010. https://core.ac.uk/download/6430808.pdf ; journal version: SIAM J. Optim. 22(2) (2012) 341–362, https://doi.org/10.1137/100802001
S. Bubeck, Convex Optimization: Algorithms and Complexity, Foundations and Trends in Machine Learning 8(3–4), 2015, §6.4. https://arxiv.org/abs/1405.4980
Yu. Nesterov, Introductory Lectures on Convex Optimization: A Basic Course, Kluwer, 2004, §2.1 (the descent lemma behind (2.3)). https://doi.org/10.1007/978-1-4419-8853-9

11 thms2 active usersReviewed

Probability Stochastic Systems·Captain: mikedeng1

Strong Convergence of an Explicit Numerical Method for SDEs with Nonglobally Lipschitz Continuous Coefficients: The Tamed Euler Scheme Converges to the SDE Solution in Uniform Lp with Order 1/2Research Paper

Motivation

Stochastic differential equations (SDEs) are simulated almost always by time-stepping schemes, and the cheapest and most widely used is the explicit Euler (Euler–Maruyama) scheme. Its strong convergence with order $\tfrac12$ is classical when the coefficients are globally Lipschitz. Many models in applications are not: Langevin dynamics with confining potentials, stochastic Ginzburg–Landau equations, and population or volatility models have drifts such as $\mu(x)=x-x^3$ that grow superlinearly.

For such drifts the explicit Euler scheme fails. Hutzenthaler, Jentzen and Kloeden showed in 2011 (Proc. R. Soc. A 467) that its absolute moments diverge to infinity whenever the drift or the diffusion grows superlinearly, so it does not converge in $L^p$ . Implicit schemes do converge (Higham, Mao and Stuart, SIAM J. Numer. Anal. 40, 2002), but each step requires solving a nonlinear equation.

The paper of this mission (Hutzenthaler, Jentzen and Kloeden, Ann. Appl. Probab. 22(4), 2012) introduces a tamed Euler scheme: an explicit method that costs the same as the Euler scheme, modified by a second-order term in the drift. It proves strong convergence with order $\tfrac12$ under a one-sided Lipschitz drift whose derivative grows at most polynomially. The scheme started a line of work on tamed and truncated methods, among them Sabanis's schemes that tame drift and diffusion together for superlinearly growing diffusion coefficients (Ann. Appl. Probab. 26(4), 2016), formalized on this platform in the series "Euler Approximations with Varying Coefficients".

Setting

Fix $T\in(0,\infty)$ and $d,m\in\mathbb N=\{1,2,\dots\}$ . Let $(\Omega,\mathcal F,\mathbb P)$ be a probability space with a filtration $(\mathcal F_t)$ , and let $W$ be an $m$ -dimensional standard $(\mathcal F_t)$ -Brownian motion. On $\mathbb R^k$ , $\|v\|$ is the Euclidean norm and $\langle v,w\rangle$ the inner product. For a matrix $A$ , $\|A\|=\sup_{\|v\|\le1}\|Av\|$ is the operator norm.

The initial value $\xi:\Omega\to\mathbb R^d$ is $\mathcal F_0$ -measurable with $\mathbb E\|\xi\|^p<\infty$ for every $p\ge1$ . The drift $\mu:\mathbb R^d\to\mathbb R^d$ is continuously differentiable and the diffusion is $\sigma:\mathbb R^d\to\mathbb R^{d\times m}$ . There is one constant $c>0$ with

\|\mu'(x)\|\le c(1+\|x\|^c),\qquad \|\sigma(x)-\sigma(y)\|\le c\|x-y\|,\qquad \langle x-y,\mu(x)-\mu(y)\rangle\le c\|x-y\|^2

for all $x,y$ . The exact solution $X$ is an adapted process with continuous paths and

X_t=\xi+\int_0^t\mu(X_s)\,ds+\int_0^t\sigma(X_s)\,dW_s,\qquad t\in[0,T],\ \mathbb P\text{-a.s.}

For $N\in\mathbb N$ , let $\Delta W^N_n=W_{(n+1)T/N}-W_{nT/N}$ . The tamed Euler scheme is $Y^N_0=\xi$ and

Y^N_{n+1}=Y^N_n+\frac{\tfrac TN\,\mu(Y^N_n)}{1+\tfrac TN\|\mu(Y^N_n)\|}+\sigma(Y^N_n)\,\Delta W^N_n,\qquad n\in\{0,\dots,N-1\}.

On $t\in[nT/N,(n+1)T/N]$ its interpolation is

\bar Y^N_t=Y^N_n+\frac{(t-nT/N)\,\mu(Y^N_n)}{1+\tfrac TN\|\mu(Y^N_n)\|}+\sigma(Y^N_n)(W_t-W_{nT/N}).

In Lean, states live in EthierKurtz.SDEState d, $\sigma$ takes values in Matrix (Fin d) (Fin m) ℝ, the standing hypotheses are the structures Coefficients and Setting, and the scheme is Y, Ybar, all in the namespace TamedEuler.Convergence.

Formalization targets

Goal: Theorem 1.1

There is a family $C_p\in[0,\infty)$ , $p\in[1,\infty)$ , such that for every solution $X$

\Big(\mathbb E\Big[\sup_{t\in[0,T]}\|X_t-\bar Y^N_t\|^p\Big]\Big)^{1/p}\le C_p\,N^{-1/2}\qquad\text{for all }N\in\mathbb N,\ p\in[1,\infty).

The constants are left unspecified; only the order $\tfrac12$ is asserted.

Milestones: Lemmas 3.1–3.10

These are the ten lemmas the proof rests on, in the paper's order:

the pathwise dominator lemma $\mathbb 1_{\Omega^N_n}\|Y^N_n\|\le D^N_n$ (3.1);
Gaussian and Brownian exponential moments (3.2, 3.3, 3.4);
moments of the dominating processes $D^N_n$ (3.5);
the decay of $\mathbb P[(\Omega^N_N)^c]$ faster than every power of $N$ (3.6);
continuous- and discrete-time Burkholder–Davis–Gundy inequalities with constant $p$ (3.7, 3.8);
the uniform moment bounds $\sup_N\sup_{n\le N}\mathbb E\|Y^N_n\|^p<\infty$ (3.9), and the same for $\mu(Y^N_n)$ and $\sigma(Y^N_n)$ (3.10).

Significance

The theorem gives an explicit scheme whose cost per step is the Euler cost and whose strong error is $O(N^{-1/2})$ in every $L^p$ , uniformly in time, for the superlinear drifts on which the explicit Euler scheme diverges. Strong error bounds of this kind are what multilevel Monte Carlo (Giles 2008) needs, as the paper notes on p. 3: its complexity analysis requires strong convergence of the coarse and fine levels, and the explicit Euler scheme supplies none for these equations. The moment bound of Lemma 3.9 answers, for the tamed scheme, the open problem Higham, Mao and Stuart formulated in 2002 on when moment bounds hold for explicit methods (quoted on p. 8).

The result is proved in the paper; it has not been formalized. The mission formalizes the statement, its ten lemmas and the scheme. The dominator construction (13)–(14) is a pathwise argument that can be checked independently of the probability. The two Burkholder–Davis–Gundy inequalities with the explicit constant $p$ are reusable tools for any numerical analysis of SDEs.

Difficulty

Once the moments of $Y^N_n$ are bounded uniformly in $N$ , the convergence proof follows the globally Lipschitz pattern: an Itô or Gronwall argument on $X-\bar Y^N$ . The difficulty is the moment bound itself. The obvious approach is a Gronwall recursion for $\mathbb E\|Y^N_{n+1}\|^p$ in terms of $\mathbb E\|Y^N_n\|^p$ . For superlinear $\mu$ it does not close: the term $\|\mu(Y^N_n)\|^2(T/N)^2$ grows faster than $\|Y^N_n\|^2$ , and for the untamed scheme the moments do in fact diverge. A newcomer's first idea, using the one-sided Lipschitz bound in expectation, fails for the same reason. The paper instead dominates the scheme pathwise on events $\Omega^N_n$ and controls the rare complement separately (Lemmas 3.1, 3.5, 3.6).

Formalization scope

The Lean development commits to the following conventions.

Normal filtration. Read as a right-continuous filtration (ℱ.IsRightContinuous) together with adaptedness to its $\mathbb P$ -completion, which is how the published SabanisEuler.Shared.IsSolution treats completeness.
Brownian motion. SabanisEuler.Shared.IsWienerMartingale, defined on $[0,\infty)$ rather than $[0,T]$ ; this is the usual extension. Time is ℝ≥0.
The solution $X$ . Any process satisfying the published IsSolution relation, whose stochastic integral is the Ethier–Kurtz Brownian Itô-integral relation. Existence and uniqueness are cited by the paper, not claimed. $\sigma$ enters IsSolution entrywise through toDiffusion. Every norm of $\sigma$ in a hypothesis or conclusion is the operator norm opNorm.
The constant $c$ . One real $c>0$ is both constant and exponent; $\|x\|^c$ and $N^{1/(2c)}$ are real powers.
Indices. $\mathbb N=\{1,2,\dots\}$ : every statement has $N\ge1$ and $n\le N$ , and $d,m\ge1$ .
Interpolation cell. In (10), $t$ uses the cell $n=\min(\lfloor tN/T\rfloor,N-1)$ , so $t=T$ uses the last cell.
(13) and (14). The supremum over $u$ in (13) is a finite maximum. The suprema in (14) are written "for all $k<n$ ", so $\Omega^N_0=\Omega$ .
Moments. Every expectation, $L^p$ norm and " $<\infty$ " is stated in $[0,\infty]$ with lintegral or eLpNorm, never as a Bochner integral that could default to $0$ .
Lemma 3.1. Stated pathwise for every $\omega$ and an arbitrary path $W$ .
Lemma 3.7. The published Itô relation demands pathwise square-integrability of the integrand. This is stronger than the printed almost-sure condition, and the deviation is disclosed.

A trivializing formalization is ruled out: the scheme is (8) exactly, with only the drift tamed and the step size $T/N$ . The error is taken against the same Brownian path that drives $X$ . The setting is satisfiable for the paper's example $\mu(x)=x-x^3$ , $\sigma(x)=x$ with $c=3$ (checked in a sorry-free verification file, given any Brownian motion). Replacing the Itô integral by a postulated operator or a pathwise integral is not admitted.

A complete development needs:

Gaussian exponential moments;
discrete and continuous BDG inequalities with explicit constants for the Ethier–Kurtz Itô integral;
an Itô-type or Gronwall argument for $X-\bar Y^N$ .

The BDG inequalities and the Gaussian moment identity are reusable beyond this mission. Contributions to any lemma, or to the stochastic-calculus infrastructure underneath, are welcome.

Selected references

M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients, Ann. Appl. Probab. 22(4), 1611–1641, 2012. https://arxiv.org/abs/1010.3756 (DOI 10.1214/11-AAP803)
M. Hutzenthaler, A. Jentzen, P. E. Kloeden, Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients, Proc. R. Soc. A 467, 1563–1576, 2011. https://doi.org/10.1098/rspa.2010.0348
D. J. Higham, X. Mao, A. M. Stuart, Strong convergence of Euler-type methods for nonlinear stochastic differential equations, SIAM J. Numer. Anal. 40(3), 1041–1063, 2002. https://doi.org/10.1137/S0036142901389530
S. Sabanis, Euler approximations with varying coefficients: the case of superlinearly growing diffusion coefficients, Ann. Appl. Probab. 26(4), 2083–2105, 2016. https://arxiv.org/abs/1308.1796
X. Mao, Stochastic Differential Equations and Their Applications, Horwood, 1997 (Theorem 2.4.1: existence and uniqueness, cited on p. 2). https://mathscinet.ams.org/mathscinet-getitem?mr=1475218
M. B. Giles, Multilevel Monte Carlo path simulation, Oper. Res. 56, 607–617, 2008. https://doi.org/10.1287/opre.1070.0496

18 thms2 active usersReviewed

Linear algebra Optimization·Captain: mikedeng1

The Primal-Dual Active Set Strategy as a Semismooth Newton Method IV: Global Convergence for Small Perturbations of M-MatricesResearch Paper

Motivation

Obstacle problems, contact problems and control-constrained optimal control problems lead, after discretization, to complementarity systems of the form "find $y \le \psi$ and $\lambda \ge 0$ with $Ay + \lambda = f$ and $\lambda_i(y_i - \psi_i) = 0$ ". The primal-dual active set strategy solves such systems by guessing which constraints are active, solving one linear system, and updating the guess. Each iteration costs a single linear solve, and in practice the method stops after few iterations.

Hintermüller, Ito and Kunisch (SIAM J. Optim. 13(3), 2002; authors' version hal-01660511v1) showed that the strategy is a semismooth Newton method for a reformulation of the complementarity system through a max-function, which gives local superlinear convergence. For convergence from arbitrary initial data the paper gives three sufficient conditions on the matrix. The first is that $A$ is an M-matrix (Theorem 3.2), a class that many discretized second-order elliptic operators belong to. The M-matrix property is not stable under small perturbations: off-diagonal entries may become positive (p. 10). Theorem 3.4, the subject of this mission, shows that global convergence nevertheless survives small perturbations of an M-matrix. The paper notes that this also makes it relevant for numerical implementations.

Setting

Fix $n \ge 0$ , vectors $f, \psi \in \mathbb{R}^n$ and a parameter $c > 0$ . Vectors are ordered componentwise. For a matrix $A \in \mathbb{R}^{n\times n}$ the problem is

Ay + \lambda = f, \qquad \lambda - \max\big(0, \lambda + c(y - \psi)\big) = 0, \tag{3.1}

with the maximum taken componentwise. The paper writes $A$ for the matrix and calligraphic $\mathcal{A}$ , $\mathcal{I}$ for the active and inactive index sets.

The primal-dual active set algorithm starts from an arbitrary pair $(y^0, \lambda^0)$ . Given $(y^k, \lambda^k)$ it sets

\mathcal{A}_k = \{ i : \lambda^k_i + c(y^k_i - \psi_i) > 0 \}, \qquad \mathcal{I}_k = \{ i : \lambda^k_i + c(y^k_i - \psi_i) \le 0 \},

and determines $(y^{k+1}, \lambda^{k+1})$ from

Ay^{k+1} + \lambda^{k+1} = f, \qquad y^{k+1} = \psi \text{ on } \mathcal{A}_k, \qquad \lambda^{k+1} = 0 \text{ on } \mathcal{I}_k .

For an index set $\mathcal{I} \subseteq \{1,\dots,n\}$ with complement $\mathcal{A}$ , $A_{\mathcal{I}}$ is the principal submatrix of $A$ on $\mathcal{I}$ and $A_{\mathcal{I}\mathcal{A}}$ the block with rows in $\mathcal{I}$ and columns in $\mathcal{A}$ . The step above is uniquely solvable exactly when $A_{\mathcal{I}_k}$ is invertible.

A square matrix $M$ is an M-matrix if it is nonsingular, $m_{ij} \le 0$ for $i \ne j$ , and $M^{-1} \ge 0$ entrywise. For a rectangular matrix $B$ , $\|B\|_1$ is the matrix norm subordinate to the one-norms on both spaces, i.e. the maximal absolute column sum $\|B\|_1 = \max_j \sum_i |b_{ij}|$ .

Formalization targets

Goal: Theorem 3.4

Let $M$ be an M-matrix. There is $\varepsilon > 0$ such that for every $K$ with $\|K\|_1 < \varepsilon$ , every $f, \psi$ and every $c > 0$ , with $A = M + K$ :

(3.1) has exactly one solution $x^* = (y^*, \lambda^*)$ ;
the algorithm is well-defined: from every $(y^k, \lambda^k)$ the step has exactly one solution;
for arbitrary initial data,

\lim_{k\to\infty} (y^k, \lambda^k) = (y^*, \lambda^*).

The threshold $\varepsilon$ is left unspecified; the theorem asserts only its existence.

Milestones (the steps of the paper's proof)

Principal submatrices $M_{\mathcal{I}}$ of an M-matrix are nonsingular M-matrices.
With $\rho \ge \|M_{\mathcal{I}}^{-1}K_{\mathcal{I}}\|_1$ for every $\mathcal{I}$ and $\rho < \tfrac12$ , every $A_{\mathcal{I}}$ is invertible and

A_{\mathcal{I}}^{-1} = \Big(I_{\mathcal{I}} + \sum_{i=1}^{\infty}\big(-M_{\mathcal{I}}^{-1}K_{\mathcal{I}}\big)^i\Big) M_{\mathcal{I}}^{-1}.

Under the same hypothesis, for $v \le 0$ , $\sum_{i\in\mathcal{I}} \big(-A_{\mathcal{I}}^{-1}v\big)_i \ge \frac{1-2\rho}{1-\rho}\,\|M_{\mathcal{I}}^{-1}v\|_1$ .
Whenever $A_{\mathcal{I}}$ is invertible, entrywise $A_{\mathcal{I}}^{-1}A_{\mathcal{I}\mathcal{A}} \le M_{\mathcal{I}}^{-1}K_{\mathcal{I}\mathcal{A}} - M_{\mathcal{I}}^{-1}K_{\mathcal{I}}(M+K)_{\mathcal{I}}^{-1}A_{\mathcal{I}\mathcal{A}}$ .
For $\|K\|_1$ small (uniformly in $f, \psi, c$ ) and $k \ge 1$ : $\sum_i y^{k+1}_i < \sum_i y^k_i$ unless $y^{k+1} = y^k$ , so $\mathcal{M}(y) = \sum_i y_i$ is a merit function.

A companion statement records the paper's remark after the theorem: for $\|K\|_1$ small, $M + K$ is a P-matrix.

Significance

The result. Theorem 3.2 makes global convergence depend on a sign pattern, which a perturbation of any size can destroy. Theorem 3.4 replaces the sign pattern by a neighbourhood: every matrix close enough to an M-matrix in $\|\cdot\|_1$ inherits unique solvability of (3.1), well-posedness of every iteration, and convergence from every starting point. As the paper observes, solvability of (3.1) for every $f$ is equivalent to $A$ being a P-matrix (Berman–Plemmons, Thm. 10.2.15), so the theorem also yields a sufficient condition for the P-matrix property.

Formalizing it. The paper's argument is a short sketch that refers back to the proof of Theorem 3.3 and asserts that " $\|K\|_1$ can be chosen sufficiently small". A formalization makes the uniformity of that choice explicit (it depends only on $M$ ), checks the strict decrease of the merit function including its equality case, and fixes the role of the first iterate, from which the sign facts used in (3.8) hold. To our knowledge none of the three global convergence theorems of the paper has a machine-checked proof.

Difficulty

The obvious argument is to transfer Theorem 3.2 by continuity. It fails because the M-matrix proof is monotonicity-based: it shows $y^* \le y^{k+1} \le y^k$ , and that ordering uses $A_{\mathcal{I}}^{-1} \ge 0$ and $A_{\mathcal{I}\mathcal{A}} \le 0$ , sign conditions that a perturbation breaks immediately. What survives is weaker: a scalar merit function $\sum_i y_i$ whose decrease must be quantified. The perturbation has to be controlled uniformly over all $2^n$ index sets, both in the multiplier term (where the decrease must remain strict whenever $\lambda^k_{\mathcal{I}} \ne 0$ ) and in the coupling term $A_{\mathcal{I}}^{-1}A_{\mathcal{I}\mathcal{A}}(y^k-\psi)_{\mathcal{A}}$ , whose positive part must be absorbed by $-\sum_{\mathcal{A}}(y^k_i - \psi_i)$ . Finally, strict decrease of a merit function does not by itself give convergence; the finiteness of the possible active sets must be combined with it.

Formalization scope

Vectors are Fin n → ℝ and matrices Matrix (Fin n) (Fin n) ℝ; index sets are Finset (Fin n) and blocks are submatrices along subtype inclusions. Conventions committed to:

$\|B\|_1$ is the maximal absolute column sum, computed as a supremum in $\mathbb{R}_{\ge0}$ , equal to $0$ for a matrix with no columns; it is not the entrywise sum or another operator norm;
the M-matrix predicate requires nonsingularity explicitly (IsUnit of the determinant), because Lean's inverse of a singular matrix is $0$ ; it is stated for any finite index type so that it applies to principal blocks;
" $\|K\|_1$ sufficiently small" is one $\varepsilon > 0$ depending only on $M$ (and $n$ ), quantified before $f$ , $\psi$ , $c$ and the initial data;
a step is a relation and a run is an infinite sequence of steps from arbitrary initial data; step (iv)'s "Stop" is not modelled;
the convergence conclusion is stated for every solution of (3.1), which by the first conclusion is the solution; convergence is in the product topology;
the merit-function milestone holds from $k = 1$ on, as in the paper.

The goal does not assume that $M + K$ is a P-matrix, that its principal blocks are invertible, that $\rho < \tfrac12$ , or that the merit function decreases: all of these are consequences. A version with $\varepsilon$ chosen after $f$ , $\psi$ , $c$ or the run, or with any of these as hypotheses, would be trivial or different and is ruled out. Well-definedness is a conclusion, so runs from every initial pair exist and the convergence conclusion is not vacuous.

In the g-estimate, the page's final $\|M^{-1}\lambda^k_{\mathcal{I}}\|_1$ is a misprint for $\|M_{\mathcal{I}}^{-1}\lambda^k_{\mathcal{I}}\|_1$ ; the statement uses $M_{\mathcal{I}}^{-1}$ .

A complete development needs: principal submatrices of M-matrices, Neumann series for matrices with the column-sum norm, the block form of one step and the identity (3.8), the sign facts $\lambda^k \le 0$ on $\mathcal{I}_k$ and $y^k \ge \psi$ on $\mathcal{A}_k$ for $k \ge 1$ , and continuity of principal minors for the P-matrix statement. The M-matrix and Neumann-series lemmas are reusable beyond this mission. Proofs of any milestone, of the companion, or of the goal through another argument are welcome.

Selected references

M. Hintermüller, K. Ito, K. Kunisch, The primal-dual active set strategy as a semismooth Newton method, SIAM J. Optim. 13(3):865–888, 2002. https://doi.org/10.1137/S1052623401383558 ; authors' version https://hal.science/hal-01660511v1
A. Berman, R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic Press, 1979 (SIAM Classics reprint 1994). https://doi.org/10.1137/1.9781611971262
R. H. W. Hoppe, Multigrid algorithms for variational inequalities, SIAM J. Numer. Anal. 24:1046–1065, 1987. https://doi.org/10.1137/0724069

10 thms2 active usersReviewed

Linear algebra Optimization·Captain: mikedeng1

The Primal-Dual Active Set Strategy as a Semismooth Newton Method II: Global Convergence and Monotonicity for M-MatricesResearch Paper

Motivation

Obstacle problems, contact problems and box-constrained quadratic programs share one finite-dimensional core: the linear complementarity problem "find $y \le \psi$ and $\lambda \ge 0$ with $Ay + \lambda = f$ and $\lambda_i(y_i - \psi_i) = 0$ ". After discretization by finite differences or piecewise linear finite elements, the matrix $A$ is typically an M-matrix, which is why this class sits at the centre of the numerical analysis of variational inequalities.

The primal-dual active set strategy solves such problems by guessing which constraints are active, solving one linear system on that guess, and updating the guess from the new primal and dual values. Hintermüller, Ito and Kunisch (SIAM J. Optim. 13 (2002) 865–888; authors' version hal-01660511) showed that this strategy is a semismooth Newton method for a nonsmooth reformulation of the problem. Local superlinear convergence then follows from Newton theory. This mission concerns the paper's global result: for M-matrices the method converges from every starting point, and its primal iterates decrease monotonically. A related global result for a purely primal active set method appears in Hoppe's 1987 analysis of multigrid methods for obstacle problems, which the paper cites as [H].

Setting

Fix $n \ge 0$ , a matrix $A \in \mathbb{R}^{n\times n}$ , vectors $f, \psi \in \mathbb{R}^n$ and a constant $c > 0$ . Vectors are compared componentwise. The paper writes $A$ for the matrix and calligraphic $\mathcal{A}$ , $\mathcal{I}$ for index sets; Lean uses A for the matrix and activeSet for $\mathcal{A}$ .

Problem (3.1). Find $(y, \lambda) \in \mathbb{R}^n \times \mathbb{R}^n$ with

Ay + \lambda = f, \qquad \lambda - \max(0, \lambda + c(y - \psi)) = 0,

the maximum taken componentwise. For $c > 0$ the second equation says exactly $y \le \psi$ , $\lambda \ge 0$ and $\lambda_i (y_i - \psi_i) = 0$ for every $i$ .

M-matrix. $A$ is an M-matrix if it is nonsingular, its off-diagonal entries satisfy $a_{ij} \le 0$ for $i \ne j$ , and its inverse is entrywise nonnegative, $A^{-1} \ge 0$ .

Blocks. For an index set $\mathcal{I}$ with complement $\mathcal{A}$ , $A_{\mathcal{I}}$ is the principal submatrix with rows and columns in $\mathcal{I}$ , and $A_{\mathcal{I}\mathcal{A}}$ the submatrix with rows in $\mathcal{I}$ and columns in $\mathcal{A}$ .

The algorithm. Start from arbitrary $(y^0, \lambda^0)$ . Given $(y^k, \lambda^k)$ , let

\mathcal{A}_k = \{ i : \lambda^k_i + c(y^k - \psi)_i > 0\}, \qquad \mathcal{I}_k = \{ i : \lambda^k_i + c(y^k - \psi)_i \le 0\},

and let $(y^{k+1}, \lambda^{k+1})$ solve

Ay^{k+1} + \lambda^{k+1} = f, \qquad y^{k+1} = \psi \text{ on } \mathcal{A}_k, \qquad \lambda^{k+1} = 0 \text{ on } \mathcal{I}_k.

A sequence obeying this rule for every $k$ is a run. The solution of (3.1) is written $x^* = (y^*, \lambda^*)$ and the iterates $x^k = (y^k, \lambda^k)$ .

Formalization targets

Goal: Theorem 3.2 (p. 7)

If $A$ is an M-matrix, then for arbitrary initial data

x^k \to x^*, \qquad y^* \le y^{k+1} \le y^k \ \text{ for all } k \ge 1, \qquad y^k \le \psi \ \text{ for all } k \ge 2.

Milestones (Appendix A, p. 20)

For every partition $\mathcal{I}, \mathcal{A}$ : $A_{\mathcal{I}}$ is regular, $A_{\mathcal{I}}^{-1} \ge 0$ and $A_{\mathcal{I}}^{-1} A_{\mathcal{I}\mathcal{A}} \le 0$ .
Complementarity (A.1): $\lambda^k_i = 0$ or $y^k_i = \psi_i$ , for all $i$ and $k \ge 1$ , for every run and every matrix.
Monotonicity: $y^{k+1} \le y^k$ for $k \ge 1$ .
Primal feasibility: $y^k \le \psi$ for $k \ge 2$ .
Lower bound: $y^* \le y^k$ for $k \ge 1$ .
Remark A.1 and dual feasibility: if $\lambda^{\bar k}_i < 0$ with $\bar k \ge 1$ , then $\lambda^k_i = 0$ and $y^k_i \le \psi_i$ for all $k \ge \bar k + 1$ ; and $\lambda^k \ge 0$ for all $k$ beyond some $k_o$ .

A companion item states that a run exists from every initial pair when $A$ is an M-matrix.

Significance

The result. Theorem 3.2 makes the active set strategy a globally convergent method for the discretized obstacle problem, with no line search, no merit function and no damping. The monotonicity statement says more: from the second iterate on, every primal iterate is feasible and lies above the solution, and the iterates descend to it. Combined with the paper's observation (p. 7) that convergent iterates of the method converge in finitely many steps, it gives finite termination from every starting point. The theorem is the model for later global results on the method under weaker hypotheses, two of which (a norm condition and small perturbations of M-matrices) appear in the same paper as Theorems 3.3 and 3.4.

Formalizing it. The theorem is proved in the paper; no machine-checked version is known. A formal proof requires the classical sign facts on principal blocks of M-matrices, which Mathlib does not contain, and a careful treatment of the iteration as a relation on pairs of vectors. The M-matrix block lemma is reusable well beyond this mission, for instance in comparison principles for discretized elliptic problems.

Difficulty

The obvious argument would bound a merit function or a contraction factor; neither is available, because the iteration map is discontinuous and the active sets change from step to step. The paper's argument is instead an order argument. Its central step is the monotonicity $y^{k+1} \le y^k$ : it requires the complementarity (A.1) of the previous iterate, which holds only from $k = 1$ on, and the sign structure of the blocks $A_{\mathcal{I}_k}^{-1}$ and $A_{\mathcal{I}_k}^{-1} A_{\mathcal{I}_k \mathcal{A}_k}$ for every possible active set. The index bounds are sharp in the sense that the statements fail at $k = 0$ (respectively $k = 1$ for feasibility), because the initial data are arbitrary. Convergence of the multipliers is a separate step: negative multiplier coordinates must be shown to vanish permanently after one more iteration.

Formalization scope

Vectors are Fin n → ℝ, matrices Matrix (Fin n) (Fin n) ℝ, the product is A *ᵥ y, and $\le$ on vectors is componentwise. Convergence of $x^k$ is in the product topology on $\mathbb{R}^n \times \mathbb{R}^n$ .
Nonsingularity of an M-matrix is IsUnit A.det. Lean's inverse of a singular matrix is the zero matrix, so without this clause every singular matrix with nonpositive off-diagonal entries would count as an M-matrix; the M-matrix definition is the page's, not an equivalent characterization.
The solution is a hypothesis IsSolution A f ψ c ystar lamstar. The paper speaks of "the unique solution"; for an M-matrix it exists and is unique by Berman–Plemmons, so the hypothesis assumes nothing the paper does not have. The uniqueness fact is not restated.
The algorithm is a relation. One step is IsStep, a run is IsRun with free initial data y 0, lam 0. The stopping option of step (iv) is not modelled. Because a universally quantified statement over runs would be empty if no run existed, the companion item run_exists states existence of a run from every initial pair.
The constant $c > 0$ is a binder of every statement that needs it; (A.1) needs neither $c > 0$ nor any hypothesis on $A$ , and has none.
Remark A.1 is printed without " $\bar k \ge 1$ "; its proof assumes it and the statement is false at $\bar k = 0$ , so the bound is kept.
A trivializing formalization is ruled out: no statement assumes monotonicity, feasibility or convergence of the iterates, and the M-matrix hypothesis cannot be met by a singular matrix.
Needed infrastructure: principal submatrices of M-matrices (Berman–Plemmons, Ch. 6), and inverse-nonnegativity arguments for block systems. Proofs of the block lemma, of (A.1) and of the order lemmas are welcome independently of the goal.

The printed closing step of the proof writes $\lambda^k = Ay^k + f$ ; step (iii) gives $\lambda^k = f - Ay^k$ . No item quotes or depends on that line.

Selected references

M. Hintermüller, K. Ito, K. Kunisch, The primal-dual active set strategy as a semismooth Newton method, SIAM J. Optim. 13(3):865–888, 2002. https://doi.org/10.1137/S1052623401383558 (authors' version: https://hal.science/hal-01660511)
A. Berman, R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic Press, 1979; reprinted as SIAM Classics in Applied Mathematics 9, 1994. https://doi.org/10.1137/1.9781611971262
R. H. W. Hoppe, Multigrid algorithms for variational inequalities, SIAM J. Numer. Anal. 24(5):1046–1065, 1987. https://doi.org/10.1137/0724069
L. Qi, J. Sun, A nonsmooth version of Newton's method, Math. Program. 58:353–367, 1993. https://doi.org/10.1007/BF01581275

9 thms2 active usersReviewed

Completed

Optimization·Captain: mikedeng1

A D.C. Optimization Algorithm for Solving the Trust-Region Subproblem IV: A Kuhn–Tucker Point with λ* + λ₁ < 0 Gives an Explicit Feasible Restart Point x̄ with f(x̄) < f(x*)Research Paper

Motivation

The trust-region subproblem — minimizing a quadratic, possibly nonconvex, over a Euclidean ball — is the inner step of every trust-region method for nonlinear optimization, and it also appears in regularized least squares and in the analysis of quadratic constraints. It is one of the few nonconvex problems whose global minimizers are characterized exactly: by the classical results of Moré and Sorensen (1983) and Gay (1981), a feasible point is a global solution if and only if it satisfies the Kuhn–Tucker conditions with a multiplier that makes the shifted matrix positive semidefinite.

Pham Dinh Tao and Le Thi Hoai An (SIAM J. Optim. 8 (1998)) solve the problem with DCA, a descent method for difference-of-convex programs. DCA converges to Kuhn–Tucker points, and a Kuhn–Tucker point need not be a global solution. Section 4.2 of the paper closes this gap: it tests a Kuhn–Tucker point for global optimality with one eigenvalue computation, and when the test fails it writes down, explicitly, a feasible point with strictly smaller objective value from which DCA can be restarted. The resulting algorithm GDCA alternates DCA runs and restarts. This mission formalizes that restart construction.

Setting

Let $A$ be a real symmetric $n\times n$ matrix, $b\in\mathbb R^n$ and $r>0$ . Write $\langle\cdot,\cdot\rangle$ for the Euclidean inner product and $\|\cdot\|$ for its norm. The trust-region subproblem (Q1) is

\min\Bigl\{\,f(x)=\tfrac12\langle x,Ax\rangle+\langle b,x\rangle \;:\; \|x\|\le r\,\Bigr\}.

A point $x^*$ is a Kuhn–Tucker point of (Q1) with multiplier $\lambda^*$ when

\lambda^*\ge0,\qquad (A+\lambda^*I)x^*=-b,\qquad \lambda^*(\|x^*\|-r)=0,\qquad \|x^*\|\le r .

In Lean, $\mathbb R^n$ is EuclideanSpace ℝ (Fin n), $A$ is a self-adjoint continuous linear operator on it, $f$ is quad A b, and the four conditions are IsKKT A b r xs lamStar.

Let $(\lambda_1,u)$ be an eigenpair of $A$ : $u\ne0$ and $Au=\lambda_1u$ . If $\lambda^*+\lambda_1<0$ then

u^T(A+\lambda^*I)u=(\lambda^*+\lambda_1)\,u^Tu<0,

so $u$ is a direction of negative curvature of the shifted matrix. For a direction $w\ne0$ and $\|x^*\|\le r$ , the restart step $\gamma$ of (23) is the nonzero root of $\|w\|^2\gamma^2+2w^Tx^*\gamma+\|x^*\|^2-r^2=0$ , so that $x^*+\gamma w$ lies on the sphere $\|x\|=r$ (Lean: restartGamma w xs r s, with a sign $s=\pm1$ choosing between the two roots when $\|x^*\|<r$ ).

Formalization targets

Goal: a Kuhn–Tucker point that fails the eigenvalue test can be improved

\text{IsKKT}(x^*,\lambda^*),\;\; Au=\lambda_1u,\;u\ne0,\;\; \lambda^*+\lambda_1<0 \;\Longrightarrow\; \exists\,\bar x:\ \|\bar x\|\le r,\ \ f(\bar x)<f(x^*).

This is the "Summing up" claim of §4.2 (p. 493). It asserts only the existence of a better feasible point; the explicit construction is spread over the milestones.

Milestones, in attack order

The multiplier (p. 491): $\lambda^*=(-\langle x^*,Ax^*\rangle-\langle x^*,b\rangle)/r^2$ .
Lemma 4.2 (p. 492): if $(A+\mu I)y=-b$ , then $f(x)=f(y)-\frac\mu2(\|x\|^2-\|y\|^2)+\frac12\langle x-y,(A+\mu I)(x-y)\rangle$ for all $x$ .
Case a) (p. 492): if $\langle b,x^*\rangle>0$ , then $\bar x=-x^*$ is feasible and $f(\bar x)=f(x^*)-2\langle b,x^*\rangle<f(x^*)$ .
Case (b.1) (p. 492): if $w^T(A+\lambda^*I)w<0$ and either $\|x^*\|<r$ , or $\|x^*\|=r$ and $w^Tx^*\ne0$ , then $\bar x=x^*+\gamma w$ with $\gamma$ from (23) satisfies $\gamma\ne0$ , $\|\bar x\|=r$ and $f(\bar x)=f(x^*)+\frac{\gamma^2}2\langle w,(A+\lambda^*I)w\rangle<f(x^*)$ .
Case (b.2), $b^Tx^*=0$ (p. 493): if $\|x^*\|=r$ , $u^Tx^*=0$ and $u^T(A+\lambda^*I)u<0$ , then $v=u+\tau x^*$ has $v^Tx^*\ne0$ and $v^T(A+\lambda^*I)v<0$ for all $\tau<0$ when $b^Tu\le0$ , and for $u^T(A+\lambda^*I)u/(2b^Tu)<\tau<0$ when $b^Tu>0$ .
Case (b.2), $b^Tx^*<0$ (p. 493): the same conclusion for $\tau_1<\tau<0$ , where $\tau_1<0$ is the smallest root of $-b^Tx^*\tau^2-2b^Tu\,\tau+u^T(A+\lambda^*I)u$ .

A companion item records the paper's remark that, for $\|x^*\|<r$ , the root of (23) with the larger $\gamma^2$ gives the smaller $f(\bar x)$ .

Significance

The restart claim turns a local method into a global one. Combined with the convergence of DCA to Kuhn–Tucker points (Theorem 4.1 of the paper) and the global-optimality characterization of (Q1), it shows that each DCA run either ends at a certified global solution or yields a feasible point that strictly improves on it, with the improvement given in closed form from one eigenvector. The construction is elementary, but its case analysis — the sign of $\langle b,x^*\rangle$ , whether the constraint is active, whether the eigenvector is orthogonal to $x^*$ — is exactly where an informal argument can go wrong, and the printed (b.2) bounds contain two errors (see Formalization scope).

The result is proved in the paper; to our knowledge none of it has a machine-checked proof. The mission produces a checked version of the case analysis with the corrected constants, and reusable facts about quadratics on the Euclidean ball (Lemma 4.2, the step to the sphere along a negative-curvature direction) that apply to any trust-region method.

Difficulty

The obvious argument — move from $x^*$ along the negative-curvature eigenvector $u$ until the boundary is reached — fails in one case. When $\|x^*\|=r$ and $u^Tx^*=0$ , the line $x^*+\gamma u$ leaves the ball for every $\gamma\ne0$ , so no step along $u$ is feasible. The paper repairs this by tilting $u$ towards $x^*$ , and the tilt must keep the curvature negative while making the direction non-orthogonal to $x^*$ ; the admissible range of the tilt depends on the signs of $b^Tx^*$ and $b^Tu$ . A second subtlety is that the identity $f(\bar x)=f(x^*)+\frac{\gamma^2}{2}\langle w,(A+\lambda^*I)w\rangle$ needs the term $\lambda^*(\|\bar x\|^2-\|x^*\|^2)$ of Lemma 4.2 to vanish, which uses complementarity when $\|x^*\|<r$ and $\|\bar x\|=r$ when the constraint is active.

Formalization scope

All statements live in namespace TaoAnDCA.Restart on EuclideanSpace ℝ (Fin n), with $A$ a continuous linear operator and symmetry the hypothesis IsSelfAdjoint A; $r>0$ throughout. The objective quad is the bare quadratic, with feasibility $\|x\|\le r$ stated separately.

Conventions and disclosed deviations from the page:

The goal assumes an arbitrary eigenpair with $\lambda^*+\lambda_1<0$ and an arbitrary Kuhn–Tucker point. The page has "the smallest eigenvalue" and "a point produced by the DCA"; neither is used, so the stated theorem is stronger.
Case (b.1) is stated for a general direction $w$ of negative curvature, since (b.2) re-enters it with a vector that is not an eigenvector; the heading " $\langle b,x^*\rangle\le0$ " of case b) is not assumed because the computation does not use it.
Cases (b.2) take a general $u$ with $u^T(A+\lambda^*I)u<0$ ; for an eigenvector $u$ with $u^Tx^*=0$ one has $b^Tu=0$ , which would make the $b^Tu>0$ branch vacuous.
The page's bound $u^T(A+\lambda^*I)u/b^Tu$ in the first (b.2) bullet lacks a factor 2 and is false for general $u$ ; the Lean uses $u^T(A+\lambda^*I)u/(2b^Tu)$ . The page's formula for $\tau_1$ has the sign of $b^Tu$ reversed; the Lean uses the smallest root, as the page's words define it. Both printed versions agree with the corrected ones when $b^Tu=0$ .

A trivializing formalization is ruled out: the goal's hypothesis is the eigenvalue test $\lambda^*+\lambda_1<0$ , never " $x^*$ is not a global solution", which would make the conclusion a restatement of the hypothesis.

No infrastructure beyond Mathlib's inner product spaces is needed. Contributions welcome: proofs of the milestones, and of the goal from them.

Selected references

Pham Dinh Tao and Le Thi Hoai An, A d.c. optimization algorithm for solving the trust-region subproblem, SIAM J. Optim. 8(2):476–505, 1998. https://doi.org/10.1137/S1052623494274313
J. J. Moré and D. C. Sorensen, Computing a trust region step, SIAM J. Sci. Stat. Comput. 4(3):553–572, 1983. https://doi.org/10.1137/0904038
D. M. Gay, Computing optimal locally constrained steps, SIAM J. Sci. Stat. Comput. 2(2):186–197, 1981. https://doi.org/10.1137/0902016

9 thms2 active usersReviewed

Convex Optimization Optimization·Captain: mikedeng1

A D.C. Optimization Algorithm for Solving the Trust-Region Subproblem III: The DCA Step x^{k+1} = P_E(x^k − (Ax^k + b)/ρ) Is a Descent Method Whose Limit Points Are Kuhn–Tucker PointsResearch Paper

Motivation

The trust-region subproblem asks for the minimum of a quadratic function over a Euclidean ball,

(Q_1)\qquad \min\Big\{\tfrac12 x^TAx + b^Tx : \|x\|\le r\Big\},

with $A$ a real symmetric $n\times n$ matrix that may be indefinite. It is solved at every iteration of a trust-region method for nonlinear optimization (Moré and Sorensen 1983; Conn, Gould and Toint 2000), so a method for it is used many thousands of times in a single run. Classical methods based on factorizations of $A + \lambda I$ do not scale to large sparse problems.

Pham Dinh Tao and Le Thi Hoai An (SIAM J. Optim. 8 (1998)) apply to $(Q_1)$ their DCA, a general algorithm for d.c. programs (minimizing a difference of two convex functions). For a suitable splitting of the objective, one DCA iteration is a matrix–vector product followed by a scaling onto the ball. This mission formalizes the convergence theorem for that iteration (Theorem 4.1 of the paper), together with the general DCA convergence results it is derived from (Theorem 3.7 and the Appendix).

Timeline (as recorded in the paper's introduction):

1978: J. F. Toland develops the duality theory of d.c. optimization (J. Math. Anal. Appl. 66 (1978)).
1988: Pham Dinh and Elbernoussi introduce the DCA based on this duality, in its subgradient form (reference [23] of the paper, Trends in Mathematical Optimization, Birkhäuser).
1994: Le Thi's thesis and the authors' Rouen research report develop the DCA's convergence theory and large-scale implementations (references [7], [27] of the paper).
1998: the present paper specializes the DCA to $(Q_1)$ , proves Theorem 4.1, and gives a restarting procedure that reaches a global solution.

Setting

Let $X = \mathbb R^n$ with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$ . For a function $\theta:X\to\mathbb R\cup\{+\infty\}$ , $\operatorname{dom}\theta = \{x:\theta(x)<+\infty\}$ , and $\Gamma_0(X)$ is the set of proper, lower semicontinuous, convex functions. For $\theta\in\Gamma_0(X)$ the conjugate is $\theta^*(y) = \sup_x\{\langle x,y\rangle-\theta(x)\}$ , and the subdifferential at $x\in\operatorname{dom}\theta$ is $\partial\theta(x) = \{y : \theta(z)\ge\theta(x)+\langle z-x,y\rangle\ \forall z\}$ . The function $\theta$ is $\rho$ -convex ( $\rho\ge0$ ) if $\theta-\frac\rho2\|\cdot\|^2$ is convex.

A d.c. program is $(P)\ \alpha = \inf\{g(x)-h(x): x\in X\}$ with $g,h\in\Gamma_0(X)$ , under the convention $+\infty-(+\infty)=+\infty$ and the standing inclusions $\operatorname{dom}g\subset\operatorname{dom}h$ , $\operatorname{dom}h^*\subset\operatorname{dom}g^*$ . Its dual objective is $h^*(y)-g^*(y)$ . The simplified DCA builds two sequences from $x^0\in\operatorname{dom}g$ :

y^k\in\partial h(x^k),\qquad x^{k+1}\in\partial g^*(y^k).

For $(Q_1)$ let $E = \{x:\|x\|\le r\}$ ( $r>0$ ), $\chi_E$ its indicator, and $f(x) = \frac12x^TAx+b^Tx+\chi_E(x)$ . Fix $\rho>0$ with $\rho I-A$ positive semidefinite and split $f = g-h$ as in (16):

g(x) = \tfrac12\rho\|x\|^2+b^Tx+\chi_E(x),\qquad h(x) = \tfrac12x^T(\rho I-A)x.

The resulting DCA starts at any $x^0\in\mathbb R^n$ and, with $w^k = (\rho I-A)x^k-b$ , sets $x^{k+1} = w^k/\rho$ if $\|w^k\|\le\rho r$ and $x^{k+1} = r\,w^k/\|w^k\|$ otherwise; equivalently $x^{k+1} = P_E\big(x^k-\frac1\rho(Ax^k+b)\big)$ . A Kuhn–Tucker point of $(Q_1)$ is an $x^*$ with some $\lambda^*\ge0$ such that $(A+\lambda^*I)x^* = -b$ , $\lambda^*(\|x^*\|-r)=0$ and $\|x^*\|\le r$ .

Formalization targets

Goal: Theorem 4.1 (p. 491)

For every run of the DCA,

f(x^{k+1})\le f(x^k)-\tfrac12(\rho+\lambda)\|x^{k+1}-x^k\|^2,\qquad \lambda = \lambda_{\min}(\rho I-A),

$f(x^k)$ decreases to a limit $\alpha'\ge\alpha$ , $\|x^{k+1}-x^k\|\to0$ , and every limit point of $\{x^k\}$ is a Kuhn–Tucker point.

Milestones

In attack order:

(P_k): $\partial g^*(y)$ is the set of minimizers of $g-\langle\cdot,y\rangle$ .
(18)–(19): the DCA step solves (18), is the projection (19), and is a simplified-DCA step for (16).
Proposition A.1(i): $(g-h)(x^{k+1})\le(h^*-g^*)(y^k)-\frac{\rho_2}2\|dx^k\|^2\le(g-h)(x^k)-\frac{\rho_1+\rho_2}2\|dx^k\|^2$ , with its equality case.
Corollary A.2(i)–(ii): the same descent with the dual decrements.
Lemma A.3: continuity of $h$ along sequences with bounded subgradients.
Theorem 3.7(i): the descent with maxima of the four weighted steps.
Theorem 3.7(iii): if $\alpha$ is finite, both objectives converge to a common $\beta\ge\alpha$ , the steps vanish, and the Fenchel–Young gaps tend to $0$ .
Theorem 3.7(iv): with bounded iterates, limit points pair into $(x^*,y^*)\in[\partial g^*(y^*)\cap\partial h^*(y^*)]\times[\partial g(x^*)\cap\partial h(x^*)]$ .
The moduli of (16): $\rho(g)=\rho$ and $\rho(h)=\rho-\lambda_n$ .
(20): $\partial h(x^*)\subset\partial g(x^*)$ if and only if $x^*$ is a Kuhn–Tucker point.

The equality case of Theorem 3.7(i) and Theorem 3.7(ii) (the dual descent) are included as further draft theorems.

Significance

Theorem 4.1 makes the DCA a descent method for $(Q_1)$ whose accumulation points satisfy the first-order conditions. That property is what the paper's restarting procedure (§4.2) builds on: at a Kuhn–Tucker point that is not a global minimizer, one eigenvector computation gives a better feasible point, and the DCA reaches a global solution after at most $2m+2$ restarts, where $m$ is the number of distinct negative eigenvalues of $A$ . Theorem 3.7 is the convergence theorem for the DCA on general d.c. programs, used throughout the later DCA literature.

All of these results are proved in the paper. As far as is known, none of them has a machine-checked proof: Mathlib has convex conjugates and subdifferentials in pieces, but no d.c. duality, no DCA and no analysis of the trust-region subproblem. This mission produces that development: a statement layer for d.c. programming with extended values, a formal proof of the general descent and criticality theorem, and its specialization to $(Q_1)$ .

Difficulty

Every inequality of the analysis runs over $\mathbb R\cup\{+\infty\}$ . The differences $g-h$ and $h^*-g^*$ need the paper's convention $+\infty-(+\infty)=+\infty$ , and the finiteness of each value along a run has to be derived from the standing inclusions. The hard step is the passage to limit points in Theorem 3.7(iv). Lower semicontinuity of $g$ , $g^*$ and $h$ gives only one-sided inequalities at a limit point, while the theorem asserts equalities of values there. The naive argument of passing to the limit in $y^k\in\partial h(x^k)$ and $x^{k+1}\in\partial g^*(y^k)$ fails, because $x^{k+1}$ and $x^k$ need not have the same limit unless $\rho_1+\rho_2>0$ . For $(Q_1)$ the specialization needs boundedness of the iterates, finiteness of $\alpha$ , and the subdifferential of the indicator of the ball. None of these is a hypothesis of Theorem 4.1, so each must be proved from the run.

Formalization scope

$X$ is EuclideanSpace ℝ (Fin n). Extended values are EReal, and every d.c. difference is dcSub g h x, equal to $+\infty$ when $g(x)=+\infty$ and to $g(x)-h(x)$ otherwise, because Mathlib's EReal has $\top-\top=\bot$ . $\Gamma_0$ , subgradients and the conjugate come from published platform definitions (IsProperClosedConvex, IsSubgradient, CondatPD.FinDim.conj). The modulus of strong convexity is never defined as a supremum. Statements take explicit constants $\rho_1,\rho_2,\rho_1^*,\rho_2^*$ with " $f_i$ is $\rho_i$ -convex", which is what the paper requires of them. Hypotheses " $\rho(g)+\rho(h)>0$ " become $\rho_1+\rho_2>0$ . The matrix $A$ is a self-adjoint continuous linear operator. " $\lambda$ is the smallest eigenvalue of $\rho I-A$ " is encoded by quantifying over every $\mu$ with $\rho I-A\succeq\mu I$ . A limit point is a cluster point. The objective $f$ is $+\infty$ outside $E$ , and the run starts from an arbitrary $x^0$ , as on the page. Statements involving $dy^{k-1}$ are for $k\ge1$ . Theorem 3.7(iv)'s limit identity is stated along a common convergent subsequence, as in its proof.

Ruled out: the goal is stated for the explicit iteration a)/b) of p. 490, never for an abstract run assumed to satisfy Theorem 3.7, and it carries no hypothesis of boundedness, finite $\alpha$ or convergence. A run predicate that no sequence satisfies would make everything vacuous, and a sanity file exhibits concrete runs.

Needed infrastructure: the Fenchel–Young equality and its equality case for $\Gamma_0$ functions in EReal, subdifferential inversion $y\in\partial g(x)\iff x\in\partial g^*(y)$ , strong-convexity versions of the subgradient inequality, and the normal cone of a Euclidean ball. These are reusable well beyond this mission and are welcome as separate contributions.

Selected references

Pham Dinh Tao and Le Thi Hoai An, A d.c. optimization algorithm for solving the trust-region subproblem, SIAM J. Optim. 8(2):476–505, 1998. https://doi.org/10.1137/S1052623494274313
J. F. Toland, Duality in nonconvex optimization, J. Math. Anal. Appl. 66(2):399–415, 1978. https://doi.org/10.1016/0022-247X(78)90243-3
J. J. Moré and D. C. Sorensen, Computing a trust region step, SIAM J. Sci. Statist. Comput. 4(3):553–572, 1983. https://doi.org/10.1137/0904038
A. R. Conn, N. I. M. Gould and Ph. L. Toint, Trust-Region Methods, SIAM, 2000. https://doi.org/10.1137/1.9780898719857
R. T. Rockafellar, Convex Analysis, Princeton University Press, 1970. https://doi.org/10.1515/9781400873173

15 thms2 active usersReviewed

Optimization·Captain: mikedeng1

Convergence Analysis of Some Algorithms for Solving Nonsmooth Equations I: At a Semismooth, Strongly BD-Regular Zero the ∂_B Newton Method Converges Superlinearly and ‖F‖ Decreases SuperlinearlyResearch Paper

Motivation

Many problems of operations research reduce to a system of nonsmooth equations $F(x)=0$ with $F:\mathbb R^n\to\mathbb R^n$ locally Lipschitz but not differentiable: the nonlinear complementarity problem $h(x)\ge 0$ , $f(x)\ge0$ , $h(x)^Tf(x)=0$ becomes $F(x)=\min\{h(x),f(x)\}=0$ componentwise, and variational inequalities, KKT systems and max-type equations admit similar reformulations. Newton's method is the natural local solver, but it needs a substitute for the Jacobian at points where $F$ is not differentiable.

Timeline:

1977. Mifflin introduces semismooth functionals.
1987–1990. Robinson proposes B-derivatives; Pang (1990) and Harker–Xiao (1990) study Newton methods based on directional derivatives (Math. Oper. Res. 15; Math. Programming 48).
1993. Qi and Sun (Math. Programming 58) prove local superlinear convergence of $x^{k+1}=x^k-V_k^{-1}F(x^k)$ , $V_k\in\partial F(x^k)$ , assuming every element of Clarke's generalized Jacobian $\partial F(x^*)$ is nonsingular.
1993. Qi (Math. Oper. Res. 18), the source of this mission, observes that this assumption already fails for $f(x)=|x|$ , where $\partial f(0)=[-1,1]$ contains $0$ , and restricts $V_k$ to the smaller B-subdifferential. Nonsingularity is then needed only on $\partial_B F(x^*)$ .

Setting

A map $F:\mathbb R^n\to\mathbb R^m$ is locally Lipschitz throughout; the norm is the Euclidean 2-norm. Let $D_F$ be the set of points where $F$ is differentiable (almost every point, by Rademacher's theorem). The B-subdifferential of $F$ at $x$ is

\partial_B F(x)=\Big\{\lim_{i\to\infty}\nabla F(x_i)\;:\;x_i\to x,\ x_i\in D_F\Big\},

and Clarke's generalized Jacobian is its convex hull, $\partial F(x)=\operatorname{conv}\partial_B F(x)$ . The directional derivative is $F'(x;h)=\lim_{t\downarrow0}(F(x+th)-F(x))/t$ . $F$ is semismooth at $x$ if $\lim Vh'$ over $V\in\partial F(x+th')$ , $h'\to h$ , $t\downarrow 0$ , exists for every $h$ . A map $F:\mathbb R^n\to\mathbb R^n$ is strongly BD-regular at $x$ if every $V\in\partial_B F(x)$ is nonsingular.

The $\partial_B$ Newton method is

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

with any choice of $V_k$ , and the directional-derivative Newton method is $x^{k+1}=x^k+d^k$ with $F(x^k)+F'(x^k;d^k)=0$ (3.11)–(3.12). In Lean, $\mathbb R^n$ is EuclideanSpace ℝ (Fin n), $\partial_B F$ is NonsmoothNewton.Shared.bJac F, $\partial F$ is clarkeJac F, semismoothness is NonsmoothNewton.Local.SemismoothAt, and a run of (3.2) is IsBDNewtonRun F x V.

Formalization targets

Goal: Theorem 3.1 (p. 234)

If $F(x^*)=0$ and $F$ is semismooth and strongly BD-regular at $x^*$ , there is $\delta>0$ such that every $V\in\partial_B F(y)$ , $\|y-x^*\|<\delta$ , is nonsingular with $\|y-V^{-1}F(y)-x^*\|<\delta$ , and every run of (3.2) from $\|x^0-x^*\|<\delta$ satisfies

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

and, when $F(x^k)\ne0$ for all $k$ ,

\lim_{k\to\infty}\frac{\|F(x^{k+1})\|}{\|F(x^k)\|}=0. \tag{1.9}

Milestones

Lemma 2.6 (p. 233): near a strongly BD-regular point, every $V\in\partial_B F(y)$ is nonsingular with $\|V^{-1}\|\le c$ , and $\|h\|\le c\|F'(y;h)\|$ where $F$ is semismooth.
Lemma 2.2 (1) ⇔ (3) (p. 230): semismoothness at $x$ is equivalent to $Vh-F'(x;h)=o(\|h\|)$ for $V\in\partial F(x+h)$ . This is the published item NonsmoothNewton.Local.theorem_2_3.
(2.2) (p. 229) at a semismooth point: $F(x+h)=F(x)+F'(x;h)+o(\|h\|)$ . This is the published item NonsmoothNewton.Local.remark_2_17.
Corollary 3.2 (p. 235): for every $\varepsilon>0$ , one $\partial_B$ Newton step near $x^*$ contracts $\|x-x^*\|$ and $\|F(x)\|$ by the factor $\varepsilon$ .

Further statements of the paper in this mission

Lemma 2.3 (degree-2 semicontinuity, p. 230); the quadratic second paragraph of Theorem 3.1; Corollary 3.3 (the directional-derivative method, p. 235); Lemma 2.2 (1) ⇔ (2); and the counterexample of Final Remark (6) (pp. 243–244), which shows that (3.12) need not be solvable near a semismooth, strongly BD-regular zero.

Significance

Theorem 3.1 is the local convergence theorem for the semismooth Newton method in the form most used later: the B-subdifferential version under strong BD-regularity is what generalized Newton methods for complementarity problems (through the min and Fischer–Burmeister reformulations) rely on, because $\partial_B F$ is finite in typical piecewise-smooth examples while $\partial F$ is a polytope that may contain singular matrices. The statement (1.9) makes the method a descent method for $\|F\|$ near the solution, which is how the paper later globalizes it (damped Newton, §4; hybrid method, §6). Corollary 3.2 is reused verbatim in the analysis of the damped method.

The results are proved in the paper; none of them is formalized. The mission produces a machine-checked version of the local theory, on top of the published B-subdifferential and semismoothness definitions of the Qi–Sun series, and makes the Clarke-to-B-subdifferential refinement available to later missions.

Difficulty

The Clarke-Jacobian argument cannot be copied. There, nonsingularity of the compact convex set $\partial F(x^*)$ propagates to a neighbourhood by upper semicontinuity of $\partial F$ . Here only the extreme set $\partial_B F(x^*)$ is assumed nonsingular, and its convex hull may contain singular matrices. Lemma 2.6 therefore has to work with limits of genuine Jacobians, and the error estimate (3.3) has to combine an estimate stated for $\partial F(x^k)$ , namely Lemma 2.2, with the bound on $V_k^{-1}$ , which exists only for $V_k\in\partial_B F(x^k)$ . A second subtlety is (1.9), which needs a lower bound $\|x^k-x^*\|\le\frac{c}{1-\varepsilon}\|F(x^k)\|$ that does not follow from superlinear convergence alone.

Formalization scope

$\mathbb R^n$ is EuclideanSpace ℝ (Fin n), with the 2-norm and the spectral operator norm, as the paper fixes ("we use 2-norms"). The dimension $n=0$ is allowed.
Standing assumption. Every statement carries LocallyLipschitz F (§1, p. 227).
Nonsingularity is the existence of a two-sided inverse $W$ (IsInverse V W), never a junk inverse. A run of (3.2) is written $V_k(x^{k+1}-x^k)=-F(x^k)$ without an inverse, so nonsingularity of $V_k$ is a conclusion.
Quantifier order. $\delta$ is chosen before the run and serves every starting point within $\delta$ and every choice $V_k\in\partial_B F(x^k)$ . Writing "for every run there is a $\delta$ " would trivialize the theorem and is ruled out.
(1.9) is stated under the paper's hypothesis $F(x^k)\ne0$ for all $k$ , since otherwise the quotient is not defined.
Disclosed deviations.
- Corollary 3.2 is stated for $\varepsilon>0$ (the page prints $\varepsilon\ge0$ , false at $0$ ).
- Corollary 3.3 adds semismoothness of $F$ on a neighbourhood of $x^*$ , which its printed proof uses through Lemma 2.6's (2.16).
- $O(\cdot)$ and "quadratic" are pinned to explicit constants.
Infrastructure needed.
- Rademacher-type facts for $\partial_B F$ : nonemptiness, compactness and closed graph. The platform already has WillR's proved lemmas NonsmoothNewton.Shared.exists_local_bJac_control, bJac_mem_of_tendsto, isCompact_bJac.
- Perturbation bounds for inverses.
- The relation between directional derivatives and $\partial_B F$ at semismooth points.
These are reusable for the other two missions of this series and for later semismooth Newton developments. Proofs of the milestones, and alternative arguments, are welcome.

Selected references

L. Qi, Convergence analysis of some algorithms for solving nonsmooth equations, Math. Oper. Res. 18(1):227–244, 1993. https://doi.org/10.1287/moor.18.1.227
L. Qi and J. Sun, A nonsmooth version of Newton's method, Math. Programming 58:353–367, 1993. https://doi.org/10.1007/BF01581275
J.-S. Pang, Newton's method for B-differentiable equations, Math. Oper. Res. 15(2):311–341, 1990. https://doi.org/10.1287/moor.15.2.311
R. Mifflin, Semismooth and semiconvex functions in constrained optimization, SIAM J. Control Optim. 15(6):959–972, 1977. https://doi.org/10.1137/0315061
A. Shapiro, On concepts of directional differentiability, J. Optim. Theory Appl. 66:477–487, 1990. https://doi.org/10.1007/BF00940933

9 thms2 active usersReviewed

Completed

Convex Optimization Optimization·Captain: mikedeng1

Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers III: Parallel Projections — the Dual Average Vanishes and ADMM Reduces to Project-and-AverageTextbook

Motivation

Many feasibility problems ask for a point satisfying several constraints at once. When each constraint set admits a manageable Euclidean projection, a parallel method can project onto all sets independently and combine the results. Chapter 5 of Boyd, Parikh, Chu, Peleato and Eckstein (2011) places this construction inside the alternating direction method of multipliers (ADMM). That connection matters when a single projection onto the intersection is difficult but projections onto its individual factors are available.

The chapter first rewrites constrained convex optimization as a split problem. It then treats feasibility for two sets and finally an arbitrary finite collection of sets. For the latter, the product of the sets captures the separate constraints, while a consensus set forces all copies of the variable to agree. These are ordinary finite dimensional convex sets, so the account applies to geometric feasibility and to constrained optimization subproblems in which projection is the computational primitive. The mission isolates the identities that let the product-space algorithm run through parallel individual projections and one average.

Setting

A Euclidean projection of a vector $v$ onto a set $C$ is a point $p\in C$ closest to $v$ in Euclidean distance. Each set $\mathcal A_i\subseteq\mathbb R^n$ , for $i=1,\ldots,N$ , is assumed nonempty, closed and convex, and $N\ge1$ . These properties give a unique projection for every input. They do not require the common intersection $\bigcap_i\mathcal A_i$ to be nonempty: the iteration and the algebraic identities still make sense when the feasibility problem has no solution.

A block vector $x=(x_1,\ldots,x_N)$ belongs to the product set $\mathcal C=\mathcal A_1\times\cdots\times\mathcal A_N$ exactly when $x_i\in\mathcal A_i$ for all $i$ . The consensus set $\mathcal D$ consists of constant block vectors $(z,\ldots,z)$ . The block average is $\bar x=N^{-1}\sum_i x_i$ . Distances in the stacked space use the ordinary Euclidean metric, whose square is $\sum_i\|x_i-y_i\|_2^2$ . A block vector need not carry a function-space norm in the formalization for this metric to be stated precisely.

Starting with arbitrary $z^0\in\mathbb R^n$ and $u^0=(u_i^0)$ , the parallel projection ADMM run has three updates: each $x_i^{k+1}$ is the projection of $z^k-u_i^k$ onto $\mathcal A_i$ ; $z^{k+1}$ is the average of $x_i^{k+1}+u_i^k$ ; and $u_i^{k+1}=u_i^k+x_i^{k+1}-z^{k+1}$ . The $u_i^k$ are scaled dual variables, formed using a penalty parameter $\rho>0$ . The source presents these equations as the specialization of scaled ADMM for the indicator functions of $\mathcal C$ and $\mathcal D$ §5.1.2, pp. 35–36.

Formalization targets

The first targets identify the projections that make the iteration parallel:

\Pi_{\mathcal C}(x)=(\Pi_{\mathcal A_1}(x_1),\ldots,\Pi_{\mathcal A_N}(x_N)),\qquad \Pi_{\mathcal D}(x)=(\bar x,\ldots,\bar x).

The goal is the full reduction of a run. Its dual average is zero after every update, and its common vector equals the primal average from the following update onward:

\bar u^{k+1}=0\quad(k\ge0),\qquad z^{k+1}=\bar x^{k+1}\quad(k\ge1).

Consequently, for every $k\ge2$ and every $i$ , the same iterates obey the simpler two-line algorithm

x_i^{k+1}=\Pi_{\mathcal A_i}(\bar x^k-u_i^k),\qquad u_i^{k+1}=u_i^k+x_i^{k+1}-\bar x^{k+1}.

If $\bar u^0=0$ , those last identities hold already for $k\ge1$ . The index qualification is necessary because $z^0$ is arbitrary, and $z^1$ still contains the initial dual average unless that average was zero. This is the precise reading of the source's “after the first step” statement §5.1.2, p. 36. The milestone list also records the indicator-function projection update on p. 33 and the averaged identities on p. 36.

Significance

The result identifies what information must pass between the parallel projection steps: each local step needs the shared average and its own dual displacement. It also fixes the relationship between the displayed three-variable ADMM iteration and the shorter project-and-average rule. Without the dual-average identity, replacing $z^k$ by $\bar x^k$ is not justified at the initial indices. The product and consensus projection identities explain why the update can be decomposed into independent local projections followed by a mean §5.1.2, p. 35.

The chapter's identities are known mathematical results; this mission asks for machine-checked proofs of their stated finite dimensional version. A completed development supplies reusable statements about nearest-point projection onto products and diagonal sets, and a precise run predicate for later analysis of parallel ADMM. The goal is a statement about every valid run rather than a claim of convergence or existence of a feasible common point. Neither convergence nor feasibility of the intersection is concluded on these pages. Those distinctions keep the mission aligned with the actual scope of Chapter 5.

Difficulty

The product projection uses a sum of squared component distances. Treating a block vector as an arbitrary function and applying its default norm would change the metric, so a projection theorem stated that way would describe a different algorithm. The consensus projection likewise depends on the Euclidean sum and the finite average, including the condition $N\ge1$ . The run reduction has a separate indexing difficulty: the vanishing average follows after one dual update, but the first simplified projection step can occur later because it uses the previous common vector. An argument that silently assumes $z^0=\bar x^0$ would erase the arbitrary initialization on the source page.

Formalization scope

Lean represents each $\mathbb R^n$ as EuclideanSpace ℝ (Fin n) and the $N$ blocks as a function on Fin N; the latter is the source's $1$ -based indexing shifted to $0$ -based indices. Squared stacked distance is defined as the sum of squared component norms. The average uses the real factor $1/N$ only under $0<N$ . Individual projections reuse the published nearest-point relation RandomGradFree.Nonsmooth.IsMetricProjection. A run records the three updates as conditions on sequences, with arbitrary $z^0,u^0$ and unused $x^0$ . The source's scaled dual variable $u$ appears directly; $\rho$ cancels from these projection updates once it is positive.

The source assumes closed convex sets; nonemptiness is stated explicitly here because a projection onto an empty set does not exist. The book prints $-z^k$ in its averaged dual update on p. 36, while its defining update on p. 35 and the next line on p. 36 require $-z^{k+1}$ . The formal statement uses the corrected index. The book also calls the individual sets $\mathcal C_i$ in its closing sentence; the sets defined for this subsection are $\mathcal A_i$ . The milestone quote preserves the printed intermediate index. The goal includes the simplified primal and dual updates, since asserting only that the dual average vanishes would omit the chapter's advertised reduction.

The development needs finite sums in Euclidean spaces, convex-set projection, and elementary properties of product and diagonal sets. Product and consensus projection results are useful beyond ADMM. Contributions proving the four milestones, the full run theorem, or equivalent projection lemmas at the same generality are within scope.

Selected references

S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein, Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers, Foundations and Trends in Machine Learning 3(1), 2011, pp. 1–122. DOI: 10.1561/2200000016.

8 thms2 active usersReviewed

Completed

Convex Optimization Optimization·Captain: mikedeng1

Convex Optimization: Algorithms and Complexity XIII: Newton's Method Converges Quadratically, ‖x_{k+1} − x‖ ≤ (M/μ)‖x_k − x‖², from ‖x₀ − x*‖ ≤ μ/(2M)Textbook

Motivation

Newton's method is the basic second-order method of continuous optimization: at the current point it replaces the objective by its second-order Taylor model and jumps to the stationary point of that model. Its defining property is speed near a nondegenerate minimum, where the error is squared at every step, so that the number of correct digits roughly doubles per iteration. This local behaviour is what makes Newton's method the inner engine of interior point methods, the polynomial-time algorithms for linear, conic and general convex programming (Nesterov and Nemirovski, 1994). In S. Bubeck's monograph Convex Optimization: Algorithms and Complexity (Foundations and Trends in Machine Learning, 2015; arXiv:1405.4980v2), §5.3.2 recalls the traditional local analysis of Newton's method, Theorem 5.3, before turning to the affine-invariant self-concordance analysis used for interior point methods. This mission formalizes that theorem and the four steps of its proof.

Setting

Let $\mathbb R^n$ carry the Euclidean norm $\|\cdot\|$ , and write $\|A\|$ for the operator norm of a linear map $A:\mathbb R^n\to\mathbb R^n$ , so that $\|Ax\|\le\|A\|\,\|x\|$ . Let $f:\mathbb R^n\to\mathbb R$ be a $C^2$ function, with gradient $\nabla f(x)\in\mathbb R^n$ and Hessian $\nabla^2 f(x)$ , a linear map $\mathbb R^n\to\mathbb R^n$ (the derivative of the gradient map). For a real number $c$ , $A\succeq cI_n$ means $\langle Av,v\rangle\ge c\|v\|^2$ for all $v\in\mathbb R^n$ .

The Hessian is $M$ -Lipschitz if $\|\nabla^2 f(x)-\nabla^2 f(y)\|\le M\|x-y\|$ for all $x,y\in\mathbb R^n$ .

Newton's method starts at $x_0\in\mathbb R^n$ and iterates, for $k\ge0$ ,

x_{k+1}=x_k-[\nabla^2 f(x_k)]^{-1}\nabla f(x_k).

A point $x^*$ is a local minimum of $f$ if $f(x^*)\le f(x)$ for all $x$ in a neighbourhood of $x^*$ ; it has strictly positive Hessian if $\nabla^2 f(x^*)\succeq\mu I_n$ for some $\mu>0$ .

Formalization targets

Goal: Theorem 5.3 (p. 320)

Assume the Hessian of $f$ is $M$ -Lipschitz, $M>0$ , and $x^*$ is a local minimum with $\nabla^2 f(x^*)\succeq\mu I_n$ , $\mu>0$ . If $\|x_0-x^*\|\le\mu/(2M)$ , then Newton's method from $x_0$ is well defined (every Hessian along the iterates is invertible, so the sequence exists and is unique) and

\|x_{k+1}-x^*\|\le\frac M\mu\,\|x_k-x^*\|^2\quad(k\ge0),\qquad x_k\to x^*.

Milestones (p. 321, the steps of the proof)

The integral formula $\int_0^1\nabla^2 f(x+sh)\,h\,ds=\nabla f(x+h)-\nabla f(x)$ .
The error representation of one Newton step, $x_{k+1}-x^*=[\nabla^2 f(x_k)]^{-1}\int_0^1[\nabla^2 f(x_k)-\nabla^2 f(x^*+s(x_k-x^*))](x_k-x^*)\,ds$ .
The Lipschitz bound $\int_0^1\|\nabla^2 f(x_k)-\nabla^2 f(x^*+s(x_k-x^*))\|\,ds\le\frac M2\|x_k-x^*\|$ .
The Hessian lower bound $\nabla^2 f(x_k)\succeq(\mu-M\|x_k-x^*\|)I_n\succeq\frac\mu2I_n$ when $\|x_k-x^*\|\le\mu/(2M)$ .

Significance

The theorem gives a quantitative basin of quadratic convergence: an explicit radius $\mu/(2M)$ , depending only on the curvature at the minimum and the Lipschitz constant of the Hessian, inside which Newton's method needs only $O(\log\log(1/\varepsilon))$ iterations to reach accuracy $\varepsilon$ . It is the classical statement whose shortcomings (dependence on a choice of norm, constants that change under linear changes of variables) motivate the self-concordance theory of the following subsections, and it is the local convergence result invoked whenever a damped or globalized Newton scheme is shown to enter its quadratic phase.

On the formal side, Mathlib has the calculus this needs (Fréchet derivatives, interval integrals of vector-valued maps, operator norms) but no convergence theorem for multivariate Newton's method for minimization. A formal proof produces reusable pieces: the integral form of the mean value theorem for gradients, the stability of a positive-definite lower bound under Lipschitz perturbations, and an inverse-operator norm bound from a quadratic-form lower bound. The result itself is classical and fully proved in the literature; what is open here is its machine-checked proof in this form.

Difficulty

The individual inequalities are short, but the argument is an induction in which well-definedness and the rate are proved together: the Hessian at $x_k$ is invertible only because $x_k$ is still in the ball of radius $\mu/(2M)$ , and $x_{k+1}$ stays in that ball only because of the rate. A proof that first assumes the sequence exists and then bounds it is circular. The proof also passes between two kinds of control on the Hessian, a lower bound on its quadratic form and an operator-norm bound on its inverse, and the second is only meaningful once invertibility is established. Finally, the integral manipulations need integrability of the maps $s\mapsto\nabla^2 f(x^*+s(x_k-x^*))(x_k-x^*)$ , which comes from the continuity of the Hessian.

Formalization scope

$\mathbb R^n$ is EuclideanSpace ℝ (Fin n). The gradient and Hessian are explicit maps $g:\mathbb R^n\to\mathbb R^n$ and $H:\mathbb R^n\to(\mathbb R^n\to_L\mathbb R^n)$ with ContDiff ℝ 2 f, HasGradientAt f (g x) x and HasFDerivAt g (H x) x at every point; the norm on $H(x)$ is Mathlib's operator norm, as on the page. $A\succeq cI_n$ is the quadratic-form inequality. A Newton run is a sequence $x:\mathbb N\to\mathbb R^n$ indexed from $0$ satisfying the linear system $\nabla^2 f(x_k)(x_k-x_{k+1})=\nabla f(x_k)$ ; no inverse of a possibly singular operator appears in any hypothesis, and "well defined" is a conclusion: a unique run exists from $x_0$ and every Hessian along it is bijective. The rate and $x_k\to x^*$ are asserted for every run. The error representation is stated with both sides multiplied by $\nabla^2 f(x_k)$ , which is equivalent to the printed form once the Hessian is invertible. Milestones 3 and 4 use only the Lipschitz property and are stated for any Lipschitz map $H$ .

Added hypothesis: $M>0$ (the radius $\mu/(2M)$ divides by $M$ ; with $M=0$ , Lean's convention $\mu/0=0$ would collapse the hypothesis to $x_0=x^*$ ). Convexity of $f$ is not assumed, as on the page; $x^*$ is a local minimum and $\nabla f(x^*)=0$ is derived, not assumed. Encoding the Newton step with Lean's inverse (which returns $0$ on singular maps), or replacing $\nabla^2 f(x^*)\succeq\mu I_n$ by mere invertibility, would change the theorem and is ruled out.

A complete development needs the fundamental theorem of calculus for $C^1$ vector-valued maps along segments, Hessian-based quadratic-form estimates, and operator-norm bounds for inverses; all are reusable for the analysis of damped Newton, cubic regularization and interior point methods. Proofs of the milestones independently of the goal are welcome.

Selected references

S. Bubeck, Convex Optimization: Algorithms and Complexity, Foundations and Trends in Machine Learning 8(3–4):231–358, 2015. arXiv:1405.4980v2, §5.3.2, Theorem 5.3, pp. 320–321.
Yu. Nesterov and A. Nemirovski, Interior-Point Polynomial Algorithms in Convex Programming, SIAM Studies in Applied Mathematics 13, 1994. doi:10.1137/1.9781611970791
Yu. Nesterov, Introductory Lectures on Convex Optimization: A Basic Course, Kluwer, 2004, Theorem 1.2.5. doi:10.1007/978-1-4419-8853-9

6 thms2 active usersReviewed

Completed

Convex Optimization Partial Differential Equations·Captain: mikedeng1

Mean Field Games: Numerical Methods for the Planning Problem II: Solutions of the Penalized Scheme Converge to a Solution of the Discrete Planning Scheme as ε → 0Research Paper

Motivation

Mean field games (MFG), introduced by Lasry and Lions and by Huang, Caines and Malhamé, describe Nash equilibria of very large populations of identical rational agents. The equilibrium is a coupled system: a backward Hamilton–Jacobi–Bellman equation for the value function $u$ of a representative agent and a forward Fokker–Planck equation for the density $m$ of the population. In the planning problem, proposed by P.-L. Lions, both the initial density $m_0$ and the final density $m_T$ are prescribed, and one asks for a cost structure under which the population moves from one to the other. This is a mean field analogue of optimal transport.

Achdou, Camilli and Capuzzo-Dolcetta (hal-00465404, SIAM J. Control Optim. 2012) propose finite-difference schemes for the planning problem. Mission I of this series treats the existence of a solution of the discrete planning scheme. Since the two boundary conditions on $m$ make the discrete system hard to solve directly, the paper also introduces a penalized scheme, in which the initial condition $M^0 = m_0$ is replaced by a penalty $U^0 = (M^0 - m_0)/\varepsilon$ ; for each $\varepsilon > 0$ this is a standard discrete MFG system with a unique solution. This mission formalizes the paper's §3.2: the penalized solutions converge, as $\varepsilon \to 0$ , to a solution of the planning scheme.

Setting

Fix integers $N_h \ge 1$ , $N_T \ge 1$ , a horizon $T > 0$ and a viscosity $\nu \ge 0$ ; put $h = 1/N_h$ and $\Delta t = T/N_T$ . The grid $\mathbb T^2_h$ consists of the points $x_{i,j}$ , $(i,j) \in (\mathbb Z/N_h\mathbb Z)^2$ (periodic indices). A grid function is a real function on $\mathbb T^2_h$ ; time levels are $n = 0, \dots, N_T$ .

$(D_1^+U)_{i,j} = (U_{i+1,j} - U_{i,j})/h$ , $(D_2^+U)_{i,j} = (U_{i,j+1} - U_{i,j})/h$ , and the discrete gradient is $[D_hU]_{i,j} = ((D_1^+U)_{i,j}, (D_1^+U)_{i-1,j}, (D_2^+U)_{i,j}, (D_2^+U)_{i,j-1}) \in \mathbb R^4$ .
$\Delta_h$ is the five-point Laplacian.
A numerical Hamiltonian $g(x_{i,j}, q)$ , $q \in \mathbb R^4$ , is monotone (G1), $C^1$ (G3), convex (G4) and coercive (G5) in $q$ .
$\mathcal B_{i,j}(U, M)$ is the discrete transport term $\mathrm{div}_h\big(M\,\nabla_q g(\cdot, [D_hU])\big)$ .
$W : \mathbb R \to \mathbb R$ is strictly convex, superlinear and $C^2$ , with $V = W'$ (hypothesis (24)).
$\mathcal K = \{M \ge 0 : h^2\sum_{i,j} M_{i,j} = 1\}$ is the set of discrete probability densities, and $m_0, m_T \in \mathcal K$ with $m_0 > 0$ .

The planning scheme (18) asks for families $(U^n, M^n)$ with, for $n < N_T$ ,

\frac{U^{n+1} - U^n}{\Delta t} - \nu\Delta_hU^{n+1} + g(x, [D_hU^{n+1}]) = V(M^n), \qquad \frac{M^{n+1} - M^n}{\Delta t} + \nu\Delta_hM^n + \mathcal B(U^{n+1}, M^n) = 0,

$M^n \in \mathcal K$ , $M^0 = m_0$ and $M^{N_T} = m_T$ . The penalized scheme (20)–(23) has the same two equations, $M^n \in \mathcal K$ for $n < N_T$ , $M^{N_T} = m_T$ , and $U^0 = (M^0 - m_0)/\varepsilon$ in place of $M^0 = m_0$ .

Formalization targets

Goal: Proposition 4

Under the hypotheses of Theorem 1, for every sequence $\varepsilon_k \downarrow 0$ , every choice of penalized solutions $(U^{\varepsilon_k}, M^{\varepsilon_k})$ , and every limit $M^{\varepsilon_k} \to M$ in $\mathcal K^{N_T+1}$ ,

\exists\, U,\ \text{a subsequence with } U^{\varepsilon_k} \to U, \quad (U, M) \text{ solves (43)–(46)}.

If $g$ is strictly convex, the whole sequence $(U^{\varepsilon_k}, M^{\varepsilon_k})$ converges to the unique solution of (43)–(46) with $\sum_{i,j} U^0_{i,j} = 0$ .

Milestones

(13): (G5) implies $\max_{i,j} g(x_{i,j}, [D_hU]_{i,j}) / \|[D_hU]\|_\infty \to +\infty$ .
Theorem 2: the penalized control problem (50), $\min \Theta^*(M,Z) + \frac{1}{2\varepsilon\Delta t}\sum(M^0 - m_0)^2$ under the discrete Fokker–Planck constraint, has a minimizer whose optimality system is (20)–(23).
Proposition 2: $\max_{i,j}|M^{\varepsilon,0}_{i,j} - (m_0)_{i,j}| \le C\varepsilon^{1/2}$ .
Proposition 3: $\max_{n,i,j}|U^{\varepsilon,n}_{i,j}| \le C$ .
Corollary 1: $\max_{i,j}|M^{\varepsilon,0}_{i,j} - (m_0)_{i,j}| \le C\varepsilon$ .
Proposition 1: uniqueness of $M$ for (18), and of $U$ normalized by $\sum U^0 = 0$ when $g$ is strictly convex.

In 3–5 the constant $C$ depends on $h$ , $\Delta t$ and the data, never on $\varepsilon$ .

Significance

Proposition 4 justifies the penalized scheme as a way to compute solutions of the discrete planning problem. The penalized system has a standard initial–terminal structure that Newton's method (the paper's §4) handles well. The planning system does not, because it prescribes two conditions on $M$ and none on $U$ . The estimates of Propositions 2–3 and Corollary 1 also give a rate: the initial mismatch of the penalized density is $O(\varepsilon)$ .

The results are proved in the paper. None of them has been machine-checked: the platform holds no formalization of mean field games, of their finite-difference schemes, or of the Lasry–Lions monotonicity argument. A complete development would yield a reusable library of the monotone finite-difference operators on the discrete torus, a formal Lasry–Lions uniqueness argument for discrete MFG systems, and a template for compactness-plus-uniqueness convergence proofs in finite dimension.

Difficulty

Compactness of the densities is free, since $\mathcal K^{N_T+1}$ is compact. Passing to the limit in the scheme needs only continuity of $g$ , $\nabla_q g$ and $V$ . The substance lies in the two uniform estimates.

Bounding $U^\varepsilon$ uniformly in $\varepsilon$ (Proposition 3) is the central difficulty. The initial condition $U^{\varepsilon,0} = (M^{\varepsilon,0} - m_0)/\varepsilon$ involves division by $\varepsilon$ , so the obvious bound $|U^{\varepsilon,0}| \le 2/(h^2\varepsilon)$ blows up. The comparison principle for the discrete HJB equation, which is the natural first idea, cannot repair this: it propagates whatever bound $U^0$ has. The HJB equation alone does not bound $U^\varepsilon$ ; the coupling with the Fokker–Planck equation, the coercivity (13), and a lower bound on $M^{\varepsilon,0}$ that is uniform in $\varepsilon$ (which is where Proposition 2 enters) all have to be used.
Proposition 2 compares the penalized problem with the planning problem through the convex duality of Theorem 2. It therefore needs a solution of the planning scheme, which is mission I's goal.
Proposition 1 needs the discrete Lasry–Lions identity, a summation by parts across the coupled system.

Formalization scope

All objects live in the namespace MFGPlanning.Penalized.

Grid points are ZMod Nh × ZMod Nh (periodicity is built in), time levels are Fin (NT + 1), and the four momentum components $q_1, \dots, q_4$ are Fin 4 indices 0, …, 3.
The data structure carries $N_h, N_T \ge 1$ , $T > 0$ and $\nu \ge 0$ .
$g$ is given only at the grid points. (G2), which only defines the continuous Hamiltonian, is not encoded.
"Coercive" in (24) is read as superlinear.
$V_h[M] = V(M_{i,j})$ , as (24) prescribes.
$\Theta^*$ and $(W+\chi)^*$ take values in EReal, so unbounded suprema are $+\infty$ .
Convergence is in the product topology, which on these finite-dimensional spaces is $\max_n \|\cdot\|_\infty$ convergence.
"The unique solution" of (43)–(46) means unique under the normalization $\sum_{i,j} U^0_{i,j} = 0$ , because $U$ is otherwise determined only up to a constant.

Three trivializing formalizations are ruled out. The constants in Propositions 2–3 and Corollary 1 are quantified before $\varepsilon$ and before the solution. The goal quantifies over every sequence $\varepsilon_k \to 0$ and assumes no convergence of $U$ and no limit solution. Its part 2 does not assume convergence of $M$ .

The existence of a solution of the planning scheme (Theorem 1) is posed in mission I and not restated here. A solver of Proposition 2 will need it once mission I's goal is proved. The existence and uniqueness of solutions of (20)–(23) are quoted by the paper from Achdou–Capuzzo-Dolcetta (2010) and are not items. Welcome contributions include discrete summation-by-parts lemmas on the torus, the discrete comparison principle for monotone schemes, and the Poincaré-type inequality $\|W\|_\infty \le c\|[D_hW]\|_\infty$ on zero-mean grid functions.

Selected references

Y. Achdou, F. Camilli, I. Capuzzo-Dolcetta, Mean field games: numerical methods for the planning problem, HAL preprint hal-00465404v1, 2010; SIAM J. Control Optim. 50 (2012). https://hal.science/hal-00465404v1, https://doi.org/10.1137/100790069
Y. Achdou, I. Capuzzo-Dolcetta, Mean field games: numerical methods, SIAM J. Numer. Anal. 48 (2010) 1136–1162. https://doi.org/10.1137/090758477
J.-M. Lasry, P.-L. Lions, Mean field games, Japanese Journal of Mathematics 2 (2007) 229–260. https://doi.org/10.1007/s11537-007-0657-8
M. Huang, R. P. Malhamé, P. E. Caines, Large population stochastic dynamic games: closed-loop McKean–Vlasov systems and the Nash certainty equivalence principle, Communications in Information and Systems 6 (2006) 221–252. https://doi.org/10.4310/CIS.2006.v6.n3.a5

11 thms2 active usersReviewed

Completed

Convex Optimization Partial Differential Equations·Captain: mikedeng1

Mean Field Games: Numerical Methods for the Planning Problem I: The Discrete Planning Scheme Has a Solution, Given by a Fenchel–Rockafellar Saddle PointResearch Paper

Motivation

Mean field games (Lasry and Lions, 2006–2007; Huang, Malhamé and Caines, 2006) model the limit of a large population of identical rational agents. Each agent solves an optimal control problem whose cost depends on the distribution $m$ of all agents, and the distribution is in turn transported by the agents' optimal feedback. In the continuous setting this gives a coupled system: a backward Hamilton–Jacobi equation for the value function $u$ and a forward Fokker–Planck equation for the density $m$ .

In the usual formulation, $m$ is prescribed at the initial time and $u$ at the final time. The planning problem, introduced by P.-L. Lions in his Collège de France lectures, prescribes instead both the initial density $m_0$ and the final density $m_T$ , and asks for a cost (through $u$ ) that steers the population from one to the other. According to the paper (§1, pp. 2–3), Lions proved existence for the continuous planning problem in mainly two cases: $\nu = 0$ with a smooth, strictly convex, superlinear Hamiltonian; and $\nu > 0$ with $H(p) = c|p|^2$ or close to it. In both cases the coupling is local and the densities are smooth and bounded away from $0$ . Existence for $\nu > 0$ and more general Hamiltonians was then open, and for sublinear $H$ , $m_0\ne m_T$ and short horizons there is no solution.

Achdou, Camilli and Capuzzo-Dolcetta (hal-00465404, 2010; SIAM J. Control Optim. 2012) introduced a finite-difference scheme for the planning problem and proved that the discrete system has a solution. The proof writes the scheme as the optimality system of a discrete optimal control problem, a Fokker–Planck equation driven by a control, and obtains a solution from a saddle point given by the Fenchel–Rockafellar duality theorem. This mission formalizes that existence result: Theorem 1 of §3.1 together with the lemmas on which its proof rests.

Setting

Fix integers $N_h, N_T \ge 1$ , a horizon $T > 0$ and a viscosity $\nu \ge 0$ ; let $h = 1/N_h$ and $\Delta t = T/N_T$ . The grid $\mathbb T^2_h$ is the periodic $N_h\times N_h$ grid on the two-dimensional torus, with points $x_{i,j}$ , $(i,j)\in(\mathbb Z/N_h)^2$ . On grid functions $U$ the scheme uses the forward differences $D_1^+$ , $D_2^+$ , the four-component discrete gradient $[D_hU]_{i,j} = ((D_1^+U)_{i,j}, (D_1^+U)_{i-1,j}, (D_2^+U)_{i,j}, (D_2^+U)_{i,j-1})$ , the five-point Laplacian $\Delta_h$ , a discrete divergence $\mathrm{div}_h$ of four-component fields, and the transport operator $\mathcal B(U, M) = \mathrm{div}_h(M\,\nabla_q g(\cdot, [D_hU]))$ .

A numerical Hamiltonian $g(x_{i,j}, q_1,q_2,q_3,q_4)$ is monotone (nonincreasing in $q_1, q_3$ , nondecreasing in $q_2, q_4$ ), $C^1$ , convex, and superlinearly coercive in the directions where monotonicity does not bound it. The coupling is local: $V = W'$ , with $W$ strictly convex, superlinear and $C^2$ . The set $\mathcal K$ consists of the discrete probability densities, $h^2\sum_{i,j}M_{i,j} = 1$ , $M\ge 0$ . The discrete planning scheme (18) asks for $(U^n, M^n)_{0\le n\le N_T}$ with

\frac{U^{n+1}-U^n}{\Delta t} - \nu\Delta_hU^{n+1} + g(x,[D_hU^{n+1}]) = V(M^n),\qquad \frac{M^{n+1}-M^n}{\Delta t} + \nu\Delta_hM^n + \mathcal B(U^{n+1},M^n) = 0,

for $0\le n<N_T$ , with $M^n\in\mathcal K$ , $M^0 = m_0$ and $M^{N_T} = m_T$ .

The duality is set up as follows. With $\chi$ the indicator of $\{m\ge0\}$ , let $\Theta(\alpha,\beta) = \sum_{n,i,j}(W+\chi)^*(\alpha^n_{i,j} + g(x_{i,j},[\beta^n]_{i,j}))$ on dual variables $(\alpha^n,\beta^n)_{1\le n\le N_T}$ . Let $\Lambda(\Psi)$ be the linear map sending $\Psi = (\Psi^n)_{0\le n\le N_T}$ to the discrete Hamilton–Jacobi operator and the discrete gradient of $\Psi^{n+1}$ . Let $\Sigma(\alpha,\beta) = \mathcal F(\Psi)$ if $(\alpha,\beta) = \Lambda(\Psi)$ with $\sum\Psi^0 = 0$ , and $+\infty$ otherwise, where $\mathcal F(\Psi) = \frac1{\Delta t}(\sum m_0\Psi^0 - \sum m_T\Psi^{N_T})$ . The Legendre–Fenchel transforms $\Theta^*$ , $\Sigma^*$ act on primal variables $(M^n, Z^n)_{0\le n<N_T}$ , where $M^n$ is paired with $\alpha^{n+1}$ (Remark 2).

Formalization targets

Goal: Theorem 1

Under (G1), (G3)–(G5), (24), $m_0, m_T\in\mathcal K$ , $m_0 > 0$ , and either $\nu > 0$ or ( $\nu = 0$ and $m_T > 0$ ):

\min_{M,Z}\ \Theta^*(M,Z) + \Sigma^*(-M,-Z) = -\min_{\alpha,\beta}\big(\Theta(\alpha,\beta) + \Sigma(\alpha,\beta)\big)

has a solution $(M,Z)$ , $(\alpha,\beta)$ with a finite common value. Moreover $(\alpha,\beta) = \Lambda(U)$ for some $U$ , and $(U, M)$ , with $M^{N_T} = m_T$ appended, solves the scheme (18), with $Z^{k,n} = M^n\,\partial_{q_k}g(x,[D_hU^{n+1}])$ .

Milestones, in attack order

§3.1, p. 7. $V$ maps $(0,\infty)$ onto $(\lambda,\infty)$ ; $(W+\chi)^*$ is finite, convex, continuous and nondecreasing, with explicit values on and off $\mathcal J_V$ .
Lemma 1. $\Theta$ is convex and continuous, $\Sigma$ is convex and l.s.c., and $\Sigma$ and $\Theta$ are both finite at some point.
Lemma 2. $\Theta^*$ and $\Sigma^*$ are convex and l.s.c., with explicit formulas.
(30). $\Sigma^*(-M,-Z)$ is $0$ on the constraint set of the control problem (26) and $+\infty$ off it.
Lemma 3. If $m_0 > 0$ , some $(M,Z)$ has $\Theta^*$ , $\Sigma^*(-M,-Z)$ finite and $\Theta^*$ finite and continuous near it.
Positivity. A discrete strong maximum principle from the proof of Theorem 1: densities in $\mathcal K$ that solve the discrete Fokker–Planck equation (42) are strictly positive before the final time.

Significance

Theorem 1 is the existence result for the finite-difference planning problem. It is used in the paper's second part (§3.2), where solutions of a penalized scheme, in which the final condition is relaxed into a penalty, are shown to converge to a solution of (18) as the penalty parameter vanishes. The penalized scheme is what is solved numerically. The discrete existence result covers general convex, monotone, coercive numerical Hamiltonians and every $\nu\ge0$ , which includes discretizations of the continuous cases that were still open when the paper was written.

The result has a published proof. What this mission adds is a machine-checked version of it: a complete discrete model of the planning problem (grid, operators, scheme, duality functionals) and the convex-analytic chain that connects it to the Fenchel–Rockafellar theorem. As far as a search of the platform shows, no part of this has been formalized; the series' second mission formalizes the convergence of the penalized scheme on the same model.

Difficulty

Existence for a coupled forward–backward nonlinear system with conditions at both ends of the time interval is not reachable by a fixed-point or time-marching argument: the Hamilton–Jacobi equation runs backward and the Fokker–Planck equation forward, and the final density is imposed rather than computed. The approach goes through duality, and three points are delicate.

All the functionals in the duality take the value $+\infty$ , and Fenchel–Rockafellar needs a qualification condition on each side. Lemma 1 gives it for the dual problem; Lemma 3 gives it for the primal one, and needs $m_0 > 0$ and the coercivity (G5).
The optimality conditions only give a complementarity system: the Hamilton–Jacobi equation holds where $M^n > 0$ and becomes an inequality where $M^n = 0$ . Recovering the scheme (18) needs the strict positivity of $M^n$ , a discrete strong maximum principle. That principle uses the monotonicity (G1) and either diffusion or a positive final density.
The bookkeeping of the time lag between primal and dual variables and of the discrete integration by parts in $\Sigma^*$ must be exact.

Formalization scope

The Lean development lives in the namespace MFGPlanning.Existence. A structure Data carries $N_h, N_T \ge 1$ , $T > 0$ , $\nu\ge0$ , $g$ at the grid points, $W$ , $m_0$ and $m_T$ . Committed conventions:

grid indices are ZMod Nh × ZMod Nh (periodic), and $d = 2$ as in the paper;
the components $q_1,\dots,q_4$ and $Z^1,\dots,Z^4$ are the Fin 4 indices 0..3;
time levels are Fin (NT+1) for $U$ and the scheme, and Fin NT for the duality variables. M k is $M^k$ and α k is $\alpha^{k+1}$ (Remark 2), and Fin.snoc M mT appends $M^{N_T} = m_T$ ;
$(W+\chi)^*$ , $\Theta$ , $\Theta^*$ , $\Sigma$ , $\Sigma^*$ are EReal-valued lattice suprema and infima, so unbounded suprema are $+\infty$ and not a junk value. Convexity of an extended-valued functional is convexity of its epigraph.

Disclosed encodings:

(G2) is not encoded, since it only defines the continuous Hamiltonian;
"coercive" in (24) is read as superlinear growth $W(m)/|m|\to\infty$ , which the paper's own consequence $V((0,\infty)) = (\lambda,\infty)$ requires;
$g$ is given only at grid points;
the optimality conditions (32)–(33) are stated through what Theorem 1 says they are equivalent to: the scheme (18) and the relation (39).

Theorem 1 is not reducible to "the scheme (18) has a solution". The goal also asserts that the primal and dual problems attain their minima, that there is no duality gap with a finite value, and that $(\alpha,\beta) = \Lambda(U)$ . The duality functionals are never real-valued suprema, which would make (30) and (31) hold or fail for junk reasons. A complete development needs finite-dimensional convex analysis on extended-valued functions: conjugates, the Fenchel–Rockafellar theorem with attainment, and subdifferential optimality conditions. These parts are reusable well beyond this mission, and proofs of them as separate theorems are welcome.

Selected references

Y. Achdou, F. Camilli, I. Capuzzo-Dolcetta, Mean field games: numerical methods for the planning problem, preprint hal-00465404v1, 2010. https://hal.science/hal-00465404 ; published in SIAM J. Control Optim. 50(1), 2012. https://doi.org/10.1137/100790069
Y. Achdou, I. Capuzzo-Dolcetta, Mean field games: numerical methods, SIAM J. Numer. Anal. 48(3), 2010. https://doi.org/10.1137/090758477
J.-M. Lasry, P.-L. Lions, Mean field games, Japanese Journal of Mathematics 2(1), 2007. https://doi.org/10.1007/s11537-007-0657-8
I. Ekeland, R. Temam, Convex Analysis and Variational Problems, North-Holland, 1976 (SIAM Classics reprint 1999). https://doi.org/10.1137/1.9781611971088

11 thms2 active usersReviewed

Completed

Linear algebra Optimization·Captain: mikedeng1

Projection-like Retractions on Matrix Manifolds V: Within Distance σ_m(X̄) of the Stiefel Manifold, the Polar Factor Is the Unique Nearest Orthonormal FrameResearch Paper

Motivation

Optimization problems with orthonormality constraints arise throughout numerical linear algebra and its applications: computing a few eigenvectors or singular vectors, Procrustes problems in statistics and shape analysis, orthogonal factor rotation, independent component analysis, and the orthogonality constraints of electronic-structure calculations. The feasible set of such a problem is the Stiefel manifold of orthonormal $m$ -frames in $\mathbb R^n$ . Riemannian optimization algorithms on this manifold (Riemannian gradient, Newton and trust-region methods; see Absil, Mahony and Sepulchre, Optimization Algorithms on Matrix Manifolds, 2008) take a step in a tangent direction and then need a retraction: a map that brings the updated point back onto the manifold while agreeing with the geometry to first order.

The most natural way to come back to a constraint set is to take the nearest point. P.-A. Absil and J. Malick, Projection-like retractions on matrix manifolds (SIAM J. Optim. 22(1), 2012; the mission follows the authors' version, HAL hal-00651608v2), show in their Proposition 3.2 that the metric projection onto a smooth submanifold yields a retraction, and then work out the projection explicitly for several matrix manifolds. For the Stiefel manifold, their Proposition 3.4 identifies the projection with a factor of the singular value decomposition, equivalently with the orthonormal factor of the polar decomposition. The paper notes that this result was mentioned without proof in Higham's survey of matrix nearness problems, and gives a proof along the lines of Horn and Johnson, Matrix Analysis, §7.4.

Section 4.5 of the same paper treats a second, "orthographic" retraction on the Stiefel manifold, which corrects a tangent step by a normal vector instead of projecting; for the orthogonal group $\mathbf O_n$ (the case $m=n$ ) it admits a closed form through a matrix square root (Proposition 4.12).

Setting

Fix natural numbers $1\le m\le n$ . Matrices are real, and $\mathbb R^{n\times m}$ carries the Frobenius norm

\|X\|^2=\sum_{i,j}X_{ij}^2=\operatorname{trace}(X^\top X).

The Stiefel manifold is

V_{n,m}=\{X\in\mathbb R^{n\times m}:\ X^\top X=I_m\},

the set of matrices with orthonormal columns; $V_{n,n}$ is the orthogonal group $\mathbf O_n$ .

The singular values of $X$ are $\sigma_1(X)\ge\sigma_2(X)\ge\dots\ge\sigma_{\min\{n,m\}}(X)\ge0$ , the square roots of the eigenvalues of $X^\top X$ . A singular value decomposition of $X$ is a factorization $X=U\Sigma V^\top$ with $U=[u_1,\dots,u_n]\in\mathbf O_n$ , $V=[v_1,\dots,v_m]\in\mathbf O_m$ , and $\Sigma\in\mathbb R^{n\times m}$ zero off its diagonal, with nonnegative nonincreasing diagonal entries. For $\bar X\in V_{n,m}$ every singular value equals $1$ ; in particular $\sigma_m(\bar X)=1$ .

A projection of $X$ onto $V_{n,m}$ is a point $Y\in V_{n,m}$ with $\|X-Y\|\le\|X-Z\|$ for all $Z\in V_{n,m}$ . A polar decomposition of $X$ is a factorization $X=WS$ with $W\in V_{n,m}$ and $S\in\mathbb R^{m\times m}$ symmetric positive definite.

Formalization targets

Goal: Proposition 3.4 (p. 10)

Let $\bar X\in V_{n,m}$ and let $X$ satisfy $\|X-\bar X\|<\sigma_m(\bar X)$ . For every singular value decomposition $X=U\Sigma V^\top$ ,

\{\,Y:\ Y\text{ is a projection of }X\text{ onto }V_{n,m}\,\}=\Big\{\sum_{i=1}^m u_iv_i^\top\Big\},

and $\sum_{i=1}^m u_iv_i^\top$ is the $W$ of the polar decomposition $X=WS$ : it is the orthonormal factor of every polar decomposition of $X$ , and a polar decomposition with this factor exists.

The statement asserts existence, uniqueness and the closed form of the projection on the whole open ball of radius $\sigma_m(\bar X)$ , for whichever singular value decomposition is supplied.

Steps of the proof (milestones)

For $Y\in V_{n,m}$ : $\|X-Y\|^2=\|X\|^2+m-2\operatorname{trace}(Y^\top X)$ .
For every $Y\in V_{n,m}$ , $\operatorname{trace}(Y^\top X)\le\sum_{i=1}^m\sigma_i$ , with equality at $Y=\sum_{i=1}^m u_iv_i^\top$ .
If $\|X-\bar X\|<\sigma_m(\bar X)$ with $\bar X\in V_{n,m}$ , then $X$ has full rank $m$ .
The polar factor of a full-rank matrix is unique (Horn and Johnson, Theorem 7.3.2).

Further items (§4.5)

(S+I)^2=I-\Omega^\top\Omega \tag{4.14}

for $X\in\mathbf O_n$ , $\Omega$ skew-symmetric and $S$ symmetric with $X+X\Omega+XS\in\mathbf O_n$ ; and Proposition 4.12,

R(X,X\Omega)=X\big(\Omega+\sqrt{I-\Omega^\top\Omega}\big),

stated as: $S_+=-I+\sqrt{I-\Omega^\top\Omega}$ is the unique symmetric correction of smallest Frobenius norm.

Significance

Proposition 3.4 makes the projective retraction on the Stiefel manifold computable from one singular value decomposition of an $n\times m$ matrix, and identifies it with the polar factor, which is also what many numerical codes already compute for re-orthonormalization. Combined with Proposition 3.2 of the paper, it yields a second-order-accurate retraction usable in any Riemannian algorithm on $V_{n,m}$ . The same statement is the orthogonal Procrustes problem in the special case of a full-rank target: the nearest orthonormal frame to $X$ .

The result is classical and proved; the paper's argument is short but cites two external facts (Weyl's perturbation bound for singular values and the uniqueness of the polar decomposition). To our knowledge no machine-checked proof of the nearest-orthonormal-frame property, of the uniqueness of the polar factor, or of the closed form of Proposition 4.12 exists in Mathlib or on this platform. A formalization produces reusable pieces: the trace inequality over the Stiefel manifold, the uniqueness of the polar decomposition, and the full-rank property near $V_{n,m}$ , all of which recur in matrix analysis and in the analysis of Riemannian algorithms.

Difficulty

Existence of a nearest point is easy (the Stiefel manifold is compact), and attainment of the trace bound is a direct computation. The content is in two places. First, the bound $\operatorname{trace}(Y^\top X)\le\sum_i\sigma_i$ for every $Y\in V_{n,m}$ requires transporting $Y$ by the orthogonal factors of the singular value decomposition and bounding diagonal entries of a matrix with orthonormal columns. Second, uniqueness does not follow from the trace argument: when $X$ is rank deficient there are many maximizers. Uniqueness needs both a perturbation bound (singular values are 1-Lipschitz in the Frobenius norm, so $\sigma_m(X)>0$ near $\bar X$ ) and the uniqueness of the polar factor, which in turn rests on the uniqueness of the positive-semidefinite square root of $X^\top X$ . Mathlib has singular values of linear maps and the spectral theorem but, at the pinned revision, neither Weyl's inequality for singular values nor the polar decomposition.

For Proposition 4.12, the printed proof compares only two solutions $S_\pm$ of (4.14), while (4.14) has other symmetric solutions (mixed signs of the square roots, or non-diagonal ones on repeated eigenvalues); minimality must be proved against all of them.

Formalization scope

Matrices are Matrix (Fin n) (Fin m) ℝ with the Frobenius norm brought in by open scoped Matrix.Norms.Frobenius. The Stiefel manifold is the set stiefel n m of matrices with Xᵀ * X = 1. Singular values are Mathlib's LinearMap.singularValues of Matrix.toEuclideanLin X, re-indexed to be 1-based as on the page (sv X m is $\sigma_m(X)$ ). A singular value decomposition is the predicate IsSVD X U S V of display (3.5), and $\sum_{i=1}^m u_iv_i^\top$ is U * E * Vᵀ with E the $n\times m$ rectangular identity (frameOfSVD U V). The projection is the set of nearest points, using the published predicate RandomGradFree.Nonsmooth.IsMetricProjection; "exists and is unique" is equality of that set with a singleton. Positive definiteness is Mathlib's Matrix.PosDef, which over $\mathbb R$ includes symmetry.

Standing assumptions and added hypotheses: $m\le n$ is the page's assumption of §3.3; $0<m$ is added so that $\sigma_m$ is meaningful. The radius is $\sigma_m(\bar X)$ with strict inequality, kept in that form although its value is $1$ . Two misprints of the page are corrected in the Lean and kept in the verbatim milestone texts: the display of Proposition 3.4 reads $P_{\mathcal R_r}(X)$ for $P_{V_{n,m}}(X)$ , and the distance identity reads $m^2$ for $m$ . The uniqueness of the polar factor is stated for positive semidefinite factors under the rank hypothesis, which contains the positive definite case of Proposition 3.4. The matrix square root of §4.5 is defined through Mathlib's spectral theorem; Proposition 4.12 is stated for every admissible symmetric correction, and is vacuous only when no symmetric correction exists.

A trivializing formalization is ruled out: the projection is not taken as a hypothesis or chosen by definition; the goal asserts that the nearest-point set equals the singleton of the explicit matrix, for every singular value decomposition supplied, so neither existence nor uniqueness can be assumed away.

Contributions welcome beyond the stated items: Weyl's inequality $|\sigma_i(X)-\sigma_i(Y)|\le\|X-Y\|$ in Mathlib's singular-value API, existence of a singular value decomposition in the matrix form (3.5), and the polar decomposition with its uniqueness, all reusable well beyond this mission.

Selected references

P.-A. Absil and J. Malick, Projection-like retractions on matrix manifolds, SIAM J. Optim. 22(1):135–158, 2012. https://doi.org/10.1137/100802529 ; authors' version https://hal.science/hal-00651608v2
P.-A. Absil, R. Mahony and R. Sepulchre, Optimization Algorithms on Matrix Manifolds, Princeton University Press, 2008. https://doi.org/10.1515/9781400830244
R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press, 1985 (cited by the paper in its 1989 printing) (Theorem 7.3.2, §7.4). https://doi.org/10.1017/CBO9780511810817
N. J. Higham, Matrix nearness problems and applications, in Applications of Matrix Theory (M. J. C. Gover and S. Barnett, eds.), Oxford University Press, 1989, pp. 1–27 (cited by the paper as [15, §4]; no DOI).
A. Edelman, T. A. Arias and S. T. Smith, The geometry of algorithms with orthogonality constraints, SIAM J. Matrix Anal. Appl. 20(2):303–353, 1998. https://doi.org/10.1137/S0895479895290954

10 thms2 active usersReviewed

Completed

Analysis Functional Analysis·Captain: mikedeng1

Theory of Reproducing Kernels VI: The Projection onto the Closed Sum of Two Subspaces as a Series in the Two ProjectionsResearch Paper

Motivation

An orthogonal projection onto a closed subspace is straightforward to describe when that subspace is given directly. It is less straightforward when the subspace is specified as the closed sum of two others: a vector can have contributions from both, and the two component projections generally do not commute. In §12 of Aronszajn's 1950 paper, this problem arises while expressing the reproducing kernel of a sum of two closed subspaces of a reproducing-kernel Hilbert space through their individual kernels. The projection formula is the analytic heart of that calculation. It also gives a series whose finite partial sums can be applied without first describing a basis for the closed sum.

Aronszajn states the result for closed subspaces of a complex Hilbert space of functions. The projection identity itself uses only Hilbert-space geometry; evaluation at a point and the reproducing property enter when the operator formula is translated into a kernel formula later in the section. This mission isolates the projection theorem so that the same formal result can be used both in that kernel calculation and in other settings with two closed subspaces.

Setting

Let $E$ be a complex Hilbert space and let $F_1,F_2$ be closed linear subspaces. The algebraic sum $F_1+F_2$ is the set of all $f_1+f_2$ with $f_i\in F_i$ . It need not be closed. Write $F'=\overline{F_1+F_2}$ for its closure and $F_0=F_1\cap F_2$ for the intersection. All four subspaces are closed except possibly the algebraic sum itself. Let $P_1,P_2,P,P_0$ be the orthogonal projections onto $F_1,F_2,F',F_0$ , respectively. Operator products denote composition, with the rightmost factor applied first; $(P_2P_1)^0$ is the identity.

The paper writes $F_1\oplus F_2$ for $F'$ , even when $F_0$ is nonzero, and uses a distinct dotted plus for the algebraic sum. These symbols can look like a direct sum in newer notation. Here the definitions of $F'$ and $F_0$ remove that ambiguity. Aronszajn assumes complex scalars throughout Part I after §1, p. 343; the statements below keep that convention. Neither a topology on an underlying set of functions nor a measure is part of these operator assertions.

Formalization targets

The goal is the strong projection series, §12, Eq. (7). For each $f\in E$ ,

Pf=P_0f+\sum_{k=1}^{\infty}\left[P_1(P_2P_1)^{k-1}+P_2(P_1P_2)^{k-1}-(P_2P_1)^k-(P_1P_2)^k\right]f.

The sum means that its finite partial sums converge in the norm of $E$ for each fixed $f$ . No operator-norm limit is claimed. The first term of the sum, at $k=1$ , includes $P_1+P_2-P_2P_1-P_1P_2$ because both zero powers are identity operators.

Three source statements are milestones. Equation (1) is a finite identity that expresses a power of $(P-P_1)(P-P_2)$ through a partial sum and a power of $P_2P_1$ . The paragraph following Eq. (4) identifies the strong limit of $(P_2P_1)^m$ as $P_0$ . The paragraph before Eq. (7) asserts that $[(P-P_1)(P-P_2)]^m$ tends strongly to zero. Together these statements specify the finite and limiting parts of the displayed target. They require no assumption that $F_0$ is zero.

Significance

The formula describes the orthogonal projection onto a closed sum in terms of projections onto its two constituents. In Aronszajn's application, an orthogonal projection on a closed subspace determines that subspace's reproducing kernel; therefore the operator expansion supplies a kernel expansion without constructing a basis of the sum. It also makes explicit why an intersection term survives when the subspaces overlap. The result is proved in the 1950 paper; this mission asks for a Lean proof of that known result, with its convergence mode and closure convention made explicit.

A machine-checked development would provide a reusable complex-Hilbert-space statement for alternating projections and closed sums. Mathlib already provides closed submodules and their orthogonal projections, so the new contribution is the finite identity and the strong-limit claims connecting them. The resulting theorem can support the paper's later kernel expression, §12, Eq. (18), once the correspondence between bounded operators and kernels from §11 is available. Equation (18) is outside this mission's target list.

Difficulty

The algebraic sum of two closed subspaces can fail to be closed, so projecting onto it directly would leave the target undefined as a Hilbert-space orthogonal projection. Passing to its closure gives the correct target but does not imply the two projections commute or that their product has a uniform contraction factor. In particular, strong convergence of the powers and of the final series does not generally improve to convergence in operator norm. The later special case in §12, where the intersection is zero and the minimal angle is positive, has a stronger uniform-convergence conclusion; that extra angle condition is absent from Eq. (7).

The finite formula also has four terms per summand whose order matters. Reversing $P_1P_2$ and $P_2P_1$ , beginning the power at $k$ rather than $k-1$ , or omitting $P_0$ changes the statement. Even a proof of pointwise convergence of each selected term would not by itself establish the convergence of the bracketed partial sums in the target.

Formalization scope

The Lean representation is an arbitrary complex Hilbert space E with two ClosedSubmodule ℂ E objects. Mathlib's Submodule.starProjection supplies $P_1$ and $P_2$ . Two small definitions name $P$ as projection onto the join of the closed submodules and $P_0$ as projection onto their meet. The closed-submodule join is the closure of the algebraic sum; the meet is the intersection. These definitions are computed from $F_1,F_2$ , not independently quantified projections. The theorem quantifies over every vector $f$ and uses Filter.Tendsto along natural-number partial sums, giving norm convergence of the resulting vectors.

The source's ambient class has a reproducing kernel, but no step in Eqs. (1)–(7) uses point evaluations. The formal statements therefore apply to any complex Hilbert space, including the paper's RKHS. The stronger scope is explicit here and does not assert a new kernel identity. The $m$ -th Lean partial sum uses indices $0,\ldots,m-1$ for the paper's $1,\ldots,m$ ; it has the same terms and is defined at $m=0$ as $P_0f$ . The finite identity is stated only for $m\ge1$ , as its $m=0$ version is false in general.

Eq. (1) on p. 375 is printed without the opening bracket of the summand (only the closing bracket appears); the formalization reads it as in Eq. (4) on p. 376, where both brackets are printed. There are two printing slips on p. 377: a running sentence drops brackets from the power of $(P-P_1)(P-P_2)$ , and it prints $P\ominus P_2$ while discussing the projection $P-P_2$ . The formalization follows the bracketed expression in Eqs. (1) and (4), using ordinary operator subtraction. The corresponding milestone preserves the printed wording. No assumption from the later angle analysis, especially $F_0=\{0\}$ , is imported into these results. The series is stated through actual partial sums; it is not defined by an arbitrary RKHS chosen to have a desired kernel.

Selected references

N. Aronszajn, Theory of Reproducing Kernels, Transactions of the American Mathematical Society 68 (1950), 337–404, §12, pp. 375–380. DOI: 10.1090/S0002-9947-1950-0051437-7.

6 thms2 active usersReviewed

Completed

Markov Chain Operations Research Probability+1·Captain: mikedeng1

Fundamentals of Queueing Theory X: Uniformization of Continuous-Time Markov ChainsTextbook

Motivation

Most Markovian queueing models have no closed-form transient solution. The M/M/1 queue already needs modified Bessel functions (Chapter 2 of the book), and a finite-capacity or multi-class model with state-dependent rates has no closed form at all. What an analyst can always write down is the system of forward equations $p'(t)=p(t)Q$ for the state probabilities. Chapter 8 of Gross, Shortle, Thompson and Harris, Fundamentals of Queueing Theory (4th ed., Wiley 2008, DOI 10.1002/9781118625651), presents two numerical techniques that turn such models into numbers: the randomization (or uniformization) method for the transient distribution of a finite continuous-time Markov chain, and the Fourier-series method for inverting a Laplace transform, as needed for the M/G/1 waiting-time transform (5.33) and the busy-period transform (5.37).

Uniformization goes back to Jensen (1953) and is the standard transient solver in performance-evaluation and reliability tools. Its appeal is that it replaces a matrix exponential, which is numerically delicate, by powers of a stochastic matrix weighted by Poisson probabilities, with an error bound that can be fixed before the computation starts (Grassmann 1977; Gross and Miller 1984). The Fourier-series method with Euler summation is due to Abate and Whitt (Abate and Whitt 1992; Abate, Choudhury and Whitt 1999).

Setting

A continuous-time Markov chain $X(t)$ on the states $\{0,1,\dots,N\}$ is described by its infinitesimal generator $Q=(q_{ij})$ : for $i\ne j$ , $q_{ij}\ge0$ is the rate of jumps from $i$ to $j$ , and the diagonal entry is $-q_i$ with

q_i=\sum_{j\ne i}q_{ij},\qquad i=0,1,\dots,N.

The transient state-probability vector $p(t)=(p_0(t),\dots,p_N(t))$ , $p_n(t)=\Pr\{X(t)=n\}$ , is the solution of the forward equations

p'(t)=p(t)Q\quad(t\ge0),

started from a given probability vector $p(0)$ . Fix a constant $\Lambda>0$ with $\Lambda\ge q_i$ for every $i$ (the book takes $\Lambda=\max_i q_i$ ) and define the uniformized matrix

\tilde P=\frac{Q}{\Lambda}+I,\qquad \tilde p_{in}=\begin{cases}q_{in}/\Lambda&(i\ne n),\\1-q_i/\Lambda&(i=n).\end{cases}

It is the transition matrix of a discrete-time chain $Y_k$ : the state of $X$ after the $k$ -th event of a Poisson process of rate $\Lambda$ that has been thinned. Write $\phi^{(k)}=p(0)\tilde P^{k}$ for its distribution after $k$ steps.

For the second half of the chapter, the Laplace transform of a real function $f$ on $[0,\infty)$ is $\bar f(s)=\int_0^\infty e^{-st}f(t)\,dt$ , and the Fourier-series approximant with parameter $A$ is

f_{A,n}(t)=\frac{e^{A/2}}{2t}\Big[\bar f\Big(\frac{A}{2t}\Big)+2\sum_{k=1}^{n}(-1)^k\,\mathrm{Re}\,\bar f\Big(\frac{A+2k\pi i}{2t}\Big)\Big],

with $f_A(t)=\lim_{n\to\infty}f_{A,n}(t)$ .

Formalization targets

Goal: the randomization formula with its truncation bound (Eqs. (8.9)–(8.12))

The forward equations have a solution, and every solution satisfies, for all $t\ge0$ ,

p(t)=\sum_{k=0}^{\infty}p(0)\tilde P^{(k)}\,\frac{e^{-\Lambda t}(\Lambda t)^k}{k!},

and whenever $\sum_{k=0}^{T}e^{-\Lambda t}(\Lambda t)^k/k!>1-\epsilon$ , every component of the sum truncated at $k=T$ is within $\epsilon$ of $p_n(t)$ .

Milestones

Eq. (8.12): $\tilde P$ has the entries above and is a stochastic matrix.
Eqs. (8.13)–(8.14): $\phi^{(k)}=\phi^{(k-1)}\tilde P$ and each $\phi^{(k)}$ is a probability vector.
p.385: $\phi=\phi\tilde P\iff0=\phi Q$ .
Eqs. (8.27)–(8.28): for bounded Lipschitz $f$ , $A>0$ and $t>0$ ,

f_A(t)-f(t)=\sum_{k=1}^{\infty}e^{-kA}f\big((2k+1)t\big),\qquad |f_A(t)-f(t)|\le\frac{Ce^{-A}}{1-e^{-A}}\ \text{ if } |f(x)|\le C \text{ for } x>3t.

The mission also contains Eqs. (8.7)–(8.8) as a further theorem, outside the milestone list: the transition probabilities satisfy $p_{in}(t)=\sum_k\tilde p^{(k)}_{in}e^{-\Lambda t}(\Lambda t)^k/k!$ , and $p_n(t)=\sum_i p_i(0)p_{in}(t)$ .

Significance

The randomization formula reduces the transient analysis of any finite Markovian queue (finite-buffer, multi-server, with balking, reneging or state-dependent rates) to repeated vector–matrix products with a sparse stochastic matrix. The truncation point is chosen from a Poisson tail alone, independently of $Q$ . Milestone 3 shows that the same matrix gives the stationary equations, so one iteration serves both transient and steady-state computation. The discretization identity (8.27) is what justifies the parameter choice in Algorithm 8.1: the error decays like $e^{-A}$ .

All of these results are classical and proved in the literature. None of them is formalized in Lean or Mathlib as far as a search of the platform and Mathlib shows. Mathlib has the matrix exponential and Poisson summation under decay hypotheses, but no continuous-time Markov chain generators, no uniformization, and no Laplace transform. This mission would add the finite-state link between generators, stochastic matrices and matrix exponentials that later chapters of queueing and reliability theory use, and a verified error formula for a numerical inversion method in wide use.

Difficulty

The book's derivation is probabilistic: it conditions on the number of events of the Poisson( $\Lambda$ ) process and thins them. A formal statement cannot rest on that picture, because $p(t)$ is defined analytically, by the forward equations. The goal therefore contains a uniqueness statement for a linear ODE on $[0,\infty)$ with one-sided derivative at $0$ , which the book never mentions. The componentwise bound then needs $\tilde P$ to be stochastic, so that every $\phi^{(k)}_n$ lies in $[0,1]$ . That is exactly where $\Lambda\ge\max_i q_i$ is used; with a smaller $\Lambda$ the matrix $\tilde P$ has negative diagonal entries and the bound fails.

For (8.27), the book gives no proof. The identity is an aliasing (Poisson-summation) formula for a periodic function assembled from the values of $f$ at all odd multiples of $t$ . The convergence of the conditionally summed series (8.24) is the delicate point: continuity of $f$ at $t$ , the book's only hypothesis, does not guarantee convergence of a Fourier series. Mathlib's Poisson summation theorems require decay of the Fourier transform that the damped, reflected function built from $f$ does not have.

Formalization scope

States are Fin (N+1); a row vector is Fin (N+1) → ℝ; $pQ$ is vecMul. A generator is a real matrix with nonnegative off-diagonal entries and diagonal $-\sum_{j\ne i}q_{ij}$ .
$p(t)$ is not defined as the series. It is any function with $p(0)=p_0$ and one-sided derivative $p(t)Q$ within $[0,\infty)$ at every $t\ge0$ . The goal also asserts that such a function exists, so it cannot hold vacuously, and it asserts the series identity for every solution. Defining $p(t)$ as the series (8.9) would make the goal a tautology and is ruled out.
$\Lambda$ is any real with $\Lambda>0$ and $\Lambda\ge q_i$ for all $i$ (the book takes equality with $\max_i q_i$ ).
The truncation bound is stated componentwise, as on p.384 ("an error bound on $p_n(t)$ of $\epsilon$ "), for an arbitrary real $\epsilon$ and truncation point $T$ .
The series (8.8), (8.9) are stated with HasSum, so convergence is part of the claim.
The Laplace transform is the Lebesgue integral over $(0,\infty)$ at a complex argument. $f_A(t)$ is the limit of the partial sums $f_{A,n}(t)$ , and the convergence is part of milestone 4.
Strengthened hypotheses in milestone 4: $f$ bounded and Lipschitz on $[0,\infty)$ replaces " $t$ is a continuity point of $f$ ", which is not sufficient for convergence.
Corrected misprints: $e^{-\lambda t}$ in (8.9) is $e^{-\Lambda t}$ ; $q_{ij}/\Lambda$ in (8.12) is $q_{in}/\Lambda$ ; $\phi(Q/\Lambda-I)$ on p.385 is $\phi(Q/\Lambda+I)$ .
Not formalized: Theorem 8.1 (Bromwich inversion) and the real form (8.21), which the book states without hypotheses on $f$ ; the limit claim $\lim_k\phi^{(k)}=\lim_t p(t)$ on p.385, which fails when $\tilde P$ is periodic; the Euler-summation approximation (8.26) and the round-off discussion, which are stated with "≈".

Useful infrastructure: the matrix exponential and its derivative (Matrix, NormedSpace.exp), uniqueness for linear ODEs (Grönwall), Fourier series on the circle, and a reusable Laplace transform file. Contributions of general lemmas on generators and stochastic matrices are welcome, as they apply to every finite Markovian model in the series.

Selected references

D. Gross, J. F. Shortle, J. M. Thompson, C. M. Harris, Fundamentals of Queueing Theory, 4th ed., Wiley, 2008, §§8.1.2–8.2. https://doi.org/10.1002/9781118625651
A. Jensen, "Markoff chains as an aid in the study of Markoff processes", Skandinavisk Aktuarietidskrift 36 (1953) 87–91.
W. K. Grassmann, "Transient solutions in Markovian queueing systems", Computers & Operations Research 4 (1977) 47–53.
D. Gross, D. R. Miller, "The randomization technique as a modeling tool and solution procedure for transient Markov processes", Operations Research 32 (1984) 343–361. https://doi.org/10.1287/opre.32.2.343
J. Abate, W. Whitt, "The Fourier-series method for inverting transforms of probability distributions", Queueing Systems 10 (1992) 5–87. https://doi.org/10.1007/BF01158520
J. Abate, G. L. Choudhury, W. Whitt, "An introduction to numerical transform inversion and its application to probability models", in W. Grassmann (ed.), Computational Probability, Kluwer, 1999, 257–323.

8 thms2 active usersReviewed

Calculus of Variations Control Theory Dynamic Programming+2·Captain: mikedeng1

Bellman's Dynamic Programming VII: Convergence of Discrete Approximations in the Calculus of VariationsTextbook

Motivation

Chapter IX of Richard Bellman's Dynamic Programming (Princeton University Press, 1957; Princeton Landmarks edition 2010, DOI 10.2307/j.ctv1nxcw0f) recasts problems of the calculus of variations with constraints as dynamic programming processes. A variational problem with an inequality constraint on the control, such as $0 \le y \le x$ , is awkward for the classical Euler–Lagrange theory: the optimal control switches between the constraint boundary and the interior, and the number and location of the switches are unknown in advance. Bellman's proposal is to replace the continuous problem by a discrete-time one and to solve that by the recurrence relations of dynamic programming, a procedure he describes as the more reliable computational one "in cases treated to date" (§ 11, p. 260).

That proposal is sound only if the discrete values converge to the continuous value as the time step goes to zero. Chapter IX states one such convergence result, Theorem 2 of § 12, and proves it through an Euler-scheme error estimate. The same chapter works through an example (§§ 10–11) whose discrete version is described by Theorem 1. Both are the subject of this mission. The convergence of time-discretized dynamic programming to the continuous value function has since become a standard topic of numerical optimal control (for example the semi-Lagrangian schemes analysed by Capuzzo-Dolcetta and Falcone), and Bellman's § 12 is an early rigorous instance of it.

Setting

A state $x(t) \ge 0$ evolves on a horizon $[0, T]$ under a control $y(t)$ with $0 \le y \le x$ , a running reward $F(x, y)$ and dynamics

\frac{dx}{dt} = G(x, y), \qquad x(0) = c .

Writing $y = \varphi x$ with a fractional control $0 \le \varphi \le 1$ turns $F$ and $G$ into $\tilde F(x, \varphi) = F(x, \varphi x)$ and $\tilde G(x, \varphi) = G(x, \varphi x)$ (phiForm). The continuous value is

f(c, T) = \sup_{\varphi} \int_0^T \tilde F(x(t), \varphi(t))\, dt,

over measurable $\varphi$ with values in $[0, 1]$ , where $x$ solves $x(t) = c + \int_0^t \tilde G(x(s), \varphi(s))\,ds$ on $[0, T]$ (IsTrajectory, contValue).

For $n = 1, 2, \dots$ the discrete problem uses step $1/n$ and $N = \lfloor Tn \rfloor$ steps (horizonSteps): for $\varphi_0, \dots, \varphi_N \in [0, 1]$ ,

x_0 = c, \quad x_{k+1} = x_k + \frac{\tilde G(x_k, \varphi_k)}{n}, \qquad J_N(\{\varphi_k\}, n) = \sum_{k=0}^{N} \frac{\tilde F(x_k, \varphi_k)}{n},

and $f(c, T, n) = \max J_N$ (eulerTraj, discretePayoff, discreteValue). A discrete control defines the step control $\varphi(t) = \varphi_k$ on $k/n \le t < (k+1)/n$ (stepControl).

The example of §§ 10–11 has a gain function $b$ with $b(0) = 0$ , $b'(0) = \infty$ , $b' > 0$ , $b'(y) \to 0$ as $y \to \infty$ , $b'' < 0$ (IsGainFunction; $b(y) = y^{1/2}$ is one), and value functions (uSeq)

u_0(c) = c, \qquad u_{N+1}(c) = \max_{0 \le v \le c}\,\bigl[c - v + u_N(c + b(v))\bigr].

Formalization targets

Goal: Chapter IX, Theorem 2 (corrected)

Assume (11): (a) $F$ and $G$ have continuous second partial derivatives; (b) $px \le G(x, y) \le qx + r$ for $x > 0$ , $0 \le y \le x$ ; (c) $G_y > 0$ throughout that region, or $G_y < 0$ throughout. Then for all $c \ge 0$ , $T > 0$ , the set of continuous payoffs is nonempty and bounded above, and

\lim_{n \to \infty} f(c, T, n) = f(c, T).

Milestones

Chapter IX, Theorem 1: the structure of $u_N$ . There are $v_N(c)$ and thresholds $c_N$ with $v_N$ decreasing in $c$ , $v_{N+1} > v_N$ , $c_N$ the unique fixed point of $v_N$ with $c_{N+1} > c_N$ , $u_N(c) = u_{N-1}(c + b(c))$ for $c \le c_N$ , $u_N(c) = c - v_N(c) + u_{N-1}(c + b(v_N(c)))$ for $c \ge c_N$ , and $u_N' \ge u_{N-1}'$ .
§ 12, Lemma (corrected): for $G$ Lipschitz on $[m, M] \times [0, 1]$ , the Euler states with a step control are within $\kappa / n$ of the exact solution on $[0, T]$ .
Eq. (12.13): $|J(\varphi) - J_N(\{\varphi_k\}, n)| \le B'/n$ for step controls.
Eq. (12.14): $f(c, T, n) \le f(c, T) + B'/n$ for all $n \ge 1$ .
Eq. (12.17): $f(c, T) \le \liminf_{n \to \infty} f(c, T, n)$ .

Milestones 3–5 are stated for $c > 0$ , as in the book's proof ("Given $c > 0$ and $T > 0$ "); the goal is stated for $c \ge 0$ , as in the theorem.

Significance

Theorem 2 says that the value of a constrained continuous-time control problem can be computed, to any accuracy, by the finite recurrence of dynamic programming on a time grid. Its upper half (12.14) gives a rate: the discrete value never exceeds the continuous one by more than $B'/n$ . The lower half (12.17) needs no rate and holds because measurable controls are approximated by step controls. Theorem 1 is a discrete counterpart of the transition curve computed in § 10: below the threshold $c_N$ the whole state is invested, above it an interior amount, and the thresholds increase with the number of remaining stages.

As far as the mission's author could establish, none of these results has a machine-checked proof. The Euler error estimate is classical and has a Gronwall-type proof; Mathlib contains Gronwall's inequality (Analysis/ODE/Gronwall) and Picard–Lindelöf for continuous right-hand sides, but no convergence theory for the value of discretized control problems. The book leaves the proof of Theorem 1 to the reader.

Difficulty

The upper bound (12.14) follows from the Lemma once the continuous and the discrete trajectories are known to stay in a common bounded strip $m \le x \le M$ ; that uniform bound is what assumption (11b) provides and must be established first. The lower bound is where the obvious argument fails: an arbitrary measurable control is not a step control on the grid $k/n$ , and passing to a step control changes the trajectory, so the payoff must be shown continuous under almost-everywhere convergence of controls. This uses the trajectory's dependence on the control in $L^1$ , not only on the initial value. Solutions exist in the Carathéodory sense only, so the integral form of the equation is required. For Theorem 1 the induction must carry concavity of $u_N$ and a strict comparison of marginal values across $N$ , and $u_N$ is not differentiable at $c = 0$ , where $b'(0) = \infty$ .

Formalization scope

States, controls and rewards are real; $F, G, b : \mathbb R \to \mathbb R$ (or $\mathbb R \to \mathbb R \to \mathbb R$ ), with the hypotheses imposed where the book imposes them. The conventions:

Misprint in (12.5). The print sets $N = [T/n]$ while using the step $1/n$ in (12.6) and the intervals $k/n \le t < (k+1)/n$ in the Lemma. Under the literal reading the discrete horizon $N/n$ tends to $0$ and Theorem 2 is false: for $F \equiv 1$ , $G(x, y) = y$ (which satisfies (11) with $p = 0$ , $q = 1$ , $r = 0$ , $G_y = 1$ ) and $T = 1$ , $f(c, 1) = 1$ but $f(c, 1, n) = 1/n$ for $n \ge 2$ . The mission uses $N = \lfloor Tn \rfloor$ and keeps the book's sum $\sum_{k=0}^{N}$ ; the range " $k = 0, \dots, n - 1$ " in (12.6) is read as $k = 0, \dots, N - 1$ .
Lemma. Its range " $0 \le t \le N$ " is read as $0 \le t \le T$ . The bounds $m \le x \le M$ are required of the solution $x(t)$ as well as of the sequence $x_k$ , since a Lipschitz condition on the strip says nothing about $G$ outside it; the proof of Theorem 2 supplies both bounds. The constant may depend on $m$ , $M$ and the Lipschitz constant.
Assumptions (11) are on the original $F(x, y)$ , $G(x, y)$ , with (11a) read as $C^2$ on $\mathbb R^2$ ; the problem is posed through $y = \varphi x$ .
"Max" over measurable controls is a supremum (the proof picks $\varepsilon$ -optimal controls). The goal asserts that the set of payoffs is nonempty and bounded above, so the Lean supremum cannot take its default value $0$ . Trajectories satisfy the integral equation with integrable right-hand side, and the reward is required integrable, so the Bochner integral's default $0$ cannot enter either.
Theorem 1. (5a) is read as non-strict monotonicity: for $N = 1$ the interior optimum solves $b'(v) = 1$ and is constant in $c$ . $v_N(c) \ge 0$ is required for $c \ge c_N$ , so that $v_N(c)$ is a feasible choice. (5f) presupposes differentiability and is stated where both derivatives exist.

A trivializing formalization is ruled out: the discrete and continuous values are defined from $F$ and $G$ by the book's recursions and integrals, not assumed as hypotheses, and the corrected step count is not a free parameter.

Needed infrastructure: Carathéodory existence and uniqueness for Lipschitz right-hand sides with measurable controls, Gronwall estimates for the Euler scheme, and approximation of measurable controls by step functions. All three are reusable beyond this mission, and contributions of any of them are welcome.

Selected references

R. Bellman, Dynamic Programming, Princeton University Press, 1957; Princeton Landmarks in Mathematics edition, 2010. Chapter IX, §§ 10–12, pp. 256–263. DOI 10.2307/j.ctv1nxcw0f
I. Capuzzo-Dolcetta, On a discrete approximation of the Hamilton–Jacobi equation of dynamic programming, Applied Mathematics and Optimization 10 (1983) 367–377. DOI 10.1007/BF01448394
M. Bardi and I. Capuzzo-Dolcetta, Optimal Control and Viscosity Solutions of Hamilton–Jacobi–Bellman Equations, Birkhäuser, 1997. DOI 10.1007/978-0-8176-4755-1

8 thms2 active usersReviewed

Completed

Convex Optimization Operations Research Optimization·Captain: mikedeng1

A New Projection Method for Variational Inequality Problems: The Hyperplane Projection Method Converges to a Solution Under Continuity and Generalized MonotonicityResearch Paper

Motivation

A variational inequality asks for a point of a convex set at which a vector field points "inward" against every feasible direction. The format covers the first-order optimality conditions of constrained optimization, nonlinear complementarity problems, traffic and economic equilibria (Wardrop, Walrasian, Nash–Cournot), and systems of nonlinear equations; see Harker and Pang's survey (Math. Programming 48, 1990) and Facchinei and Pang's monograph (Springer, 2003).

When the map has no special structure (not strongly monotone, not Lipschitz with known constant, not affine) and the feasible set is a general closed convex set, the practical algorithms are projection methods. The oldest is Korpelevich's extragradient method (1976). Without a known Lipschitz constant, extragradient-type methods need a linesearch in which every trial point costs one projection onto the feasible set, and projection onto a general convex set is itself an optimization problem.

Solodov and Svaiter (SIAM J. Control Optim. 37 (1999) 765–776) proposed a method that spends exactly two projections per iteration, whatever the linesearch does, and proved global convergence under only continuity of the map and a generalized monotonicity condition weaker than pseudomonotonicity. The method, often called the hyperplane projection method, is a standard reference point for later projection and extragradient-type algorithms.

Timeline:

1976: Korpelevich, extragradient method, Lipschitz monotone maps.
1987–1994: Khobotov (1987), Iusem (1994) and others: extragradient variants with Armijo-type stepsize rules, which need one projection per trial step.
1997: Iusem and Svaiter, a separating-hyperplane variant of extragradient for monotone maps (reference [9] of the paper).
1999: Solodov and Svaiter, Algorithm 2.1: two projections per iteration, convergence under condition (1.2) below.

Setting

Work in $\mathbb{R}^n$ with the Euclidean inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$ . Let $C \subseteq \mathbb{R}^n$ be closed and convex and $F : \mathbb{R}^n \to \mathbb{R}^n$ continuous. The problem $\mathrm{VI}(F, C)$ is to find $x^*$ with

x^* \in C, \qquad \langle F(x^*), x - x^*\rangle \ge 0 \quad \text{for all } x \in C. \tag{1.1}

Its solution set is $S$ . The projection onto a nonempty closed convex set $K$ is $P_K[x] := \arg\min_{y \in K}\|y - x\|$ . The projected residual is $r(x) := x - P_C[x - F(x)]$ ; its zeros are exactly the points of $S$ .

Condition (1.2) requires, for every $x^* \in S$ ,

\langle F(x), x - x^*\rangle \ge 0 \qquad \text{for all } x \in C. \tag{1.2}

It holds when $F$ is monotone or pseudomonotone, and in cases where $F$ is neither.

Algorithm 2.1. Fix $\gamma, \sigma \in (0,1)$ and $x^0 \in C$ . Given $x^i$ : if $r(x^i) = 0$ , stop. Otherwise let $k_i$ be the smallest nonnegative integer $k$ with

\langle F(x^i - \gamma^k r(x^i)), r(x^i)\rangle \ge \sigma\|r(x^i)\|^2, \tag{2.1}

set $\eta_i = \gamma^{k_i}$ , $z^i = x^i - \eta_i r(x^i)$ , $H_i = \{x \mid \langle F(z^i), x - z^i\rangle \le 0\}$ , and

x^{i+1} = P_{C \cap H_i}[x^i].

The hyperplane $\partial H_i$ separates $x^i$ from $S$ .

Formalization targets

Goal: Theorem 2.1

If $C$ is closed and convex, $F$ is continuous, $S \ne \emptyset$ and (1.2) holds, then every sequence generated by Algorithm 2.1 converges to a single point of $S$ :

\exists\, \hat x \in S:\quad x^i \to \hat x \quad (i \to \infty).

The theorem fixes no rate and no constant; it asserts convergence of the whole sequence, not only of a subsequence.

Milestones (in attack order)

Lemma 2.1 (p. 768): for nonempty closed convex $B$ , $\langle x - P_B[x], z - P_B[x]\rangle \le 0$ for $z \in B$ , and $\|P_B[x]-P_B[y]\|^2 \le \|x-y\|^2 - \|P_B[x]-x+y-P_B[y]\|^2$ .
Residual characterization (p. 767): $x \in S \iff r(x) = 0$ .
(2.5) (p. 769): $\langle F(x), r(x)\rangle \ge \|r(x)\|^2$ for $x \in C$ .
Linesearch well-definedness (p. 769): for $x \in C$ with $r(x) \ne 0$ , some $k$ satisfies (2.1).
Lemma 2.2 (p. 768): $x^{i+1} = P_{C\cap H_i}[\bar x^i]$ with $\bar x^i = P_{H_i}[x^i]$ .
(2.6) (pp. 769–770): $\|x^{i+1}-x^*\|^2 \le \|x^i-x^*\|^2 - \|x^{i+1}-\bar x^i\|^2 - \big(\eta_i\sigma/\|F(z^i)\|\big)^2\|r(x^i)\|^4$ for every $x^* \in S$ .
(2.8) (p. 770): $\eta_i\|r(x^i)\| \to 0$ .

Significance

Theorem 2.1 gives global convergence of a projection method for variational inequalities with no Lipschitz constant, no monotonicity and no knowledge of the problem beyond continuity and (1.2), at a fixed cost of two projections per iteration. Condition (1.2) covers pseudomonotone maps, which arise as gradients of pseudoconvex functions and in equilibrium models where monotonicity fails. The separating-hyperplane-and-project template of the proof is reused throughout the later literature on projection, proximal and hybrid methods for monotone inclusions.

The result has been proved since 1999. To the best of available knowledge no machine-checked proof of it, or of any convergence theorem for a projection method for variational inequalities, exists in Lean or Mathlib. This mission produces the statement and the supporting layer: a Euclidean projection onto closed convex sets with its standard inequalities, variational inequality solution sets, the projected residual, and a formal model of an Armijo-type linesearch algorithm with termination.

Difficulty

The Fejér-type inequality (2.6) quickly gives bounded iterates and $\eta_i\|r(x^i)\| \to 0$ . The obvious next step, concluding $r(x^i) \to 0$ , fails: nothing prevents the stepsizes $\eta_i$ from tending to zero, and in that regime the product going to zero says nothing about the residual. This regime is where the minimality of $k_i$ and the continuity of $F$ enter, and it is the step a naive formalization (for instance one that drops minimality, or fixes the stepsize) cannot reach. A second subtlety is that (1.2) is needed at an accumulation point that is only known to lie in $S$ at the end of the argument, which is why the condition must hold for every $x^* \in S$ . Finally, subsequential convergence must be upgraded to convergence of the whole sequence to one solution; convergence of a subsequence, or of the distance to $S$ , is strictly weaker.

Formalization scope

$\mathbb{R}^n$ is EuclideanSpace ℝ (Fin n) (not Fin n → ℝ, whose norm is the sup norm). The accumulation-point step needs finite dimension; no Hilbert-space generalization is intended.
Projection encoding. projOnto K x is a nearest point of $K$ to $x$ when one exists, chosen by Classical.choose, and the junk value $x$ otherwise. On nonempty closed convex sets it is exactly $P_K[x]$ ; the paper only projects onto such sets ( $C$ , $H_i$ , $C \cap H_i$ ), so the junk value is never reached under the hypotheses.
Stopping-rule encoding. A run is a sequence x : ℕ → ℝⁿ with Armijo indices k : ℕ → ℕ (predicate IsAlg21Run). If $r(x^i) = 0$ the method has stopped and the run stalls, $x^{i+1} = x^i$ ; otherwise $k_i$ is the least index satisfying (2.1) and $x^{i+1} = P_{C \cap H_i}[x^i]$ . A stalled point is a solution, so finitely terminating runs are included in the goal.
Parameters $\gamma, \sigma$ are real with $0 < \gamma < 1$ , $0 < \sigma < 1$ , universally quantified; $n$ , $C$ , $F$ and $x^0 \in C$ are arbitrary.
(2.6) is stated for one generic step ( $x \in C$ , $r(x)\ne 0$ , $k$ satisfying (2.1)) rather than along a run; it is the same inequality with $x^i, k_i$ abstracted.
Trivializing formalizations are ruled out: condition (1.2) is quantified over every solution and every $x \in C$ (not replaced by monotonicity or an existential), $k_i$ is the least index satisfying (2.1), the update projects $x^i$ onto $C \cap H_i$ (not onto $C$ alone), the stopped case is pinned down by the stall encoding, and the conclusion is convergence of the whole sequence to one solution, not $r(x^i) \to 0$ or $\operatorname{dist}(x^i, S) \to 0$ .
Needed infrastructure: existence, uniqueness and variational characterization of the projection (Mathlib has exists_norm_eq_iInf_of_complete_convex and norm_eq_iInf_iff_real_inner_le_zero), firm nonexpansiveness, the explicit projection onto a halfspace, and a bounded-sequence subsequence argument in $\mathbb{R}^n$ . The projection lemmas are reusable for any projection-type method; contributions proving them as standalone lemmas are welcome.

Selected references

M. V. Solodov and B. F. Svaiter, A New Projection Method for Variational Inequality Problems, SIAM J. Control Optim. 37(3), 765–776, 1999. https://doi.org/10.1137/S0363012997317475
G. M. Korpelevich, The extragradient method for finding saddle points and other problems, Matecon 12, 747–756, 1976.
A. N. Iusem and B. F. Svaiter, A variant of Korpelevich's method for variational inequalities with a new search strategy, Optimization 42, 309–321, 1997. https://doi.org/10.1080/02331939708844365
P. T. Harker and J.-S. Pang, Finite-dimensional variational inequality and nonlinear complementarity problems: a survey of theory, algorithms and applications, Math. Programming 48, 161–220, 1990. https://doi.org/10.1007/BF01582255
F. Facchinei and J.-S. Pang, Finite-Dimensional Variational Inequalities and Complementarity Problems, Springer, 2003. https://doi.org/10.1007/b97543

16 thms2 active usersReviewed

Completed

Linear algebra Operations Research Optimization·Captain: mikedeng1

Conjugate Gradient Methods with Inexact Searches: The Self-Scaled Direction Is a Multiple of Beale's Restart DirectionResearch Paper

Motivation

Conjugate gradient methods minimize a smooth function $f:\mathbb R^n\to\mathbb R$ using only gradients and a handful of stored vectors. This makes them the standard choice when $n$ is too large for Newton or quasi-Newton methods, which store an $n\times n$ matrix. On a strictly convex quadratic with exact line searches the classical method of Hestenes and Stiefel terminates in at most $n$ steps. On general functions, and with the inexact line searches used in practice, its behaviour is much less clear.

D. F. Shanno's 1978 paper in Mathematics of Operations Research (doi:10.1287/moor.3.3.244) links conjugate gradient methods to quasi-Newton methods. It writes the search direction as $-\hat H g$ , where $\hat H$ is a positive definite approximation of the inverse Hessian that is never stored. The resulting "memoryless" BFGS directions give descent without exact line searches. The paper's new algorithm uses two BFGS updates: one from the last restart and one from the current step. Its first update is scaled by the Oren–Spedicato factor $\gamma_t$ . Shanno and Phua's CONMIN code implements the algorithm, and the memoryless BFGS direction is the one-pair case of the later limited-memory BFGS methods.

1952: Hestenes and Stiefel, linear conjugate gradients.
1964: Fletcher and Reeves, nonlinear conjugate gradients.
1969: Polak and Ribière, a second nonlinear variant.
1972: Beale, a restart procedure that keeps the computed direction $d_t$ .
1977: Powell's restart criterion (Powell 1977).
1978: Shanno's reformulation as memoryless and two-update quasi-Newton methods (this paper).

Setting

Vectors are columns in $\mathbb R^n$ . A prime denotes transpose: $u'v$ is the inner product and $uv'$ the outer product. An iterative method produces points $x_k$ , steps $p_k = x_{k+1}-x_k = \alpha_k d_k$ along search directions $d_k$ , gradients $g_k = \nabla f(x_k)$ , and gradient changes $y_k = g_{k+1}-g_k$ . A line search is exact when $p_k'g_{k+1} = 0$ .

The BFGS update of a matrix $H$ with the pair $(p,y)$ is

H^+ = H - \frac{p\,y'H + H y\,p'}{p'y} + \left(1+\frac{y'Hy}{p'y}\right)\frac{pp'}{p'y}.

A restart cycle begins at iteration $t$ . At a later iteration $k>t$ , Shanno's self-scaled restart matrix is

\hat H_k = \gamma_t\left(I - \frac{p_ty_t'+y_tp_t'}{p_t'y_t} + \frac{y_t'y_t}{p_t'y_t}\frac{p_tp_t'}{p_t'y_t}\right) + \frac{p_tp_t'}{p_t'y_t}, \qquad \gamma_t = \frac{p_t'y_t}{y_t'y_t}.

The matrix $\hat H_{k+1}$ is its BFGS update with $(p_k,y_k)$ , and the self-scaled two-update direction is $d_{k+1} = -\hat H_{k+1}g_{k+1}$ . The unscaled variant uses the BFGS update of $I$ in place of the first matrix.

Beale's direction is

d_{k+1} = -g_{k+1} + \frac{y_k'g_{k+1}}{d_k'y_k}d_k + \frac{y_t'g_{k+1}}{d_t'y_t}d_t.

The quadratic case has gradient $g(x) = Ax + c$ with $A$ symmetric positive definite.

Formalization targets

Goal: reduction of the self-scaled method to Beale's method

Let $A$ be symmetric positive definite, $g_i = Ax_i + c$ and $t<k$ . Assume that for $t\le i\le k$ we have $x_{i+1} = x_i + p_i$ , $p_i = \alpha_i d_i$ and $p_i'g_{i+1}=0$ , that $p_t'Ap_i = 0$ for $t<i\le k$ , and that $p_t, p_k \ne 0$ . Then

-\hat H_{k+1}g_{k+1} = \gamma_t\left(-g_{k+1} + \frac{y_k'g_{k+1}}{d_k'y_k}d_k + \frac{y_t'g_{k+1}}{d_t'y_t}d_t\right).

This is the paper's claim that "for $f(x)$ quadratic with exact searches each of the above methods reduces exactly to Beale's method defined by (28)", with the conclusion (44). The scale is exactly $\gamma_t$ .

Companion statements

The unscaled two-update direction equals Beale's direction exactly.
Both two-update directions are descent directions, $g_{k+1}'d_{k+1} < 0$ , whenever $p_t'y_t > 0$ and $p_k'y_k > 0$ . No exact search is needed.

Milestones on the path

(34): the expansion of $-\hat H_{k+1}g_{k+1}$ .
(40): its form under an exact search.
(41): gradients along a run on a quadratic.
$p_t'g_{k+1} = 0$ .
(38), corrected by a factor 2: the action of the self-scaled restart matrix.
(42): its form when $p_t'g_{k+1} = 0$ .
(43): the direction after substitution.

Significance

The result. The reduction says that the new algorithm reproduces Beale's restarted conjugate gradient directions on a quadratic with exact line searches. So it keeps the finite-termination and rate-of-convergence properties behind Beale's restart. Away from that setting it behaves as a quasi-Newton method, whose directions are descent directions under any line search with $p'y>0$ . The two regimes are what justify relaxing the line search, which the paper's computations exploit. The scale $\gamma_t$ changes only the length of the step, not its direction.

Formalizing it. The claim is proved in the paper by a short computation, and no machine-checked version is known. This mission produces:

a checked statement of the claim with every hypothesis explicit, including the conjugacy the proof takes as known;
a corrected version of display (38), which is misprinted;
reusable definitions of the additive BFGS update and of Beale's direction.

Difficulty

The obvious attempt is to expand both BFGS updates symbolically and compare with Beale's formula. This fails without two facts that are not algebraic identities. The first is that the restart step stays orthogonal to every later gradient, $p_t'g_{k+1}=0$ . It needs the affine gradient of a quadratic, the exact search at the restart step, and conjugacy of $p_t$ with all later steps. The second is $y_k'p_t = 0$ , which again comes from conjugacy. Beale's formula also has to be matched in its $d$ -form: the coefficient $\frac{y'g}{d'y}d$ equals $\frac{y'g}{p'y}p$ only when the step length is nonzero. The descent statements need a different argument: the BFGS update of a positive definite matrix with $p'y>0$ must be shown to remain positive definite, and this has to be done twice.

Formalization scope

Vectors are Fin n → ℝ and matrices Matrix (Fin n) (Fin n) ℝ. The inner product $u'v$ is u ⬝ᵥ v, the outer product $uv'$ is vecMulVec u v, and $Hv$ is H *ᵥ v. The quadratic enters only through its gradient A *ᵥ x + c with A.PosDef. The paper's (4) is the case $c = -A\hat x$ . Iterates, steps, directions and gradients are sequences indexed by ℕ. Division is Lean's total division. Every statement that divides therefore carries hypotheses making its denominators nonzero: $p_t \ne 0$ and $p_k\ne0$ in the quadratic statements, and $p_t'y_t\ne 0$ or $p'y>0$ in the generic ones.

The conjugacy $p_t'Ap_i=0$ for $t<i\le k$ is a hypothesis, exactly as the paper's proof uses it. It is not derived from a full run of Beale's algorithm. The range $k\le t+n-1$ of Beale's formula is not assumed.

Several formalizations would make the claim easier than the paper's, and none of them is used:

defining the direction by the expanded formula (34), or the restart matrix by (38);
adding orthogonality or conjugacy hypotheses beyond those listed;
concluding only that the two directions are parallel;
dropping the restart term of Beale's direction;
allowing a zero denominator.

The generic milestones — (34), (40), (38), (42), (43) and the descent statements — are statements about arbitrary vectors and matrices and are reusable for any BFGS-based method. Proofs of any milestone, and alternative derivations of the goal, are welcome.

Selected references

D. F. Shanno, Conjugate Gradient Methods with Inexact Searches, Mathematics of Operations Research 3(3) (1978) 244–256. https://doi.org/10.1287/moor.3.3.244
E. M. L. Beale, A derivation of conjugate gradients, in F. A. Lootsma (ed.), Numerical Methods for Nonlinear Optimization, Academic Press, 1972, 39–43.
M. J. D. Powell, Restart procedures for the conjugate gradient method, Mathematical Programming 12 (1977) 241–254. https://doi.org/10.1007/BF01593790
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
S. S. Oren and E. Spedicato, Optimal conditioning of self-scaling variable metric algorithms, Mathematical Programming 10 (1976) 70–90. https://doi.org/10.1007/BF01580654

18 thms2 active usersReviewed

Completed

Convex Optimization Optimization·Captain: mikedeng1

The Relaxation Method for Linear Inequalities III: Reflexion in a Closed Bounded Convex Set Terminates or Ends in OscillationResearch Paper

Motivation

The relaxation method for a system of linear inequalities, introduced by Agmon and by Motzkin and Schoenberg in back-to-back papers of the Canadian Journal of Mathematics (1954), solves $\sum_j a_{ij}x_j + b_i \ge 0$ by repeatedly moving a point towards, or across, the most violated half-space. It is the ancestor of the perceptron algorithm, of Kaczmarz-type projection methods, and of the method of alternating projections, all of which are still used in feasibility problems, tomography and machine learning.

Motzkin and Schoenberg's paper ends (Part IV, §§9–10) by asking whether the behaviour of the reflexion process — the relaxation step with factor $\lambda = 2$ — survives when the finite family of half-spaces is replaced by an infinite one. They answer this for one natural infinite family: all supporting half-spaces of a closed bounded convex set. This mission formalizes that answer, Theorem 3 of the paper.

Timeline of the thread this mission belongs to:

1922 — Fejér observes that a sequence approaching every point of a set monotonically has useful convergence properties (Fejér, Math. Annalen 85, 1922).
1954 — Agmon proves convergence of the relaxation method for $0 < \lambda < 2$ (Agmon, Canad. J. Math. 6, 1954, pp. 382–392); Motzkin and Schoenberg prove finite termination of the reflexion method ( $\lambda = 2$ ) for full-dimensional solution polytopes (Theorem 1), the oscillation behaviour in lower dimension (Theorem 2), and the convex-body version (Theorem 3) (Motzkin–Schoenberg 1954).

Setting

Let $E_n$ be $n$ -dimensional Euclidean space and let $A \subseteq E_n$ be a nonempty, closed, bounded, convex set. Its dimension $r$ is the dimension of its affine span $L_r$ , the smallest flat containing $A$ .

For $p \notin A$ let $q$ be the point of $A$ nearest to $p$ ; it exists and is unique. The image of $p$ with respect to $A$ is

p_1 = F(p) = p + 2(q - p), \qquad (3.1)

the reflexion of $p$ through $q$ . The reflexion process starts at $p_0 \notin A$ and sets $p_{\nu+1} = F(p_\nu)$ as long as $p_\nu \notin A$ (3.2). Either the process terminates with some $p_N \in A$ , or it produces an infinite sequence of points outside $A$ .

The family $F$ of supporting half-spaces of $A$ consists of the closed half-spaces $H \supseteq A$ whose bounding hyperplane touches $A$ . The paper observes that the half-space $H_0$ through $q$ normal to $pq$ is the member of $F$ farthest from $p$ , so (3.1) is exactly the reflexion step of the relaxation method applied to the infinite family $F$ .

A sequence $\{q_\nu\}$ of points outside $A$ is Fejér-monotone with respect to $A$ if $q_\nu \ne q_{\nu+1}$ and $|q_{\nu+1} - a| \le |q_\nu - a|$ for all $a \in A$ . Two points $u, v$ are symmetric with respect to a flat $L$ if their midpoint lies in $L$ and $u - v$ is orthogonal to $L$ .

Formalization targets

Goal: Theorem 3 (p. 402)

For every nonempty closed bounded convex $A \subseteq E_n$ with affine span $L_r$ , every $p_0 \notin A$ and every run $\{p_\nu\}$ of the reflexion process:

\textbf{Case 1: } r = n \implies \exists N,\ p_N \in A.

\textbf{Case 2: } r < n \implies \begin{cases} p_0 \in L_r \implies \exists N,\ p_N \in A,\\[2pt] p_0 \notin L_r \implies p_\nu \notin A\ \forall \nu, \text{ and } \exists \nu_0, u \ne v \text{ symmetric w.r.t. } L_r:\ \{p_\nu, p_{\nu+1}\} = \{u, v\}\ \forall \nu > \nu_0. \end{cases}

No number of steps is fixed: termination is finite but not uniformly bounded.

Milestones (attack order)

§9 — the half-space $H_0$ belongs to $F$ , maximizes $\operatorname{dist}(p, H)$ over $F$ , and every farthest member of $F$ yields the step (3.1).
§10 — an infinite run of the reflexion process is Fejér-monotone with respect to $A$ .
Lemma 1, Case 1 — a sequence Fejér-monotone with respect to a set of dimension $n$ converges to a point.
§10 — the limit of an infinite run lies in $A$ and on its boundary.
Theorem 3, Case 1 — if $r = n$ the process always terminates.
§10 — orthogonal projection on a flat $L \supseteq A$ commutes with the image map, and each step keeps the distance to $L$ while switching sides.

Significance

The result. Theorem 3 shows that the dichotomy proved in the paper for finitely many half-spaces — finite termination when the target is full-dimensional, eventual oscillation otherwise — holds for the infinite family of all supporting half-spaces of a convex body. Case 1 says that reflecting through the nearest point, a method that uses no information about $A$ beyond metric projection, reaches a full-dimensional convex body in finitely many steps from any start. Case 2 says that for a lower-dimensional body the process detects this: the iterates settle into a two-cycle whose midpoint is a point of $A$ , so a solution can be read off.

Formalizing it. The theorem is proved in the paper (1954), with a short proof of Case 1 that relies on geometric intuition about normal cones near a boundary point. To our knowledge no machine-checked proof exists. A complete development produces a verified finite-termination theorem for a projection method on general convex bodies, together with reusable facts about Fejér-monotone sequences and metric projections that recur throughout the analysis of projection algorithms.

Difficulty

The obvious argument for Case 1 is: the run is Fejér-monotone, hence converges (Lemma 1), and its limit $a$ lies on the boundary of $A$ ; then derive a contradiction. For a polytope (Theorem 1) the contradiction comes from finiteness: eventually every reflexion is in one of finitely many hyperplanes through $a$ , which keeps the iterates on a sphere around $a$ . For a convex body there are infinitely many supporting hyperplanes near $a$ , and the iterates are reflected in a different one at every step; no finiteness argument is available. What replaces finiteness is an argument about how the supporting hyperplanes of $A$ at boundary points near $a$ are oriented, and the paper gives it only as an informal geometric sketch. Making this step rigorous is the main work of the mission.

Numerical experiments show that termination can take tens of thousands of steps near sharp corners of a polygon, so no bound on the number of steps in terms of the distance of $p_0$ to $A$ alone can be expected.

Formalization scope

Space. $E_n$ is EuclideanSpace ℝ (Fin n). $A$ is a Set with four explicit hypotheses: A.Nonempty, IsClosed A, Convex ℝ A, Bornology.IsBounded A. Nonemptiness is implicit in the paper ("of dimension $r$ ", "the point of $A$ nearest to $p$ ").
The image as a relation. IsImage A p p₁ holds when $p_1 = p + 2(q - p)$ for a nearest point $q$ of $A$ ; the nearest point is not chosen by a function. For closed convex nonempty $A$ this relation is a function on $E_n$ . A run is any sequence with IsImage A (p ν) (p (ν+1)) whenever p ν ∉ A; values after entering $A$ are unconstrained. Every statement quantifies over every run from every $p_0 \notin A$ ; a formalization that only asserts the existence of some terminating run is ruled out.
Dimension. $L_r$ is affineSpan ℝ A; $r = n$ is affineSpan ℝ A = ⊤, $r < n$ is affineSpan ℝ A ≠ ⊤. No separate natural number $r$ is introduced.
Symmetry. IsSymmetricWrt L u v: midpoint ℝ u v ∈ L and u - v orthogonal to L.direction. For $r < n - 1$ this is a point reflection through the foot of the perpendicular, not a reflection in a hyperplane. The goal also requires $u \ne v$ and strict alternation $p_\nu \leftrightarrow p_{\nu+1}$ for $\nu > \nu_0$ (strict inequality, as printed).
Boundary is frontier A; distances to sets are Metric.infDist.
Generalizations recorded. Lemma 1, Case 1 is stated for an arbitrary set $A \subseteq E_n$ with full affine span (the paper states it for the polytope of (1.4) and applies it in §10 to a convex body). The projection milestone is stated for any flat $L \supseteq A$ , not only after the paper's reduction to $r = n - 1$ . The §9 milestone expresses "the reflexion process with respect to $F$ amounts to (3.1)" as: every farthest member of $F$ , with any nearest point on it, gives the same step.
Duplication. Case 1 appears both as a milestone and inside the goal, because the paper's proof of Case 2 cites Case 1. Solvers will need to transfer Case 1 from $E_n$ to the flat $L_r$ (an isometric copy of $E_r$ ).
Infrastructure. Mathlib's metric projection onto complete convex sets (exists_norm_eq_iInf_of_complete_convex, norm_eq_iInf_iff_real_inner_le_zero) and EuclideanGeometry.orthogonalProjection cover the basic geometry. Normal cones of convex sets and their upper semicontinuity are not in Mathlib; a contribution there is reusable well beyond this mission. Proofs of any milestone, and alternative proofs of Case 1, are welcome.

Selected references

T. S. Motzkin and I. J. Schoenberg, The relaxation method for linear inequalities, Canadian Journal of Mathematics 6 (1954), 393–404. https://doi.org/10.4153/CJM-1954-038-x
S. Agmon, The relaxation method for linear inequalities, Canadian Journal of Mathematics 6 (1954), 382–392 (the companion paper in the same issue).
L. Fejér, Über die Lage der Nullstellen von Polynomen, die aus Minimumforderungen gewisser Art entspringen, Mathematische Annalen 85 (1922), 41–48.

11 thms2 active usersReviewed

Completed

Convex Optimization Operations Research Optimization·Captain: mikedeng1

Globally Convergent Type-I Anderson Acceleration for Nonsmooth Fixed-Point Iterations: The Stabilized Type-I Anderson Acceleration Converges to a Fixed Point of Every Nonexpansive MapResearch Paper

Motivation

Many first-order methods in optimization are fixed-point iterations $x^{k+1}=f(x^k)$ of a nonexpansive map $f:\mathbb R^n\to\mathbb R^n$ : proximal gradient descent, projected gradient descent, alternating projections, ISTA and Douglas–Rachford splitting (with the conic solver SCS as an instance) all have this form (Zhang, O'Donoghue, Boyd 2020, §4.2). The averaged (Krasnosel'skiĭ–Mann) iteration $x^{k+1}=(1-\alpha)x^k+\alpha f(x^k)$ converges to a fixed point whenever one exists, but it is often slow in its terminal phase. Anderson acceleration (AA) combines the last few iterates through a quasi-Newton update of an approximate Jacobian; it is used for electronic-structure computations and, in the stabilized form studied here, in the solver SCS 2.1. Its type-I variant (AA-I) is often faster in practice than the more studied type-II variant, but it can be numerically unstable, and no global convergence guarantee existed for it on nonsmooth problems.

Timeline (as recounted in the paper's related-work section):

1965. Anderson introduces the method for nonlinear integral equations (J. ACM 12).
1978. Gay and Schnabel prove local Q-superlinear convergence of a full-memory AA-I-type method (Broyden with projected updates), assuming $f$ continuously differentiable near the solution.
2009. Fang and Saad connect AA with multisecant Broyden methods and distinguish the two types (Numer. Linear Algebra Appl. 16).
2011. Walker and Ni show the essential equivalence of full-memory AA with GMRES for affine $f$ (SIAM J. Numer. Anal. 49); Rohwedder and Schneider prove local Q-linear convergence of limited-memory AA-II for differentiable $f$ .
2015. Toth and Kelley give local linear convergence of AA-II for contractive $f$ (SIAM J. Numer. Anal. 53).
2020. Zhang, O'Donoghue and Boyd introduce a stabilized AA-I method (Powell-type regularization, restart checking, safeguarding) and prove global convergence for every nonexpansive map with a fixed point, without differentiability (SIAM J. Optim. 30(4), 3170–3197; longer version arXiv:1808.03971).

Setting

Let $f:\mathbb R^n\to\mathbb R^n$ satisfy $\|f(x)-f(y)\|_2\le\|x-y\|_2$ for all $x,y$ (the Euclidean norm), and assume the solution set $X=\{x^\star\mid x^\star=f(x^\star)\}$ is nonempty. The residual is $g(x)=x-f(x)$ , $g_k=g(x^k)$ , and the averaged operator is $f_\alpha(x)=(1-\alpha)x+\alpha f(x)$ .

Powell's weight. For $\bar\theta\in(0,1)$ , $\phi_{\bar\theta}(\eta)=1$ if $|\eta|\ge\bar\theta$ and $\phi_{\bar\theta}(\eta)=(1-\operatorname{sign}(\eta)\bar\theta)/(1-\eta)$ otherwise, with $\operatorname{sign}(0)=1$ .

Window matrices. For vectors $s_0,\dots,s_{m_k-1}$ and $y_0,\dots,y_{m_k-1}$ , let $\hat s_i$ be their unnormalized Gram–Schmidt orthogonalization (3.2), $B^0=I$ , and

B^{i+1}=B^i+\frac{(\tilde y_i-B^is_i)\hat s_i^T}{\hat s_i^Ts_i},\qquad \tilde y_i=\theta^iy_i+(1-\theta^i)B^is_i,\qquad \theta^i=\phi_{\bar\theta}\Bigl(\frac{\hat s_i^T(B^i)^{-1}y_i}{\|\hat s_i\|_2^2}\Bigr).

The matrix norm $\|\cdot\|_2$ is the induced $\ell_2$ operator norm.

Algorithm 3.1 (AA-I-S-m). With parameters $\bar\theta,\tau,\alpha\in(0,1)$ , $D,\epsilon>0$ and max-memory $m\ge1$ : start from $H_0=I$ , $m_0=0$ , $n_{AA}=0$ , $\bar U=\|g_0\|_2$ and $x^1=\tilde x^1=f_\alpha(x^0)$ . At iteration $k\ge1$ , set $s_{k-1}=\tilde x^k-x^{k-1}$ and $y_{k-1}=g(\tilde x^k)-g(x^{k-1})$ , orthogonalize $s_{k-1}$ against the current window, and restart the window (memory back to $1$ , $H_{k-1}$ replaced by $I$ ) if the memory would exceed $m$ or $\|\hat s_{k-1}\|_2<\tau\|s_{k-1}\|_2$ . Then apply one Powell-regularized rank-one update to obtain $H_k$ and the trial point $\tilde x^{k+1}=x^k-H_kg_k$ . The trial point is accepted if $\|g_k\|_2\le D\bar U(n_{AA}+1)^{-(1+\epsilon)}$ (and $n_{AA}$ increases); otherwise $x^{k+1}=f_\alpha(x^k)$ .

Formalization targets

Goal: Theorem 4.1

\text{for every run of Algorithm 3.1:}\qquad \lim_{k\to\infty}x^k=x^\star\quad\text{for some } x^\star=f(x^\star).

The only hypotheses are nonexpansiveness of $f$ , $X\ne\emptyset$ and the parameter ranges. The limit is not specified: it depends on $x^0$ and the parameters.

Milestones

Lemma 3.2. Well-defined updates give $|\det B^{m_k}|\ge\bar\theta^{m_k}>0$ .
Lemma 3.3. If $\|y_i\|_2\le2\|s_i\|_2$ , $\|\hat s_i\|_2\ge\tau\|s_i\|_2$ and $m_k\le m$ , then $\|B^{m_k}\|_2\le3((1+\bar\theta+\tau)/\tau)^m-2$ .
Corollary 3.4. $\|H_k\|_2\le(3((1+\bar\theta+\tau)/\tau)^m-2)^{n-1}/\bar\theta^m$ (3.8).
Corollary 3.5. Along the algorithm, unless a solution is hit, (3.8) holds and $\mathrm{cond}(H_k)\le(3((1+\bar\theta+\tau)/\tau)^m-2)^n/\bar\theta^m$ .
Eq. (4.3). $\|x^k-y\|_2\le\|x^0-y\|_2+CD\bar U\sum_{i\ge0}(i+1)^{-(1+\epsilon)}$ for every $y\in X$ .
Eq. (4.6). $\lim_k\|g_k\|_2=0$ .
Eq. (4.7). $\|x^{k+1}-y\|_2^2\le\|x^k-y\|_2^2+\epsilon_k$ with $\epsilon_k\ge0$ summable.
§4.1, Step 2. $\|x^k-y\|_2$ converges for every $y\in X$ .

Significance

The theorem places a quasi-Newton acceleration scheme under the same hypotheses as the plain averaged iteration: nonexpansiveness and existence of a fixed point. It therefore applies at once to the nonexpansive examples of §4.2 of the paper (proximal gradient, projected gradient, alternating projections, ISTA, Douglas–Rachford splitting and SCS), with no smoothness, strong convexity or local assumption. The matrix bounds of Lemma 3.3 and Corollaries 3.4–3.5 are also of independent use: they give explicit, iteration-independent control of the approximate inverse Jacobians of a limited-memory type-I method, an assumption that other globalization frameworks (for example SuperMann) have to impose.

The result is proved in the paper; to our knowledge none of it is machine-checked. Mathlib has neither the Krasnosel'skiĭ–Mann iteration, nor Fejér monotonicity, nor any Anderson-type method. A formalization provides these pieces, checks the index bookkeeping of the restart and safeguard steps, and makes explicit the one place where the printed algorithm and the analysis disagree (line 9 at a window start; see below).

Difficulty

The accelerated step $\tilde x^{k+1}=x^k-H_kg_k$ need not decrease the distance to the solution set, so the Fejér argument for averaged iterations does not apply to it directly. The naive fix, bounding $\|H_kg_k\|_2$ , requires a bound on $\|H_k\|_2$ that holds uniformly along the run; without the Powell regularization $B_k$ can be singular, and without the restart rule $\|B_k\|_2$ can grow without bound as the window becomes nearly linearly dependent. The determinant and norm bounds of Lemmas 3.2–3.3 have to be established for arbitrary windows and then connected to the run of the algorithm, where the window, its orthogonalization and the matrices are defined by an intertwined recursion with resets. The safeguard then turns the uniform bound into a summable perturbation of a Fejér-monotone sequence, and the final step needs an Opial-type argument to pass from convergence of distances to convergence of the iterates.

Formalization scope

$\mathbb R^n$ is EuclideanSpace ℝ (Fin n), so all vector norms are Euclidean; matrices are continuous linear maps with the operator norm; $u\hat s^T$ is rankOne ℝ u ŝ; $\det$ is LinearMap.det; inverses are Ring.inverse (zero on singular maps, excluded by Lemma 3.2). The orthogonalization (3.2) is InnerProductSpace.gramSchmidt. Nonexpansive is LipschitzWith 1 f. The algorithm is the predicate IsAAISRun, a deterministic recursion: every run is determined by $f$ , the parameters and $x^0$ , and a sorry-free check that runs exist (for $f=\mathrm{id}$ ) was built locally. The reset of $H_{k-1}$ in line 8 is local to its iteration. $\mathrm{sign}(0)=1$ is encoded explicitly. Constants are the paper's explicit expressions; no milestone replaces them with an existential constant.

Line 9. The paper prints $\tilde y_{k-1}=\theta_{k-1}y_{k-1}-(1-\theta_{k-1})g_{k-1}$ , derived from (3.3) through $B_{k-1}s_{k-1}=-g_{k-1}$ , which fails when the window starts afresh ( $m_k=1$ ). The mission uses (3.3) there: $\tilde y_{k-1}=\theta_{k-1}y_{k-1}+(1-\theta_{k-1})s_{k-1}$ when $m_k=1$ , and the printed formula when $m_k\ge2$ .

Milestones of §4.1 and Corollary 3.5 carry the hypothesis $f(x^k)\ne x^k$ for all $k$ , as the paper does ("we temporarily assume for simplicity that a solution to (1.1) is not found in finite steps"); the goal does not. A goal proved from an unsatisfiable run predicate, or one that assumes a bound on $\|H_k\|_2$ , contractivity of $f$ , or that the accelerated step is never taken, is not this theorem. Division by zero in Lean can only occur once a fixed point has been reached, after which all later iterates coincide.

Useful reusable infrastructure includes the Krasnosel'skiĭ–Mann inequality $\|f_\alpha(x)-y\|_2^2\le\|x-y\|_2^2-\alpha(1-\alpha)\|g(x)\|_2^2$ , quasi-Fejér monotone sequences and their convergence, and determinant and norm identities for rank-one updates (the matrix determinant lemma and Sherman–Morrison). Proofs of the window lemmas, of the run invariants (the window matrices coincide with the algorithm's $H_k^{-1}$ ) and of the convergence steps are all welcome.

Selected references

J. Zhang, B. O'Donoghue, S. Boyd, Globally Convergent Type-I Anderson Acceleration for Nonsmooth Fixed-Point Iterations, SIAM J. Optim. 30(4), 3170–3197, 2020. https://doi.org/10.1137/18M1232772
J. Zhang, B. O'Donoghue, S. Boyd, longer version, 2018. https://arxiv.org/abs/1808.03971
D. M. Gay, R. B. Schnabel, Solving systems of nonlinear equations by Broyden's method with projected updates, in Nonlinear Programming 3, Academic Press, 245–281, 1978 (reference [20] of the paper). https://doi.org/10.1137/18M1232772
D. G. Anderson, Iterative procedures for nonlinear integral equations, J. ACM 12(4), 547–560, 1965. https://doi.org/10.1145/321296.321305
H. Fang, Y. Saad, Two classes of multisecant methods for nonlinear acceleration, Numer. Linear Algebra Appl. 16(3), 197–221, 2009. https://doi.org/10.1002/nla.617
H. F. Walker, P. Ni, Anderson acceleration for fixed-point iterations, SIAM J. Numer. Anal. 49(4), 1715–1735, 2011. https://doi.org/10.1137/10078356X
A. Toth, C. T. Kelley, Convergence analysis for Anderson acceleration, SIAM J. Numer. Anal. 53(2), 805–819, 2015. https://doi.org/10.1137/130919398
H. H. Bauschke, P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer, 2011. https://doi.org/10.1007/978-1-4419-9467-7

12 thms2 active usersReviewed

Numerical Analysis

Missions

Proximal Alternating Linearized Minimization for Nonconvex and Nonsmooth Problems II: The Proximal Map of the Nonnegative Sparsity Constraint Is Hard Thresholding of the Positive PartResearch Paper

Motivation

Setting

Formalization targets

Goal: Proposition 4.1 (Proximal map formula), p. 28

Milestones (in the order the proof uses them)

Significance

Difficulty

Formalization scope

Selected references

Sparsity Constrained Nonlinear Optimization: Optimality Conditions and Algorithms II: Iterative Hard Thresholding for Sparse Least Squares Converges to an L-Stationary PointResearch Paper

Motivation

Setting

Formalization targets

Goal: Theorem 3.2

Milestones, in the order the proof uses them

Significance

Difficulty

Formalization scope

Selected references

Proximal Alternating Linearized Minimization for Nonconvex and Nonsmooth Problems I: Every Bounded PALM Sequence of a KL Function Has Finite Length and Converges to a Critical PointResearch Paper

Motivation

Setting

Formalization targets

Goal: Theorem 3.1 (p. 21)

Milestones, in proof order

Significance

Difficulty

Formalization scope

Selected references

Sparsity Constrained Nonlinear Optimization: Optimality Conditions and Algorithms III: Accumulation Points of the Greedy Sparse-Simplex Method Are Coordinate-Wise MinimaResearch Paper

Motivation

Setting

Formalization targets

Goal: Theorem 3.3

Milestones

Significance

Difficulty

Formalization scope

Selected references

Efficiency of Coordinate Descent Methods on Huge-Scale Optimization Problems 4: Accelerated Coordinate Descent ACDM Has Expected Error at Most (n/(k + 1))²·(2‖x₀ − x*‖₁² + (f(x₀) − f*)/n²)Research Paper

Motivation

Setting

Formalization targets

Goal: Theorem 6, (5.3)

Milestones

Significance

Difficulty

Formalization scope

Selected references

Sparsity Constrained Nonlinear Optimization: Optimality Conditions and Algorithms I: Every Coordinate-Wise Minimum Is an L2(f)-Stationary PointResearch Paper

Motivation

Setting

Formalization targets

Goal: Theorem 2.4

Milestones

Significance

Difficulty

Formalization scope

Selected references

Efficiency of Coordinate Descent Methods on Huge-Scale Optimization Problems 1: Random Block Coordinate Descent RCDM(α, x₀) Has Expected Error at Most 2·S_α·R²_{1−α}(x₀)/(k + 4)Research Paper

Motivation

Setting

Formalization targets

Goal: Theorem 1

Milestones

Significance

Difficulty

Formalization scope

Selected references

Strong Convergence of an Explicit Numerical Method for SDEs with Nonglobally Lipschitz Continuous Coefficients: The Tamed Euler Scheme Converges to the SDE Solution in Uniform Lp with Order 1/2Research Paper

Motivation

Setting

Formalization targets

Goal: Theorem 1.1

Milestones: Lemmas 3.1–3.10

Significance

Difficulty

Efficiency of Coordinate Descent Methods on Huge-Scale Optimization Problems 4: Accelerated Coordinate Descent ACDM Has Expected Error at Most (n/(k + 1))²·(2‖x₀ − x‖₁² + (f(x₀) − f)/n²)Research Paper

Convex Optimization: Algorithms and Complexity XIII: Newton's Method Converges Quadratically, ‖x_{k+1} − x‖ ≤ (M/μ)‖x_k − x‖², from ‖x₀ − x*‖ ≤ μ/(2M)Textbook