Variable metric inexact line-search based methods for nonsmooth optimization

Introduction

When in particular $f_{1}$ reduces to the indicator function of a convex set $\Omega$, i.e. $f_{1}=\iota_{\Omega}$ with

a simple and well studied algorithm for the solution of (1) is the gradient projection (GP) method, which is particularly appealing for large scale problems. In the last years, several variants of such method have been proposed , with the aim to accelerate the convergence which, for the basic implementation, can be very slow. In particular, reliable acceleration techniques have been proposed for the so called gradient projection method with line–search along the feasible direction [6, Chapter 2], whose iteration consists in

where ${y}^{(k)}$ is the Euclidean projection of the point ${x}^{(k)}-\nabla f_{0}({x}^{(k)})$ onto the feasible set $\Omega$ and ${\lambda^{(k)}}\in$ is a steplength parameter ensuring the sufficient decrease of the objective function. Typically, ${\lambda^{(k)}}$ is determined by means of a backtracking loop until an Armijo-type inequality is satisfied. Variants of the basic scheme are obtained by introducing a further variable stepsize parameter ${\alpha_{k}}$, which controls the step along the gradient, in combination with a variable choice of the underlying metric. In practice, the point ${y}^{(k)}$ can be defined as

In this paper we generalize the GP scheme (2)–(3), by introducing the concept of descent direction for the case where $f_{1}$ is a general convex function and we propose a suitable variant of the Armijo rule for the nonsmooth problem (1). In particular, we focus on the case when the descent direction has the form ${y}^{(k)}-{x}^{(k)}$, with

Formally, the scheme (2)-(4) is a forward–backward (or proximal gradient) method depending on the parameters ${\lambda^{(k)}}$, $\sigma^{(k)}$.

Forward–backward algorithms based on a variable metric have been recently studied also in for the convex case and in for the nonconvex case under the Kurdyka-Łojasiewicz assumption (see also ). Even if our scheme is formally very similar to those in , the involved parameters have a substantially different meaning. In our case, the theoretical convergence is ensured by the Armijo parameter ${\lambda^{(k)}}$ in combination with the descent direction properties; this results in an almost complete freedom to choose the other algorithm parameters (e.g. ${\alpha_{k}}$ and $D_{k}$), without necessarily relating them to the Lipschitz constant of $\nabla f_{0}$ (actually, our analysis, except the convergence rate estimate, is performed without this assumption). We believe that this is also one of the main strength of our method, since acceleration techniques based on suitable choices of ${\alpha_{k}}$ and $D_{k}$, originally proposed for smooth optimization, can be adopted, leading to an improvement of the practical performances. The other crucial ingredient of our method is the inexact computation of the minimizer in (4): this issue has been considered in several papers in the context of proximal and proximal gradient methods (see for example and references therein). The approach we follow in this paper is more similar to the one proposed in and has the advantage to provide an implementable condition for the approximate computation of the proximal point. Moreover, we also generalize the ideas proposed in for the inexact computation of the projection onto a convex set. Finally, we also mention the papers for the use of non Euclidean distances in the context of forward–backward and proximal methods.

The paper is organized as follows: some background material is collected in Section 2, while the concept of descent direction for problem (1) is presented and developed in Section 3. In Section 4, the modified Armijo rule is discussed. Then, a general convergence result for line–search descent algorithms based on this rule is proved, in the nonconvex case. Two different inexactness criteria, called of $\epsilon$-type and $\eta$-type are proposed in Sections 4.2 and 4.3, and the related implementation is discussed in Sections 5.1 and 5.4. Section 4.5 deals with the convex case, where the convergence of an $\epsilon$-approximation based algorithm is proved and the related convergence rate is analyzed. The results of a numerical experience on a total variation based image restoration problem are presented in Section 6 while our conclusions are given in Section 7.

Definitions and basic properties

When $f$ is smooth at ${x}$, then $f^{\prime}({x};{d})=\nabla f({x})^{T}{d}$. When $f$ is convex, its directional derivative has the following property.

The following proposition states a useful property of the conjugate.

A family of descent directions

Thanks to Theorem 2.1, the previous definition is well posed. In this section we define a family of descent directions for problem (1). To this end, we define the following set of non–negative functions.

