# A General Proximal Alternating Minimization Method with Application to Nonconvex Nonsmooth 1D Total Variation Denoising

Authors: Xiaoya Zhang; Tao Sun; Lizhi Cheng
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

We deal with a class of problems whose objective functions are compositions of nonconvex nonsmooth functions, which has a wide range of applications in signal/image processing. We introduce a new auxiliary variable, and an efficient general proximal alternating minimization algorithm is proposed. This method solves a class of nonconvex nonsmooth problems through alternating minimization. We give a brilliant systematic analysis to guarantee the convergence of the algorithm. Simulation results and the comparison with two other existing algorithms for 1D total variation denoising validate the efficiency of the proposed approach. The algorithm does contribute to the analysis and applications of a wide class of nonconvex nonsmooth problems.


---

## Body

## 1. Introduction

In the past few years, increasing attentions have been paid to convex optimization problems, which consist of minimizing a sum of convex or smooth functions [1–3]. Each of the objective functions enjoys several appreciated properties, like strong convexity, Lipschitz continuity, or other convex conditions. These properties can usually lead to great advantage in computing. Meanwhile works on such convex problems have provided a sound theoretic foundation. Both of the theoretic and computing advantages created many benefits to practical use. It is particularly noticeable in signal/image processing, machine learning, computer vision, and so forth. However, what deserves the special attention is the fact that the convex or smooth models are always approximations of nonconvex nonsmooth problems. For example, nonconvexl0-norm in sparse recovery problems is routinely relaxed as convexl1-norm, and many related works were developed [4, 5]. Although the difference between the nonconvex nonsmooth problem and its approximations vanishes in certain case, it is nonnegligible sometimes, like the problem in paper [6]. On the other hand, excellent numerical performances of various nonconvex nonsmooth algorithms inspire researchers to continue on their directions to the nonconvex methodology.Nonconvex and nonsmooth convex optimization problems are ubiquitous in different disciplines, including signal denoising [7], image deconvolution [8, 9], or other ill-posed inverse problems [10], to name a few. In this paper, we aim at solving the following generic nonconvex nonsmooth optimization problem, formulated in real Hilbert spaces X andUmm=1,2,…,M, for someM∈N+:(1)x^∈arg⁡minx⁡Px≔fx+∑m=1MhmLmx,where (i) f(·):X→R∪{+∞}andhm(·):Um→R∪{+∞}are proper lower semicontinuous functions; (ii) the operatorsLm:X→Umare linear; and (iii) the set of minimizers is supposed to be nonempty.It is quite meaningful to find a common convergent point in the optimal set of sums of simple functions [2, 11]. Insightful studies for nonconvex problems are presented in [12, 13]: if nonconvex structured functions of the type L=f(x)+Q(x,y)+g(y)has the Kurdyka-Lojasiewicz (K-L) property, then each bounded sequence generated by a proximal alternating minimization algorithm converges to a critical point ofL. Our main work is mainly based on this convergence result and introduces a generic proximal minimization method for the general model (1). Note that if some of the functions in (1) are zeros and linear operatorLmis the identityId, the model would reduce to common one in compressed sensing [14]. The model we consider is the generalization of many application problems, such as the common “lasso” problem [15], and composition in papers [2, 16]. Here, we provide a few of the typical examples.Example 1 (1D total variation minimization [7]).
In this application, we need to solve the following denoising problem:(2)x^∈arg⁡minx∈X⁡x-y22+λDx0,wherey∈X is the input signal, X=Rn,z0 counts the number of nonzero elements of z, and Dx≔(x2-x1,x3-x2,…,xn-xn-1)is defined as derivation of original signalx.Noise removal is the basis and requisite of other subsequential applications and algorithm dealing with total variation (TV) regularizer; this regularizer is of great importance since it can efficiently deal with noisy signals which have sparse derivatives (or gradients), for instance, piecewise constant (PWC) signal that has flat sections with a number of abrupt jumps. The 1D total variation minimization can be extended to related 2-dimension image restoration.Example 2 (group lasso [17]).
In this application, one needs to solve(3)x^∈arg⁡minx∈RN⁡12Ax-y22+λ∑i=1nxi2,wherexT=x1T,x2T,…,xnT,xi∈Rpiare the decision block variables and∑i=1npi=Nwithpiis the corresponding block size.In the past few years, researches on structural sparse signal recovery have been very popular and group lasso is typical one of those important problems. It attracts many attentions in face recognition, multiband signal processing, and other machine learning problems. The general case is also applied to many other kinds of structural sparse recovery problems, likel12-minimization [18] and block sparse recovery [19].Example 3 (image deconvolution [20]).
In this application, one needs to solve(4)X^∈arg⁡minX∈X⁡AX-YF2+λTVX,where X=Rm×nand XF=∑i∑jXi,j21/2. The discrete total variation, denoted byTVin (4), is defined as follows. We definer×rmatrix: (5)Dr≔1-1-11⋱⋱-11, and the discrete gradient operatorD:X→X2is defined as (6)DXi,j=dxi,j,dyi,j,dx=DmX,dy=XDnT. Then we have TV(X)=DX1,2=∑i∑j|dxi,j|2+|dyi,j|2.The concept of deconvolution finds lots of applications in signal processing and image processing [21–23]. In this paper, it would just be considered as a specific case to problem (1), although paying attention to this problem (4) with other details is equally important.The main difficulty in solving (1) lies in thatxis coupled byLm. In order to surmount the computation barrier, we introduce a new auxiliary variable and split the problem into two sequences of subproblems, minutely described in the next section. Then our problem is an extension of problem given in paper [12]. The paper aims at giving a generic proximal alternating minimization method for a class of nonconvex nonsmooth problem (1), to be applied in many fields. The motivation is introduction of auxiliary variable and splitting the original problem into two kinds of easier nonconvex nonsmooth subproblems. Recent studies often give the regularizationhma reasonable assumption; namely, the proximal map ofhmis easy to calculate. Then convergence analysis can be extended by the context of the present work [12]. In the last section, we show application to nonconvex nonsmooth 1D total variation signal denoising.


