# Sparse Signal Inversion with Impulsive Noise by Dual Spectral Projected Gradient Method

Authors: Liang Ding; Hui Zhao; Yixin Dou
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We consider sparse signal inversion with impulsive noise. There are three major ingredients. The first is regularizing properties; we discuss convergence rate of regularized solutions. The second is devoted to the numerical solutions. It is challenging due to the fact that both fidelity and regularization term lack differentiability. Moreover, for ill-conditioned problems, sparsity regularization is often unstable. We propose a novel dual spectral projected gradient (DSPG) method which combines the dual problem of multiparameter regularization with spectral projection gradient method to solve the nonsmoothl1+l1 optimization functional. We show that one can overcome the nondifferentiability and instability by adding a smooth l2 regularization term to the original optimization functional. The advantage of the proposed functional is that its convex duality reduced to a constraint smooth functional. Moreover, it is stable even for ill-conditioned problems. Spectral projected gradient algorithm is used to compute the minimizers and we prove the convergence. The third is numerical simulation. Some experiments are performed, using compressed sensing and image inpainting, to demonstrate the efficiency of the proposed approach.


---

## Body

## 1. Introduction

In the present manuscript we are concerned with ill-posed linear operator equation:(1)Ax=y,where x is sparse with respect to an orthonormal basis and A:D(A)⊂X→Y is a bounded linear operator. In practice, exact data y are not known precisely, but only an approximation yδ with (2)y-yδ≤δis available. We call yδ the noisy data and δ the noise level. It is well known that the conventional method for solving (1) is sparsity regularization, which provides an efficient way to extract the essential features of sparse solutions compared with oversmoothed classical Tikhonov regularization.In the past ten years, sparsity regularization has certainly become an important concept in inverse problems. The theory of sparse recovery has largely been driven by the needs of applications in compressed sensing [1, 2], bioluminescence tomography [3], seismic tomography [4], parameter identification [5], and so forth. For accounts of the regularizing properties and computational techniques in sparsity regularization we refer the reader to [5–7] and the references given there. In general, sparsity regularization is given by(3)minxAx-yδl22+αxw,pp,where xw,pp=∑γωγφγ,xp(1≤p≤2), α is the regularization parameter balancing the fidelity Ax-yδl22 and regularization term xw,pp. The functional in (3) is not convex if p<1, it is challenging to investigate the regularizing properties and numerical computing method of minimizers. Limited work has been done for p<1; we refer the reader to [8–11] for a recent account of the theory. In this paper, we will focus our main attention on the situation of p=1.The aim of this paper is to consider a regularization functional of the form(4)minxAx-yδl1+α∑γωγφγ,x.We call (4) l1+l1 problem. A main motivation to investigate the l1+l1 problem is that noisy data yδ often contain impulsive noise. For Gaussian noise, l2 fidelity is a natural choice. However, a typical nondifferentiable fidelity used in application involving impulsive noise is the l1 fidelity, which is more robust than l2 fidelity [12].l
1 fidelity technique was motivated by [13] for signal processing and later had attracted considerable attention in image processing, especially in image denoising. In image processing, l1 fidelity typically combines with TV regularization term. For details of l1+TV problems, we refer the reader to [14–17]. Nowadays l1 fidelity has received growing interest in the inverse problems where solutions are sparse with respect to an orthonormal basis. Minimizers of cost functions involving l1 fidelity combined with sparsity regularization have been studied. We refer the reader to [18–23] and the references given there.Forl1 fidelity, regularizing properties must be handled specifically. Burger and Osher [24] proved convergence rate when the regularization parameter is chosen arbitrarily fixed but sufficiently small; the authors call this phenomenon exact penalization. The most curious situation is that, in [25], the optimal convergence rate requires that regularization parameter must be equal to a fixed value. In [26], Grasmair et al. proved that source condition together with finite basis injectivity (FBI) condition is weaker than the restricted isometry property and obtained linear convergence O(δ). Flemming and Hegland discussed convergence rates in l1-regularization when the basis is not smooth enough [27]. König et al. obtained high order polynomial rates of convergence in the special corrupted domain, even though the underlying problem is exponentially ill-posed in the classical sense [28].Thoughl1 fidelity is robust, more researchers prefer to use l2 fidelity because of its differentiability. Hence a key issue for the l1 fidelity is the numerical computing methods. In the past few years, numerous algorithms have been systematically proposed for the l1+TV problems. On the other hand, in spite of growing interests in the l1 fidelity, we can indicate limited work has been done for numerical methods of l1+l1 problems. For sparsity regularization, the popular algorithms, for example, homotopy (LARS) method [29], iteratively reweighted least squares (IRLS) method [30], and iterative thresholding algorithm [31, 32], cannot be directly applied to l1+l1 problem due to the fact that both fidelity and regularization term lack differentiability. There are only a few papers, in which numerical algorithms for l1+l1 problems have been discussed systematically. References [20, 22] proposed the new model by statistics method. Reference [21] discussed the prior sparse representation and the data-fidelity term and proposed the two-phase approach. In [23], authors propose a robust bisparsity model (RBSM) to effectively exploit the prior knowledge about the similarities and the distinctions of signals. However, the researchers often devoted them to compressive sensing problem, where random matrices are well-conditioned. For ill-conditioned problems, these methods are often unstable [33] [Chap. 5]. Moreover, the researchers assume that the solution is sparse itself, which is different from the general assumption that the solution is sparse with respect to an orthonormal basis. In [19], Yang and Zhang proposed a Primal Dual-Interior Point Methods (PD-IPM) for EIT problem, which is efficient at dealing with the nondifferentiability. However, they did not give the convergence proof. Yang and Zhang reformulated the l1+l1 problem into the basis pursuit model which can be solved effectively by ADM method [18]. It is a competitive method compared with other algorithms for compressive sensing. In [34], Xiao et al. applied ADM method to l1+l1 problem directly. Numerical results illustrated that the proposed algorithm performs better than Yall1 [18].In this paper, we investigate regularizing properties ofl1+l1 problems. The convergence rates of xαδ-x†l2 can be shown to be O(δ1-ε) and O(δ) by a priori and a posteriori choice with the additional source condition and FBI condition. As mentioned above, dual is a conventional technique to solve the Tikhonov regularization with l1 fidelity. However, there are some limitations to this approach to l1+l1 problems due to the fact that it is difficult to obtain the dual formulation of l1+l1 problems. Inspired by [14] and multiregularization theory [35–37], a smooth l2 term is added to original functional of regularization. The dual problem of this new cost function is reduced to a constraint smooth functional. Moreover, the smooth l2 regularization term can improve the stability. The conventional method to solve the constraint smooth functional is projected gradient algorithm. We use spectral projected gradient (SPG) method to seek for the minimizers of the constraint functional and prove the convergence.An outline of this paper is as follows. We devote Section2 to a discussion of regularizing properties, including well-posedness and convergence rate. The convergence rate for a priori and a posteriori choice is established. In Section 3, inspired by the theory of multiregularization, we construct a new functional which is equal to a constraint smooth functional according to Frechel duality. In Section 4, the spectral projected gradient method is applied to compute the minimizers. Section 5 provides a detailed exposition of multiparameter choice rule based on the balancing principle. Numerical experiments involving compressed sensing and image inpainting are presented in Section 6, showing that our proposed approaches are robust and efficient.