$d_{{\sigma}}({z},{x})$ is continuous in $(\sigma,{z},{x})$;

$d_{{\sigma}}({z},{x})$ is smooth w.r.t. ${z}\in\Omega$;

$d_{{\sigma}}({z},{x})$ is strongly convex w.r.t. ${z}$:

where $m>0$ does not depend on $\sigma$ or ${x}$ (here $\nabla_{1}$ denotes the gradient with respect to the first argument of a function);

$d_{{\sigma}}({z},{x})=0$ if and only if ${z}={x}$ (which implies that $\nabla_{1}d_{{\sigma}}({x},{x})=0$ for all ${x}\in\Omega$).

where $h_{{\sigma}}^{\prime}({z},{x};{d})$ denotes the directional derivative of $h_{{\sigma}}(\ \cdot\ ,{x})$ at the point ${z}$ with respect to $d$. From assumption $({\mathcal{D}}_{3})$, it follows that $h_{{\sigma}}(\cdot,{x})$ is strongly convex and admits a unique minimum point for any ${x}\in\Omega$.

Now we introduce the following operator ${p}:\Omega_{0}\rightarrow\Omega$ associated to any function $h_{{\sigma}}$ of the form (11)

When $d_{{\sigma}}$ is chosen as in (10), the operator (13) becomes

Under assumption $({\mathcal{D}}_{3})$, one can show that ${p}({x};h_{{\sigma}})$ depends continuously on $({x},\sigma)$.

Let $d_{{\sigma}}\in{\mathcal{D}}(\Omega,S)$ and $h_{{\sigma}}$ be defined as in (11). Then ${p}({x};h_{{\sigma}})$ depends continuously on $({x},\sigma)$.

Assumption $({\mathcal{D}}_{3})$ expressed in ${y}$ and $u$ gives:

Let ${y}_{1}={p}({x}_{1};h_{\sigma_{1}})$ and ${y}_{2}={p}({x}_{2};h_{\sigma_{2}})$. Adding the previous inequality for ${y}={y}_{1}$ (resp. ${y}={y}_{2}$) and choosing $u={y}_{2}$ (resp. $u={y}_{1}$), one finds:

the stationarity condition (8) can be reformulated in terms of fixed points of the operator ${p}(\ \cdot\ ;h_{{\sigma}})$;

To this purpose, we collect in the following proposition some properties of the function $h_{{\sigma}}$ and the associated operator ${p}(\ \cdot\ ;h_{{\sigma}})$.

Proof. (a) is a direct consequence of definition (14) and condition $({\mathcal{D}}_{3})$ on $d_{{\sigma}}$.

Conversely, let ${x}\in\Omega$ be a stationary point of (1) and assume by contradiction that ${x}\neq{y}$. Then, by Proposition 3.2 (d) we obtain $f^{\prime}({x},{y}-{x})<0$, which contradicts the stationarity assumption on ${x}$.

(b) $\Longleftrightarrow$ (c) See Proposition 3.2 (c). $\square$

A line–search algorithm based on a modified Armijo rule

In this section we consider the modified Armijo rule described in Algorithm LS, which is a generalization of the one in . Indeed the rule proposed in is recovered when $d_{{\sigma}}$ is chosen as in (10) and $\gamma\in[0,1)$. In the following we will prove that Algorithm LS is well defined and classical properties of the Armijo condition still hold for this modified case.

Here and in the following we will define the function $h_{{\sigma}}(\cdot,\cdot)$ as in (11) and, for sake of simplicity, we will make the following assumption

where $\lambda^{(k)}$ and ${d}^{(k)}$ are computed with Algorithm LS, then we have

where the second inequality is obtained by means of the Jensen inequality applied to the convex function $f_{1}$. Taking limits on the right hand side for $j\rightarrow\infty$ we obtain

where the second inequality follows from the non–negativity of $d_{{\sigma}}\in{\mathcal{D}}(\Omega,S)$ and the last one from (18). Since $0<\beta<1$, this is an absurdum.