## 2. Algorithm

In order to reduce computation complexity caused by composite of nonconvex functionhm·and operatorLm, we introduce a sequence of auxiliary variablesθmm=1M. Then, the problem in (1) can be represented equivalently as follows: for eachm∈{1,…,M}(7)x^,θ^m∈arg⁡minx,θm⁡Pmx,θm≔fx+Qmx,θm+Hθ,whereθmm=1Mare represented as a concise wholeθand (8)Qx,θ≔∑m=1Mτm2θm-Lmx22;Qmx,θm≔τm2θm-Lmx22;Hθ≔∑m=1Mhmθm. The last term∑m=1Mτm/2θm-Lmx22 ensures the high similarity betweenθmandLmx(m=1,…,M). And this quadratic function minimization can be easily solved. Hence, the original complex composite is split into two simpler objectives.Apparently, now this form could be solved by a series of alternating proximal minimization methods [12]: for eachm∈{1,…,M},(9a)xk+1∈arg⁡minu⁡fu+Qu,θk+12ζku-xk22;(9b)θmk+1∈arg⁡minvm⁡Qmxk+1,vm+hmvm+12μkvm-θmk22.We then make the following assumptions about (9a) and (9b): HinfX×∏m=1MUm⁡fx+Hθ+Qx,θ>-∞;fx+Hθ0+Qx,θ0is  proper;for  some  positive0<r-<r+,the  stepsize  sequences  satisfyr-<ζk,μk<r+∀k≥0.The assumption is weakly required and can be easily meet. In the first section, the given three examples all satisfy the first two conditions, and the third condition holds when proper parameters are set in practical tests. Besides, each of the proximal terms1/2ζku-xk22and1/2μkvm-θmk22 in (9a) and (9b) is used to meet the decrease condition of objective function. If we omit them, the algorithm can also perform well but can not build direct general convergence theories. At that time, it reduces to alternating projection minimization method.


