Convex-Concave Backtracking for Inertial Bregman Proximal Gradient Algorithms in Non-Convex Optimization

TL;DR

This paper introduces a novel double convex-concave backtracking method for inertial Bregman proximal gradient algorithms, enabling adaptive step size and extrapolation control, with proven convergence and improved performance in non-convex problems.

Contribution

It proposes a new double backtracking procedure controlling both step size and extrapolation in inertial Bregman proximal methods, with theoretical convergence guarantees.

Findings

Convergence to critical points is proven for the proposed algorithms.

Numerical experiments show improved performance in image processing and machine learning tasks.

The method effectively balances inertial step and backtracking for non-convex optimization.

Abstract

Backtracking line-search is an old yet powerful strategy for finding a better step sizes to be used in proximal gradient algorithms. The main principle is to locally find a simple convex upper bound of the objective function, which in turn controls the step size that is used. In case of inertial proximal gradient algorithms, the situation becomes much more difficult and usually leads to very restrictive rules on the extrapolation parameter. In this paper, we show that the extrapolation parameter can be controlled by locally finding also a simple concave lower bound of the objective function. This gives rise to a double convex-concave backtracking procedure which allows for an adaptive choice of both the step size and extrapolation parameters. We apply this procedure to the class of inertial Bregman proximal gradient methods, and prove that any sequence generated by these algorithms converges globally to a critical point of the function at hand. Numerical experiments on a number of challenging non-convex problems in image processing and machine learning were conducted and show the power of combining inertial step and double backtracking strategy in achieving improved performances.

Equations

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

$$y^k$$

$$x^{k+1}$$

$$x^{k+1}={\mathrm{argmin}}_{x\in {}^d}\left\{g\left(y^k\right)+⟨\nabla g\left(y^k\right),x-y^k⟩+\frac{{\bar{L}}_k}{2}{\Vert x-y^k\Vert }^2\right\}.$$

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

$$-\frac{{\underline{L}}_k}{2}{\Vert x-y^k\Vert }^2\le g\left(x\right)-g\left(y^k\right)-⟨\nabla g\left(y^k\right),x-y^k⟩\le \frac{{\bar{L}}_k}{2}{\Vert x-y^k\Vert }^2,$$

$$\overline{\gamma }<\sqrt{\frac{\bar{L}}{\underline{L}+\bar{L}}}\left(=\frac{1}{\sqrt{2}}\text{for}\bar{L}=\underline{L}\right),$$

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

$$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⟩.$$