## 2. Regularization Properties

### 2.1. Notation and Assumptions

For the approximate solutions ofAx=y, we consider the minimization of the regularization functional (5)Jα,δx=Ax-yδl1+αRx,x∈domA∩domR,+∞,x∉domA∩domR,where Rx≔∑γωγφγ,x and the subdifferential of R(x) at x is denoted by ∂R(x)⊂X. All along this paper, X denotes a Hilbert space which is a subspace of l2 space and 〈·,·〉 denotes the inner product; Y denotes a Banach space which is a subspace of l1 space. A:dom(A)⊆X→Y is a bounded linear operator and dom(A)∩dom(R)≠∅. (φγ)γ∈Λ⊂X is an orthonormal basis where Λ is some countable index set. From now, we denote (6)xγ=x,φγ,xl2=∑γφγ,x21/2=∑γxγ21/2.In order to prove convergence rate results we denote byxαδ the minimizer of the regularization functional Jα,δ(x) for every α>0 and use the following definition of R(x)-minimum norm solution.Definition 1.
An elementx† is called a R(x)-minimum norm solution of linear problem Ax=y if (7)Ax†=y,Rx†=minRx∣Kx=y.We define the sparsity as follows.Definition 2.
x
∈
X is sparse with respect to {φγ}γ∈Λ in the sense that supp(x)≔{γ∈Λ:〈φγ,x〉≠0} is finite. If supp(x)0=s for some s∈N, the x∈X is called s-sparse.For a subset of indicesJ⊂Λ, we denote by(8)J=γ∈Λ:φγ,ξ≥ωmin>0,W=spanφγ:γ∈J,where ξ∈∂R(x†). We denote by Jc the complement of J. In addition, let PJx≔∑γ∈Jφγ,xφγ and PJcx≔∑γ∉Jφγ,xφγ.Source condition or variational inequality are necessary for analysis of convergence rate. With these conditions one can obtain convergence rateo(δ1/2), even saturation convergence rate o(δ2/3). But for sparsity regularization, in order to obtain convergence rate o(δ1-ϵ)(0<ϵ<1), even linear convergence o(δ), one needs the following assumptions.Assumption 3.
Assume that the following hold:(i)
Source condition: there existsη∈Y∗ and some ξ∈∂R(x†) satisfying ξ=A∗η.(ii)
Finite basis injectivity (FBI) condition: for every finite setΛ⊆N, the restriction of A to {φγ:γ∈Λ} is injective.The FBI condition was originally introduced in [8]. In [26], Grasmair et al. used a slightly weaker condition based on FBI condition and showed that Assumption 3 is in some sense the weakest possible condition that guarantees the linear convergence rate. Actually, if x† is the unique R(x)-minimizing solution, Assumption 3 is satisfied [26]. For 0≤p<1, x† being the unique R(x)-minimizing solution can not guarantee that Assumption 3 is satisfied. Even if Assumption 3 is satisfied for 0≤p<1, we do not know if linear convergence can be obtained due to fact that traditional proof for linear convergence needs the convexity of R(X).


### 2.2. Well-Posedness and Convergence Rate

In this subsection, well-posedness and convergence rate of the regularization method are given.Proposition 4 (well-posedness).
For everyα>0 and yδ∈Y, minimizing Jα,δ(x) is well-defined, stable, and convergent.Proof.
InJα,δ(x), A is a linear bounded operator, X is Hilbert space, and R(x) is sequentially lower semicontinuous. It is easy to verify Assumption 3.13 [38]; then proof is along the lines of the proofs of Theorem 3.22 (existence), 3.23 (stablity), and 3.26 (convergence) in [38].We now turn to convergence rate. Forl2+l1 problem, Lorenz proved convergence rate O(δ1/2) with source condition and FBI condition [7]. Linear convergence rate O(δ) was improved by Grasmair et al. for nonlinear equation with additional two nonlinear conditions [8]. In [26], Grasmair et al. proved that source condition together with FBI condition is weaker than the restricted isometry property and obtained linear convergence O(δ). Flemming discussed convergence rates in l1-regularization when the basis is not smooth enough [27].Remark 5.
We note that Bregman distance cannot be used as an error measure in this section due to the fact thatRx≔∑γωγφγ,x fails to be strictly convex. We refer reader to [8, 39] for details. In this section we use l2 norm as error measure.The next results show a convergence rate for a priori parameter choice rule.Lemma 6.
Assume that Assumption3 holds. Then there exists a constant C>0 such that (9)x-x†l2≤CAx-Ax†l1+1+CAPJcxl2for all x∈X.Proof.
Sinceξγ=φγ,ξ∈X, we see that (10)limγ→∞ξγ=0.Together with the definition of J, (10) implies that J is a finite set. From Assumption 3, the restriction A to every x∈W is injective. Then for some constant C>0(11)xl2=A-1Axl2≤A-1Axl1≤CAxl1for all x∈W. By the definition of J and x†∈X we have that φγ,x†=0 for every γ∉J and x†∈W. From (11), it is easy to show that(12)PJx-x†l2≤CAPJx-Ax†l1.Then we conclude that (13)x-x†l2=PJx-x†+PJcxl2≤CAPJx-Ax†l1+PJcxl2=CAx-x†-APJcxl1+PJcxl2≤CAx-x†l1+CAPJcxl2+PJcxl2≤CAx-Ax†l1+CA+1PJcxl2,and the proof is completed.Lemma 7.
Assume there existβ1,β2>0 such that(14)Rx-Rx†≥β1x-x†l2-β2Ax-Ax†l1for all x∈dom(A). If αβ2<1, then (15)xαδ-x†l2≤1+αβ2αβ1δ.Proof.
Sincexαδ minimize Jα,δ(x), then(16)Axαδ-yδl1+αRxαδ≤Ax†-yδl1+αRx†.Together with condition (14) this implies the inequality (17)δ≥Ax†-yδl1≥Axαδ-yδl1+αRxαδ-Rx†≥Axαδ-yδl1+αβ1xαδ-x†l2-αβ2Axαδ-Ax†l1≥Axαδ-yδl1+αβ1xαδ-x†l2-αβ2Axαδ-yδl1-αβ2δ≥αβ1xαδ-x†l2-αβ2δ,which proves Lemma 7.Theorem 8.
Assume that Assumption3 holds. Then for the choice α=O(δϵ)(0<ϵ<1), we deduce that (18)xαδ-x†l2=Oδ1-ϵ.Proof.
Letm≔max⁡{φγ,ξ∣γ∉J}; we can estimate (19)PJcxl2=∑γ∉Jφγ,x21/2≤∑γ∉Jφγ,x≤1ωmin-m∑γ∉Jωγ-mφγ,x≤1ωmin-m∑γ∉Jωγφγ,x-φγ,ξφγ,x≤1ωmin-m∑γ∉Jωγφγ,x-ωγφγ,x†-φγ,ξφγ,x-x†=1ωmin-mRx-Rx†-ξ,x-x†=1ωmin-mRx-Rx†+η,Ax†-Ax=1ωmin-mRx-Rx†+ηl∞Ax-Ax†l1.From Lemma 6, we obtain (20)x-x†l2≤CAx-Ax†l1+CA+11ωmin-mRx-Rx†+ηl∞Ax-Ax†l1;that is,(21)Rx-Rx†≥ωmin-mCA+1x-x†l2-ηl∞+Cωmin-mCA+1Ax-Ax†l1.Let(22)β1=ωmin-mCA+1,β2=ηl∞+Cωmin-mCA+1,αβ2<1.According to Lemma 7, it follows that (23)xαδ-x†l2≤1+αβ2αβ1δ.The assertion follows from α=O(δϵ).Next we turn our attention to a posteriori parameter choice which is so-called discrepancy principle due to Morozov [40]. The regularization parameter defined via the discrepancy principle is (24)α≔supα>0∣Axαδ-yδl1=τδfor some τ≥1. Morozov discrepancy principle is a widely used and easy numerical implementation rule. Furthermore, for τ=1, (4) is equivalent to a constrained minimization problem [26] (25)Rx†⟶minsubjecttoAxαδ-yδl1≤τδ,which is widely used in compressive sensing. The next results show a convergence rate by the Morozov discrepancy principle.Theorem 9.
Assume that Assumption3 holds and that the regularization parameter is determined by (24). Then (26)xαδ-x†l2=Oδ.Proof.
Sincexαδ is a minimizer of Jα,δ(x), the inequality (27)Axαδ-yδl1+αRxαδ≤Ax†-yδl1+αRx†≤δ+αRx†holds. This together with Axαδ-yδl1=τδτ≥1 implies that (28)Rxαδ≤Rx†.Proceeding as in the proof of Theorem 8, we obtain that (29)PJcxαδl2≤1ωmin-mRxαδ-Rx†+ηl∞Axαδ-Ax†l1≤ηl∞Axαδ-Ax†l1ωmin-m.From Lemma 6, it follows that (30)xαδ-x†l2≤CAxαδ-Ax†l1+CA+1PJcxαδl2≤CAxαδ-Ax†l1+CA+1ηl∞Axαδ-Ax†l1ωmin-m≤C+ηl∞CA+1ωmin-mAxαδ-Ax†l1≤C+ηl∞CA+1ωmin-mAxαδ-yδl1+Ax†-yδl1≤C+ηl∞CA+1ωmin-mτ+1δ.The proof is completed.Remark 10.
Linear convergence rateO(δ) can also be obtained when the a priori choice and the Morozov discrepancy principle are applied to Kxαδ-yδlpp(1<p≤2) fidelity. However, one must let α=O(δ) in the a priori parameter choice rule.


