First Order Methods beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems

TL;DR

This paper extends first order methods to nonconvex, nonsmooth problems without requiring Lipschitz gradient continuity, introducing smooth adaptable functions and applying them to quadratic inverse problems with convergence guarantees.

Contribution

It generalizes the descent lemma and develops a Bregman-based proximal gradient method for nonconvex problems, broadening applicability beyond convex settings.

Findings

Established a full extended descent lemma for nonconvex functions

Proposed a globally convergent Bregman proximal gradient method

Applied framework to quadratic inverse problems with new convergence results

Abstract

We focus on nonconvex and nonsmooth minimization problems with a composite objective, where the differentiable part of the objective is freed from the usual and restrictive global Lipschitz gradient continuity assumption. This longstanding smoothness restriction is pervasive in first order methods (FOM), and was recently circumvent for convex composite optimization by Bauschke, Bolte and Teboulle, through a simple and elegant framework which captures, all at once, the geometry of the function and of the feasible set. Building on this work, we tackle genuine nonconvex problems. We first complement and extend their approach to derive a full extended descent lemma by introducing the notion of smooth adaptable functions. We then consider a Bregman-based proximal gradient methods for the nonconvex composite model with smooth adaptable functions, which is proven to globally converge to a critical point under natural assumptions on the problem's data. To illustrate the power and potential of our general framework and results, we consider a broad class of quadratic inverse problems with sparsity constraints which arises in many fundamental applications, and we apply our approach to derive new globally convergent schemes for this class.

Equations

$$(\mathcal{P})\inf \left\{f\left(x\right)+g\left(x\right):x\in \overline{C}\right\},$$

$$D_h\left(x,y\right):=h\left(x\right)-\left[h\left(y\right)+⟨\nabla h\left(y\right),x-y⟩\right].$$

$$h\text{is convex if and only if}{D}_h(x,y)\ge 0,\forall x\in \mathrm{dom}h,y\in \mathrm{int}\mathrm{dom}h.$$

$$D_h\left(x,z\right)-D_h\left(x,y\right)-D_h\left(y,z\right)=⟨\nabla h\left(y\right)-\nabla h\left(z\right),x-y⟩.$$

$$D_{\alpha h_1+\beta h_2}\left(x,y\right)=\alpha D_{h_1}\left(x,y\right)+\beta D_{h_2}\left(x,y\right),$$

$$\left|g\left(x\right)-g\left(y\right)-⟨\nabla g\left(y\right),x-y⟩\right|\le L{D}_h\left(x,y\right),\forall x,y\in \mathrm{int}\mathrm{dom}h.$$

$$\left|D_g\left(x,y\right)\right|\le L{D}_h\left(x,y\right),\forall x,y\in \mathrm{int}\mathrm{dom}h.$$

$$\left|D_g\left(x,y\right)\right|\equiv \left|g\left(x\right)-g\left(y\right)-⟨\nabla g\left(y\right),x-y⟩\right|\le \frac{L}{2}{\Vert x-y\Vert }^2,\forall x,y\in {\mathbb{R}}^d.$$

$$\exists L>0\text{such that}Lh-g\text{is convex on}\mathrm{int}\mathrm{dom}h⟺D_g\left(x,y\right)\le L{D}_h\left(x,y\right),\forall x,y\in \mathrm{int}\mathrm{dom}h.\text{vspace}-0.1in$$

$$Lh-g=L\left(h-\omega \right)-\left(g-L\omega \right):=L\bar{h}-\bar{g},$$

$$\inf \left\{\Psi \left(x\right)\equiv f\left(x\right)+g\left(x\right):x\in \overline{C}\right\},$$

$$T_{\lambda }\left(x\right)$$

$$=\mathrm{argmin}\left\{f\left(u\right)+⟨\nabla g\left(x\right),u-x⟩+\frac{1}{\lambda }{D}_h\left(u,x\right):u\in {\mathbb{R}}^d\right\},$$

$$\underset{\Vert u\Vert \to \infty }{\lim }\frac{h\left(u\right)+\lambda f\left(u\right)}{\Vert u\Vert }=\infty .$$

$$T_{\lambda }(x)\subset C,\forall x\in C.$$

$${\partial }^{\infty }f\left(x\right)\cap \left(-{\partial }^{\infty }h\left(x\right)\right)=\left\{0\right\},\forall x\in {\mathbb{R}}^d.$$

$${\Psi }_h\left(u\right):=\lambda f\left(u\right)+\lambda ⟨\nabla g\left(x\right),u-x⟩+D_h\left(u,x\right),$$

$${\Psi }_h\left(u\right)=\lambda f(u)+h(u)+⟨\gamma ,u⟩+\rho ,$$

$${\Psi }_h\left(u\right)\ge \lambda f\left(u\right)+h\left(u\right)-\Vert \gamma \Vert \cdot \Vert u\Vert -\left|\rho \right|=\Vert u\Vert \left(\frac{\lambda f\left(u\right)+h\left(u\right)}{\Vert u\Vert }-\Vert \gamma \Vert -\frac{\left|\rho \right|}{\Vert u\Vert }\right).$$

$$x^k\in \mathrm{argmin}\left\{f\left(x\right)+⟨x-x^{k-1},\nabla g\left(x^{k-1}\right)⟩+\frac{1}{\lambda }{D}_h\left(x,x^{k-1}\right):x\in \overline{C}\right\}.$$

$$x^{+}\in \mathrm{argmin}\left\{f\left(u\right)+⟨u-x,\nabla g\left(x\right)⟩+\frac{1}{\lambda }{D}_h\left(u,x\right):u\in \overline{C}\right\},(\lambda >0),$$

$$\lambda \Psi \left(x^{+}\right)\le \lambda \Psi \left(x\right)-\left(1-\lambda L\right)D_h\left(x^{+},x\right).$$

$$f\left(x^{+}\right)+⟨x^{+}-x,\nabla g\left(x\right)⟩+\frac{1}{\lambda }{D}_h\left(x^{+},x\right)\le f\left(x\right).$$

$$g\left(x^{+}\right)+f\left(x^{+}\right)$$

$$\le g\left(x\right)+L{D}_h\left(x^{+},x\right)+f\left(x\right)-\frac{1}{\lambda }{D}_h\left(x^{+},x\right)$$

$$=g\left(x\right)+f\left(x\right)-\left(\frac{1}{\lambda }-L\right)D_h\left(x^{+},x\right),$$

$$\lambda \left(\Psi \left(x^{+}\right)-\Psi \left(u\right)\right)\le D_h\left(u,x\right)-D_h\left(u,x^{+}\right)-\left(1-\lambda L\right)D_h\left(x^{+},x\right)-\lambda D_g(u,x),\forall u\in \mathrm{dom}h.$$

$$\lambda \left(\Psi \left(x^{+}\right)-\Psi \left(u\right)\right)\le \left(1+\lambda L\right)D_h\left(u,x\right)-\left(1-\lambda L\right)D_h\left(x^{+},x\right)-D_h\left(u,x^{+}\right),\forall u\in \mathrm{dom}h.$$

$$\lambda \left(\Psi \left(x^k\right)-\Psi \left(x^{k-1}\right)\right)\le -\left(1-\lambda L\right)D_h\left(x^k,x^{k-1}\right).$$

$$\sum _{k=1}^n{D}_h\left(x^k,x^{k-1}\right)\le \frac{\lambda \left(\Psi \left(x^0\right)-\Psi \left(x^n\right)\right)}{1-\lambda L}\le \frac{\lambda \left(\Psi \left(x^0\right)-{\Psi }_{\ast }\right)}{1-\lambda L},$$

$$n\underset{1\le k\le n}{\min }{D}_h\left(x^k,x^{k-1}\right)\le \sum _{k=1}^n{D}_h\left(x^k,x^{k-1}\right)\le \frac{\lambda \left(\Psi \left(x^0\right)-{\Psi }_{\ast }\right)}{1-\lambda L},$$

