Optimization

Motivation

Accelerated gradient methods use information from earlier iterates to reduce the objective gap of a smooth convex optimization problem faster than ordinary gradient descent. Their behavior depends on the precise momentum rule. Attouch, Chbani, Fadili, and Riahi study a discrete rule with an additional difference of consecutive gradients, motivated by an inertial differential equation with Hessian driven damping. This mission concerns their inertial gradient algorithm with Hessian damping, called IGAHD, and the rates proved for it in Theorem 6 of the source paper. The result matters when one wants a rate for objective values together with control of the gradients seen during the run.

The paper appeared as a 2020 preprint and subsequently in Mathematical Programming. The mission follows the pinned arXiv version, whose theorem and equation numbers supply all citations here. The proximal and strongly convex algorithms elsewhere in the paper have separate missions; their dynamics and rate statements differ from IGAHD.

Setting

Let $H$ be a real Hilbert space and let $f:H\to\mathbb R$ be convex and differentiable. Assume that its gradient $\nabla f$ is $L$-Lipschitz: $\|\nabla f(u)-\nabla f(v)\|\le L\|u-v\|$ for every $u,v\in H$, where $L>0$. Choose a minimizer $x^*$, so $f(x^*)=\min_H f$. The positive step size $s$ satisfies $s\le1/L$. The momentum parameter $\alpha$ is at least $3$, and the damping parameter $\beta$ lies in $[0,2\sqrt{s})$.

Given $x_0,x_1\in H$, the algorithm takes the following two steps for every integer $k\ge1$:

$$\begin{aligned} y_k={}&x_k+\left(1-\frac{\alpha}{k}\right)(x_k-x_{k-1}) -\beta\sqrt{s}\bigl(\nabla f(x_k)-\nabla f(x_{k-1})\bigr) -\frac{\beta\sqrt{s}}{k}\nabla f(x_{k-1}),\\ x_{k+1}={}&y_k-s\nabla f(y_k). \end{aligned}$$

The paper sets $t_{k+1}=k/(\alpha-1)$ and $v_k=x_{k-1}-x^*+t_k(x_k-x_{k-1}+\beta\sqrt{s}\nabla f(x_{k-1}))$. Its Lyapunov quantity is $E_k=t_k^2(f(x_k)-f(x^*))+(2s)^{-1}\|v_k\|^2$ for $k\ge1$. These expressions are fixed by the source formulas (14)–(15), p. 15, not chosen to simplify a later proof.

Formalization targets

The central target is the full rate statement of Theorem 6. In the range where the paper's monotonicity argument applies, $E_k$ is eventually nonincreasing. For every IGAHD run,

$$f(x_k)-\min_H f=O(k^{-2})\quad(k\to\infty).$$

When $\beta>0$, the theorem also gives both weighted gradient summability claims:

$$\sum_{k\ge0} k^2\|\nabla f(y_k)\|^2<\infty, \qquad \sum_{k\ge0} k^2\|\nabla f(x_k)\|^2<\infty.$$

The milestone list follows claims the authors state while proving that target. It begins with the extended descent inequality (Lemma 1, (28)), records the identities for $t_k$ and $v_k$, and then captures the displayed energy and $\varepsilon$ inequalities. Remark 4 asks for the $O(k^{-3})$ bound on the smallest squared gradient norm among the first $k$ iterates; Remark 5 asks for the $O(k^{-2})$ objective gap at $y_k$. Those are companion statements and are not substituted for the main theorem.

Significance

The rate bounds give three kinds of information about one trajectory. The objective gap controls the values attained at the iterates. The two summable series imply that the gradients at both the extrapolated points and the iterates become small in a weighted aggregate sense when damping is positive. The paper then uses this series estimate to derive the best-iterate estimate of Remark 4 and the extrapolated-value estimate of Remark 5, all on p. 18.

The result is proved in the cited paper; the open work here is its Lean proof and the supporting smooth convex analysis. The development makes the precise discrete algorithm, the Lyapunov quantity, and the descent inequality reusable as mathematical interfaces. A complete formal proof would also make explicit where the author's finite-index monotonicity wording needs correction.

Difficulty

The usual smooth descent estimate alone gives a decrease after one gradient step, but the extrapolated point $y_k$ also carries a momentum term and a gradient difference. Its objective value is therefore not the same quantity as $f(x_k)$, and a direct comparison of successive objective values does not give the claimed summability. The source uses a stronger descent inequality with a squared difference of gradients. The energy change then contains a quadratic expression in $\nabla f(y_k)$ and $\nabla f(x_k)$ whose sign depends on $\beta<2\sqrt{s}$ and, at first, on how large $k$ is. Establishing that sign uniformly over the relevant tail of the run is the central obstruction.

Formalization scope

The Lean space is a complete real inner product space. Convexity is ConvexOn ℝ Set.univ f; differentiability is explicit because the phrase “whose gradient is $L$-Lipschitz” presupposes a gradient. The existence of a minimizer is represented by a selected $x^*$ with $f(x^*)\le f(z)$ for every $z$. Positive $L$ and positive $s$ make $s\le1/L$ meaningful. The algorithm is a predicate on sequences $x,y:\mathbb N\to H$, imposed from $k=1$ so no coefficient is divided by zero. The first two $x$ terms are free, and $y_0$ is unused. The $O$ statements use the filter at infinity, so finite initial terms do not affect them; summability uses the ordinary series over natural indices.

The printed claim that $(E_k)$ is nonincreasing from its first term is false for the paper's own definitions. On $H=\mathbb R$ with $f(u)=u^2/2$, $L=s=1$, $\alpha=3$, $\beta=0$, $x_0=0$, and $x_1=1$, one gets $E_1=0$ and $E_2=1/8$. The formal goal keeps eventual monotonicity, as supported by the proof on pp. 16–18. This correction is disclosed rather than treated as a new hypothesis. No energy decrease, gradient bound, or conclusion is assumed in the run predicate. The milestone definitions also specify $B_k$ exactly as on p. 17. Contributions to the goal, to the extended descent lemma, or to the explicit energy estimates are all useful; the definition layer and smooth convex lemmas can serve later optimization formalizations.

Selected references

• H. Attouch, Z. Chbani, J. Fadili, and H. Riahi, First-order optimization algorithms via inertial systems with Hessian driven damping, Mathematical Programming, 2020; arXiv:1907.10536v2, especially §3.2, pp. 15–18, and Appendix A.1, p. 34.

Motivation

Accelerated first-order methods for convex minimization, such as Nesterov's method and Polyak's heavy ball method, are now routinely studied through their continuous-time limits: second-order differential equations whose trajectories descend toward a minimizer. In the time-continuous setting the convergence analysis is a Lyapunov argument, and a good Lyapunov function for the ODE usually suggests one for the algorithm obtained by discretizing it.

