# A Smoothedl0-Norm and l1-Norm Regularization Algorithm for Computed Tomography

Authors: Jiehua Zhu; Xiezhang Li
Journal: Journal of Applied Mathematics
 (2019)

---

## Abstract

The nonmonotone alternating direction algorithm (NADA) was recently proposed for effectively solving a class of equality-constrained nonsmooth optimization problems and applied to the total variation minimization in image reconstruction, but the reconstructed images suffer from the artifacts. Though by thel0-norm regularization the edge can be effectively retained, the problem is NP hard. The smoothed l0-norm approximates the l0-norm as a limit of smooth convex functions and provides a smooth measure of sparsity in applications. The smoothed l0-norm regularization has been an attractive research topic in sparse image and signal recovery. In this paper, we present a combined smoothed l0-norm and l1-norm regularization algorithm using the NADA for image reconstruction in computed tomography. We resolve the computation challenge resulting from the smoothed l0-norm minimization. The numerical experiments demonstrate that the proposed algorithm improves the quality of the reconstructed images with the same cost of CPU time and reduces the computation time significantly while maintaining the same image quality compared with the l1-norm regularization in absence of the smoothed l0-norm.


---

## Body

## 1. Introduction

Statistical and iterative reconstruction algorithms in computed tomography (CT) are widely applied since they yield more accurate results than analytic approaches for low-dose and limited-view reconstruction. These algorithms involve solving a linear system:(1)Φf+e=u,where Φ is an m×n2 projection matrix, f∈Rn2 represents a 2D n×n image to be reconstructed, e is an additive noise with e2≤ε for some known ε>0, and u∈Rm is the noisy projection data. For limited-view reconstruction, the underdetermined system (m≪n2) has infinitely many solutions. An optimal solution representing the original image as well as possible is sought by iteration methods.The theory of compressed sensing [1–3] has recently shown that signals and images that have sparse representations in some orthonormal bases can be reconstructed at high quality from much less data than what the Nyquist sampling theory [4] requires. In many cases in tomography, images can be approximately modeled to be essentially piecewise constant so the gradients are sparse. With the gradient operator as a sparse transform, compressed sensing provides a novel framework for CT image reconstruction [1, 3, 5]. The l0-norm of a vector, the number of its nonzero elements, is an appropriate measurement for the sparsity of a vector. So the l0-norm regularization might be an approach for the reconstruction of an image with sparse gradient. However, the l0-norm is a nonconvex function and the l0-norm regularization problem is NP hard [6]. A pseudoinverse transform of the discrete gradient and the iterative hard thresholding algorithm are used to address the issues of the l0-norm [7]. The smoothed l0-norm (sl0-norm) provides a smooth measure of sparsity and is applied in compressed sensing MRI imaging [8]. The sl0-norm is used to find the jointly sparse representation via the low-resolution image [9]. sl0-norm regularization model is proposed for sparse-view X-ray CT reconstruction [10]. Two new smoothed functions approximating the l0-norm are proposed in the mechanism of reweighted regularization [11, 12]. An edge-preserving image reconstruction method based on l0-regularized gradient prior for limited-angle CT is investigated [13]. The local and global minimizers of l0 gradient regularized model with box constraints for image restoration are discussed [14]. A new l0 regularization and wavelet tight framelets to suppress the slope artifacts in the limited-angle X-ray CT reconstruction are addressed [15]. An image model based on l0 and l2 regularizations for the limited-angle CT is proposed, and the existence of a solution and a local convergence analysis under certain conditions are proved [16].Thel1-norm regularization, known as total variation (TV) minimization [17], has been widely used in CT reconstruction. Under certain conditions, the l1-norm minimization has a unique sparsest solution which is a good approximation to the original image [18]. However, the l1-norm regularization provides less sparsity representation than the l0-norm regularization and may lose some detailed features and contrast. In addition, the discrete gradient transformation ∇f does not satisfy the restricted isometry property required by the compressed sensing theory.Researchers in this area have been seeking for other efficient and stable algorithms inspired by the compressed sensing theory, for example, solving anl1-norm minimization with an l2-norm constraint using the alternating direction method (ADM) [19, 20] and the nonmonotone alternating direction method (NADA) [21, 22]. Regularization including multiple norms is developed to address the drawbacks of the l1-norm and l0-norm. A combined l1-norm and sl0-norm regularization minimization with an l2-norm constraint using SART algorithm and the gradient decent method is proposed for sparse-view CT image reconstruction in [23]. It is observed in the paper that the convergence is slow and the computation is time consuming because of the alternative minimization of the l1-norm and sl0-norm.The purposes of this paper are multifold:(i) Presenting a combinedsl0-norm and l1-norm regularization algorithm for image reconstruction in CT: The proposed algorithm is unique in the following aspects: A new parameter is introduced to balance the sl0-norm and l1-norm terms; two-norm regularization is combined in one Lagrangian object function to be minimized(ii) Adopting a newly developed alternating direction method NADA to efficiently solve minimization(iii) Resolving the computation challenge problem caused from thesl0-norm minimizationThe numerical experiments demonstrate that the proposed algorithm improves the quality of the recovered images for the same cost of CPU time and reduces the computation time significantly while maintaining the same image quality compared with thel1-norm regularization in absence of the sl0-norm.The rest of the paper is organized as follows. Section2 includes the background and notations. In Section 3, a combined sl0-norm and l1-norm regularization algorithm with the NADA is developed. Numerical experiments with the Shepp-Logan phantom and a cardiac image are presented in Section 4. Finally, Section 5 concludes the paper by discussions and conclusions.


