Splitting methods with variable metric for KL functions

Nonconvex and nonsmooth optimization ; Kurdyka-Łojasiewicz inequality ; Descent methods ; Convergence rates ; Variable metric ; Gauss-Seidel method ; Newton-like method.

The second and third authors are partly supported by Conicyt Anillo Project ACT-1106, ECOS-Conicyt Project C13E03 and Millenium Nucleus ICM/FIC P10-024F. The third author is also partly supported by FONDECYT Grant 1140829 and Basal Project CMM Universidad de Chile.

P. Frankel & G. Garrigos Institut de Mathématiques et Modélisation de Montpellier, UMR 5149 CNRS. Université Montpellier 2, Place Eugène Bataillon, 34095 Montpellier cedex 5, France. Email: p.frankel30@orange.fr, guillaume.garrigos@gmail.com

G. Garrigos & J. Peypouquet Departamento de Matemática & AM2V. Universidad Técnica Federico Santa María, Avenida España 1680, Valparaíso, Chile. Email: guillaume.garrigos@gmail.com, juan.peypouquet@usm.cl

Introduction

In this paper we present a class of numerical methods to find critical points for a class of nonsmooth and nonconvex functions defined on a Hilbert space. Our analysis relies on the Kurdyka-Łojasiewicz (KŁ) inequality, initially formulated by Łojasiewicz for analytic functions in finite dimension , and later extended to nonsmooth functions in more general spaces . Gradient-like systems governed by potentials satisfying this KŁ inequality enjoy good asymptotic properties: under a compactness assumption, the corresponding trajectories have finite length and converge strongly to equilibria or critical points. These ideas were used in to study nonlinear first-order evolution equations (see also ). Second-order systems were considered in and a Schrödinger equation in .

The convergence analysis of algorithms in this context is more recent. See for gradient-related methods, for the proximal point algorithm and for a nonsmooth subgradient-oriented descent method. The celebrated Forward-Backward algorithm, a splitting method exploiting the nonsmooth/smooth structure of the objective function, has been studied in , and extended in to take in account a variable metric. Another splitting approach comes from Gauss-Seidel-like methods, which apply to functions with separated variables, and consist in doing a descent method relatively to each (block of) variables alternatively. See for a proximal alternating method, and for a variable-metric version. Recent papers propose to combine these two splitting approaches in order to exploit both the smooth/nonsmooth character and the separated structure of the function.

Most of the algorithms studied in the aforementioned papers share the same asymptotic behavior: under a compactness assumption, the sequences generated converge strongly to critical points, and the affine interpolations have finite length. This is not surprising since the algorithms described in together with the ones of (without extrapolation step) fall into the general convergence result for abstract descent methods of Attouch, Bolte and Svaiter . Besides, these methods essentially share the same hypotheses on the parameters with the abstract method of : the step sizes (resp. the eigenvalues of the matrices underlying the metric) are required to remain in a compact subinterval of the positive numbers. Moreover they have little flexibility regarding the presence of computational errors. To our knowledge vanishing step sizes (resp. unbounded eigenvalues) or sufficiently general errors have never been treated in the KŁ context.

Another interesting aspect is that the convergence rate of several of these methods are essentially the same, and depend on the KŁ inequality rather than the nature of the algorithm. Therefore, it seems reasonable to consider the existence of an abstract convergence rate result for general descent methods.

We present now the structure of the paper and underline its main contributions: in Section 2 we recall some definitions, well-known facts, and set the notation. Section 3 contains the main theoretical results of the paper. More precisely, in Subsection 3.1, we present an abstract inexact descent method, which is inspired by but extending their setting in order to account for additive computational errors and more versatility in the choice of the parameters. The strong convergence of the iterates with a finite-length condition, and a capture property are proved under certain hypotheses. Since the proofs are very close to those of , most arguments are given in Appendix A.1. Then, in Subsection 3.2 we prove new and interesting general convergence rates. They are similar to the ones obtained in . Surprisingly, an explicit form of the algorithm terminates in a finite number of iterations in several cases. A link with convergence rates for some continuous-time dynamical systems is also given. Sections 4 and 5 contain the main practical contributions. In Section 4, we present a particular instance of the model, which provides further insight into a large class of known methods and present some innovative variants. More exactly, we revisit the Alternating Forward-Backward methods, already considered in , but allowing inexact computation of the iterates and a dynamic choice of metric. This setting includes also the generalized Levenberg-Marquardt algorithm, a Newton-like method adapted for nonconvex and nonsmooth functions. In Section 5, we briefly describe an instance of this algorithm to produce a new method for the sparse and low-rank matrix decomposition. Finally, some perspectives are discussed in Section 6.

