# A Preconditioning Technique for First-Order Primal-Dual Splitting Method in Convex Optimization

Authors: Meng Wen; Jigen Peng; Yuchao Tang; Chuanxi Zhu; Shigang Yue
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

We introduce a preconditioning technique for the first-order primal-dual splitting method. The primal-dual splitting method offers a very general framework for solving a large class of optimization problems arising in image processing. The key idea of the preconditioning technique is that the constant iterative parameters are updated self-adaptively in the iteration process. We also give a simple and easy way to choose the diagonal preconditioners while the convergence of the iterative algorithm is maintained. The efficiency of the proposed method is demonstrated on an image denoising problem. Numerical results show that the preconditioned iterative algorithm performs better than the original one.


---

## Body

## 1. Introduction

Many real-world application problems arising in signal and image processing [1–3], machine learning [4–6], and medical image reconstruction [7, 8] can be modeled as solving some convex optimization problems (maybe nonsmooth). In recent years, the optimization of the sum of two convex functions has received much attention, which takes the form of(1)minx∈Rnfx+hBx,where f∈Γ0(Rn), h∈Γ0(Rm), and B:Rn→Rm is a linear transformation matrix. Here Γ0(Rn) denotes the set of proper, lower semicontinuous convex functions from Rn to (-∞,+∞]. The functions f and g in (1) usually denote the data error term and the regularization term, respectively.Chambolle and Pock [9] proposed a general primal-dual method to solve problem (1). Under the assumption that the functions f and h have a closed-form solution of proximity operator, they proved the convergence of the proposed iterative algorithms. He and Yuan [10] have shown that the primal-dual method of Chambolle and Pock [9] is equivalent to the proximity point algorithm (PPA). Thus, the original convergence analysis can be easily obtained by the well-known PPA method. In order to accelerate the primal-dual method [9], Pock and Chambolle introduced in [11] a preconditioning technique for the primal-dual iterative algorithm for which the constant iterative parameters were replaced by some precondition iteration matrices. Then, the convergence of the preconditioned iterative algorithm followed directly from the PPA method. They also gave a practical way to choose the precondition iteration matrix. As an application, Sidky et al. [12] applied the primal-dual method of Chambolle and Pock [9] to solve a variety of application problems arising in medical image reconstruction. The primal-dual method and the corresponding precondition method presented the performance very well. Some related works can also be found in [13, 14].If the objective function in (1) has a differentiable term with Lipschitz continuous gradient, such as the least squares loss function, the primal-dual method introduced by Chambolle and Pock [9] did not use the gradient operator of the function. In order to solve such a more general problem, Combettes and Pesquet [15], Condat [1], and Vũ [16] introduced a primal-dual method. For instance, Condat [1] considered an optimization problem involving the sum of three convex functions, with a smooth function with Lipschitz continuous gradient, a nonsmooth proximable function, and linear composite functions. The problem is presented below:(2)minx∈Rnfx+gx+hBx,where f∈Γ0(Rn), g∈Γ0(Rn), h∈Γ0(Rm), f:Rn→R is differentiable, and its gradient ∇f is Lipschitz continuous with Lipschitz constant β>0; B:Rn→Rm is a linear transformation matrix. To solve problem (2), he proposed a primal-dual splitting method and proved the convergence of the new iterative algorithm in an infinite-dimensional Hilbert space based on the fixed point theory of nonexpansive mappings.The purpose of this paper is to introduce a preconditioned primal-dual splitting method to solve problem (2). The advantage of the preconditioning technique is that the iterative parameters will be updated self-adaptively. Furthermore, we give a family of preconditioners which are restricted to diagonal matrices and guarantee the convergence of the algorithm. To illustrate the efficiency of the proposed method, we compare it with the original method on image denoising problem. In addition, Combettes et al. [17] recently proposed a variable metric primal-dual method; they also referred to the preconditioning technique. But our algorithm is different from the method proposed by Combettes et al. In our algorithm, we use two different metrics: one for the primal variable; one for the dual variable(s). The other difference is in [17], which is a smooth term in the dual problem. When it is zero, the convergence conditions are stronger than the ones presented here. What is more, we explain how to choose the variable metric.The rest of the paper is organized as follows. In Section2, we briefly review the primal-dual splitting method proposed by Condat [1]. In Section 3, we present the preconditioning technique and provide a practical way to choose the iterative parameters matrices. In Section 4, we make several experiments in image denoising problems. Finally, we make a brief conclusion on this paper.


