All missions

Motivation

For a real symmetric matrix, the largest eigenvalue is the maximum of the quadratic form $\mathbf x^\top A\mathbf x$ on the unit sphere, and it is computable to any accuracy in polynomial time. Tensor eigenvalues, introduced independently by Lim (2005) and Qi (2005), extend this picture to higher order: they are the stationary values of the cubic form $\mathcal A(\mathbf x,\mathbf x,\mathbf x)$ on the unit sphere. They appear in multilinear statistics, signal processing, hypergraph theory and polynomial optimization, and for these applications the symmetric case is the natural one.

Hillar and Lim, in Most Tensor Problems Are NP-Hard (J. ACM 2013), show that the eigenvalue problem stays NP-hard when the tensor is required to be symmetric (their Theorem 9.3). The proof is short once one ingredient is available: a theorem of Nesterov (2003), in a form corrected by Hillar and Lim, that expresses the stability number of a graph as the maximum of an explicit cubic polynomial over a sphere. This mission formalizes that chain.

Setting

Let $G=(V,E)$ be a simple graph on $V=\{1,\dots,v\}$, $v\ge1$. A set of vertices is stable if no two of its members are adjacent; the stability number $\alpha(G)$ is the size of a largest stable set, so $1\le\alpha(G)\le v$.

Put $n=v+v(v-1)/2$ and write a vector of $\mathbb R^n$ as $\mathbf z=(\mathbf x,\mathbf y)$, with $\mathbf x=(x_1,\dots,x_v)$ indexed by vertices and $\mathbf y=(y_{ij})_{i<j}$ indexed by pairs of vertices. The tensor $\mathcal S_G=[\![s_{abc}]\!]\in\mathbb R^{n\times n\times n}$ has entry $1$ at the six orderings of the triple of coordinates $(x_i,x_j,y_{ij})$ for every non-adjacent pair $i<j$, and $0$ elsewhere. It is symmetric: its entries do not change under any permutation of the three indices.

For a real 3-tensor $\mathcal A=[\![a_{abc}]\!]$ the cubic form is $\mathcal A(\mathbf z,\mathbf z,\mathbf z)=\sum_{a,b,c}a_{abc}z_az_bz_c$, and a real $\lambda$ is a unit $\ell^2$-eigenvalue of $\mathcal A$ if

$$\sum_{a,b}a_{abc}z_az_b=\lambda z_c\quad\text{for all }c,\qquad\text{for some }\mathbf z\text{ with }\|\mathbf z\|_2=1 .$$

Finally, $\lambda_l=2\sqrt{\tfrac23\left(1-\tfrac1l\right)}$ for $l\ge1$. In Lean these are TensorNP.SymEigen.stabTensor, cubicForm, IsUnitL2Eigenvalue and lambdaL, and $\alpha(G)$ is Mathlib's SimpleGraph.indepNum.

Formalization targets

Goal: Theorem 9.3, core

$$\alpha(G)=\max\{\,l\in\{1,\dots,v\}:\ \lambda_l\text{ is a unit }\ell^2\text{-eigenvalue of }\mathcal S_G\,\}.$$

That is, $\lambda_{\alpha(G)}$ is an eigenvalue of $\mathcal S_G$, and no $\lambda_l$ with $\alpha(G)<l\le v$ is. This is the statement on which the paper's $v$-query reduction rests; nothing is claimed about $l<\alpha(G)$.

Milestones, in the paper's order

1. Theorem 9.2 (Nesterov). With $\mathbb S^{n-1}=\{(\mathbf x,\mathbf y):\|\mathbf x\|_2^2+\|\mathbf y\|_2^2=1\}$,

$$\sqrt{1-\frac1{\alpha(G)}}=3\sqrt{\frac32}\cdot\max_{(\mathbf x,\mathbf y)\in\mathbb S^{n-1}}\ \sum_{i<j,\ \{i,j\}\notin E}x_ix_jy_{ij}.$$

2. The tensor. $\mathcal S_G$ is symmetric and $\mathcal S_G(\mathbf z,\mathbf z,\mathbf z)=6\sum_{i<j,\{i,j\}\notin E}x_ix_jy_{ij}$.
3. The maximum is an eigenvalue. $\max_{\|\mathbf z\|_2=1}\mathcal S_G(\mathbf z,\mathbf z,\mathbf z)$ is attained and is a unit $\ell^2$-eigenvalue of $\mathcal S_G$.
4. Its value. $\max_{\|\mathbf z\|_2=1}\mathcal S_G(\mathbf z,\mathbf z,\mathbf z)=\lambda_{\alpha(G)}$.
5. Further result: Theorem 6.9 (Banach), (23). For every symmetric $\mathcal S\in\mathbb R^{n\times n\times n}$,

$$\sup_{\mathbf x,\mathbf y,\mathbf z\ne\mathbf 0}\frac{|\mathcal S(\mathbf x,\mathbf y,\mathbf z)|}{\|\mathbf x\|_2\|\mathbf y\|_2\|\mathbf z\|_2}=\sup_{\mathbf x\ne\mathbf 0}\frac{|\mathcal S(\mathbf x,\mathbf x,\mathbf x)|}{\|\mathbf x\|_2^3},$$