## 2.1. Notation and Assumptions

For the approximate solutions ofAx=y, we consider the minimization of the regularization functional (5)Jα,δx=Ax-yδl1+αRx,x∈domA∩domR,+∞,x∉domA∩domR,where Rx≔∑γωγφγ,x and the subdifferential of R(x) at x is denoted by ∂R(x)⊂X. All along this paper, X denotes a Hilbert space which is a subspace of l2 space and 〈·,·〉 denotes the inner product; Y denotes a Banach space which is a subspace of l1 space. A:dom(A)⊆X→Y is a bounded linear operator and dom(A)∩dom(R)≠∅. (φγ)γ∈Λ⊂X is an orthonormal basis where Λ is some countable index set. From now, we denote (6)xγ=x,φγ,xl2=∑γφγ,x21/2=∑γxγ21/2.In order to prove convergence rate results we denote byxαδ the minimizer of the regularization functional Jα,δ(x) for every α>0 and use the following definition of R(x)-minimum norm solution.Definition 1.
An elementx† is called a R(x)-minimum norm solution of linear problem Ax=y if (7)Ax†=y,Rx†=minRx∣Kx=y.We define the sparsity as follows.Definition 2.
x
∈
X is sparse with respect to {φγ}γ∈Λ in the sense that supp(x)≔{γ∈Λ:〈φγ,x〉≠0} is finite. If supp(x)0=s for some s∈N, the x∈X is called s-sparse.For a subset of indicesJ⊂Λ, we denote by(8)J=γ∈Λ:φγ,ξ≥ωmin>0,W=spanφγ:γ∈J,where ξ∈∂R(x†). We denote by Jc the complement of J. In addition, let PJx≔∑γ∈Jφγ,xφγ and PJcx≔∑γ∉Jφγ,xφγ.Source condition or variational inequality are necessary for analysis of convergence rate. With these conditions one can obtain convergence rateo(δ1/2), even saturation convergence rate o(δ2/3). But for sparsity regularization, in order to obtain convergence rate o(δ1-ϵ)(0<ϵ<1), even linear convergence o(δ), one needs the following assumptions.Assumption 3.
Assume that the following hold:(i)
Source condition: there existsη∈Y∗ and some ξ∈∂R(x†) satisfying ξ=A∗η.(ii)
Finite basis injectivity (FBI) condition: for every finite setΛ⊆N, the restriction of A to {φγ:γ∈Λ} is injective.The FBI condition was originally introduced in [8]. In [26], Grasmair et al. used a slightly weaker condition based on FBI condition and showed that Assumption 3 is in some sense the weakest possible condition that guarantees the linear convergence rate. Actually, if x† is the unique R(x)-minimizing solution, Assumption 3 is satisfied [26]. For 0≤p<1, x† being the unique R(x)-minimizing solution can not guarantee that Assumption 3 is satisfied. Even if Assumption 3 is satisfied for 0≤p<1, we do not know if linear convergence can be obtained due to fact that traditional proof for linear convergence needs the convexity of R(X).


## 2.2. Well-Posedness and Convergence Rate