$$(\mathcal{P})\inf \left\{\Psi \equiv f\left(x\right)+g\left(x\right):x\in \overline{C}\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.$$

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

$${\mathrm{prox}}_f^h\left(x\right)\in \mathrm{argmin}\left\{f\left(u\right)+D_h\left(u,x\right):u\in {}^d\right\},\forall x\in \mathrm{int}\mathrm{dom}h.$$

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

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

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

$$y^k$$

$$x^{k+1}$$

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

$$y^k$$

$$\left(\delta -\epsilon \right)D_h\left(x^{k-1},x^k\right)\ge \left(1+{\underline{L}}_k{\tau }_{k-1}\right)D_h\left(x^k,y^k\right)$$

$$g\left(x^k\right)\ge g\left(y^k\right)+⟨\nabla g\left(y^k\right),x^k-y^k⟩-{\underline{L}}_k{D}_h\left(x^k,y^k\right).$$

$$x^{k+1}\in {\mathrm{argmin}}_u\left\{f\left(u\right)+⟨\nabla g\left(y^k\right),u-y^k⟩+\frac{1}{{\tau }_k}{D}_h\left(u,y^k\right)\right\}$$

$$g\left(x^{k+1}\right)\le g\left(y^k\right)+⟨\nabla g\left(y^k\right),x^{k+1}-y^k⟩+{\bar{L}}_k{D}_h\left(x^{k+1},y^k\right).$$

$$0\le {\gamma }_k\le \overline{\gamma },\overline{\gamma }<\sqrt{\frac{{\bar{L}}_{k-1}}{{\bar{L}}_{k-1}+{\underline{L}}_k}},$$

$$\alpha \left(h\right):=\inf \left\{\frac{D_h\left(x,y\right)}{D_h\left(y,x\right)}:x,y\in \mathrm{int}\mathrm{dom}h,x\ne y\right\}\in \left[0,1\right].$$

$$\alpha \left(h\right)D_h\left(x,y\right)\le D_h\left(y,x\right)\le \alpha {\left(h\right)}^{-1}{D}_h\left(x,y\right),$$

$$\left(\delta -\epsilon \right)D_h\left(x^{k-1},x^k\right)\ge \left(1+{\underline{L}}_k{\tau }_{k-1}\right)D_h\left(x^k,y^k\right),$$

$$D_h\left(x_1,y\right)\le \kappa D_h\left(x_2,x_1\right),\forall \gamma \in \left[0,{\gamma }^{\ast }\right].$$

$$D_h\left(y,x_2\right)$$

$$=D_h\left(y,x_1\right)+D_h\left(x_1,x_2\right)+\gamma ⟨\nabla h\left(x_1\right)-\nabla h\left(x_2\right),x_1-x_2⟩$$

$$=D_h\left(y,x_1\right)+D_h\left(x_1,x_2\right)+\gamma \left(D_h\left(x_1,x_2\right)+D_h\left(x_2,x_1\right)\right).$$

$$D_h\left(y,x_2\right)\le \frac{1}{\alpha \left(h\right)}\left[D_h\left(x_1,y\right)+\left(\gamma \alpha \left(h\right)+1+\gamma \right)D_h\left(x_2,x_1\right)\right].$$

Full text

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

Keywords: Composite minimization, proximal gradient algorithms, inertial methods, convex-concave backtracking, non-Euclidean distances, Bregman distance, global convergence, Kurdyka-Łojasiewicz property.

1 Introduction

In this work we are interested in tackling non-convex additive composite minimization problems, which include the sum of two extended-valued functions: a non-smooth function denoted by $f$ (possibly non-convex) and a smooth function denoted by $g$ (possibly non-convex). More precisely, we consider problems of the following form

[TABLE]

where $\overline{C}$ is a nonempty, closed and convex set in d. We will give a more precise statement in Section 2 about the involved functions and set. There is a tremendous number of applications in machine learning, computer vision, statistics, and many more, that can be formulated in this framework.

Motivated by challenging applications as illustrated in Section 6, we consider here an instance of problem $(\mathcal{P})$, where the smooth function $g$ has a gradient that is not necessarily globally Lipschitz continuous. The restrictive assumption of having Lipschitz continuous gradient can be replaced with a certain convexity condition, which was proposed and developed first in [5] for problems $(\mathcal{P})$ with convex functions, and recently extended to the non-convex setting in [12]. More details on these recent developments will be given below in Section 2.

This convexity condition easily yields an approximation of the objective function at hand by a convex function from above (majorant) and a concave function from below (minorant). In the traditional setting, where the gradient of the smooth function $g$ is Lipschitz continuous, the majorant and the minorant are quadratic functions. In this case, it is well-known that the tightness of the quadratic approximations is directly related to restrictions on the step size to be used in the algorithm. The same relation is true for the convexity condition. In addition to their global existence, these approximations can be locally improved by backtracking (line search) strategies and it is well-known that tight approximations are advantageous, as we explain below in more detail.

Interestingly, while the step size is usually restricted by the quality of the majorant, the extrapolation (also known as inertia or over-relaxation) parameter is also affected by the quality of the minorant. This observation suggests to adapt the majorant and the minorant independently. In this paper we propose an efficient backtracking strategy that locally determines a tight majorant and minorant to exploit as much information as possible from the objective function, to be used in the proposed algorithm. This leads to a highly efficient algorithm, which is able to detect “the degree of local convexity” of the objective function (see Section 3 for details). As the backtracking procedure seeks for tight convex majorants and concave minorants, our idea is to combine it with an inertial step. We propose an inertial version of the Bregman Proximal Gradient (BPG) algorithm, which uses a convex-concave backtracking procedure to dynamically adjust the step size and the extrapolation parameter. Therefore, we call our algorithm Convex-Concave Inertial BPG (CoCaIn BPG in short). We prove a global convergence result of this algorithm (see Section 3.2 for an overview of the results and Section 5 for the details) to critical points of the objective function. The efficiency, which we demonstrate on several practical applications, comes from combining the inertial step with the novel convex-concave backtracking strategy, which fully exploits the power of tight local approximations in achieving large step sizes and large extrapolation parameters that can be used at the same time.

Before concluding this section, we would like to give the reader a first intuition about the convex-concave backtracking strategy on a simple instance of problem $(\mathcal{P})$.

A simple illustrative example.

In the following, we consider the following particular instance of problem $(\mathcal{P})$: $C={}^{d}$, $f\equiv 0$ and the gradient of $g$ is $L$-Lipschitz continuous. Even in this simpler setting, the convex-concave backtracking strategy is novel.

In this smooth and non-convex setting, an update step of a classical inertial based gradient method, starting with some $x^{0}\in{}^{d}$, reads as follows

[TABLE]

where $\gamma_{k}\in\left[0,1\right]$, $k\in\mathbb{N}$, is an extrapolation parameter and $\bar{L}_{k}>0$. If $g$ is convex and the extrapolation parameter $\gamma_{k}$ is carefully chosen, this recovers the popular Nesterov Accelerated Gradient method [32] (for $f\neq 0$, again in the convex setting, see [7]). It is well-known that the gradient step above, can be equivalently written as follows

[TABLE]

For a proper $\bar{L}_{k}$, the function to be minimized above is a convex quadratic majorant of the function $g$ (due to the classical Descent Lemma), which is a property that is also crucial for the convergence analysis of the algorithm. Classically, $\bar{L}_{k}\geq L$, $k\in\mathbb{N}$, is a sufficient condition to guarantee the existence of a quadratic majorant. However, locally, i.e., between the points $y^{k}$ and $x^{k+1}$, the parameter $\bar{L}_{k}$ may be significantly smaller than the global Lipschitz constant $L$ (which will immediately affect the step size of the algorithm). More precisely, note that the Descent Lemma,

[TABLE]

actually guarantees the existence of a quadratic minorant and a quadratic majorant that are determined by the same (global) parameter $L$. However, only the majorant limits the step size that is used in the algorithm. As shown in Figure 1, tighter approximations can be computed if the parameters of the minorant and the majorant are allowed to differ:

[TABLE]

i.e., the minorant parameter $\underline{L}_{k}$ could be different from the majorant parameter $\bar{L}_{k}$.

While the step size of the algorithm only depends on the majorant parameter $\bar{L}_{k}$, the extrapolation parameter $\gamma_{k}$ also depends on the minorant parameter $\underline{L}_{k}$. When $\bar{L}_{k}=\bar{L}$ and $\underline{L}_{k}=\underline{L}$, for all $k\in\mathbb{N}$, it was established in [46] that for any $0\leq\gamma_{k}\leq\overline{\gamma}$, when

[TABLE]

the generated sequence converges linearly (under certain error bound condition).

If the minorant parameter $\underline{L}_{k}$ is close to [math], which means that the function $g$ is “locally convex”, the extrapolation parameter $\gamma_{k}$ can be taken close to $1$, which makes the algorithm we present “similar” to an Accelerated Gradient method in the non-convex setting.

Below, we will show that using the minorant and the majorant in a local fashion (instead of their global counterparts) is very useful in developing the inertial Bregman Proximal Gradient method.

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

2 The Bregman Framework

In this section we will first recall the definition of Bregman distance, which stands at the heart of our developments. It was introduced in [14] and popularized by [17]. Based on that we will shortly review the recent concept of smooth adaptable functions, which in some sense extends and generalizes the class of smooth functions with globally Lipschitz continuous gradient. Then, we will provide the basic and essential ingredients to deal with the Bregman Proximal Gradient method.

We begin with the notion of kernel generating distance functions, which was recently stated in [12] (in this respect see also [4]).

Definition 2.1**.**

(Kernel Generating Distance) Let $C$ be a nonempty, convex and open subset of d. Associated with $C$, a function $h:{}^{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)$, the Bregman distance that is associated to $h$, is a proximity measure $D_{h}:\mathrm{dom}\,{h}\times\mathrm{int\,dom}\,{h}\to{}_{+}$ which is defined by

[TABLE]

This object is not a distance according to the classical definition (for example, it is not symmetric in general). However, the Bregman distance between two points is nonnegative if and only if the function $h$ is convex. If $h$ is known to be strictly convex, we have that $D_{h}\left(x,y\right)=0$ if and only if $x=y$. The classic example of a Bregman distance is the squared Euclidean distance, which is generated by $h(x)=\left\|{x}\right\|^{2}$. For more examples, results and applications of Bregman distances, see [18, 43, 22, 6, 44] and references therein.

An important property that is always crucial when dealing with Bregman distances is the well-known three-points identity [20, Lemma 3.1]: for any $y,z\in\mathrm{int\,dom}\,{h}$ and $x\in\mathrm{dom}\,{h}$,

[TABLE]

We conclude this part by restating our optimization model

[TABLE]

and making the first connection to the Bregman framework. One important feature of using Bregman distances in optimization algorithms is the ability of relate the constraint set $C$ to a certain kernel generating distances function $h\in\mathcal{G}(C)$. From now on, we make the following assumption.

Assumption A**.**

• $\rm{(i)}$

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

• $\rm{(ii)}$

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

• $\rm{(iii)}$

$g:{}^{d}\rightarrow\left(-\infty,+\infty\right]$ is a proper and lower semicontinuous function (possibly non-convex) with $\mathrm{dom}\,{h}\subset\mathrm{dom}\,{g}$, which is continuously differentiable on $C$.

• $\rm{(iv)}$

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

2.1 Smooth Adaptable Functions

One goal of this work is to deal with the non-convex optimization model ($\mathcal{P}$) where the gradient of the smooth function $g$ is not globally Lipschitz. Recently, Bauschke, Bolte and Teboulle [5], observed that the property of having a Lipschitz continuous gradient can be interpreted equivalently as a certain convexity condition on the function itself. This opens the gate for generalizing known results in the convex setting. It was extended to the non-convex setting in [12] with the concept of smooth adaptable functions given below.

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

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

The convexity requirement of $Lh+g$ can be written with respect to a different parameter $\ell\leq L$, which is key to the proposed double backtracking procedure to be developed in Section 3.1. In this section, for the sake of simplicity, we use $\ell=L$.

The optimization model ($\mathcal{P}$) appears with a smooth term in the objective function which is very common in many fields of applications. A crucial pillar in designing and analyzing algorithms for tackling this model, is usually based on the fact that the smooth part in the objective function has a Lipschitz continuous gradient. This property, via the well-known Descent Lemma, guarantees us that a lower and an upper quadratic approximation exist. For $L$-smooth adaptable functions, we will use the following extended version of the Descent Lemma (see [12, Lemma 2.1, p. 2134]).

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

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

[TABLE]

Remark 2.1** (Invariance to Strong Convexity).**

We would like to note that the $L$-smooth adaptable property is invariant when $h$ is additionaly assumed to be $\sigma$-strongly convex. Indeed, as described in [12], since convexity of $g$ is not needed, we can define $\omega(x):=\left(\sigma_{1}/2\right)\left\|{x}\right\|^{2}$, and then for any $0<\sigma_{1}<\sigma$, we have

[TABLE]

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

2.2 The Bregman Proximal Gradient Algorithm

In this section we review the basic notations and results needed to study Bregman based optimization methods. We first recall the definition of the Bregman proximal mapping [43], which is associated with a proper and lower semi-continuous function $f:{}^{d}\rightarrow\left(-\infty,+\infty\right]$, and is defined by

[TABLE]

With $h\equiv\left(1/2\right)\left\|{\cdot}\right\|^{2}$, the above boils down to the classical set-valued Moreau proximal mapping introduced in [29]. We refer the reader to the recent survey paper [44], and references therein. Here, we will focus on the Bregman proximal gradient mapping, which will take a central role in the algorithm to be developed in the next section. Given $x\in\mathrm{int\,dom}\,{h}$ and a step size parameter $\tau>0$, the Bregman proximal gradient mapping is defined by

[TABLE]

where the second equality follows from the fact that $\mathrm{dom}\,{h}\subset\overline{C}$. Note that here with $h\equiv\left(1/2\right)\left\|{\cdot}\right\|^{2}$, the above recovers the classical proximal gradient mapping. Since $f$ could be non-convex, the mapping $T_{\tau}$ is not, in general, single-valued. This mapping emerges from the usual approach, which consists of linearizing the differentiable function $g$ around a point $x$ and regularizing it with a proximal distance from that point. Similar to [12], the following assumption guarantees that the Bregman proximal gradient mapping is well-defined.

Assumption B**.**

• $\rm{(i)}$

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

[TABLE]

• $\rm{(ii)}$

For all $x\in C$, we have $T_{\tau}\left(x\right)\subset C$.

Assumption B(i) is a standard coercivity condition, which is for instance automatically satisfied when $\overline{C}$ is compact. On the other hand, Assumption B(ii) can be shown to hold under a classical constraint qualification condition. It also holds automatically when $f$ is convex or when $C={}^{d}$. The following result from [12], ensures that the Bregman proximal gradient mapping is well-defined.

Lemma 2.2** (Well-Posedness of $T_{\tau}$).**

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

3 The Inertial Bregman Proximal Gradient Method

Our proposed algorithm belongs to the class of inertial based optimization methods. The most well-known method in this class is the so-called Heavy-ball method, which was introduced by Polyak [40] to minimize convex and smooth functions. A popular variant of the method, when applied to the additive composite model ($\mathcal{P}$) with $C={}^{d}$, takes the following form. Start with any $x^{0}=x^{1}\in{}^{d}$, and generate iteratively a sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ via

[TABLE]

where $\gamma_{k}\in\left[0,1\right]$ is an extrapolation parameter and $\tau_{k}>0$ is a step size paramter. In [37], an inertial proximal gradient algorithm, called iPiano, was proposed555With a small modification that the proximity term is centered around the extrapolated point $y^{k}$, while the gradient of $g$ is evaluated at $x^{k}$.. It was shown that under Assumption A, if $f$ is convex and $g$ has a globally Lipschitz continuous gradient, the sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ converges globally to a critical point (in this setting, under additional error-bound condition, a linear rate of convergence was proved in [46]). The case where also the function $f$ is not necessarily convex was treated in [13, 34]. Two years later, in [39] a block version of the method, called iPALM was proposed and analyzed in the fully non-convex setting, i.e., both $f$ and $g$ are non-convex. In this case, a global convergence result to critical points was also established. A unified analysis was presented in [36].

In this work we propose a Bregman variant of the method mentioned above (see steps (3.1) and (3.2)), which also handle the two involved parameter $\gamma_{k}$ and $\tau_{k}$, $k\in\mathbb{N}$, in a dynamic fashion. To this end we incorporate into our basic steps two routines aiming at controlling and updating these parameters.

3.1 The Convex-Concave Backtracking Procedure

As we already illustrated on a simple example in the introduction, the origin of this procedure comes from the fact that for smooth adaptable functions we can build lower and upper approximations as given in Lemma 2.1:

[TABLE]

Even though the existence of the parameters $\underline{L}$ and $\bar{L}$ could be globally guaranteed, in practice it is often difficult or computationally expensive to evaluate them. In such cases it is recommended to apply a backtracking procedure that can locally verify the validity of the inequalities given in (3.3). However, in most cases only the upper approximation and the corresponding parameter $\bar{L}$ are used. Here, we will develop a double backtracking procedure that locally verifies both the lower and the upper approximations, in order to better control and update the extrapolation parameter $\gamma_{k}$ and the step size parameter $\tau_{k}$ at each iteration $k\in\mathbb{N}$. To the best of our knowledge, this is the first attempt to use the lower approximation in algorithms for tackling non-convex problems. It should be noted that in the case that $g$ is convex we have by definition $\underline{L}=0$, or even a convex quadratic lower approximation can be found when $g$ is strongly convex (see [44] for a discussion and references about a strong convexity property with respect to a Bregman distance). Based on the concepts described above, we will make the following additional assumptions on the involved functions.

Assumption C**.**

• $\rm{(i)}$

The function $h:{}^{d}\rightarrow\left(-\infty,+\infty\right]$ is $\sigma$-strongly convex on $C$.

• $\rm{(ii)}$

The pair of functions $(g,h)$ is $L$-smooth adaptable on $C$.

• $\rm{(iii)}$

There exists $\alpha\in\real$ such that $f\left(\cdot\right)-\left(\alpha/2\right)\left\|{\cdot}\right\|^{2}$ is convex666Such functions are called semi-convex with modulus $\alpha$ (see [34, 35])..

A few comments on the assumption above are now in order. The first item is related to Remark 2.1, which says that the smooth adaptable property is invariant to strongly convex kernel generating distance functions $h$. The third assumption allows us to deal with non-convex functions $f$ since $\alpha$ could be negative. See Section 6 for examples of functions that satisfy all these assumptions. Now we are ready to present our algorithm, which is called Convex-Concave Inertial (CoCaIn) Bregman Proximal Gradient.

The two input parameters $\delta$ and $\varepsilon$ are free to be chosen by the user. As we will see later the parameter $\varepsilon$ measures the descent to be achieved at each iteration of the algorithm.

The steps (3.4) and (3.7) are the classical steps of the inertial proximal gradient method, while here since we are dealing with the Bregman variant, it must be guaranteed that the auxiliary vector $y^{k}$ as defined in (3.4) belongs to $\mathrm{int\,dom}\,{h}$. Otherwise the Bregman proximal gradient step (3.7) is not defined (see Section 2.2). Even though, in general, it is not easy to guarantee that, in our case this will not be an issue. Indeed, in order to derive global convergence results of Bregman based algorithms in the non-convex setting an essential assumption seems to be that the kernel generating distance function $h$ has a full domain, i.e., $\mathrm{dom}\,{h}={}^{d}$ (see, for instance, [12] for more details about this limitation). The steps (3.6) and (3.8) implement the double backtracking procedure (see Section 5.4). The step (3.5) is designed to control the extrapolation parameter $\gamma_{k}$, $k\in\mathbb{N}$, and should be validated at each iteration. However, a natural question would be if such a parameter always exists? We postpone the positive answer to this question, to Section 4, and conclude this section with a list of our theoretical contributions.

3.2 Summary of the Convergence Results

Before we proceed with the well-posedness of CoCaIn BPG and the convergence analysis, we provide here a brief summary of our results.

• •

We show the well-posedness of CoCaIn BPG, in the sense that, one can always find $\gamma_{k}$ such that (3.5) is satisfied for all $k\in\mathbb{N}$ (see Lemma 4.1). Moreover, we show that it suffices to know the Bregman symmetric coefficient $\alpha\left(h\right)$ (Definition 4.1), in order to estimate the extrapolation parameter $\gamma_{k}$, $k\in\mathbb{N}$.

• •

In the Euclidean setting, i.e., when $h=\left(1/2\right)\left\|{\cdot}\right\|^{2}$, we provide an explicit formula for the maximal extrapolation parameter

[TABLE]

which uses the majorant parameter $\bar{L}_{k-1}$ from the previous iterate, which is a key for the efficient implementation of the proposed convex-concave backtracking procedure. When $\bar{L}_{k-1}=\underline{L}_{k}$, we easily recover that $\overline{\gamma}<1/\sqrt{2}$.

• •

Stability and convergence of the objective function values of CoCaIn BPG, which relies on finding an appropriate sequence of Lyapunov functions that enjoys a sufficient descent property (see Proposition 5.1).

• •

Global convergence of a sequence generated by the CoCaIn BPG method to critical points of the objective function $\Psi$ (see Theorem 5.2). This result relies on the concept of Gradient-like Descent Sequences (see Definition 5.1 below).

4 Well-Posedness of CoCaIn BPG

Now, we would like to verify the well-posedness of the CoCaIn BPG algorithm. An important tool in achieving our goal is the recently introduced symmetry coefficient of a Bregman distance, which measures the lack of symmetry in $D_{h}\left(\cdot,\cdot\right)$, see [5].

Definition 4.1** (Symmetry Coefficient).**

Given $h\in\mathcal{G}(C)$, its symmetry coefficient is defined by

[TABLE]

An important and immediate consequence of this definition is the fact that for all $x,y\in\mathrm{int\,dom}\,{h}$ we have

[TABLE]

where we have adopted the convention that $0^{-1}=+\infty$ and $+\infty\times r=+\infty$ for all $r\geq 0$. Clearly, the closer is $\alpha\left(h\right)$ to $1$, the more symmetric $D_{h}$ is with perfect symmetry when $\alpha\left(h\right)=1$ (which holds if and only if $h=\left\|{\cdot}\right\|^{2}$).

To this end, we need to convince the reader about the existence of $\gamma_{k}$, $k\in\mathbb{N}$, which satisfies (3.5), i.e., that

[TABLE]

holds true. The following result provides a positive answer to the existences question and information on the relevant extrapolation parameters that satisfy this inequality.

Lemma 4.1** (General Extrapolation Behavior).**

Given $h\in\mathcal{G}(C)$ with $\alpha\left(h\right)>0$. Let $x_{1},x_{2},y\in\mathrm{int\,dom}\,{h}$ and $y:=x_{1}+\gamma\left(x_{1}-x_{2}\right)$ with $\gamma\geq 0$. Then, for a given $\kappa>0$, there exists $\gamma^{\ast}>0$ such that

[TABLE]

Proof.

From the three points identity (see (2.1)) we have

[TABLE]

Now, from (4.1), we obtain that

[TABLE]

On the other hand, since $x_{1}=\left(y+\gamma x_{2}\right)/\left(1+\gamma\right)$, we can use the fact that $u\rightarrow D_{h}\left(u,v\right)$, for a fixed $v\in\mathrm{int\,dom}\,{h}$, is a convex function and therefore

[TABLE]

where the last inequality follows from (4.1). By combining the last two inequalities we derive that

[TABLE]

and, by re-arranging we have

[TABLE]

First, it is easy to verify that for $\gamma<\alpha\left(h\right)^{2}/\left(1-\alpha\left(h\right)^{2}\right)$, the denominator is positive. In addition, to find $\gamma$ such that

[TABLE]

we will use simple algebraic manipulations. Indeed, by re-arranging we have

[TABLE]

Since $\alpha\left(h\right)^{2}\leq 1$, it follows that $b>0$. We also have that $\Delta=b^{2}+4a\alpha\left(h\right)^{2}\kappa>0$, and thus there exists a positive root denoted by $\gamma^{\ast}$. Therefore, for any $\gamma\in\left[0,\gamma^{\ast}\right]$, the desired result follows. ∎

Remark 4.1**.**

Note that in the above lemma, $\gamma^{\ast}$ depends only on the symmetry coefficient $\alpha\left(h\right)$. Therefore, for the Euclidean distance with $\alpha\left(h\right)=1$, this implies that,

[TABLE]

However, for the Euclidean distance, the expression in (4.2), can be simplified significantly. Indeed, since we take $h=(1/2)\left\|{\cdot}\right\|^{2}$, then using the fact that $y^{k}-x^{k}=\gamma_{k}\left(x^{k}-x^{k-1}\right)$ we obtain that $\gamma_{k}\leq\sqrt{\kappa}$. In the case of CoCaIn BPG, we have the following restriction on the maximal extrapolation parameter that can be used

[TABLE]

A related bound also appeared in [46] as we discussed in the introduction. When, the values of $\underline{L}_{k}$ and $\bar{L}_{k-1}$ are almost equal and $\delta-\varepsilon\approx 1$, then it is possible to choose the inertial parameter $\gamma_{k}$ such that $\gamma_{k}\approx 1/\sqrt{2}$. We discuss more about bounds of $\gamma_{k}$, $k\in\mathbb{N}$, in Section 5.3.

5 Convergence Analysis of CoCaIn BPG

Before we proceed to the convergence analysis, we need the following technical lemma.

Lemma 5.1** (Function Descent Property).**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a sequence generated by CoCaIn BPG. Then, for all $k\in\mathbb{N}$, we have

[TABLE]

Proof.

Fix $k\geq 1$. From the convexity of $f\left(\cdot\right)-\left(\alpha/2\right)\left\|{\cdot}\right\|^{2}$, which holds thanks to Assumption C(iii), we obtain from the sub-gradient inequality [41, Example 8.8 and Proposition 8.12] that

[TABLE]

where $\xi^{k+1}\in\partial f\left(x^{k+1}\right)$. By rearranging the inequality we obtain

[TABLE]

From the optimality condition of step (3.7), we have that

[TABLE]

which combined with (5.1) yields that

[TABLE]

where the last equality follows from the three-points identity (see (2.1)). On the other hand, using the lower approximation given in (3.6) and the upper approximation given in (3.8), we have that

[TABLE]

Combining the last two inequalities and using the fact that $\tau_{k}^{-1}\geq{\bar{L}_{k}}$, implies that

[TABLE]

which completes the proof. ∎

Since we are dealing with inertial based methods, which belong to the class of non-descent methods, we can not expect to use classical convergence techniques for non-convex problems (see below for more information about it). In order to overcome the lack of descent, we will use the Lyapunov technique, which involves the construction of a sequence of new functions, which will be used to “better” measure the progress of the algorithm, where by progress we mean a decrement in the Lyapunov function values. In several cases a trivial Lyapunov function would be to use the function itself, however in the case of non-descent methods, it is not a good choice, since it does not capture well the behavior of the iterates. The behavior of two subsequent iterates must be taken into consideration along with the function, as observed in [37, 42].

5.1 Lyapunov Function Descent Property of CoCaIn BPG

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a sequence generated by CoCaIn BPG. We define, at iterate $k\in\mathbb{N}$, the following Lyapunov function

[TABLE]

This Lyapunov function involves two terms: (i) the term $\tau_{k-1}\left(\Psi\left(x^{k}\right)-v(\mathcal{P})\right)$, which measures the progress in original function values $\Psi$ with respect to the global optimal value of problem ($\mathcal{P}$) and (ii) the term given by $\delta D_{h}\left(x^{k-1},x^{k}\right)$, which ensures that the iterates stay close enough, with respect to the Bregman distance. Before we motivate further the usage of this Lyapunov function, we show its descent property.

Proposition 5.1**.**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a sequence generated by CoCaIn BPG. Then, for all $k\in\mathbb{N}$, we have

[TABLE]

Proof.

Multiplying (5.1) with $\tau_{k}$, we obtain

[TABLE]

By the definition of the Lyapunov function $\Phi_{\delta}^{k}$ and the fact that $\tau_{k}\leq\tau_{k-1}$ we have

[TABLE]

With $1-\delta>0$ and the strong convexity of $h\left(\cdot\right)$, that follows from Assumption C(i), we obtain

[TABLE]

where the last inequality holds, since $\tau_{k}^{-1}\geq{\bar{L}}_{k}$ and ${\bar{L}}_{k}\geq-\alpha/\left(1-\delta\right)\sigma$. Next, we observe that

[TABLE]

where the first inequality is due to the step (3.5) of the algorithm and the second inequality is due to fact that $\tau_{k}\leq\tau_{k-1}$. By rearranging we obtain,

[TABLE]

thus completing the proof. ∎

Proposition 5.2**.**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a sequence generated by CoCaIn BPG. Then, the following assertions hold:

• $\rm{(i)}$

The sequence $\left\{\Phi_{\delta}^{k+1}\left(x^{k+1},x^{k}\right)\right\}_{k\in\mathbb{N}}$ is nonincreasing.

• $\rm{(ii)}$

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

• $\rm{(iii)}$

$\min_{1\leq k\leq n}D_{h}\left(x^{k-1},x^{k}\right)\leq\Phi_{\delta}^{1}\left(x^{1},x^{0}\right)/\left(\varepsilon n\right)$.

Proof.

• $\rm{(i)}$

This follows trivially from Proposition 5.1, since $\varepsilon>0$.

• $\rm{(ii)}$

Let $n$ be a positive integer. Summing (5.4) from $k=1$ to $n$ we get

[TABLE]

since $\Phi_{\delta}^{n+1}\left(x^{n+1},x^{n}\right)\geq 0$. Taking the limit as $n\rightarrow\infty$, we obtain the first desired assertion, from which we immediately deduce that $\left\{D_{h}\left(x^{k-1},x^{k}\right)\right\}_{k\in\mathbb{N}}$ converges to zero.

• $\rm{(iii)}$

From (5.5) we also obtain,

[TABLE]

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

∎

In order to proceed with the global convergence analysis of CoCaIn BPG, we will need throughout the rest of this section, to additionally assume the following.

Assumption D**.**

• $\rm{(i)}$

$\mathrm{dom}\,{h}={}^{d}$.

• $\rm{(ii)}$

$\nabla h$ and $\nabla g$ are Lipschitz continuous on any bounded subset of d.

5.2 Global Convergence for CoCaIn BPG

In this subsection we show the global convergence result of CoCaIn BPG. The goal is to show that the whole sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$, that is generated by CoCaIn BPG, converges to a critical point. To this end, we denote the set of critical points by

[TABLE]

Note that, such a set is well-defined due to Fermat’s rule [41, Theorem 10.1, p. 422] and due to the concept of limiting subdifferential.

From now on we will make the following assumption regarding the sequence of majorant parameters $\left\{\bar{L}_{k}\right\}_{k\in\mathbb{N}}$: there exists an integer $K\in\mathbb{N}$ such that $\bar{L}_{k}=\bar{L}$ for all $k\geq K$ ($K$ can be as large as the user wishes). It should be noted that thanks to Assumption C(ii) and Lemma 2.1, there exists a global majorant parameter $\bar{L}$ such that (3.8) holds true for all $k\in\mathbb{N}$. On the other hand, since in anyway we require that the parameters do not decrease between two successive iterations, it makes sense that at some point we will stop changing them and continue with a fixed value. However, it is very important not using the global parameter $\bar{L}$ right from the beginning since in practice the parameter $\bar{L}_{k}$ determined by (3.8) might be much smaller (especially in early stages of the algorithm).

In the second phase of the algorithm, i.e., when $k\geq K$, it also makes sense to assume that $\tau_{k}=\tau$ for all $k\geq K$ where $\tau\leq\bar{L}^{-1}$. This immediately suggests that our Lyapunov function can also be simplified. More precisely, we define the following new Lyapunov function:

[TABLE]

where $\delta_{1}=\delta/\tau$.

The global convergence result is based on showing that CoCaIn BPG generates a gradient-like descent sequence according to Definition 5.1 (see below). This involves three properties which need to be verified: “sufficient descent condition”, “relative error condition” and “continuity condition”. Such a convergence analysis is based on a recent technique, which was initiated by Attouch and Bolte [1], and later on was simplified and unified in [11]. A more general framework was proposed in [36].

The main tool that stands behind this technique is the Kurdyka-Łojasiewicz (KL) property [26, 27] (see [8] for the non-smooth case), which is properly defined in the appendix. This property has been used in several recent works that deal with non-convex optimization problems (see [1, 3, 11] for early foundational works). For more details and information on the KL property, we refer the reader to the following papers [8, 1, 10, 2, 3, 11, 36] and references therein.

Verifying that a given function satisfies the KL property could be difficult, however in their seminal work [8], Bolte, Daniilidis and Lewis prove that any proper, lower semicontinuous and semi-algebraic function satisfies the KL property on its domain. This important result makes this proof technique very powerful, since we are familiar with many semi-algebraic functions that appear very often in applications. In fact, the same result holds for (possibly non-smooth) functions that are definable in an o-minimal structure [8, 9]. For examples and more details about the relations between KL and other important notions, see [8, 10] and references therein.

In order to derive the global convergence of our algorithm we follow this proof technique that we shortly recall now. For the interested readers we refer to [12, Appendix 6, p. 2147], where a short and self-contained summary of this proof methodology can be found. It should be noted again that here we consider a modification, which fits non-descent methods like CoCaIn BPG.

Definition 5.1** (Gradient-like Descent Sequence).**

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

• $\rm{(C1)}$

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

[TABLE]

• $\rm{(C2)}$

Relative error condition. There exist an integer $K\in\mathbb{N}$ and a positive scalar $\rho_{2}$ such that

[TABLE]

• $\rm{(C3)}$

Continuity condition. 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}}\Psi\left(x^{k}\right)\leq\Psi\left(\overline{x}\right)$.