## 3. Convergence Results

In fact, now the problem is a proximal alternating minimization case, whose global convergence has been detailedly analyzed in paper [12]. The paper mainly concentrates on theory analysis of problemL(x,y)=f(x)+Q(x,y)+g(y)with the following form:(10)xk+1∈arg⁡minx⁡Lx,yk+12ζkx-xk22,x∈Rm,yk+1∈arg⁡miny⁡Lxk+1,y+12μky-yk22,y∈Rn.And ifLhas the Kurdyka-Lojasiewicz (K-L) property, then each bounded sequence generated by the above algorithm converges to a critical point ofL. Even convergence rate of the algorithm can be computed, which depends on the geometrical properties of the functionLaround its critical points. It is remarkable that assumption with K-L property can be verified in many common functions.The convergence difference between algorithm of paper [12] and ours is that our minimization objective is not two variables(x,y)butxand a sequence of variablesθmm=1M. In this section, our work is to give similar consequence of algorithms (9a) and (9b).


### 3.1. Preliminary

Definition 4 (subdifferentials [24, 25]).
Letg:RN→(-∞,+∞]be a proper and lower semicontinuous function. (1)
For a givenx∈dom⁡(g), the Fréchet subdifferential ofgat x, written as∂^g(x), is the set of all vectorsu∈RNwhich satisfy (11)limy≠x⁡infy→x⁡Jy-Jx-u,y-xy-x2≥0.Whenx∉dom⁡(g), we set∂^J(x)=∅.(2)
The “limiting” subdifferential, or simply the subdifferential, ofJatx∈RN, written as∂g(x), is defined through the following closure process: (12)∂gx≔u∈RN:∃xk⟶x,gxk⟶gx,uk∈∂^gxk⟶uask⟶∞.A necessary condition forz∈X to be a minimizer ofgis(13)0∈∂gx.A point that satisfies (13) is called limiting-critical or simply critical point. The set of critical points ofgis denoted bycritg.Being givenx0,θ0∈X×∏m=1MUm, recall that sequence generated by (9a) and (9b) is of the form(xk,θk)→(xk+1,θk)→(xk+1,θk+1). According to the basic properties ofP, we can deduce a few important points.Corollary 5.
Assume sequences(xk,θk)are generated by (9a) and (9b) under assumption H, and then they are well defined. Moreover, consider the following:(i)
The following estimate holds:(14)Pxk,θk+12ζkxk-xk-122+∑m=1M12μkθmk-θmk-122⩽Pxk-1,θk-1,∀k≥1;henceLxk,θkdose not increase.(ii)
(15)∑k=1∞xk-xk-122+∑m=1Mθmk-θmk-122<∞;hencelimk→∞⁡xk-xk-12+∑m=1Mθmk-θmk-12=0.(iii)
Fork≥1, we have(16)∑m=1M∇xQmxk,θmk-∇xQmxk,θmk-1,0-1ζk-1xk-xk-1,∏m=1M1μk-1θmk-θmk-1∈∂Pxk,θk,where above0is a multivector with the same dimension of∏m=1Mθm.
Besides, for all bounded subsequencexkj,θkjof(xk,θk),dist⁡(0,∂P(xkj,θkj))→0,kj→+∞.Corollary 6.
Assume thatH hold. Letxk,θkbe a sequence complying with (9a) and (9b). Letωx0,θ0denote the set (possibly empty) of its limit points. Then (i)
if(xk,θk)is bounded, thenω(x0,θ0)is a nonempty compact connected set anddist⁡((xk,θk),ω(x0,θ0))→0ask→+∞,(ii)
ω(x0,θ0)⊂critP,(iii)
Pis finite and constant onω(x0,y0), equal toinfk∈N⁡P(xk,yk)=limk→+∞⁡P(xk,θk).The above proposition gives some convergence results about sequences generated by (7), (9a), and (9b). Point (ii) guarantees that all limiting points produced by (9a) and (9b) must be limiting-critical points. And (iii) gives the point that objectivePconverges to the finite and constant.