## 2. Background and Notations

The TV minimization of an imagef for solving system (1) is mathematically represented as [17](2)minfTVsubjecttoΦf+e=u, where the total variation fTV is the l1-norm of the magnitude of the discrete gradients:(3)fTV=∇f1=∑i,jfi+1,j-fi,j2+fi,j+1-fi,j2.Then a 2D image f with sparse gradients and noisy measurements in CT can be reconstructed by solving the following l1-norm minimization with an l2-norm constraint:(4)minf∇f1subjecttoΦf-u2≤ε,which can be implemented [24] by(5)minfμ2Φf-u22+∇f1.Here μ>0 is a regularization parameter which controls the trade-off between sparsity of gradients and accuracy of approximation.For convenience, we denoteDpf as the forward difference of f at a pixel p in both horizontal and vertical directions; i.e.,(6)Dpf=Di,jf=fi+1,j-fi,j,fi,j+1-fi,jT∈R2,forp=i-1n+j. Then ∇f1=∑p=1n2‍Dpf2 and minimization (5) is rewritten as(7)minfμ2Φf-u22+∑p=1n2Dpf2.With the Lagrangian function [19, 20](8)Lf,v,λ=μ2Φf-u22+∑p=1n2vp2-λpTDpf-vp+β2Dpf-vp22,minimization (7) is converted to(9)minfLf,v,λ,which can be solved by the ADM [19–22]. The ADM is a variant of the classic augmented Lagrangian method for optimization. It has been applied to solve different types of l1-minimization problems for sparse solution recovery. The kth step of the ADM for solving (9) involves the procedures(10)fk+1=argminfLf,vk,λk;(11)vk+1=argminvLfk+1,v,λk;(12)λk+1=λk+βDfk+1-vk+1.Minimization (10) and minimization (11) can be successfully solved using the framework of the NADA with a nonmonotone line search scheme recently developed in [21, 22]. The idea is outlined as follows.(i) Choose a directiondk=-∂/∂fLfk,vk,λk.(ii) Select a step sizesk uniformly bounded above such that L(fk+skdk,vk,λk)≤Ck-skδdk22, for 0<δ<1.(iii) Setfk+1=fk+skdk.(iv) ComputeCk+1 such that L(fk+1,vk,λk)≤Ck+1≤Ck.Combined with the solution of (11), the sequence {L(fk,vk,λk)} is bounded above by a monotonically nonincreasing sequence {Ck} though {L(fk,vk,λk)} itself is not decreasing. The advantage of the NADA lies in the fact that L(f,v,λ) is not required to be differentiable while reducing the computation complexity and improving the efficiency. The convergence analysis of the NADA for solving a general model including minimization (11) can be found in [21].Thesl0-norm approximates the l0-norm by a smooth function [8]. The sl0-norm of a vector x∈Rn2, denoted by xs0, is defined as(13)xs0=n2-∑j=1n2limσ→0exp⁡-xj22σ2.It is easy to see that xs0=x0.A new approach for solving a more general minimization than the regularization minimization in [23] is proposed in the next section.