Based on Definition 5.1 and the KL property, the following global convergence result holds true. We provide its proof in the appendix.

Theorem 5.1** (Global Convergence).**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a bounded gradient-like descent sequence for minimizing $\Psi_{\delta_{1}}$. If $\Psi$ 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}\,{\Psi}$.

Now, in a sequence of lemmas, we prove that CoCaIn BPG generates a gradient-like descent sequence for minimizing $\Psi_{\delta_{1}}$. In order to prove condition (C1), we first note that Proposition 5.2 is also valid for the new Lyapunov function $\Psi_{\delta_{1}}$ as recorded now (for the sake of simplicity we omit the exact details of the proof, which is almost identical to the proof above).

Proposition 5.3**.**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a sequence generated by CoCaIn BPG. Then, the following assertions hold:

• $\rm{(i)}$

The sequence $\left\{\Psi_{\delta_{1}}\left(x^{k+1},x^{k}\right)\right\}_{k\in\mathbb{N}}$ is nonincreasing, converging and condition (C1) of Definition 5.1 holds true.

• $\rm{(ii)}$

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

• $\rm{(iii)}$

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

Now we can prove the following result, which means that condition (C2) holds true.

Proposition 5.4**.**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a bounded sequence generated by CoCaIn BPG. Then, there exist $w^{k+1}\in\partial\Psi_{\delta_{1}}\left(x^{k+1},x^{k}\right)$ and a positive scalar $\rho_{2}$ such that

