Numerical Analysis

Motivation

Many problems in statistics, signal processing and machine learning minimize a smooth loss plus a nonsmooth regularizer over a set with manifold structure: sparse principal component analysis over the Stiefel manifold, sparse blind deconvolution over the sphere, and clustering and dictionary-learning models with orthogonality constraints. In Euclidean space the standard method for such composite problems is the proximal gradient method, whose $O(1/k)$ rate for convex problems is classical (Beck–Teboulle 2009). Huang and Wei (arXiv:1909.06065v4, Mathematical Programming 2021) propose a Riemannian proximal gradient method (RPG) that solves its proximal subproblem on the tangent space and returns to the manifold through a retraction. They prove three things about it: global convergence, an $O(1/k)$ rate under a convexity notion adapted to retractions, and local rates under a Riemannian Kurdyka–Łojasiewicz property. This mission covers the second result, Theorem 3.2.

Earlier Riemannian proximal methods either required the exponential map and parallel transport, or solved a Euclidean proximal problem in the ambient space (Chen, Ma, So, Zhang 2020, ManPG). Theorem 3.2 is stated for a general retraction, with an explicit constant.

Setting

Let $\mathcal M$ be a finite-dimensional Riemannian manifold, with inner product $\langle\cdot,\cdot\rangle_x$ and norm $\|\cdot\|_x$ on each tangent space $T_x\mathcal M$. A retraction is a family of smooth maps $R_x:T_x\mathcal M\to\mathcal M$ with $R_x(0_x)=x$ and $\mathrm DR_x(0_x)=\mathrm{id}$. The problem is

$$\min_{x\in\mathcal M} F(x)=f(x)+g(x),$$

where $f$ is differentiable with Riemannian gradient $\operatorname{grad} f$, and $g$ is continuous and possibly nonsmooth.

For a function $h$ on a real inner-product space $V$, the Clarke generalized directional derivative is $h^\circ(\eta;v)=\limsup_{\xi\to\eta,\,t\downarrow0}(h(\xi+tv)-h(\xi))/t$. The generalized subdifferential is $\partial h(\eta)=\{\zeta\mid\langle\zeta,v\rangle\le h^\circ(\eta;v)\ \forall v\}$. A point $\eta$ is stationary for $h$ if $0\in\partial h(\eta)$.

Given a constant $\tilde L$, the RPG method (Algorithm 1) produces iterates $x_k$ and steps $\eta^*_{x_k}\in T_{x_k}\mathcal M$ as follows. Set

$$\ell_{x_k}(\eta)=\langle\operatorname{grad} f(x_k),\eta\rangle_{x_k}+\frac{\tilde L}{2}\|\eta\|_{x_k}^2+g(R_{x_k}(\eta)).$$

Choose a stationary point $\eta^*_{x_k}$ of $\ell_{x_k}$ on $T_{x_k}\mathcal M$ with $\ell_{x_k}(0)\ge\ell_{x_k}(\eta^*_{x_k})$, and set $x_{k+1}=R_{x_k}(\eta^*_{x_k})$.

The assumptions involve the sublevel set $\Omega_{x_0}=\{x\mid F(x)\le F(x_0)\}$ and two notions defined with respect to $R$ on a set $\mathcal N$:

• $L$-retraction-smoothness: $h(R_x(\eta))\le h(x)+\langle\operatorname{grad} h(x),\eta\rangle_x+\frac L2\|\eta\|_x^2$;
• retraction-convexity: $q_x=h\circ R_x$ satisfies $q_x(\eta)\ge q_x(\xi)+\langle\zeta,\eta-\xi\rangle_x$, where $\zeta$ is the gradient of $q_x$ at $\xi$, or any Riemannian subgradient of $q_x$ at $\xi$ when $h$ is nonsmooth.

Assumption 3.1 asks that $F$ be bounded below and $\Omega_{x_0}$ compact. Assumption 3.3 provides an open $\Omega\supseteq\Omega_{x_0}$ on which $f$ is $L$-retraction-smooth and retraction-convex and $g$ is retraction-convex. Assumption 3.4 bounds how far the retraction is from Euclidean: for all $x,y,z\in\Omega$,

$$\bigl|\|R_x^{-1}(z)-R_x^{-1}(y)\|_x^2-\|R_y^{-1}(z)\|_y^2\bigr|\le\kappa_\Omega\|R_x^{-1}(y)\|_x^2 .$$

Finally, $\beta=(\tilde L-L)/2$.

Formalization targets

Goal: Theorem 3.2

Let $\tilde L>\max(L,0)$, suppose Assumptions 3.1, 3.3 and 3.4 hold, and let $x_*$ be any accumulation point of $\{x_k\}$. Then for every $k\ge1$

$$F(x_k)-F(x_*)\le\frac1k\left(\frac{\tilde L}{2}\|R_{x_0}^{-1}(x_*)\|_{x_0}^2+\frac{\tilde L\kappa_\Omega}{2\beta}\bigl(F(x_0)-F(x_*)\bigr)\right).\tag{3.16}$$

The constants are those of the paper.

Milestones

1. Lemma 3.1 (descent): $F(x_k)-F(x_{k+1})\ge\beta\|\eta^*_{x_k}\|_{x_k}^2$.
2. Lemma 3.4 (Riemannian three-point inequality): for one step $z=R_x(\eta^*_x)$ and any $y=R_x(\xi_x)\in\Omega$, $F(z)\le F(y)+\frac{\tilde L}{2}(\|\xi_x\|_x^2-\|\xi_x-\eta^*_x\|_x^2)$.
3. (3.17): the one-step estimate obtained from Lemma 3.4 with $y=x_*$ and from Assumption 3.4.
4. (3.18): the one-step estimate with $\|\eta^*_{x_k}\|^2$ replaced by the decrease of $F$.
5. (3.19): the averaged estimate over the first $k$ steps.

Significance

Theorem 3.2 extends the $O(1/k)$ function-value rate of the proximal gradient method to manifolds and general retractions. The constant is explicit: it combines the initial distance $\|R_{x_0}^{-1}(x_*)\|$ with a correction proportional to $\kappa_\Omega$, which measures how far the retraction is from Euclidean. In the Euclidean case $R_x(\eta)=x+\eta$ one can take $\kappa_\Omega=0$, and the bound reduces to the classical one. The rate concerns $F(x_k)-F(x_*)$, which the iteration bound of Theorem 3.1 (the norm of $\eta^*_{x_k}$) does not control.

The result is proved in the paper. None of it is formalized, and no machine-checked version of the Euclidean proximal-gradient rate is on the platform either. A formalization would supply the first checked version of Clarke stationarity on tangent spaces, of the retraction-based convexity and smoothness notions, and of a rate proof for a nonsmooth Riemannian algorithm. These pieces carry over to other retraction-based methods, such as the accelerated variants of §4 of the paper.

Difficulty

The arithmetic of the rate proof is a telescoping sum, so the work lies elsewhere. First, Lemma 3.4 needs a usable form of the stationarity condition $0\in\partial\ell_x(\eta^*_x)$: a subgradient of $g\circ R_x$ at $\eta^*_x$ equal to $-(\operatorname{grad} f(x)+\tilde L\eta^*_x)$. This requires a sum rule for the Clarke subdifferential of a $C^1$ function plus a locally Lipschitz one, on a tangent space whose norm comes from the Riemannian metric rather than from the model space. Second, the gradient of $f\circ R_x$ at $0_x$ has to be identified with $\operatorname{grad} f(x)$ through $\mathrm DR_x(0_x)=\mathrm{id}$. In Mathlib this mixes manifold derivatives (mfderiv) with Fréchet derivatives on a space carrying two equivalent norms. Third, it has to be shown that the iterates and the accumulation point remain in $\Omega$, so that Assumptions 3.3 and 3.4 apply at every step.