## 3.sl0-Norm and l1-Norm Regularization Algorithm

In this section we presentsl0-norm and l1-norm regularization algorithm for CT reconstruction using the NADA. Consider the following constrained minimization problem for image reconstruction in CT:(14)minf∇f1+α∇fs0subjecttoΦf-u22≤ε,where the parameter α is used to balance the two terms of the minimization. While the l1-norm is used for the approximation accuracy, the sl0-norm is introduced in consideration of the sparsity of gradients. It is remarked that minimization (14) is an extension of minimization (4) (α=0) and the minimization (α=1) in [23].


### 3.1. New Approach

The term∇fs0 in (14) can be approximated by a smooth function(15)sσ∇f≔n2-∑p=1n2exp⁡-∇fp22σ2=n2-∑p=1n2exp⁡-Dpf222σ2for a small positive parameter σ. So minimization (14) can be approximately obtained by solving(16)minf∇f1+αsσ∇fsubjecttoΦf-u22≤ε,smallσ>0.To deal with the constraint Φf-u22≤ε in (16), we rewrite (16) with a positive parameter μ as(17)minfμ2Φf-u22+∇f1+αsσ∇f,forsmallσ>0, or(18)minfμ2Φf-u22+∑p=1n2vp2+αsσvp2,forsmallσ>0,subjecttoDpf=vp,1≤p≤n2.Using a positive parameter β and Lagrange vectors λp∈R2, 1≤p≤n2, we define a Lagrangian function (19)Lαf,v,λ=μ2Φf-u22+∑p=1n2vp2-λpTDpf-vp+β2Dpf-vp22+αsσvp2.Finally, minimization (18) is converted to the following minimization:(20)minfLαf,v,λ.It is obvious that Lα(f,v,λ) is an extension of L(f,v,λ) in (8). Consequently the NADA is adopted to solve minimization (20). The kth iteration for solving (20) includes the following steps:(21)fk+1=argminfLαf,vk,λk,(22)vk+1=argminvLαfk+1,v,λk,(23)λk+1=argminλLαfk+1,vk+1,λ.It is noted that the objective functionLα(f,v,λ) in (19) has an extra term αsσ(vp), depending only on the vector v, compared with L(f,v,λ) in (8). So the methods for the calculation of fk+1 in (10) and λk+1 in (12) can be applied to the calculation of fk+1 in (21) and λk+1 in (23). For example, in (23), λk+1=λk+β(Dfk+1-vk+1). However, finding vk+1 in (22) is a challenging problem.


### 3.2. Calculation ofvpk+1