[TABLE]

Proof.

Fix $k\geq K$. By the definition of the Lyapunov function $\Psi_{\delta_{1}}\left(\cdot,\cdot\right)$ we obtain that

[TABLE]

Writing the optimality condition of the optimization problem which defines $x^{k+1}$ (see (3.7) and recall that for $k\geq K$, we have that $\tau_{k}=\tau$) yields that

[TABLE]

Therefore

[TABLE]

and by defining

[TABLE]

and $w_{2}^{k+1}\equiv\delta_{1}\left(\nabla h\left(x^{k}\right)-\nabla h\left(x^{k+1}\right)\right)$ we obviously obtain that $w^{k+1}\in\partial\Psi_{\delta_{1}}\left(x^{k+1},x^{k}\right)$ where $w^{k+1}=\left(w_{1}^{k+1},w_{2}^{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 d (see Assumption D(ii)), there exists $M>0$ such that

[TABLE]

where the last inequality follows also from the fact that $\left\|{\nabla^{2}h\left(x^{k+1}\right)}\right\|\leq M$, since $\nabla h$ is Lipschitz continuous on bounded subsets of d. Using step (3.4) we obtain that

[TABLE]

where we have used the fact that $\gamma_{k}\leq 1$, $k\in\mathbb{N}$. Since, we also have that

[TABLE]

the desired result is proved and condition (C2) also holds true. ∎

Now we are left with showing that CoCaIn BPG generates a sequence that satisfies condition (C3).

Proposition 5.5**.**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a bounded sequence generated by CoCaIn BPG. Let $x^{\ast}$ be a limit point of a subsequence $\left\{x^{k}\right\}_{k\in{\cal K}}$, then $\limsup_{k\in{\cal K}\subset\mathbb{N}}\Psi\left(x^{k}\right)\leq\Psi\left(x^{\ast}\right)$.

Proof.

Consider a subsequence $\left\{x^{n_{k}}\right\}_{k\in\mathbb{N}}$ which converges to $x^{\ast}$ (there exists such a subsequence since the sequence $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ is assumed to be bounded). Using Proposition 5.3(ii) and the strong convexity of $h\left(\cdot\right)$, 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}$. From the definition of $y^{k}$, see (3.4), it also follows that $\left\{y^{n_{k}-1}\right\}_{k\in\mathbb{N}}$ also converges to $x^{\ast}$. In addition, since $h$ is continuously differentiable on d we have that $\lim_{k\rightarrow\infty}D_{h}\left(x^{\ast},y^{{n_{k}}-1}\right)=0$. Now, from (3.7), it follows (after some simplifications), for all $k\geq K$, that

[TABLE]

Substituting $k$ by $n_{k}$ and letting $k\rightarrow\infty$, we obtain from the fact that $g$ is continuously differentiable on d, that

[TABLE]

Using this, and recalling that here $g$ is continuous, we obtain that $\limsup_{k\in{\cal K}\subset\mathbb{N}}\Psi\left(x^{n_{k}}\right)\leq\Psi\left(x^{\ast}\right)$, where ${\cal K}=\left\{n_{k}:\,k\geq K\right\}$. ∎

The global convergence of CoCaIn BPG now easily follows from our general result on gradient-like descent sequences (see Theorem 5.1)

Theorem 5.2** (Global Convergence of CoCaIn BPG).**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a bounded sequence generated by CoCaIn BPG. If $f$ and $g$ satisfy 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}\,{\Psi}$.