both sides finite. The paper uses it in §10 to transfer Theorem 9.3 to the spectral norm and the best rank-one approximation of symmetric tensors.

Significance

The result. Theorem 9.3 places symmetric tensor eigenvalue among the NP-hard problems of Hillar and Lim's Table I, and through Banach's theorem (Theorem 10.2 of the paper) the same holds for the largest singular value, the spectral norm and the best rank-one approximation of a symmetric 3-tensor. These are the quantities computed by higher-order power methods and by moment-based estimators, so the result marks where the matrix theory stops transferring.

Formalizing it. All results here are proved in the literature; none is machine-checked. Nesterov's theorem, which rests on the Motzkin–Straus theorem for the complement graph, is not on the platform or in Mathlib, and its formal statement settles the constant that the original paper got wrong by a factor $1/\sqrt2$ (Hillar–Lim, footnote 12). Banach's theorem on symmetric multilinear forms is likewise absent from Mathlib and is useful well beyond this mission.

Difficulty

The tensor identity (milestone 2) is bookkeeping. The analytic content sits in two places. First, Nesterov's theorem: the obvious route, fixing $\mathbf x$ and maximizing over $\mathbf y$ by Cauchy–Schwarz, reduces the cubic problem to a quartic one in $\mathbf x$, whose value still has to be identified with $\alpha(G)$; that identification is the Motzkin–Straus theorem, which is not available in Lean and is itself a nontrivial optimization over the simplex. Second, the passage from a maximum to an eigenvalue needs a Lagrange multiplier argument on the sphere, and the symmetry of $\mathcal S_G$ is what makes the gradient of the cubic form match the eigen-equation. Part (b) of the goal needs a further, unprinted observation relating every unit eigenpair to the cubic form. Banach's theorem has no short proof; the classical arguments go through polarization estimates or a variational argument specific to symmetric forms.

Formalization scope

• Tensors and indices. A 3-tensor is a function ι → ι → ι → ℝ on a finite index type. For $\mathcal S_G$ the index type is the disjoint union of the vertices Fin v and the pairs $\{(i,j): i<j\}$; this is $\mathbb R^n$ with coordinates labelled by name instead of by the paper's enumeration $\varphi(i,j)=(i-1)v-i(i-1)/2+j-i$. Vertices are $0,\dots,v-1$.
• Normalization. Eigenvalues are taken with a unit eigenvector. The paper's Problem 9.1 and Definition 5.1 say only $\mathbf x\ne\mathbf 0$; since $(\lambda,\mathbf x)\mapsto(t\lambda,t\mathbf x)$ preserves the eigen-equation, that reading would make every nonzero real an eigenvalue of $\mathcal S_G$ and the goal false. The unit sphere (13) is where the paper derives the notion.
• Explicit readings. The paper's index range "$1\le i<j<k\le v$" in the definition of $s_{ijk}$ is read as $k\le n$ (taken literally it gives $\mathcal S_G=0$). Every "max" is stated as a greatest element, so attainment is claimed. In (23) the suprema are real suprema over nonzero vectors with boundedness part of the conclusion. $v\ge1$ is assumed wherever the paper uses $\alpha(G)\in\{1,\dots,v\}$.
• Not formalized. The oracle procedure of the proof of Theorem 9.3, the polynomial size of $\mathcal S_G$, the input model of Problem 9.1 with quadratic irrationalities (Remark 9.4), and the NP-completeness of the stability number (Karp 1972, via $\alpha(G)=\omega(\overline G)$). The goal is the mathematical core of the Turing reduction, not a complexity statement.
• Trivializations ruled out. Without the unit normalization the eigenvalue statements are trivial or false; a goal stating only (a), or only an inequality between the maximum and $\lambda_{\alpha(G)}$, would not determine $\alpha(G)$; both are excluded by the statements above.
• Infrastructure. Everything rests on Mathlib (SimpleGraph.indepNum, real analysis on finite-dimensional spaces, Lagrange multipliers). A Motzkin–Straus theorem for SimpleGraph and Banach's theorem for symmetric trilinear forms are both reusable beyond this mission; contributions of either, or of a general "maximum of a symmetric cubic form on the sphere is an eigenvalue" lemma, are welcome.

Selected references

• C. J. Hillar, L.-H. Lim, Most Tensor Problems Are NP-Hard, Journal of the ACM 60(6), Article 45, 2013. https://doi.org/10.1145/2512329 (arXiv:0911.1393v5: https://arxiv.org/abs/0911.1393)
• Yu. Nesterov, Random walk in a simplex and quadratic optimization over convex polytopes, CORE Discussion Paper 2003/71, 2003. http://edoc.bib.ucl.ac.be:83/archive/00000238/01/dp2003-71.pdf
• E. de Klerk, The complexity of optimizing over a simplex, hypercube or sphere: a short survey, Central European Journal of Operations Research 16, 111–125, 2008. https://doi.org/10.1007/s10100-007-0052-9
• T. S. Motzkin, E. G. Straus, Maxima for graphs and a new proof of a theorem of Turán, Canadian Journal of Mathematics 17, 533–540, 1965. https://doi.org/10.4153/CJM-1965-053-6
• S. Banach, Über homogene Polynome in ($L^2$), Studia Mathematica 7(1), 36–44, 1938 (reference [Banach 1938] of Hillar–Lim).
• L.-H. Lim, Singular values and eigenvalues of tensors: a variational approach, Proc. IEEE CAMSAP 1, 129–132, 2005. https://arxiv.org/abs/math/0607648
• L. Qi, Eigenvalues of a real supersymmetric tensor, Journal of Symbolic Computation 40, 1302–1324, 2005. https://doi.org/10.1016/j.jsc.2005.05.007
• R. M. Karp, Reducibility among combinatorial problems, in Complexity of Computer Computations, 85–103, 1972. https://doi.org/10.1007/978-1-4684-2001-2_9