In this subsection we address how to calculatevpk+1 in (22) for 1≤p≤n2. Actually, for each p, we need to solve for vp from the following equation:(24)∂Lα∂vp=vpvp2+λp-βDpf+βvp+αvpσ2exp⁡-vp222σ2=0.It follows from (24) that(25)vp1+1βvp2+αβσ2exp⁡-vp222σ2=Dpf-λpβ. If vp2 is calculated then(26)vp=vp2Dpf-λp/βDpf-λp/β2.Thus, it suffices to determine vp2 from the following equation:(27)vp21+αβσ2exp⁡-vp222σ2=Dpf-λpβ2-1β.Denote c=Dpf-λp/β2-1/β. If c≤0 we set vp=0. If c>0, it follows from (27) that vp2 is a positive real solution of the equation(28)hz≔z+αzβσ2exp⁡-z22σ2-c=0.It is obvious that there is a positive zero z of h(z) between 0 and c and z→c as σ→0.Equation (28) can be rewritten as(29)exp⁡-z22σ2=βσ2c-zαzorz22σ2=ln⁡αzβσ2c-z. By substituting z=w, the above equation becomes(30)sw≔w-2σ2ln⁡αwβσ2c-w=0,on0,c2.Thus, we have shown the following lemma.Lemma 1.
Leth(z) and s(w) be given in (28) and (30), respectively. Then w is a zero of s(w) on (0,c2) if and only if z=w is a zero of h(z) on (0,c).Before we develop a procedure to calculate a positive root ofs(w)=0, we will list some properties of s(w) below.It is known that(31)s′w=1-σ2cwc-w=-ww+cw-σ2cwc-w. With z=w, the sign of s′(w) on (0,c2) is the same as that of the function(32)qz≔-z3+cz2-σ2con (0,c). By calculus, it is easy to see that q(z) is increasing on (0,2/3c) and decreasing on (2/3c,c) since q′(z)=-z(3z-2c). Thus, q(z) achieves a local maximum 4/27c3-σ2c at z=2/3c. Therefore, if 4/27c3-σ2c>0 (or σ<2/33c) then the equation q(z)=0 has two positive roots z1,z2 with z1<z2. It follows that q(z)<0 on (0,z1)∪(z2,c) and q(z)>0 on (z1,z2). On the other hand, if σ>2/33c, then q(z)<0 on (0,c). We have the following monotonicity of s(w).Property 2.
Ifσ<2/33c then s(w) is increasing on (z12,z22) and decreasing on (0,z12) and (z22,c2), where z1 and z2 are two positive zeros of q(z). If σ≥2/33c then s(w) is always decreasing on (0,c2).The second derivative(33)s′′w=σ2cw2c-w2c-32w indicates the following features of s(w).Property 3.
s ′ ′ ( w ) > 0 on (0,4/9c2) and s′′(w)<0 on (4/9c2,c2). The graph of s(w) has an inflection point (4/9c2, 4/9c2-2σ2ln⁡2α/βσ2). The function s(w) achieves a local minimum value at w=z12 and a local maximum value at w=z22. In addition, we have limw→0+s(w)=∞ and limw→c2-s(w)=-∞.Finally, we can determine the locations of zeros ofs(w) by convergent Newton’s iterations.Theorem 4.
Lets(w) and q(z) be, respectively, given in (30) and (32), where z=w. Let z2 be a larger positive root of q(z) on (0,c) for σ<2/33c. Then a positive root of s(w)=0 on (0,c2) can be calculated by Newton’s method. The convergence of Newton’s method is guaranteed if an initial guess w0 is selected as follows.
Case 1. (i) σ<2/33c, s(z22)>0 or (ii) σ≥2/33c, 2/9c2-σ2ln⁡2α/βσ2>0.
Choosew0<c2 and near c2 with s(w0)<0.
Case 2. (i) σ<2/33c, s(z22)<0 or (ii) σ≥2/33c, 2/9c2-σ2ln⁡2α/βσ2<0.
Choosew0=(βσ2c/α+βσ2)2.Proof.
We identify an interval containing a positive root ofs(w)=0 on which s′(w) and s′′(w) do not change signs for each subcase. Then the convergence of Newton’s method with an initial guess above is guaranteed.
Case 1. For subcase (i), there is a positive solution of s(w)=0 on (z22,c2) on which both s′(w) and s′′(w) are negative. For subcase (ii), the inflection point P of the graph of s(w) is on the upper half plane. So a positive zero is on (4/9c2,c2) on which both s′(w) and s′′(w) are negative. So one should choose w0<c2 and near c2 with s(w0)<0 for the convergence of Newton’s method.
Case 2. For subcase (i), there is a positive solution of s(w)=0 on (0,z12) on which s′(w)<0 and s′′(w)>0. For subcase (ii), the inflection point P is on the lower half plane. So a positive zero is on (0,4/9c2) on which s′(w)<0 and s′′(w)>0. One could choose w0=(βσ2c/α+βσ2)2 for the convergence of Newton iteration since βσ2c/α+βσ2 is a lower bound of positive solutions of h(z)=0 on (0,c) in (28).Thus,vpk+12, 1≤p≤n2, can be calculated from Lemma 1 and Theorem 4. Then vk+1 in (26) can be obtained.