Before we conclude this section, we provide a simplified variant of CoCaIn BPG.

5.3 CoCaIn BPG Without Backtracking

Note that CoCaIn BPG uses a local estimate of the minorant and majorant parameters $\underline{L}_{k}$ and $\bar{L}_{k}$, $k\in\mathbb{N}$, determined by the backtracking steps (3.6) and (3.8), respectively. However, when the global parameter $L$ is known (guaranteed in Assumption C(ii)), we can skip the backtracking steps, and provide a simplified variant of CoCaIn BPG.

For the inertial step (5.9), when $h=\left(1/2\right)\left\|{\cdot}\right\|^{2}$ we can obtain that

[TABLE]

with $\bar{L}=\underline{L}$. Using Remark 4.1, if $\delta-\epsilon\approx 1$, one could choose the extrapolation parameter as follows $\gamma_{k}\approx 1/\sqrt{2}$. However, in general, the closed form expression for $\gamma_{k}$ is difficult to obtain, for which backtracking line-search strategy can be used. Recently, in [31] the authors showed a technique to obtain closed form inertia for general Bregman distances. We use their technique later in the context of Quadratic inverse problems to propose a new variant of CoCaIn BPG with closed form inertia.

5.4 Implementing the Double Backtracking Procedure

The update steps of CoCaIn BPG are based on the double backtracking strategy (see steps (3.6) and (3.8)). Here, we describe some implementation details of these two steps. Note that the inner loops for finding the minorant and the majorant parameters $\underline{L}_{k}$ and $\bar{L}_{k}$, $k\in\mathbb{N}$, are implemented in a sequential fashion. By this, we mean that at iteration $k\in\mathbb{N}$ we first execute the steps (3.4), (3.5) and (3.6) in order to compute an appropriate $y^{k}$, only then we proceed to steps (3.7) and (3.8) in order to compute $x^{k+1}$. Note that the fact that the sequence $\left\{\bar{L}_{k}\right\}_{k\in\mathbb{N}}$ does not decrease is crucial in order to decouple the steps (3.4) and (3.7). More precisely, we now describe the backtracking procedure to find $\underline{L}_{k}$. Let $\underline{\nu}>1$ be a scaling parameter and arbitrarily initialize $\underline{L}_{k,0}>0$. Then, we find the smallest $\underline{L}_{k}\in\left\{\underline{\nu}^{0}\underline{L}_{k,0},\underline{\nu}^{1}\underline{L}_{k,0},\underline{\nu}^{2}\underline{L}_{k,0},\ldots\right\}$ that satisfies (3.6) and such that $\gamma_{k}\geq 0$ satisfies