Preliminaries

If $x^{k}\overset{f}{\longrightarrow}x$ and $\liminf\limits_{n\to+\infty}\|\partial f(x^{k})\|_{-}=0$, then $0\in\partial f(x)$.

2. The Kurdyka-Łojasiewicz property

holds for all $x$ in the strict local upper level set

Semi-algebraic and bounded sub-analytic functions in finite dimension satisfy a KŁ inequality (), as well as some, but not all, convex functions (see for details and a counterexample). See , and the references therein, for more information in the general context of o-minimal functions. See for characterizations in infinite-dimensional Hilbert spaces.

3. Proximal operator in a given metric

Observe that $\mbox{prox}_{f}^{A}(x)\neq\emptyset$ if $f$ is weakly lower-semicontinuous and bounded from below (see [29, Theorem 3.2.5]), which holds in many relevant applications. If $f$ is the indicator function of a set, then $\mbox{prox}_{f}^{A}(x)$ is the nearest point mapping relatively to the metric induced by $A$.

Convergence of an abstract inexact descent method

$a_{k}\geq\underline{a}>0$ for all $k\geq 0$.

In Section 4, we complement this axiomatic description of descent methods by providing a large class of implementable algorithms that produce sequences verifying hypotheses $\mathbf{H}_{1}$, $\mathbf{H}_{2}$ and $\mathbf{H}_{3}$. A simple example is:

If $f$ is differentiable, a gradient-related method (see ) is an algorithms where each iteration has the form $x^{k+1}=x^{k}+\lambda_{k}d^{k}$, where $\lambda_{k}>0$ and $d^{k}$ agrees with the steepest descent direction $-\nabla f(x^{k})$ in the sense that $\langle d^{k},\nabla f(x^{k})\rangle+C\|d^{k}\|^{2}\leq 0$ and $\|\nabla f(x^{k})+d^{k}\|\leq C\|d^{k}\|+e_{k}$, with $C>0$ and $\lim_{k\to\infty}e_{k}=0$. If $\nabla f$ is Lipschitz-continuous, it is easy to find conditions on the sequence $(\lambda_{k})$ to verify hypotheses $\mathbf{H}_{1}$, $\mathbf{H}_{2}$ and $\mathbf{H}_{3}$.

Sequences generated by the procedure described above converge strongly to critical points of $f$ and the piecewise linear curve obtained by interpolation has finite length. More precisely, we have:

It is possible in Theorem 1 to drop the $f$-precompactness assumption and obtain a capture result, near a global minimum of $f$. To simplify the notation, for $x^{*}\in H$, $\eta\in]0,+\infty]$ and $\delta>0$, define the relaxed local upper level set by

As mentioned in , Theorem 2 admits a more general formulation, for instance, if $x^{*}$ is a local minimum of $f$ where a growth assumption is locally satisfied (see [17, Remark 2.11]).

The proofs of Theorems 1 and 2 follow the arguments in [17, Subsection 2.3], adapted to the presence of errors and the variability of the parameters. They are given in Appendix A.1 for the reader’s convenience.

2. Rates of Convergence

We assume that $\mathbf{H}_{1}$, $\mathbf{H}_{2}$ and $\mathbf{H}_{3}$ hold, and for simplicity and precision, we restrict ourselves to the case where $\varepsilon_{k}\equiv 0$. Suppose that $x^{k}$ $f$-converges to a point $x^{*}$ where $f$ has the KŁ property. We study three types of convergence rate results, depending on the nature of the desingularizing function $\varphi$:

Theorem 3 establishes the relationship between the distance to the limit $\|x^{k}-x^{*}\|$ and the gap $f(x^{k})-f(x^{*})$, for a generic desingularizing function. It is similar to the result in for the proximal method in the convex case.

Theorem 4 gives explicit convergence rates in terms of the parameters $-$ both for the distance and the gap $-$ when the desingularizing function is of the form $\varphi(t)=\frac{C}{\theta}t^{\theta}$ with $C>0$ and $\theta\in]0,1]$. Several results obtained in the literature for various methods are recovered.

Finally, Theorem 5 provides convergence rates when $\mathbf{H}_{2}$ is replaced by a slightly different hypothesis that holds for certain explicit schemes, namely gradient-related methods. This result is valid for a generic desingularizing function $\varphi$. However, when $\varphi$ is of the form $\varphi(t)=\frac{C}{\theta}t^{\theta}$ ($C>0$, $\theta\in]0,1]$) the prediction is considerably better than the one provided by Theorem 4.

for all $k\geq K$. Summing this inequality for $k=K,\dots,N$, we obtain

Using the triangle inequality and passing to the limit, we get