where the first inequality follows from the non–negativity of $d_{{\sigma}}$, the second one is obtained by adding and subtracting $f_{0}({x}^{(k)})$ and the last one is a consequence of $f({x}^{(k+1)})\leq f({x}^{(k)})$.

Let us show that the only limit point of $\Delta^{(k)}$ is zero. We observe that from (18) and (19) we obtain

for all sufficiently large $k\in\bar{K}$. Repeating the same arguments employed in the first part of the proof, we obtain

Summing the previous inequality for $k=0,...,j$ gives

Proposition 4.1 allows the convergence analysis of a wide class of descent methods based on the Armijo condition (16). The crucial ingredients of these methods are

Then $\bar{x}$ is a stationary point for problem (1).

Proof. First, we notice that Algorithm LS is well defined, since (18) holds. We observe that, since $h_{\sigma^{(k)}}$ is strongly convex with modulus of convexity $m$ and ${y}^{(k)}$ is its minimum point, we have

Thus we can apply Proposition 4.1 and obtain

Combining the previous equality with (15) and (25) yields

Since $h_{\sigma^{(k)}}({y}^{(k)},{x}^{(k)})\leq 0$, this implies $\lim_{k\rightarrow\infty,k\in K}h_{\sigma^{(k)}}({y}^{(k)},{x}^{(k)})=0$. Expressing inequality (26) for ${z}={x}^{(k)}$, we can write

2 ϵitalic-ϵ\epsilon- approximations

In this section we will assume that $d_{{\sigma}}$ has the form (10) and, in this case, we will describe a sufficient condition for (25).

We recall that $h_{{\sigma}}(\ \cdot\ ,{x})$ is strongly convex with modulus $m=2/(\alpha\mu)$ and ${y}$ is its minimizer. This yields

where the rightmost inequality follows from (31) with $w={y}$. $\square$ The previous result combined with Theorem 4.1 directly implies the following Corollary.

3 η𝜂\eta-approximations

A different approach to define a suitable approximation of the operator (13) is based on the following definition.

for some $\eta\in(0,1]$. This idea of inexactness was introduced first in to approximate the projection operator onto a convex set in the context of scaled gradient projection methods for smooth optimization. Clearly, if

Proof. We set ${y}^{(k)}={p}({x}^{(k)};h_{\sigma^{(k)}})$ and we first observe that $\gamma\leq 1$ and (35) imply

Since $h_{\sigma^{(k)}}(\cdot,{x}^{(k)})$ is strongly convex with modulus of convexity $m$, and ${y}^{(k)}$ is the minimizer of $h_{\sigma^{(k)}}(\cdot,{x}^{(k)})$, we can write

which, since $h_{\sigma^{(k)}}({y}^{(k)},{x}^{(k)})\leq 0$, implies

Invoking again the strong convexity of $h_{\sigma^{(k)}}(\ \cdot\ ,{x}^{(k)})$, we obtain

with, together with (38) gives $\lim_{k\to\infty,k\in K}\|{y}^{(k)}-{x}^{(k)}\|^{2}=0$. Thus, $\bar{y}=\bar{x}$ and by Proposition 3.3, we conclude that $\bar{{x}}$ is stationary. $\square$

4 Remarks

Different notions of inexactness have been proposed in the literature (see and references therein), especially in the context of proximal point methods, with the aim of approximating the resolvent operator, and some of them could be considered also in our framework. A synthetic description of possible inexactness notions and their relationships is given in Figure 1.

It is difficult to insert the inexactness criterion (34) in the scheme in Figure 1, since the shape of $P_{\eta}$ in (34) depends on ${x}$, while the implications in Figure 1 are independent of ${x}$. In general, we observe that from inequality (37) and by definition of $\epsilon$-subdifferential we have

5 Convergence analysis in the convex case with ϵitalic-ϵ\epsilon-approximations

Proof. First of all we recall the basic norm equality

which, recalling that ${z}^{(k)}={x}^{(k)}-{\alpha_{k}}D_{k}^{-1}\nabla f_{0}({x}^{(k)})$, writes also as

For $w=\hat{x}$, the previous inequality gives