[TABLE]

We can now describe the procedure to find $\bar{L}_{k}$. Let $\overline{\nu}>1$ and initialize $\bar{L}_{k,0}:=\bar{L}_{k-1}$, then we take the smallest $\bar{L}_{k}\in\left\{\overline{\nu}^{0}\bar{L}_{k,0},\overline{\nu}^{1}\bar{L}_{k,0},\overline{\nu}^{2}\bar{L}_{k,0},\ldots\right\}$ that satisfies (3.8). Therefore, $\left\{\bar{L}_{k}\right\}_{k\in\mathbb{N}}$ is monotonically non-decreasing. Note, however, we do not require any monotonicity of the sequence $\left\{\underline{L}_{k}\right\}_{k\in\mathbb{N}}$.

The double backtracking strategy preserves the sign of $\underline{L}_{k}$, however, only $-\underline{L}_{k}\leq\bar{L}_{k}$ is required. Changing the sign of $\underline{L}_{k}$ when the function is locally strongly convex might lead to additional acceleration. However, we leave this kind of adaptation for future work.

6 Numerical Experiments

Our goal in this section is to illustrate the performance of CoCaIn BPG in various situations. We start with minimization of univariate functions, which emphasizes the power of incorporating inertial terms into the BPG algorithm and using the double backtracking procedure. Then we provide some insights on the following practical applications: Quadratic Inverse Problems in Phase Retrieval and Non-convex Robust Denoising with Non-convex Total Variation Regularization. More recently, the efficiency of CoCaIn BPG is also demonstrated in related work for Matrix Factorization [30] and Deep Linear Neural Networks [31].

6.1 Finding Global Minima of Univariate Functions

We begin with two examples of minimizing univariate non-convex functions, which shed some light on the two main features of our algorithm: (i) inertial term and (ii) double backtracking procedure. We consider unconstrained minimization of functions $g:\real\rightarrow\real$, with Lipschitz continuous gradient, i.e., model ($\mathcal{P}$) with $d=1$, $f\equiv 0$ and $C=\real$. The two functions are: $g\left(x\right)=\log\left(1+x^{2}\right)$ and $g\left(x\right)=\left(1+e^{x}\right)^{-1}$. We compare three methods: CoCaIn BPG with $h=\left(1/2\right)\left\|{\cdot}\right\|^{2}$ and refer to it as CoCaIn with Euclidean distance, classical Gradient Descent (GD) method with backtracking (which is actually CoCaIn with Euclidean distance and with $\gamma_{k}=0$ for all $k\in\mathbb{N}$), and iPiano777In this particular case, the method coincides with the Heavy-ball method [40]. [37] (with the inertial parameter set to $0.7$). When using a backtracking procedure in GD and iPiano methods, we mean that only the majorant parameter is varied. We use the same initialization for all the algorithms and report the performance in Figure 2.

In the second experiment, we illustrate the robustness of CoCaIn BPG to local minima and critical points. We consider the non-smooth and non-convex function $\Psi\left(x\right)=\left|x\right|+\sin\left(x\right)+\cos\left(x\right)$, with many critical points as shown in the center plot of Figure 3, and set $f\left(x\right)=\left|x\right|$ and $g\left(x\right)=\sin\left(x\right)+\cos\left(x\right)$ (which is obviously a non-convex function with Lipschitz continuous gradient). Here again we take $h=\left(1/2\right)\left\|{\cdot}\right\|^{2}$. In order to apply CoCaIn BPG, the main computational step is of the following form:

[TABLE]

which results in the following update step

[TABLE]

We compare CoCaIn BPG with Euclidean distance to the classical Proximal Gradient (PG) method with backtracking (CoCaIn BPG with Euclidean distance and $\gamma_{k}=0$, $k\in\mathbb{N}$), and iPiano. As mentioned in the first experiment, when using a backtracking procedure in PG and iPiano methods we mean that only the majorant parameter is varied.