### 3.3. Algorithm

Based on the above discussion, the calculation process of the proposedsl0+l1 regularization algorithm is summarized in Algorithm 1.Algorithm 1:Calculation process of sl0+l1 regularization algorithm (sl0-norm and l1-norm regularization algorithm).
inputΦ,u,α,ε initializeβ,μ,f0,v0,λ0,maxit,σ,δ,ratio,tol whilek<maxit 1. updatefk+1=argminfLα(f,vk,λk) by the NADA 1.1. decent directiondk=-∂/∂fLα(fk,vk,λk) 1.2. step sizesk such that Lα(fk+skdk,vk,λk)≤Ck-skδdk22 1.3. updatefk+1=fk+skdk 1.4. computeCk+1 such that Lα(fk+1,vk,λk)≤Ck+1≤Ck 2. updatevk+1=argminvLα(fk+1,v,λk) forp=1 to n2 2.1. setc=Dpfk+1-λpk/β2-1/β 2.2. ifc≤0 then vpk+1=0 else select an initial guessw0 from Theorem 4 find a zerow of s(w) in (30) by Newton’s iterations vpk+1=wDpfk+1-λpk/βDpfk+1-λpk/β2 end if end for 3. updateλk+1=λk+β(Dfk+1-vk+1) 4. iferror<tol then exit 5. updateσ=σ×ratio;k=k+1 end while outputfk+1 and error


## 3.1. New Approach

The term∇fs0 in (14) can be approximated by a smooth function(15)sσ∇f≔n2-∑p=1n2exp⁡-∇fp22σ2=n2-∑p=1n2exp⁡-Dpf222σ2for a small positive parameter σ. So minimization (14) can be approximately obtained by solving(16)minf∇f1+αsσ∇fsubjecttoΦf-u22≤ε,smallσ>0.To deal with the constraint Φf-u22≤ε in (16), we rewrite (16) with a positive parameter μ as(17)minfμ2Φf-u22+∇f1+αsσ∇f,forsmallσ>0, or(18)minfμ2Φf-u22+∑p=1n2vp2+αsσvp2,forsmallσ>0,subjecttoDpf=vp,1≤p≤n2.Using a positive parameter β and Lagrange vectors λp∈R2, 1≤p≤n2, we define a Lagrangian function (19)Lαf,v,λ=μ2Φf-u22+∑p=1n2vp2-λpTDpf-vp+β2Dpf-vp22+αsσvp2.Finally, minimization (18) is converted to the following minimization:(20)minfLαf,v,λ.It is obvious that Lα(f,v,λ) is an extension of L(f,v,λ) in (8). Consequently the NADA is adopted to solve minimization (20). The kth iteration for solving (20) includes the following steps:(21)fk+1=argminfLαf,vk,λk,(22)vk+1=argminvLαfk+1,v,λk,(23)λk+1=argminλLαfk+1,vk+1,λ.It is noted that the objective functionLα(f,v,λ) in (19) has an extra term αsσ(vp), depending only on the vector v, compared with L(f,v,λ) in (8). So the methods for the calculation of fk+1 in (10) and λk+1 in (12) can be applied to the calculation of fk+1 in (21) and λk+1 in (23). For example, in (23), λk+1=λk+β(Dfk+1-vk+1). However, finding vk+1 in (22) is a challenging problem.


## 3.2. Calculation ofvpk+1