where the second inequality is obtained adding and subtracting $f_{1}({x}^{(k)})$ and by the convexity of $f_{0}$, the third one from the fact that $\hat{x}$ is a minimum point and the last one by definition of ${x}^{(k+1)}$. By equality (40) with $a={x}^{(k+1)}$, $b={x}^{(k)}$, $c=\hat{x}$, $D=D_{k}$ we obtain

Proof. By substituting (45) in (44) we obtain

The rest of the proof follows exactly from the same arguments employed in Theorem 4.3. $\square$ We will show in Section 5.4 how the conditions (32) and (45) can be satisfied in practice.

Assumption (H2) is analogous to the one proposed in . A special case of it consists in the following

which implies ${D_{k+1}}\preceq\mu_{k}\mu_{k+1}D_{k}$. Moreover, $\mu_{k}\mu_{k+1}$ can be written as $\mu_{k}\mu_{k+1}=1+\zeta_{k}$, where $\zeta_{k}=\sqrt{(1+\xi_{k})(1+\xi_{k+1})}-1$. Since $\lim_{x\to 0}\sqrt{1+x}/x=1/2$, it follows that $\sum_{k=0}^{\infty}\xi_{k}$ and $\sum_{k=0}^{\infty}\zeta_{k}$ have the same behaviour. Then, we can conclude that (H2’) implies (H2).

We also observe that, employing the same arguments above, we can also prove that $\mu_{k+1}\mu_{k}{D_{k+1}}\succeq{D_{k}}$, and, as a consequence, (H2’) also implies that $(1+\zeta_{k}){D_{k+1}}\succeq{D_{k}}$ with $\sum_{k=0}^{\infty}\zeta_{k}<\infty$. In practice, (H2’) says that the scaling matrices have to converge to the identity matrix at a certain rate, while (H2) implies the convergence to some symmetric positive definite matrix (see Lemma 2.3 in ).

5.2 Convergence rate analysis

where the last inequality follows from (5). $\square$

Proof. In view of (45)–(46), setting $a=\alpha_{\max}\mu$, one obtains

If $\nabla f_{0}$ is Lipschitz continuous on $\Omega$ with Lipschitz constant $L$, then from the descent lemma [6, p.667] we have

where $\lambda\in$. By combining inequalities (48) and (49) we further obtain

The previous inequality ensures that the Armijo condition

Assume that the hypotheses of Theorem 4.4 hold and, in addition, that the gradient of $f_{0}$ is Lipschitz continuous on $\Omega$. Let $f^{*}$ be the optimal function value for problem (1). Then, we have

Proof. If we do not neglect the term $f({x}^{(k)})-f(\hat{x})=f({x}^{(k)})-f^{*}$ in (41) and in all the subsequent inequalities, instead of (43) we obtain

where we set $\zeta=1+\max_{k}\zeta_{k}$, $a=2\lambda_{\min}\alpha_{\min}$, where $\lambda_{\min}$ is defined in Proposition 4.2. Summing the previous inequality from 0 to $k$ gives

Practical computation of η𝜂\eta- and ϵitalic-ϵ\epsilon- approximations

Moreover, for all sufficiently large $l$ we have $a_{l}>h_{\sigma^{(k)}}({y}^{(k)},{x}^{(k)})/\eta$ which, together with (54) gives

2 Composition with a linear operator

In this section we assume that $f_{1}(x)$ is given by

with ${z}^{(k)}={x}^{(k)}-{\alpha_{k}}D_{k}^{-1}\nabla f_{0}({x}^{(k)})$. The dual problem is obtained by computing the minimum of the primal–dual function with respect to ${y}$, which is given by ${y}={z}^{(k)}-{\alpha_{k}}D_{k}^{-1}A^{T}v$, and substituting it in (57), obtaining the explicit expression of the dual function

By definition of the primal–dual and dual functions, the following inequalities hold

is satisfied, i.e. (55) with $a_{l}=\Psi_{\sigma^{(k)}}(v^{(l)},{x}^{(k)})$.

For example, one can apply a forward–backward method , called also ISTA or its accelerated version (FISTA, ) to the dual problem. As an alternative, also the saddle point problem

3 Preserving feasibility

4 Computing ϵitalic-ϵ\epsilon-approximations