Attouch, Chbani, Fadili and Riahi (arXiv:1907.10536v2, Math. Program. 2020) study inertial dynamics with Hessian-driven damping: besides the viscous friction term $\gamma\dot x$, the equation contains a term $\beta\nabla^2 f(x(t))\dot x(t)$, which equals $\beta\frac{d}{dt}\nabla f(x(t))$. This damping is known to reduce the oscillations of the heavy ball trajectory, and because it is a time derivative of the gradient, explicit and implicit time discretizations of it produce first-order algorithms with no Hessian evaluation. This mission covers the strongly convex continuous-time result of the paper, Theorem 7 (§4.1, pp. 19–21). For $\beta=0$ it reduces to the heavy ball equation with friction $2\sqrt\mu$, whose exponential rate is due to Siegel (arXiv:1903.05671, Theorem 2.2), as the paper's Remark 7 records. For $\beta>0$, Remark 7 points to a result on a related but different dynamical system by Wilson, Recht and Jordan (arXiv:1611.02635, Theorem 1), whose rate is slightly worse. The algorithms (IPAHD-SC) and (IGAHD-SC) of §5 of the paper are the discrete counterparts of this theorem.

Setting

Let $\mathcal H$ be a real Hilbert space with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$. Let $f:\mathcal H\to\mathbb R$ be a convex function of class $\mathcal C^2$, with gradient $\nabla f$ and Hessian $\nabla^2 f$. For $\mu>0$, $f$ is $\mu$-strongly convex (Definition 1) if $f-\frac\mu2\|\cdot\|^2$ is convex. Such an $f$ is assumed to have a minimizer $x^\star$, and $\min_{\mathcal H} f=f(x^\star)$.

Fix $t_0>0$ and a real constant $\beta$. A solution trajectory of

$$\ddot x(t)+2\sqrt{\mu}\,\dot x(t)+\beta\nabla^2 f(x(t))\dot x(t)+\nabla f(x(t))=0 \tag{19}$$