In this subsection we address how to calculatevpk+1 in (22) for 1≤p≤n2. Actually, for each p, we need to solve for vp from the following equation:(24)∂Lα∂vp=vpvp2+λp-βDpf+βvp+αvpσ2exp⁡-vp222σ2=0.It follows from (24) that(25)vp1+1βvp2+αβσ2exp⁡-vp222σ2=Dpf-λpβ. If vp2 is calculated then(26)vp=vp2Dpf-λp/βDpf-λp/β2.Thus, it suffices to determine vp2 from the following equation:(27)vp21+αβσ2exp⁡-vp222σ2=Dpf-λpβ2-1β.Denote c=Dpf-λp/β2-1/β. If c≤0 we set vp=0. If c>0, it follows from (27) that vp2 is a positive real solution of the equation(28)hz≔z+αzβσ2exp⁡-z22σ2-c=0.It is obvious that there is a positive zero z of h(z) between 0 and c and z→c as σ→0.Equation (28) can be rewritten as(29)exp⁡-z22σ2=βσ2c-zαzorz22σ2=ln⁡αzβσ2c-z. By substituting z=w, the above equation becomes(30)sw≔w-2σ2ln⁡αwβσ2c-w=0,on0,c2.Thus, we have shown the following lemma.Lemma 1.
Leth(z) and s(w) be given in (28) and (30), respectively. Then w is a zero of s(w) on (0,c2) if and only if z=w is a zero of h(z) on (0,c).Before we develop a procedure to calculate a positive root ofs(w)=0, we will list some properties of s(w) below.It is known that(31)s′w=1-σ2cwc-w=-ww+cw-σ2cwc-w. With z=w, the sign of s′(w) on (0,c2) is the same as that of the function(32)qz≔-z3+cz2-σ2con (0,c). By calculus, it is easy to see that q(z) is increasing on (0,2/3c) and decreasing on (2/3c,c) since q′(z)=-z(3z-2c). Thus, q(z) achieves a local maximum 4/27c3-σ2c at z=2/3c. Therefore, if 4/27c3-σ2c>0 (or σ<2/33c) then the equation q(z)=0 has two positive roots z1,z2 with z1<z2. It follows that q(z)<0 on (0,z1)∪(z2,c) and q(z)>0 on (z1,z2). On the other hand, if σ>2/33c, then q(z)<0 on (0,c). We have the following monotonicity of s(w).Property 2.
Ifσ<2/33c then s(w) is increasing on (z12,z22) and decreasing on (0,z12) and (z22,c2), where z1 and z2 are two positive zeros of q(z). If σ≥2/33c then s(w) is always decreasing on (0,c2).The second derivative(33)s′′w=σ2cw2c-w2c-32w indicates the following features of s(w).Property 3.
s ′ ′ ( w ) > 0 on (0,4/9c2) and s′′(w)<0 on (4/9c2,c2). The graph of s(w) has an inflection point (4/9c2, 4/9c2-2σ2ln⁡2α/βσ2). The function s(w) achieves a local minimum value at w=z12 and a local maximum value at w=z22. In addition, we have limw→0+s(w)=∞ and limw→c2-s(w)=-∞.Finally, we can determine the locations of zeros ofs(w) by convergent Newton’s iterations.Theorem 4.
Lets(w) and q(z) be, respectively, given in (30) and (32), where z=w. Let z2 be a larger positive root of q(z) on (0,c) for σ<2/33c. Then a positive root of s(w)=0 on (0,c2) can be calculated by Newton’s method. The convergence of Newton’s method is guaranteed if an initial guess w0 is selected as follows.
Case 1. (i) σ<2/33c, s(z22)>0 or (ii) σ≥2/33c, 2/9c2-σ2ln⁡2α/βσ2>0.
Choosew0<c2 and near c2 with s(w0)<0.
Case 2. (i) σ<2/33c, s(z22)<0 or (ii) σ≥2/33c, 2/9c2-σ2ln⁡2α/βσ2<0.
Choosew0=(βσ2c/α+βσ2)2.Proof.
We identify an interval containing a positive root ofs(w)=0 on which s′(w) and s′′(w) do not change signs for each subcase. Then the convergence of Newton’s method with an initial guess above is guaranteed.
Case 1. For subcase (i), there is a positive solution of s(w)=0 on (z22,c2) on which both s′(w) and s′′(w) are negative. For subcase (ii), the inflection point P of the graph of s(w) is on the upper half plane. So a positive zero is on (4/9c2,c2) on which both s′(w) and s′′(w) are negative. So one should choose w0<c2 and near c2 with s(w0)<0 for the convergence of Newton’s method.
Case 2. For subcase (i), there is a positive solution of s(w)=0 on (0,z12) on which s′(w)<0 and s′′(w)>0. For subcase (ii), the inflection point P is on the lower half plane. So a positive zero is on (0,4/9c2) on which s′(w)<0 and s′′(w)>0. One could choose w0=(βσ2c/α+βσ2)2 for the convergence of Newton iteration since βσ2c/α+βσ2 is a lower bound of positive solutions of h(z)=0 on (0,c) in (28).Thus,vpk+12, 1≤p≤n2, can be calculated from Lemma 1 and Theorem 4. Then vk+1 in (26) can be obtained.