Formalization scope

Theorem numbers and pages are those of arXiv:1909.06065v4. $\mathcal M$ is a smooth manifold modelled on a finite-dimensional real inner-product space $E$, with Mathlib's RiemannianBundle on the tangent bundle. All norms and inner products on $T_x\mathcal M$ are the Riemannian ones. Retractions and gradients are the published RiemOpt.BFGS.IsRetraction and RiemOpt.FR.IsGradient. The Clarke derivative is an EReal-valued limit superior along $\mathcal N(\eta)\times\mathcal N_{>0}(0)$. Algorithm 1 is modelled as a run (sequences $x_k$, $\eta^*_{x_k}$ with the three defining properties), and every statement holds for every run. Indices start at $0$, accumulation points are cluster points of the sequence, and the bounds divided by $k$ are stated for $k\ge1$.

The formalization commits to the following readings:

1. $L$-retraction-smoothness is required for every tangent vector at the points of the set. This is a strengthening of Definition 3.1, which only covers $\eta$ with $R_x(\eta)$ in the set, and it is needed because Lemma 3.1 fails under the literal reading.
2. In retraction-convexity the vector $\zeta$ depends on $\xi$. It is the gradient of $f\circ R_x$, respectively any Clarke subgradient of $g\circ R_x$, as the sentence after (3.10) states.
3. $\kappa_\Omega$ is a single constant.
4. $R^{-1}$ is data with $R_x(R_x^{-1}(y))=y$ on $\Omega$, together with the identification $R_{x_k}^{-1}(x_{k+1})=\eta^*_{x_k}$ that the proof uses.
5. Local Lipschitz continuity of $g\circ R_x$ is assumed, as p. 5 presupposes when it defines the subdifferential.

A single $\zeta$ for all $\eta,\xi$ in (3.10) would make retraction-convexity force $q_x$ to be affine, and the theorem would only cover affine problems. That reading is excluded. Statements with the case $k=0$, where Lean's $1/0=0$ makes (3.16) false, are excluded as well.

Contributions welcome: a Clarke sum rule for $C^1$ plus locally Lipschitz functions, the identification of the Fréchet derivative of $f\circ R_x$ at $0$ with the Riemannian gradient, and the five milestones in order. The Clarke layer and the derivative identification are reusable beyond this mission.

Selected references

• W. Huang, K. Wei, Riemannian proximal gradient methods, Mathematical Programming, 2021; extended version arXiv:1909.06065v4. https://arxiv.org/abs/1909.06065v4, https://doi.org/10.1007/s10107-021-01632-3
• A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imaging Sciences, 2009. https://doi.org/10.1137/080716542
• S. Chen, S. Ma, A. M.-C. So, T. Zhang, Proximal gradient method for nonsmooth optimization over the Stiefel manifold, SIAM J. Optimization, 2020. https://doi.org/10.1137/18M122457X
• F. H. Clarke, Optimization and Nonsmooth Analysis, Wiley, 1983. https://doi.org/10.1137/1.9781611971309
• P.-A. Absil, R. Mahony, R. Sepulchre, Optimization Algorithms on Matrix Manifolds, Princeton University Press, 2008. https://doi.org/10.1515/9781400830244

Motivation

Many problems in statistics, signal processing and machine learning minimize a smooth loss plus a nonsmooth regularizer over a set that is a manifold rather than a vector space: sparse principal component analysis over the Stiefel manifold, sparse blind deconvolution over the sphere, low-rank recovery over fixed-rank matrices. In Euclidean space the standard tool for $\min f + g$ with $f$ smooth and $g$ nonsmooth is the proximal gradient method, whose convergence theory (stationarity of limit points, an $O(1/\epsilon^2)$ bound on the number of iterations to reach an $\epsilon$-small step, an $O(1/k)$ rate under convexity) is classical (Beck–Teboulle 2009).

Huang and Wei (arXiv:1909.06065, Mathematical Programming 2021) proposed the Riemannian proximal gradient method (RPG): at each iterate the proximal subproblem is posed on the tangent space, the nonsmooth term is pulled back through a retraction, and the next iterate is obtained by retracting the computed step. Unlike the earlier ManPG method of Chen, Ma, So and Zhang (2020), which uses the Euclidean proximal mapping of an embedding space, RPG is defined intrinsically on any manifold with a retraction. Its global convergence theorem is the subject of this mission.

Setting

Let $\mathcal M$ be a finite-dimensional manifold with a smooth Riemannian metric $\langle\cdot,\cdot\rangle_x$ on each tangent space $T_x\mathcal M$, with norm $\|\cdot\|_x$. A retraction is a smooth map $R : T\mathcal M\to\mathcal M$, $(x,\eta)\mapsto R_x(\eta)$, with $R_x(0_x) = x$ and $DR_x(0_x) = \mathrm{id}$. The vector transport by differentiated retraction is $\mathcal T_{\eta_x} = DR_x(\eta_x) : T_x\mathcal M\to T_{R_x(\eta_x)}\mathcal M$; $\mathcal T^\sharp_{\eta_x}$ denotes its adjoint for the two metrics and $\mathcal T^{-\sharp}_{\eta_x}$ the adjoint of its inverse.

The problem is

$$\min_{x\in\mathcal M} F(x) = f(x) + g(x),$$

with $f$ continuously differentiable, Riemannian gradient $\operatorname{grad} f$, and $g$ continuous, possibly nonsmooth. For $h : \mathcal M\to\mathbb R$ the Riemannian Clarke subdifferential is defined through the pull-back $\hat h_x = h\circ R_x$ on the Hilbert space $T_x\mathcal M$:

$$\hat h_x^\circ(\eta; v) = \limsup_{\xi\to\eta,\ t\downarrow 0}\frac{\hat h_x(\xi+tv)-\hat h_x(\xi)}{t},\qquad \partial h(x) = \{\zeta\in T_x\mathcal M \mid \langle\zeta,v\rangle_x\le\hat h_x^\circ(0_x;v)\ \forall v\}.$$

A point $x$ is stationary if $0\in\partial F(x)$.

Algorithm 1 (RPG). Given $x_0$ and a constant $\tilde L > L$, for $k = 0,1,\dots$ let

$$\ell_{x_k}(\eta) = \langle\operatorname{grad} f(x_k),\eta\rangle_{x_k} + \frac{\tilde L}{2}\|\eta\|_{x_k}^2 + g(R_{x_k}(\eta)),$$

choose $\eta^*_{x_k}\in T_{x_k}\mathcal M$ that is a stationary point of $\ell_{x_k}$ with $\ell_{x_k}(0)\ge\ell_{x_k}(\eta^*_{x_k})$ (condition (3.1)), and set $x_{k+1} = R_{x_k}(\eta^*_{x_k})$.

Assumption 3.1. $F$ is bounded below and $\Omega_{x_0} = \{x \mid F(x)\le F(x_0)\}$ is compact. Assumption 3.2. $f$ is $L$-retraction-smooth on $\Omega_{x_0}$: $f(R_x(\eta))\le f(x) + \langle\operatorname{grad} f(x),\eta\rangle_x + \frac L2\|\eta\|_x^2$ (inequality (3.2)).