## 2. A Primal-Dual Splitting Method of Condat [1]

In this section, we briefly review the primal-dual splitting method introduced by Condat [1]. First, we introduce some definitions and notations. Let H be a real Hilbert space, with its inner product ·,· and norm ∥·∥1/2. We denote by Γ0(H) the set of proper, lower semicontinuous, convex functions from H to (-∞,+∞]. Let f∈Γ0(H); its Fenchel conjugate f∗∈Γ0(H) is defined by f∗(u)=supx〈x,u〉-f(x) and its proximity operator by proxλfu=argminx1/2x-u2+λf(x); λ>0 is a positive constant. We define the subdifferential of f as the set-valued operator ∂f(u)={v∈H:f(u′)≥f(u)+〈v,u′-u〉,∀u′∈H}. If f is differential at u, then ∂f(u)={∇f(u)}. Let H1 and H2 be two real Hilbert spaces and B:H1→H2 be a bounded linear operator with adjoint B∗ and induced norm (3)B=supx≤1Bx,x∈H1<+∞.Condat [1] considered solving the following optimization problem:(4)Findx^∈argminxfx+gx+hBx,where f:H1→R is convex, differentiable on H1, and its gradient ∇f is β-Lipschitz continuous, for some β>0; that is, ∇f(x)-∇f(y)≤βx-y, ∀x,y∈H1. g∈Γ0(H1) and h∈Γ0(H2), and their proximity operator have a closed-form solution. B:H1→H2 is a bounded linear operator with adjoint B∗. The dual formulation of the primal problem (4) is(5)Findy^∈argminyf+g∗-B∗y+h∗y.The corresponding saddle point of primal problem (4) and dual problem (5) is as follows:(6)Findx^,y^∈arg⁡minx⁡maxyfx+gx+Bx,y-h∗y.The pair (x^,y^) can be found via the following monotone variational inclusion:(7)00∈∂gx^+B∗y^+∇fx^-Bx^+∂h∗y^,where ∂g and ∂h∗ are the subdifferential of g and h∗, respectively.Condat [1] proposed the following iterative algorithm to solve problem (4).Algorithm 1 (primal-dual splitting method (PDS) for solving problem (4)).
Choose the proximal parametersτ>0,  σ>0 and the relaxation parameters {ρn}. Give an initial value (x0,y0)∈H1×H2. For n=0,1,2,…, the iterative {xn} and {yn} are updated by(1)
x~n+1=proxτgxn-τ(∇f(xn)+ef,n)-τB∗yn+eg,n;(2)
y~n+1=proxσh∗yn+σB(2x~n+1-xn)+eh,n;(3)
xn+1,yn+1=ρnx~n+1,y~n+1+1-ρnxn,yn,where the error terms ef,n∈H1, eg,n∈H1, and eh,n∈H2 model the inexact computation of the operators ∇f, proxτg, and proxσh, respectively.
If some stopping criteria have been reached, then the algorithm stops.The convergence of Algorithm1 was ensured by the following theorem.Theorem 2 (see [1]).
Letτ>0, σ>0, and the sequences {ρn}, {ef,n}, {eg,n}, and {eh,n} be the parameters of Algorithm 1. Let β be the Lipschitz constant of ∇f. Suppose that β>0 and the following conditions hold:(i)
1/τ-σB2>β/2;(ii)
{ρn}⊂(0,δ), where δ=2-β/21/τ-σB2-1;(iii)
∑ρn(δ-ρn)=+∞;(iv)
∑ρnef,n<+∞, ∑ρneg,n<+∞, and ∑ρneh,n<+∞. Then the sequences {xn} and {yn} generated by Algorithm 1 converge weakly to x^ of (4) and y^ of (5), respectively.Remark 3.
(1) Due to Moreau’s identityproxσh∗(u)=u-σprox1/σh(1/σu), the proximity operator proxσh∗ can be computed from prox1/σh.
(2) The role of primal variablex and dual variable y in Algorithm 1 can be exchanged. The convergence of the new iterative algorithm can also be ensured, according to the same Theorem 2. In practice, the performance of the two iterative algorithms is nearly the same.