Theorem 4 below is qualitatively analogous to the results in : convergence in a finite number of steps if $\theta=1$, exponential convergence if $\theta\in[\frac{1}{2},1[$ and polynomial convergence if $\theta\in]0,\frac{1}{2}[$. In the general convex case, finite-time termination of the proximal point algorithm was already proved in and (see also ).

Assume $\varphi(t)=\frac{C}{\theta}t^{\theta}$ for some $C>0$, $\theta\in]0,1]$.

$f(x^{k})-f(x^{*})=O\left(\exp\left(-c{\sum\limits_{n=k_{0}}^{k-1}b_{n+1}}\right)\right)$, and

$\|x^{*}-x^{k}\|=O\left(\exp\left(-\dfrac{c}{2}{\sum\limits_{n=k_{0}}^{k-2}b_{n+1}}\right)\right)$.

$f(x^{k})-f(x^{*})=O\left(\left({\sum\limits_{n=k_{0}}^{k-1}b_{n+1}}\right)^{\frac{-1}{1-2\theta}}\right)$, and

$\|x^{*}-x^{k}\|=O\left(\left({\sum\limits_{n=k_{0}}^{k-2}b_{n+1}}\right)^{\frac{-\theta}{1-2\theta}}\right)$.

for each $k\geq k_{0}$. Let us now consider different cases for $\theta$:

Subcase $\theta\in[\frac{1}{2},1[$: Since $r_{k}\to 0$ and $0<2-2\theta\leq 1$, we may assume, by enlarging $k_{0}$ if necessary, that $r_{k+1}^{2-2\theta}\geq r_{k+1}$ for all $k\geq k_{0}$. Inequality (7) implies $(r_{k}-r_{k+1})\geq\beta_{k+1}r_{k+1}$ or, equivalently, $r_{k+1}\leq r_{k}\left(\dfrac{1}{1+\beta_{k+1}}\right)$ for all $k\geq k_{0}$. By induction, we obtain

for all $k\geq k_{0}$. But $\ln\left(\dfrac{1}{1+\beta_{n+1}}\right)\leq\dfrac{-\beta_{n+1}}{1+\beta_{n+1}}\leq\dfrac{-1}{1+\bar{b}}\beta_{n+1}$, and so

Subcase $\theta\in]0,\frac{1}{2}[$: Recall from inequality (7) that $r_{k+1}^{2\theta-2}(r_{k}-r_{k+1})\geq\beta_{k+1}$. Set $\phi(t):=\frac{C}{1-2\theta}t^{2\theta-1}$. Then $\phi^{\prime}(t)=-Ct^{2\theta-2}$, and

On the one hand, if we suppose that $r_{k+1}^{2\theta-2}\leq 2r_{k}^{2\theta-2}$, then

On the other hand, suppose that $r_{k+1}^{2\theta-2}>2r_{k}^{2\theta-2}$. Since $2\theta-2<2\theta-1<0$, we have $\frac{2\theta-1}{2\theta-2}>0$. Thus $r_{k+1}^{2\theta-1}>qr_{k}^{2\theta-1}$, where $q:=2^{\frac{2\theta-1}{2\theta-2}}>1$. Therefore,

with $C^{\prime}:=\frac{C}{1-2\theta}(q-1)r_{k_{0}}^{2\theta-1}>0$. Since $\beta_{k+1}\leq\frac{\bar{b}m}{C^{2}}$, we can write

Setting $c:=\min\{\frac{C}{2},\frac{C^{\prime}C^{2}}{\bar{b}m}\}>0$ we can write $\phi(r_{k+1})-\phi(r_{k})\geq c\beta_{k+1}$ for all $k\geq k_{0}$. This implies

which is precisely $r_{k+1}\leq D\left(\sum\limits_{n=k_{0}}^{k}b_{n+1}\right)^{\frac{-1}{1-2\theta}}$ with $D=\left(\frac{cm(1-2\theta)}{C^{3}}\right)^{\frac{-1}{1-2\theta}}$. As before, Theorem 3 gives the second part. ∎

2.3. Sharper results for gradient-related methods

Convergence rates for the continuous-time gradient system

where $f$ is some integral functional, are given in . For any $\varphi$, [34, Theorem 2.7] states that

$f(x^{k})-f(x^{*})=O\left(\Phi^{-1}(t-\hat{t})\right)$, and

$\|x^{*}-x^{k}\|_{L^{2}(\Omega)}=O\left(\varphi\circ\Phi^{-1}(t-\hat{t})\right)$,

$f(x^{k})-f(x^{*})=O\left(\Phi^{-1}\left(m\sum\limits_{n=k_{0}}^{k-1}b_{n+1}\right)\right)$, and

$\|x^{*}-x^{k}\|=O\left(\varphi\circ\Phi^{-1}\left(m\sum\limits_{n=k_{0}}^{k-1}b_{n+1}\right)\right)$.

To see this, let $k_{0}$ be large enough to have $x^{k}\in\Gamma_{\eta}(x^{*},\delta)$ where the KŁ inequality holds for all $k\geq k_{0}$. We apply successively $\mathbf{H}_{1}$, $\mathbf{H}_{2}^{\prime}$, the KŁ inequality and $\mathbf{H}_{3}$ to obtain

Let $\Phi$ be a primitive of $-(\varphi^{\prime})^{2}$. Then

because $\varphi^{\prime}$ is decreasing. Therefore,

as claimed. Now let us analyze the two cases:

For the second one, since $\varphi$ is concave and differentiable, we have

by $\mathbf{H}_{1}$. The KŁ property and $\mathbf{H}_{2}^{\prime}$ then give

Descent methods with errors and variable metric

As stressed in , the abstract scheme developed in Section 3 covers, among others, the gradient-related methods (a wide variety of schemes based on the gradient method sketched in ), the proximal algorithm (introduced in and further developed in ), and the forward-backward algorithm (a combination of the preceding, see ). This last one is a splitting method, used to solve structured optimization problems with the following form

It satisfies $\mathbf{H}_{1}$, $\mathbf{H}_{2}$ and $\mathbf{H}_{3}$ (see [17, Theorem 5.1]) and falls into the setting of Theorem 1. We shall extend this class of algorithms in different directions:

Alternative choice of metric for the ambient space, which may vary at each step (see and the references therein). Considering metrics induced by a sequence $(A_{k})\subset\mathcal{S}_{\text{\tiny++}}(H)$, the forward-backward method becomes

(recall Subsection 2.3). Indeed, (13) can be rewritten as

At each step, an approximation of $f$, replacing its smooth part $h$ by a quadratic model, is minimized. See for a similar algorithm called Variable Metric Forward-Backward, and for an approach considering more general models. Note that when $A_{k}=\frac{1}{\lambda_{k}}id_{H}$ one recovers (12). Allowing variable metric can improve convergence rates, help to implicitly deal with certain constraints, or compensate the effect of ill-conditioning. Rather than simply giving a convergence result for a general choice of $A_{k}$, we handle, in Subsection 4.3, a detailed method to select these operators, using second-order information.

where $g_{1},g_{2}$ are nonsmooth proper l.s.c functions and $h$ is differentiable with Lipschitz gradient. One approach is the regularized Gauss-Seidel method, which exploits the fact that the variables are separated in the nonsmooth part of $f$, as considered in . It consists in minimizing alternatively a regularized version of $f$ with respect to each variable. In other words, it is an alternating proximal algorithm, of the form:

But this algorithm does not exploit the smooth nature of $h$. An alternative is to use an alternating minimization method which can deal with the nonsmooth character, while it benefits from the smooth features. An Alternating Forward-Backward Method considering variable metrics is presented below. A constant-metric version, namely the Proximal Alternating Linearized Minimization Algorithm, can be found in . A forthcoming paper deals with the same algorithm, called Block Coordinate Variable Metric Forward-Backward, with a non-cyclic way of selecting the variables to minimize. Nevertheless, our setting differs from the aforementioned works in the following ways:

We allow more flexibility in the choice of parameters, accounting, in particular, for vanishing step sizes or unbounded eigenvalues for the metrics.

Convergence of this method with errors is given in Theorem 6.

Let $H_{1},\dots,H_{p}$ be Hilbert spaces, each $H_{i}$ provided with its own inner product $\langle\cdot,\cdot\rangle_{H_{i}}$ and norm $\|\cdot\|_{H_{i}}$. If there is no ambiguity, we will just note $\|x_{i}\|$ instead of $\|x_{i}\|_{H_{i}}$. Set $H=\prod\limits_{i=1}^{p}H_{i}$ and endow it with the inner product $\langle\cdot,\cdot\rangle={\sum\limits_{i=1}^{p}\langle\cdot,\cdot\rangle_{H_{i}}}$ and the associated norm $\|\cdot\|=\sqrt{\langle\cdot,\cdot\rangle}$. Consider the problem

has a $L$-Lipschitz continuous gradient. We shall present an algorithm that generates sequences converging to critical points of $f$. The sequences will be updated cyclically, meaning that given $(x_{1}^{k},...,x_{p}^{k})$, we start by updating the first variable $x_{1}^{k}$ into $x_{1}^{k+1}$, and then we consider $(x_{i}^{k+1},x_{2}^{k},...,x_{p}^{k})$ to update the second variable, and so on. In order to have concise and clear notations, throughout this section we shall denote:

Observe that $X_{1}^{k}=X^{k}$ and that we can write $X_{p+1}^{k}=X^{k+1}$.

We shall consider some hypotheses on the operators $A_{i,k}$. Define $\alpha_{k}=\min\limits_{i=1..p}\alpha(A_{i,k})$ and $\beta_{k}:=\max\limits_{i=1..p}\VERT A_{i,k}\VERT$, which give bounds on the spectral values of $(A_{i,k})_{i=1..p}$. We make the following assumptions:

Here $\mathbf{HP}_{1}$ is a bound on the spectral values by the Lipschitz constant of the gradient of $h$, in order to enforce the descent property of the sequence. For operators of the form $\frac{1}{\lambda_{i,k}}id_{H_{i}}$, we recover the classical bound $L\lambda_{i,k}\leq L\bar{\lambda}<1$. In , the authors prove that, with an additional convexity assumption on the $g_{i}$’s, and boundedness of the parameters, one can consider $L\lambda_{i,k}\leq L\bar{\lambda}<2$. Item $\mathbf{HP}_{2}$ states that the spectral values may diverge, but not too fast. Finally, $\mathbf{HP}_{3}$ can be seen as an hypothesis on the variations of the extreme spectral values of the chosen operators. It clearly holds for instance if $\beta_{k}$ is bounded. It is also sufficient to assume that the condition numbers

are bounded, with also $\min\left\{\frac{\alpha_{k}}{\alpha_{k+1}},\frac{\beta_{k}}{\beta_{k+1}}\right\}$ remaining bounded.

Even if $\nabla h$ is globally Lipschitz continuous, $L$ is not the Lipschitz constant of $\nabla h$ but a common Lipschitz constant for the functions defined in (18). As a consequence the partial gradients $\nabla_{i}h$ are $\sqrt{p}L$-Lipschitz continuous while $\nabla h$ is $pL$-Lipschitz. This allows us to have a better bound in $\mathbf{HP}_{1}$ which is of particular importance in the applications (see Section 5). In , the authors give a more precise analysis: at each substep $X_{i}^{k}$ of the algorithm, they consider $L_{i,k}$ as the Lipschitz constant of the gradient of $x\in H_{i}\mapsto h(x_{1}^{k+1},...,x_{i-1}^{k+1},x,x_{i+1}^{k},...,x_{p}^{k})$. Then they take step sizes equal to $\lambda_{i,k}=\frac{\epsilon_{i}}{L_{i,k}}$ where $\epsilon_{i}<1$ is a fixed non-negative constant. This approach can be related to the one in . However, they suppose a priori that the values $L_{i,k}$ remain bounded. It would be interesting to know if it is possible to combine the two approaches (a variable Lipschitz constant and vanishing step sizes).

2. The AFB method with errors

In order to allow for approximate computation of the descent direction or the proximal mapping, we go further by considering an inexact AFB method. We introduce the sequences $(r_{i}^{k})$ and $(s_{i}^{k})$ for $i\in\{1,...,p\}$ which correspond respectively to errors arising at the explicit and implicit steps relatively to the variable $x_{i}$. The AFB method with Errors is computed from an initial $(x_{1}^{0},...,x_{p}^{0})\in H$ by

We do specific hypothesis on the errors in view to guarantee the convergence of the method. Observe in particular that we do not assume a priori that the errors converge to zero:

This AFB algorithm (with errors) is related to the abstract descent method studied in Section 3. This is stated in the next proposition, whose proof is left in Appendix A.2.

Any sequence $Y^{k}=(y_{1}^{k},...,y_{p}^{k})$ generated by the AFB algorithm with errors satisfies $\mathbf{H}_{1}$, $\mathbf{H}_{2}$ and $\mathbf{H}_{3}$.

Given this result, one could directly apply Theorem 1 to obtain convergence of the sequence $(Y^{k})$ to a critical point of $f$. But this result would suffer from some drawbacks. First, we are expecting that $(X^{k})$ converges to a critical point, not $(Y^{k})$. So we should make the assumption that the errors $S^{k}:=X^{k}-Y^{k}$ tend to zero. Moreover we would suppose that $(Y^{k})$ is $f$-precompact, while we may only have an access to $(X^{k})$. To handle this, we make the link between the asymptotic behaviour of $(Y^{k})$ and $(X^{k})$:

For any sequence generated by the AFB method with errors:

If $(Y^{k})$ has finite length, then so does $(X^{k})$.

$(Y^{k})$ is precompact if and only if $(f(Y^{k}))$ is bounded from below and $(X^{k})$ is precompact.

Item 1 comes directly from $\text{{HE}}_{1}$. To prove item 2, we use Proposition 1: from $\text{{H}}_{1}$ and $\text{{H}}_{3}(i)$ we have that

hence $(f(Y^{k}))$ is a decreasing sequence. Then we can sum inequality (21) to obtain that

An other disadvantage to the direct application of Theorem 1 is that it asks the $f$-precompactness of $(Y^{k})$. In some cases, precompactness of a sequence can be deduced using compact embeddings between Hilbert spaces. Sequences remaining in a sublevel set of an inf-compact function $f$ are also precompact. However, $f$-precompactness is harder to obtain without further continuity assumption on $f$. Actually, both limit and $f$-limit points coincide whenever the parameters are bounded:

If either $\beta_{k}\leq\bar{\beta}$ or $f$ is continuous on its domain, then $(Y^{k})$ is $f$-precompact if and only if it is precompact.

and the latter implies (using Cauchy-Schwartz and $\VERT A_{i,k}\VERT\leq\bar{\beta}$):

Now recall that $y_{i}^{k+1}=y_{i}^{k_{n}}$ tends to $y_{i}^{\infty}$ while $r_{i}^{k}+s_{i}^{k}$ goes to zero (see Proposition 2). Observe also that $\nabla_{i}h(Y_{i}^{k}+S_{i}^{k})$ is bounded since it converges to $\nabla_{i}h(Y^{\infty})$. Moreover, $\|y_{i}^{\infty}-y_{i}^{k}\|$ goes also to zero since we have

As a direct consequence of Propositions 1, 2, 3 together with Theorem 1, we finally get our convergence result for the AFB algorithm with errors. It extends the results of (when taking a cyclic permutation on the variables) in two directions: the functions $g_{i}$ need not be continuous on their domain, or the step sizes can tend to 0.

Let $f$ be a KŁ function. Let $(Y^{k})$ be a precompact sequence generated by the AFB algorithm with errors, with (HP) and (HE) satisfied. Suppose that either $\beta_{k}$ remains bounded, or that $f$ is continuous on its domain. Hence, the sequence $(X^{k})$ has finite length and converges toward a critical point of $f$.

An analog of the capture result in Theorem 2 can also be deduced:

Suppose that the KŁ property holds in a global minimum $X^{*}$ of $f$. Let $(X^{k})$ be a sequence generated by the AFB algorithm with errors, satisfying (HP) and (HE) with $\mu_{k}\equiv 0$. Hence, there exist $\gamma>0$ and $\eta>0$ such that if $X^{0}\in\underline{\Gamma}_{\eta}(X^{*},\gamma)$, then $(X^{k})$ has finite length and converges to a global minimum of $f$.

To prove this theorem, it suffices to use $Y^{0}=X^{0}$, and to see at the end of the proof of Proposition 1 that $\mu_{k}=0$ iff $\epsilon_{k}=0$, where $\epsilon_{k}$ is the parameter involved in $\mathbf{H}_{3}$. Then, apply Theorem 2 together with Propositions 1 and 2.

3. Variable metric: towards generalized Newton methods

gives the minimum of $h$ over $C$ in one single step. For a general function $h$, (24) reduces to the minimisation over $C$ of a quadratic model of $h$, as stressed in (14). One can see on this example that computing the proximal operator relatively to the metric $A_{n}$ used in the explicit step (and not the ambient metric !) is of crucial importance in this method.

The spirit here is to use second-order information from $h$ in order to improve the convergence of the method. In the unconstrained case, a popular choice of metric is given by Newton-like methods, where the metric at step $k$ is induced by (an approximation of) the Hessian $\nabla^{2}h(x^{k})$. Since it is often impossible to know in advance whether or not the Hessian is uniformly elliptic at each $x^{k}$, a positive definite approximation has to be chosen.

We detail here a natural way to chose this positive definite $A_{k}\sim\nabla^{2}h(x^{k})$ in closed loop, and show that this method remains in the setting of Theorem 6. Since it generalizes the Levenberg-Marquardt method used in the convex case (see ) we will refer to the Generalized Levenberg-Marquardt method for this way of designing $A_{k}$. One of the interesting aspect of the method is that such a matrix can be defined even if $h$ is only $C^{1,1}$ and not $C^{2}$, since the differentiability of $\nabla h$ is not necessary in Theorem 6. Another interesting aspect is that the splitting approach led us to solve constrained minimization problems with a Newton-projected approach.

A globalized version of the method can be considered by taking step sizes ensuring descent. Then the following convergence result holds:

where $A_{k}$ is selected with the Generalized Levenberg-Marquardt process detailed above, and the stepsizes $\lambda_{k}$ satisfy:

Then the sequence has finite length and is converging to a critical point of $f$.

Start by observing that $\mbox{\rm proj\,}_{C}^{A_{k}}=\mbox{\rm proj\,}_{C}^{\lambda_{k}^{-1}A_{k}}$, so the algorithm falls in the setting of the AFB algorithm. According with the previous notations, $\nabla h$ being $L$-Lipschitz continuous implies that the sequence $(\VERT H_{k}\VERT)$ is bounded by $L$, and so $(\VERT P_{k}\VERT)$ remains bounded by $2L$. To conclude through Theorem 6 we just need to check the hypotheses (HP) on the parameters $\frac{1}{\lambda_{k}}A_{k}$. We have here $\alpha_{k}=\alpha(\frac{1}{\lambda_{k}}A_{k})\geq\varepsilon{\lambda_{k}}^{-1}\geq\varepsilon\bar{\lambda}^{-1}>L$ and $\beta_{k}=\VERT\frac{1}{\lambda_{k}}A^{k}\VERT\leq(2L+\epsilon)\lambda_{k}^{-1}.$ Thus $\mathbf{HP}_{1}$ is satisfied, while items $\mathbf{HP}_{2}$ and $\mathbf{HP}_{3}$ follows directly from the hypotheses made on $(\lambda_{k})$. Since the indicator function $\delta_{C}$ is continuous on its domain, the hypotheses of Theorem 6 are satisfied.∎

This extends, in a way, results from the convex setting to the nonconvex one, enforcing moreover the strong convergence (see [42, Theorem 7.1]).

A drawback of this method is that the Hessian increases the complexity of implementation since a matrix must be inverted in the explicit step. An alternative is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) update scheme (see ,), using only first-order information to compute the inverse of the Hessian. On the other hand, the implicit step gains also in complexity since one must project onto a constraint relatively to a given metric, which is nontrivial even for simple constraints. For linear constraints, a particular second-order model of the Hessian can be taken in order to reduce the implicit step in a trivial orthogonal projection step (see ).

Newton-like methods are expected to have good convergence rates in exchange for a more expensive implementation. An interesting question is whether one can obtain convergence rates beyond the results in Subsection 3.2, by exploiting, not only the KŁ nature of the function, but also the specific properties of the matrices selected by the Generalized Levenberg-Marquardt process.

Applications

The framework presented in this paper is suitable for the numerical resolution of a wide variety of structured problems. Consider for instance the problems arising in image processing and data compression, which are generally semi-algebraic by nature . Indeed, they generally involve the semi-algebraic counting norm $\|x\|_{0}:=\sharp\{i\ |\ x_{i}\neq 0\}$, whose proximal operator (the hard shrinkage operator, see ) is easily implementable. Feasability problems with semi-algebraic (eventually nonconvex) constraints are also well suited for the AFB method (see ). The search for equilibria of nonlinear partial differential equations has already been tackled using the KŁ inequality . It should now be improved by using splitting methods more adapted to the structure of the problem. Let us end by discussing in some detail the sparse and low-rank matrix decomposition, for which the AFB method is particularly well adapted, in view of its structure.

The KŁ framework is well adapted to the original nonconvex (but semialgebraic!) problem and offers convergent numerical methods. Moreover, the AFB method is well suited for its structure in separated variables involving smooth and nonsmooth parts. It leads to an Alternating Averaged Projected Method: given $(X_{0},Y_{0})$, take $(\lambda_{k})$, $(\mu_{k})$ with $0<\underline{\tau}\leq\lambda_{k},\ \mu_{k}\leq\bar{\tau}<1$. For $k\geq 0$, define

Projection onto $\{\mbox{\rm rank\,}\cdot\leq r\}$ can be done using the Singular Value Decomposition (see Eckart-Young’s Theorem). To project onto $\{\|\cdot\|_{0}\leq s\}$, one simply sets all the coefficients to zero, except for the $s$ largest ones (in absolute value). Theorem 7 guarantees convergence to the solution for sufficiently close initialization. This example illustrate the discussion in Remark 2: here we have $L=1$, while if one consider the Lipschitz constant of the gradient of $(X,Y)\mapsto\frac{1}{2}\|A-X-Y\|_{F}^{2}$, we would have had $L=2$, that is a strictly smaller upper bound for the parameters.

Concluding Remarks

We have given a unified way to handle various recent descent algorithms, and derived general convergence rate results in the KŁ framework. These are applicable to potential future numerical methods. Some improvements have been explored, and a novel projected Newton-like method has been proposed.

A challenging task is to extend the present convergence analysis to algorithms that do not satisfy the sufficient decrease condition $\mathbf{H}_{1}$. This will allow to consider acceleration schemes like the ones studied in , or primal-dual methods based on a Lagrangian approach. A recent preprint seems to be an interesting first attempt in this direction.

Finally, it is worth mentioning that the results in Section 3 remain true in the more general context of a normed space, adapting the definition of subdifferential and lazy slope in an obvious manner.

Acknowledgements : The authors thank H. Attouch for useful remarks. They would also like to thank the anonymous reviewer for his careful reading and constructive comments.

Appendix A Appendix

The argument is a straightforward adaptation of the ideas in the proof of [17, Lemma 2.6]. One first proves:

$\mathbf{S}(x^{*},\delta,\rho)$: There exist $\delta>\rho>0$ such that

The initial point $x^{0}$ belongs to $\Gamma_{\eta}(x^{*},\rho)$ and

The basic asymptotic properties are given by the following result:

Let $\mathbf{H}_{1}$, $\mathbf{H}_{2}$, $\mathbf{H}_{3}$ and $\mathbf{S}(x^{*},\delta,\rho)$ hold. Then $x^{k}\in\underline{\Gamma}_{\eta}(x^{*},\rho)$ for all $k$ and converges to some $\overline{x}$ lying in the closed ball $\overline{B(x^{*},\rho)}$. Moreover $\sum_{k=1}^{\infty}\|x^{k+1}-x^{k}\|<\infty$, $\liminf_{k\to\infty}\|\partial f(x^{k})\|_{-}=0$, and $f(\overline{x})\leq\lim_{k\to\infty}f(x^{k})=f(x^{*})$.

We are now in position to complete the proofs of Theorems 1 and 2.

which implies part i) of $\mathbf{S}(x^{*},\delta,\rho)$. Now take $K\geq K_{0}$ such that