Formalization targets

Goal: Theorem 3.1

For every run of Algorithm 1: if $\eta^*_{x_k} = 0$ then $x_k$ is stationary. Under Assumptions 3.1 and 3.2, $\{x_k\}$ has an accumulation point, every accumulation point $x_*$ is stationary, and for every $\epsilon > 0$ some iterate with index at most

$$\frac{F(x_0)-F(x_*)}{\beta\epsilon^2},\qquad\beta = \frac{\tilde L - L}{2},$$

satisfies $\|\eta^*_{x_k}\|_{x_k}\le\epsilon$.

Milestones, in the order the proof uses them

• Lemma 3.1 (descent): $F(x_k) - F(x_{k+1})\ge\beta\|\eta^*_{x_k}\|_{x_k}^2$.
• (3.6): $\lim_k\|\eta^*_{x_k}\|_{x_k} = 0$.
• Lemma 3.2: for a continuous vector field $\xi$, $\|\xi_y - \mathcal T^{-\sharp}_{\eta_x}\xi_x\|_y\to 0$ as $\eta_x\to 0$ on the tangent bundle, $y = R_x(\eta_x)$.
• (3.8): at every iteration some $\zeta\in\partial g(x_{k+1})$ satisfies $\operatorname{grad} f(x_k)+\tilde L\eta^*_{x_k}+\mathcal T^\sharp_{R_{\eta^*_{x_k}}}\zeta = 0$.
• (3.9): where the transport is invertible, $\operatorname{grad} f(x_{k+1}) - \mathcal T^{-\sharp}(\operatorname{grad} f(x_k)+\tilde L\eta^*_{x_k}) = \operatorname{grad} f(x_{k+1})+\zeta\in\partial F(x_{k+1})$.

Significance

Theorem 3.1 is the basic convergence guarantee of RPG: the method is a descent method without any convexity of $g$, its limit points satisfy the first-order necessary condition $0\in\partial F$, and the number of iterations needed to make the step small is explicit, with the same $O(1/\epsilon^2)$ dependence as in the Euclidean case. It underlies the further analyses of the same paper (the $O(1/k)$ rate under retraction convexity and the Kurdyka–Łojasiewicz analysis) and of later Riemannian proximal methods.

The theorem is proved on paper; no machine-checked version exists. Formalizing it requires a Clarke subdifferential on tangent spaces, its calculus (sum rule with a $C^1$ function, chain rule through a local diffeomorphism) and its closedness along convergent sequences on a manifold. None of these is in Mathlib. They are reusable well beyond this paper by any Riemannian nonsmooth optimization method.

Difficulty

The descent inequality and the iteration bound are elementary once Lemma 3.1 is available. The difficulty is the stationarity of accumulation points. The subproblem's optimality condition (3.8) lives in the tangent space at $x_k$ and involves a subgradient of $g$ at $x_{k+1}$; to pass to the limit, the gradients and subgradients must be compared across different tangent spaces, which is what the inverse-adjoint transport and Lemma 3.2 do. This needs joint continuity of the retraction and its differential on the tangent bundle and continuity of the metric, and then a closedness property of the Riemannian Clarke subdifferential that the paper cites from Hosseini, Huang and Yousefpour (2018, Theorem 2.2(c)) and that a complete development must prove. A pointwise argument at a fixed base point does not suffice, because the base points $x_{k_j}$ move.

Formalization scope

Theorem numbers and pages are those of arXiv:1909.06065v4. The manifold is modelled on a finite-dimensional real inner-product space $E$ and carries a RiemannianBundle with a smooth metric (IsContMDiffRiemannianBundle). Retraction, differentiated-retraction transport and Riemannian gradient are the published definitions RiemOpt.BFGS.Setting and RiemOpt.FR.Setting; joint smoothness of $R$ on the tangent bundle is added, as on p. 4. The Clarke limit superior is taken in EReal along $\mathcal N(\eta)\times\mathcal N_{>0}(0)$, so it never takes a default value. Accumulation points are MapClusterPt, and "in at most $N$ iterations" is an index $k\le N$. A run of Algorithm 1 is any sequence of iterates and steps satisfying (3.1); the step is a selection, not a function of the iterate.

Two hypotheses are explicit in the Lean statements. Each pull-back $g\circ R_x$ is locally Lipschitz, which p. 5 requires for the subdifferential to be defined. $f$ is $C^1$, because the proof applies Lemma 3.2 to $\operatorname{grad} f$. Assumption 3.2 is stated in the form the proof of Lemma 3.1 uses: (3.2) at every $x\in\Omega_{x_0}$ for every tangent vector. Definition 3.1 read literally only constrains $\eta$ with $R_x(\eta)\in\Omega_{x_0}$, and under that reading Lemma 3.1 fails: a nonconvex $g$ can send $x_1$ far outside $\Omega_{x_0}$. The paper's own sufficient condition gives the strengthened form.

The formalization must not trivialize stationarity. The subdifferential is defined from the Clarke derivative of the pull-back alone. It does not involve $\mathcal T^{-\sharp}$ or ContinuousLinearMap.inverse, so $0\in\partial F(x)$ is a genuine condition. The run hypotheses are satisfiable, for instance on $\mathcal M = \mathbb R$ with $R_x(\eta) = x+\eta$, $f(x) = x^2$, $g = 0$, $x_k = 0$, $\eta_k = 0$.

Welcome contributions: the Clarke calculus on a real Hilbert space (sum rule, chain rule through a $C^1$ map, upper semicontinuity), local invertibility of $DR_x(\eta)$ near the zero section, and the closedness of the Riemannian subdifferential on the tangent bundle.

Selected references

• W. Huang, K. Wei, Riemannian proximal gradient methods, Mathematical Programming, 2021. arXiv:1909.06065v4. https://arxiv.org/abs/1909.06065, https://doi.org/10.1007/s10107-021-01632-3
• S. Hosseini, W. Huang, R. Yousefpour, Line search algorithms for locally Lipschitz functions on Riemannian manifolds, SIAM Journal on Optimization 28(1):596–619, 2018 (reference [29] of Huang–Wei; source of the Riemannian Clarke subdifferential and of its closedness, Theorem 2.2(c)). https://epubs.siam.org/toc/sjope8/28/1
• P.-A. Absil, R. Mahony, R. Sepulchre, Optimization Algorithms on Matrix Manifolds, Princeton University Press, 2008. https://doi.org/10.1515/9781400830244
• 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 (2020). https://doi.org/10.1137/18M122457X
• A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM Journal on Imaging Sciences 2 (2009). https://doi.org/10.1137/080716542
• F. H. Clarke, Optimization and Nonsmooth Analysis, Wiley, 1983. https://doi.org/10.1137/1.9781611971309

Why Sinkhorn's algorithm needs a convergence rate

Sinkhorn's algorithm is the standard method for computing entropically regularized optimal transport between two histograms. It is used throughout machine learning, statistics, imaging and graphics because each iteration is two matrix–vector products followed by entrywise divisions, which parallelize on GPUs and differentiate easily. The same iteration, under the names matrix scaling, RAS or iterative proportional fitting, has been used since the 1930s in transportation planning, economics and contingency-table analysis. Whether the iteration converges, and how fast, determines how many iterations a practitioner must run and which stopping criteria are reliable.

This mission formalizes the global convergence analysis of Sinkhorn's algorithm as presented in Chapter 4 (§4.2, Remarks 4.12–4.14) of G. Peyré and M. Cuturi, Computational Optimal Transport (2019).

