Table Of ContentIEEESIGNALPROCESSINGLETTERS(PREPRINT) 1
Total Variation Denoising via the Moreau Envelope
Ivan Selesnick
Abstract—Total variation denoising is a nonlinear filtering A. Relation to Prior Work
method well suited for the estimation of piecewise-constant
signals observed in additive white Gaussian noise. The method Numerous non-convex penalties and algorithms have been
is defined by the minimization of a particular non-differentiable proposed to outperform (cid:96) -norm regularization for the esti-
1
convex cost function. This paper describes a generalization of mation of sparse signals e.g., [8], [10], [12], [13], [15], [33],
this cost function that can yield more accurate estimation of
[34], [39], [43], [51], [53]. However, few of these methods
piecewiseconstantsignals.Thenewcostfunctioninvolvesanon-
7 maintain the convexity of the cost function. The prescription
convex penalty (regularizer) designed to maintain the convexity
1 of the cost function. The new penalty is based on the Moreau of non-convex penalties maintaining cost function convexity
0
envelope. The proposed total variation denoising method can be was pioneered by Blake, Zisserman, and Nikolova [6], [36],
2
implemented using forward-backward splitting. [39], [40], and further developed in Refs. [3], [4], [14], [20],
n
[28], [30], [32], [42], [46], [48]. These works rely on the
a
I. INTRODUCTION presence of both strongly and weakly convex terms, which
J
Total variation (TV) denoising is a nonlinear filtering is also exploited in [35].
2
method based on the assumption that the underlying signal is Theproposedpenaltyisexpressedasadifferentiableconvex
] piecewise constant (equivalently, the derivative of the under- function subtracted from the standard penalty (i.e., (cid:96) norm).
C 1
lying signal is sparse) [45]. Such signals arise in geoscience, Previous works also use this idea [41], [42], [47]. But the dif-
O
biophysics, and other areas [31]. The TV denoising technique ferentiable convex functions used therein are either separable
h. is also used in conjunction with other methods in order to [41], [42] or sums of bivariate functions [47].
t process more general types of signals [20], [23], [24], [26]. In parallel with the submission of this paper, Carlsson has
a
m Total variation denoising is prototypical of methods based alsoproposedusingMoreauenvelopestoprescribenon-trivial
on sparse signal models. It is defined by the minimization convexcostfunctions[9].Whiletheapproachin[9]startswith
[
of a convex cost function comprising a quadratic data fidelity a given non-convex cost function (e.g., with the (cid:96) pseudo-
0
1 termandanon-differentiableconvexpenaltyterm.Thepenalty norm penalty) and seeks the convex envelope, our approach
9v term is the composition of a linear operator and the (cid:96)1 norm. startswiththe(cid:96)1-normpenaltyandseeksaclassofconvexity-
3 Although the (cid:96)1 norm stands out as the convex penalty that preserving penalties.
4 most effectively induces sparsity [27], non-convex penalties Some forms of generalized TV are based on infimal con-
0 can lead to more accurate estimation of the underlying signal volution (related to the Moreau envelope) [5], [7], [11],
0 [37], [38], [40], [44], [50]. [49]. But these works propose convex penalties suitable for
1. Afewrecentpapersconsidertheprescriptionofnon-convex non-piecewise-constant signals, while we propose non-convex
0 penalties that maintain the convexity of the TV denoising penalties suitable for piecewise-constant signals.
7 cost function [1], [30], [32], [48]. (The motivation for this is
1 to leverage the benefits of both non-convex penalization and
:
v convexoptimization,e.g.,toaccuratelyestimatetheamplitude II. TOTALVARIATIONDENOISING
i of jump discontinuities while guaranteeing the uniqueness of
X Definition 1. Given y ∈ RN and λ > 0, total variation
thesolution.)Thepenaltiesconsideredintheseworksaresep-
r denoising is defined as
a arable (additive). But non-separable penalties can outperform
separable penalties in this context. This is because preserving
tvd(y;λ)=arg min(cid:8)1(cid:107)y−x(cid:107)2+λ(cid:107)Dx(cid:107) (cid:9) (1)
the convexity of the cost function is a severely limiting x∈RN 2 2 1
requirement. Non-separable penalties can more successfully
=prox (y) (2)
λ(cid:107)D·(cid:107)
meet this requirement because they are more general than 1
separable penalties [47]. where D is the (N −1)×N matrix
This paper proposes a non-separable non-convex penalty
for total variation denoising that generalizes the standard −1 1
penalty and maintains the convexity of the cost function to be −1 1
menivneilmopizee,dc.a1nTmheorneeawccpuernataelltyy,ewsthimicahteisthbeaasemdplointudtheesoMfojuremapu D = ... ... . (3)
−1 1
discontinuities in an underlying piecewise constant signal.
TheauthoriswiththeTandonSchoolofEngineering,NewYorkUniversity, As indicated in (2), TV denoising is the proximity operator
NewYork,USA.Email:[email protected]. [16] of the function x (cid:55)→ λ(cid:107)Dx(cid:107) . It is convenient that TV
This work was supported by NSF under grant 1525398 and ONR under 1
denoising can be calculated exactly in finite-time [17], [18],
grantN00014-15-1-2314.
1Softwareisavailableathttp://eeweb.poly.edu/iselesni/mtvd [22], [29].
2 IEEESIGNALPROCESSINGLETTERS(PREPRINT)
III. MOREAUENVELOPE TheproposedpenaltyisupperboundedbythestandardTV
penalty, which is recovered as a special case.
Before we define the non-differentiable non-convex penalty
inSec.IV,wefirstdefineadifferentiableconvexfunction.We Proposition 5. Let α(cid:62)0. The penalty ψ satisfies
α
use the Moreau envelope from convex analysis [2].
ψ (x)=(cid:107)Dx(cid:107) , ∀x∈RN (12)
0 1
Definition 2. Let α(cid:62)0. We define S : RN →R as
α and
S (x)= min(cid:8)(cid:107)Dv(cid:107) + α(cid:107)x−v(cid:107)2(cid:9) (4) 0(cid:54)ψ (x)(cid:54)(cid:107)Dx(cid:107) , ∀x∈RN. (13)
α v∈RN 1 2 2 α 1
Proof: It follows from (5) and (7).
where D is the first-order difference matrix (3).
When a convex function is subtracted from another convex
If α>0, then S is the Moreau envelope of index α−1 of function [as in (11)], the resulting function may well be
α
the function x(cid:55)→(cid:107)Dx(cid:107) . negative on part of its domain. Inequality (13) states that the
1
proposed penalty ψ avoids this fate. This is relevant because
Proposition 1. The function S can be calculated by α
α the penalty function should be non-negative.
S (x)=0 (5) Figures in the supplemental material show examples of the
0
proposed penalty ψ and the function S .
S (x)=(cid:107)Dtvd(x;1/α)(cid:107) α α
α 1
+ α(cid:107)x−tvd(x;1/α)(cid:107)2, α>0. (6) V. ENHANCEDTVDENOISING
2 2
Wedefine‘Moreau-enhanced’TVdenoising.Ifα>0,then
Proof: For α =0: Setting α =0 and v =0 in (4) gives
(5).Forα>0:BythedefinitionofTVdenoising,thev ∈RN the proposed penalty penalizes large amplitude values of Dx
less than the (cid:96) norm does (i.e., ψ (x)(cid:54)(cid:107)Dx(cid:107) ), hence it is
minimizing the function in (4) is the TV denoising of x, i.e., 1 α 1
less likely to underestimate jump discontinuities.
vopt =tvd(x,1/α).
Definition 4. Given y ∈ RN, λ > 0, and α (cid:62) 0, we define
Proposition 2. Let α(cid:62)0. The function S satisfies
α Moreau-enhanced total variation denoising as
0(cid:54)Sα(x)(cid:54)(cid:107)Dx(cid:107)1, ∀x∈RN. (7) mtvd(y;λ,α)=arg min(cid:8)1(cid:107)y−x(cid:107)2+λψ (x)(cid:9) (14)
Proof: From (4), we have S (x)(cid:54)(cid:107)Dv(cid:107) +(α/2)(cid:107)x− x∈RN 2 2 α
α 1 where ψ is given by (11).
v(cid:107)2 for all v ∈ RN. In particular, v = x leads to S (x) (cid:54) α
2 α
(cid:107)Dx(cid:107) . Also, S (x) (cid:62) 0 since S (x) is defined as the The parameter α controls the non-convexity of the penalty.
1 α α
minimum of a non-negative function. Ifα=0,thenthepenaltyisconvexandMoreau-enhancedTV
denoising reduces to TV denoising. Greater values of α make
Proposition 3. Let α (cid:62) 0. The function S is convex and
α the penalty more non-convex. What is the greatest value of
differentiable.
α that maintains convexity of the cost function? The critical
Proof: It follows from Proposition 12.15 in Ref. [2]. value is given by Theorem 1.
Proposition 4. Let α(cid:62)0. The gradient of Sα is given by Theorem 1. Let λ>0 and α(cid:62)0. Define Fα: RN →R as
∇S0(x)=0 (8) Fα(x)= 21(cid:107)y−x(cid:107)22+λψα(x) (15)
(cid:0) (cid:1)
∇S (x)=α x−tvd(x;1/α) , α>0 (9) where ψ is given by (11). If
α α
where tvd denotes total variation denoising (1). 0(cid:54)α(cid:54)1/λ (16)
Proof: Since S is the Moreau envelope of index α−1 thenFα isconvex.If0(cid:54)α<1/λthenFα isstronglyconvex.
α
of the function x (cid:55)→ (cid:107)Dx(cid:107)1 when α > 0, it follows by Proof: We write the cost function as
Proposition 12.29 in Ref. [2] that
F (x)= 1(cid:107)y−x(cid:107)2+λ(cid:107)Dx(cid:107) −λS (x) (17)
(cid:0) (cid:1) α 2 2 1 α
∇Sα(x)=α x−prox(1/α)(cid:107)D·(cid:107) (x) . (10) = 1(cid:107)y−x(cid:107)2+λ(cid:107)Dx(cid:107)
1 2 2 1
This proximity operator is TV denoising, giving (9). −λ min(cid:8)(cid:107)Dv(cid:107) + α(cid:107)x−v(cid:107)2(cid:9) (18)
v∈RN 1 2 2
= max(cid:8)1(cid:107)y−x(cid:107)2+λ(cid:107)Dx(cid:107)
IV. NON-CONVEXPENALTY v∈RN 2 2 1
To strongly induce sparsity of Dx, we define a non-convex −λ(cid:107)Dv(cid:107) − λα(cid:107)x−v(cid:107)2(cid:9) (19)
1 2 2
generalization of the standard TV penalty. The new penalty is = max(cid:8)1(1−λα)(cid:107)x(cid:107)2+λ(cid:107)Dx(cid:107) +g(x,v)(cid:9) (20)
defined by subtracting a differentiable convex function from v∈RN 2 2 1
the standard penalty. = 1(1−λα)(cid:107)x(cid:107)2+λ(cid:107)Dx(cid:107) + maxg(x,v) (21)
Definition 3. Let α(cid:62)0. We define the penalty ψ : RN →R 2 2 1 v∈RN
α
where g(x,v) is affine in x. The last term is convex as it is
as
the point-wise maximum of a set of convex functions. Hence,
ψ (x)=(cid:107)Dx(cid:107) −S (x) (11)
α 1 α F is a convex function if 1−λα (cid:62) 0. If 1−λα > 0, then
α
where D is the matrix (3) and S is defined by (4). F is strongly convex (and strictly convex).
α α
3
VI. ALGORITHM In summary, the proposed Moreau-enhanced TV denoising
Proposition 6. Let y ∈ RN, λ > 0, and 0 < α < 1/λ. Then method comprises the steps:
x(k) produced by the iteration 1) Set the regularization parameter λ (λ>0).
z(k) =y+λα(cid:0)x(k)−tvd(x(k);1/α)(cid:1) (22a) 2) Set the non-convexity parameter α (0(cid:54)α<1/λ).
3) Initialize x(0) =0.
x(k+1) =tvd(z(k);λ). (22b) 4) Run iteration (22) until convergence.
converges to the solution of the Moreau-enhanced TV denois-
ing problem (14).
VII. OPTIMALITYCONDITION
Proof: If the cost function (15) is strongly convex, then
To avoid terminating the iterative algorithm too early, it is
the minimizer can be calculated using the forward-backward
useful to verify convergence using an optimality condition.
splitting (FBS) algorithm [2], [16]. This algorithm minimizes
a function of the form Proposition 7. Let y ∈RN, λ>0, and 0<α<1/λ. If x is
a solution to (14), then
F(x)=f (x)+f (x) (23)
1 2
(cid:2) (cid:0) (cid:1)(cid:3)
where both f1 and f2 are convex and ∇f1 is additionally C (x−y)/λ+α(tvd(x;1/α)−x) n ∈sign([Dx]n) (32)
Lipschitz continuous. The FBS algorithm is given by
for n=0,...,N −1, where C ∈R(N−1)×N is given by
z(k) =x(k)−µ(cid:2)∇f (x(k))(cid:3) (24a)
1 (cid:40)
x(k+1) =argmxin(cid:8)12(cid:107)z(k)−x(cid:107)22+µf2(x)(cid:9) (24b) Cm,n = 10,, mm(cid:62)<nn, i.e., [Cx]n = (cid:88) xm (33)
m(cid:54)n
where 0 < µ < 2/ρ and ρ is the Lipschitz constant of ∇f .
1
The iterates x(k) converge to a minimizer of F. and sign is the set-valued signum function
To apply the FBS algorithm to the proposed cost function
(15), we write it as {−1}, t<0
F (x)= 1(cid:107)y−x(cid:107)2+λψ (x), (25) sign(t)= [−1,1], t=0 (34)
α 2 2 α {1}, t>0.
= 1(cid:107)y−x(cid:107)2+λ(cid:107)Dx(cid:107) −λS (x) (26)
2 2 1 α
=f1(x)+f2(x) (27) According to (32), if x ∈ RN is a minimizer, then the
points ([Dx] ,u )∈R2 must lie on the graph of the signum
where n n
function, where u denotes the value on the left-hand side
n
f1(x)= 12(cid:107)y−x(cid:107)22−λSα(x) (28a) of (32). Hence, the optimality condition can be depicted as a
f (x)=λ(cid:107)Dx(cid:107) . (28b) scatter plot. Figures in the supplemental material show how
2 1
the points in the scatter plot converge to the signum function
The gradient of f is given by
1 as the algorithm (22) progresses.
∇f (x)=x−y−λ∇S (x) (29) Proof of Proposition 7: A vector x minimizes a convex
1 α
(cid:0) (cid:1) function F if 0 ∈ ∂F(x) where ∂F(x) is the subdifferential
=x−y−λα x−tvd(x;1/α) (30)
ofF atx.Thesubdifferentialofthecostfunction(15)isgiven
usingProposition4.SubtractingS fromf doesnotincrease by
α 1
theLipschitzconstantof∇f ,thevalueofwhichis1.Hence,
1
we may set 0<µ<2. ∂Fα(x)=x−y−λ∇Sα(x)+∂(λ(cid:107)D·(cid:107)1)(x) (35)
Using (28), the FBS algorithm (24) becomes
which can be written as
z(k) =x(k)−µ(cid:2)x(k)−y
−λα(cid:0)x(k)−tvd(x(k);1/α)(cid:1)(cid:3) (31a) ∂F (x)={x−y−λ∇S (x)+λDTu
α α
x(k+1) =argmxin(cid:8)12(cid:107)z(k)−x(cid:107)22+µλ(cid:107)Dx(cid:107)1(cid:9). (31b) :un ∈sign([Dx]n), u∈RN−1}. (36)
Note that (31b) is TV denoising (1). Using the value µ = Hence, the condition 0∈∂F (x) can be written as
α
1 gives iteration (22). (Experimentally, we found this value
yields fast convergence.) (y−x)/λ+∇S (x)
α
Each iteration of (22) entails solving two standard TV ∈{DTu:u ∈sign([Dx] ), u∈RN−1}. (37)
n n
denoising problems. In this work, we calculate TV denoising
usingthefastexactClanguageprogrambyCondat[17].Like LetC beamatrixofsize(N−1)×N suchthatCDT =−I,
the iterative shrinkage/thresholding algorithm (ISTA) [19],
e.g.,(33).Itfollowsthatthecondition0∈∂F (x)impliesthat
α
[25], algorithm (22) can be accelerated in various ways.
We suggest not setting α too close to the critical value 1/λ (cid:2)C(cid:0)(x−y)/λ−∇S (x)(cid:1)(cid:3) ∈sign([Dx] ) (38)
α n n
because the FBS algorithm generally converges faster when
the cost function is more strongly convex (α<1). for n=0,...,N −1. Using Proposition 4 gives (32).
4 IEEESIGNALPROCESSINGLETTERS(PREPRINT)
(a) Noisy signal (σ = 0.50) TV denoising
6
L1 norm
4 0.5 Scalar log
Scalar MC
E
S 0.4 Proposed
2 M
R
e 0.3
0 g
a
er
v 0.2
−2 A
0.1
0 50 100 150 200 250
0
(b) Standard convex penalty (L1 norm) 0.2 0.4 0.6 0.8 1
6 Noise standard deviation (σ)
RMSE = 0.318, MAE = 0.225
4 Fig. 2. TV denoising using four penalties: RMSE as a function of noise
level.
2
0 Figure 1 shows the result of TV denoising with three
different penalties. In each case, a convex cost function is
−2 minimized.Figure1(b)showstheresultusingstandardTVde-
noising (i.e., using the (cid:96) -norm). This denoised signal consis-
0 50 100 150 200 250 1
tently underestimates the amplitudes of jump discontinuities,
(c) Scalar MC penalty especially those occurring near other jump discontinuities of
6 opposite sign. Figure 1(c) shows the result using a separable
RMSE = 0.240, MAE = 0.172
non-convex penalty [48]. This method can use any non-
4
convex scalar penalty satisfying aprescribed set of properties.
Here we use the minimax-concave (MC) penalty [3], [52]
2
with non-convexity parameter set to maintain cost function
0 convexity. This result significantly improves the root-mean-
squareerror(RMSE)andmean-absolute-deviation(MAE),but
−2 still underestimates the amplitudes of jump discontinuities.
Moreau-enhancedTVdenoising,showninFig.1(d),further
0 50 100 150 200 250
reduces the RMSE and MAE and more accurately estimates
the amplitudes of jump discontinuities. The proposed non-
(d) Proposed penalty
6 separable non-convex penalty avoids the consistent underes-
RMSE = 0.200, MAE = 0.130 timation of discontinuities seen in Figs. 1(b) and 1(c).
4 To further compare the denoising capability of the consid-
ered penalties, we calculate the average RMSE as a function
2
of the noise level. We let the noise standard deviation span
the interval 0.2(cid:54)σ (cid:54)1.0. For each σ value, we calculate the
0
average RMSE of 100 noise realizations. Figure 2 shows that
−2 the proposed penalty yields the lowest average RMSE for all
σ (cid:62) 0.4. However, at low noise levels, separable convexity-
0 50 100 150 200 250
preserving penalties [48] perform better than the proposed
non-separable convexity-preserving penalty.
Fig.1. Totalvariationdenoisingusingthreedifferentpenalties.(Thedashed
lineisthetruenoise-freesignal.) IX. CONCLUSION
This paper demonstrates the use of the Moreau envelope to
define a non-separable non-convex TV denoising penalty that
VIII. EXAMPLE
maintains the convexity of the TV denoising cost function.
This example applies TV denoising to the noisy piecewise ThebasicideaistosubtractfromaconvexpenaltyitsMoreau
constant signal shown in Fig. 1(a). This is the ‘blocks’ signal envelope. This idea should also be useful for other problems,
(length N = 256) generated by the Wavelab [21] function e.g., analysis tight-frame denoising [41].
MakeSignal with additive white Gaussian noise (σ =0.5). Separableconvexity-preservingpenalties[48]outperformed
√
We settheregularization parameterto λ= Nσ/4 following theproposedoneatlownoiselevelsintheexample.Itisyetto
adiscussioninRef.[22].ForMoreau-enhancedTVdenoising, be determined if a more general class of convexity-preserving
we set the non-convexity parameter to α=0.7/λ. penalties can outperform both across all noise levels.
5
REFERENCES [28] W.He,Y.Ding,Y.Zi,andI.W.Selesnick.Sparsity-basedalgorithmfor
detecting faults in rotating machines. Mechanical Systems and Signal
Processing,72-73:46–64,May2016.
[1] F. Astrom and C. Schnorr. On coupled regularization for non-convex
variational image enhancement. In IAPR Asian Conf. on Pattern [29] N.A.Johnson. Adynamicprogrammingalgorithmforthefusedlasso
Recognition(ACPR),pages786–790,November2015. andL0-segmentation. J.Computat.Graph.Stat.,22(2):246–260,2013.
[30] A.Lanza,S.Morigi,andF.Sgallari. Conveximagedenoisingvianon-
[2] H. H. Bauschke and P. L. Combettes. Convex Analysis and Monotone
convex regularization with parameter selection. J. Math. Imaging and
OperatorTheoryinHilbertSpaces. Springer,2011.
[3] ˙I.Bayram. Penaltyfunctionsderivedfrommonotonemappings. IEEE Vision,pages1–26,2016.
[31] M.A.LittleandN.S.Jones.Generalizedmethodsandsolversfornoise
SignalProcessingLetters,22(3):265–269,March2015.
removal from piecewise constant signals: Part I – background theory.
[4] I. Bayram. On the convergence of the iterative shrinkage/thresholding
Proc.R.Soc.A,467:3088–3114,2011.
algorithmwithaweaklyconvexpenalty. IEEETrans.SignalProcess.,
[32] M. Malek-Mohammadi, C. R. Rojas, and B. Wahlberg. A class
64(6):1597–1608,March2016.
of nonconvex penalties preserving overall convexity in optimization-
[5] S.BeckerandP.L.Combettes. Analgorithmforsplittingparallelsums
basedmeanfiltering. IEEETrans.SignalProcess.,64(24):6650–6664,
of linearly composed monotone operators, with applications to signal
December2016.
recovery. J.NonlinearandConvexAnalysis,15(1):137–159,2014.
[33] Y. Marnissi, A. Benazza-Benyahia, E. Chouzenoux, and J.-C. Pesquet.
[6] A.BlakeandA.Zisserman. VisualReconstruction. MITPress,1987.
Generalizedmultivariateexponentialpowerpriorforwavelet-basedmul-
[7] M.Burger,K.Papafitsoros,E.Papoutsellis,andC.-B.Scho¨nlieb.Infimal tichannelimagerestoration. InProc.IEEEInt.Conf.ImageProcessing
convolution regularisation functionals of BV and Lp spaces. J. Math. (ICIP),pages2402–2406,September2013.
ImagingandVision,55(3):343–369,2016. [34] H. Mohimani, M. Babaie-Zadeh, and C. Jutten. A fast approach for
[8] E. J. Cande`s, M. B. Wakin, and S. Boyd. Enhancing sparsity by overcompletesparsedecompositionbasedonsmoothedl0norm. IEEE
reweighted l1 minimization. J. Fourier Anal. Appl., 14(5):877–905, Trans.SignalProcess.,57(1):289–301,January2009.
December2008. [35] T. Mo¨llenhoff, E. Strekalovskiy, M. Moeller, and D. Cremers. The
[9] M. Carlsson. On convexification/optimization of functionals including primal-dual hybrid gradient method for semiconvex splittings. SIAM
anl2-misfitterm. https://arxiv.org/abs/1609.09378,September2016. J.Imag.Sci.,8(2):827–857,2015.
[10] M.CastellaandJ.-C.Pesquet. OptimizationofaGeman-McClurelike [36] M. Nikolova. Estimation of binary images by minimizing convex
criterionforsparsesignaldeconvolution.InIEEEInt.WorkshopComput. criteria. InProc.IEEEInt.Conf.ImageProcessing(ICIP),pages108–
Adv.Multi-SensorAdaptiveProc.,pages309–312,December2015. 112vol.2,1998.
[11] A. Chambolle and P.-L. Lions. Image recovery via total variation [37] M.Nikolova.Localstronghomogeneityofaregularizedestimator.SIAM
minimization and related problems. Numerische Mathematik, 76:167– J.Appl.Math.,61(2):633–658,2000.
188,1997. [38] M.Nikolova. Analysisoftherecoveryofedgesinimagesandsignals
[12] R.Chartrand. Shrinkagemappingsandtheirinducedpenaltyfunctions. byminimizingnonconvexregularizedleast-squares. MultiscaleModel.
InProc.IEEEInt.Conf.Acoust.,Speech,SignalProcessing(ICASSP), Simul.,4(3):960–991,2005.
pages1026–1029,May2014. [39] M. Nikolova. Energy minimization methods. In O. Scherzer, editor,
[13] L. Chen and Y. Gu. The convergence guarantees of a non-convex HandbookofMathematicalMethodsinImaging,chapter5,pages138–
approachforsparserecovery.IEEETrans.SignalProcess.,62(15):3754– 186.Springer,2011.
3767,August2014. [40] M. Nikolova, M. K. Ng, and C.-P. Tam. Fast nonconvex nonsmooth
[14] P.-Y. Chen and I. W. Selesnick. Group-sparse signal denoising: Non- minimization methods for image restoration and reconstruction. IEEE
convex regularization, convex optimization. IEEE Trans. Signal Pro- Trans.ImageProcess.,19(12):3073–3088,December2010.
cess.,62(13):3464–3478,July2014. [41] A. Parekh and I. W. Selesnick. Convex denoising using non-convex
[15] E. Chouzenoux, A. Jezierska, J. Pesquet, and H. Talbot. A majorize- tightframeregularization.IEEESignalProcessingLetters,22(10):1786–
minimizesubspaceapproachfor(cid:96)2−(cid:96)0 imageregularization. SIAMJ. 1790,October2015.
Imag.Sci.,6(1):563–591,2013. [42] A.ParekhandI.W.Selesnick. Enhancedlow-rankmatrixapproxima-
[16] P.L.CombettesandJ.-C.Pesquet. Proximalsplittingmethodsinsignal tion. IEEESignalProcessingLetters,23(4):493–497,April2016.
processing.InH.H.Bauschkeetal.,editors,Fixed-PointAlgorithmsfor [43] J. Portilla and L. Mancera. L0-based sparse approximation: two
InverseProblemsinScienceandEngineering,pages185–212.Springer- alternative methods and some applications. In Proceedings of SPIE,
Verlag,2011. volume6701(WaveletsXII),SanDiego,CA,USA,2007.
[44] P. Rodriguez and B. Wohlberg. Efficient minimization method for a
[17] L.Condat. Adirectalgorithmfor1-Dtotalvariationdenoising. IEEE
generalized total variation functional. IEEE Trans. Image Process.,
SignalProcessingLetters,20(11):1054–1057,November2013.
18(2):322–332,February2009.
[18] J. Darbon and M. Sigelle. Image restoration with discrete constrained
[45] L.Rudin,S.Osher,andE.Fatemi.Nonlineartotalvariationbasednoise
totalvariationPartI:Fastandexactoptimization.J.Math.Imagingand
removalalgorithms. PhysicaD,60:259–268,1992.
Vision,26(3):261–276,2006.
[46] I.W.SelesnickandI.Bayram. Sparsesignalestimationbymaximally
[19] I.Daubechies,M.Defrise,andC.DeMol. Aniterativethresholdingal-
sparseconvexoptimization. IEEETrans.SignalProcess.,62(5):1078–
gorithmforlinearinverseproblemswithasparsityconstraint.Commun.
1092,March2014.
PureAppl.Math,57(11):1413–1457,2004.
[47] I. W. Selesnick and I. Bayram. Enhanced sparsity by non-separable
[20] Y. Ding and I. W. Selesnick. Artifact-free wavelet denoising: Non-
regularization. IEEE Trans. Signal Process., 64(9):2298–2313, May
convex sparse regularization, convex optimization. IEEE Signal Pro-
2016.
cessingLetters,22(9):1364–1368,September2015.
[48] I.W.Selesnick,A.Parekh,andI.Bayram. Convex1-Dtotalvariation
[21] D. Donoho, A. Maleki, and M. Shahram. Wavelab 850, 2005.
denoising with non-convex regularization. IEEE Signal Processing
http://www-stat.stanford.edu/%7Ewavelab/.
Letters,22(2):141–144,February2015.
[22] L. Du¨mbgen and A. Kovac. Extensions of smoothing via taut strings. [49] S.Setzer,G.Steidl,andT.Teuber. Infimalconvolutionregularizations
Electron.J.Statist.,3:41–75,2009. with discrete l1-type functionals. Commun. Math. Sci., 9(3):797–827,
[23] S.DurandandJ.Froment. Reconstructionofwaveletcoefficientsusing 2011.
total variation minimization. SIAM J. Sci. Comput., 24(5):1754–1767, [50] M. Storath, A. Weinmann, and L. Demaret. Jump-sparse and sparse
2003. recovery using Potts functionals. IEEE Trans. Signal Process.,
[24] G.R.Easley,D.Labate,andF.Colonna. Shearlet-basedtotalvariation 62(14):3654–3666,July2014.
diffusion for denoising. IEEE Trans. Image Process., 18(2):260–268, [51] D. P. Wipf, B. D. Rao, and S. Nagarajan. Latent variable Bayesian
February2009. modelsforpromotingsparsity.IEEETrans.Inform.Theory,57(9):6236–
[25] M. Figueiredo and R. Nowak. An EM algorithm for wavelet-based 6255,September2011.
imagerestoration. IEEETrans.ImageProcess.,12(8):906–916,August [52] C.-H.Zhang.Nearlyunbiasedvariableselectionunderminimaxconcave
2003. penalty. TheAnnalsofStatistics,pages894–942,2010.
[26] A.GholamiandS.M.Hosseini. AbalancedcombinationofTikhonov [53] H. Zou and R. Li. One-step sparse estimates in nonconcave penalized
andtotalvariationregularizationsforreconstructionofpiecewise-smooth likelihoodmodels. Ann.Statist.,36(4):1509–1533,2008.
signals. SignalProcessing,93(7):1945–1960,2013.
[27] T. Hastie, R. Tibshirani, and M. Wainwright. Statistical learning with
sparsity:thelassoandgeneralizations. CRCPress,2015.
6 IEEESIGNALPROCESSINGLETTERS(PREPRINT)
X. SUPPLEMENTALFIGURES
To gain intuition about the proposed penalty function and
how it induces sparsity of Dx while maintaining convexity of
the cost function, a few illustrations are useful.
Figure 3 illustrates the proposed penalty, its sparsity-
inducing behavior, and its relationship to the differentiable
||Dx||
convex function S . Figure 4 illustrates how the proposed 1
α
4
penalty is able to maintain the convexity of the cost function.
Figure3showstheproposedpenalty ψ definedin(11)for
α
α=1andN =2.Itcanbeseenthatthepenaltyapproximates 3
the standard TV penalty (cid:107)D·(cid:107) for signals x for which Dx
1
is approximately zero. But it increases more slowly than the 2
standardTVpenaltyas(cid:107)Dx(cid:107)→∞.Inthatsense,itpenalizes
large values less than the standard TV penalty. 1
As shown in Fig. 3, the proposed penalty is expressed as
the standard TV penalty minus the differentiable convex non- 0
2
negative function S . Since S is flat around the null space
α α
of D, the penalty ψ approximates the standard TV penalty 0
around the null spaαce of D. In addition, since Sα is non- x2 −2 −2 −1 x 0 1 2
negative, the penalty ψ lies below the standard TV penalty. 1
α
Figure 4 shows the differentiable part of the cost function
S (x)
Fα for α = 1, λ = 1, and N = 2. The differentiable part is α
given by f in (28a). The total cost function is obtained by 4
1
adding the standard TV penalty to f , see (23). Hence, F is
1 α
convex if the differentiable part f1 is convex. As can be seen 3
in Fig. 4, the function f is convex. We note that the function
1
f1 in this figure is not strongly convex. This is because we 2
have used α=1/λ. If α<1/λ, then the function f will be
1
strongly convex (and hence F will also be strongly convex
α 1
and have a unique minimizer). We recommend α<1/λ.
Figure5showsthedifferentiablepartf ofthecostfunction
1 0
F for α = 2, λ = 1, and N = 2. Here, the function f is 2
α 1
non-convex because α > 1/λ which violates the convexity 0
condition. x −2 −2 −1 0 1 2
In order to simplify the illustration, we have set y = 0 in 2 x
1
Fig. 4. In practice y (cid:54)=0. But the only difference between the
cases y = 0 and y (cid:54)= 0 is an additive affine function which ψ (x) = ||Dx|| − S (x)
α 1 α
does not alter the convexity properties of the function.
1.5
InpracticeweareinterestedinthecaseN (cid:29)2,i.e.,signals
much longer than two samples. However, in order to illustrate
the functions, we are limited to the case of N = 2. We note 1
thatthecaseofN =2doesnotfullyillustratethebehaviorof
the proposed penalty. In particular, when N = 2 the penalty
is simply a linear transformation of a scalar function, which 0.5
does not convey the non-separable behavior of the penalty for
N >2.
0
A separate document has additional supplemental figures
2
illustrating the convergence of the iterative algorithm (22).
0
These figures show the optimality condition as a scatter plot.
The points in the scatter plot converge to the signum function x2 −2 −2 −1 x 0 1 2
1
as the algorithm converges.
Fig.3. Penaltyψα withα=1forN =2.
7
0.5 ||x||2 0.5 ||x||2
2 2
4 4
3 3
2 2
1 1
0 0
2 2
0 0
x −2 −2 −1 0 1 2 x −2 −2 −1 0 1 2
2 x 2 x
1 1
S (x) S (x)
α α
4 4
3 3
2 2
1 1
0 0
2 2
0 0
x −2 −2 −1 0 1 2 x −2 −2 −1 0 1 2
2 x 2 x
1 1
f (x) = 0.5 ||x||2 − λ S (x) f (x) = 0.5 ||x||2 − λ S (x)
1 2 α 1 2 α
4 4
3 3
2
2
1
1
0
0
2 2
0 0
x −2 −2 −1 0 1 2 x −2 −2 −1 0 1 2
2 x 2 x
1 1
Fig. 4. Differentiable part f1 of cost function with λ = 1, α = 1 for Fig. 5. Differentiable part f1 of cost function with λ = 1, α = 2 for
N =2.Thefunctionisconvexasα(cid:54)1/λ. N =2.Thefunctionisnon-convexasα>1/λ.