Motivation

Interior-point methods for conic optimization are best understood on symmetric cones: the nonnegative orthant, the second-order cone and the cone of positive semidefinite matrices. On these cones the barrier functions are self-scaled, and Nesterov and Todd showed that every primal-dual pair $(x, s)$ then has a unique scaling point $w$ with $s = F''(w)x$ and $F'(x) = F''(w)F_*'(s)$. The resulting primal-dual algorithms, implemented in solvers such as SeDuMi and MOSEK, combine long steps with good worst-case complexity.

Many models in statistics, geometric programming, entropy maximization and logistic regression need the exponential cone, which is not symmetric. Dahl and Andersen (Math. Program. 194, 2022) designed the primal-dual algorithm behind MOSEK's exponential-cone support. To explain why they cannot simply transplant the Nesterov–Todd construction, they use a property that sits between "self-scaled" and "arbitrary": negative curvature of the barrier. Barriers of symmetric cones have it; a barrier is self-scaled if and only if both it and its conjugate have it (Nesterov–Tunçel, 2016); and barriers with negative curvature still admit a unique scaling point satisfying one secant equation. Dahl and Andersen observe in §5 that the exponential-cone barrier lies outside this class, by an explicit witness. This mission formalizes that observation and the derivative formulas of their Appendix A it rests on.

Timeline:

• 1997–1998: Nesterov and Todd introduce self-scaled barriers and their primal-dual scalings (Math. Oper. Res. 22; SIAM J. Optim. 8).
• 2001: Tunçel generalizes primal-dual scalings to arbitrary convex cones (Found. Comput. Math. 1).
• 2009: Chares studies the exponential cone and its 3-self-concordant barrier (thesis, Université catholique de Louvain).
• 2016: Nesterov and Tunçel characterize self-scaled barriers through negative curvature of the barrier and its conjugate.
• 2022: Dahl and Andersen give a primal-dual algorithm for the exponential cone with Tunçel-type scalings, and note that its barrier does not have negative curvature.

Setting

Points of $\mathbb{R}^3$ are written $x = (x_1, x_2, x_3)$. The exponential cone is the closure

$$K_{\exp} = \operatorname{cl}\{x \in \mathbb{R}^3 \mid x_1 \ge x_2 \exp(x_3/x_2),\ x_2 > 0\},$$

whose interior is $\{x \mid x_2 > 0,\ x_1 > x_2\exp(x_3/x_2)\}$. Its standard barrier is

$$F(x) = -\log\bigl(x_2\log(x_1/x_2) - x_3\bigr) - \log x_1 - \log x_2, \qquad (2)$$

defined and smooth on $\operatorname{int}(K_{\exp})$. Appendix A of the paper writes $F = g + h$ with

$$\psi(x) = x_2\log(x_1/x_2) - x_3,\qquad g(x) = -\log\psi(x),\qquad h(x) = -\log x_1 - \log x_2 .$$

For a smooth $F$ on an open set, $F''(x)$ is its Hessian, a symmetric linear map of $\mathbb{R}^n$, and the third directional derivative is the symmetric linear map

$$F'''(x)[u] = \frac{d}{dt}F''(x + tu)\Big|_{t=0}, \qquad F'''(x)[u, v, v] = \langle F'''(x)[u]\,v, v\rangle .$$

A barrier $F$ of a cone $K$ has negative curvature if for all $x \in \operatorname{int}(K)$ and for $u \in K$

$$F'''(x)[u] \preceq 0,$$

where $\preceq$ is the Loewner order: $\langle F'''(x)[u]\,v, v\rangle \le 0$ for every $v$. The direction $u$ ranges over the closed cone $K$, including its boundary.

Formalization targets

Goal: the exponential-cone barrier does not have negative curvature

$$\neg\Bigl(\forall x \in \operatorname{int}(K_{\exp}),\ \forall u \in K_{\exp},\ \forall v \in \mathbb{R}^3:\ \langle F'''(x)[u]\,v, v\rangle \le 0\Bigr).$$

The goal asserts only the qualitative failure, so it does not depend on the particular witness.

Milestones