Proof. From the definition of the primal–dual gap, a simple computation shows that

where the last inequality follows from Proposition 2.1. Thus, if (62) holds, the previous inequality yields

Rearranging terms, the previous inequality writes also as

5 Equivalence between η𝜂\eta and ϵitalic-ϵ\epsilon approximations

Numerical illustration

In order to validate the proposed approach, we consider a relevant image restoration problem, whose variational formulation consists in minimizing the sum of a discrepancy functional plus a regularization term. Following the Bayesian paradigm, when the noise affecting the data is of Poisson type, a typical choice for measuring the discrepancy of a given image ${x}$ from the observed data $b$ is the following Kullback-Leibler divergence

In our experiments we assume that $H$ corresponds to a convolution operator associated to a Gaussian kernel, with reflective boundary conditions, so that the matrix-vector products involving $H$ can be performed via the Discrete Cosine Transform .

We implement our inexact algorithm, which is summarized in Algorithm VMILA, in Matlab environment with the following settings:

Step 1, metric selection: the scaling matrix $D_{k}$ is chosen mimicking the split-gradient idea . In particular, at each outer iteration it is defined as the diagonal matrix with positive entries as follows

where $\mu_{k}=\sqrt{1+10^{10}/k^{2}}$, so that assumption (H2’) is satisfied. We choose a large initial range for the scaling matrix selection to allow more freedom of choice at the first iterates, where the benefits of the scaling matrix are more relevant .

Step 1, steplength selection: the parameter $\alpha_{k}$ is chosen by the same strategy used e.g. in , and its value is constrained in the interval $[\alpha_{\min},\alpha_{\max}]$ with $\alpha_{\min}=10^{-5}$, $\alpha_{\max}=10^{2}$.

Other parameters setting: the line–search parameters $\delta,\beta,\gamma$ have been set respectively equal to $0.5,10^{-4},1$.

All the following results have been obtained on a PC equipped by an Intel Core i7-2620M processor with CPU at 2.70GHz and 8GB of RAM, running Windows 7 OS and MATLAB Version 7 (R2010b).

We investigate first the impact of the inexactness parameter $\eta$ choice on the overall method. In Figure 3 the relative decrease of the objective function values in the first 500 iterates is reported with respect to both the iteration number (first row) and the computational time, in seconds (second row). It can be observed that a higher precision can accelerate the progress toward the solution, but this usually results in a very large number of inner iterations and, consequently, it is extremely time consuming (for example, for the test problem cameraman with $\eta=10^{-6},10^{-2},5\cdot 10^{-1}$ the mean number of inner iterations per outer iteration is 28, 54, 409, respectively). This is typical of inexact algorithms based on the iterative solution of an inner subproblem. We find that a good balance between convergence speed and computational cost is obtained by allowing a relatively large tolerance, corresponding to $\eta=10^{-6}$.

Conclusions and future work

In this paper we presented and analyzed an inexact variable metric forward–backward method based on an Armijo–type line–search along a suitable descent direction. The inexactness of the method relies in the possibility of using an approximation of the proximal operator, while the underlying metric may change at each iterations and also non Euclidean metrics are allowed. We performed the convergence analysis of the method, obtaining results in both the nonconvex and convex cases and providing also a convergence rate estimate in the latter one. The main strengths of the method are listed below.

The convergence is ensured by a line–search procedure, which does not depend on any user supplied parameter (actually the constants $\gamma,\beta,\delta$ have to be chosen, but the behaviour of the whole algorithm is not sensitive to these choices). On the other side, the “free” parameter $\sigma$ in (13) could be exploited to accelerate the convergence speed.

The possibility of using at each iterate an approximation of $p({x}^{(k)};h_{{\sigma}})$ makes the method well suited for the solution of a wide variety of structured problems.

The numerical results on a large scale convex problems shows that the performances of the inexact method are promising and comparable with those of a state-of-the-art method.

Future work will be addressed especially to deepen the theoretical and numerical analysis in the nonconvex case, investigating the possibility to obtain convergence results stronger than the ones stated in Theorems 4.1 and 4.2, at least for some classes of nonconvex functions (e.g. Kurdyka-Łojasiewicz functions).

References