Timeline. Sinkhorn (1964) proved that a positive matrix can be scaled to be doubly stochastic, and Sinkhorn and Knopp (1967) extended this to prescribed marginals. Birkhoff (1957) and Samelson (1957) independently showed that a positive linear map is a strict contraction of the positive cone for Hilbert's projective metric, with an explicit ratio. Franklin and Lorenz (1989) combined the two to prove that Sinkhorn's scalings converge linearly in Hilbert's metric at rate $\lambda(K)^2$. Linial, Samorodnitsky and Wigderson (1998) related the iteration to the permanent.

Setting

Fix $n, m \ge 1$, a matrix $K \in \mathbb{R}^{n\times m}$ with positive entries (in the book, the Gibbs kernel $K_{i,j} = e^{-C_{i,j}/\varepsilon}$), and vectors $a \in \mathbb{R}^n$, $b \in \mathbb{R}^m$ with positive entries. Products $\odot$ and quotients of vectors are entrywise. The regularized transport plan has the form $P = \mathrm{diag}(u)K\mathrm{diag}(v)$, where the scalings $u, v$ solve

$$u \odot (Kv) = a, \qquad v \odot (K^\top u) = b. \tag{4.14}$$

Sinkhorn's algorithm alternately enforces each equation:

$$v^{(0)} = \mathbb{1}_m,\qquad u^{(\ell+1)} = \frac{a}{Kv^{(\ell)}},\qquad v^{(\ell+1)} = \frac{b}{K^\top u^{(\ell+1)}}. \tag{4.15}$$

For vectors with positive entries, Hilbert's projective metric is