1. The witness: $\hat x := (1, e^{-2}, 0) \in \operatorname{int}(K_{\exp})$ and $\hat u := (1, 0, 0) \in K_{\exp}\setminus\operatorname{int}(K_{\exp})$ (§5, p. 356).
2. First-order derivatives (A.1, p. 367): $\psi'(x) = (x_2/x_1,\ \log(x_1/x_2) - 1,\ -1)$, $h'(x) = -(1/x_1, 1/x_2, 0)$, $g'(x) = -\psi'(x)/\psi(x)$, $F' = g' + h'$ on $\operatorname{int}(K_{\exp})$.
3. Second-order derivatives (A.2, p. 367):

$$F''(x) = -\frac{1}{\psi(x)}\Bigl(\psi''(x) - \frac{\psi'(x)\psi'(x)^T}{\psi(x)}\Bigr) + h''(x).$$

4. Third-order directional derivatives (A.3 and (33), p. 368): the closed form of $F'''(x)[u]$ in terms of $\psi$, its derivatives and $h'''(x)[u]$.
5. The matrix at the witness (§5, p. 356):

$$F'''(\hat x)[\hat u] = \begin{bmatrix} -4 & e^2/2 & e^2/2 \\ e^2/2 & 0 & 0 \\ e^2/2 & 0 & -e^4/4 \end{bmatrix}.$$

6. The positive value (§5, p. 356): $F'''(\hat x)[\hat u, \hat v, \hat v] = 4$ for $\hat v := (1, 8e^{-2}, 4e^{-2})$.

Significance

The result. Negative curvature is the property that gives symmetric-cone barriers their long-step Hessian estimates and a unique scaling point with one secant equation. Its failure for the exponential-cone barrier means that none of these consequences can be invoked for the exponential cone. An algorithm for this cone has to obtain its scalings by other means; Dahl and Andersen take them from Tunçel's framework, via the quasi-Newton secant updates of §5. The witness also shows that the failure occurs already for a direction on the boundary of the cone, at an interior point with simple coordinates.

Formalizing it. The claim is proved in the paper by a computation "using the expressions in the appendix", with the details left to the reader. The appendix formulas for the gradient, Hessian and third directional derivative of the exponential-cone barrier are the ones an implementation evaluates, including the higher-order corrector term $-\tfrac12F'''(x)[u, v]$. A machine-checked version certifies these closed forms, including the zero padding of the $2\times2$ "leading parts" in (33), and gives a reusable definition of negative curvature that applies to any cone and barrier. As far as is known, none of these statements has a machine-checked proof.

Difficulty

The obstacle is computational, not conceptual. The barrier is a composition of logarithms with a function $\psi$ that is itself a perspective of the logarithm, so its third derivative is a sum of terms with denominators $\psi(x)^k$, $x_1^k$ and $x_2^k$. Stating it as a derivative of the Hessian map requires showing that the Hessian is differentiable on the open interior, which in turn requires a description of that interior, a set defined as the interior of a closure. The membership claims about $\hat u$ are not definitional either: $\hat u$ lies on the boundary, reached only as a limit of points of the generating set, and its non-interiority needs a neighbourhood argument. Finally, the cone condition on $u$ matters: replacing $u \in K_{\exp}$ by $u \in \mathbb{R}^3$ makes the claim nearly trivial, because $F'''(x)[-u] = -F'''(x)[u]$.

Formalization scope

• Points of $\mathbb{R}^3$ are EuclideanSpace ℝ (Fin 3); the paper's $(x_1, x_2, x_3)$ are (x 0, x 1, x 2), and vectors are written !₂[·, ·, ·].
• $K_{\exp}$ is the closure of the printed set, not its interior or an equivalent description. $F$, $\psi$, $g$, $h$ are defined on all of $\mathbb{R}^3$ by their formulas with Real.log; every statement evaluates them, and their derivatives, only at points of $\operatorname{int}(K_{\exp})$, where all logarithms are the genuine ones. The appendix formulas carry the hypothesis $x \in \operatorname{int}(K_{\exp})$.
• $F''$ is the published SelfScaledIPM.ShortStep.hess (derivative of the gradient), and $F'''(x)[u]$ is fderiv ℝ (hess F) x u. The Loewner order is stated through the quadratic form for every $v$. Matrices act through the standard basis (Matrix.toEuclideanCLM).
• Disclosed deviations from the page: milestone 1 asserts $\hat x \in \operatorname{int}(K_{\exp})$, stronger than the printed $\hat x \in K_{\exp}$ and the form the definition needs. In A.3 the $2\times2$ leading parts $\hat h'''$, $\hat\psi'''$ of (33) are padded with a zero third row and column, and $g''(x)$ is written out by the A.2 formula. Only the first display of A.2 is stated, not the factorization $F''(x) = R(x)R(x)^T$.
• A trivializing formalization would let $u$ range over all of $\mathbb{R}^3$, or read $\preceq 0$ entrywise; both are excluded, and $u$ ranges over the closed cone exactly as in the paper.
• The claim that $F$ is a 3-self-concordant barrier (quoted from Chares) is not part of the mission. The definition of negative curvature applies to any cone in any dimension and can be reused for other barriers. Contributions of calculus lemmas for logarithmic barriers on open sets (differentiability of the Hessian map, interiors of closures of epigraph-type sets) are welcome.

Selected references