In this subsection, well-posedness and convergence rate of the regularization method are given.Proposition 4 (well-posedness).
For everyα>0 and yδ∈Y, minimizing Jα,δ(x) is well-defined, stable, and convergent.Proof.
InJα,δ(x), A is a linear bounded operator, X is Hilbert space, and R(x) is sequentially lower semicontinuous. It is easy to verify Assumption 3.13 [38]; then proof is along the lines of the proofs of Theorem 3.22 (existence), 3.23 (stablity), and 3.26 (convergence) in [38].We now turn to convergence rate. Forl2+l1 problem, Lorenz proved convergence rate O(δ1/2) with source condition and FBI condition [7]. Linear convergence rate O(δ) was improved by Grasmair et al. for nonlinear equation with additional two nonlinear conditions [8]. In [26], Grasmair et al. proved that source condition together with FBI condition is weaker than the restricted isometry property and obtained linear convergence O(δ). Flemming discussed convergence rates in l1-regularization when the basis is not smooth enough [27].Remark 5.
We note that Bregman distance cannot be used as an error measure in this section due to the fact thatRx≔∑γωγφγ,x fails to be strictly convex. We refer reader to [8, 39] for details. In this section we use l2 norm as error measure.The next results show a convergence rate for a priori parameter choice rule.Lemma 6.
Assume that Assumption3 holds. Then there exists a constant C>0 such that (9)x-x†l2≤CAx-Ax†l1+1+CAPJcxl2for all x∈X.Proof.
Sinceξγ=φγ,ξ∈X, we see that (10)limγ→∞ξγ=0.Together with the definition of J, (10) implies that J is a finite set. From Assumption 3, the restriction A to every x∈W is injective. Then for some constant C>0(11)xl2=A-1Axl2≤A-1Axl1≤CAxl1for all x∈W. By the definition of J and x†∈X we have that φγ,x†=0 for every γ∉J and x†∈W. From (11), it is easy to show that(12)PJx-x†l2≤CAPJx-Ax†l1.Then we conclude that (13)x-x†l2=PJx-x†+PJcxl2≤CAPJx-Ax†l1+PJcxl2=CAx-x†-APJcxl1+PJcxl2≤CAx-x†l1+CAPJcxl2+PJcxl2≤CAx-Ax†l1+CA+1PJcxl2,and the proof is completed.Lemma 7.
Assume there existβ1,β2>0 such that(14)Rx-Rx†≥β1x-x†l2-β2Ax-Ax†l1for all x∈dom(A). If αβ2<1, then (15)xαδ-x†l2≤1+αβ2αβ1δ.Proof.
Sincexαδ minimize Jα,δ(x), then(16)Axαδ-yδl1+αRxαδ≤Ax†-yδl1+αRx†.Together with condition (14) this implies the inequality (17)δ≥Ax†-yδl1≥Axαδ-yδl1+αRxαδ-Rx†≥Axαδ-yδl1+αβ1xαδ-x†l2-αβ2Axαδ-Ax†l1≥Axαδ-yδl1+αβ1xαδ-x†l2-αβ2Axαδ-yδl1-αβ2δ≥αβ1xαδ-x†l2-αβ2δ,which proves Lemma 7.Theorem 8.
Assume that Assumption3 holds. Then for the choice α=O(δϵ)(0<ϵ<1), we deduce that (18)xαδ-x†l2=Oδ1-ϵ.Proof.
Letm≔max⁡{φγ,ξ∣γ∉J}; we can estimate (19)PJcxl2=∑γ∉Jφγ,x21/2≤∑γ∉Jφγ,x≤1ωmin-m∑γ∉Jωγ-mφγ,x≤1ωmin-m∑γ∉Jωγφγ,x-φγ,ξφγ,x≤1ωmin-m∑γ∉Jωγφγ,x-ωγφγ,x†-φγ,ξφγ,x-x†=1ωmin-mRx-Rx†-ξ,x-x†=1ωmin-mRx-Rx†+η,Ax†-Ax=1ωmin-mRx-Rx†+ηl∞Ax-Ax†l1.From Lemma 6, we obtain (20)x-x†l2≤CAx-Ax†l1+CA+11ωmin-mRx-Rx†+ηl∞Ax-Ax†l1;that is,(21)Rx-Rx†≥ωmin-mCA+1x-x†l2-ηl∞+Cωmin-mCA+1Ax-Ax†l1.Let(22)β1=ωmin-mCA+1,β2=ηl∞+Cωmin-mCA+1,αβ2<1.According to Lemma 7, it follows that (23)xαδ-x†l2≤1+αβ2αβ1δ.The assertion follows from α=O(δϵ).Next we turn our attention to a posteriori parameter choice which is so-called discrepancy principle due to Morozov [40]. The regularization parameter defined via the discrepancy principle is (24)α≔supα>0∣Axαδ-yδl1=τδfor some τ≥1. Morozov discrepancy principle is a widely used and easy numerical implementation rule. Furthermore, for τ=1, (4) is equivalent to a constrained minimization problem [26] (25)Rx†⟶minsubjecttoAxαδ-yδl1≤τδ,which is widely used in compressive sensing. The next results show a convergence rate by the Morozov discrepancy principle.Theorem 9.
Assume that Assumption3 holds and that the regularization parameter is determined by (24). Then (26)xαδ-x†l2=Oδ.Proof.
Sincexαδ is a minimizer of Jα,δ(x), the inequality (27)Axαδ-yδl1+αRxαδ≤Ax†-yδl1+αRx†≤δ+αRx†holds. This together with Axαδ-yδl1=τδτ≥1 implies that (28)Rxαδ≤Rx†.Proceeding as in the proof of Theorem 8, we obtain that (29)PJcxαδl2≤1ωmin-mRxαδ-Rx†+ηl∞Axαδ-Ax†l1≤ηl∞Axαδ-Ax†l1ωmin-m.From Lemma 6, it follows that (30)xαδ-x†l2≤CAxαδ-Ax†l1+CA+1PJcxαδl2≤CAxαδ-Ax†l1+CA+1ηl∞Axαδ-Ax†l1ωmin-m≤C+ηl∞CA+1ωmin-mAxαδ-Ax†l1≤C+ηl∞CA+1ωmin-mAxαδ-yδl1+Ax†-yδl1≤C+ηl∞CA+1ωmin-mτ+1δ.The proof is completed.Remark 10.
Linear convergence rateO(δ) can also be obtained when the a priori choice and the Morozov discrepancy principle are applied to Kxαδ-yδlpp(1<p≤2) fidelity. However, one must let α=O(δ) in the a priori parameter choice rule.


## 3. Dual Problem