$$\underset{1\le k\le n}{\min }{\mathrm{dist}}^2\left(x^{k-1},T_{\lambda }\left(x^{k-1}\right)\right)\le \underset{1\le k\le n}{\min }{\Vert x^k-x^{k-1}\Vert }^2\le \frac{1}{n}\cdot \frac{\lambda \left(\Psi \left(x^0\right)-{\Psi }_{\ast }\right)}{\sigma \left(1-\lambda L\right)}.$$

Full text

2010 Mathematics Subject Classification: Primary 90C25; Secondary 26B25, 49M27, 52A41, 65K05.

Keywords: Composite nonconvex nonsmooth minimization, proximal-gradient algorithms, descent lemma, Non Euclidean distances, Bregman distance, global convergence, Kudyka-Łosiajewicz proeprty, semi-algebraic functions, quadratic inverse problems, phase retrieval.

1 Introduction

The gradient method, forged by Cauchy about 170 years ago, is still at the heart of many fundamental schemes in modern computational optimization through many variants and relatives known as First Order Methods (FOM). These are currently the leading algorithms for solving large scale problems to medium accuracy. The focus and intensity of today research is mostly in the convex case: it goes from complexity to decomposition methods, from Lagrangian approaches to stochastic variants. For an appetizer, see for instance [15, 22, 24] and references therein, as well as the recent comprehensive text of Bertsekas [8] with many relevant up-to-date and annotated sources and references.

A crucial and standard assumption common to almost all FOM is that the gradient of the smooth part in a given objective function has to be globally Lipschitz continuous on the entire underlying space. This is a very restrictive assumption which has only been circumvented by line-search approaches and/or quite complex inner loops which are unavoidably distorting the efficiency and the complexity of the initial method.

Recently, Bauschke, Bolte and Teboulle (BBT) [4] solved this longstanding issue and dealt with non globally Lipschitz continuous gradient. Their framework is simple and far reaching, it is based on adapting the geometry to the objective through the Bregman distance paradigm. It allows them to derive a new Descent Lemma, whereby the usual upper quadratic approximation of a smooth function is replaced by a more general proximity measure which captures, all at once, the geometry of the function and the one of the feasible set. The corresponding FOM derived in [4] come with guaranteed complexity estimates and pointwise global convergence results for convex composite minimization. Very recently, the publication [26] developed independently similar ideas for analyzing the convergence of a Bregman proximal gradient applied to the convex composite model in Banach spaces.

This work is a major departure from the current convex composite model and hence with far reaching consequences. Our main goal is to extend the BBT framework to the nonconvex setting. We complement and extend their approach to study nonconvex composite minimization problems. These consist in minimizing the sum of two nonconvex functions: an extended valued function $f$ and a continuously differentiable function $g$

[TABLE]

where $\overline{C}$ denotes the closure of $C$ which is a nonempty, convex and open set in $\mathbb{R}^{d}$ (see Section 3 for a more precise statement).

Thus, here not only the functions $f$ and $g$ are not convex, but we also consider problems where $g$ does not admit a globally Lipschitz continuous gradient. The usual restrictive requirement of Lipschitz continuity of the gradient of $g$ in problem $(\mathcal{P})$ is not needed, and replaced by a convexity condition which adapts to the geometry of the function $g$. To solve problem $(\mathcal{P})$, we focus on a proximal-based gradient method which involves the non Euclidean distances of Bregman type, and which covers standard proximal-based gradient methods. We derive a class of proximal-based gradient algorithms which are proven to globally converge to a critical point of $(\mathcal{P})$ when $C=\mathbb{R}^{d}$ and when the data is semi-algebraic. Note that one of the earliest work describing and analyzing the classical proximal gradient method in the nonconvex setting goes back to the work of Fukushima and Milne [17]. See also [3, 11] for a full convergence analysis in the semi-algebraic setting, but which also imposed the usual restrictive global Lipschitz continuity of the gradient of $g$. These results are significantly improved in the present paper within the non Euclidean proximal framework.

To illustrate the power of our general framework and results, we consider a broad class of quadratic inverse problems with sparsity constraints which arise in many fundamental applications (see, for instance, [7] and reference therein). We apply our approach to derive new and simple provably convergent schemes, which to the best of our knowledge are the first globally convergent algorithms for this important class of problems.

Outline of the paper. The paper is organized as follows. We first complement and extend the BBT approach to derive a full Descent Lemma by introducing the notion of smooth adaptable functions, see Section 2. The following section presents the mathematical tools necessary to handle the nonconvex setting. It describes the problem and the corresponding proximal- based gradient method for the nonconvex composite model, freed from the usual global Lipschitz gradient continuity restriction. The analysis is developed in Section 4 where the resulting scheme is proven to globally converges to a critical point under natural assumptions on the problem’s data when $C=\mathbb{R}^{d}$. We demonstrate the potential of our framework by showing in Section 5 how it can be successfully applied to a broad class of quadratic inverse problems with sparsity constraints, resulting in computationally simple and explicit iterative formulas. To make the paper self contained, we provide an appendix (see Appendix 6) which includes the relevant material and results for the convergence analysis of algorithm in the semi-algebraic setting.

Notation. We use standard notation and concepts which, unless otherwise specified can all be found in [23].

2 Smooth Adaptable Functions and a Descent Lemma

We begin by defining the notion of smooth adaptable functions. This notion is motivated by the recent work [4] in which a Lipschitz-like/Convexity condition was introduced. This condition allows to lift the usual smoothness of the gradient of a given convex function, and derive a one sided Descent Lemma, whereby the usual upper quadratic approximation of $C^{1}$ functions is given in terms of a more general proximity measure. Here, we extend and complement this notion by accommodating differentiable functions which are not necessarily convex, and we also derive a natural two sided Descent Lemma.

2.1 Preliminaries on Proximity Measures

We first introduce our setting with some notations and definitions that will be used throughout the paper.

Definition 2.1**.**

(Kernel Generating Distance) Let $C$ be a nonempty, convex and open subset of $\mathbb{R}^{d}$. Associated with $C$, a function $h:\mathbb{R}^{d}\rightarrow\left(-\infty,+\infty\right]$ is called a kernel generating distance if it satisfies the following:

• $\rm{(i)}$

$h$ is proper, lower semicontinuous and convex, with $\mathrm{dom}\,{h}\subset\overline{C}$ and $\mathrm{dom}\,{\partial}h=C$.

• $\rm{(ii)}$

$h$ is $C^{1}$ on $\mathrm{int\,dom}\,{h}\equiv C$.

We denote the class of kernel generating distances by $\mathcal{G}(C)$.

Given $h\in\mathcal{G}(C)$, define the proximity measure $D_{h}:\mathrm{dom}\,{h}\times\mathrm{int\,dom}\,{h}\to\mathbb{R}_{+}$

[TABLE]

The proximity measure $D_{h}$ is the so-called Bregman Distance [12]. It measures the proximity between $x$ and $y$. Indeed, thanks to the gradient inequality, one has

[TABLE]

If in addition $h$ is strictly convex, equality holds if and only if $x=y$. However, note that $D_{h}$ is not symmetric in general, unless $h$ is the energy kernel $h=\left(1/2\right)\left\|{\cdot}\right\|^{2}$, which yields the classical squared Euclidean distance. For early foundation papers and key results on Bregman distances, associated proximal-based algorithms, as well as many examples of kernels $h$, generating Bregman distances, we refer the reader to [13, 25, 14, 16, 5].

Note that the structural form of $D_{h}$ is also useful when $h$ is not convex. It measures the discrepancy or error, between the value of $h$ at a given point $x\in\mathrm{dom}\,{h}$ from its linear approximation around $y\in\mathrm{int\,dom}\,{h}$. In that case, obviously, the distance like property $D_{h}\left(x,y\right)\geq 0$ and equal to zero if and only if $x=y$, is no longer valid. However, $D_{h}$ still enjoys two simple, but remarkable properties, which follows from elementary algebra:

• •