• J. Dahl, E. D. Andersen, A primal-dual interior-point algorithm for nonsymmetric exponential-cone optimization, Mathematical Programming 194 (2022), 341–370. https://doi.org/10.1007/s10107-021-01631-4
• Yu. Nesterov, M. J. Todd, Self-scaled barriers and interior-point methods for convex programming, Mathematics of Operations Research 22 (1997), 1–42. https://doi.org/10.1287/moor.22.1.1
• Yu. Nesterov, M. J. Todd, Primal-dual interior-point methods for self-scaled cones, SIAM Journal on Optimization 8 (1998), 324–364. https://doi.org/10.1137/S1052623495290209
• L. Tunçel, Generalization of primal-dual interior-point methods to convex optimization problems in conic form, Foundations of Computational Mathematics 1 (2001), 229–254. https://doi.org/10.1007/s102080010008
• Yu. Nesterov, L. Tunçel, Local superlinear convergence of polynomial-time interior-point methods for hyperbolicity cone optimization problems, SIAM Journal on Optimization 26 (2016), 139–170. https://doi.org/10.1137/140978296
• P. R. Chares, Cones and interior-point algorithms for structured convex optimization involving powers and exponentials, PhD thesis, Université catholique de Louvain, 2009.

Motivation

A matrix Gaussian series $\sum_k \gamma_k A_k$, a sum of fixed self-adjoint matrices each multiplied by an independent standard normal variable, is the simplest sum of independent random matrices. It already appears in the analysis of random designs, in compressed sensing, in randomized numerical linear algebra and in quantum information, wherever one needs to know how large the spectrum of a random linear combination of matrices can be. Its Rademacher analogue $\sum_k \varepsilon_k A_k$, with random signs, plays the same role in symmetrization arguments.

For real coefficients $a_k$ the answer is classical: $\sum_k \gamma_k a_k$ is normal with variance $\sum_k a_k^2$, so its upper tail is at most $e^{-t^2/2\sigma^2}$. Tropp's paper (arXiv:1004.4389v7) shows that the same bound holds for the largest eigenvalue of a matrix series, with an explicit dimensional factor and with the variance parameter $\|\sum_k A_k^2\|$.

Timeline. Noncommutative Khintchine inequalities (Lust-Piquard 1986; Lust-Piquard and Pisier 1991) control the moments of matrix Rademacher and Gaussian series and imply tail bounds of this type. Ahlswede and Winter (2002, arXiv:quant-ph/0012127) introduced the matrix Laplace transform method, whose bounds depend on the sum of the norms $\sum_k\|A_k^2\|$. Oliveira (2010, arXiv:1004.3821) obtained bounds of the same form as Theorem 4.1. Tropp (2012) derived the result from a master inequality built on Lieb's concavity theorem (1973), which produces the norm of the sum, $\|\sum_k A_k^2\|$, as variance parameter. The paper itself notes that Theorem 4.1 and Corollary 4.2 are not new results; the contribution is the route through the master bound.

Setting

Fix a dimension $d\ge1$. Matrices are $d\times d$ complex arrays; a matrix is self-adjoint if it equals its conjugate transpose. For a self-adjoint $A$, $\lambda_{\max}(A)$ is its largest eigenvalue, and $\|A\|$ is the spectral norm, the operator norm on $\mathbb C^d$ with the Euclidean norm. Functions of a self-adjoint matrix are defined through its eigen-decomposition: $e^A$, $\log A$ (for positive-definite $A$), $\cosh A$. The semidefinite order $A\preccurlyeq H$ means that $H-A$ is positive semidefinite.

A random matrix is a measurable map from a probability space $(\Omega,\mathcal F,\mathbb P)$ into the $d\times d$ complex matrices; its expectation $\mathbb E X$ is taken entry by entry. A Rademacher variable takes the values $\pm1$ with probability $1/2$ each; a standard normal variable has law $N(0,1)$.

The data of the mission are fixed self-adjoint matrices $A_1,\dots,A_n$ and independent scalar variables $\xi_1,\dots,\xi_n$, all standard normal or all Rademacher. The variance parameter is

$$\sigma^2 := \Big\|\sum_k A_k^2\Big\|.$$

Formalization targets

Goal: Theorem 4.1 (p. 14)

For all $t\ge0$,

$$\mathbb P\Big\{\lambda_{\max}\Big(\sum_k\xi_kA_k\Big)\ge t\Big\}\le d\,e^{-t^2/2\sigma^2},\qquad \mathbb P\Big\{\Big\|\sum_k\xi_kA_k\Big\|\ge t\Big\}\le 2d\,e^{-t^2/2\sigma^2}.$$

Milestones

1. Theorem 3.6 (master tail bound, p. 12): for independent random self-adjoint $X_k$ and every $t$,

$$\mathbb P\Big\{\lambda_{\max}\Big(\sum_kX_k\Big)\ge t\Big\}\le\inf_{\theta>0}e^{-\theta t}\operatorname{tr}\exp\Big(\sum_k\log\mathbb Ee^{\theta X_k}\Big).$$