## 3. A Preconditioned Primal-Dual Splitting Method for Solving (4)

In this section, we give a precondition version of Algorithm1. The main idea is motivated by the work of Pock and Chambolle [11]. The iterative parameters in Algorithm 1 are replaced by some symmetrical positive matrices. First, we give the detailed iterative algorithm below. Then we will prove its convergence.Algorithm 4 (preconditioned primal-dual splitting method (PPDS) for solving problem (4)).
Choose the symmetrical positive define matricesT, Σ and the relaxation parameters {ρn}. Give an initial value (x0,y0)∈H1×H2. For n=0,1,2,…, the iterative {xn} and {yn} are updated by (1)
x~n+1=proxTgxn-T(∇f(xn)+ef,n)-TB∗yn+eg,n;(2)
y~n+1=proxΣh∗yn+ΣB(2x~n+1-xn)+eh,n;(3)
(xn+1,yn+1)=ρn(x~n+1,y~n+1)+(1-ρn)(xn,yn),where the error terms ef,n∈H1, eg,n∈H1, and eh,n∈H2 model the inexact computation of the operators ∇f, proxTg, and proxΣh, respectively.
If some stopping criteria have been reached, then the algorithm stops.Theorem 5.
LetE be a symmetrical positive matrix and the sequences {ρn}, {ef,n}, {eg,n}, and {eh,n} be the parameters of Algorithm 4. Let β be the Lipschitz constant of ∇f. Suppose that β>0 and the following conditions hold:(i)
T-1-ΣB2>E/2, E≥β.(ii)
{ρn}⊂(0,δ), where δ=2-E/2T-1-ΣB2-1.(iii)
∑ρn(δ-ρn)=+∞.(iv)
∑ρnef,n<+∞, ∑ρneg,n<+∞, and ∑ρneh,n<+∞. Then the sequences {xn} and {yn} generated by Algorithm 4 converge weakly to x^ of (4) and y^ of (5), respectively.Proof.
The main idea of the proof is based on the method of Theorem2. First, we give some definitions and notations. Let Z=H1×H2 with inner product 〈·,·〉I: (8)z,z′I=x,x′+y,y′,∀z=x,y,z′=x′,y′∈Z.Then, Z is a real Hilbert space with the defined inner product.
Definezn=(xn,yn). It follows from the proof of Theorem 2; we know that the iterative sequences (1)–(3) in Algorithm 4 could be rewritten as follows: (9)zn+1=ρnI+P-1∘A-1zn-P-1∘Bzn-e2,n+e1,n+1-ρnzn,where P=T-1-B∗-BΣ-1, A=∂gB∗-B∂H∗, B=∇f0, e1,n=(eg,n,eh,n), and  e2,n=P-1(ef,n,-2Beg,n).
Notice that, from condition (i), it is easy to check thatP is bounded, self-adjoint, and strictly positive; that is, 〈z,Pz〉I>0, for every z≠0. Hence, we can define another inner product 〈·,·〉P and ·P=〈·,·〉P1/2 in Z as (10)z,z′P=z,Pz′I,for everyz,z′∈Z.Define M2=P-1∘B and k=T-1-ΣB2/E; we will prove that kM2 is 1-Lipschitz continuous. In fact, for any z,z′∈Z, we have(11)M2z-M2z′P2=P-1∘Bz-P-1∘Bz′,Bz-Bz′I=Σ-1T-1Σ-1-B∗B-1∇fx-∇fx′,∇fx∇fx′≤Σ-1T-1Σ-1-B2-1∇fx-∇fx′2≤β2Σ-1T-1Σ-1-B2-1x-x′2≤Ekx-x′2.Let us define a linear operator Q:(x,y)→((T-1-ΣB2)x,0). Then, P-Q is positive definite. So,(12)Ekx-x′2≤z-z′,Qz-z′I≤z-z′,Pz-z′=z-z′P2.Substituting (12) into (11), we have (13)kM2z-M2z′P≤z-z′P.The rest of the convergence proof follows the same argument of Theorem 2. So we omitted it here. This completes the proof.As a matter of fact,T and Σ in Theorem 5 could be any symmetric, positive definite maps. In order to ensure that the proximity operator of G and F∗ has a closed-form solution, it is sufficient to choose diagonal matrices for both of them. In the following, we give a practical way to choose the symmetric positive matrices and ensure the convergence of Theorem 5. To facilitate our proof, we need the following lemmas, which were obtained by [11].Lemma 6 (see [11]).
LetK be a well-defined matrix and T and Σ be symmetric positive definite maps satisfying Σ-1/2KT-1/22<1. Define the matrix M as (14)M≔T-KT-KΣ.Then the matrix M is symmetric and positive definite.Therefore, we are in the position to give our way to choose the matricesT and Σ accordingly.Lemma 7.
Assume thatE is a diagonal matrix with E=diag(β,β,…,β), where β is the Lipschitz constant of ∇f. Fix s∈[0,2]; let diagonal matrices T=diag(τ1,…,τn) and Σ=diag(σ1,…,σm) with (15)τj=1β/2+∑i=1mBi,j2-s,j=1,2,…,n,σi=1∑j=1nBi,js,i=1,2,…,m,and then it holds that(16)T-1-12E≻0,Σ≻0,(17)T-1-ΣB2>E2.Proof.
Conclusion (16) follows from the definition of diagonal matrices T and Σ. In fact, (18)1τj-12β=∑i=1mBi,j2-s>0,∀j=1,2,…,n,σi=1∑j=1nBi,js>0,∀i=1,2,…,m.We will prove (17). It is easy to see that the proof of (17) is equivalent to the positive definite map of the following matrix M: (19)M≔T-1-12E-BT-BΣ-1.From Lemma 6, it is sufficient to prove(20)Σ1/2BT-1-12E-1/2<1,which ensure that the matrix M is positive definite. For any j=1,2,…,n, let us define τj~=∑i=1m|Bi,j|2-s. We have(21)Σ1/2BT-1-12E-1/2x2=Σ1/2BT-1-12E-1/2x,Σ1/2BT-1-12E-1/2x=∑i=1mσi∑j=1nBi,jτj~-1/2xj2≤∑i=1mσi∑j=1nBi,js/2Bi,j1-s/2τj~-1/2xj2due to the Cauchy-Schwartz inequality≤∑i=1mσi∑j=1nBi,js∑j=1nBi,j2-sτj~-1xj2notice the condition ofσi=∑i=1m∑j=1nBi,j2-sτj~-1xj2=∑j=1n∑i=1mBi,j2-sτj~-1xj2from the definition ofτj~=x.By the definition of the operator norm, we obtain (22)Σ1/2BT-1-12E-1/22≤1.Finally, the strict inequality of (20) can be obtained from the above proof process with one of the above inequalities becoming strictly smaller.