It follows that $\|x^{k+1}-x^{*}\|\leq\|x^{k+1}-x^{k}\|+\|x^{k}-x^{*}\|<\sqrt{\frac{\eta}{\underline{a}}}+\rho<\frac{4}{3}\rho<\delta$, and so $x^{k+1}\in\underline{\Gamma}_{\eta}(x^{*},\delta)$. Finally, we have

A.2. Proof of Proposition 1

Since $X_{i}^{k}=Y_{i}^{k}+S_{i}^{k}$, we can rewrite the algorithm as

We start by showing that $\mathbf{H}_{1}$ is satisfied.

Let $i=1..p$ be fixed. Using the definition of the proximal operator $\mbox{prox}_{g_{i}}^{A_{i,k}}$ in (27) and developing the squared norms gives

Using $\mathbf{HE}_{3}$ in (28), the latter results in

where $\rho A_{i,k}-Lid_{H_{i}}$ remains coercive, since $\rho\alpha_{k}>L$. Using successively the Cauchy-Schwartz inequality, the Lipschitz property of $\nabla_{i}h$ (see Remark 2) and $\mathbf{HE}_{1}$, one gets

Inserting this estimation in (32) we deduce that

We can now conclude by summing all these inequalities for $i=1,\dots,p$:

so $\mathbf{H}_{1}$ is fulfilled with $a_{k}=\frac{\rho\alpha_{k}-L(\sigma\frac{\sqrt{p}}{p}+1)}{2}$. To prove $\mathbf{H}_{2}$, fix $i=1,\dots,p$ and use Fermat’s first order condition in (27) to get:

Define $w_{i}^{k+1}:=\nabla_{i}h(Y^{k})-\nabla_{i}h(Y_{i}^{k}+S_{i}^{k})-A_{i,k}(y_{i}^{k+1}-y_{i}^{k})+A_{i,k}(r_{i}^{k}+s_{i}^{k})$ which lies in $\partial g_{i}(y_{i}^{k+1})+\nabla_{i}h(Y^{k+1})$, by (35). The triangle inequality gives

where we use the error estimations from (HE)

and the $\sqrt{p}L$-Lipschitz continuity of $\nabla_{i}h$:

Define now $W^{k+1}:=(w_{1}^{k+1},...,w_{p}^{k+1})\in\partial f(Y^{k+1})$ (recall the definition of $w_{i}^{k+1}$). Then through the sum over $i=1..p$ of inequality (39) we have (using $\sqrt{p}\leq p\leq p^{2}$)

Hence $\mathbf{H}_{2}$ is verified with $b_{k+1}=\frac{1}{p^{2}(1+\sigma)(\beta_{k}+L)}$ and $\epsilon_{k+1}=\frac{\beta_{k}\mu_{k}}{p(1+\sigma)(\beta_{k}+L)}$.

References