The three point identity [14, Lemma 3.1] For any $y,z\in\mathrm{int\,dom}\,{h}$ and $x\in\mathrm{dom}\,{h}$,

[TABLE]

• •

Linear Additivity For any $\alpha,\beta\in\mathbb{R}$, and any functions $h_{1}$ and $h_{2}$ we have:

[TABLE]

for all couple $(x,y)\in\left(\mathrm{dom}\,{h}_{1}\cap\mathrm{dom}\,{h}_{2}\right)^{2}$ such that both $h_{1}$ and $h_{2}$ are differentiable at $y$.

2.2 Smooth Adaptivity and an Extended Descent Lemma

Throughout the paper we will work with a pair of functions $(g,h)$ satisfying:

• $\rm{(i)}$

$h\in\mathcal{G}(C)$.

• $\rm{(ii)}$

$g:\mathbb{R}^{d}\rightarrow\left(-\infty,+\infty\right]$ is a proper and lower semicontinuous function with $\mathrm{dom}\,{h}\subset\mathrm{dom}\,{g}$, which is continuously differentiable on $C=\mathrm{int\,dom}\,{h}$

Definition 2.2** (L-smooth adaptable).**

A pair $(g,h)$ is called $L$-smooth adaptable on $C$ if there exists $L>0$ such that $Lh-g$ and $Lh+g$ are convex on $C$.

Note that the above definition holds for any convex function $h$ which is $C^{1}$ on any open subset of $\mathrm{dom}\,{h}$. The additional properties required in the class $\mathcal{G}(C)$ are not necessary. Only the convexity of $Lh\pm g$ plays a central role. However, for the sake of consistency with the algorithmic development and results that follow, we will always use $h\in\mathcal{G}(C)$.

The above definition allows for immediately obtaining the promised two-sided Descent Lemma which naturally complements and extends the one derived in [4, Lemma 1, p. 333].

Lemma 2.1** (Full Extended Descent Lemma).**

The pair of functions $(g,h)$ is $L$-smooth adaptable on $C$ if and only if:

[TABLE]

Proof.

By using Definition 2.2, the pair $(g,h)$ is $L$-smooth adaptable on $C$ if and only if $Lh-g$ and $Lh+g$ are convex on $\mathrm{int\,dom}\,{h}$. Therefore, thanks to the definition of the Bregman distance and its linear additivity property, this holds if and only if $LD_{h}\left(x,y\right)-D_{g}\left(x,y\right)=D_{Lh-g}\left(x,y\right)\geq 0$, and likewise for $D_{Lh+g}\left(x,y\right)\geq 0$, from which the result immediately follows. ∎

Clearly, in the setting of [4], that is, when the function $g$ is also assumed to be convex, the condition $Lh+g$ is convex, trivially holds, and hence redundant. In this case, the Nolips Descent Lemma given in [4, Lemma 1, p .333] is recovered, i.e., $D_{g}\left(x,y\right)\leq LD_{h}\left(x,y\right)$, though here $g$ needs not be convex.

Using the structural definition of $D_{g}$ (even though $D_{g}$ is not necessarily a proximity measure à la Bregman, since here $g$ is not convex) the full Descent Lemma reads compactly:

[TABLE]

In the special case when the set $C$ is the whole space $\mathbb{R}^{d}$, and $h=\left(1/2\right)\left\|{\cdot}\right\|^{2}$, the classical Descent Lemma for a function $g$ with an $L$-smooth gradient on $\mathbb{R}^{d}$, which provides lower and upper quadratic approximation for $g$, is recovered, i.e.,

[TABLE]

The new Descent Lemma allows for more general lower and upper approximations for $g$ by exploiting the geometry of the set $C$ through the use of the kernel function $h\in\mathcal{G}(C)$ and its associated Bregman proximity measure $D_{h}$.

Remark 2.1**.**

The convexity requirement of $Lh+g$ (and therefore also the left-hand side inequality in (2.1)) can be written with respect to a different parameter $\ell\leq L$. For simplicity, we have used here the same parameter $L$ which appears in the condition $Lh-g$.

Similarly to [4, Proposition 1, p. 334], it is easy to see that the convexity condition on $Lh-g$ (and likewise on $Lh+g$), admit various alternative reformulations which can facilitate its checking. In particular, when both $g$ and $h$ are twice continuously differentiable on $C$, the usual Hessian test can be used, see the applications section for interesting cases.

For the purposes of this paper, and throughout the rest (unless otherwise specified) it will be enough to consider only the condition $Lh-g$ convex, and its corresponding one-sided Descent Lemma. Thus, for convenience we adopt the following short hand terminology for such a pair $(g,h)$:

L-smad holds on $C$ (SMooth ADaptable) if and only if

$\exists\;L>0\;\text{such that}\;Lh-g\;\text{is convex on}\;\mathrm{int\,dom}\,{h}\;\Longleftrightarrow\;D_{g}\left(x,y\right)\leq LD_{h}\left(x,y\right),\;\forall\;x,y\in\mathrm{int\,dom}\,{h}.\vspace{-0.1in}$

We end this section by observing that the L-smad property is invariant when $h$ is assumed $\sigma$-strongly convex, where $\sigma>0$ stands for the strong convexity modulus, namely $h-\left(\sigma/2\right)\left\|{x}\right\|^{2}$ is convex. Indeed, since no convexity is needed/assumed for $g$, we obviously have with $\omega\left(\cdot\right):=\left(1/2\right)\left\|{\cdot}\right\|^{2}$:

[TABLE]

namely, L-smad holds on $C$ for the new pair $\left({\bar{g}},{\bar{h}}\right)$.

3 The Problem and a Bregman Proximal Gradient Algorithm

Our aim is to solve the nonconvex problem $(\mathcal{P})$ with $C:=\mathrm{int\,dom}\,{h}$, that is,

[TABLE]

where the following assumptions on the problem’s data are made throughout the paper.

Assumption A**.**

• $\rm{(i)}$

$h\in\mathcal{G}(C)$ with $\overline{C}=\overline{\mathrm{dom}\,{h}}$.

• $\rm{(ii)}$

$f:\mathbb{R}^{d}\rightarrow\left(-\infty,+\infty\right]$ is a proper and lower semicontinuous function with $\mathrm{dom}\,{f}\cap C\neq\emptyset$.

• $\rm{(iii)}$

$g:\mathbb{R}^{d}\rightarrow\left(-\infty,+\infty\right]$ is a proper and lower semicontinuous function with $\mathrm{dom}\,{h}\subset\mathrm{dom}\,{g}$, which is $C^{1}$ on $C$.

• $\rm{(iv)}$

$v(\mathcal{P}):=\inf\left\{\Psi\left(x\right):\;x\in\overline{C}\right\}>-\infty$.

3.1 A Bregman Proximal Gradient Algorithm

To present the algorithm for tackling problem $(\mathcal{P})$, we define for all $x\in\mathrm{int\,dom}\,{h}$ and any step-size $\lambda>0$, the Bregman proximal gradient mapping

[TABLE]

where the second inequality follows from the fact that $\mathrm{dom}\,{h}\subset\overline{C}$. Note that, since $f$ is nonconvex, the mapping $T_{\lambda}$ is not in general single-valued. This map emerges from the usual approach which consists of linearizing the differential part $g$ around $x$, and regularize it with a proximal distance from that point. Clearly, with $h\equiv\left(1/2\right)\left\|{\cdot}\right\|^{2}$, the above boils down to the classical set-valued proximal gradient mapping.

We now briefly discuss the properties of the Bregman proximal mapping $T_{\lambda}$ in the nonconvex setting. For that purpose, we make the following two additional assumptions.

Assumption B**.**

The function $h+\lambda f$ is supercoercive for all $\lambda>0$, that is,

[TABLE]

Assumption C**.**

For all $x\in C$ we have

[TABLE]

Assumption B is a quite standard coercivity condition, it is for instance automatically satisfied when $\overline{C}$ is compact. On the other hand Assumption C can be shown to hold under a classical qualification condition which ensures the use of a partial sum rule [23, Corollary 10.9, p. 430] and which is formulated through the use of horizon subgradients