2. Corollary 3.7 (p. 12): if $\mathbb Ee^{\theta X_k}\preccurlyeq e^{g(\theta)A_k}$ with $g\ge0$, then the tail is at most $d\inf_{\theta>0}e^{-\theta t+g(\theta)\rho}$ with $\rho=\lambda_{\max}(\sum_kA_k)$.
3. Display (2.4) (p. 8): $\cosh(A)\preccurlyeq e^{A^2/2}$ for every self-adjoint $A$.
4. Lemma 4.3 (p. 15): $\mathbb Ee^{\varepsilon\theta A}\preccurlyeq e^{\theta^2A^2/2}$ for Rademacher $\varepsilon$ and $\mathbb Ee^{\gamma\theta A}=e^{\theta^2A^2/2}$ for standard normal $\gamma$.

Further item: Corollary 4.2 (p. 15)

For fixed $d_1\times d_2$ matrices $B_k$ and $\sigma^2=\max\{\|\sum_kB_kB_k^*\|,\|\sum_kB_k^*B_k\|\}$,

$$\mathbb P\Big\{\Big\|\sum_k\xi_kB_k\Big\|\ge t\Big\}\le(d_1+d_2)\,e^{-t^2/2\sigma^2}.$$

Significance

The result. Theorem 4.1 is the template for every inequality in the paper: a scalar tail bound carried over to matrices at the cost of a factor $d$, with a variance parameter that is the norm of a sum rather than a sum of norms. The difference matters: $\sum_k\|A_k^2\|$ can exceed $\|\sum_kA_k^2\|$ by a factor of $d$, and §4 of the paper shows that both the variance parameter and the dimensional factor in (4.3) are needed in general. Corollary 4.2 transfers the bound to the largest singular value of a rectangular series, which covers Gaussian matrices with nonuniform variances (§4.3).

Formalizing it. The theorem is proved in the paper; the work here is to formalize that proof. The milestones are reusable beyond this mission: Theorem 3.6 and Corollary 3.7 are the common entry point of the Chernoff, Bernstein and Azuma missions of this series, and the Rademacher half of Lemma 4.3 is used again in the matrix Azuma mission. No machine-checked proof of these matrix inequalities was found in Mathlib or on the platform when this mission was drafted.

Difficulty

The scalar argument bounds $\mathbb E e^{\theta\sum_k X_k}$ by the product of the individual mgfs, using $e^{a+b}=e^ae^b$. For matrices this identity fails when the summands do not commute, and the trace inequality that replaces it for two matrices (Golden–Thompson) has no analogue for three or more. Controlling the trace mgf of a sum of many noncommuting independent matrices therefore requires a genuinely matrix-analytic input with no scalar counterpart; the naive product bound is not available. A second difficulty is infrastructural: matrix functions, the semidefinite order and expectations of random matrices must be combined in a single setting.

Formalization scope

Matrices are Matrix (Fin d) (Fin d) ℂ with [NeZero d]; self-adjointness is IsHermitian; $e^A$, $\log A$, $\cosh A$ are Mathlib's continuous functional calculus cfc; $\preccurlyeq$ is Mathlib's Loewner order under MatrixOrder; $\lambda_{\max}$ is the supremum of the real spectrum; $\|\cdot\|$ is the operator norm on $\ell_2$. Expectation is entrywise, and the paper's standing regularity assumption is made explicit as entrywise integrability of every matrix whose expectation is taken. The infimum over $\theta>0$ is stated pointwise: the bound is asserted for every $\theta>0$ at which the mgfs exist, which is equivalent to the infimum form. Independence is iIndepFun of the scalar variables, and "Gaussian or Rademacher" is the disjunction of two hypotheses, never a mixture of laws.

The variance parameter is the norm of the sum $\|\sum_kA_k^2\|$; a formalization with $\sum_k\|A_k^2\|$, with a constant different from $1/2$, or with the dimensional factor omitted is a different theorem. At $\sigma^2=0$, Lean's division by zero makes the bound the trivial $d$; this is the only degenerate case and it is valid.

The development needs: the spectral calculus for Hermitian matrices, monotonicity of the trace exponential, operator monotonicity of the matrix logarithm, Lieb's theorem, and moments of the Gaussian. These are reusable well beyond this paper; proofs of the milestones, of supporting matrix-analysis lemmas, and alternative routes (for instance through noncommutative Khintchine inequalities) are all welcome.

Selected references

• J. A. Tropp, User-Friendly Tail Bounds for Sums of Random Matrices, Found. Comput. Math. 12 (2012) 389–434; cited from arXiv:1004.4389v7. https://arxiv.org/abs/1004.4389, https://doi.org/10.1007/s10208-011-9099-z
• R. Ahlswede and A. Winter, Strong converse for identification via quantum channels, IEEE Trans. Inform. Theory 48 (2002) 569–579. https://arxiv.org/abs/quant-ph/0012127
• R. I. Oliveira, Sums of random Hermitian matrices and an inequality by Rudelson, Electron. Commun. Probab. 15 (2010) 203–212. https://arxiv.org/abs/1004.3821
• E. H. Lieb, Convex trace functions and the Wigner–Yanase–Dyson conjecture, Adv. Math. 11 (1973) 267–288. https://doi.org/10.1016/0001-8708(73)90011-X
• F. Lust-Piquard and G. Pisier, Non commutative Khintchine and Paley inequalities, Ark. Mat. 29 (1991) 241–260. https://doi.org/10.1007/BF02384340