### 3.2. Convergence to a Critical Point

This part gives more precise convergence analysis about the proximal algorithms (9a) and (9b).Letf:X→R∪{+∞}be a proper lower semicontinuous function. For-∞<η1<η2≤+∞, let us set[η1<f<η2]=x∈X:1<fx<2. Then we give an important definition in the optimization theory.Definition 7 (K-L property [12]).
The functionfis said to have the K-L property atx¯∈dom⁡∂fif there existη∈(0,+∞], a neighborhoodUofx¯,and a continuous concave functionϕ:[0,η)→R+such that (i)
ϕ(0)=0;(ii)
ϕisC1on(0,η);(iii)
for alls∈(0,η), ϕ′(s)>0;(iv)
for allx∈U∩[f(x¯)<f<f(x¯)+η], the Kurdyka-Lojasiewicz inequality holds: (17)ϕfx-fx¯dist⁡0,∂fx≥1.If we justify that a function has K-L property, we should estimateη,U,ϕ. Many convex functions, for instance, satisfy the above property withU=Rnandη=+∞. Besides, a lot of nonconvex examples are also given in paper [12].Below, we will give convergence analysis to critical point.Theorem 8 (convergence).
Assume thatPsatisfies H and has the Kurdyka-Lojasiewicz property at each point of the domain off. Then (i)
either(xk,θk)2tends to infinity,(ii)
orxk-xk-1,θk-θk-1isl1, that is, (18)∑k=1+∞xk-xk-12+∑m=1Mθmk-θmk-12<+∞, and, as a consequence,xk,θkconverges to a critical point ofP.The above theorem’s proof is based on the same analysis process in paper [12], so here we just present the convergence results but omit their proofs.


## 3.1. Preliminary

Definition 4 (subdifferentials [24, 25]).
Letg:RN→(-∞,+∞]be a proper and lower semicontinuous function. (1)
For a givenx∈dom⁡(g), the Fréchet subdifferential ofgat x, written as∂^g(x), is the set of all vectorsu∈RNwhich satisfy (11)limy≠x⁡infy→x⁡Jy-Jx-u,y-xy-x2≥0.Whenx∉dom⁡(g), we set∂^J(x)=∅.(2)
The “limiting” subdifferential, or simply the subdifferential, ofJatx∈RN, written as∂g(x), is defined through the following closure process: (12)∂gx≔u∈RN:∃xk⟶x,gxk⟶gx,uk∈∂^gxk⟶uask⟶∞.A necessary condition forz∈X to be a minimizer ofgis(13)0∈∂gx.A point that satisfies (13) is called limiting-critical or simply critical point. The set of critical points ofgis denoted bycritg.Being givenx0,θ0∈X×∏m=1MUm, recall that sequence generated by (9a) and (9b) is of the form(xk,θk)→(xk+1,θk)→(xk+1,θk+1). According to the basic properties ofP, we can deduce a few important points.Corollary 5.
Assume sequences(xk,θk)are generated by (9a) and (9b) under assumption H, and then they are well defined. Moreover, consider the following:(i)
The following estimate holds:(14)Pxk,θk+12ζkxk-xk-122+∑m=1M12μkθmk-θmk-122⩽Pxk-1,θk-1,∀k≥1;henceLxk,θkdose not increase.(ii)
(15)∑k=1∞xk-xk-122+∑m=1Mθmk-θmk-122<∞;hencelimk→∞⁡xk-xk-12+∑m=1Mθmk-θmk-12=0.(iii)
Fork≥1, we have(16)∑m=1M∇xQmxk,θmk-∇xQmxk,θmk-1,0-1ζk-1xk-xk-1,∏m=1M1μk-1θmk-θmk-1∈∂Pxk,θk,where above0is a multivector with the same dimension of∏m=1Mθm.
Besides, for all bounded subsequencexkj,θkjof(xk,θk),dist⁡(0,∂P(xkj,θkj))→0,kj→+∞.Corollary 6.
Assume thatH hold. Letxk,θkbe a sequence complying with (9a) and (9b). Letωx0,θ0denote the set (possibly empty) of its limit points. Then (i)
if(xk,θk)is bounded, thenω(x0,θ0)is a nonempty compact connected set anddist⁡((xk,θk),ω(x0,θ0))→0ask→+∞,(ii)
ω(x0,θ0)⊂critP,(iii)
Pis finite and constant onω(x0,y0), equal toinfk∈N⁡P(xk,yk)=limk→+∞⁡P(xk,θk).The above proposition gives some convergence results about sequences generated by (7), (9a), and (9b). Point (ii) guarantees that all limiting points produced by (9a) and (9b) must be limiting-critical points. And (iii) gives the point that objectivePconverges to the finite and constant.