[TABLE]

It also holds automatically when $f$ is convex or when $C=\mathbb{R}^{d}$. Another approach to warrant Assumption C is to consider extending the notion of prox-boundedness as defined in [23, Chapter 1, p. 19-20] to Bregman proximal-based map and their envelopes. However, to keep our presentation simple and transparent, these technical issues will not be pursued here, since they are irrelevant in the context of this paper.

We have the following basic result.

Lemma 3.1** (Well-Posedness of $T_{\lambda}$).**

Suppose that Assumptions A, B and C hold, and let $x\in\mathrm{int\,dom}\,{h}$. Then, the set $T_{\lambda}\left(x\right)$ is a nonempty and compact subset of $C$.

Proof.

Fix any $x\in\mathrm{int\,dom}\,{h}$, and $\lambda>0$. For any $u\in\mathbb{R}^{d}$, we define the function

[TABLE]

so that $T_{\lambda}\left(x\right)=\mathrm{argmin}_{u\in\mathbb{R}^{d}}\Psi_{h}\left(u\right)$. Using the definition of the Bregman distance, note that $\Psi_{h}$ can be rewritten as:

[TABLE]

where $\gamma:=\lambda\nabla g\left(x\right)-\nabla h\left(x\right)\in\mathbb{R}^{d}$ and $\rho:=-h\left(x\right)-\left\langle{\gamma,x}\right\rangle\in\mathbb{R}$ are constant quantities. We now show that $\Psi_{h}$ is level bounded on $\mathbb{R}^{d}$, i.e., $\lim_{\left\|{u}\right\|\rightarrow\infty}\Psi_{h}\left(u\right)=\infty$. Using the Cauchy-Schwarz inequality in the above definition of $\Psi_{h}$ we obtain

[TABLE]

Passing to the limit as $\left\|{u}\right\|\rightarrow\infty$, the supercoercivity of $h+\lambda f$ gives $\lim_{\left\|{u}\right\|\rightarrow\infty}\Psi_{h}\left(u\right)=\infty$. Therefore, since $\Psi_{h}$ is also proper and lower semicontinuous, invoking the modern form of Weierstrass’ theorem (see, e.g., [23, Theorem 1.9, p. 11]), it follows that the value $\inf_{\mathbb{R}^{d}}\Psi_{h}$ is finite, and the set $\mathrm{argmin}_{u\in\mathbb{R}^{d}}\Psi_{h}\left(u\right)\equiv T_{\lambda}\left(x\right)$ is nonempty and compact. ∎

Remark 3.1** (The case when $f$ is convex).**

First note that in this case problem ($\mathcal{P}$) remains a nonconvex composite minimization which consists in minimizing the sum of a nonsmooth and convex function $f$ with a nonconvex and continuously differentiable function $g$. However, in this case, under our Assumption A, the function $\Psi_{h}$ which is proper and lower semicontinuous is now convex. It can then be shown (see [4, Lemma 2, p. 336]), that in this case the mapping $T_{\lambda}$ is also single-valued from $\mathrm{int\,dom}\,{h}$ to $\mathrm{int\,dom}\,{h}$.

We are now ready to describe our algorithm for solving the nonconvex composite problem $(\mathcal{P})$.

Bregman Proximal Gradient - BPG

Input. A function $h\in\mathcal{G}(C)$ with $C=\mathrm{int\,dom}\,{h}$ such that L-smad holds on $C$.

Initialization. $x^{0}\in\mathrm{int\,dom}\,{h}$ and let $\lambda>0$.

General Step. For $k=1,2,\ldots$, compute

$x^{k}\in\mathrm{argmin}\left\{f\left(x\right)+\left\langle{x-x^{k-1},\nabla g\left(x^{k-1}\right)}\right\rangle+\frac{1}{\lambda}D_{h}\left(x,x^{k-1}\right):\;x\in\overline{C}\right\}.$

(3.4)

Under our standing Assumptions A, B and C the algorithm is well-defined by Lemma 3.1. In the next section we analyze its properties, and establish a global convergence result to a critical point of $\Psi$.

4 Convergence Analysis of BPG

Throughout this section, we take the following as our blanket assumption

• $\rm{(i)}$

L-smad holds on $C$.

• $\rm{(ii)}$

Assumptions A, B and C hold.

4.1 Properties of the Algorithm

Using the Descent Lemma (see Lemma 2.1), we easily obtain the following key estimation for the composite objective function $\Psi$, which will play an essential role to derive our main convergence results.

Lemma 4.1** (Sufficient Decrease Property).**

For any $x\in\mathrm{int\,dom}\,{h}$, and any $x^{+}\in\mathrm{int\,dom}\,{h}$ defined by

[TABLE]

we have

[TABLE]

In particular, with $0<\lambda L<1$, the sufficient decrease in the composite objective function value $\Psi$ is ensured.

Proof.

By the definition of global optimality for (4.1), we obtain by taking $u=x\in\mathrm{int\,dom}\,{h}$, that

[TABLE]

Invoking the Descent Lemma (see Lemma 2.1) for $g$ we then get with the above,

[TABLE]

and hence with $\Psi=f+g$, the desired inequality for $\Psi$ is proved, and the last statement immediately follows with $0<\lambda L<1$. ∎

Remark 4.1**.**

• $\rm{(i)}$

When $f$ is assumed convex, using the global optimality condition for the corresponding convex problem, which defines $x^{+}$ through (4.1), followed by using the three points identity property of a Bregman distance, and the L-smad property, one easily see that for any $x\in\mathrm{int\,dom}\,{h}$, the inequality (4.2) can be strengthened and reads

[TABLE]

When $g$ is also convex, then the last term $-\lambda D_{g}\left(u,x\right)\leq 0$, and the above result recovers the key estimation result proven in [4, Lemma 5, p. 340].

• $\rm{(ii)}$

When $g$ is nonconvex and the Full Descent Lemma (see Lemma 2.1) holds, i.e., there exists $L>0$ with $Lh+g$ convex, then we have $-D_{g}\left(u,x\right)\leq LD_{h}\left(u,x\right)$ and hence the above inequality reduces to

[TABLE]

The result given in Lemma 4.1 is valid for any $\lambda>0$, and as seen, imposing the condition $0<\lambda L<1$ ensures a sufficient decrease in the composite objective function value $\Psi$. This fact yields the following result.

Proposition 4.1**.**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a sequence generated by BPG with $0<\lambda L<1$. Then the following assertions hold:

• $\rm{(i)}$

The sequence $\left\{\Psi\left(x^{k}\right)\right\}_{k\in\mathbb{N}}$ is nonincreasing.

• $\rm{(ii)}$

$\sum_{k=1}^{\infty}D_{h}\left(x^{k},x^{k-1}\right)<\infty$, and hence the sequence $\left\{D_{h}\left(x^{k},x^{k-1}\right)\right\}_{k\in\mathbb{N}}$ converges to zero.

• $\rm{(iii)}$

$\min_{1\leq k\leq n}D_{h}\left(x^{k},x^{k-1}\right)\leq\frac{\lambda}{n}\left(\frac{\Psi\left(x^{0}\right)-\Psi_{\ast}}{1-\lambda L}\right),$* where $\Psi_{\ast}=v(\mathcal{P})>-\infty$ (by Assumption A(iv)).*

Proof.

• $\rm{(i)}$

Fix $k\geq 1$. Under our assumptions and using the iterative step (3.4), we can apply Lemma 4.1 with $x=x^{k-1}$ and $x^{+}=x^{k}$, to obtain

[TABLE]

Thus, with $0<\lambda L<1$, we immediately obtain that the sequence $\left\{\Psi\left(x^{k}\right)\right\}_{k\in\mathbb{N}}$ is nonincreasing.

• $\rm{(ii)}$

Let $n$ be a positive integer. Summing the above inequality from $k=1$ to $n$ we get

[TABLE]