$$d_{\mathcal H}(u,u') = \log\max_{i,j}\frac{u_i u'_j}{u_j u'_i},$$

which vanishes exactly when $u$ and $u'$ are proportional, and equals the variation seminorm $\|\log u - \log u'\|_{\mathrm{var}}$ with $\|f\|_{\mathrm{var}} = \max_i f_i - \min_i f_i$. Birkhoff's constants of $K$ are

$$\eta(K) = \max_{i,j,k,\ell}\frac{K_{i,k}K_{j,\ell}}{K_{j,k}K_{i,\ell}},\qquad \lambda(K) = \frac{\sqrt{\eta(K)}-1}{\sqrt{\eta(K)}+1}.$$

Write $P^{(\ell)} = \mathrm{diag}(u^{(\ell)})K\mathrm{diag}(v^{(\ell)})$ and $P^\star = \mathrm{diag}(u^\star)K\mathrm{diag}(v^\star)$ for a positive solution $(u^\star, v^\star)$ of (4.14).

Formalization targets

Goal: Theorem 4.2 (Franklin–Lorenz)

For every $\ell$,

$$d_{\mathcal H}(v^{(\ell)}, v^\star) \le \lambda(K)^{2\ell}\, d_{\mathcal H}(v^{(0)}, v^\star),\qquad d_{\mathcal H}(u^{(\ell+1)}, u^\star) \le \lambda(K)^{2\ell+1}\, d_{\mathcal H}(v^{(0)}, v^\star),$$

so both distances tend to $0$; for $\ell \ge 1$,

$$d_{\mathcal H}(u^{(\ell)}, u^\star) \le \frac{d_{\mathcal H}(P^{(\ell)}\mathbb{1}_m, a)}{1-\lambda(K)^2},\qquad \max_{i,j}\big|\log P^{(\ell)}_{i,j} - \log P^\star_{i,j}\big| \le d_{\mathcal H}(u^{(\ell)},u^\star) + d_{\mathcal H}(v^{(\ell)},v^\star),$$

and for $\ell \ge 0$ the mirrored bound for $v^{(\ell)}$ with the half-step coupling $\mathrm{diag}(u^{(\ell+1)})K\mathrm{diag}(v^{(\ell)})$. The book writes the rate as $O(\lambda(K)^{2\ell})$; the constants stated are those its proof yields.

Milestones, in attack order

1. (4.21): $d_{\mathcal H}(u,u') = \|\log u - \log u'\|_{\mathrm{var}}$.
2. Remark 4.12: $d_{\mathcal H}$ is a distance on rays (triangle inequality; zero exactly on proportional vectors).
3. Invariances: $d_{\mathcal H}(v,v') = d_{\mathcal H}(v/v',\mathbb{1}) = d_{\mathcal H}(\mathbb{1}/v,\mathbb{1}/v')$ and invariance under positive diagonal scaling.
4. Theorem 4.1 (Birkhoff–Hopf): $d_{\mathcal H}(Kv,Kv') \le \lambda(K)\,d_{\mathcal H}(v,v')$ and $\lambda(K) < 1$.
5. One half-step contracts by $\lambda(K)$: $d_{\mathcal H}(u^{(\ell+1)},u^\star) \le \lambda(K)\,d_{\mathcal H}(v^{(\ell)},v^\star)$, and the $v$-analogue.
6. The display behind (4.23): $d_{\mathcal H}(u^{(\ell)},u^\star) \le d_{\mathcal H}(a, u^{(\ell)}\odot Kv^{(\ell)}) + \lambda(K)^2 d_{\mathcal H}(u^{(\ell)},u^\star)$, and its mirror.

Significance

The result. Theorem 4.2 gives a global linear rate that holds from the first iteration, independent of the starting point, with a ratio computable from $K$ alone. Through (4.21) the rate transfers to the dual potentials $\varepsilon\log u^{(\ell)}$, and the bounds (4.23) turn the marginal violations $P^{(\ell)}\mathbb{1}_m$ versus $a$ — which are computed anyway — into certified stopping criteria. Theorem 4.1 alone also yields a quantitative Perron–Frobenius theorem (Remark 4.13).

Formalizing it. All results here are classical and proved in the literature (Birkhoff 1957; Franklin–Lorenz 1989); Theorem 4.1 is quoted by the book without proof. To our knowledge neither Hilbert's projective metric on the positive cone, the Birkhoff–Hopf theorem, nor the convergence rate of Sinkhorn's algorithm has a machine-checked proof in Mathlib. The mission produces a reusable development of $d_{\mathcal H}$ and Birkhoff's contraction theorem, and a certified rate and stopping rule for the most widely used optimal-transport solver.

Difficulty

Milestones 1–3 and 5–6 are finite manipulations of maxima of ratios. The central difficulty is Theorem 4.1. The obvious argument — bound $(Kv)_i/(Kv')_i$ entrywise by $\max_k v_k/v'_k$ — only shows that $K$ is non-expansive ($\lambda \le 1$); it loses the strict contraction entirely. The sharp ratio $(\sqrt{\eta}-1)/(\sqrt{\eta}+1)$ depends on how $K$ acts on the two vectors jointly, which the entrywise argument does not see. The book quotes the theorem without proof, so the argument has to be supplied from Birkhoff's original paper or later expositions.

Formalization scope

Lean namespace CompOT.Sinkhorn. Indices are $0$-based (Fin n, Fin m); vectors are functions Fin n → ℝ, matrices Matrix (Fin n) (Fin m) ℝ, and $\odot$ and $/$ are the pointwise operations on functions. $d_{\mathcal H}$, $\eta$ and $\|\cdot\|_{\mathrm{var}}$ use real iSup/iInf over the finite index sets. Every theorem assumes $0 < n$, $0 < m$ and positive entries for $K$, $a$, $b$ and every vector passed to $d_{\mathcal H}$, which excludes Lean's junk values for division by zero and logarithms of non-positive numbers. Sinkhorn's algorithm is a run predicate (IsSinkhornRun) stating the recursion (4.15) with $v^{(0)} = \mathbb{1}_m$; $u^{(0)}$ is not an iterate and is unconstrained. The solution $(u^\star, v^\star)$ is any positive solution of (4.14), taken as a hypothesis (it exists iff $\sum_i a_i = \sum_j b_j$); the book's histograms are probability vectors, but only positivity is used. $P^\star$ in (4.24) is $\mathrm{diag}(u^\star)K\mathrm{diag}(v^\star)$, the solution of the regularized problem (4.2) by Proposition 4.3.

Deviations from the print, all disclosed in the items: the $O(\cdot)$ of (4.22) is replaced by the proof's explicit constants; the limit $(u^{(\ell)}, v^{(\ell)}) \to (u^\star, v^\star)$ is stated as $d_{\mathcal H} \to 0$, since scalings are determined only up to $(ru, v/r)$; the second inequality of (4.23) is false as printed (its right side is $0$ for $\ell \ge 1$) and is stated with the half-step coupling the mirrored proof gives.

A definition of $d_{\mathcal H}$ without the logarithm, or as a maximum of single-entry ratios, would make the statements different theorems; the definition here is the book's, and (4.21) is a milestone checking it against the variation seminorm. A run predicate not tied to $K$, $a$, $b$, or an unsatisfiable one, is ruled out: the predicate is the recursion itself and is inhabited for every $K$, $a$, $b$.

Contributions welcome: the metric properties of $d_{\mathcal H}$ (reusable for Perron–Frobenius and nonlinear consensus), a proof of the Birkhoff–Hopf theorem, and positivity lemmas for Sinkhorn iterates.

Selected references

• G. Peyré, M. Cuturi, Computational Optimal Transport, Foundations and Trends in Machine Learning 11(5–6):355–607, 2019. https://doi.org/10.1561/2200000073 (§4.2, pp. 431–443)
• J. Franklin, J. Lorenz, On the scaling of multidimensional matrices, Linear Algebra and its Applications 114/115:717–735, 1989. https://doi.org/10.1016/0024-3795(89)90490-4
• G. Birkhoff, Extensions of Jentzsch's theorem, Transactions of the AMS 85(1):219–227, 1957. https://doi.org/10.1090/S0002-9947-1957-0087058-6
• H. Samelson, On the Perron–Frobenius theorem, Michigan Mathematical Journal 4(1):57–59, 1957. https://doi.org/10.1307/mmj/1028990177
• R. Sinkhorn, A relationship between arbitrary positive matrices and doubly stochastic matrices, Annals of Mathematical Statistics 35(2):876–879, 1964. https://doi.org/10.1214/aoms/1177703591
• N. Linial, A. Samorodnitsky, A. Wigderson, A deterministic strongly polynomial algorithm for matrix scaling and approximate permanents, STOC 1998. https://doi.org/10.1145/276698.276880

Motivation

First-order methods use function gradients to minimize an objective without solving a second-order subproblem at each iteration. Inertial methods also use the previous iterate, which can improve the convergence rate but complicates the choice of momentum and step size. Attouch, Chbani, Fadili, and Riahi study a family of such methods derived from continuous dynamical systems with Hessian driven damping. Their explicit strongly convex scheme, IGAHD-SC, uses a difference of consecutive gradients in addition to momentum and the current gradient. The difference represents the damping contribution without evaluating a Hessian matrix. Attouch et al., §5.2, pp. 28–29.

The question in this mission is quantitative. If the objective has strong curvature and a Lipschitz gradient, do the iterates, objective values, and a discounted sum of squared gradients decrease at the geometric rates claimed in Theorem 11? This is the final independent convergence theorem of the paper's series of discrete algorithms. Attouch et al., Theorem 11, p. 29.

Setting

Let $H$ be a real Hilbert space and $f:H\to\mathbb R$ a continuously differentiable objective. For $\mu>0$, the paper calls $f$ $\mu$-strongly convex when $z\mapsto f(z)-\frac\mu2\|z\|^2$ is convex. Let $x^*$ be a minimizer of $f$. The paper assumes that the gradient is $L$-Lipschitz for a positive constant $L$: $\|\nabla f(u)-\nabla f(w)\|\le L\|u-w\|$ for every $u,w\in H$. Strong convexity makes the minimizer unique. Attouch et al., Definition 1, p. 19; Theorem 11, p. 29.

Choose parameters $\beta\ge0$ and $s>0$, and write $r=\sqrt{\mu s}$. Starting from arbitrary $x_0,x_1\in H$, IGAHD-SC determines $x_{k+1}$ for $k\ge1$ by

$$x_{k+1}=x_k+\frac{1-r}{1+r}(x_k-x_{k-1}) -\frac{\beta\sqrt s}{1+r}\bigl(\nabla f(x_k)-\nabla f(x_{k-1})\bigr) -\frac{s}{1+r}\nabla f(x_k).$$

The two differences have different roles: $x_k-x_{k-1}$ is the inertial part, and $\nabla f(x_k)-\nabla f(x_{k-1})$ is the Hessian damping term in this explicit discretization. The algorithm is exactly the displayed scheme on p. 28, with its two starting points free. Attouch et al., (25) and (IGAHD-SC), p. 28.

Formalization targets

Theorem 11 assumes $\beta\le1/\sqrt\mu$ and the two smoothness bounds in (26). Define

$$q=\frac{1}{1+\frac12\sqrt{\mu s}},\qquad \theta=\frac{1}{1+\sqrt{\mu s}}.$$

Both constants lie strictly between zero and one. The mission's goal is the paper's three-part conclusion:

$$\|x_k-x^*\|=O(q^{k/2}),\qquad f(x_k)-f(x^*)=O(q^k),$$

and

$$\theta^k\sum_{j=0}^{k-2}\theta^{-j}\|\nabla f(x_j)\|^2=O(q^k) \quad (k\to\infty).$$

The milestones record the auxiliary vector and energy from the proof, their velocity and strong-convexity estimates, the gradient majorization, the one-step energy inequality, and the geometric energy bound. These are claims displayed by the authors on pp. 30–32, in the order in which they lead to the goal. Attouch et al., proof of Theorem 11, pp. 30–32.

Significance

The first two bounds give geometric decay of distance to the optimizer and of objective error for the specified parameter regime. The third controls a discounted accumulation of squared gradients. Together they describe more than convergence of the objective alone: the same run has explicit asymptotic rates for its points and gradient behavior. The rate constants depend on the strong-convexity parameter, step size, and allowed damping range. Attouch et al., Theorem 11, p. 29.

This result is proved in the cited paper; the open work here is its Lean formalization. A completed development would include a reusable representation of the explicit inertial run, a precise strong-convexity interface matching Definition 1, and the energy estimates that support the theorem. The formal statements in this proposal are proof targets, not claims that already have machine-checked proofs. The two comparison methods mentioned in Remark 10 use different discretizations or hypotheses, so their rate theorems cannot replace this target. Attouch et al., Remark 10, p. 29.

Difficulty

Ordinary smooth gradient descent is controlled by the current gradient and one step difference. Here the recurrence also contains the preceding gradient and the preceding step. A direct estimate of $f(x_{k+1})-f(x^*)$ does not isolate a contracting quantity. The paper therefore tracks an auxiliary vector and an energy involving both the objective gap and that vector; the precise coefficients in the two bounds of (26) are needed for the final one-step inequality. The proof's intermediate display on p. 31 also changes a gradient index between lines, so each formalized displayed claim needs scrutiny against the recurrence and the final p. 32 inequality. Attouch et al., proof of Theorem 11, pp. 30–32.

Formalization scope

Lean represents $H$ as a complete real inner-product space, $\nabla f$ by Mathlib's gradient, and the run as a predicate on a sequence $x:\mathbb N\to H$ obeying (IGAHD-SC) for every $k\ge1$. The source's arbitrary initial pair is preserved. A witness $x^*$ and the hypothesis that it minimizes $f$ express the nonempty argmin assumption without taking an infimum of real values. The strong-convexity definition is exactly Definition 1. Positive $s$, $\mu$, and $L$ make every square root and step-size denominator meaningful. The first bound in (26) is written $8\beta L\le\sqrt\mu$, which retains the paper's intended $\beta=0$ case; the paper later uses $0\le\beta$ explicitly, so it is included. Attouch et al., pp. 19, 28–32.

The rates use =O[Filter.atTop]. Lean writes $q^{k/2}$ as $(\sqrt q)^k$, and the sum through $k-2$ as Finset.range (k - 1). For $k=0,1$ the finite range is empty; the asymptotic claim concerns large $k$. The source prints $p$ in the sum limit and $j$ in its summand; both are read as $j$. The goal itself mentions the iterates, objective values, and gradients, without assuming an energy decay bound. A complete proof can contribute reusable inequalities for strongly convex functions with Lipschitz gradients and geometric recurrences, alongside the paper-specific run and energy lemmas.

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.

Motivation

Many first-order optimization methods are discretizations of second-order differential equations. In the continuous system the velocity carries the momentum, and the damping terms decide how fast the trajectory settles at a minimizer. Attouch, Chbani, Fadili and Riahi study inertial systems with Hessian driven damping: a friction term proportional to $\nabla^2 f(x(t))\dot x(t)$, which damps the oscillations that plain heavy-ball and Nesterov-type dynamics show. Because $\nabla^2 f(x(t))\dot x(t)$ is the time derivative of $\nabla f(x(t))$, a discretization can use differences of gradients and never evaluate a Hessian. Attouch et al., §1, pp. 1–3.

For a $\mu$-strongly convex objective, the paper shows that the dynamic

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

makes $f(x(t))-\min f$ decay like $e^{-\frac{\sqrt\mu}{2}t}$ (Theorem 7, p. 19). Section 5.1 asks whether that rate survives an implicit (proximal) discretization. The answer is Theorem 9: the algorithm converges linearly, and the gradients also go to zero at a geometric rate. The authors remark that they know of no earlier result of this kind for such a proximal algorithm. Attouch et al., Theorem 9 and Remark 9, p. 23.

Setting

Let $H$ be a real Hilbert space with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$, and let $f:H\to\mathbb R$ be convex and continuously differentiable, with gradient $\nabla f$. For $\mu>0$, $f$ is $\mu$-strongly convex when $z\mapsto f(z)-\frac\mu2\|z\|^2$ is convex (Definition 1, p. 19). Such an $f$ has at most one minimizer; write $x^\star$ for it.

For $\gamma>0$ and $y\in H$, the proximal point $\operatorname{prox}_{\gamma f}(y)$ is the minimizer of $z\mapsto\gamma f(z)+\frac12\|z-y\|^2$.

Fix $\beta\ge0$ and a step $s>0$ (the paper writes $s=h^2$ with $h$ the step size). Starting from arbitrary $x_0,x_1\in H$, the inertial proximal algorithm with Hessian damping for strongly convex functions, (IPAHD-SC), produces for $k\ge1$

$$\begin{cases} y_k=x_k+\Big(1-\frac{2\sqrt{\mu s}}{1+2\sqrt{\mu s}}\Big)(x_k-x_{k-1})+\beta\sqrt s\Big(1-\frac{2\sqrt{\mu s}}{1+2\sqrt{\mu s}}\Big)\nabla f(x_k),\\[4pt] x_{k+1}=\operatorname{prox}_{\frac{\beta\sqrt s+s}{1+2\sqrt{\mu s}}f}(y_k). \end{cases}$$

The first line is an extrapolation with momentum and an explicit gradient correction; the second is one implicit (proximal) step. Attouch et al., (IPAHD-SC), p. 23.

The rates are

$$q=\frac{1}{1+\frac12\sqrt{\mu s}},\qquad \theta=\frac{1}{1+\sqrt{\mu s}},$$

both in $]0,1[$, and the constant of the theorem is

$$E_1=f(x_1)-f(x^\star)+\tfrac12\Big\|\sqrt\mu(x_1-x^\star)+\tfrac{1}{\sqrt s}(x_1-x_0)+\beta\nabla f(x_1)\Big\|^2.$$

Formalization targets

Goal: Theorem 9

Assume $0\le\beta\le\frac{1}{2\sqrt\mu}$ and $\sqrt s\le\beta$. Then every sequence generated by (IPAHD-SC) satisfies, for every $k\ge1$,

$$\frac\mu2\|x_k-x^\star\|^2\le f(x_k)-\min_H f\le E_1q^{k-1},$$

and

$$\theta^k\sum_{j=0}^{k-2}\theta^{-j}\|\nabla f(x_j)\|^2=\mathcal O(q^k)\qquad(k\to+\infty).$$

Milestones

The milestone targets are the paper's displayed intermediate claims, in order:

1. the equivalent second-order form (20) of the algorithm;
2. the increment of the auxiliary vector $v_k=\sqrt\mu(x_k-x^\star)+\frac{1}{\sqrt s}(x_k-x_{k-1})+\beta\nabla f(x_k)$;
3. the two "elementary algebra" quadratic-form inequalities, Term 1 and Term 2;
4. the one-step inequality $\frac{1}{\sqrt s}(E_{k+1}-E_k)+\frac12\sqrt\mu E_{k+1}\le0$ for $E_k=f(x_k)-f(x^\star)+\frac12\|v_k\|^2$;
5. the energy decay (21), $E_k\le E_1q^{k-1}$;
6. the recursive inequality (22) for $Z_k=2\beta(f(x_k)-f(x^\star))+\sqrt\mu\|x_k-x^\star\|^2$;
7. the discounted gradient-sum bound (24).

Significance

Theorem 9 gives linear convergence of a proximal algorithm for a smooth strongly convex objective with no Lipschitz assumption on $\nabla f$ and no upper bound on the step in terms of a smoothness constant: the step condition $\sqrt s\le\beta\le\frac1{2\sqrt\mu}$ involves only the strong convexity modulus. It also gives a geometric rate for a weighted sum of squared gradients, which the authors present as new for proximal algorithms of this type. In §5.1.2 the same estimate, applied to the Moreau envelope, gives linear convergence for nonsmooth strongly convex objectives (Theorem 10), and §5.2 shows the explicit gradient counterpart (Theorem 11).

The result is proved in the paper. A formal development would certify the exact constants ($q$, $\theta$, and the factor $4E_1/\sqrt\mu$ in (24)) and the index bookkeeping of the gradient sums, where the printed proof contains small slips. It would also provide reusable proximal optimality and strong convexity results for later formalizations.

Difficulty

The continuous-time estimate of Theorem 7 does not transfer term by term to the discrete energy $E_k$. Its difference contains cross terms such as $\beta\mu\langle\nabla f(x_{k+1}),x_{k+1}-x^\star\rangle$ and $\sqrt\mu\langle\nabla f(x_{k+1}),x_{k+1}-x_k\rangle$, whose signs are not automatic. The upper bound on $\beta$ and the lower bound on $\sqrt s$ must both survive the resulting estimates. The gradient conclusion introduces a separate weighted sum, with different index ranges in Theorem 9 and the final displayed line of its proof. The algorithm begins at $k=1$ and uses two free initial points, $x_0$ and $x_1$.

Formalization scope

The Lean development works on a real Hilbert space H (NormedAddCommGroup, InnerProductSpace ℝ, CompleteSpace), with $\nabla f$ Mathlib's gradient f. Conventions:

• Hypotheses. $f$ is ConvexOn ℝ Set.univ f and ContDiff ℝ 1 f, and $\mu$-strongly convex in the literal sense of Definition 1 (HessianDamping.DINSC.IsStronglyConvex, convexity of $f-\frac\mu2\|\cdot\|^2$), with $0<\mu$. The standing assumption $\operatorname{argmin}f\neq\emptyset$ is a point $x^\star$ with $f(x^\star)\le f(y)$ for every $y$, so $\min_H f$ is $f(x^\star)$. The step $s$ is positive, as the page implies ($s=h^2$, $h>0$). The parameter hypotheses are $0\le\beta\le\frac1{2\sqrt\mu}$ and $\sqrt s\le\beta$, as printed.
• Algorithm. The proximal step is the published predicate GoldenRatioVI.Shared.IsProxPoint: $x_{k+1}$ minimizes $z\mapsto\gamma f(z)+\frac12\|z-y_k\|^2$. No proximal function is defined. A run imposes the step for every $k\ge1$; $x_0$ and $x_1$ are free.
• Rates and sums. $q$ and $\theta$ are written with Real.sqrt (μ * s). $\theta^{-j}$ is an integer power. $\sum_{j=0}^{k-2}$ is a sum over Finset.range (k - 1). $\mathcal O(q^k)$ is =O[Filter.atTop] with the constant allowed to depend on all the data. The side claims $0<q<1$ and $0<\theta<1$ are a conjunct of the goal.
• Not trivial. The goal does not assume the energy decay, and it does not state the result for a single parameter choice or a vacuous hypothesis set: the hypotheses hold, for example, for $f=\frac12\|\cdot\|^2$, $\mu=1$, $s=\frac1{16}$, $\beta=\frac14$. The goal also does not mention $v_k$, $E_k$ or $Z_k$; only the printed constant $E_1$ appears.

Infrastructure needed: the first-order optimality condition of a proximal step for a differentiable function, and the first-order characterization of strong convexity, $f(y)\ge f(x)+\langle\nabla f(x),y-x\rangle+\frac\mu2\|y-x\|^2$, from Definition 1. Both are reusable well beyond this mission, as is the iteration of linear recursive inequalities $a_k\le qa_{k-1}+b_k$. Contributions of either kind, and proofs of the individual milestones, are welcome.

Selected references

• H. Attouch, Z. Chbani, J. Fadili, H. Riahi, First-order optimization algorithms via inertial systems with Hessian driven damping, Mathematical Programming, 2020; preprint arXiv:1907.10536v2. https://arxiv.org/abs/1907.10536
• H. Attouch, J. Peypouquet, P. Redont, Fast convex minimization via inertial dynamics with Hessian driven damping, Journal of Differential Equations 261(10) (2016), 5734–5783. https://doi.org/10.1016/j.jde.2016.08.020
• B. T. Polyak, Some methods of speeding up the convergence of iteration methods, USSR Computational Mathematics and Mathematical Physics 4 (1964), 1–17. https://doi.org/10.1016/0041-5553(64)90137-5

Motivation

The ℓ1-regularized least squares problem (the Lasso, or basis pursuit denoising) is the model problem of compressed sensing and sparse inverse problems: given a matrix $A\in\mathbb R^{m\times n}$, data $d\in\mathbb R^m$ and a weight $\mu>0$, one seeks a vector $x$ that fits $Ax\approx d$ while having few nonzero entries. The standard first-order method for it is the iterative shrinkage thresholding algorithm (ISTA), a forward-backward splitting scheme whose backward step is componentwise soft thresholding. Each iteration costs two matrix–vector products, which is why ISTA and its variants are used on problems far too large for interior-point methods.

The classical complexity guarantee for ISTA is $f(x_k)-\min f=O(1/k)$, the rate of the proximal gradient method on any convex composite problem (Beck–Teboulle 2009). ISTA itself goes back to Daubechies–Defrise–De Mol 2004, and its link with forward-backward splitting to Combettes–Wajs 2005. An asymptotic linear rate was known, with an explanation through partial smoothness (Liang–Fadili–Peyré), but there was no global, explicit, non-asymptotic bound. Bolte, Nguyen, Peypouquet and Suter (arXiv:1510.08234v3, Math. Program. 2017) derive such a bound from two ingredients: an error bound for the Lasso objective on $\ell_1$ balls, with a constant written through a Hoffman constant, and their general theory turning a Kurdyka–Łojasiewicz (KL) inequality into complexity bounds for descent methods.

Setting

Write $\mathbb R^n$ with the Euclidean norm $\|\cdot\|$, and $\|x\|_1=\sum_i|x_i|$. For a matrix $M$, $\|M\|$ is its spectral norm, the operator norm from $\mathbb R^n$ to $\mathbb R^m$ with Euclidean norms. The objective is

$$f(x)=\mu\|x\|_1+\tfrac12\|Ax-d\|^2,\qquad \mu>0,$$

with $\min f=\inf_x f(x)$ and $S=\operatorname{argmin} f$. Put $g=\mu\|\cdot\|_1$ and $h=\tfrac12\|A\cdot-d\|^2$. Then $h$ is convex with gradient $A^T(Ax-d)$, Lipschitz with constant $L=\|A^TA\|$.

Given stepsizes $0<\lambda^-\le\lambda_k\le\lambda^+$ with $\lambda^+L<2$ and any $x_0$, ISTA is

$$x_{k+1}=\operatorname{prox}_{\lambda_k\mu\|\cdot\|_1}\big(x_k-\lambda_k(A^TAx_k-A^Td)\big),\qquad k\ge0,$$

where $\operatorname{prox}_\varphi(v)$ is the minimizer of $\varphi(y)+\tfrac12\|y-v\|^2$. More generally, in a real Hilbert space $H$ with $g$ proper, lower semicontinuous and convex and $h$ convex with $L$-Lipschitz gradient, the forward-backward method takes $x_{k+1}$ to be a minimizer of $g(z)+\langle\nabla h(x_k),z-x_k\rangle+\frac1{2\lambda_k}\|z-x_k\|^2$.

A sequence is a subgradient descent sequence for a convex $f$ with constants $a,b>0$ if $x_0\in\operatorname{dom} f$ and, for every $k\ge1$,

• (H1) $f(x_k)+a\|x_k-x_{k-1}\|^2\le f(x_{k-1})$;
• (H2) some $\omega_k\in\partial f(x_k)$ has $\|\omega_k\|\le b\|x_k-x_{k-1}\|$.

With $\min f=0$, $f$ has the KL property on a set $X\cap[0<f<\bar r]$ with desingularizing function $\varphi$ if $\varphi'(f(x))\|v\|\ge1$ for every such $x$ and every $v\in\partial f(x)$.

A number $\nu\ge0$ is a Hoffman constant for a pair of matrices $(A',E)$ if, for every right-hand side with $X=\{A'x\le a'\}$ and $Y=\{Ex=e\}$ intersecting, $\operatorname{dist}(x,X\cap Y)\le\nu\|Ex-e\|$ for all $x\in X$.

Formalization targets

Goal: Theorem 25

Put $a=\frac1{\lambda^+}-\frac L2$ and $b=\frac1{\lambda^-}+L$, and $R=\max\big(\frac{f(x_0)}\mu,1+\frac{\|d\|^2}{2\mu}\big)$. Let $\gamma_R>0$ be such that $f(x)-\min f\ge2\gamma_R\operatorname{dist}^2(x,S)$ whenever $\|x\|_1\le R$. Then ISTA converges to some $x^*\in S$ and

$$f(x_k)-\min f\le\frac{f(x_0)-\min f}{q^k}\quad(k\ge0),\qquad \|x_k-x^*\|\le C\,\frac{\sqrt{f(x_0)-\min f}}{q^{(k-1)/2}}\quad(k\ge1),$$ $$q=1+\frac{2a\gamma_R}{b^2},\qquad C=\frac1{\sqrt a}\Bigg(1+\frac1{ab^{-2}\gamma_R\sqrt{1+\frac1{2ab^{-2}\gamma_R}}}\Bigg).$$

Milestones

1. Lemma 10, (9)–(10): for $R>\|d\|^2/(2\mu)$ the error bound holds with $\gamma_R=1/\big(4\nu^2(1+\mu R+(R\|A\|+\|d\|)(4R\|A\|+\|d\|))\big)$. Here $\nu$ is a Hoffman constant of an explicit pair of matrices built from $A$ and $\mu$.
2. Lemma 10, KL consequence: under the error bound, $f-\min f$ has the KL property on the $\ell_1$ ball with $\varphi(s)=\sqrt{2s/\gamma_R}$.
3. Proposition 13: the forward-backward method satisfies (H1) and (H2) with the constants $a,b$ above.
4. §5.3: every ISTA iterate satisfies $\|x_k\|_1\le R$.
5. Corollary 20 in the stable-set form of Corollary 19: a subgradient descent sequence that stays in a set $X$, on which the KL inequality holds with $\psi(s)=\frac\ell2s^2$, satisfies $f(x_k)\le f(x_0)/(1+2a\sigma)^k$ and the companion bound on $\|x_k-x^*\|$, $\sigma=\ell b^{-2}$.

Significance

The result. The theorem replaces the sublinear $O(1/k)$ guarantee of ISTA by a global linear rate with explicit constants depending only on $\mu$, $A$, $d$, the stepsizes, $x_0$ and a Hoffman constant. It also bounds the distance of the iterates to the limit, not only the objective gap, and needs no uniqueness of the minimizer. Through $\gamma_R$, it shows that conditioning of the objective near its minimizers, rather than strong convexity, governs the complexity of the method. The general machinery behind it (Theorem 16 of the paper, a separate mission of this series) applies to any method producing subgradient descent sequences.

Formalizing it. The result is proved on paper, and no machine-checked proof of it or of its ingredients is known. A formalization would give a verified linear-rate certificate for a widely deployed algorithm. It would also yield reusable lemmas: the descent and relative-error properties of the forward-backward method in a Hilbert space, and quadratic growth for the Lasso on $\ell_1$ balls. Lemma 10 rests on an external result (Beck–Shtern, Lemma 2.5) that a complete development must also formalize.

Difficulty

The Lasso objective is not strongly convex when $A$ has a nontrivial kernel, and its minimizer need not be unique. The textbook argument for linear convergence of proximal gradient methods, which uses strong convexity, therefore does not apply. Quadratic growth $f-\min f\ge2\gamma\operatorname{dist}^2(\cdot,S)$ is the substitute, but it holds only on bounded sets, with a constant degrading in the radius. The analysis must therefore show that the iterates stay in a fixed $\ell_1$ ball and keep track of how the growth constant enters the rate. Deriving the explicit $\gamma_R$ requires recasting the problem as a quadratic program over a polyhedron with $2^n+1$ facets and invoking a Hoffman-type error bound for it.

Formalization scope

Citations refer to the arXiv preprint arXiv:1510.08234v3 (20 July 2016), with its numbering and pages; the journal version is not used. $\mathbb R^n$ is EuclideanSpace ℝ (Fin n) and matrices act through Matrix.toEuclideanLin. Spectral norms are operator norms of the associated continuous linear maps, never Mathlib's default entrywise matrix norm. The data vector, called $b$ in §3.2.1 of the paper and $d$ in §5.3, is $d$ throughout; $b$ is the constant of (H2).

Functions with values in $(-\infty,+\infty]$ are EReal-valued. The convex subdifferential, Γ₀, $\|\cdot\|_1$, the prox predicate and the class $\mathcal K(0,\bar r)$ of desingularizing functions are reused from published platform definitions. The stepsize condition $\lambda^+<2/L$ is written $\lambda^+L<2$, so that $A=0$ (where $L=0$) is covered. $\min f$ is the infimum of $f$, and $q^{(k-1)/2}$ is a real power.

In the goal, $\gamma_R$ is any positive constant satisfying the error bound on the ball of radius $R$; it is not fixed to the value (10) and not existentially quantified. The constants $a,b,q,C,R$ are fixed formulas, never existentials. A statement with an unspecified constant, with $\gamma_R$ chosen after the sequence, or with the stepsize hypothesis written as $\lambda^+<2/L$ (unsatisfiable or vacuous at $L=0$) would trivialize the goal and is ruled out.

The subgradient-descent and KL notions duplicate, in the sub-namespace ErrBoundCplx.ISTA, those of the series' first mission, because draft items cannot import each other. Proofs of any milestone are welcome, in particular Proposition 13 and the KL consequence of Lemma 10, which are self-contained. A formalization of Hoffman's error bound and of Beck–Shtern's Lemma 2.5 is the main infrastructure Lemma 10 needs.

Selected references

• J. Bolte, T. P. Nguyen, J. Peypouquet, B. W. Suter, From error bounds to the complexity of first-order descent methods for convex functions, Math. Program. 165 (2017) 471–507; arXiv:1510.08234v3. https://arxiv.org/abs/1510.08234v3
• A. Beck, S. Shtern, Linearly convergent away-step conditional gradient for non-strongly convex functions, preprint (cited as [10], Lemma 2.5). https://arxiv.org/abs/1504.05002
• A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imaging Sci. 2 (2009) 183–202. https://doi.org/10.1137/080716542
• A. J. Hoffman, On approximate solutions of systems of linear inequalities, J. Res. Natl. Bur. Stand. 49 (1952) 263–265. https://doi.org/10.6028/jres.049.027
• I. Daubechies, M. Defrise, C. De Mol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Comm. Pure Appl. Math. 57 (2004) 1413–1457. https://doi.org/10.1002/cpa.20042
• J. Liang, J. Fadili, G. Peyré, Local linear convergence of forward-backward under partial smoothness, preprint. https://arxiv.org/abs/1407.5611
• P. L. Combettes, V. R. Wajs, Signal recovery by proximal forward-backward splitting, Multiscale Model. Simul. 4 (2005) 1168–1200. https://doi.org/10.1137/050626090
• H. Attouch, J. Bolte, B. F. Svaiter, Convergence of descent methods for semi-algebraic and tame problems, Math. Program. 137 (2013) 91–129. https://doi.org/10.1007/s10107-011-0484-9