## 3.3. Algorithm

Based on the above discussion, the calculation process of the proposedsl0+l1 regularization algorithm is summarized in Algorithm 1.Algorithm 1:Calculation process of sl0+l1 regularization algorithm (sl0-norm and l1-norm regularization algorithm).
inputΦ,u,α,ε initializeβ,μ,f0,v0,λ0,maxit,σ,δ,ratio,tol whilek<maxit 1. updatefk+1=argminfLα(f,vk,λk) by the NADA 1.1. decent directiondk=-∂/∂fLα(fk,vk,λk) 1.2. step sizesk such that Lα(fk+skdk,vk,λk)≤Ck-skδdk22 1.3. updatefk+1=fk+skdk 1.4. computeCk+1 such that Lα(fk+1,vk,λk)≤Ck+1≤Ck 2. updatevk+1=argminvLα(fk+1,v,λk) forp=1 to n2 2.1. setc=Dpfk+1-λpk/β2-1/β 2.2. ifc≤0 then vpk+1=0 else select an initial guessw0 from Theorem 4 find a zerow of s(w) in (30) by Newton’s iterations vpk+1=wDpfk+1-λpk/βDpfk+1-λpk/β2 end if end for 3. updateλk+1=λk+β(Dfk+1-vk+1) 4. iferror<tol then exit 5. updateσ=σ×ratio;k=k+1 end while outputfk+1 and error


## 4. Numerical Experiments and Results