As shown in Figure 3, CoCaIn BPG achieves the global minimum, whereas the PG with backtracking gets stuck in a local minimum. We performed the same experiment starting at $100$ equidistant points sampled from the interval $\left[-15,15\right]$. The average final function value for CoCaIn was $2.75$, whereas for PG method with backtracking it was $3.21$ and for the iPiano it was $3.37$. This means that CoCaIn BPG reaches the global minimum from $52$ points, PG method with backtracking achieves the global minimum only from $27$ points and iPiano from $39$ points. Hence, the behavior illustrated in Figure 3 is not due to the choice of initialization, but rather due to additional features of the CoCaIn BPG algorithm. This illustrates the great power of using double backtracking procedure in minimizing univariate non-convex functions.

6.2 Escaping Spurious Stationary Points

Here, we provide evidence that CoCaIn BPG can escape spurious stationary points in minimizing non-convex functions of two variables. Let $b_{i}\in\real$, $i=1,2,\ldots,m$, be samples of a noisy signal with additive Gaussian noise. A very common task in signal processing is to recover the true data. However, due to the noise, data can be prone to several outliers. In such cases, a robust loss [23] is used. Moreover, prior information about the data, can be embedded through a regularizing term (for instance, a sparsity promoting regularizer). Given $\lambda,\rho>0$, we consider minimization of

[TABLE]

with

[TABLE]

The function $f$ is a non-convex sparsity promoting regularizer (also known as the log-sum penalty term [16, 33]) and the function $g$ is a robust loss. For illustration purposes, we consider a simple instance of problem (6.3) where $m=2$, $\lambda=0.5$ and $\rho=100$. For minimizing this function we set $C={}^{2}$ and $h\left(x\right):=\left(1/2\right)\left(x_{1}^{2}+x_{2}^{2}\right)$ to be used in the CoCaIn BPG method.

Before presenting the numerical results, we would like to note that in this example, the function $f\left(x\right)-\left(\alpha/2\right)h\left(x\right)$ is convex for any $\alpha\leq-1$ and $Lh-g$ is convex for all $L\geq 100$. Each iteration of CoCaIn BPG would require to compute the Bregman proximal gradient mapping, which in this case reduces to the classical proximal gradient mapping (due to the choice of $h$). Note that due to the separability of the functions $f$ and $g$, the needed minimization problem can be split into two individual minimizations with respect to $x_{1}$ and $x_{2}$. These two optimization problems (after simple manipulations) reduces to computation of the proximal mapping of the univariate function $\log\left(1+\left|x\right|\right)$. A closed form formula can be found in [24] and reads as follows:

[TABLE]

where

[TABLE]

with $\left[x\right]_{+}:=\max\left\{0,x\right\}$.

Now we can apply CoCaIn BPG method and the function behavior is described in Figure 4.

The performance of CoCaIn BPG is illustrated in Figure 5, which shows that CoCaIn BPG can indeed escape spurious critical points to reach the global minimum.

6.3 Quadratic Inverse Problems in Phase retrieval

Phase retrieval has been an active research topic for several years [15, 45, 21, 28]. It gained a lot of attention from the optimization community, due to resulting hard non-convex problems [12, 21, 19]. The phase retrieval problem can be described as follows. Given sampling vectors $a_{i}\in{}^{d}$, $i=1,2,\ldots,m$, and measurements $b_{i}>0$, we seek to find a vector $x\in{}^{d}$ such that the following system of quadratic equations is approximately satisfied,

[TABLE]

One typical way to tackle this system is by solving an optimization problem that seeks to minimize a certain error/noise measure in accomodating the equations. The objective function also depends on the type of noise [19] in the system (for instance, Gaussian or Poisson noise). We assume additive Gaussian noise and the squared error measure

[TABLE]

with

[TABLE]

The function $f$ acts as a regularizing term and is used to incorporate certain prior information on the wished solution. We conduct experiments with two options of regularizing functions: (i) squared $\ell_{2}$-norm, $f\left(x\right)=\left(\lambda/2\right)\left\|{x}\right\|^{2}$ and (ii) $\ell_{1}$-norm, $f\left(x\right)=\lambda\left\|{x}\right\|_{1}$. When applying here the CoCaIn BPG method we use the following kernel generating distance function

[TABLE]

We obviously have that $\mathrm{dom}\,{h}={}^{d}$ and we record below a result [12, Lemma 5.1, p. 2143], which shows that the pair $\left(g,h\right)$ satisfies the L-smad property (see Definition 2.2).

Lemma 6.1**.**

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

[TABLE]

the function $Lh-g$ is convex on d.

By the design of CoCaIn BPG algorithm, the inertial parameter $\gamma_{k}$ must satisfy (3.5). However, this involves backtracking over $\gamma_{k}$, which can computationally expensive for high dimensional problems. To this regard, following [31], we propose closed form expression for $\gamma_{k}$ which satisfies (3.5). We also illustrate with our numerical experiments, that CoCaIn BPG variant with closed form inertia is competitive to our main algorithm CoCaIn BPG.

Lemma 6.2** (Closed form inertia).**

For $h$ defined in (6.6), we obtain the following gradient

[TABLE]

and for any $a\in\mathbb{R}^{d}$, we have

[TABLE]

Proof.

Consider the expansion at $x+a$ till second order terms, we thus have

[TABLE]

The first order terms result in (6.7) and we also have

[TABLE]

where the inequality follows due to Cauchy-Schwarz inequality. ∎

Lemma 6.3** ([31]).**

Let $h\in\mathcal{G}(C)$ be twice continuously differentiable on $C$. Then, the following identity holds

[TABLE]

Proposition 6.1**.**

Denote $\Delta_{k}:={x^{k}-x^{k-1}}$, for any $k\geq 1$ the following holds

[TABLE]

Proof.

We use the strategy from [31, Lemma 15]. From Lemma 6.3, we have

[TABLE]

where in the last step we used the upper bound (6.8) from Lemma 6.2. Also, we used the following inequality

[TABLE]

where in the last step we used $\gamma_{k}^{2}\leq 1$ and $(t_{1}+(1-t_{1})t)^{2}\leq 1$. With $\int_{0}^{1}(1-t)dt=\frac{1}{2}$ the result follows.∎

Therefore, in this case, Assumptions A, B, C and D are valid. We now discuss the update step of CoCaIn BPG, which requires the solution of the following subproblem

[TABLE]

Following [12], we provide closed form formulas for these optimization problems when $f$ is either the squared $\ell_{2}$-norm or the $\ell_{1}$-norm.

$\ell_{1}$-norm.

Here we use the following closed form solution, derived in [12, Proposition 5.1, p. 2145]. First, we define the soft-thresholding operator with respect to the parameter $\theta>0$, as follows

[TABLE]

where all operations are applied coordinate-wise. Then the closed form solution of problem (6.9) is given by

[TABLE]

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

[TABLE]

Squared $\ell_{2}$-norm.

Using similar arguments as of [12, Proposition 5.1, p. 2145], we can easily derive that the solution of problem (6.9) is given by

[TABLE]

where $t^{\ast}$ is the unique real root of the following cubic equation

[TABLE]

We illustrate, in Figure 6, the performance of CoCaIn BPG and CoCaIn BPG with closed form inertia (CoCaIn BPG CFI), compared with two other algorithms: (i) the Bregman Proximal Gradient Method with backtracking (denoted by BPG-WB) using the same kernel generating distance function (which is exactly CoCaIn BPG with $\gamma_{k}=0$ for all $k\in\mathbb{N}$) and (ii) the Inexact Bregman Proximal Minimization Line Search Algorithm (denoted by IBPM-LS) of [38]. We also compare with the Bregman Proximal Gradient (BPG) method of [12] without backtracking and with the parameter $L$ as derived in Lemma 6.1.

6.4 Non-convex Robust Denoising with Non-convex TV Regularization

We consider the problem of image denoising of a given image $b\in\mathbb{R}^{M\times N}$, where $M,N\in\mathbb{N}$. The goal is to obtain the true image, denoted by $x\in{}^{M\times N}$. However, in real world applications, it is possible that the measurements are noisy with outliers. The standard routine to deal with outliers is to use robust loss function. The basic idea is to heavily penalize small errors and reasonably penalize large errors. This is done to ensure that the predicted data $x$, is not influenced significantly by outliers. We consider a fully non-convex formulation of the problem, which includes a non-convex loss function along with a non-convex regularization.