Motivation

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

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

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

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

Setting

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

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

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

The problem is

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

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

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

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

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

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

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

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

Formalization targets

Goal: Theorem 3.2

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

Milestones, in the order the proof uses them

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

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

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

Significance

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

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

Difficulty

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

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

Formalization scope

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

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

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

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

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

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

Selected references

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

Motivation

Coordinate descent methods update one block of variables at a time. They are the method of choice when the number of variables is so large that even one full gradient is too expensive, while a single partial derivative is cheap: support vector machines, $\ell_1$-regularized regression, and large sparse least-squares problems are the standard examples. Yu. Nesterov's paper Efficiency of coordinate descent methods on huge-scale optimization problems (CORE Discussion Paper 2010/2; journal version SIAM J. Optim. 22 (2012), doi:10.1137/100802001) gave the first global efficiency estimates for randomized coordinate descent, and these estimates started a large literature on random coordinate and block methods.

Many applications carry simple constraints on each block: bounds on each variable, a simplex or a ball per block. Section 4 of the paper shows that randomized coordinate descent handles such block-separable constraints at no loss: a projected block step, chosen uniformly at random, enjoys the same $O(n/k)$ expected rate as the unconstrained method, and a linear rate under strong convexity. This mission formalizes that result, Theorem 5.

Setting

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

The objective $f:\mathbb R^N\to\mathbb R$ is convex and differentiable, and its gradient is coordinate-wise Lipschitz (2.2): writing $f'_i(x)=U_i^T\nabla f(x)$ for the $i$-th block of the gradient, there are constants $L_i>0$ with

$$\|f'_i(x+U_ih)-f'_i(x)\|^*_{(i)}\le L_i\|h\|_{(i)}\quad\text{for all }x,\ i,\ h .$$

The weighted norm is $\|x\|_1^2=\sum_iL_i\|x^{(i)}\|_{(i)}^2$ (3.5). The function $f$ is strongly convex in $\|\cdot\|_1$ with constant $\sigma>0$ if $f(y)\ge f(x)+\langle\nabla f(x),y-x\rangle+\frac\sigma2\|y-x\|_1^2$ for all $x,y$ (3.1).

The problem is (4.1), $\min_{x\in Q}f(x)$ with $Q=Q_1\times\cdots\times Q_n$ and each $Q_i\subseteq\mathbb R^{n_i}$ nonempty, closed and convex. Its optimal value is $f^*$, attained at an optimal solution $x_*\in Q$. The constrained coordinate update (4.2) is