In this section, letωγ in (4) take the same value; that is, ωγ=μ>0 for all γ∈Γ. It is reasonable because convergence can be obtained when δ/α→0 [6]. Let α≔αμ=αωγ; then (4) is equivalent to (31)minxAx-yδl1+α∑γφγ,x.Let u=(x1,x2,…,xγ,…)∈l2, where xγ=〈φγ,x〉. In addition, we denote by D:l2→l2 a dictionary satisfied with u=Dx and u†=Dx†. For example, in the field of wavelet transform, D is a wavelet decomposition operator and DT is a wavelet reconstruction operator [41, 42]. Let K=A∘DT; then (4) is equivalent to (32)P:minuJαu=Ku-yδl1+αul1.Dual is a conventional method for solving Tikhonov regularization with l1 fidelity. However, there are some limitations to this approach to solve (32). The main difficulty is that both the l1 fidelity and the l1 regularization term are nondifferentiable. Moreover, for ill-conditioned problems, sparsity regularization is often unstable. We add the smooth penalty β/2ul22 to (32) to construct the following functional: (33)Pβ:minuJα,βu=Ku-yδl1+αul1+β2ul22.The advantage of problem (33) in place of (32) is that the dual problem of (33) is a constraint smooth functional and projected gradient algorithm can be used to compute minimizers. Moreover, the regularization effect of l1 penalty is weak and l2 penalty can improve the stability of (32).Next we will investigate the convergence of the minimizers to the functionalPβ as β tending to zero.Theorem 11.
Letα be fixed and {βn} be a sequence tending to zero. Then the minimizers uαβn to the functional Jα,βn(u) have a subsequence converging to u∗ being a minimizer of the functional Jα(u).Proof.
By the definition ofuαβn, we have (34)Jα,βnuαβn=Kuαβn-yδl1+αuαβnl1+βn2uαβnl22≤Ku†-yδl1+αu†l1+βn2u†l22≤δ+αu†l1+βn2u†l22,and hence (35)uαβnl1≤δα+u†l1+βn2αu†l22.Since α is a fixed value and βn→0, there exists a constant C>0 such that (36)uαβnl2≤uαβnl1≤δα+u†l1+βn2αu†l22<C.It follows that {uαβn} is uniformly bounded. Therefore, a subsequence {uαβm} of {uαβn} and u∗ exist such that (37)uαβm⟶u∗.By the weak lower semicontinuity of norm, we obtain that (38)Ku∗-yδl1≤liminfKuαβm-yδl1,u∗l1≤liminfuαβml1,and hence that (39)Ku∗-yδl1+αu∗l1≤liminfKuαβm-yδl1+αuαβml1+βm2uαβml22≤limsupKuαβm-yδl1+αuαβml1+βm2uαβml22≤limsupKu-yδl1+αul1+βm2ul22=Ku-yδl1+αul1for all u∈D(K). Therefore, u∗ is a minimizer of the functional Jα(u) and the proof is completed.Next we consider the dual problem ofPβ. We will show that the constraint smooth minimization problems Pβ∗(40)Pβ∗:min12βK∗pl22-p,yδ,s.t.pl∞≤minα,1,are the dual problem of Pβ.Theorem 12.
P
β
∗ is the dual problem of Pβ. The solutions uαβ of Pβ and pαβ of Pβ∗ have the following relation: (41)uαβi=Sα/βK∗pαβiβ,i=1,2,…,Kuαβ-yδ,p-pαβ≥0for all p∈l2 with pl∞≤min{α,1}, where(42)Sα/βx=x-αβ,ifx≥αβ,0,ifx<αβ,x+αβ,ifx≤-αβis the soft threshold function.Proof.
Let(43)Fu=u-yδl1;Ru=β2ul22+αul1,then problem Pβ can be rewritten as (44)infu∈l2FKu+Ru.Let us denote by F∗ and R∗ the conjugate function of F and R. By the Fenchel duality [43] [Chap 3, Chap 10], it follows that (45)R∗p=12βpl22,ifpl∞≤α,∞,else,F∗p=p,yδl2,ifpl∞≤1,∞,else.In Jα,δ(x), the functionals F and R are proper, convex lower semicontinuous. By the Fenchel duality theorem [43](Chap 3, Prop. 2.4, Prop. 4.1, Rem. 4.2), it is easy to show that (46)infu∈l2FKu+Ru=supp∈l2-F∗-p-R∗K∗p;that is,(47)infu∈l2FKu+Ru=-infp∈l2F∗-p+R∗K∗p.Since uαβ and pαβ are minimizers of Pβ and Pβ∗, we have that(48)FKuαβ+Ruαβ+F∗-pαβ+R∗K∗pαβ=0.Moreover, the extremality condition (48) is equivalent to the Kuhn-Tucker conditions(49)-pαβ∈∂FKuαβ,K∗pαβ∈∂Ruαβ.By the definition of the subgradient and -pαβ∈∂F(Kuαβ) in (49), it follows that (50)Kuαβ∈∂F∗-pαβ,Kuαβ-yδ∈∂I-pαβl∞≤1.Then we obtain (51)Kuαβ-yδ,p-pαβ≥0.By the definition of the subgradient and K∗pαβ∈∂R(uαβ) in (49), we conclude that (52)K∗pαβi=βuαβi+αsignuαβi,i=1,2,….For (uαβ)i>0, this leads to (uαβ)i=(K∗(pαβ))i/β-α/β; we must impose in this case that K∗pαβi/β>α/β. For (uαβ)i<0, this leads to (uαβ)i=K∗pαβi/β+α/β, valid only when K∗pαβi/β<-α/β. When |(K∗(pαβ))i/β|<α/β, we put (uαβ)i=0. Summarizing (53)uαβi=Sα/βK∗pαβiβ,where Sα/β is soft threshold function which is defined by (54)Sα/βx=x-αβ,ifx≥αβ,0,ifx<αβ,x+αβ,ifx≤-αβ,the proof is completed.Remark 13.
In [39] (chap 1), numerical experiments showed that the choice of p slightly larger than 1 gives even better results than p=1. It is necessary to discuss the functional (55)minuJαu=Ku-yδl1+αpulpp,where 1<p<2; we call (55) l1+lp problem. We do not have to add additional l2 penalty to (55); its dual problem can be obtained by Fenchel dual theory directly.


## 4. Computation of Minimizers

As mentioned above, the functional equation (33) can be transformed using Fenchel duality into a smooth functional with a box constraint, which can be solved effectively using projected gradient-type methods. The spectral projected gradient (SPG) method [44] has been proved to be effective for the minimization of differentiable functional with box constraints. The function f(p) is defined by (56)fp=12βK∗pl22-p,yδ;then (57)∇fp=1βKK∗p-yδ.By setting (58)Bc≔p∈U∗∣pl∞≤c,wherec=minα,1,we denote by PB(c)(z) the orthogonal projection on the ball B(c)(59)PBcxi,j=cxi,jmaxc,xi,j.Given initial value x0 and a constraint on step length αk, that is, 0<αmin≤αk≤αmax, let parameter γ∈(0,1), and 0<σ1<σ2<1, M>1. We use the convergence criteria given by (60)Pcpk+∇fpk-pkl2≤ϵ(see Algorithm 1). We have the following convergence results.Algorithm 1:Spectral projected gradient method for Pβ∗.
(1) Setγ, ϵ, σ1, σ2, αmin, αmax>0, M≥0, k=1; choose p1∈l2, α1
(2)dk=Pcpk-αk∇fpk-pk
(3) Selectλk using line research Algorithm 2
(4) computepk+1=pk+λkdk, sk=pk+1-pk, qk=∇f(pk+1)-∇f(pk), βk=〈sk,qk〉
(5)
if  βk≤0  then
Set αk+1=αmax
Else
Set αk+1=min⁡αmax,maxαmin,sk,sk/βk
end if
(6)k=k+1
UntilPc(xk+∇f(xk))-xk2≤ϵ or k=kmaxAlgorithm 2:Line research.
(1) Setλ=1
(2) Computep-=pk+λdk
(3)
if
f
p
-
≤
m
a
x
0
≤
j
≤
min
⁡
k
,
M
-
1
f
(
p
k
-
j
)
+
γ
λ
d
k
,
∇
f
p
k
then
Let λk=λ
Else
Let λ∈[σ1λ,σ2λ], and go to (2)
end ifProposition 14.
Forγ∈(0,1) and 0<σ1<σ2<1, algorithm SPG is well-defined; the sequence pk generated by algorithm SPG converges to a stationary point of the functional Pαβ.Proof.
Sincef(p) is convex and has continuous derivatives, the convergence follows from standard results (refer to Theorem 2.1 in [44]).