In this section, the proposed combinedsl0-norm and l1-norm regularization is compared with the regularization without sl0-norm for its performance. Both regularization algorithms are implemented using the NADA. The MATLAB code is developed based on the software package TVAL3 [22], and the numerical experiments are conducted on an Intel Core i7 3.40 GHz PC. The 2D Shepp-Logan phantom and a cardiac image [25] of size 128×128 are tested. In each test, a random matrix Φ∈Rm×n2 (m≈0.3n2) is used to create the same system u=Φf+e for two regularization algorithms, where the noise e= 0.02∗mean(Φf)∗randn(m). The parameters are taken as α=1, β=23, μ=25 for the Shepp-Logan phantom and α=1, β=24, μ=27 for the cardiac image, respectively. The values of the parameter σ in the sl0-norm are decreasing at a ratio of 0.9 with a starting value σ0=0.1.Experiments are conducted to compare the reconstruction by the two algorithms after the same number of iterations. The original/reconstructed Shepp-Logan phantom images and the original/reconstructed cardiac images after same numbers of iterations are shown in Figure1. The reconstructed images show that the proposed sl0+l1-norm regularization produces better images while taking about 75% of the CPU time.Figure 1
Original and reconstructed images. First row: Shepp-Logan phantom after 150 iterations. Second row: cardiac image after 80 iterations.The quality of images is evaluated using the relative errorf-freconF/fF of the reconstructed image frecon in Frobenius norm. The root-mean-square error (RMSE), the normalized root-mean-square deviation (NRMSD), the normalized mean absolute deviation (NMAD), and the structural similarity index (SSIM) are also measured to reflect different aspects of the quality of the recovered images. RMSE evaluates the reconstruction quality on a pixel-by-pixel basis. NRMSD emphasizes large errors in a few pixels of the recovered image. NMAD focuses on small errors in the recovered image. SSIM compares the quality of the images using the original phantom image as a reference. A greater value of an SSIM indicates the better image quality. The experimental results from 100 tests are summarized in Table 1. The data in Table 1 indicates that the proposed sl0+l1-norm regularization provides a better accuracy after the same iteration number.Table 1
Comparison of two algorithms after same number of iterations.
(a) Shepp-Logan Phantom, after 70 iterations Time(s) Error RMSE NRMSD NMAD SSIM l 1 regularization 9.5 0.121 0.030 0.031 0.124 0.826 s l 0 + l 1 regularization 6.3 0.067 0.016 0.017 0.068 0.912  (b) Cardiac Image, after 50 iterations Time(s) Error RMSE NRMSD NMAD SSIM l 1 regularization 5.1 0.131 0.047 0.033 0.152 0.696 s l 0 + l 1 regularization 4.9 0.116 0.042 0.029 0.131 0.735Experiments are also conducted to compare iteration numbers and CPU time by the two algorithms when the same relative error is achieved and thus the image quality is same. The tolerance of the relative error to terminate the iteration is selected as 0.05. The average iteration numbers and CPU time from 100 tests of thesl0+l1-norm regularization and l1-norm only regularization are listed in Table 2. The comparison indicates that the CPU time is significantly reduced in the proposed algorithm while producing the same image quality. The graph of relative errors vs. the number of iterations for the two algorithms, shown in Figure 2, indicates that the sl0+l1-norm regularization requires less iterations to achieve the same accuracy.Table 2
Comparison of two algorithms with the same tolerance of relative error 0.05.
Phantom l 1 Regularization s l 0 + l 1 Regularization CPU Time Saved No. Iter. Time(s) No. Iter. Time(s) Shepp-Logan 129 17.4 86 7.4 58%  Cardiac 194 24.9 178 19.4 22%Figure 2
Relative error vs. iteration number in reconstruction by two algorithms to achieve the same tolerance of relative error 0.05. (a) Shepp-Logan phantom; (b) cardiac image.
(a) (b)Both types of the numerical experiments demonstrate that the proposedsl0+l1-norm regularization is superior to the regular l1-norm regularization minimization.


## 5. Discussions and Conclusions

The application of thesl0-norm in CT reconstruction has become one of active research topics recently. A combined sl0-norm and l1-norm regularization algorithm is proposed in this paper. The sl0-norm of a vector x is the limit of a smooth convex function sσ(x). In iterations of the proposed algorithm, the parameter σ is getting smaller and closer to zero. The smooth function sσ(∇f) in the objective function Lα(f,v,λ) approaches ∇f0, and the sparsity of the recovered image is fulfilled quickly. So the proposed algorithm improves the quality of the reconstructed images and the efficiency in terms of the iteration number and CPU time even with the extra computation from the added sl0-norm term. The numerical results demonstrate that the proposed algorithm improves the l1-norm regularization with the NADA in reducing the computation time significantly. The effects of the weight parameter α of the sl0-norm and the sl0-norm parameter σ are also tested. From our experiments, the parameter α from 0.5 to 2 and the initial value σ0 from 0.05 to 0.9 almost do not affect the efficiency of the algorithm. The decreasing ratio of parameter σ between 0.8 and 0.95 is a good choice, but a ratio smaller than 0.7 is not recommended.There are some aspects to be further investigated. It is a challenging problem to select good values of other parameters such asβ and μ in the Lagrangian function Lα(f,v,λ). The impact of these parameters on the performance of the algorithm will be further investigated.



---

*Source: 8398035-2019-06-02.xml*