is a twice differentiable map $x:[t_0,+\infty[\to\mathcal H$ satisfying the equation at every $t\ge t_0$. This is the dynamic $(\mathrm{DIN})_{\gamma,\beta}$ of the paper with the constant viscous damping $\gamma=2\sqrt\mu$.

The proof uses two auxiliary functions. The first is the Lyapunov function

$$\mathcal E(t)=f(x(t))-\min_{\mathcal H}f+\tfrac12\big\|\sqrt\mu(x(t)-x^\star)+\dot x(t)+\beta\nabla f(x(t))\big\|^2 .$$

The second is $Z(t)=2\beta(f(x(t))-f(x^\star))+\sqrt\mu\|x(t)-x^\star\|^2$. Both appear in the milestones under the names lyap and Zfun.

Formalization targets

Goal: Theorem 7 (p. 19)

Assume $0\le\beta\le\frac1{2\sqrt\mu}$. Then three statements hold. First, for all $t\ge t_0$,

$$\frac\mu2\|x(t)-x^\star\|^2\le f(x(t))-\min_{\mathcal H}f\le C e^{-\frac{\sqrt\mu}{2}(t-t_0)},\qquad C:=f(x(t_0))-\min_{\mathcal H}f+\mu\|x(t_0)-x^\star\|^2+\|\dot x(t_0)+\beta\nabla f(x(t_0))\|^2 .$$

Second, there is $C_1>0$ such that for all $t\ge t_0$

$$e^{-\sqrt\mu t}\int_{t_0}^t e^{\sqrt\mu s}\|\nabla f(x(s))\|^2\,ds\le C_1e^{-\frac{\sqrt\mu}{2}t}.$$

Third,

$$\int_{t_0}^\infty e^{\frac{\sqrt\mu}{2}t}\|\dot x(t)\|^2\,dt<+\infty .$$

Milestones

The milestones are the displayed steps of the paper's proof, in order:

1. the strong convexity inequality $\langle\nabla f(y),y-x^\star\rangle\ge f(y)-f(x^\star)+\frac\mu2\|y-x^\star\|^2$;
2. the first differential inequality for $\mathcal E$;
3. the nonnegativity of $\frac\mu4X^2+\frac{\beta}{2\sqrt\mu}Y^2-\beta\sqrt\mu XY$;
4. $\frac{d}{dt}\mathcal E+\frac{\sqrt\mu}2\mathcal E+\frac{\sqrt\mu}2\|\dot x\|^2\le0$;
5. $\mathcal E(t)\le\mathcal E(t_0)e^{-\frac{\sqrt\mu}2(t-t_0)}$;
6. the differential inequality for $Z$.

Companion: the case $\beta=0$

The last sentence of Theorem 7 states that for $\beta=0$, $f(x(t))-\min f=\mathcal O(e^{-\sqrt\mu t})$. The paper gives no proof of it and attributes it to Siegel. It is a separate item, listed as the last milestone, and is not part of the goal.

Significance

Theorem 7 gives exponential decay, at the rate $e^{-\frac{\sqrt\mu}2 t}$, of the values and of the distance to the minimizer for the whole range $0\le\beta\le1/(2\sqrt\mu)$ of Hessian damping. It also gives an averaged exponential decay of $\|\nabla f(x(t))\|^2$, which the authors describe as new in the literature, and weighted integrability of the squared velocity. These estimates guide the design of the strongly convex algorithms of §5, whose linear rates $q=1/(1+\frac12\sqrt{\mu s})$ mirror the continuous rate.

Formalizing it builds machine-checked continuous-time Lyapunov arguments in a Hilbert space: differentiation of energy functions along a trajectory that is only one-sidedly differentiable at $t_0$, the chain rule through $\nabla f$ (the Hessian term), and Grönwall-type integration of differential inequalities. To our knowledge none of Theorem 7 has a machine-checked proof; the paper's proof is complete apart from the gap at $\beta=0$ in part (ii) noted under Difficulty.

Difficulty

The algebra of the first differential inequality is a direct but long expansion, and the quadratic form is where the upper bound on $\beta$ enters. The main analytic steps are elsewhere. First, the energy must be shown differentiable along the trajectory, which requires the derivative of $t\mapsto\nabla f(x(t))$, a chain rule through the Fréchet derivative of the gradient. Second, the differential inequalities must be integrated on $[t_0,+\infty[$ with one-sided derivatives at $t_0$. Third, part (ii) must be derived from the $Z$-inequality. The paper's derivation divides by $\beta^2$, so it covers only $\beta>0$. At $\beta=0$ part (ii) is still true, but it needs a separate argument, for example local Lipschitz continuity of $\nabla f$ near $x^\star$ together with the convergence $x(t)\to x^\star$ from part (i). A proof that only follows the printed text does not cover the case $\beta = 0$.

Formalization scope

• Space and function. $\mathcal H$ is a real inner product space that is complete. $f$ is ConvexOn ℝ univ and ContDiff ℝ 2, $\nabla f$ is Mathlib's gradient f, and the Hessian applied to $v$ at $x$ is fderiv ℝ (gradient f) x v.
• Strong convexity is Definition 1 verbatim, with $\mu>0$ as a separate hypothesis. Convexity of $f$ is kept from the standing hypothesis (H) of p. 2 even though strong convexity implies it.
• The minimizer. $\operatorname{argmin}f\neq\emptyset$ is a point xstar with $f(x^\star)\le f(y)$ for all $y$. $\min f$ is written f xstar.
• The trajectory is a triple $x,\dot x,\ddot x:\mathbb R\to\mathcal H$. At every $t\ge t_0$, $x$ has derivative $\dot x(t)$ and $\dot x$ has derivative $\ddot x(t)$ within $[t_0,+\infty[$, and (19) holds. Values before $t_0$ play no role, and $\dot x(t_0)$ in the constant $C$ is the right derivative. $\beta$ is a real constant and $t_0>0$, as in (H).
• Derivatives in milestones are stated as the existence of a number that is the one-sided derivative and satisfies the inequality. Such a number is unique.
• Integrals. Part (ii) is an interval integral of a function continuous on $[t_0,t]$. The "Moreover" clause is stated as integrability on $[t_0,+\infty[$ of the nonnegative integrand. The constant $C_1$ is existentially quantified after the trajectory, so it may depend on it, as on the page.
• Ruling out trivial readings. The goal never mentions $\mathcal E$ or $Z$: its conclusions are about $f(x(t))$, $\|x(t)-x^\star\|$, $\nabla f(x(t))$ and $\dot x(t)$ only. The solution predicate is satisfiable: the constant trajectory at $x^\star$ satisfies it.
• Infrastructure needed: the chain rule for $\nabla f\circ x$ with one-sided derivatives, differentiation of $\|\cdot\|^2$ and $\langle\cdot,\cdot\rangle$ along curves, and a Grönwall-type lemma $u'\le-cu\Rightarrow u(t)\le u(t_0)e^{-c(t-t_0)}$ on $[t_0,\infty)$ with one-sided derivatives. These are reusable for the other continuous-time missions of the series. Contributions of these general lemmas are welcome as separate theorems.

Selected references

• H. Attouch, Z. Chbani, J. Fadili, H. Riahi, First-order optimization algorithms via inertial systems with Hessian driven damping, Math. Program. (2020); arXiv:1907.10536v2. https://arxiv.org/abs/1907.10536
• J. W. Siegel, Accelerated first-order methods: differential equations and Lyapunov functions, 2019. https://arxiv.org/abs/1903.05671
• A. C. Wilson, B. Recht, M. I. Jordan, A Lyapunov analysis of momentum methods in optimization, J. Mach. Learn. Res. 22 (2021); arXiv:1611.02635. https://arxiv.org/abs/1611.02635
• B. T. Polyak, Some methods of speeding up the convergence of iteration methods, USSR Comput. Math. Math. Phys. 4 (1964). https://doi.org/10.1016/0041-5553(64)90137-5

Motivation

A decision maker who chooses $x$ to minimize an expected cost $c(x,\mathbb P)=\mathbb E_{\mathbb P}[\gamma(x,\xi)]$ rarely knows the distribution $\mathbb P$ of the random parameter $\xi$. What is available is a sample $\xi_1,\dots,\xi_T$. Any rule that turns the sample into a cost estimate can disappoint: the realized expected cost may exceed the estimate. In a cost-minimization context such underestimates are more harmful than overestimates, which is why distributionally robust optimization replaces the unknown $\mathbb P$ by a worst case over a set of distributions consistent with the data.

Van Parys, Mohajerin Esfahani and Kuhn (arXiv:1704.04118, Management Science 2021) asked which such rule is best. They define a data-driven predictor to be least conservative among all predictors whose probability of disappointment decays exponentially at a prescribed rate $r$, and they prove that a single predictor, the worst-case expected cost over a relative-entropy ball around the empirical distribution, is the unique strong solution. Their main results are for a finite set of outcomes $\Xi=\{1,\dots,d\}$. Section 5 extends the predictor result to an arbitrary compact $\Xi\subseteq\mathbb R^d$; this mission formalizes that extension. The version used throughout is the arXiv preprint arXiv:1704.04118v3 (22 December 2019), and every index below is the preprint's.

Setting

Fix compact sets $X\subseteq\mathbb R^n$ (decisions) and $\Xi\subseteq\mathbb R^d$ (outcomes), and a cost $\gamma:X\times\Xi\to\mathbb R$ that is jointly continuous.

• The model class $\mathcal P$ is the set of Borel probability distributions on $\Xi$, with the topology of weak convergence. $X\times\mathcal P$ carries the product topology.
• The model-based predictor is $c(x,\mathbb P)=\int_\Xi\gamma(x,\xi)\,\mathrm d\mathbb P(\xi)$.
• The relative entropy of $\mathbb P'$ with respect to $\mathbb P$ (Definition 8) is $I(\mathbb P',\mathbb P)=\int_\Xi\log\frac{\mathrm d\mathbb P'}{\mathrm d\mathbb P}\,\mathrm d\mathbb P'$ if $\mathbb P'\ll\mathbb P$, and $+\infty$ otherwise.
• Given independent samples $\xi_1,\dots,\xi_T$ from $\mathbb P$, the empirical distribution is $\hat{\mathbb P}_T=\frac1T\sum_{t=1}^T\delta_{\xi_t}$, and $\mathbb P^\infty$ denotes the law of the sample path.
• A data-driven predictor is a continuous function $\hat c:X\times\mathcal P\to\mathbb R$; the estimate after $T$ samples is $\hat c(x,\hat{\mathbb P}_T)$. Its out-of-sample disappointment is $\mathbb P^\infty\big(c(x,\mathbb P)>\hat c(x,\hat{\mathbb P}_T)\big)$.
• Problem (5): $\hat c$ is feasible if, for all $x\in X$ and $\mathbb P\in\mathcal P$,

$$\limsup_{T\to\infty}\frac1T\log\mathbb P^\infty\big(c(x,\mathbb P)>\hat c(x,\hat{\mathbb P}_T)\big)\le -r,$$

and strongly optimal if it is feasible and $\hat c(x,\mathbb P')\le\hat c'(x,\mathbb P')$ for all $(x,\mathbb P')$ and every feasible $\hat c'$.

• The distributionally robust predictor is $\hat c_r(x,\mathbb P')=\sup\{c(x,\mathbb P):\mathbb P\in\mathcal P,\ I(\mathbb P',\mathbb P)\le r\}$.
• $\bar\gamma(x)=\max_{\xi\in\Xi}\gamma(x,\xi)$ is the worst-case cost and $\Xi^\star(x)$ the set of its maximizers.

Formalization targets

Goal: Theorem 10 (p. 26)

$$r\ge0\ \Longrightarrow\ \hat c_r\text{ is feasible in (5)},\qquad r>0\ \Longrightarrow\ \hat c_r\text{ is strongly optimal in (5)}.$$

Both halves are in the goal. The case $r=0$ of feasibility is the statement that $\hat c_0=c$ is continuous.

Milestones

• Lemma 1 (p. 23): $c$ is continuous on $X\times\mathcal P$.
• Lemma 2 (33), p. 29: $\hat c_r$ equals a supremum over pairs $(\mathbb P_c,p)$ with $\mathbb P'\ll p\,\mathbb P_c\ll\mathbb P'$, where the mass $1-p$ is placed at cost $\bar\gamma(x)$.
• Lemma 3 (p. 31): the perturbed problem (34), $\hat c_{r,\epsilon}$, satisfies $\hat c_r\le\hat c_{r,\epsilon}\le\hat c_r+\epsilon$.
• Lemma 4 (35), p. 31:

$$\hat c_{r,\epsilon}(x,\mathbb P')\le\min_{\alpha\ge\bar\gamma(x)+\epsilon}\ \alpha-e^{-r}\exp\Big(\int_\Xi\log(\alpha-\gamma(x,\xi))\,\mathrm d\mathbb P'(\xi)\Big),$$

with a minimizer $\alpha^\star\le(\bar\gamma(x)+\epsilon-e^{-r}c(x,\mathbb P'))/(1-e^{-r})$, and with equality for $\epsilon>0$.

• Proposition 6 (p. 25): $\hat c_r$ is continuous on $X\times\mathcal P$ for $r\ge0$.
• Theorem 9 (24a)/(24b), p. 25: the weak large deviation principle for $\hat{\mathbb P}_T$ in the weak topology, with the closure $\operatorname{cl}\mathcal D$ in the upper bound and the interior in the lower bound.
• (41), p. 35: if $\mathbb P_0(\Xi^\star(x))<1$ and $c(x,\mathbb P_0)\ge\hat c_r(x,\mathbb P')$, then $I(\mathbb P',\mathbb P_0)\ge r$.
• Proposition 1(ii), (iii) in the general setting (p. 24): $I$ is jointly convex and jointly lower semicontinuous on $\mathcal P\times\mathcal P$.

Proposition 5, the dual representation (23) of $\hat c_r$ itself, is included as a further item; the proof of Theorem 10 goes through Lemmas 3–4 and does not use it.

Significance

Theorem 10 says that, on any compact outcome space, no continuous predictor is less conservative than $\hat c_r$ while keeping its disappointment probability below $e^{-rT}$ to first order in the exponent. Relative-entropy distributionally robust optimization is therefore not one choice among many ambiguity sets but the optimal one for this criterion, and the dual representation turns its evaluation into a one-dimensional convex problem. The paper also shows that the continuous-state analogue of its prescriptor result (Theorem 11) only holds up to an arbitrarily small shift, so the predictor statement formalized here is the clean part of the extension.

The results are proved in the paper; none is formalized. A complete development yields a formal Sanov-type weak large deviation principle for empirical measures in the weak topology, joint convexity and lower semicontinuity of the Kullback–Leibler divergence for probability measures on a compact metric space, and a measure-theoretic dual representation of a relative-entropy worst-case expectation. Mathlib has the Kullback–Leibler divergence klDiv and the weak topology on probability measures, but none of these three results.

Difficulty

The finite-state proof of feasibility bounds the disappointment through the large deviations upper bound over the disappointment set $\mathcal D$ itself. On a continuous space that bound only holds over the closure $\operatorname{cl}\mathcal D$ in the weak topology (24a), and the paper notes that this invalidates the finite-state argument (p. 25). The rate must then be controlled on a larger, closed set of empirical distributions, where the strict inequality defining a disappointment is lost; this is where the mass that $\mathbb P_0$ puts on the worst-case scenarios $\Xi^\star(x)$ starts to matter, as the hypothesis of (41) shows.

Continuity of $\hat c_r$ in the weak topology is the second obstacle. The supremum runs over an infinite-dimensional set of distributions, a worst-case distribution may put mass outside the support of $\mathbb P'$, and the integrals involved can diverge or take the value $-\infty$. Weak convergence controls integrals of bounded continuous functions only, and $\log(\alpha-\gamma)$ is unbounded below at $\alpha=\bar\gamma(x)$.

Formalization scope

• $\mathcal P$ is ProbabilityMeasure ↥Ξ for the Borel σ-algebra of the subtype, with Mathlib's weak-convergence topology. $X$ and $\Xi$ are compact subsets of Euclidean spaces, $\gamma:↥X\to↥\Xi\to\mathbb R$ is jointly continuous; these standing assumptions of §2 (p. 5) and §5 (p. 23) appear as hypotheses of every theorem. No other hypothesis is added. Nonemptiness of $\Xi$ is automatic wherever a distribution on $\Xi$ is in scope.
• $I$ is Mathlib's InformationTheory.klDiv, valued in $[0,\infty]$. It equals Definition 8, including the value $+\infty$ when $\mathbb P'\not\ll\mathbb P$ or when the log-likelihood ratio is not integrable.
• The empirical distribution is defined for $T\ge1$. The probability of an event about $\hat{\mathbb P}_T$ is the outer measure under the product measure $\mathbb P^{\otimes T}$, and is $0$ at $T=0$; all statements are asymptotic in $T$.
• Decay rates are stated without logarithms: $\limsup\frac1T\log p_T\le-s$ becomes "for every $r'<s$, eventually $p_T\le e^{-r'T}$", which is equivalent and handles $p_T=0$. A naive Real.log encoding would be wrong, since Lean sets $\log0=0$. Predictors must be continuous for the weak topology; dropping continuity, or replacing the weak topology by a finer one, changes the problem.
• The constraint $\int\log(\frac1p\frac{\mathrm d\mathbb P'}{\mathrm d\mathbb P_c})\,\mathrm d\mathbb P'\le r$ in (33)–(34) carries an explicit integrability requirement, so that a divergent integral fails the constraint instead of evaluating to $0$. The geometric mean $\exp\int\log(\alpha-\gamma)\,\mathrm d\mathbb P'$ follows the paper's convention $\log0=-\infty$; it is defined as $\inf_{\delta>0}\exp\int\log(\alpha+\delta-\gamma)\,\mathrm d\mathbb P'$, which is the same value for $\alpha\ge\bar\gamma(x)$.
• At the end of the proof of Theorem 10 the page writes "when $\epsilon>0$" for strong optimality; the theorem statement says $r>0$, which is what is formalized.

Contributions are welcome at every level: a Sanov upper and lower bound for empirical measures on compact metric spaces, convexity and lower semicontinuity of klDiv, and the measure-theoretic duality of Lemma 4 are each reusable beyond this mission.

Selected references

• B. P. G. Van Parys, P. Mohajerin Esfahani, D. Kuhn, From Data to Decisions: Distributionally Robust Optimization is Optimal, Management Science 67(6), 2021. Preprint arXiv:1704.04118v3. https://arxiv.org/abs/1704.04118
• I. Csiszár, A simple proof of Sanov's theorem, Bulletin of the Brazilian Mathematical Society 37(4), 2006. https://doi.org/10.1007/s00574-006-0023-1
• T. van Erven, P. Harremoës, Rényi divergence and Kullback–Leibler divergence, IEEE Transactions on Information Theory 60(7), 2014. https://doi.org/10.1109/TIT.2014.2320500
• A. Dembo, O. Zeitouni, Large Deviations Techniques and Applications, 2nd ed., Springer, 1998. https://doi.org/10.1007/978-1-4612-5320-4

Motivation

A classifier trained from finitely many examples can perform well on those examples while responding poorly to changes in the feature distribution or to mislabeled observations. Distributionally robust learning addresses this by choosing a classifier against every probability distribution within a specified distance of the empirical distribution. The resulting optimization appears to range over an infinite-dimensional family of distributions. Shafieezadeh-Abadeh, Kuhn, and Mohajerin Esfahani show that, for linear classification with a convex Lipschitz loss, this problem has an exact finite convex reformulation. That result is Theorem 3.11(ii) of the 2019 arXiv preprint, which is the source used here.

The authors' transport metric allows a data point's features to move and its binary label to flip, with separate prices for those two changes. This matters when robustness to inaccurate labels is part of the model. A label-preserving metric cannot express that choice. The paper builds on its Section 2 learning model and the transport duality developed in Appendix A.1; the present mission isolates the classification result and the intermediate statements used for it. The preprint's Corollary 3.14 specializes the result to logloss, while its Corollary 3.12 treats hinge loss. A logloss-specific reformulation and a label-preserving classification result already have proved statements on Prove2Me; the general convex-Lipschitz, finite label-cost theorem remains the target of this mission.

Setting

Let $V$ be a finite-dimensional real vector space with a chosen norm $\|\cdot\|$. Features are $x\in V$ and labels are $y\in\{-1,+1\}$. A linear classifier is a continuous linear functional $w:V\to\mathbb R$, with score $\langle w,x\rangle=w(x)$. Its dual norm $\|w\|_*$ is the operator norm induced by $\|\cdot\|$. Given $N\ge1$ training pairs $(\hat x_i,\hat y_i)$, their empirical distribution is $\hat P_N=N^{-1}\sum_{i=1}^N\delta_{(\hat x_i,\hat y_i)}$.

For $\kappa>0$, the transport cost between pairs is

$$d((x,y),(x',y'))=\|x-x'\|+\kappa\mathbf1_{y\ne y'}.$$

The Wasserstein distance $W(Q,P)$ is the least expected transport cost among joint distributions whose two marginals are $Q$ and $P$. For a radius $\rho\ge0$, the Wasserstein ball $\mathbb B_\rho(\hat P_N)$ contains probability distributions $Q$ satisfying $W(Q,\hat P_N)\le\rho$. The loss is $L(y\langle w,x\rangle)$, where $L:\mathbb R\to[0,\infty)$ is convex and Lipschitz continuous. The Lipschitz modulus $\operatorname{lip}(L)$ is the smallest Lipschitz bound for $L$, defined as the supremum of its difference quotients. These are the objects and standing conventions of Sections 1.1, 2, and 3.2 of the preprint.

Formalization targets

The goal is Theorem 3.11(ii): for every fixed $w$, the worst-case expected loss equals the value of program (18),

$$\sup_{Q\in\mathbb B_\rho(\hat P_N)}\mathbb E^Q[L(y\langle w,x\rangle)] = \inf_{\lambda,s}\left\{\lambda\rho+\frac1N\sum_{i=1}^{N}s_i: \begin{array}{l} L(\hat y_i\langle w,\hat x_i\rangle)\le s_i,\\ L(-\hat y_i\langle w,\hat x_i\rangle)-\kappa\lambda\le s_i\quad(i=1,\ldots,N),\\ \operatorname{lip}(L)\|w\|_*\le\lambda \end{array}\right\}.$$

Taking the infimum of both sides over $w$ gives equality between problem (4) and the full program (18). The fixed-$w$ identity is stated explicitly because it is the stronger assertion established in the proof. The milestone list follows the paper's appendix: Lemma A.1 converts the robust expectation into a scalar transport-price infimum; the label split isolates the two possible labels; Lemma A.3 evaluates the feature supremum; the ensuing display gives the program in terms of conjugate slopes; and the final identity identifies those slopes with $\operatorname{lip}(L)$. Each milestone is tied to its printed page in the preprint.

Significance

The theorem replaces a search over probability distributions by a finite convex program with one scalar $\lambda$ and one slack variable $s_i$ per sample, in addition to the classifier $w$. It makes the price of label changes explicit through $\kappa\lambda$ and the price of feature changes through $\operatorname{lip}(L)\|w\|_*$. Consequently, different convex Lipschitz classification losses can be treated within the same transport model rather than by deriving a new distributional optimization problem for each loss. The preprint applies the result to hinge, smoothed hinge, and logloss examples in the paragraphs following Theorem 3.11.

The result is proved in the cited paper. Its general statement and the paper's appendix milestones are not yet machine-checked in this mission. A verified development would connect the published definitions of the Wasserstein ball and Lipschitz modulus to an exact program-value theorem that later work can import. The loss-specific proved statements already available on Prove2Me provide useful comparisons, but their narrower losses or label-preserving transport costs do not replace this theorem.

Difficulty

The visible finite constraints contain no probability distributions, while the left-hand side optimizes over all probability measures in a Wasserstein ball. Establishing that the two values coincide requires control of the transport budget and the inner worst-case loss at the same time. The straightforward comparison of feasible values gives only one inequality. It does not show that every improvement allowed by a distribution can be represented by the scalar transport price and the samplewise constraints. The case $\rho=0$ also needs care: Lemma A.1's strong-duality proof cites a result for strictly positive radius, whereas the theorem itself allows zero radius and states an infimum rather than an attained minimum.

Formalization scope

Lean represents $V$ by an abstract finite-dimensional real normed space with its Borel measurable structure. This retains the arbitrary norm of $\mathbb R^n$ in the paper. Labels are Bool, interpreted as $+1$ and $-1$; the published DRLogReg_Reformulation_Core definition supplies their sign map, the metric, Wasserstein distance, ball, and empirical distribution. The published WassersteinDRO_Duality_lipschitzModulus definition supplies the exact Lipschitz modulus. The mission's own definitions add the general loss, the feasible set of (18), the fixed-classifier program value, and the effective domain $\Theta=\{\theta:L^*(\theta)<\infty\}$ of the convex conjugate. The same objects and conventions are used across all milestones.

All potentially infinite robust values are represented in $[0,\infty]$; the unbounded feature suprema in Lemma A.3 and the label split use extended reals. The standing loss condition $L\ge0$ comes from the paper's Section 2.1 loss codomain. It also makes the lower Lebesgue integral and the finite program's extended-nonnegative objective faithful. The sample size is positive, $\kappa>0$, and the goal allows $\rho=0$. Lemma A.1 separately assumes $\rho>0$, matching the strong-duality condition used on page 29. The printed extra index $j\in[J]$ in the page-35 classification display is carried in its verbatim milestone provenance but not in the Lean statement, because no $j$ exists in part (ii). The formalization does not obtain the theorem by restricting $L$ to a constant loss, by taking a singleton feature space, or by substituting an arbitrary Lipschitz bound for its modulus. Contributions that establish transport duality, the conjugate-domain slope identity, or the zero-radius case are directly reusable.

Selected references

• S. Shafieezadeh-Abadeh, D. Kuhn, and P. Mohajerin Esfahani, Regularization via Mass Transportation, arXiv preprint arXiv:1710.10016v3, 2019. Preprint.
• S. Shafieezadeh-Abadeh, P. Mohajerin Esfahani, and D. Kuhn, Distributionally Robust Logistic Regression, Advances in Neural Information Processing Systems 28, 2015. Published Prove2Me formalization.

Motivation

Many problems in statistics and signal processing ask for an orthonormal frame $X\in\mathbb R^{n\times r}$ that is both a good fit to the data and sparse or otherwise structured. Sparse principal component analysis, compressed modes in physics, unsupervised feature selection and sparse blind deconvolution all take the form

$$\min_{X}\ F(X)=f(X)+h(X)\quad\text{s.t.}\quad X^\top X=I_r,$$

with $f$ smooth and $h$ convex but nonsmooth, typically an $\ell_1$ penalty. Before the work of Chen, Ma, So and Zhang, algorithms for this class fell into three groups: subgradient-oriented methods, which solve a quadratic program over sampled subgradients in every step; proximal point algorithms, whose subproblems are as hard as the original problem; and operator-splitting methods such as SOC and MADMM. The manifold proximal gradient method (ManPG) carries the Euclidean proximal gradient method over to the Stiefel manifold. The paper shows that ManPG keeps the global convergence and the $O(\varepsilon^{-2})$ iteration complexity that Riemannian gradient descent has for smooth objectives (Boumal, Absil and Cartis, 2019).

Setting

The feasible set is the Stiefel manifold

$$\mathcal M=\mathrm{St}(n,r)=\{X\in\mathbb R^{n\times r}: X^\top X=I_r\},\qquad r\le n.$$

The ambient space $\mathbb R^{n\times r}$ carries the Frobenius inner product $\langle A,B\rangle=\operatorname{Tr}(A^\top B)$ and norm $\|A\|_F$. The tangent space at $X\in\mathcal M$ is $T_X\mathcal M=\{V: V^\top X+X^\top V=0\}$. The orthogonal projection onto it is $\operatorname{Proj}_{T_X\mathcal M}(Y)=(I_n-XX^\top)Y+\tfrac12X(X^\top Y-Y^\top X)$.

The standing assumptions on $F=f+h$ are imposed on all of $\mathbb R^{n\times r}$:

• (i) $f$ is differentiable, and its gradient $\nabla f$ is $L$-Lipschitz;
• (ii) $h$ is convex and $L_h$-Lipschitz.

A point $X\in\mathcal M$ is stationary (Definition 3.3) if

$$0\in\operatorname{grad}f(X)+\operatorname{Proj}_{T_X\mathcal M}\partial h(X),\qquad \operatorname{grad}f(X)=\operatorname{Proj}_{T_X\mathcal M}\nabla f(X),$$

where $\partial h$ is the convex subdifferential.

A retraction $\mathrm{Retr}_X(\xi)$ maps tangent vectors at $X$ back to $\mathcal M$, with $\mathrm{Retr}_X(0)=X$. On the compact manifold $\mathcal M$ every retraction satisfies, for some $M_1,M_2>0$ and all $X\in\mathcal M$, $\xi\in T_X\mathcal M$ (Fact 3.6),

$$\|\mathrm{Retr}_X(\xi)-X\|_F\le M_1\|\xi\|_F,\qquad \|\mathrm{Retr}_X(\xi)-(X+\xi)\|_F\le M_2\|\xi\|_F^2.$$

ManPG (Algorithm 1) takes $X_0\in\mathcal M$, $\gamma\in(0,1)$ and $t>0$. In iteration $k$ it solves the convex subproblem

$$V_k=\operatorname*{argmin}_{V\in T_{X_k}\mathcal M}\ \langle\nabla f(X_k),V\rangle+\frac1{2t}\|V\|_F^2+h(X_k+V).\tag{4.3}$$

It then backtracks $\alpha=1,\gamma,\gamma^2,\dots$ until $F(\mathrm{Retr}_{X_k}(\alpha V_k))\le F(X_k)-\alpha\|V_k\|_F^2/(2t)$, and sets $X_{k+1}=\mathrm{Retr}_{X_k}(\alpha_kV_k)$.

A point $X\in\mathcal M$ is $\varepsilon$-stationary (Definition 5.4) if the solution of (4.3) at $X$ with $t=1/L$ has $\|V\|_F\le\varepsilon/L$. The proof of Lemma 5.2 uses a bound $G$ on $\|\nabla f\|_F$ over $\mathcal M$ and the constants

$$c_0=M_2G+\frac{LM_1^2}{2},\qquad \bar\alpha=\frac{1}{2(c_0+L_hM_2)t}.$$

Formalization targets

Goal: Theorem 5.5

For every run of ManPG:

1. every limit point $\bar X$ of $\{X_k\}$ is stationary;
2. with $t=1/L$, for every $\varepsilon>0$ some

$$k\le\left\lceil\frac{2L\,(F(X_0)-F^*)}{\gamma\bar\alpha\varepsilon^2}\right\rceil$$

has $\|V_k\|_F\le\varepsilon/L$ and $X_k$ is $\varepsilon$-stationary.

Here $F^*$ is the optimal value of (1.1). The constant is the paper's, kept exactly.

Milestones, in the order the proof uses them

• Lemma 5.1: $g(\alpha V_k)-g(0)\le\frac{(\alpha-2)\alpha}{2t}\|V_k\|_F^2$ for $\alpha\in[0,1]$, where $g$ is the objective of (4.3).
• (5.5): $f(\mathrm{Retr}_{X}(\alpha V))-f(X)\le\alpha\langle\nabla f(X),V\rangle+c_0\alpha^2\|V\|_F^2$.
• Lemma 5.2, split into three parts:
  • sufficient decrease for every $0<\alpha\le\min\{1,\bar\alpha\}$;
  • well-definedness of the line search: a run exists from every $X_0$;
  • the run form $\gamma\min\{1,\bar\alpha\}\le\alpha_k\le1$, with $F(X_{k+1})-F(X_k)\le-\frac{\alpha_k}{2t}\|V_k\|_F^2$.
• Lemma 5.3: $V_k=0$ implies $X_k$ is stationary.
• Two displays of the proof of Theorem 5.5:
  • $\|V_k\|_F^2\to0$;
  • the telescoped chain $F(X_0)-F^*>\frac{tK\varepsilon^2}{2}\gamma\bar\alpha$ when $\|V_k\|_F>\varepsilon/L$ for all $k<K$.

Significance

Theorem 5.5 is the paper's main result. It certifies ManPG as a first-order method for nonsmooth problems on the Stiefel manifold with the same $O(\varepsilon^{-2})$ complexity as Riemannian gradient descent in the smooth case. When $h=0$ the bound reduces to the Boumal–Absil–Cartis rate (Remark 5.6). Later proximal methods on matrix manifolds build on this analysis. Lemmas 5.1–5.3 are the reusable steps: descent of the subproblem, sufficient decrease after retraction, and the stopping test as a stationarity certificate.

The result is proved in the paper; it is not formalized anywhere. No proof assistant library contains the convergence theory of a proximal method on a matrix manifold. The mission produces a machine-checked version of the full chain from the subproblem to the iteration bound. Along the way it builds infrastructure that Mathlib lacks: the Frobenius geometry of $\mathbb R^{n\times r}$ as a trace formula, the tangent space and projection of the Stiefel manifold, and optimality conditions for convex problems restricted to a subspace.

Difficulty

The bookkeeping inequalities (Lemma 5.1, (5.5), the decrease estimate, the telescoping) are short computations once the right facts are available. The difficulty lies elsewhere.

• Lemma 5.3 and the limit-point claim need the optimality condition of the subproblem (4.3): a subdifferential sum rule for a convex function restricted to the linear subspace $T_X\mathcal M$, followed by the projection formula. The paper cites this condition rather than proving it.
• The limit-point claim also needs a closedness argument. Subgradients of the Lipschitz function $h$ are bounded by $L_h$, so along a convergent subsequence they have a limit point that is again a subgradient. The projection $X\mapsto\operatorname{Proj}_{T_X\mathcal M}$ must be shown continuous. The paper's proof leaves both steps implicit.
• The well-definedness of the line search requires existence of a minimizer of (4.3), by strong convexity and coercivity on a subspace.
• The descent lemma for $f$ has to be derived from the derivative along lines and the Lipschitz gradient.

A natural first idea is to restate everything for the Riemannian gradient and reuse smooth Riemannian results. That route does not work: $h$ is nonsmooth, and the argument runs in the ambient space through the Euclidean subdifferential.

Formalization scope

• Representation.
  • Matrices are Matrix (Fin n) (Fin r) ℝ.
  • The inner product is trace (Aᵀ * B) and the norm is its square root. No Mathlib matrix norm instance is used.
  • $\mathrm{St}(n,r)$ is the published definition ProjLikeRetr.Stiefel.stiefel.
• Assumptions.
  • $f$ comes with a gradient map gradf, tied to it by derivatives along every line.
  • $\nabla f$ is $L$-Lipschitz in the Frobenius norm, with $L>0$.
  • $h$ is convex and $L_h$-Lipschitz on the whole ambient space.
  • $r\le n$.
• Retraction. The retraction is any map sending tangent vectors at points of $\mathcal M$ into $\mathcal M$, with $\mathrm{Retr}_X(0)=X$ and the bounds (3.1)–(3.2). By Fact 3.6 this covers every retraction; smoothness is never used.
• Constants.
  • $G$, $M_1$, $M_2$, $L$, $L_h$ are parameters, each with its defining hypothesis.
  • $c_0$ and $\bar\alpha$ are the explicit formulas above.
  • $F^*$ is an attained lower bound of $F$ on $\mathcal M$.
• The run. A run of Algorithm 1 is a predicate on sequences $(X_k,V_k,\alpha_k)$. Each $V_k$ minimizes (4.3), $\alpha_k=\gamma^j$ for the least passing $j$, and $X_{k+1}=\mathrm{Retr}_{X_k}(\alpha_kV_k)$.
• The goal's conventions.
  • A limit point is a cluster point of the sequence.
  • "Returns an $\varepsilon$-stationary point within $N$ iterations" means some $k\le N$ passes the stopping test $\|V_k\|_F\le\varepsilon/L$, with $N=\lceil\cdot\rceil$ as Nat.ceil.

Two encodings would trivialize the theorem, and the mission excludes both.

• An unsatisfiable run predicate would make the goal vacuous. The well-definedness milestone asserts that a run exists from every $X_0\in\mathcal M$.
• A projection or stationarity condition that holds trivially would make the first claim empty. Stationarity uses the explicit projection formula and the convex subdifferential on all of $\mathbb R^{n\times r}$.

Lemma 5.2 as printed fixes $\alpha\le\min\{1,\bar\alpha\}$ and asserts the decrease for the accepted iterate. What its proof establishes is three separate facts: decrease at every trial stepsize $\alpha\le\min\{1,\bar\alpha\}$, termination of the line search, and $\alpha_k\ge\gamma\min\{1,\bar\alpha\}$ along a run. The mission states these three facts.

Contributions are welcome at every level: proofs of the milestones, and general lemmas such as the descent lemma for $C^{1,1}$ functions on $\mathbb R^{n\times r}$, the optimality condition of a strongly convex problem on a subspace, and properties of the Stiefel tangent projection. Theorem numbers and pages refer to arXiv:1811.00980v2, the version read for this mission.

Selected references

• S. Chen, S. Ma, A. M.-C. So, T. Zhang, Proximal gradient method for nonsmooth optimization over the Stiefel manifold, SIAM Journal on Optimization 30(1), 2020. https://doi.org/10.1137/18m122457x (arXiv:1811.00980v2: https://arxiv.org/abs/1811.00980)
• N. Boumal, P.-A. Absil, C. Cartis, Global rates of convergence for nonconvex optimization on manifolds, IMA Journal of Numerical Analysis 39(1), 2019. https://doi.org/10.1093/imanum/drx080
• P.-A. Absil, R. Mahony, R. Sepulchre, Optimization Algorithms on Matrix Manifolds, Princeton University Press, 2008. https://doi.org/10.1515/9781400830244

Motivation

Best subset selection asks for the linear model that fits a response best using only a few of many candidate predictors. In its penalized form it is the least-squares problem with an $L_0$ penalty, the number of nonzero coefficients. The problem is NP-hard. Mixed-integer optimization can solve it to certified optimality at moderate sizes (Bertsimas, King, Mazumder 2016), but statisticians who fit whole regularization paths with tens of thousands of predictors need something faster. Hazimeh and Mazumder (arXiv:1803.01454, Operations Research 2020) proposed cyclic coordinate descent for the $L_0$-penalized problem, the method that glmnet made standard for the lasso. Their algorithm underlies the L0Learn software package.

Coordinate descent on a convex smooth-plus-separable objective converges by classical arguments (Tseng 2001). Those arguments do not cover the $L_0$ penalty, which is discontinuous. The paper adds occasional spacer steps to plain cyclic coordinate descent and proves that the modified algorithm converges to a coordinate-wise minimum. This mission formalizes that convergence theorem.

Setting

Let $X\in\mathbb R^{n\times p}$ be a model matrix whose columns $X_1,\dots,X_p$ have unit Euclidean norm, and let $y\in\mathbb R^n$. For $\beta\in\mathbb R^p$ let

$$f(\beta)=\tfrac12\|y-X\beta\|^2+\lambda_1\|\beta\|_1+\lambda_2\|\beta\|_2^2,\qquad F(\beta)=f(\beta)+\lambda_0\|\beta\|_0,$$

with $\lambda_0>0$ and $\lambda_1,\lambda_2\ge0$; $\|\beta\|_0$ counts the nonzero coordinates. Problem (2) is $\min_\beta F(\beta)$. Its three named cases are (L0) ($\lambda_1=\lambda_2=0$), (L0L1) ($\lambda_1>0=\lambda_2$) and (L0L2) ($\lambda_1=0<\lambda_2$). $\mathrm{Supp}(\beta)=\{i:\beta_i\ne0\}$.

Fix all coordinates but the $i$-th. With unit-norm columns, minimizing $F$ over $\beta_i$ amounts to minimizing a scalar function of $\tilde\beta_i=\langle y-\sum_{j\ne i}X_j\beta_j,X_i\rangle$. The minimizer set is $\tilde T(\tilde\beta_i,\lambda_0,\lambda_1,\lambda_2)$ (7). The thresholding operator $T$ (12) picks one element of it:

$$T(a,\lambda_0,\lambda_1,\lambda_2)=\mathrm{sign}(a)\frac{|a|-\lambda_1}{1+2\lambda_2}\ \text{ if }\ \frac{|a|-\lambda_1}{1+2\lambda_2}\ge\sqrt{\frac{2\lambda_0}{1+2\lambda_2}},\qquad 0\ \text{ otherwise.}$$

A coordinate-wise (CW) minimum is a $\beta$ such that every $\beta_i$ minimizes $F$ in coordinate $i$ with the others fixed (Definition 2). A stationary solution is a $\beta$ whose lower directional derivative $F'(\beta;d)=\liminf_{\alpha\downarrow0}(F(\beta+\alpha d)-F(\beta))/\alpha$ is nonnegative in every direction $d$ (Definition 1).

Algorithm 1 (CDSS) takes an initial point $\beta^0$ and a positive integer $C$. It sweeps the coordinates cyclically. Each non-spacer step sets the current coordinate to $T(\tilde\beta_i,\lambda_0,\lambda_1,\lambda_2)$ and increments a counter for the resulting support. When the counter of a support $S$ reaches $Cp$, a spacer step follows: one sequential pass over the coordinates of $S$ with the $L_0$ term switched off, $\beta_i\leftarrow T(\tilde\beta_i,0,\lambda_1,\lambda_2)$. The counter of $S$ is then reset. The iterates $\beta^k$ count both kinds of step.

For (L0) and (L0L1) the paper adds Assumption 1: every $\min\{n,p\}$ columns of $X$ are linearly independent. It also adds Assumption 2, which constrains the initial point when $p>n$: $F(\beta^0)\le\lambda_0 n$ for (L0), and $F(\beta^0)\le\min_\beta f(\beta)+\lambda_0 n$ for (L0L1).

Formalization targets

Goal: Theorem 2 (p. 12)

For (L0), (L0L1) and (L0L2), under Assumptions 1–2 in the first two cases, there are $m$ and a support $S$ with

$$\mathrm{Supp}(\beta^k)=S\ \ (k\ge m),\qquad \beta^k\to B,\qquad B\text{ a CW minimum},\qquad \mathrm{Supp}(B)=S .$$

The two parts share $S$: the support stabilizes, and the limit is a CW minimum on that support.

Milestones

• Lemma 5, monotone descent: $F(\beta^{k+1})\le F(\beta^k)$ and $F(\beta^k)\downarrow F^*\ge0$.
• Lemma 6: $\|\beta^k\|_0\le\min\{n,p\}$ for (L0), (L0L1).
• Lemma 8: $\{\beta^k\}$ is bounded.
• Lemma 1: stationarity is equivalent to $\nabla_S f(\beta)=0$ on the support $S$.
• Lemma 9: along a support $S$ recurring in non-spacer steps, the spacer iterates eventually keep support $S$. Every subsequence with support $S$ converges to the unique minimizer $\beta^*$ of $f$ on $S$, and $|\beta^*_j|\ge\sqrt{2\lambda_0/(1+2\lambda_2)}$ on $S$.
• Lemma 10: the support of every limit point occurs infinitely often.
• Lemma 2: closed form of $\tilde T$.
• Lemma 11: for two limit points whose supports differ by one index $j$, $|B^{(2)}_j|$ exceeds $\sqrt{2\lambda_0/(1+2\lambda_2)}$, or equals it, according to whether $X_j$ is correlated with the smaller support.
• Lemma 12: zeroing a nonzero coordinate decreases $F$ by at least $\frac{1+2\lambda_2}{2}\big(|\beta^k_j|-\sqrt{2\lambda_0/(1+2\lambda_2)}\big)^2$.
• Lemma 3: the coordinatewise characterization (8) of CW minima, with the condition $|\tilde\beta_i|>\lambda_1$ of (5) on the support (see Formalization scope).

Significance

Theorem 2 guarantees that exact cyclic coordinate descent on Problem (2), with spacer steps, has a limit, and that the limit is a CW minimum. That class is strictly smaller than the stationary solutions and the fixed points of iterative hard thresholding ((4), p. 5). It also serves as the inner solver of the paper's local combinatorial search (Algorithm 2, Theorem 4). So the guarantee transfers to the partial-swap-inescapable minima studied in the companion mission. Support stabilization (part 1) is also what reduces the asymptotic analysis to smooth strongly convex coordinate descent on a fixed support (Theorem 3).

The result is proved in the paper; no machine-checked version exists. A formalization would check the convergence argument for a discontinuous objective, where the bookkeeping of supports, counters and step indices is easy to get wrong. It would also provide a reusable Lean model of $L_0$-penalized least squares and its thresholding operators.

Difficulty

The standard proof for cyclic coordinate descent needs a uniform decrease per step or a continuous objective. Neither holds here. A coordinate may switch between zero and a nonzero value infinitely often while $F$ decreases by arbitrarily small amounts, so the descent property alone gives neither support stabilization nor convergence of the iterates. Coordinates dropped infinitely often and coordinates added infinitely often must both be excluded. The tie-breaking rule in $T$ matters here: with the other choice on ties the algorithm is different. The exclusion of coordinates added infinitely often is only sketched in the paper (p. 43), and a complete proof has to supply it.

Formalization scope

Vectors are Fin p → ℝ with 0-based indices. Every norm is written as an explicit sum; Mathlib's sup norm on function types is never used. Unit-norm columns, $\lambda_0>0$, $\lambda_1,\lambda_2\ge0$ and $p\ge1$ are hypotheses. The case $\lambda_1,\lambda_2>0$, which the paper never names, is excluded from the statements about the algorithm's limits. $\mathrm{sign}$ is Real.sign. $T$ chooses the nonzero solution on ties, as in (12). The lower directional derivative takes values in EReal, since it is $+\infty$ in directions leaving the support. Lemma 3 is stated with the condition $|\tilde\beta_i|>\lambda_1$ on the support, taken from (5). The printed condition (8) omits it, and without it the equivalence fails for $\lambda_1>0$. Indices and pages throughout are those of arXiv:1803.01454v3.

Algorithm 1 is a deterministic state machine: iterate, cyclic pointer, the counter array Count, and a flag for a pending spacer step. The iterates are therefore defined from $\beta^0$ and $C$, not quantified over. The pointer advances only on non-spacer steps, following the pseudocode box; the prose formula $i=(k+1)\bmod(p+1)$ on p. 10 contradicts it. A statement about an arbitrary sequence satisfying a step predicate, or one that adds a step size or a sufficient-decrease hypothesis, would describe a different algorithm and is ruled out.

Useful infrastructure includes the scalar thresholding lemmas and coercivity of $f$ under Assumption 1 on supports of size at most $\min\{n,p\}$. Contributions of proofs for any milestone are welcome, as are lemmas on the state machine (descent of each step type, the counter bookkeeping).

Selected references

• H. Hazimeh, R. Mazumder, Fast Best Subset Selection: Coordinate Descent and Local Combinatorial Optimization Algorithms, Operations Research 68(5), 2020; arXiv:1803.01454v3. https://arxiv.org/abs/1803.01454
• D. Bertsimas, A. King, R. Mazumder, Best subset selection via a modern optimization lens, Annals of Statistics 44(2), 2016. https://doi.org/10.1214/15-AOS1388
• P. Tseng, Convergence of a block coordinate descent method for nondifferentiable minimization, J. Optim. Theory Appl. 109, 2001. https://doi.org/10.1023/A:1017501703105
• 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
• J. Friedman, T. Hastie, R. Tibshirani, Regularization paths for generalized linear models via coordinate descent, J. Stat. Softw. 33(1), 2010. https://doi.org/10.18637/jss.v033.i01