## 4. Applications

In this section, we present an application of our proposed iterative algorithm. We aim at solving the following constrained total variation (TV) denoising problem:(23)minx∈C12x-u22+λxTV,where u∈Rn is a noisy image which was contaminated by Gaussian noise, λ>0 is the regularization parameter, and C is a closed convex set representing the prior information of the denoised image. By using the indicator function, constrained (TV) denoising problem (23) could be formulated by the following unconstrained optimization problem as follows:(24)minx12x-u22+λxTV+ιCx,where the indicator function ιC∈Γ0(Rn):x∈Rn↦{0,ifx∈C;+∞,otherwise}, since the total variation term xTV can be represented by a combination of convex function φ and linear transformation matrix B; that is, xTV=φ(Bx). See, for example, [18]. Then, optimization problem (24) is actually a special case of the general optimization problem (2) with f(x)=1/2x-u22, g(x)=ιC(x), and h(x)=λφ(x). Note that the gradient of the function f(x) is ∇f(x)=x-u and the Lipschitz constant β=1.If the constrained setC=Rn, then the constrained TV denoising problem reduces to the unconstrained TV denoising problem:(25)minx12x-u22+λxTV.The above TV denoising problem is often referred to as the ROF denoising model which was first introduced in computer vision by Rudin et al. [19].