where $\Psi_{\ast}=v(\mathcal{P})>-\infty$. Taking the limit as $n\rightarrow\infty$, we obtain the first desired assertion (ii), from which we immediately deduce that $\left\{D_{h}\left(x^{k},x^{k-1}\right)\right\}_{k\in\mathbb{N}}$ converges to zero.

• $\rm{(iii)}$

From (4.5) we also obtain,

[TABLE]

which after division by $n$ yields the desired result.

∎

Clearly, with $h$ being the energy kernel on $C\equiv\mathbb{R}^{d}$, the above proposition yields all the classical results for the nonconvex proximal gradient method (e.g., set $\lambda L=1/2$). In particular, in this case one obtains that the corresponding classical gradient mapping , defined by $G\left(x^{k-1}\right):=\left\|{x^{k-1}-T_{\lambda}\left(x^{k-1}\right)}\right\|$ (see [6, Theorem 2.3, p. 61]) converges to zero at a rate of $O(1/\sqrt{n})$.

Finally, note that thanks to the invariance of $L$-smad (cf. Section 2.2), we can assume that $h$ is $\sigma$-strongly convex on $C$. In that case, for any nonempty set $S\subset\mathbb{R}^{d}$, let $\mathrm{dist}(y,S):=\inf\left\{\left\|{u-y}\right\|:\;u\in S\right\}$. Since $x^{k}\in T_{\lambda}\left(x^{k-1}\right)$, using the $\sigma$-strong convexity of $h$ combined with Proposition 4.1(iii) we immediately obtain the following rate of convergence result for two successive iterates:

[TABLE]

4.2 Global Convergence of BPG

Throughout this section we consider problem $(\mathcal{P})$ defined on $C\equiv\mathbb{R}^{d}$, namely

[TABLE]

Here, throughout this subsection, we additionally assume that

Assumption D**.**

• $\rm{(i)}$

$\mathrm{dom}\,{h}=\mathbb{R}^{d}$ and $h$ is $\sigma$-strongly convex on $\mathbb{R}^{d}$.

• $\rm{(ii)}$

$\nabla h$ and $\nabla g$ are Lipschitz continuous on any bounded subset of $\mathbb{R}^{d}$.

Noting that Assumption C is automatically fulfilled since $C=\mathbb{R}^{d}$. As usual we use the concept of limiting subdifferential for $f$ so that thanks to Fermat’s rule [23, Theorem 10.1, p. 422], the set of critical points of $\Psi$ is given by:

[TABLE]

To prove the global convergence of the sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ generated by BPG to a critical point of $\Psi$, we apply the methodology developed in [11]. For the reader’s convenience, and to make this paper self-contained, the main tools and relevant proofs are given in Appendix 6. First, we describe three key ingredients of the methodology [11].

Definition 4.1** (Gradient-like descent sequence).**

Let $F:\mathbb{R}^{d}\to(-\infty,\infty]$ be a proper and lower semicontinuous function. A sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is called a gradient-like descent sequence for $F$ if the following three conditions hold:

• $\rm{(C1)}$

Sufficient decrease property. There exists a positive scalar $\rho_{1}$ such that

[TABLE]

• $\rm{(C2)}$

A subgradient lower bound for the iterates gap. There exists a positive scalar $\rho_{2}$ such that

[TABLE]

• $\rm{(C3)}$

Let $\overline{x}$ be a limit point of a subsequence $\left\{x^{k}\right\}_{k\in{\cal K}}$, then $\limsup_{k\in{\cal K}\subset\mathbb{N}}F\left(x^{k}\right)\leq F\left(\overline{x}\right)$.

The two conditions (C1) and (C2) are typical for any descent type algorithm (see, e.g., [1]). They are also the main tools to prove subsequential convergence, as recorded below. The condition (C3) is a minimal and very weak continuity requirement, which, in particular, holds when $F$ is continuous. We denote by $\omega\left(x^{0}\right)$ the set of all limit points of $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$.

Lemma 4.2** (Abstract subsequence convergence).**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a bounded gradient-like descent sequence for $F$. Then, $\omega\left(x^{0}\right)$ is a nonempty and compact subset of $\mathrm{crit}\,{F}$, and we have

[TABLE]

In addition, the function $F$ is finite and constant on $\omega\left(x^{0}\right)$.

Proof.

See Appendix 6. ∎

Remark 4.2** (Boundedness of Sequences).**

Observe that the mere coercivity of $F$ ensures that any gradient-like sequence is bounded.

The above, together with the so-called nonsmooth Kurdyka-Łojasiewicz (KL) property (see [9] and Appendix 6 for details) allows us to establish our main convergence results of BPG.

Theorem 4.1** (Convergence of BPG).**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a sequence generated by BPG which is assumed to be bounded and let $0<\lambda L<1$. The following assertions hold.

• $\rm{(i)}$

Subsequential convergence.* Any limit point of the sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is a critical point of $\Psi$.*

• $\rm{(ii)}$

Global convergence.* Suppose that $\Psi$ satisfies the KL property on $\mathrm{dom}\,{\Psi}$. Then, the sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ has finite length and converges to a critical point $x^{\ast}$ of $\Psi$.*

Proof.

(i) We first prove that the sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is a gradient-like descent sequence; namely it satisfies the three conditions (C1), (C2) and (C3) as described in Definition 6.1, and the result will then be established by invoking Lemma 4.2.

From Proposition 4.1 (see (4.4)) we have that

[TABLE]

where the last inequality follows from the $\sigma$-strong convexity of $h$ (see Assumption D(i)). This proves that condition (C1) holds true.

Writing the optimality condition of the optimization problem which defines $x^{k+1}$ (see (3.4)) yields that

[TABLE]

Therefore, by defining

[TABLE]

we obviously obtain that $w^{k+1}\in\partial\Psi\left(x^{k+1}\right)$. Since $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is a bounded sequence and both $\nabla h$ and $\nabla g$ are Lipschitz continuous on bounded subsets of $\mathbb{R}^{d}$ (see Assumption D(ii)), there exists $M>0$ such that

[TABLE]

This proves that condition (C2) also holds true.

Consider a subsequence $\left\{x^{n_{k}}\right\}_{n\in\mathbb{N}}$ which converges to some point $x^{\ast}$ (there exists such a subsequence since the sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is assumed to be bounded). Using (4.7) and Proposition 4.1(iii) we obtain that $\lim_{k\rightarrow\infty}\left\|{x^{k}-x^{k-1}}\right\|=0$. Therefore, the sequence $\left\{x^{n_{k}-1}\right\}_{k\in\mathbb{N}}$ also converges to $x^{\ast}$. In addition, since $h$ is continuously differentiable on $\mathbb{R}^{d}$ we have that $\lim_{k\rightarrow\infty}D_{h}\left(x^{\ast},x^{k-1}\right)=0$. Now, from (3.4), it follows that

[TABLE]

that is,

[TABLE]

Substituting $k$ by $n_{k}$ and letting $k\rightarrow\infty$, we obtain from Proposition 4.1(iii) and the facts mentioned above, that

[TABLE]

This proves condition (C3), and thanks to Lemma 4.2, the first item of the theorem follows. Moreover, we also obtain that $\Psi$ is constant on $\omega(x^{0})$.

(ii) Since $\Psi$ is assumed to have the KL property, and we have just shown that $\Psi$ is constant on $\omega(x^{0})$, we can then apply Theorem 6.2 given in Appendix 6, with $F:=\Psi$. This proves that the sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ has finite length, and is globally convergent to a critical point $x^{\ast}$ of $\Psi$. ∎

We end by noting that, when the data $f$ and $g$ are semi-algebraic, the convergence rate of the generated sequence described in Theorem 4.1 follows by invoking Theorem 6.3, see the Appendix 6 for details.

5 Application to Quadratic Inverse Problems

5.1 Motivation and Problem’s Statement