## 5. Choice of Parameterα and β

The solutionuαβ of Pβ converges theoretically to the solution uαδ of P as β→0. Obviously smaller β is better. However, 1/2βK∗pl22 tends to infinite as β→0. Moreover, the smaller β weakens the regularization effect of the l2 penalty, which leads to instability. If the solution is sparse, we can say that the nonzero coefficients are impulsive parts and the zero coefficients are smooth parts. Multiparameter regularization (61)Pαβ:minu∈l2Jα,βu=Ku-yδl1+αul1+β2ul22is a conventional method for ill-posed problems if the solutions have a number of different structures. In [35, 37], numerical experiments show that multiparameter regularization can effectively recover the different part of solutions. If the solutions contain only a single structure, multiparameter regularization also has better performance. Hence it does not necessarily require the parameter β tending to zero when we choose the parameters α and β. A conventional method for the choice of parameters α and β is multiparameter regularization choice principle, for example, discrepancy principle [35] and balance principle [37]. In this paper, we use multiparameter balance principle for the choice of regularization parameters α and β.Balance principle is to compute minimizers of the function(62)Φγα,β=cγinfKu-yδl1+αul1+β/2ul22γ+2αβ,where cγ is a fixed constant, and (62) is equivalent to [37] (63)γαuα,βδl1=γβuα,βδl22=Kuα,βδ-yδl1.A fixed point algorithm for α and β according to (63) is as in Algorithm 3.Algorithm 3:Fixed point algorithm for α and β.
(1) Chooseα0 and β0, and let k=0
(2)uk+1=argminxKu-yδl1+αul1+β2ul22
(3)αk+1=11+γKuk+1-yδl1+βuk+1l22uk+1l1
β
k
+
1
=
1
1
+
γ
K
u
k
+
1
-
y
δ
l
1
+
β
u
k
+
1
l
1
u
k
+
1
l
2
2
Until a stopping criterion satisfied.


## 6. Numerical Simulations

In this section, we present some numerical experiments to illustrate the efficiency of the proposed method. In Section6.1, numerical experiments involve compressive sensing. We aim to demonstrate that l1 fidelity is more stable than the l2 fidelity and is capable of handing both impulsive and Gaussian noises. In Section 6.2, we first compare the performance of the DSPG method with the alternating direction method (ADM) and TNIP method by well-conditioned compressive sensing problems. In the second example, we discuss an ill-posed problem where the condition number of linear operator A is 255; we aim to demonstrate that the proposed method is stable even with large condition numbers. In Section 6.3, we discuss the image inpainting where images are sparse with respect to the Daubechies wavelets. For image inpainting, the linear operator A is moderate ill-condition and the condition number is around 4000. In order to compare the restoration results, the quality of the computed solution x is measured by relative error Rerr and PSNR which are, respectively, defined by (64)Rerrx=x-x†x†×100%,PSNRyδ=-20log10x-x†n.All experiments were performed under Windows 7 and Matlab R2010a on HP ProBook 4431s with Intel Core i5 2410M CPU 2.30 GHz 2.30 GHz and 4 GB of memory.


### 6.1. Comparison withl1 and l2 Fidelity

This example involves compressive sensing problemAx=yδ, where matrix A80×200 is random Gaussian and yδ=Ax†+δ is the observed data containing white noise or impulsive noise. The true solution x† is 16-sparse with respect to natural basis of l2 space which is defined by (65)φγ=eγ=0,0,…,0,1︸γ,0,….White noise is generated such that data yδ attains a desired SNR, which is defined by (66)SNRyδ=20log10yδ-Eyδδ.The impulsive noise is measured by relative error, which is defined by (67)Rerrδ=y-yδy×100%.Figure 1 presents the restoration results by l1 and l2 fidelities. The left column describes the data with different white noise levels. The SNR of data are 30 dB, 20 dB, 10 dB, and 5 dB. The right column is restoration results; the regularization parameters (α,β) are (2.37e-3,1.26e-3), (6.32e-2,2.63e-2), (4.35e-1,1.68e-1), and (9.05e-1,7.45e-1). As can be seen from Figure 1, the quality of restoration by l1 and l2 fidelities is similar if the data are contaminated with low noise levels. If data contain high noise levels, the performance of l2 fidelity is better. However, we can not recover high-quality solutions with high noise levels. In Figure 2, the left column describes the data which are contaminated by different impulsive noise levels. Rerr(δ) of noise level are 3%, 7%, 15%, and 22%. The value of impulsive noise is ±1 at random positions and 0 at other positions. The right column contains restoration results according to different noise levels. The regularization parameters (α,β) are (2.43e-2,1.58e-2), (8.29e-2,4.63e-2), (2.12e-1,1.06e-1), and (4.25e-1,2.46e-1). As can be seen from Figure 2, the l1 fidelity is more stable for impulsive noise and always offers high-quality restoration even with poor data. In contrast to the l2 fidelity, the quality of restoration results by the l2 fidelity is always poor.Figure 1
Restorations ofl1 and l2 fidelities with different Gaussian noise levels.Figure 2
Restorations ofl1 and l2 fidelities with different impulsive noise levels.


### 6.2. Comparison of DSPG with ADM-l1, ADM-l2, and TNIP