### 4.1. Numerical Experiments

In the following, we present some preliminary numerical results and show the efficiency of our proposed methods. All the experiments are run on a personal Lenovo computer with Pentium(R) Dual-Core CPU @ 2.8 GHz and RAM 4 GB.For all the tested iterative algorithms the stopping criterion is(26)xk+1-xkxk≤ϵ,where ϵ>0 is a given small constant, or the maximum iteration numbers 4000 are reached. The reconstructed image is evaluated in terms of the signal to noise (SNR) defined by(27)SNR=10log10xori2xrec-xori2,where xori is the original clear image and xrec is the reconstructed image. The reconstructed time is denoted by “T(s)” and the iteration number “Iter” is recorded when the stopping criteria satisfied. The tested images are the well-known “Barbara,” “Boat,” and “Lena”. These images have the same size of 512×512 and are displayed in Figure 1.Figure 1
Original images. (a) “Barbara”; (b) “Boat”; (c) “Lena.”
(a) (b) (c)Experiment 1.
We present how the iterative parameters are chosen. First, we choose different combination of parametersσ and τ and then apply Algorithm 1 to solve image denoising problem (23). We choose “Barbara” as the test image and add it by random Gaussian noise with zero mean and standard deviation 0.05. For convenience, we define τmax=2/β, σmax=1/B21/τ-β/2. Here, the Lipschitz constant β=1 and B2=8. The numerical results are reported in Table 1.
Meanwhile, we plot the SNR versus the iteration numbers in Figure2. We can see from Figure 2 that when the iterative parameter τ=0.9τmax and τ=0.5τmax, the greater σ, the faster the convergence speed. For the small τ, it has no apparent difference by the choice of parameter σ. So we select τ=0.1τmax and σ=0.1σmax for Algorithm 1 for the comparative experiments.
For our proposed Algorithm4, the corresponding preconditioned iterative matrices are chosen according to Lemma 7. There is only one parameter s that needs to be set, and the experiments results are reported in Table 2. Figure 3 shows the SNR value versus the iteration number. We can see from it that the performance of the three choices of s is nearly the same. So we choose s=1 for Algorithm 4 in the following test.Table 1
Numerical results obtained by Algorithm1.
Constraint setC Parameterτ= Parameterσ= ϵ = 1 0 - 3 T   ( s ) I t e r S N R x ∣ x ≥ 0 0.9 τ m a x σ m a x 2.17 59 23.56 0.5 ∗ σ m a x 1.29 37 22.76 0.1 ∗ σ m a x 1.22 35 20.93 0.5 τ m a x σ m a x 0.46 13 23.61 0.5 ∗ σ m a x 0.60 17 23.59 0.1 ∗ σ m a x 1.02 29 23.13 0.1 τ m a x σ m a x 0.88 25 23.57 0.5 ∗ σ m a x 0.88 25 23.58 0.1 ∗ σ m a x 0.92 26 23.59  x∣0≤x≤1 0.9 τ m a x σ m a x 2.17 59 23.56 0.5 ∗ σ m a x 1.33 37 22.76 0.1 ∗ σ m a x 1.21 34 20.91 0.5 τ m a x σ m a x 0.48 13 23.62 0.5 ∗ σ m a x 0.62 17 23.60 0.1 ∗ σ m a x 1.07 29 23.13 0.1 τ m a x σ m a x 0.91 25 23.57 0.5 ∗ σ m a x 0.91 25 23.58 0.1 ∗ σ m a x 0.94 26 23.59Table 2
Numerical results obtained by Algorithm4.
Constraint setC Parameters= ϵ = 1 0 - 3 T (s) Iter SNR x ∣ x ≥ 0 0.5 0.80 23 23.58 1 0.81 23 23.58 1.5 0.81 23 23.58  x ∣ 0 ≤ x ≤ 1 0.5 0.90 23 23.58 1 0.82 23 23.58 1.5 0.82 23 23.58Figure 2
(a) The relation between SNR and iteration number with the nonnegative constraint, whereτ=0.9τmax in (a1), τ=0.5τmax in (a2), and τ=0.1τmax in (a3). (b) The relation between SNR and iteration number with the bounded constraint, where τ=0.9τmax in (b1), τ=0.5τmax in (b2), and τ=0.1τmax in (b3).
(a) (b)Figure 3
(a) The relation between SNR and iteration number with the nonnegative constraint for different choice ofs. (b) The relation between SNR and iteration number with the bounded constraint for different choice of s.
(a) (b)Experiment 2.
We show the performance comparison between Algorithm1 and our proposed Algorithm 4. The constraint set C is set as nonnegative set; that is, C=x∣x≥0. To perform fair comparison, we add each of these images with random Gaussian noise with 0 mean value and different level of standard variation σ. The regularization parameter λ=0.02 when the noise level σ≤0.05 and λ=0.05 when σ>0.05.
From Table3, we can see that our proposed Algorithm 4 converges faster than Algorithm 1 in terms of iteration numbers and iteration time in CPU time. For the large noise level, Algorithm 4 reaches higher SNR value than Algorithm 1 with less iteration numbers. Both of the algorithms get the cleared image, which have nearly the same SNR value finally. To visualize the reconstructed images, we present the noised image and the final denoising image in Figures 4, 5, and 6, respectively.Table 3
Performance comparison between Algorithms1 and 4.
Image Noise level Methods ϵ = 10 - 3 ϵ = 10 - 4 T   ( s ) Iter SNR T   ( s ) I t e r S N R “Barbara” 0.01 Algorithm1 0.87 26 26.06 2.18 64 26.07 Algorithm4 0.80 23 26.03 1.17 34 26.08 0.05 Algorithm1 0.88 26 23.59 1.86 55 23.60 Algorithm4 0.80 23 23.58 1.12 32 23.60 0.1 Algorithm1 1.37 32 19.67 3.20 94 19.71 Algorithm4 0.81 23 19.70 1.26 36 19.71  “Boat” 0.01 Algorithm1 0.88 26 28.63 2.06 62 28.62 Algorithm4 0.81 23 28.56 1.18 34 28.63 0.05 Algorithm1 1.02 26 25.70 1.82 54 25.72 Algorithm4 0.89 23 25.70 1.12 32 25.72 0.1 Algorithm1 1.03 31 22.39 3.05 92 22.45 Algorithm4 0.79 23 22.44 1.23 36 22.45  “Lena” 0.01 Algorithm1 0.93 26 30.55 2.19 65 30.48 Algorithm4 0.79 23 30.39 1.21 35 30.48 0.05 Algorithm1 0.87 26 26.47 1.92 57 26.51 Algorithm4 0.79 23 26.48 1.09 32 26.50 0.1 Algorithm1 1.05 31 23.21 3.23 96 23.34 Algorithm4 0.80 23 23.33 1.26 37 23.34Figure 4
(a1)–(a3) are the noisy image of “Barbara” with noise levels0.01, 0.05, and 0.1, respectively; (b1)–(b3) are the denoised image by Algorithm 1; (c1)–(c3) are the denoised image by our proposed Algorithm 4.Figure 5
(a1)–(a3) are the noisy image of “Boat” with noise levels0.01, 0.05, and 0.1, respectively; (b1)–(b3) are the denoised image by Algorithm 1; (c1)–(c3) are the denoised image by our proposed Algorithm 4.Figure 6
(a1)–(a3) are the noisy image of “Lena” with noise levels0.01, 0.05, and 0.1, respectively; (b1)–(b3) are the denoised image by Algorithm 1; (c1)–(c3) are the denoised image by our proposed Algorithm 4.