This section illustrates the potential of our approach and results. To this end, we consider the broad class of problems which consists of solving approximately a system of quadratic equations and can be described as follows. Given a finite number of symmetric matrices $A_{i}\in\mathbb{R}^{d\times d}$, $i=1,2,\ldots,m$, describing the model, and a vector of possibly noisy measurements $b\in\mathbb{R}^{m}$, the goal is to find $x\in\mathbb{R}^{d}$, that solves the following system

[TABLE]

This type of problems is a natural extension of the classical linear inverse problem which arise in the broad area of signal recovery, whereby linear measurements are replaced by quadratic measurements, see [7] and references therein. It also includes the well-known and fundamental class of phase retrieval problems as a special case, which has been, and is still intensively studied in the literature; see [20] for a very recent review on this problem and references therein.

Commonly, the system under consideration is underdetermined, thus some prior information on the original signal is incorporated into the model by means of some regularizer represented by a function $f$, possibly nonconvex, nonsmooth, and extended valued (to accommodate constraints). Adopting the usual least-squares model to measure the error, the problem can then be reformulated as the nonconvex optimization problem

[TABLE]

where $\theta>0$ plays the role of a penalty parameter controlling the trade-off between matching the data fidelity criteria versus its regularizer $f$.

Clearly, the nonconvex function $g:\mathbb{R}^{d}\rightarrow\left(-\infty,+\infty\right]$ defined by $g\left(x\right):=\left(1/4\right)\sum_{i=1}^{m}\left(x^{T}A_{i}x-b_{i}\right)^{2}$ is continuously differentiable on $\mathbb{R}^{d}$, and it does not admit a global Lipschitz continuous gradient, thus precluding the use of the classical proximal gradient method.

In the following, we demonstrate the power of our approach, illustrating how new and simple globally convergent schemes can be derived for the broad class of problems (QIP). In particular, we design two new algorithms: one for the unconstrained $\ell_{1}$ regularized (QIP) model, and the other for the sparsity $\ell_{0}$ constrained (QIP); both models being known to promote sparse solutions. Our algorithms are given through an explicit closed form formula, and are thus straightforward to implement. To the best of our knowledge, they appear to be the first globally convergent schemes for this class of problems.

5.2 Two Algorithms for Sparse Quadratic Inverse Problems

Throughout this section, the underlying space is $C\equiv\mathbb{R}^{d}$. Given the nonconvex function $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$ defined by

[TABLE]

we consider solving the nonconvex model (QIP) in the following two cases:

• (a)

A convex $\ell_{1}$-norm regularization. $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ with $f\left(x\right)=\left\|{x}\right\|_{1}$.

• (b)

A nonconvex $\ell_{0}$-ball constraint. $f:\mathbb{R}^{d}\rightarrow\left(-\infty,+\infty\right]$ with $f\left(x\right)=\delta_{\mathbb{B}_{0}^{s}}\left(x\right)$, the indicator of the $\ell_{0}$-ball, where for a positive integer $s<d$,

[TABLE]

and $\left\|{x}\right\|_{0}$ is the “$\ell_{0}$ norm” which counts the number of nonzero components in $x$.

To apply our approach on the (QIP) model in both cases (a) and (b), we first need to identify a suitable function $h\in\mathcal{G}(\mathbb{R}^{d})$ such that L-smad holds for the pair $\left(g,h\right)$. Here, we use $h:\mathbb{R}^{d}\rightarrow\mathbb{R}$ given by

[TABLE]

Equipped with this $h$, we now show that L-smad holds, i.e., there exists $L>0$ such that $Lh-g$ is convex on $\mathbb{R}^{d}$.

Lemma 5.1**.**

Let $g$ and $h$ as defined above. Then, for any $L$ satisfying

[TABLE]

the function $Lh-g$ is convex on $\mathbb{R}^{d}$.

Proof.

Let $x\in\mathbb{R}^{d}$. Since $g$ and $h$ are $C^{2}$ on $\mathbb{R}^{d}$, in order to warrant the convexity of $Lh-g$ it is sufficient to find $L>0$ such that $L\lambda_{\min}\left(\nabla^{2}h\left(x\right)\right)\geq\lambda_{\max}\left(\nabla^{2}g\left(x\right)\right)$, where $\lambda_{\min}(M)$ and $\lambda_{\max}(M)$ denote the minimal and maximal eigenvalue of a matrix $M$, respectively.

By a straightforward computation we obtain that

[TABLE]

and due to the well-known fact that $\lambda_{\max}(M)\leq\left\|{M}\right\|$ we can write

[TABLE]

On the other hand, since $\nabla^{2}h\left(x\right)\succeq\left(\left\|{x}\right\|_{2}^{2}+1\right)I_{d}$, we obtain that $\lambda_{\min}\left(\nabla^{2}h\left(x\right)\right)\geq\left\|{x}\right\|_{2}^{2}+1$. Therefore, taking $L\geq\sum_{i=1}^{m}\left(3\left\|{A_{i}}\right\|^{2}+\left\|{A_{i}}\right\|\left|b_{i}\right|\right)$ yields

[TABLE]

which proves the desired result. ∎

In order to apply our results (cf. Section 4.2) for both scenarios (a) and (b), we observe that $h$ given above is $1$-strongly convex on $\mathbb{R}^{d}$ and it easy to see that AssumptionsA-D hold. Furthermore, the function $g$ is a real polynomial, hence semi-algebraic, and both functions $\left\|{x}\right\|_{0}$ and $\left\|{x}\right\|_{1}$ is also semi-algebraic (see, [11, Appendix 5, p. 490]). Therefore, since the addition of semi-algebraic functions results in a semi-algebraic function, it follows that for both models (a) and (b), the objective $\Psi$ is semi-algebraic, and BPG can be applied on the (QIP) model in each case, to produce a globally convergent sequence which converges to a critical point of $\Psi$.

To apply BPG, the main computational step require us to compute the Bregman proximal gradient map (see (3.1)):

[TABLE]

We now show that for both scenarios (a) and (b), this step yields an explicit closed form formula.

Before doing so, we introduce some convenient notations and recall some well known operators that will be used in the rest of this section. Let $\lambda>0$ and fix any $x\in\mathbb{R}^{d}$. Define

[TABLE]

For the pair $(g,h)$ given above, a direct computation of their gradients gives $p_{\lambda}\left(x\right)$. Now, omitting constant terms in the above expression defining $T_{\lambda}$, we can write

[TABLE]

We now recall two very well-known operators that will be used to compute the resulting $T_{\lambda}$.

• Soft-thresholding (with parameter $\tau$). For any $y\in\mathbb{R}^{d}$,

[TABLE]

with the absolute value understood to be component-wise.

• Hard-thresholding (with parameter $\tau$). For any $y\in\mathbb{R}^{d}$,

[TABLE]

where we assumed, without the loss of generality, that the vector $y$ is sorted in a decreasing order according to the absolute value of its entries (ties may be broken arbitrarily), i.e., $|y_{1}|\geq|y_{2}|\geq\cdots\geq|y_{n}|$.

We are now ready to establish the promised formula of $T_{\lambda}$ for the two cases (a) and (b).

Proposition 5.1** (Bregman Proximal Formula for the $\ell_{1}$-Norm Regularization).**

Let $f=\left\|{\cdot}\right\|_{1}$, and for $x\in\mathbb{R}^{d}$, let $v(x):=\mathcal{S}_{\lambda\theta}\left(p_{\lambda}\left(x\right)\right)$. Then, $x^{+}=T_{\lambda}\left(x\right)$ is given by

[TABLE]

where $t^{\ast}$ is the unique positive real root of

[TABLE]

which admits an explicit formula.555The general solution of a cubic equation was apparently discovered first by Scipione del Ferro and published by Gerolamo Cardano in his Ars magna in 1545, see https://www.britannica.com/topic/cubic-equation.

Proof.

From the formulation of $T_{\lambda}$ given in (5.2), we have

[TABLE]

where we recall that $p=\lambda\nabla g\left(x\right)-\nabla h\left(x\right)$. By the first order global optimality condition for the strongly convex problem (5.5) we obtain for each $i=1,2,\ldots,d$,

[TABLE]