## 3.2. Convergence to a Critical Point

This part gives more precise convergence analysis about the proximal algorithms (9a) and (9b).Letf:X→R∪{+∞}be a proper lower semicontinuous function. For-∞<η1<η2≤+∞, let us set[η1<f<η2]=x∈X:1<fx<2. Then we give an important definition in the optimization theory.Definition 7 (K-L property [12]).
The functionfis said to have the K-L property atx¯∈dom⁡∂fif there existη∈(0,+∞], a neighborhoodUofx¯,and a continuous concave functionϕ:[0,η)→R+such that (i)
ϕ(0)=0;(ii)
ϕisC1on(0,η);(iii)
for alls∈(0,η), ϕ′(s)>0;(iv)
for allx∈U∩[f(x¯)<f<f(x¯)+η], the Kurdyka-Lojasiewicz inequality holds: (17)ϕfx-fx¯dist⁡0,∂fx≥1.If we justify that a function has K-L property, we should estimateη,U,ϕ. Many convex functions, for instance, satisfy the above property withU=Rnandη=+∞. Besides, a lot of nonconvex examples are also given in paper [12].Below, we will give convergence analysis to critical point.Theorem 8 (convergence).
Assume thatPsatisfies H and has the Kurdyka-Lojasiewicz property at each point of the domain off. Then (i)
either(xk,θk)2tends to infinity,(ii)
orxk-xk-1,θk-θk-1isl1, that is, (18)∑k=1+∞xk-xk-12+∑m=1Mθmk-θmk-12<+∞, and, as a consequence,xk,θkconverges to a critical point ofP.The above theorem’s proof is based on the same analysis process in paper [12], so here we just present the convergence results but omit their proofs.


## 4. Application to 1D TV Denoising