In this subsection, we present the comparison results of DSPG with ADM-l1, ADM-l2, and TNIP. ADM-l1 [18] is an efficient alternating direction method for l1+l1 problem. ADM-l2 [18] and TNIP [45] (truncated Newton interior point) are for l2+l1 problem. In the first experiment, we use random Gaussian matrix Am×n, where sampling length is m and signal length n=200. The condition number of random Gaussian matrix Am×n is around 5. The signal is p-sparsity. We add 1% impulsive noise to true signal. For each fixed pair (m,p), we take 100 iterations, where m/n=0.5, 0.4, 0.3, 0.2 and 0.1, p/m=0.1 and 0.2. Figure 3 presents comparison results for four different iterative algorithms. The left column describes convergence rates of Rerr(x) for p/m=0.1. The right column describes convergence rates of Rerr(x) for p/m=0.2. m/n increase from top to bottom row. As can be seen from Figure 3, ADM-l1 converges faster than DSPG and the accuracy is better than DSPG when m/n=0.4. With m/n decreasing, DSPG method performs more competitively. The accuracy of DSPG method is even better than ADM-l1 when m/n=0.1. Though the optimal relative error of ADM-l1 is better than DSPG method, the corresponding optimal iteration number or stopping tolerance of ADM-l1 is difficult to estimate in practice. The other l2 fidelity algorithms, especially TNIP method, converge obviously faster than ADM-l1 and DSPG method. However, the accuracy of the two algorithms is worse than ADM-l1 and DSPG method no matter what value the m/n is.Figure 3
Comparisons of DSPG with ADM-l1, ADM-l2, and TNIP.Next, in order to test the stability of the DSPG method for ill-conditioned problems, we use matrixAn×n (n=200) whose condition number is 255. This problem was discussed by Lorenz in [8] where the ill-conditioned matrix is generated by Matlab code: A = tril(ones(200)). The signal is p-sparsity where p/n=0.1 and 0.2. We add 1% impulsive noise to data. As can be seen from Figure 4, DSPG method converges obviously faster than the ADM-l1 method. The relative error of DSPG method is also better than ADM-l1 method. It is shown that DSPG method is stable even for large condition number matrices. In Table 1, data contain impulsive noise with corruption Rerr(δ)=0.1%, 0.3%, 0.5%, 1%, 3%, 5%, 10%. As can be seen from Table 1, the quality of restoration improves as noise level δ decreases. Theoretically, ADM-l1 and DSPG methods are adept to process impulsive noise. However, ADM-l1 method is sensitive to noise when the operators are ill-conditioned. In this case, ADM-l1 cannot obtain reasonable restoration. ADM-l2 and DSPG methods are more stable to noise level δ even if matrix A has large condition numbers. For low noise levels, ADM-l2 and DSPG methods have advantage over the other two methods. For high noise levels, the restoration results of ADM-l1 method are very close to ADM-l2 and TNIP methods. Restoration results of the DSPG method are obviously better than the other three methods.Table 1
Comparisons of different algorithms.
Impulsive noiseδ
TNIP
AMD-l2
AMD-l1
DSPG
Rerror
PSNR
Rerror
PSNR
Rerror
PSNR
Rerror
PSNR
0.1%
0.610%
54.49
0.399%
57.97
0.907%
50.84
0.387%
57.28
0.3%
1.412%
47.01
0.559%
55.04
1.813%
44.83
0.548%
55.36
0.5%
2.295%
42.78
1.078%
49.34
2.719%
41.31
0.956%
50.14
1.0%
3.187%
39.92
3.178%
39.93
3.619%
38.82
2.305%
42.65
3.0%
4.993%
36.03
4.956%
36.12
5.542%
35.31
3.063%
40.15
5.0%
8.608%
31.30
8.602%
31.30
9.101%
30.87
5.217%
36.25
10%
17.66%
25.06
17.66%
25.08
18.11%
24.84
12.64%
29.21Figure 4
Comparisons of DSPG with ADM-l1, ADM-l2, and TNIP.
(a)
p
/
n
=
0.1
(b)
p
/
n
=
0.2


### 6.3. Image Inpainting

We consider 2D image inpainting problems. The image is Barbara (n=512; cf. Figure 5). We randomly remove 40% pixels of Barbara to create an incomplete image. In this case, the image inpainting is a moderate ill-posed problem. The condition number of the matrix is around 4000. For our purpose, we make use of Daubechies 4 wavelet basis as a dictionary. We use four scales, for a total of 8192×512 wavelet and scaling coefficients (cf. Figure 6). As seen from Figure 6, the representation of the image with respect to Daubechies 4 basis is sparse. We add salt-and-pepper noise by Matlab code: imnoise(image, ‘salt & pepper’, d). In this example, d=0.05. The restoration results are shown in Figure 5(d). Restoration results of four images with different noise levels are given in Table 2. Restoration results show that if images have a sparse representation with respect to an orthogonal basis, DSPG method can obtain satisfactory results even if the image inpainting are moderate ill-posed problems.Table 2
Restoration of four images with different salt-and-pepper noise levels.
Noise level
Lena
Boat
Barbara
Goldhill
Cameraman
Peppers
Mandrill
Pirate
0.05
26.21
25.16
22.16
25.58
25.12
26.25
26.68
25.09
0.08
24.87
23.88
21.52
25.17
25.03
26.11
26.16
25.01
0.10
22.82
21.56
20.35
24.07
24.23
25.01
25.32
24.58
0.15
20.42
19.55
19.23
22.74
22.42
23.16
23.85
22.75Figure 5
Restoration of Barbara with salt-and-pepper noise,d=0.2.
(a)
True image
(b)
Painted image
(c)
Painted image with 0.2 salt-and-pepper noise
(d)
RestorationFigure 6
Wavelet and scale coefficients.


## 6.1. Comparison withl1 and l2 Fidelity

This example involves compressive sensing problemAx=yδ, where matrix A80×200 is random Gaussian and yδ=Ax†+δ is the observed data containing white noise or impulsive noise. The true solution x† is 16-sparse with respect to natural basis of l2 space which is defined by (65)φγ=eγ=0,0,…,0,1︸γ,0,….White noise is generated such that data yδ attains a desired SNR, which is defined by (66)SNRyδ=20log10yδ-Eyδδ.The impulsive noise is measured by relative error, which is defined by (67)Rerrδ=y-yδy×100%.Figure 1 presents the restoration results by l1 and l2 fidelities. The left column describes the data with different white noise levels. The SNR of data are 30 dB, 20 dB, 10 dB, and 5 dB. The right column is restoration results; the regularization parameters (α,β) are (2.37e-3,1.26e-3), (6.32e-2,2.63e-2), (4.35e-1,1.68e-1), and (9.05e-1,7.45e-1). As can be seen from Figure 1, the quality of restoration by l1 and l2 fidelities is similar if the data are contaminated with low noise levels. If data contain high noise levels, the performance of l2 fidelity is better. However, we can not recover high-quality solutions with high noise levels. In Figure 2, the left column describes the data which are contaminated by different impulsive noise levels. Rerr(δ) of noise level are 3%, 7%, 15%, and 22%. The value of impulsive noise is ±1 at random positions and 0 at other positions. The right column contains restoration results according to different noise levels. The regularization parameters (α,β) are (2.43e-2,1.58e-2), (8.29e-2,4.63e-2), (2.12e-1,1.06e-1), and (4.25e-1,2.46e-1). As can be seen from Figure 2, the l1 fidelity is more stable for impulsive noise and always offers high-quality restoration even with poor data. In contrast to the l2 fidelity, the quality of restoration results by the l2 fidelity is always poor.Figure 1
Restorations ofl1 and l2 fidelities with different Gaussian noise levels.Figure 2
Restorations ofl1 and l2 fidelities with different impulsive noise levels.


## 6.2. Comparison of DSPG with ADM-l1, ADM-l2, and TNIP