where $\gamma_{i}=\mathrm{sgn}(x_{i}^{+})$ if $x_{i}^{+}\neq 0$, and $\gamma_{i}\in\left[-1,1\right]$ if $x_{i}^{+}=0$. For each $i=1,2,\ldots,d$, set $v_{i}=p_{i}+\lambda\theta\gamma_{i}$. Then, from the above equation, $x_{i}^{+}=-tv_{i}$ with $t>0$, and

[TABLE]

Now consider the following two cases for any $i=1,2,\ldots,d$:

• •

If $v_{i}=0$, then $x_{i}^{+}=0$, and $\gamma_{i}\in\left[-1,1\right]$, i.e., we obtain $x_{i}^{+}=0$ whenever $\left|p_{i}\right|\leq\lambda\theta$.

• •

If $v_{i}\neq 0$, then $x_{i}^{+}\neq 0$, and $\gamma_{i}=\mathrm{sgn}(x_{i}^{+})=\mathrm{sgn}(-tv_{i})=-\mathrm{sgn}(v_{i})$. Therefore, using $\gamma_{i}=\left(v_{i}-p_{i}\right)/\left(\lambda\theta\right)$, the latter relation reads $\lambda\theta\mathrm{sgn}(v_{i})+v_{i}-p_{i}=0$, which is nothing else but the following optimality condition for $v$:

[TABLE]

where the last equality is from the definition of the soft threshold operator (see (5.3).

To summarize, in all cases, we have thus obtained that $x^{+}=-t^{\ast}\mathcal{S}_{\lambda\theta}\left(p\right)$, where $t^{\ast}$ is the positive root of the cubic equation $t^{3}\left\|{v}\right\|_{2}^{2}+t-1=0$, which is given by a closed form expression. ∎

Now, we turn to the $\ell_{0}$-norm constrained version. First we recall the following useful result established in [21, Proposition 4.3 p. 79].

Lemma 5.2**.**

Given $0\neq a\in\mathbb{R}^{d}$ and a positive integer $s<d$, then

[TABLE]

with optimal solution obtained at $z^{\ast}=\mathcal{H}_{s}\left(a\right)/\left\|{\mathcal{H}_{s}\left(a\right)}\right\|_{2}$.

Proposition 5.2** (Bregman Proximal Formula for the $\ell_{0}$-Ball Constraint).**

Let $f=\delta_{\mathbb{B}_{0}^{s}}$ and $x\in\mathbb{R}^{d}$. Then, $x^{+}=T_{\lambda}\left(x\right)$ is given by

[TABLE]

where $\sqrt{t^{\ast}}\equiv\eta^{\ast}$ is the unique nonnegative real root of the cubic equation

[TABLE]

Proof.

We use the same notations as in Proposition 5.1. Here $x^{+}=T_{\lambda}\left(x\right)$ consists in finding $x^{+}$ which solves the nonconvex problem:

[TABLE]

Introducing the new variable $t:=\left\|{u}\right\|_{2}^{2}\geq 0$, the latter is equivalent to solving

[TABLE]

Invoking Lemma 5.2 for the inner minimization with respect to $u$, one easily obtains:

[TABLE]

with solution, $u^{\ast}=\sqrt{t}\left\|{\mathcal{H}_{s}\left(p\right)}\right\|_{2}^{-1}\mathcal{H}_{s}\left(p\right)$ where we used the fact that for the hard threshold operator: $\mathcal{H}_{s}\left(p\right)=\mathcal{H}_{s}\left(-p\right)$). Therefore, problem (5.7) reduces to solve

[TABLE]

and it follows that $x^{+}=\sqrt{t^{\ast}}\left\|{\mathcal{H}_{s}\left(p\right)}\right\|_{2}^{-1}\mathcal{H}_{s}\left(p\right)$, where $t^{\ast}$ is an optimal solution of (5.8), thus proving the first part of the proposition. Now, if $\left\|{\mathcal{H}_{s}\left(p\right)}\right\|_{2}=0$ then clearly $t^{\ast}=0$. Otherwise, with $\left\|{\mathcal{H}_{s}\left(p\right)}\right\|_{2}>0$, since $\lim_{t\rightarrow 0^{+}}\psi^{\prime}\left(t\right)=-\infty$ and the above objective $\psi$ is strongly convex over $\mathbb{R}_{+}$, the problem (5.8) admits a unique minimizer $t^{\ast}>0$ which satisfies the optimality condition:

[TABLE]

Thus, in either cases, the unique optimal solution $t^{\ast}$ for problem (5.8) is the non-negative real root of the last equation, which can be rewritten as (5.6), and hence the desired result is proved. ∎

Equipped with Lemma 5.1 together with Propositions 5.1 and 5.2, thanks to Theorem 4.1, we have now all the ingredients to formulate simple explicit algorithms in terms of the problem’s data for solving a broad class of sparse quadratic inverse problems (SQIP). A thorough computational study of (SQIP), and some other related extensions will be developed in a separate paper.

6 Appendix: Global Convergence for KL Functions

The objective of this appendix is twofold. First, to make the paper self-contained, and second to outline here the general convergence mechanism as recently described in [11] on the so-called PALM algorithm, so that it can be used and applied to any given algorithm as described below. To this end, let $F:\mathbb{R}^{d}\rightarrow\left(-\infty,+\infty\right]$ be a proper and lower semicontinuous function which is bounded from below and consider the problem

[TABLE]

Consider a generic algorithm $\mathcal{A}$ which generates a sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ via the following:

[TABLE]

The main goal is to prove that the whole sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$, generated by the algorithm $\mathcal{A}$, converges to a critical point of $F$. For that purpose we first outline the three key ingredients of the forthcoming methodology.666To keep this appendix as a completely independent unit, we repeat here a definition and results already stated in Section 4.

Definition 6.1** (Gradient-like descent sequence).**

A sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is called a gradient-like descent sequence for $F$ if the following three conditions hold:

• $\rm{(C1)}$

Sufficient decrease property. There exists a positive scalar $\rho_{1}$ such that

[TABLE]

• $\rm{(C2)}$

A subgradient lower bound for the iterates gap. There exists a positive scalar $\rho_{2}$ such that

[TABLE]

• $\rm{(C3)}$

Let $\overline{x}$ be a limit point of a subsequence $\left\{x^{k}\right\}_{k\in{\cal K}}$, then $\limsup_{k\in{\cal K}\subset\mathbb{N}}F\left(x^{k}\right)\leq F\left(\overline{x}\right)$.

The two conditions (C1) and (C2) are typical for any descent type algorithm (see, e.g., [1]), and basic to prove subsequential convergence. The condition (C3) is a minimal and weak requirement, which, in particular, holds when $F$ is continuous.

The set of all limit points of $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is defined by

[TABLE]

The next result establish the promised subsequential convergence.

Lemma 6.1** (Subsequence Convergence).**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a gradient-like descent sequence for $F$ which is assumed to be bounded. Then, $\omega\left(x^{0}\right)$ is a nonempty and compact subset of $\mathrm{crit}\,{F}$, and we have

[TABLE]

In addition, the objective function $F$ is finite and constant on $\omega\left(x^{0}\right)$.

Proof.

Since $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is bounded there is $x^{\ast}\in\mathbb{R}^{d}$ and a subsequence $\left\{x^{k_{q}}\right\}_{q\in\mathbb{N}}$ such that $x^{k_{q}}\rightarrow x^{\ast}$ as $q\rightarrow\infty$ and hence $\omega\left(x^{0}\right)$ is nonempty. Moreover, the set $\omega\left(x^{0}\right)$ is compact since it can be viewed as an intersection of compact sets. Now, from condition (C3) and the lower semicontinuity of $F$, we obtain

[TABLE]

