Variable Metric Forward-Backward Algorithm for Composite Minimization Problems
Audrey Repetti, Yves Wiaux
TL;DR
This paper introduces a variable metric forward-backward algorithm for minimizing sums of differentiable and nonconvex nonsmooth functions, especially when the nonsmooth part is a composite of a concave and convex function, with proven convergence and practical advantages.
Contribution
It proposes a novel C2FB algorithm that handles composite nonconvex functions without explicit proximity computation, extending reweighting methods with convergence guarantees.
Findings
Requires fewer iterations to converge in image processing tasks.
Achieves better critical points compared to traditional methods.
Converges even with inexact proximity operator computations.
Abstract
We present a forward-backward-based algorithm to minimize a sum of a differentiable function and a nonsmooth function, both being possibly nonconvex. The main contribution of this work is to consider the challenging case where the nonsmooth function corresponds to a sum of non-convex functions, resulting from composition between a strictly increasing, concave, differentiable function and a convex nonsmooth function. The proposed variable metric Composite Function Forward-Backward algorithm (C2FB) circumvents the explicit, and often challenging, computation of the proximity operator of the composite functions through a majorize-minimize approach. Precisely, each composite function is majorized using a linear approximation of the differentiable function, which allows one to apply the proximity step only to the sum of the nonsmooth functions. We prove the convergence of the algorithm…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Variable Metric Forward-Backward Algorithm
for Composite Minimization Problems
Audrey Repetti*†⋆* and Yves Wiaux⋆
† School of Mathematics and Computer Sciences, Heriot-Watt University, Edinburgh, UK
⋆ School of Engineering and Physical Sciences, Heriot-Watt University, Edinburgh, UK
{a.repetti, y.wiaux}@hw.ac.uk
Abstract
We present a forward-backward-based algorithm to minimize a sum of a differentiable function and a nonsmooth function, both being possibly nonconvex. The main contribution of this work is to consider the challenging case where the nonsmooth function corresponds to a sum of non-convex functions, resulting from composition between a strictly increasing, concave, differentiable function and a convex nonsmooth function. The proposed variable metric Composite Function Forward-Backward algorithm (C2FB) circumvents the explicit, and often challenging, computation of the proximity operator of the composite functions through a majorize-minimize approach. Precisely, each composite function is majorized using a linear approximation of the differentiable function, which allows one to apply the proximity step only to the sum of the nonsmooth functions. We prove the convergence of the algorithm iterates to a critical point of the objective function leveraging the Kurdyka-Łojasiewicz inequality. The convergence is guaranteed even if the proximity operators are computed inexactly, considering relative errors. We show that the proposed approach is a generalization of reweighting methods, with convergence guarantees. In particular, applied to the log-sum function, our algorithm reduces to a generalized version of the celebrated reweighted method. Finally, we show through simulations on an image processing problem that the proposed C2FB algorithm necessitates less iterations to converge and leads to better critical points compared with traditional reweighting methods and classic forward-backward algorithms.
Keywords. Nonconvex optimization, nonsmooth optimization, proximity operator, composite minimization problem, majorize-minimize method, forward-backward algorithm, reweighting algorithm, inverse problems
MSC. 90C26, 90C59, 65K10, 49M27, 68W25, 68U10, 94A08
1 Introduction
In this work we consider optimization problems of the form
[TABLE]
where is a Lipschitz differentiable function with constant , and is a sum of composite functions as follows
[TABLE]
Composite optimization problems of this form have been studied extensively during the last decades, notably in, e.g., [14, 17, 27, 29, 36, 47, 48, 57, 59]. In general, all these works rely on the same strategy, consisting in iteratively minimizing approximations to the objective function .
1.1 Related work
A first notable example is when , is the identity function and is a proper, lower semi-continuous function whose proximity operator can be computed. In this context, a common approach to solve (1.1) is the forward-backward (FB) algorithm [18, 25, 38, 54], which alternates between a gradient step on and a proximity step on . Precisely, at each iteration , given the current iterate , the next iterate is defined as
[TABLE]
where denotes the proximity operator111The definition of the proximity operator and all the other mathematical definitions and notation used throughout this paper will be given in section 2. of , and . The convergence of the iterates of the FB algorithm to a minimizer of has been established when both and are convex, and choosing (see e.g. [13, 25]). This result has been extended in [3] to the case when both and are nonconvex. Assuming that , the authors have proved the convergence of to a critical point of , using the Kurdyka-Łojasiewicz (KL) inequality [7, 35, 39]. However, as many first-order minimization methods, it may suffer from slow convergence [18]. Accelerated versions of the FB algorithm have been proposed in the literature, mainly based either on Nesterov’ accelerations [4, 37] (also called subspace accelerations, or FISTA) or on preconditioning strategies [20, 21, 30]. While the former uses information from the previous iterates, the second aims to improve the step-size at each iteration by introducing a symmetric positive definite (SPD) matrix , leading to the following variable metric forward-backward algorithm (VMFB):
[TABLE]
When is chosen to be equal to the identity matrix , then the basic FB algorithm is recovered. However, as shown in [20, 49], wiser choices of can drastically accelerate the convergence of the iterates, and even outperform FISTA. In particular, in [20], the authors proposed to choose the preconditioning matrices using a majorize-minimize (MM) approach [33, 40, 53, 58]. Indeed, one can notice that (1.4) can be equivalently rewritten as an MM algorithm:
[TABLE]
The above equation also shows that the VMFB (and FB) algorithm can be interpreted as the proximal regularization of linearized at the current iterate [4, 9]. When , using the descent lemma, it is straightforward to notice that is a majorant function of at in the sense that, for every , and . In [20], the authors proposed to choose to define a more accurate majorant function at each iteration, and proved the convergence of sequences generated by (1.4) to a critical point of using the KL inequality (when is assumed to be convex). However, to the best of our knowledge, there exists no version of the FB algorithm with convergence guarantees able to solve the general composite problem (1.1)-(1.2) for more general choices of , in particular when the proximity operator of cannot be computed.
During the last years, many optimization methods, mainly based on MM strategies, arose to solve the general composite problem in the case where , , and , i.e.
[TABLE]
In particular, we can distinguish two main approaches. The first one consists in majorizing the outer function in (1.6). This approach has been investigated in [31, 44]. Precisely, in [44], the authors assume that where and, for every , is convex, and is coordinate-wise non-decreasing. They also do not require to be differentiable, but to be proper, lower semi-continuous and convex. In this context, the authors propose to solve (1.1)-(1.6) by defining
[TABLE]
where is a convex, proper, component-wise non-decreasing majorant function of at . Using the KL inequality, the authors show that converges to a critical point of . However, problem (1.7) needs to be solved exactly at each iteration, which may not be possible in practice. Note that convergence of general MM algorithms in a non-convex setting (not necessarily for composite functions) has also been investigated in [8], but it necessitates as well problem (1.7) to be solved exactly at each iteration to ensure the convergence of . In [31], the authors propose a similar approach, using the same update as in (1.7), but under different assumptions, and with a particular form for the majorant function . Precisely, for every , the authors propose to choose the majorant function as follows
[TABLE]
where and is the Bregman distance relative to a Legendre function [51]. In particular, when is the usual Euclidean norm squared, then (1.8) corresponds to a quadratic majorant of at . In this work, the authors do not assume to be differentiable, but only proper and lower semi-continuous. The function is assumed to be continuously differentiable, and to be differentiable such that is convex for some constant , and is locally Lipschitz continuous on (for strongly convex given by the Bregman distance in (1.8)). Under these technical assumptions and using the KL inequality, the convergence of is then guaranteed. However, both functions and necessitate to be differentiable, and problem (1.7) must be solved accurately at each iteration (often using sub-iterations). The second approach to solve the full composite problem (1.1)-(1.6) consists in using Taylor-like models, and has been investigated in [28, 45]. In both the works, the authors propose to investigate a general problem, consisting in minimizing the function on . In particular, in [28], the authors propose to define the next iterate as
[TABLE]
where is a model function for at , in the sense that there exists a growth function such that, for every , . In [45], the authors propose a modified version of (1.9), adding a Bregman distance to that changes at each iteration. In both the cases, sub-iterations are needed at each iteration. In addition, both the two works have similar convergence guarantees, obtained without the need of the KL inequality. It is interesting to note that, in [45], the authors pointed out that, at each iteration, the function must be solved accurately to reach asymptotic convergence.
It is worth mentioning that other types of methods have been proposed to minimise composite functions. In particular one can mention the Gauss-Newton methods [46]. These methods have been originally proposed to minimize functionals of the particular fom of , by solving iteratively approximated problems where is linearized. This approach can also be used to minimize the sum of and an additional term, however traditional Gauss-Newton methods may be unstable when the additional term is non-smooth [55]. To circumvent this issue, a Gaussian-Newton method with proximal linearization have been proposed in [34]. Nevertheless this approach cannot handle more general composite functions than choosing . Note that Gauss-Newton methods have also been used to minimize more general composite functions (see e.g. [15]) for convex and continuously differentiable, but without additional term.
In the current work, we propose a novel algorithm merging the structure of the VMFB algorithm with the MM strategy consisting in approximating the composite function of interest.
1.2 Proposed approach
In this work, we will aim to solve problem (1.1)-(1.2) where, for every , is convex, proper, lower semi-continuous and Lipschitz continuous on its domain, and is a concave, strictly increasing and differentiable function, such that is Lipschitz-continuous on the domain of , where denotes the first derivative of .
To solve problem (1.1), we propose to use a VMFB-based algorithm. Basically, as in (1.4), the Lipschitz differentiable function is handled through a gradient step, while the non-smooth term is handled using a proximal step. However, due to the composite form of in (1.2), its the proximity operator might not be computable, either efficiently, or at all. To overcome this difficulty, we propose to replace at each iteration , the function by an approximation denoted by . Precisely, this approximation is chosen to be a majorant function of at :
[TABLE]
such that, for every , is a majorant function of at , i.e.
[TABLE]
and is obtained by taking, for every , the tangent of the concave differentiable function at :
[TABLE]
Then, at each iteration , the proposed VMFB algorithm with approximated proximal step reads
[TABLE]
where , and is an SPD matrix. As recalled in (1.5), the classic FB (and VMFB) algorithm to solve (1.1) can be seen as the proximal regularization of linearized at the current iterate . Interestingly, in the proposed algorithm (1.13), the second function is also (partially) linearized through the linearization of the functions .
To avoid computing the approximated function at each iteration, we propose to fix it for a given finite number of iterations. Precisely, at each iteration , the majorant function is computed using the current iterate , and kept fixed for VMFB iterations, applied to . Then, the proposed variable metric composite function forward-backward (C2FB) algorithm to solve problem (1.1)-(1.2) is given by
[TABLE]
where, for every , and every , , and is an SPD matrix. One can notice that in the case when for every , algorithm (1.14) can be interpreted as an MM algorithm of the same flavour as those proposed in [28, 31, 44, 45]. Indeed, in this case, at each iteration , we have
[TABLE]
In this work, we prove the convergence of sequences generated by the C2FB algorithm given in (1.14) to a critical point of , using the KL inequality. Subsequently, according to the remark above, we show that algorithm (1.15) converges to a critical point of , if each sub-problem is solved using VMFB iterates, independently of the number of iterations towards solving each sub-problem.
The remainder of the paper is organized as follows. We introduce our notation and give useful definitions of non-convex optimization in section 2. The proposed C2FB method, including an inexact version allowing the proximity operator to be computed inexactly, are given in section 3. In section 4, we investigate the asymptotic behaviour of the proposed method, and we give the main convergence result of this work. Particular cases of the proposed approach, including reweighting algorithms, are described in section 5. Finally, simulation results on a small image restoration problem are provided in section 6.
2 Optimization background
In this section we give the definitions and notation used throughout the paper. For additional definitions on non-convex optimization, we refer the reader, e.g., to [52].
2.1 Analysis notation
Definition 2.1
Let .
- (i)
The level set of at height is defined as . 2. (ii)
The domain of is defined as . 3. (iii)
The function is proper if its domain is non-empty.
Definition 2.2
Let be a subset of and . The distance from to is defined by . If , then .
Definition 2.3
Let and be two SPD matrices. We denote by the Loewner partial ordering on , defined as, for every , . The weighted norm associated with is defined as, for every , .
2.2 Proximity operator
The proximity operator, relative to a metric, is defined as follows (see e.g. [32, Sec. XV.4] and [3, 20, 24]).
Definition 2.4
Let be a proper, lower-semicontinuous function. Let be an SPD matrix, and let . The proximity operator at of relative to the metric induced by is given by .
Remark 2.5
- (i)
In the definition of the proximity operator, since is coercive and is proper and lower-semicontinuous, if is bounded from below by an affine function, then, for every , is a non-empty set. 2. (ii)
If is convex, for every , is unique. In addition, if , then is the proximity operator originally defined in [42].
2.3 Sub-gradients
Definition 2.6
Let and . The Fréchet sub-differential of at is denoted by , and is given by
[TABLE]
If , then .
The limiting sub-differential of at is denoted by , and is given by
[TABLE]
Remark 2.7
- (i)
An equivalent definition of (2.1) is given by [52, Def. 8.3]: \widehat{\partial}f(\overline{x})=\Big{\{}\widehat{v}(\overline{x})\in\mathbb{R}^{N}|(\forall x\in\mathbb{R}^{N})\;f(x)\geqslant f(\overline{x})+{\left\langle{{\widehat{v}(\overline{x})}\mid{x-\overline{x}}}\right\rangle}+o(|x-\overline{x}|)\Big{\}} 2. (ii)
A necessary condition for to be a minimizer of is that is a critical point of , i.e. . If is convex, this condition is also sufficient.
The next results are basic chain rules that will be used throughout the paper.
Proposition 2.8** **(Basic chain rules [52])
- (i)
Let be a differentiable function and , then we have . 2. (ii)
*Let, for every , be a convex, proper and lower-semicontinuous function. Then we have \partial\Big{(}\sum_{p=1}^{P}\psi_{p}\Big{)}=\sum_{p=1}^{P}\partial\psi_{p}. *
The following closedness property (see e.g [52, Thm. 8.6]) will be useful to establish the convergence of the proposed C2FB algorithm to a critical point of . This property is often used to investigate convergence of sequences in a nonconvex setting (see e.g. [3, 10, 11, 43]).
Theorem 2.9
Let \big{(}x_{k},t(x_{k})\big{)}_{k\in\mathbb{N}} be a sequence belonging to . If \big{(}x_{k},t(x_{k})\big{)}_{k\in\mathbb{N}} converges to \big{(}x^{\star},t^{\star}\big{)}, and converges to , then \big{(}x^{\star},t^{\star}\big{)}\in\text{\rm graph}\,\partial f.
The following proposition will be used to link the sub-gradients of the objective function and its majorant. This result is similar to the one presented in [44, Lemma 1], with slightly different conditions (in particular, we do not assume that the outer function has locally Lipschitz continuous gradient). To the best of our knowledge, this sub-gradient calculation rule has not been presented elsewhere.
Proposition 2.10
Let be a proper function which is continuous on its domain, and let be a concave, strictly increasing and differentiable function. We further assume that is continuous on its domain. Then, for every , we have .
Proof. Firstly, let us show the first inclusion, i.e., , for . Let v(\overline{x})\in\partial\big{(}\phi\circ\psi\big{)}(\overline{x}). According to Definition 2.6, there exists \big{(}x_{k},v(x_{k})\big{)}_{k\in\mathbb{N}} converging to \big{(}\overline{x},v(\overline{x})\big{)}, such that and, for every , v(x_{k})\in\widehat{\partial}\big{(}\phi\circ\psi\big{)}(x_{k}). Thus, according to Remark 2.7(i), for every , we have
[TABLE]
Since is a concave and differentiable function, we have, for every , . Let and , then
[TABLE]
Combining the last inequality with (2.3) leads to
[TABLE]
Since is a strictly increasing function, for every , . Then , and the last inequality is equivalent to
[TABLE]
Then, by definition of , we have \big{(}(\phi^{\prime}\circ\psi)(x_{k})\big{)}^{-1}v(x_{k})\in\widehat{\partial}\psi(x_{k}). In addition, since and is a continuous function, we have . Finally, since , is continuous on , and , we have \big{(}(\phi^{\prime}\circ\psi)(x_{k})\big{)}^{-1}v(x_{k})\to\big{(}(\phi^{\prime}\circ\psi)(\overline{x})\big{)}^{-1}v(\overline{x}). Therefore, using Definition 2.6, we can conclude that \big{(}(\phi^{\prime}\circ\psi(\overline{x})\big{)}^{-1}v(\overline{x})\in\partial\psi(\overline{x}), i.e. .
We will now show the second inclusion, i.e. , for . Let , with . Then, according to Definition 2.6, there exists a sequence \big{(}x_{k},r(x_{k})\big{)}_{k\in\mathbb{N}} converging to , such that and . According to [52, Prop. 8.5], for every , on a neighbourhood of , there exists a differentiable function such that , , and , for every , . Let . Since is differentiable, the function is also differentiable, and we have
[TABLE]
In addition, we have and, since is a continuous and strictly increasing function, for every , . Then, using again [52, Prop. 8.5], we deduce that . In addition, since and is a continuous function, we have . Finally, since and is continuous, we have . In addition, since , we have . According to Definition 2.6, we can deduce that .
3 Proposed optimization method and assumptions
Before giving the assumptions necessary to prove the convergence of the proposed method, we would emphasize that the C2FB algorithm described in (1.14) can be rewritten using the proximity operator of instead of the proximity operator of .
Let, for every and ,
[TABLE]
Then, according to (1.12), for every , we have , where . Consequently, algorithm (1.14) can be reformulated as follows:
[TABLE]
Using the definition of the proximity operator, we can deduce that, for every and we have .
We can observe two particular cases of algorithm (3.2). On the one hand, in the particular case when, for every , , then algorithm (3.2) reads
[TABLE]
where is given by (3.1). Algorithm (3.3) requires to redefine the majorant function at each iteration , while in algorithm (3.2) the majorant function is fixed for a finite number of iterations . On the other hand, as emphasized in the introduction, under technical assumptions, in the limit case when, for every , , according to [20] the sequence converges to a critical point of . In other words, each inner-loop in algorithm (3.2) corresponds to the VMFB algorithm as defined in [20], for minimizing , leading to algorithm (1.15).
3.1 Assumptions
In the remainder of this work, we will focus on functions and satisfying the following assumptions. Examples of functions satisfying the needed assumptions are described in section 5.
Assumption 3.1
- (i)
The function is differentiable, and has a -Lipschitzian gradient, with , i.e., for every , . 2. (ii)
For every , is a convex, proper and lower-semicontinuous function. Moreover, it is Lipschitz-continuous on its domain. 3. (iii)
For every , the function is a concave and strictly increasing function (i.e. for every ). Moreover, it is differentiable on . 4. (iv)
For every , is locally Lipschitz continuous on its domain. 5. (v)
The function is coercive, i.e. .
The remark below provides comments on 3.1.
Remark 3.2
- (i)
According to 3.1(i)-(iii), is continuous on . In addition, according to 3.1(v), we can deduce that, for every , is compact ([22, Prop. 11.12]). 2. (ii)
According to 3.1(ii), there exists such that, for every , , for every , with . This assumption is satisfied for simple choices of (see section 5 for examples). 3. (iii)
According to 3.1(ii)-(iii), is continuous on its domain. 4. (iv)
According to 3.1(ii)-(iii) remark 3.2-(iii), for every and , the parameter introduced in (3.1) is well defined. 5. (v)
3.1(iv) holds if, for every , is on its domain. 6. (vi)
3.1(iv) holds if and only if, for every , the function is Lipschitz continuous on every compact subset of . Thus, under 3.1(iv), there exists such that, for every , and (x^{\prime},x^{\prime\prime})\in\big{(}\text{lev}_{\leqslant f(x)}f\big{)}^{2}, .
For every , the SPD matrices are used in practice to accelerate the convergence of usual FB methods. They are chosen using the method proposed in [20, 21, 30], leveraging an MM approach. We define them as follows:
Assumption 3.3
Let, for every , be generated by algorithm (1.14).
- (i)
For every , we have
[TABLE] 2. (ii)
There exists such that, for every and , .
Remark 3.4
According to [20, Lem. 2.1], under Assumption 3.1(i), Assumption 3.3 is trivially satisfied when choosing, for every and for every , and .
The two last assumptions are made to ensure that the step-sizes are bounded, and the inner-iteration numbers are finite.
Assumption 3.5
There exists such that, for every and for every we have .
Assumption 3.6
There exists such that, for every , .
3.2 Inexact algorithm
The proximity operator relative to an arbitrary metric may not have a closed form expression. This is also true for some evolved functions, even when the preconditioning operators are diagonal matrices. For instance when, for every , is a composition between an -norm and a non-orthogonal matrix (e.g. to promote sparsity in a redundant dictionary), the computation of the proximity operator of is done iteratively [23]. To circumvent this difficulty, we propose an inexact version of the proposed C2FB algorithm given in (3.2):
[TABLE]
In algorithm (3.5), for every and , is an SDP matrix satisfying 3.3, and (see proposition 2.8). Thus, finding is equivalent to finding, for every , .
Under our assumptions, algorithm (3.5) can be viewed as an inexact version of algorithm (3.2) (or equivalently algorithm (1.14)), where, at each iteration , the proximity operator of (or equivalently ) can be computed inexactly (i.e. using sub-iterations). The first inequality in algorithm (3.5) is called sufficient-decrease condition, while the second inequality is the inexact optimality condition. These two conditions allow to handle possible errors arising when the proximity operator is computed approximately, and are often referred to relative errors. Although they are common for FB algorithms in a nonconvex context (see e.g. [3, 20, 21]), they are more of theoretical interest, showing that the algorithm is robust with respect to inexact computations, than of practical use. Note that implementable inexactness conditions for the FB algorithm have been proposed in other works (see e.g. [56] in the convex context and [11] in the non-convex context).
We now show formally that, in particular, sequence and, for every , generated by algorithm (3.2) satisfy the sufficient decrease condition and the inexact optimality condition of algorithm (3.5). This will show that algorithm (3.5) can be viewed as an inexact version of algorithm (3.2).
Let and . Using definition 2.4, we have
[TABLE]
According to 3.5 the first condition in algorithm (3.5) (i.e. sufficient-decrease condition) is obtained with .
The second condition in algorithm (3.5) (i.e. inexact optimality condition) is obtained combining algorithm (3.2) with 3.3(ii) and 3.5. Indeed, using the variational characterization of the proximity operator, we have, for every and ,
[TABLE]
Therefore, there exists such that v_{k}(\widetilde{x}_{k,i+1})=\gamma_{k,i}^{-1}A_{k,i}\big{(}\widetilde{x}_{k,i}-\widetilde{x}_{k,i+1}\big{)}-\nabla h(\widetilde{x}_{k,i}). Then, using 3.3(ii) and 3.5, we obtain
[TABLE]
and the inexact optimality condition is obtained with .
4 Convergence analysis
4.1 Descent properties
In this section, we provide convergence properties on for generated by algorithm (3.5). We also investigate the behaviour of the approximate objective function at the current iteration , defined by
[TABLE]
As a first step, for fixed, we focus on the sub-iterations and we investigate the behaviour of . In particular, in the following lemma we show that is non-increasing.
Lemma 4.1
Let and, for every , be generated by algorithm (3.5). Let and . Under Assumptions 3.1, 3.3, 3.5 and 3.6, we have
[TABLE]
where , and hence
[TABLE]
Proof. Let and, for every , be generated by algorithm (3.5). Using the sufficient decrease condition in algorithm (3.5) (i.e. the first inequality), we have, for every and ,
[TABLE]
Let . Summing on , we obtain
[TABLE]
[TABLE]
Summing the last inequality on , we obtain
[TABLE]
which, combined with (4.5), leads to
[TABLE]
Then, using 3.3(ii), there exists such that (4.2) is satisfied.
Equation (4.3) is then obtained by adding the term \sum_{p=1}^{P}\Big{(}(\phi_{p}\circ\psi_{p})(x_{k})-\lambda_{p,k}\psi_{p}(x_{k})\Big{)}, constant over , on both sides of (4.8), and using the definitions of function (see (1.10)-(1.12)), of function (see (3.1)) and of function (see (4.1)).
In the following proposition, we investigate the behaviour of the global sequences \big{(}f(x_{k})\big{)}_{k\in\mathbb{N}} and \big{(}\|x_{k+1}-x_{k}\|^{2}\big{)}_{k\in\mathbb{N}}, where is generated by algorithm (3.5).
Proposition 4.2
Let and, for every , be generated by algorithm (3.5). Under Assumptions 3.1, 3.3, 3.5 and 3.6, the following assertions hold:
- (i)
For every , we have
[TABLE]
where is given in lemma 4.1, and .
In addition, we have
[TABLE] 2. (ii)
\big{(}f(x_{k})\big{)}_{k\in\mathbb{N}}* is a converging non-increasing sequence.* 3. (iii)
We have and , and thus and .
Proof.
- (i)
The first inequality (4.9) is a direct consequence of lemma 4.1. On the one hand, the right-sided inequality in (4.9) is obtained by choosing , and in (4.3). Indeed, noticing that , and that , we have:
[TABLE]
On the other hand, the left-sided inequality in (4.9) is obtained using the majoration condition of function (see (1.10)):
[TABLE]
The second inequality (4.10) is obtained using Jensen’s inequality, as
[TABLE]
The result is then obtained using 3.6. 2. (ii)
We deduce directly from (4.9) that \big{(}f(x_{k})\big{)}_{k\in\mathbb{N}} is a non-increasing sequence. In addition, combining the fact that, for every , (according to the definition of the proximity operator) with remark 3.2(i), we deduce that belongs to a compact subset of . Thus, being lower bounded (according to 3.1 is continuous on its domain and coercive), \big{(}f(x_{k})\big{)}_{k\in\mathbb{N}} converges to a real number . 3. (iii)
According to eq. 4.9, we have \|\chi_{k}\|^{2}\leqslant\overline{\alpha}^{-1}\Big{(}f(x_{k})-f(x_{k+1})\Big{)}. Let be a positive integer. Summing the last inequality from to we have \sum_{k=0}^{K-1}\|\chi_{k}\|^{2}\leqslant\overline{\alpha}^{-1}\Big{(}f(x_{0})-f(x_{K})\Big{)}. Since is a non-increasing sequence converging to , we thus obtain \sum_{k=0}^{K-1}\|\chi_{k}\|^{2}\leqslant\overline{\alpha}^{-1}\Big{(}f(x_{0})-\xi\Big{)}. The result is then obtained by taking the limit for .
The same arguments can be applied for , leveraging (4.10).
Finally, the next proposition will be useful to show that the sequences generated by algorithm (3.5) are approaching the set of critical points of the global objective function .
Proposition 4.3
Let and, for every , be generated by algorithm (3.5). Under Assumptions 3.1, 3.3, 3.5 and 3.6, we have, for every ,
[TABLE]
where , is defined in proposition 4.2, and
[TABLE]
with, for every , provided by algorithm (3.5).
Proof. Firstly, using the chain rule given in proposition 2.8(i), we have . In addition, using proposition 2.10 and proposition 2.8(ii), we notice that where . Thus, we can deduce that as defined in (4.15) belongs to .
We now show inequality (4.14). Let be defined as in algorithm (3.5), i.e. , where, for every , . Using the triangular inequality and remark 3.2(ii), we have
[TABLE]
Let, for every , be defined as in algorithm (3.5). Then, using the fact that and Jensen’s inequality, we have
[TABLE]
where the last majoration is obtained using 3.1(i), 3.3(ii) and the second inequality condition in algorithm (3.5). Then, by definition of , we obtain
[TABLE]
In addition, according to 3.1(iv) and remark 3.2(vi), since belongs to the compact set , for every , we have . Thus, combining this last inequality with (4.16) and (4.18), we obtain . Noticing that , we obtain \|t(x_{k+1})\|\leqslant\big{(}\sqrt{2(\mu^{2}+\overline{\nu}\beta^{2})}+\upsilon\eta\big{)}\|\chi_{k}\|.
4.2 Convergence results
Before giving our main convergence result, we need to introduce our last assumption, concerning the Kurdyka-Łojasiewicz inequality [2, 5, 6, 7, 35, 39]. To this aim, we introduce an additional notation: for , we denote by the class of concave and continuous functions such that , is on , and for every .
Assumption 4.4
The function satisfies the Kurdyka-Łojasiewicz inequality, i.e., for every , there exist , a neighbourhood of , and a function , such that for every satisfying , we have
[TABLE]
As emphasized, e.g., in [2, 3, 9], the KL inequality is satisfied for a wide class of functions, and in particular by sub-analytic, log-exp ans semi-algebraic functions222A function is semi-algebraic if its graph is a finite union and intersection of sets defined by a finite number of polynomial equalities and inequalities..
In order to use 4.4 to investigate the convergence of sequences generated by algorithm (3.5), we need to use the following property, initially introduced in [1, Lemma 1] and [9, Lemma 6].
Lemma 4.5** **(Uniformized KL property)
Let be a compact subset of . Let be a proper and lower-semicontinuous function, constant on and satisfying the KL inequality on . Then, there exists , and such that, for every , and for every satisfying
[TABLE]
we have \varphi^{\prime}\big{(}f(x)-f(\overline{x})\big{)}\operatorname{dist}\big{(}0,\partial f(x)\big{)}\geqslant 1.
We will now investigate the properties of the limit point set. Let be a sequence generated by algorithm (3.5), and be a starting point. The associated set of limit points is defined as \Omega_{x_{0}}=\big{\{}\overline{x}\in\mathbb{R}^{N}\,|\,\exists\text{ an increasing sequence of} \text{ integers }(k_{j})_{k\in\mathbb{N}}\text{ such that }x_{k_{j}}\underset{j\to\infty}{\to}\overline{x}\big{\}}.
Remark 4.6
According to [9, Lemma 5] and proposition 4.2(iii), we have
- (i)
\displaystyle\lim_{k\to\infty}\text{dist }\Big{(}x_{k},\Omega_{x_{0}}\Big{)}=0, 2. (ii)
is a non-empty, compact and connected set.
In addition, following the same arguments as in [9, Lemma 5], we have the following results.
Lemma 4.7** **(Properties of the limit point set)
Let be a sequence generated by algorithm (3.5) with starting point . The following assertions hold under Assumptions 3.1, 3.3, 3.5 and 3.6.
- (i)
If , then is a critical point of , i.e. . 2. (ii)
The objective function is finite and constant on .
Proof.
- (i)
Let be a limit point of . On the one hand, since for every , , and is continuous on (see remark 3.2(i)), there exists a subsequence converging to such that as . In addition, since is a convergent and non-increasing sequence (according to proposition 4.2(ii)), we have as . On the other hand, using proposition 4.2(iii) and proposition 4.3, we have and as . So, using the closedness property of given in theorem 2.9, we have , and thus is a critical point of . 2. (ii)
On the one hand, according to proposition 4.2(ii), there exist such that converges to as . On the other hand, let , i.e. there exists a subsequence of converging to as . As seen in the proof of statement (i), we have as , thus . Therefore, the restriction of to is constant, and equal to .
The next theorem is our main convergence result, analysing the convergence of the sequences generated by the C2FB algorithm given in (3.5).
Theorem 4.8
Let be a sequence generated by algorithm (3.5). Under Assumptions 3.1, 3.3, 3.5, 3.6 and 4.4, the following holds.
- (i)
* has a finite length, i.e. .* 2. (ii)
* converges to a critical point of . *
Proof.
- (i)
First, in order to apply the result given in lemma 4.5, we need to show that there exists such that, for every , satisfies condition (4.20) for . According to remark 4.6(ii), is non-empty. So there exists a convergent subsequence , i.e. there exists such that as . In addition, using the same arguments as in the proof of lemma 4.7, we can deduce that . On the one hand, being a convergent and a non-increasing sequence, we have, for every ,
[TABLE]
On the other hand, according to remark 4.6(i), for every , there exists such that, for every , \operatorname{dist}\big{(}x_{k},\Omega_{x_{0}}\big{)}<\kappa. So, there exists such that, for every , satisfies condition (4.20). According to remark 4.6(ii), is compact, and using lemma 4.7(ii), is constant on , thus we can apply lemma 4.5. Then, for every , we have
[TABLE]
Combining this inequality with proposition 4.3, we obtain
[TABLE]
In addition, since is concave, we have, for every , . Then, by taking and , we obtain
[TABLE]
Let \Delta_{k}=\varphi\Big{(}f(x_{k})-f(x^{\star})\Big{)}-\varphi\Big{(}f(x_{k+1})-f(x^{\star})\Big{)}. Combining the last two inequalities leads to . In addition, using proposition 4.2(i), we obtain . Thus,
[TABLE]
where the last majoration is obtained using the fact that, for every , .
Summing the last inequality from to , we obtain
[TABLE]
In addition, we have
[TABLE]
where the last inequality is obtained using the fact that, for every , . Therefore, combining ((i)) and ((i)), we obtain \sum_{k=k^{\star}+1}^{K}\|\chi_{k}\|\leqslant\overline{\alpha}^{-1}\overline{\beta}\varphi\Big{(}f(x_{k^{*}+1})-f(x^{\star})\Big{)}+\|\chi_{k^{*}}\|. Since , , using Jensen’s inequality as in ((i)) and 3.6 we have . So we conclude that has finite length. 2. (ii)
Since has finite length, it is a Cauchy sequence, hence converging to a point . Then the result follows directly from lemma 4.7(i).
Remark 4.9
The proof of theorem 4.8 is very similar to the one of [9, Theorem 1]. However, the later cannot directly be applied to algorithm (3.5) due to the sub-iterations on , where the majorant function of is fixed. Precisely, in (4.23), the bound used on the sub-gradient of at depends on the sub-iterates (see proposition 4.3). This prevents the use of the proof of [9, Theorem 1] that only relies on the main sequence .
5 Particular cases of the proposed method
In this section we describe particular cases for the proposed C2FB algorithm, paying attention to the assumptions necessary to ensure that the convergence of Theorem 4.8 holds.
5.1 Variable Metric Forward Backward algorithm
The proposed algorithm boils down to the VMFB developed in [20] when choosing, for every , , , , and being a proper and Lipschitz-continuous function. In this case, we have and the global minimization problem is of the form
[TABLE]
In this particular case, we have, for every , , and consequently, algorithms (1.14) reduces to
[TABLE]
We would emphasize that the assumptions on to ensure convergence of (5.2) are different than in [20]. On the one hand, in [20] the function only needs to be continuous on its domain (while in Assumption (ii) we assume that is Lipschitz-continuous on its domain). Other works present convergence results for algorithm (5.2) when is non-convex either in the case without variable metric [3] or for alternating minimization [21]. On the other hand, in [20] the step-sizes are chosen such that, there exists such that, for every , , while in Assumption 3.5 it is assumed that .
5.2 Reweighted algorithms
The second example is interesting in computational imaging since it is related to iteratively reweighted algorithms [44].
We consider problem (1.1)-(1.2) with, for every , satisfying Assumption 3.1 and, for every , , where and . In this case, we have
[TABLE]
and, for every and , q_{p}(x,x_{k})=\theta\log(\psi_{p}(x_{k})+\varepsilon)+\lambda_{p,k}\big{(}\psi_{p}(x)-\psi_{p}(x_{k})\big{)}, with \lambda_{p,k}=\theta{\big{(}\psi_{p}(x_{k})+\varepsilon\big{)}^{-1}}.
In this context, the proposed C2FB algorithm given in (3.2) reduces to
[TABLE]
5.2.1 Reweighted algorithm for log-sum penalization
Let us consider the following particular case. For every ,
[TABLE]
where is a linear operator (e.g. wavelet transform [41]), and denotes the -th component of its argument. Then we have, for every , \lambda_{p,k}={\theta}{\big{(}|[Wx_{k}]^{(p)}|+\varepsilon\big{)}^{-1}}. In this context, the function corresponds to a log-sum penalization composed with a linear operator:
[TABLE]
and algorithm (5.4) reduces to the re-weighted algorithm, initially proposed in [16] for , where each sub-problem is solved using a VMFB algorithm. Let be the operator giving the diagonal matrix whose diagonal elements are given by its argument. The resulting algorithm reads
[TABLE]
Note that (5.6) with is known as the log-sum penalization and its proximity operator has an explicit form [19]. Let A=\operatorname{Diag}\big{(}(a^{(n)})_{1\leqslant n\leqslant N}\big{)}\in[0,+\infty[^{N\times N} be a diagonal SDP matrix. Then we have
[TABLE]
Consequently, in this context, one can apply the classic VMFB algorithm, for diagonal matrices , directly to minimize .
5.2.2 Cauchy penalization
Similarly to the reweighting algorithm described above, another particular case is when choosing, for every and , \psi_{p}(x)=\big{(}[Wx]^{(p)}\big{)}^{2}, where is a linear operator. In this case, corresponds to the Cauchy penalization composed with a linear operator:
[TABLE]
and, for every and , the weights in algorithm (5.4) are given by \lambda_{p,k}=\theta{\big{(}([Wx]^{(p)})^{2}+\varepsilon\big{)}^{-1}}.
5.3 Nonconvex norms
In this section we show that our method can be used to handle nonconvex -norms, where , defined by
[TABLE]
The proximity operator (relative to the Euclidean norm) of , relative to the metric induced by the diagonal matrix A=\text{Diag}\big{(}(a^{(n)})_{1\leqslant n\leqslant N}\big{)}, has an explicit formula, however it necessitates to find roots of polynomial equations [12]:
[TABLE]
where is such that . In practice is found approximately using the Newton method [19].
The proposed approach cannot handle directly the function, and we need to introduce an approximation of it. Indeed, to define the norm as a composition of functions, we take, for every , to be the function defined in (5.5). In this case, to obtain the function defined in (5.10), we should choose, for every and , , where is a regularization parameter. However this function is not differentiable at [math], hence violating 3.1(iii). Instead we propose to choose \phi_{p}(u)=\theta\big{(}(u+\varepsilon)^{\rho}-\varepsilon^{\rho}\big{)}, where . Then, is chosen to be a slightly modified version of the nonconvex -norm, denoted by , composed with a linear operator :
[TABLE]
In this context, for every , and , we have q_{p}(x,x_{k})=\theta\Big{(}\big{(}|[Wx_{k}]^{(p)}|+\varepsilon\big{)}^{\rho}-\varepsilon^{\rho}\Big{)}+\lambda_{p,k}\Big{(}|[Wx]^{(p)}|-|[Wx_{k}]^{(p)}|\Big{)}, where \lambda_{p,k}=\theta\rho\big{(}|[Wx_{k}]^{(p)}|+\varepsilon\big{)}^{\rho-1}. Therefore, the proposed algorithm given in (3.2) reduces to
[TABLE]
6 Simulations
Many imaging problems such as reconstruction, restoration, inpainting, etc., can be formulated as inverse problems. In this context, the objective is to find an estimate of an original unknown image from degraded observations , given by , where is a linear observation operator, and is a realization of an additive independent identically distributed (i.i.d.) random noise. When the random noise is normally distributed, with zero-mean, a common approach to find is to define it as the minimizer of a penalized least-squares criterion, i.e. solve (1.1) with, for every , .
We consider a restoration example in image processing. Precisely, we choose to be the image jetplane of size shown in Figure 1, and to model a blurring operator. In this case, is implemented as a convolution operator such that the image is convolved with a motion blurring kernel of length 5 and angle . The noisy observation is then obtained building as a realization of an i.i.d. Gaussian variable with zero-mean and standard deviation . In our simulations we will consider two different noise levels, defined through the input signal-to-noise ratio (iSNR): . Precisely, we will consider the cases when iSNR dB and iSNR dB. In these cases, the SNR of the observed images are equal approximately to dB and dB, respectively. For each case, we run simulations for 50 realizations of random noise. In Figure 1 are shown the blurred image and an example of a noisy observation , when iSNR dB.
For the matrices , since the preconditioning scheme has been already widely discussed in the literature (see e.g. [20, 21, 24, 50]), we will not investigate this functionality in our simulations. In our simulations, we propose to use the same approach as proposed in [20]. In this context, the matrices are diagonal, fixed over iterations, and can be seen as diagonal approximation of the Hessian of .
We will consider two different regularization terms: the log-sum penalization and the penalization described in sections 5.2.1 and 5.3, respectively. For the two considered penalization terms we will provide reconstruction results obtained using a classic VMFB algorithm [20]333Note that the VMFB algorithm does not have convergence guaranties when considering to be a non-convex function [20], unlike the FB algorithm [3]. Nevertheless, the convergence guaranties can be deduced from [21], for standard numerical implementations of the log function. and the proposed C2FB algorithm. Since C2FB is proved to converge for any value of , we will consider different values of , keeping it fixed over iterations. As explained with equation (1.15), when is chosen large enough (i.e. when the inner-loop has converged), C2FB has a similar behaviour as state-of-the-art methods to minimize composite functions (see e.g. [28, 31, 44, 45]). For both the methods we will investigate the convergence behaviour, the reconstruction quality of the estimate, and the convergence speed. For the convergence behaviour, we will compare, for each experiment, the value of the objective function at convergence. Precisely, we will evaluate C(x^{\star}_{VMFB},x^{\star}_{C2FB}):=\Big{(}f(x^{\star}_{VMFB})-f(x^{\star}_{C2FB})\Big{)}/|f(x^{\star}_{VMFB})|, where and are the estimates obtained with VMFB and C2FB, respectively. If this criteria is positive, it means that our method reached a better critical point than VMFB, i.e. . To evaluate the reconstruction quality of the estimate, we use the signal-to-noise ratio (SNR), which is defined, for an image as \text{SNR}=10\log_{10}\Big{(}\|\overline{x}\|^{2}/\|\overline{x}-x\|^{2}\Big{)}. Finally, we will investigate the convergence speed of both the methods, in terms of number of iterations needed to reach convergence. We consider that both the algorithms have converged when the following stopping criteria are fulfilled:
[TABLE]
where is the sequence generated either by VMFB or by the proposed C2FB. It is important to emphasize that, when the computation of the proximity operator does not require sub-iterations, one iteration of VMFB has similar computational cost as one inner-iteration of C2FB. For the three above-mentioned evaluation criteria, we will show that the proposed approach leads to better results than the VMFB algorithm. In addition, choosing different values of (number of inner-iterations) will allow us to show that there is an optimal value for in terms of convergence speed and reconstruction quality, suggesting that there is no need for reaching convergence in the inner-iterations before recomputing the majorant function.
6.1 Log-sum penalization
In this section we present the simulation results obtained when solving the problem described above, using the log-sum regularization described in section 5.2.1. More precisely, we propose to
[TABLE]
where models the Db8 wavelet transform [26] with decomposition levels, and and are chosen to maximize the reconstruction quality. In our experiments we have (resp. ) for iSNR dB (resp. dB). When (6.2) is solved with the VMFB algorithm, we use the proximity operator of the log-sum function given in (5.8). When solved by our C2FB algorithm, we use algorithm (5.7), with .
Results are given in Figure 2, considering (top row) iSNR dB and (bottom row) iSNR dB. The left plots give as a function of . For the middle and right plots, the blue curves are obtained with C2FB, and the red curves are obtained considering VMFB, computing exactly the proximity operator. The middle plots give the SNR values as a function of ; and the right plots show the total number of iterations needed to reach convergence as a function of . For C2FB, the total number of iterations is given by , where is the number of outer iterations computed in algorithm (5.7) to satisfy the stopping criteria (6.1). For VMFB, the curves are constant as there is no inner-loop in the algorithm. For both the methods, the continuous lines represent the average values, and the dotted lines show the associated results within 1 standard deviation around the mean.
For both the considered noise levels, we observe that increases with and is always positive, showing that C2FB provides a better critical point. In addition, C2FB leads to better reconstruction results in terms of SNR than VMFB. When iSNR dB (resp. iSNR dB), the SNR obtained with C2FB is dB (resp. SNR dB) when (resp. ), while VMFB leads to results with SNR dB (resp. SNR dB). In addition, C2FB necessitates less global iterations to reach convergence, and we observe that there is an optimal value for in terms of total iteration number. Precisely, when iSNR dB (resp. iSNR dB), the optimal value is around (resp. ), for a total number of iterations of (resp. ). In comparison, VMFB necessitates iterations to converge in both the cases. This observation suggests that reducing the number of iterations in the inner-loop (instead of reaching convergence in each inner-loop before re-computing the weights, as suggested in classical reweighting algorithms [16]) can accelerate the convergence of the reweighting algorithm without altering the reconstruction quality.
6.2 penalization
In this section we present the simulation results obtained when solving the problem described at the beginning of the section, using the regularization described in section 5.3, with . More precisely, we propose to
[TABLE]
where is the approximation of the norm defined in (5.12), models the Db8 wavelet transform [26], with decomposition levels, and and are chosen to maximize the reconstruction quality. In our experiments we have (resp. ) for iSNR dB (resp. dB). Problem (6.3) can be solved using our C2FB algorithm, with algorithm (5.13). In our experiments, we run C2FB for different number of inner-iterations . In the limit case when , then reduces to the exact norm, and problem (6.3) can be solved using the VMFB algorithm, where the proximity operator of the norm is given by (5.11). Note that even if the proximity operator of the norm has an explicit formula, in practice it requires to compute sub-iterations.
Results are given in Figure 3, with (top row) iSNR dB and (bottom row) iSNR dB. The left plots show as a function of . For the middle and right plots, the blue curves are obtained with C2FB, and the red curves are obtained using VMFB where the proximity operator is computed using equation (5.11) combined with a Newton method. The central (resp. middle) plots show SNR values (resp. total iteration number needed to reach convergence) as a function of .
In this experiment, the conclusions are slightly different depending on the considered noise level. For the noise level corresponding to iSNR dB (top row), the criteria is decreasing when increases, and is positive for all values of , showing that the proposed C2FB algorithm provides a better critical point than VMFB independently from the number of inner-iterations. In particular, the best results are obtained when taking . This conclusion is also true when observing the middle and right plots. For the SNR (see middle plots), C2FB leads to results with SNR dB for , decreasing to SNR dB when . In comparison, the SNR for the estimate obtained with VMFB is equal to dB. On the right plots, we observe that the total iteration number needed to reach convergence increases with for C2FB, starting at for and finishing at for . For comparison, VMFB necessitates iterations to reach convergence. It is worth noticing that for this penalization function, unlike our method, the VMFB algorithm necessitates sub-iterations to compute the proximity operator. For the noise level corresponding to iSNR dB (bottom row), the criteria decreases when increases, and is positive for . This shows that the critical point obtained with C2FB is better than the one obtained with VMFB only for a small number of inner-iterations . For the reconstruction quality (see middle plots), the results obtained with C2FB have a higher SNR, equal to dB for , and decreasing to dB for . For comparison, the SNR of the estimate obtained with VMFB is equal to dB. Finally, for the computational cost (see right plots), we observe that the total number of iterations needed by C2FB to reach convergence, is decreasing for and increasing for . When (i.e. corresponding to the best reconstruction, both in terms of objective value and reconstruction quality), C2FB needs iterations to reach convergence. This total iteration number drops to for . For comparison, VMFB requires iterations to reach convergence (without counting the sub-iterations to compute the proximity operator).
We can conclude that the proposed approach provides a good alternative to the classic VMFB method, in terms of both quality reconstruction and convergence speed. In addition, the proposed C2FB algorithm outperforms as well state-of-the-art methods to minimize composite functions, obtained when (see [28, 31, 44, 45]).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] H. Attouch and J. Bolte. On the convergence of the proximal algorithm for nonsmooth functions involving analytic features. Math. Program. , 116:5–16, 2009.
- 2[2] H. Attouch, J. Bolte, P. Redont, and A. Soubeyran. Proximal alternating minimization and projection methods for nonconvex problems: an approach based on the kurdyka-łojasiewicz inequality. Math. Oper. Res. , 35:438–457, 2010.
- 3[3] H. Attouch, J. Bolte, and B. F. Svaiter. Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods. Math. Program. , 137:91–129, Feb. 2011.
- 4[4] Amir Beck and Marc Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM journal on imaging sciences , 2(1):183–202, 2009.
- 5[5] J. Bolte, A. Daniilidis, and A. Lewis. The Łojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems. SIAM J. Optim. , 17:1205–1223, 2006.
- 6[6] J. Bolte, A. Daniilidis, A. Lewis, and M. Shiota. Clarke subgradients of stratifiable functions. SIAM J. Optim. , 18(2):556–572, 2007.
- 7[7] J. Bolte, A. Daniilidis, O. Ley, and L. Mazet. Characterizations of Łojasiewicz inequalities: subgradient flows, talweg, convexity. Trans. Amer. Math. Soc. , 362(6):3319–3363, 2010.
- 8[8] J. Bolte and E. Pauwels. Majorization-minimization procedures and convergence of sqp methods for semi-algebraic and tame programs. Math. Oper. Res. , 41:442–465, 2016.