$$u^{(i)}(x)=\arg\min_{u\in Q_i}\Big[\langle f'_i(x),u-x^{(i)}\rangle+\tfrac{L_i}2\|u-x^{(i)}\|_{(i)}^2\Big],\qquad V_i(x)=x+U_i\big(u^{(i)}(x)-x^{(i)}\big).$$

The uniform coordinate descent method UCDM$(x_0)$ (4.5) draws $i_k$ uniformly from $\{1,\dots,n\}$ and sets $x_{k+1}=V_{i_k}(x_k)$. Its expected objective is $\varphi_k=\mathbb E f(x_k)$, the expectation over $i_0,\dots,i_{k-1}$. The level-set radius $R_1(x_0)$ is the largest $\|\cdot\|_1$-distance between a feasible $x$ with $f(x)\le f(x_0)$ and an optimal solution.

Formalization targets

Goal: Theorem 5

For every $k\ge0$,

$$\varphi_k-f^*\le\frac n{n+k}\Big[\tfrac12R_1^2(x_0)+f(x_0)-f^*\Big],$$

and, if $f$ is strongly convex in $\|\cdot\|_1$ with constant $\sigma>0$,

$$\varphi_k-f^*\le\Big(1-\frac{2\sigma}{n(1+\sigma)}\Big)^k\Big(\tfrac12R_1^2(x_0)+f(x_0)-f^*\Big).\qquad(4.6)$$

Milestones

1. (4.3): the optimality condition $\langle f'_i(x)+L_iB_i(u^{(i)}(x)-x^{(i)}),u-u^{(i)}(x)\rangle\ge0$ for all $u\in Q_i$.
2. (4.4): for feasible $x$, $f(x)-f(V_i(x))\ge\frac{L_i}2\|u^{(i)}(x)-x^{(i)}\|_{(i)}^2$.
3. (4.7): with $r_k=\|x_k-x_*\|_1$, $\ \mathbb E_{i_k}\big(\frac12r_{k+1}^2+f(x_{k+1})-f^*\big)\le\frac12r_k^2+f(x_k)-f^*+\frac1n\langle\nabla f(x_k),x_*-x_k\rangle$.
4. Footnote 2: under (2.2), the strong convexity constant satisfies $\sigma\le1$.

Significance

The result. Theorem 5 shows that the cost of a constrained problem with separable constraints is, up to the projection onto each $Q_i$, the cost of the unconstrained one: $O(n/\epsilon)$ block steps for accuracy $\epsilon$, and $O\big(\frac n\sigma\log\frac1\epsilon\big)$ under strong convexity. The method needs no knowledge of $\sigma$; the same iteration attains whichever rate applies. The Lyapunov function $\frac12\|x-x_*\|_1^2+f(x)-f^*$ of (4.7) needs no level-set argument, so it applies wherever the projected step is cheap.

Formalizing it. The theorem is proved in the paper; no machine-checked proof of it, or of any rate for a randomized constrained coordinate method, appears in the Prove2Me catalogue or in Mathlib. The Prove2Me catalogue holds Bubeck's scalar, unconstrained random coordinate descent results (ConvexOptAlg.CoordDescent.*) and the deterministic Tseng–Yun coordinate gradient descent; neither covers a projected block step. A formalization produces the first verified rate for a randomized method with constraints, and a reusable model of block-structured spaces, block projections and expectations over random index sequences.

Difficulty

The unconstrained analysis of the paper (Theorem 1) bounds the distance to the optimum by a level-set radius and uses the decrease $f(x)-f(T_i(x))\ge(\|f'_i(x)\|^*_{(i)})^2/(2L_i)$ of a gradient step. Neither carries over: with constraints the gradient at the optimum need not vanish, so the decrease of a projected step is not controlled by the gradient norm, and the unconstrained recursion $\varphi_k-\varphi_{k+1}\ge(\varphi_k-f^*)^2/C$ is not available. A different potential is needed, and its one-step estimate (4.7) is the central milestone. For the linear rate, the contraction factor must be shown nonnegative without assuming $\sigma\le1$, which is where footnote 2 enters.

Formalization scope

$\mathbb R^N$ is the dependent product Blocks E $=\prod_{i:\mathrm{Fin}\,n}E_i$, each $E_i$ a finite-dimensional real inner product space whose inner product is the paper's $\langle B_i\cdot,\cdot\rangle$. Indices are 0-based. The dual norm is the operator norm of linear functionals, and $f'_i(x)$ is the Fréchet derivative of $f$ composed with the block inclusion. Random draws are explicit index sequences, so $\varphi_k$ is the finite average of $f(x_k)$ over all $n^k$ sequences of draws; no measure theory is involved.

Explicit choices where the paper is silent:

• $n\ge1$, every block nontrivial ($n_i\ge1$), and $L_i>0$ are hypotheses.
• Every $Q_i$ is nonempty, and $x_0\in Q$. The paper leaves both implicit; without $x_0\in Q$ the iterates are infeasible and $\varphi_k-f^*$ can be negative.
• $u^{(i)}(x)$ is a choice among the minimizers of the block subproblem. Under the hypotheses the minimizer exists and is unique, so the choice is the paper's update.
• $R_1(x_0)$ is never computed: the theorem takes any $R$ with $\|x-x_*\|_1\le R$ for every feasible $x$ with $f(x)\le f(x_0)$ and every optimal solution $x_*$. This is equivalent to the paper's bound when $R_1(x_0)$ is finite and avoids the default value of a supremum. The level set is taken inside $Q$, the reading of the p. 7 definition for problem (4.1).
• $f^*$ is the constrained optimum $f(x_*)$ for an optimal solution $x_*\in Q$, not an unconstrained minimum.
• Strong convexity is (3.1) on all of $\mathbb R^N$ in $\|\cdot\|_1$. The linear rate does not assume $\sigma\le1$; footnote 2 derives it.
• Printed typos are corrected in the Lean: $U_i^T$ in (4.2) is $U_i$; $u^i$ in (4.3) is $u^{(i)}$.

A formalization in which the update is an arbitrary point of $Q_i$ or an unprojected gradient step, the draws are not uniform, or $f^*$ is the unconstrained minimum, proves a different theorem and is ruled out by the definitions above.

Needed infrastructure: existence and the variational characterization of the minimizer of a strongly convex quadratic over a closed convex set (Mathlib has projections onto closed convex sets in Hilbert spaces), the block descent inequality (2.3) obtained from (2.2), and manipulation of the finite expectations over index sequences (conditioning on the last draw). The block model, the weighted norm and the expectation are shared with the other missions of this series. Proofs of the milestones, and of (2.3) as a supporting lemma, are welcome.

Selected references

• Yu. Nesterov, Efficiency of coordinate descent methods on huge-scale optimization problems, CORE Discussion Paper 2010/2, Université catholique de Louvain, 2010. https://core.ac.uk/download/6430808.pdf
• Yu. Nesterov, Efficiency of coordinate descent methods on huge-scale optimization problems, SIAM Journal on Optimization 22(2) (2012), 341–362. https://doi.org/10.1137/100802001
• Yu. Nesterov, Introductory Lectures on Convex Optimization: A Basic Course, Kluwer, 2004 (Section 2.1, cited in the proof). https://doi.org/10.1007/978-1-4419-8853-9
• P. Tseng and S. Yun, A coordinate gradient descent method for nonsmooth separable minimization, Mathematical Programming 117 (2009), 387–423. https://doi.org/10.1007/s10107-007-0170-0