On the other hand, from conditions (C1) and (C2), we know that there is $w^{k}\in\partial F\left(x^{k}\right)$, $k\in\mathbb{N}$, such that $w^{k}\rightarrow 0$ as $k\rightarrow\infty$. The closedness property777Let $\left\{\left(x^{k},u^{k}\right)\right\}_{k\in\mathbb{N}}$ be a sequence in $\mathrm{graph}\,{\left(\partial F\right)}$ that converges to $\left(x,u\right)$ as $k\rightarrow\infty$. By the very definition of $\partial F\left(x\right)$, if $F(x^{k})$ converges to $F(x)$ as $k\rightarrow\infty$, then $\left(x,u\right)\in\mathrm{graph}\,{(\partial F)}$. of $\partial F$ implies thus that $0\in\partial F\left(x^{\ast}\right)$. This proves that $x^{\ast}$ is a critical point of $F$, and hence (6.1) is valid.

To complete the proof, let $\lim_{{k}\rightarrow{\infty}}F\left(x^{k}\right)=l\in\mathbb{R}$. Then $\left\{F\left(x^{k_{q}}\right)\right\}_{q\in\mathbb{N}}$ converges to $l$ and from (6.2) we have $F\left(x^{\ast}\right)=l$. Hence the restriction of $F$ to $\omega\left(x^{0}\right)$ equals $l$. ∎

To achieve our main goal, i.e., to establish global convergence of the whole sequence, we need an additional assumption on the class of functions $F$: it must satisfy the so-called nonsmooth Kurdyka-Łojasiewicz (KL) property [9](see [18, 19] for smooth cases). We refer the reader to [10] for an in depth study of the class of KL functions, as well as references therein. We provide now the formal definition of the KL property and two important results.

Denote $[\alpha<F<\beta]:=\left\{x\in\mathbb{R}^{d}:\;\alpha<F\left(x\right)<\beta\right\}$. Let $\eta>0$, and set

[TABLE]

Definition 6.2** (The nonsmooth KL property).**

A proper and lower semicontinuous function $F:\mathbb{R}^{d}\rightarrow\left(-\infty,+\infty\right]$ has the Kurdyka-Łojasiewicz (KL) property locally at $\overline{u}\in\mathrm{dom}\,{F}$ if there exist $\eta>0$, $\varphi\in\Phi_{\eta}$, and a neighborhood $U\left(\overline{u}\right)$ such that

[TABLE]

for all $u\in U\left(\overline{u}\right)\cap\left[F\left(\overline{u}\right)<F\left(u\right)<F\left(\overline{u}\right)+\eta\right]$.

Verifying the KL property of a given function might often be a difficult task. However, thanks to a fundamental result established in [9], it holds for the broad class of semi-algebraic functions.

Theorem 6.1**.**

Let $F:\mathbb{R}^{d}\rightarrow\left(-\infty,+\infty\right]$ be a proper and lower semicontinuous function. If $F$ is semi-algebraic then it satisfies the KL property at any point of $\mathrm{dom}\,{F}$.

Our last ingredient is a key uniformization of the KL property proven in [11, Lemma 6, p. 478], which we record below.

Lemma 6.2** (Uniformized KL property).**

Let $\Omega$ be a compact set and let $F:\mathbb{R}^{d}\rightarrow\left(-\infty,+\infty\right]$ be a proper and lower semicontinuous function. Assume that $F$ is constant on $\Omega$ and satisfies the KL property at each point of $\Omega$. Then, there exist $\varepsilon>0$, $\eta>0$ and $\varphi\in\Phi_{\eta}$ such that for all $\overline{x}$ in $\Omega$ one has,

[TABLE]

and all $x\in\left\{x\in\mathbb{R}^{d}:\;\mathrm{dist}\left(x,\Omega\right)<\varepsilon\right\}\cap\left[F\left(\overline{x}\right)<F\left(x\right)<F\left(\overline{x}\right)+\eta\right]$.

We can now conveniently summarize the convergence results of [11] through the following abstract convergence result.

Theorem 6.2** (Global convergence).**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a bounded gradient-like descent sequence for $F$. If $F$ satisfies the KL property, then the sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ has finite length, i.e., $\sum_{k=1}^{\infty}\left\|{x^{k+1}-x^{k}}\right\|<\infty$ and it converges to $x^{\ast}\in\mathrm{crit}\,{F}$.

Proof.

Since $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is bounded there exists a subsequence $\left\{x^{k_{q}}\right\}_{q\in\mathbb{N}}$ such that $x^{k_{q}}\rightarrow\overline{x}$ as $q\rightarrow\infty$. In a similar way as in Lemma 4.2 we get that

[TABLE]

If there exists an integer $\bar{k}$ for which $F(x^{\bar{k}})=F\left(\overline{x}\right)$ then condition (C1) would imply that $x^{\bar{k}+1}=x^{\bar{k}}$. A trivial induction show then that the sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is stationary and the announced results are obvious. Since $\left\{F\left(x^{k}\right)\right\}_{k\in\mathbb{N}}$ is a nonincreasing sequence, it is clear from (6.4) that $F\left(\overline{x}\right)<F\left(x^{k}\right)$ for all $k>0$. Again from (6.4) for any $\eta>0$ there exists a nonnegative integer $k_{0}$ such that $F\left(x^{k}\right)<F\left(\overline{x}\right)+\eta$ for all $k>k_{0}$. From Lemma 4.2 we know that $\lim_{k\rightarrow\infty}\mathrm{dist}\left(x^{k},\omega\left(x^{0}\right)\right)=0$. This means that for any $\varepsilon>0$ there exists a positive integer $k_{1}$ such that $\mathrm{dist}\left(x^{k},\omega\left(x^{0}\right)\right)<\varepsilon$ for all $k>k_{1}$.

From Lemma 4.2, we know that $\omega\left(x^{0}\right)$ is nonempty and compact, the function $F$ is finite and constant on $\omega\left(x^{0}\right)$. Hence, we can apply the Uniformization Lemma with $\Omega=\omega\left(x^{0}\right)$. Therefore, for any $k>l:=\max\left\{k_{0},k_{1}\right\}$, we have

[TABLE]

This makes sense since we know that $F\left(x^{k}\right)>F\left(\overline{x}\right)$ for any $k>l$. From condition (C2) we get that

[TABLE]

For convenience, we define for all $p,q\in\mathbb{N}$ and $\overline{x}$ the following quantity

[TABLE]

From the concavity of $\varphi$ we get that

[TABLE]

Combining condition (C1) with (6.6) and (6.7) yields, for any $k>l$, that

[TABLE]

Using the fact that $2\sqrt{\alpha\beta}\leq\alpha+\beta$ for all $\alpha,\beta\geq 0$, we infer from the later inequality that

[TABLE]

Summing up (6.9) for $i=l+1,\ldots,k$ yields

[TABLE]

where the last inequality follows from the fact that $\Delta_{p,q}+\Delta_{q,r}=\Delta_{p,r}$ for all $p,q,r\in\mathbb{N}$. Since $\varphi\geq 0$, recalling the definition of $\Delta_{l+1,k+1}$, we thus have for any $k>l$ that

[TABLE]

which implies that $\sum_{k=1}^{\infty}\left\|{x^{k+1}-x^{k}}\right\|<\infty$, i.e., $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is a Cauchy sequence and hence together with Lemma 4.2, we obtain the global convergence to a critical point is established. ∎

Finally, we record a generic rate of convergence result (see [1]) which is also valid under the same assumptions as Theorem 6.2.

Theorem 6.3** (Rate of convergence).**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a bounded gradient-like descent sequence for $F$. Assume that $F$ satisfies the KL property where the desingularizing function $\varphi$ of $F$ is of the following form

[TABLE]

Let $\overline{x}$ be the limit point of $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$. Then the following estimations hold:

• $\rm{(i)}$

If $\theta=0$ then the sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ converges to $\overline{x}$ in a finite number of steps.

• $\rm{(ii)}$

If $\theta\in\left(0,1/2\right]$ then there exist $\omega>0$ and $\tau\in\left[0,1\right)$ such that

[TABLE]

• $\rm{(iii)}$

If $\theta\in\left(1/2,1\right)$ then there exist $\omega>0$ such that

[TABLE]