In this subsection, we present the comparison results of DSPG with ADM-l1, ADM-l2, and TNIP. ADM-l1 [18] is an efficient alternating direction method for l1+l1 problem. ADM-l2 [18] and TNIP [45] (truncated Newton interior point) are for l2+l1 problem. In the first experiment, we use random Gaussian matrix Am×n, where sampling length is m and signal length n=200. The condition number of random Gaussian matrix Am×n is around 5. The signal is p-sparsity. We add 1% impulsive noise to true signal. For each fixed pair (m,p), we take 100 iterations, where m/n=0.5, 0.4, 0.3, 0.2 and 0.1, p/m=0.1 and 0.2. Figure 3 presents comparison results for four different iterative algorithms. The left column describes convergence rates of Rerr(x) for p/m=0.1. The right column describes convergence rates of Rerr(x) for p/m=0.2. m/n increase from top to bottom row. As can be seen from Figure 3, ADM-l1 converges faster than DSPG and the accuracy is better than DSPG when m/n=0.4. With m/n decreasing, DSPG method performs more competitively. The accuracy of DSPG method is even better than ADM-l1 when m/n=0.1. Though the optimal relative error of ADM-l1 is better than DSPG method, the corresponding optimal iteration number or stopping tolerance of ADM-l1 is difficult to estimate in practice. The other l2 fidelity algorithms, especially TNIP method, converge obviously faster than ADM-l1 and DSPG method. However, the accuracy of the two algorithms is worse than ADM-l1 and DSPG method no matter what value the m/n is.Figure 3
Comparisons of DSPG with ADM-l1, ADM-l2, and TNIP.Next, in order to test the stability of the DSPG method for ill-conditioned problems, we use matrixAn×n (n=200) whose condition number is 255. This problem was discussed by Lorenz in [8] where the ill-conditioned matrix is generated by Matlab code: A = tril(ones(200)). The signal is p-sparsity where p/n=0.1 and 0.2. We add 1% impulsive noise to data. As can be seen from Figure 4, DSPG method converges obviously faster than the ADM-l1 method. The relative error of DSPG method is also better than ADM-l1 method. It is shown that DSPG method is stable even for large condition number matrices. In Table 1, data contain impulsive noise with corruption Rerr(δ)=0.1%, 0.3%, 0.5%, 1%, 3%, 5%, 10%. As can be seen from Table 1, the quality of restoration improves as noise level δ decreases. Theoretically, ADM-l1 and DSPG methods are adept to process impulsive noise. However, ADM-l1 method is sensitive to noise when the operators are ill-conditioned. In this case, ADM-l1 cannot obtain reasonable restoration. ADM-l2 and DSPG methods are more stable to noise level δ even if matrix A has large condition numbers. For low noise levels, ADM-l2 and DSPG methods have advantage over the other two methods. For high noise levels, the restoration results of ADM-l1 method are very close to ADM-l2 and TNIP methods. Restoration results of the DSPG method are obviously better than the other three methods.Table 1
Comparisons of different algorithms.
Impulsive noiseδ
TNIP
AMD-l2
AMD-l1
DSPG
Rerror
PSNR
Rerror
PSNR
Rerror
PSNR
Rerror
PSNR
0.1%
0.610%
54.49
0.399%
57.97
0.907%
50.84
0.387%
57.28
0.3%
1.412%
47.01
0.559%
55.04
1.813%
44.83
0.548%
55.36
0.5%
2.295%
42.78
1.078%
49.34
2.719%
41.31
0.956%
50.14
1.0%
3.187%
39.92
3.178%
39.93
3.619%
38.82
2.305%
42.65
3.0%
4.993%
36.03
4.956%
36.12
5.542%
35.31
3.063%
40.15
5.0%
8.608%
31.30
8.602%
31.30
9.101%
30.87
5.217%
36.25
10%
17.66%
25.06
17.66%
25.08
18.11%
24.84
12.64%
29.21Figure 4
Comparisons of DSPG with ADM-l1, ADM-l2, and TNIP.
(a)
p
/
n
=
0.1
(b)
p
/
n
=
0.2


## 6.3. Image Inpainting

We consider 2D image inpainting problems. The image is Barbara (n=512; cf. Figure 5). We randomly remove 40% pixels of Barbara to create an incomplete image. In this case, the image inpainting is a moderate ill-posed problem. The condition number of the matrix is around 4000. For our purpose, we make use of Daubechies 4 wavelet basis as a dictionary. We use four scales, for a total of 8192×512 wavelet and scaling coefficients (cf. Figure 6). As seen from Figure 6, the representation of the image with respect to Daubechies 4 basis is sparse. We add salt-and-pepper noise by Matlab code: imnoise(image, ‘salt & pepper’, d). In this example, d=0.05. The restoration results are shown in Figure 5(d). Restoration results of four images with different noise levels are given in Table 2. Restoration results show that if images have a sparse representation with respect to an orthogonal basis, DSPG method can obtain satisfactory results even if the image inpainting are moderate ill-posed problems.Table 2
Restoration of four images with different salt-and-pepper noise levels.
Noise level
Lena
Boat
Barbara
Goldhill
Cameraman
Peppers
Mandrill
Pirate
0.05
26.21
25.16
22.16
25.58
25.12
26.25
26.68
25.09
0.08
24.87
23.88
21.52
25.17
25.03
26.11
26.16
25.01
0.10
22.82
21.56
20.35
24.07
24.23
25.01
25.32
24.58
0.15
20.42
19.55
19.23
22.74
22.42
23.16
23.85
22.75Figure 5
Restoration of Barbara with salt-and-pepper noise,d=0.2.
(a)
True image
(b)
Painted image
(c)
Painted image with 0.2 salt-and-pepper noise
(d)
RestorationFigure 6
Wavelet and scale coefficients.


## 7. Conclusion

With source and FBI conditions, we have proved thatl1+l1 regularization method yields convergence rates of order δ1-ϵ and δ. For numerical solutions, we have proposed a novel DSPG approach for sparsity regularization with l1 fidelity. Numerical results indicate that the proposed DSPG algorithm performs competitively with several state-of-art algorithms such as ADM method. On various classes of test problems with different condition numbers, the proposed DSPG method has exhibited the following advantages: (i) compared with other algorithms, the dual problem of our methods is more simple which is a box-constraint smooth functional and can be solved effectively by projected methods; (ii) for ill-conditioned problems, DSPG method is more stable with respect to noise level. We remark that for well-conditioned problems, for example, compressive sensing, optimal accuracy of ADM-l1 is better than DSPG method. However, the optimal accuracy of ADM-l1 method is strongly dependant on stopping tolerance values which can be difficult to estimate in practice. To the best of the author’s knowledge, multiparameter regularization is first used to obtain the dual formulation of l1+l1 problems. One can also try to use multiparameter regularization strategy to solve dual problem of the other nonsmooth functionals, for example, l1+TV problems.



---

*Source: 2736306-2017-09-14.xml*