In practical scientific and engineering contexts, noise removal is the basis and requisite of other subsequential applications. It has received extensive attentions. A range of computational algorithms have been proposed to solve the denoising problem [26–28]. Among these solvers, total variation (TV) regularizer is of great importance since it can efficiently deal with noisy signals that have sparse derivatives (or gradients). For instance, piecewise constant (PWC) [29] signal with noise, whose derivative is sparse relative to signal dimension, could be denoised by powerful TV denoising method.In 1D TV denoising problem [10], one needs to solve model (2). TV denoising minimizes a composite of two parts. The first part is to keep the error, between the observed data and the original, as small as possible. The second is devoted to minimizing the sparsity of the gradients. Usually, denoising model is defined as one combination of a quadratic data fidelity term and a convex regularization term or a differential regularization, for example, convex but nonsmooth problem [30](19)arg⁡minx∈Rn⁡y-x22+λDx1:x∈Rnor differential but nonconvex problem [7](20)arg⁡minx∈Rn⁡y-x22+λ∑i1alog⁡1+aDxi:x∈Rn,a>0,wherex1=∑i=1nxi. Exact solution of the above two types can be obtained by very fast direct algorithms [7, 30]. In fact, convexl1-norm is the replacement of nonconvexl0-norm in (19) since convex optimization techniques have been deeply studied. The latter, like logarithmic penalty and arctangent penalty, can be solved by MM update iteration, in which total objective function (including data fidelity and regularization terms) should meet strictly convex condition.In this test, we apply our algorithms (9a) and (9b) to this example. Auxiliary variableθis introduced to reduce complexity of the compositionDx0. Then (2) is represented as(21)arg⁡minx∈Rn,θ∈Rn-1⁡x-y22+τ2Dx-θ22+λθ0.Apparently this problem satisfies the convergence conditions [12]. Concrete steps by algorithms (9a) and (9b) are shown in (22a)xk+1∈arg⁡minu∈Rn⁡u-y22+τ2Du-θk22+12ζku-xk22;(22b)θk+1∈arg⁡minv∈Rn-1⁡λv0+τ2Dxk+1-v22+12μkv-θk22.In fact, when testsζkandμkcan be set as very large constants, the last proximal terms1/2ζku-xk22and 1/2μkv-θk22in each computing step could be negligible. Hence, the computation of (22a) and (22b) is as follows.Computation of (22a). The former step (22a) is a quadratic function and could be computed through many techniques, like gradient descent.Computation of (22b). Apparently, the latter (22b) could be rewritten as a proximal operator of functiongθ [12]; that is,proxλ/τθ0⁡Dxk+1. Consider the proximal operatorproxλ/τ·0⁡(u). Whenn=1, l0norm is reduced to·0, where one easily establishes that (23)proxλ/τ·0⁡u=u,ifu>2λτ;0,u,ifu=2λτ;0,otherwise.Whennis arbitrary, trivial algebraic manipulations are given, withu=(u1,u2,…,un)∈Rn: (24)proxλ/τ·0⁡ui=proxλ/τ·0⁡ui, and thusproxλ/τθ0⁡Dxk+1is a perfectly known object.Total variation denoising examples with three convex and nonconvex regularization instances (the two others are convex and nonconvex but smooth algorithms in [7, 30], resp.) are figured in Figure 1. Original piece signal datax0∈Rnwith lengthn=256is obtained with MakeSignal in paper [7]. The noisy data is obtained using additive white Gaussian noise (AWGN) (σ=0.5). For both convex and nonconvex cases, we setλ=nσ/4, consistent with the range suggested in [30] for standard (convex) TV denoising and nonconvexity parameter is set to its maximal value,a=1/(4λ)default in [7]. These settings could lead to the best denoising result in their papers. All the other settings are consistent with paper [7]. The maximum iteration numbers are all500. All the codes are tested in the same computer.Figure 1
Total variation denoising with nonconvex nonsmoothl0penalty (RMSE = 0.1709), compared with convexl1penalty [21] (RMSE = 0.2720) and smooth log penalty [7] (RMSE = 0.2611).
(a) (b) (c) (d)According to the comparison between our algorithm for TV-l0norm and the proposed algorithms in papers [7, 30], our algorithm has better result with smaller Root-Mean-Square-Error (RMSE), whereRMSE(x)≔(1/n)∑i=1nxi-xi022. Referring to Figure 1, the best RMSE results for 1D TV denoising with convexl1penalty [21] and smoothlogpenalty [7] are 0.2720 and 0.2611, respectively. And ours is 0.1709, much better than the convex and smooth cases.


## 5. Conclusion

Nonconvex nonsmooth algorithm finds many interesting applications in many fields. In this paper, we give a general proximal alternating minimization method for a kind of nonconvex nonsmooth problems with complex composition. It has concise form, good theory results, and promising numerical result. For specific 1D standard TV denoising problem, the improvement is more dramatic compared to the existing algorithms [7, 30]. Besides, our algorithm works on other nonconvex nonsmooth problems, such as block sparse recovery, group lasso, and image deconvolution, of which the examples are just too numerous to mention.



---

*Source: 5053434-2016-11-03.xml*