We need the following technical details to provide the full problem statement. The spatial finite difference operator is given by

[TABLE]

where $i\in[M]$ and $j\in[N]$. The horizontal spatial finite differences are given by $(\mathcal{D}x)^{1}_{i,j}:=x_{i+1,j}-x_{i,j}$ for all $i<M$ and [math] otherwise. The vertical spatial finite differences are given by $(\mathcal{D}x)^{2}_{i,j}:=x_{i,j+1}-x_{i,j}$ for all $j<N$ and [math] otherwise.

The problem involves the following functions

[TABLE]

where $\lambda,\rho>0$. The function $f$ is non-smooth non-convex and $g$ is smooth non-convex. The function $g$ is a non-convex variant of the popular Total Variation (TV) regularizer, which is used to prefer smooth signals while preserving sharp changes in the signal (such as edges of images). For an overview on non-convex regularizations we refer the reader to [33, 47]. Consider $h\left(x\right)=\left(1/2\right)\left\|{x}\right\|^{2}_{F}$. It is easy to prove the convexity of $f\left(x\right)-\left(\alpha/2\right)\left\|{x}\right\|^{2}_{F}$, by checking that its right derivative is monotonically increasing [25, Theorem 6.4], for all $\alpha\leq-1$. The function $Lh-g$ is convex for $L\geq 16\lambda\rho$.

Due to separability of the function $f$, we can split the computation of the corresponding Bregman Proximal Gradient mapping, into the following separable subproblems

[TABLE]

which as discussed in Section 6.2, can be reduced to the computation of the proximal mapping of the function $\log\left(1+\left|x-b\right|\right)$.

We consider two additional experimental settings apart from our main setting given by (6.12) and (6.13). Firstly, we use the $\ell_{2}$-norm based data term with the same regularization as in (6.13). Secondly, we use the squared $\ell_{1}$-norm based data term with regularization as in (6.13). We use the good image given in Figure 7(a) and add severe noise randomly of $10^{5}$ magnitude. We illustrate the robustness of the model given by (6.12) and (6.13) to such outliers. The reconstructed image from $\ell_{2}$-norm based data penalty term is given in Figure 7(c) and the reconstructed image from $\ell_{1}$-norm based data penalty term is given in Figure 7(d), after applying CoCaIn BPG. Clearly the $\ell_{1}$-norm based data penalty is better than $\ell_{2}$-norm based data penalty term, which is due to the robustness properties of $\ell_{1}$-norm. However, even using $\ell_{1}$-norm is not enough in the presence of severe outliers, the robustness properties are not so significant. This is mitigated by our setting, where the reconstructed image is given in Figure 7(e). In our setting, the data term in (6.12) is very robust to outliers. In all the settings, we used $\lambda=10$ and $\rho=1$. The convergence plots for the experiments with (6.12) and (6.13) are given in Figure 7(f) and 7(g). Note that CoCaIn BPG CFI uses the closed form inertia with Euclidean distance. BPG-WB and BPG are same as in earlier experiments. IBPM-LS is a general purpose line-search algorithm for nonconvex nonsmooth problems proposed in [38]. Even though, IBPM-LS is general, BPG based methods are much faster. The comparisons also illustrate that CoCaIn BPG is better in terms of convergence with respect to iterations and competitive with respect to time. CoCaIn BPG CFI performs very similar to CoCaIn BPG and as anticipated the time plots illustrate that CoCaIn BPG CFI is slightly faster than CoCaIn BPG.

7 Acknowledgments

Mahesh Chandra Mukkamala and Peter Ochs acknowledge funding by the German Research Foundation (DFG Grant OC 150/1-1). Thomas Pock acknowledges support by the ERC starting grant HOMOVIS, No. 640156.

8 Appendix: Proof of Theorem 5.1

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

[TABLE]

We first prove the following result.

Lemma 8.1**.**

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

[TABLE]

In addition, the objective function $\Psi$ 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{}^{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 conditions (C1) and (C3), and the lower semicontinuity of $\Psi$ (which follows from the lower semi-continuity of $f$ and $g$, see Assumption A), we obtain

[TABLE]

and therefore

[TABLE]

On the other hand, from conditions (C1) and (C2), we know that there is $w^{k+1}\in\partial\Psi_{\delta_{1}}\left(x^{k+1},x^{k}\right)$, $k\in\mathbb{N}$, such that $w^{k+1}\rightarrow{\bf 0}$ as $k\rightarrow\infty$. The closedness property of $\partial\Psi_{\delta_{1}}$ implies thus that ${\bf 0}\in\partial\Psi_{\delta_{1}}\left(x^{\ast},x^{\ast}\right)=\left(\partial\Psi\left(x^{\ast}\right),{\bf 0}\right)$. This proves that $x^{\ast}$ is a critical point of $\Psi$, and hence (8.1) is valid.

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

We recall now the definition of the Kurdyka-Łojasiewicz (KL) property [26, 27] and [8] (for the non-smooth case). Denote $[\alpha<F<\beta]:=\left\{x\in{}^{d}:\;\alpha<F\left(x\right)<\beta\right\}$. Let $\eta>0$, and set

[TABLE]

Definition 8.1** (The Non-smooth KL Property).**

A proper and lower semicontinuous function $F:{}^{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]$.

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

Lemma 8.2** (Uniformized KL Property).**

Let $\Omega$ be a compact set and let $F:{}^{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 $\tilde{\varepsilon}>0$, $\eta>0$ and $\varphi\in\Phi_{\eta}$ such that for all $\overline{x}$ in $\Omega$ one has,

[TABLE]

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

We can now restate and prove Theorem 5.1.

Theorem 8.1**.**

Let $\{{x}^{{k}}\}_{{k}\in\mathbb{N}}$ be a bounded gradient-like descent sequence for minimizing $\Psi_{\delta_{1}}$. If $\Psi$ and $h$ satisfy 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}\,{\Psi}$.

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 8.1 we get that

[TABLE]

If there exists an integer $\bar{k}$ for which $\Psi_{\delta_{1}}\left(x^{\bar{k}+1},x^{\bar{k}}\right)=\Psi\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\{\Psi_{\delta_{1}}\left(x^{k+1},x^{k}\right)\right\}_{k\in\mathbb{N}}$ is a nonincreasing sequence, it is clear from (8.4) that $\Psi\left(\overline{x}\right)<\Psi_{\delta_{1}}\left(x^{k+1},x^{k}\right)$ for all $k>0$. Again from (8.4) for any $\eta>0$ there exists a nonnegative integer $k_{0}$ such that $\Psi_{\delta_{1}}\left(x^{k+1},x^{k}\right)<\Psi\left(\overline{x}\right)+\eta$ for all $k>k_{0}$. From Lemma 8.1 we know that $\lim_{k\rightarrow\infty}\mathrm{dist}\left(x^{k},\omega\left(x^{0}\right)\right)=0$. This means that for any $\tilde{\varepsilon}>0$ there exists a positive integer $k_{1}$ such that $\mathrm{dist}\left(x^{k},\omega\left(x^{0}\right)\right)<\tilde{\varepsilon}$ for all $k>k_{1}$.

From Lemma 8.1 applied to $\Psi_{\delta_{1}}$, we know that $\omega\left(x^{0}\right)$ is nonempty and compact and that the function $\Psi$ is finite and constant on $\omega\left(x^{0}\right)$. Hence, we can apply the Uniformization Lemma 8.2 applied to $\Psi_{\delta_{1}}$, which satisfies the KL property since $\Psi$ and $h$ do, with $\Omega=\omega\left(x^{0}\right)$. Therefore, for any $k\geq l:=\max\left\{k_{0},k_{1}\right\}+1$, we have

[TABLE]

This makes sense since we know that $\Psi_{\delta_{1}}\left(x^{k},x^{k-1}\right)>\Psi\left(\overline{x}\right)$ for any $k>l$. Combining (8.5) with condition (C2), see Proposition 5.4, 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 (8.6) and (8.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]

and thus

[TABLE]

Summing up (8.8) for $i=l+2,\ldots,k$ yields

[TABLE]

where the last equality 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+2,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 8.1, we obtain the global convergence to a critical point. ∎