## 4.1. Numerical Experiments

In the following, we present some preliminary numerical results and show the efficiency of our proposed methods. All the experiments are run on a personal Lenovo computer with Pentium(R) Dual-Core CPU @ 2.8 GHz and RAM 4 GB.For all the tested iterative algorithms the stopping criterion is(26)xk+1-xkxk≤ϵ,where ϵ>0 is a given small constant, or the maximum iteration numbers 4000 are reached. The reconstructed image is evaluated in terms of the signal to noise (SNR) defined by(27)SNR=10log10xori2xrec-xori2,where xori is the original clear image and xrec is the reconstructed image. The reconstructed time is denoted by “T(s)” and the iteration number “Iter” is recorded when the stopping criteria satisfied. The tested images are the well-known “Barbara,” “Boat,” and “Lena”. These images have the same size of 512×512 and are displayed in Figure 1.Figure 1
Original images. (a) “Barbara”; (b) “Boat”; (c) “Lena.”
(a) (b) (c)Experiment 1.
We present how the iterative parameters are chosen. First, we choose different combination of parametersσ and τ and then apply Algorithm 1 to solve image denoising problem (23). We choose “Barbara” as the test image and add it by random Gaussian noise with zero mean and standard deviation 0.05. For convenience, we define τmax=2/β, σmax=1/B21/τ-β/2. Here, the Lipschitz constant β=1 and B2=8. The numerical results are reported in Table 1.
Meanwhile, we plot the SNR versus the iteration numbers in Figure2. We can see from Figure 2 that when the iterative parameter τ=0.9τmax and τ=0.5τmax, the greater σ, the faster the convergence speed. For the small τ, it has no apparent difference by the choice of parameter σ. So we select τ=0.1τmax and σ=0.1σmax for Algorithm 1 for the comparative experiments.
For our proposed Algorithm4, the corresponding preconditioned iterative matrices are chosen according to Lemma 7. There is only one parameter s that needs to be set, and the experiments results are reported in Table 2. Figure 3 shows the SNR value versus the iteration number. We can see from it that the performance of the three choices of s is nearly the same. So we choose s=1 for Algorithm 4 in the following test.Table 1
Numerical results obtained by Algorithm1.
Constraint setC Parameterτ= Parameterσ= ϵ = 1 0 - 3 T   ( s ) I t e r S N R x ∣ x ≥ 0 0.9 τ m a x σ m a x 2.17 59 23.56 0.5 ∗ σ m a x 1.29 37 22.76 0.1 ∗ σ m a x 1.22 35 20.93 0.5 τ m a x σ m a x 0.46 13 23.61 0.5 ∗ σ m a x 0.60 17 23.59 0.1 ∗ σ m a x 1.02 29 23.13 0.1 τ m a x σ m a x 0.88 25 23.57 0.5 ∗ σ m a x 0.88 25 23.58 0.1 ∗ σ m a x 0.92 26 23.59  x∣0≤x≤1 0.9 τ m a x σ m a x 2.17 59 23.56 0.5 ∗ σ m a x 1.33 37 22.76 0.1 ∗ σ m a x 1.21 34 20.91 0.5 τ m a x σ m a x 0.48 13 23.62 0.5 ∗ σ m a x 0.62 17 23.60 0.1 ∗ σ m a x 1.07 29 23.13 0.1 τ m a x σ m a x 0.91 25 23.57 0.5 ∗ σ m a x 0.91 25 23.58 0.1 ∗ σ m a x 0.94 26 23.59Table 2
Numerical results obtained by Algorithm4.
Constraint setC Parameters= ϵ = 1 0 - 3 T (s) Iter SNR x ∣ x ≥ 0 0.5 0.80 23 23.58 1 0.81 23 23.58 1.5 0.81 23 23.58  x ∣ 0 ≤ x ≤ 1 0.5 0.90 23 23.58 1 0.82 23 23.58 1.5 0.82 23 23.58Figure 2
(a) The relation between SNR and iteration number with the nonnegative constraint, whereτ=0.9τmax in (a1), τ=0.5τmax in (a2), and τ=0.1τmax in (a3). (b) The relation between SNR and iteration number with the bounded constraint, where τ=0.9τmax in (b1), τ=0.5τmax in (b2), and τ=0.1τmax in (b3).
(a) (b)Figure 3
(a) The relation between SNR and iteration number with the nonnegative constraint for different choice ofs. (b) The relation between SNR and iteration number with the bounded constraint for different choice of s.
(a) (b)Experiment 2.
We show the performance comparison between Algorithm1 and our proposed Algorithm 4. The constraint set C is set as nonnegative set; that is, C=x∣x≥0. To perform fair comparison, we add each of these images with random Gaussian noise with 0 mean value and different level of standard variation σ. The regularization parameter λ=0.02 when the noise level σ≤0.05 and λ=0.05 when σ>0.05.
From Table3, we can see that our proposed Algorithm 4 converges faster than Algorithm 1 in terms of iteration numbers and iteration time in CPU time. For the large noise level, Algorithm 4 reaches higher SNR value than Algorithm 1 with less iteration numbers. Both of the algorithms get the cleared image, which have nearly the same SNR value finally. To visualize the reconstructed images, we present the noised image and the final denoising image in Figures 4, 5, and 6, respectively.Table 3
Performance comparison between Algorithms1 and 4.
Image Noise level Methods ϵ = 10 - 3 ϵ = 10 - 4 T   ( s ) Iter SNR T   ( s ) I t e r S N R “Barbara” 0.01 Algorithm1 0.87 26 26.06 2.18 64 26.07 Algorithm4 0.80 23 26.03 1.17 34 26.08 0.05 Algorithm1 0.88 26 23.59 1.86 55 23.60 Algorithm4 0.80 23 23.58 1.12 32 23.60 0.1 Algorithm1 1.37 32 19.67 3.20 94 19.71 Algorithm4 0.81 23 19.70 1.26 36 19.71  “Boat” 0.01 Algorithm1 0.88 26 28.63 2.06 62 28.62 Algorithm4 0.81 23 28.56 1.18 34 28.63 0.05 Algorithm1 1.02 26 25.70 1.82 54 25.72 Algorithm4 0.89 23 25.70 1.12 32 25.72 0.1 Algorithm1 1.03 31 22.39 3.05 92 22.45 Algorithm4 0.79 23 22.44 1.23 36 22.45  “Lena” 0.01 Algorithm1 0.93 26 30.55 2.19 65 30.48 Algorithm4 0.79 23 30.39 1.21 35 30.48 0.05 Algorithm1 0.87 26 26.47 1.92 57 26.51 Algorithm4 0.79 23 26.48 1.09 32 26.50 0.1 Algorithm1 1.05 31 23.21 3.23 96 23.34 Algorithm4 0.80 23 23.33 1.26 37 23.34Figure 4
(a1)–(a3) are the noisy image of “Barbara” with noise levels0.01, 0.05, and 0.1, respectively; (b1)–(b3) are the denoised image by Algorithm 1; (c1)–(c3) are the denoised image by our proposed Algorithm 4.Figure 5
(a1)–(a3) are the noisy image of “Boat” with noise levels0.01, 0.05, and 0.1, respectively; (b1)–(b3) are the denoised image by Algorithm 1; (c1)–(c3) are the denoised image by our proposed Algorithm 4.Figure 6
(a1)–(a3) are the noisy image of “Lena” with noise levels0.01, 0.05, and 0.1, respectively; (b1)–(b3) are the denoised image by Algorithm 1; (c1)–(c3) are the denoised image by our proposed Algorithm 4.


## 5. Conclusion

In this paper, we have studied the general optimization problem with the sum of three convex functions which is composed of a differential function with Lipschitz continuous gradient, a proximable function, and a linear composition function. Many interesting problems arising in image restoration and image reconstruction are special case of this problem. Inspired by the preconditioning technique proposed by Pock and Chambolle, we have introduced a primal-dual splitting algorithm with self-adaptive step-size to solve such problem. We also proposed a practical way to choose these step-sizes with a proof of convergence. Numerical results on image denoising problem showed that the precondition iterative algorithm performs better than the original one with constant step-size.



---

*Source: 3694525-2017-01-03.xml*