Accelerated Chambolle Projection Algorithms for Image Restoration

: In this paper, the accelerated Chambolle projection algorithms based on Frank–Wolfe have been proposed. For solving the image restoration under the additive Gaussian noise, the Chambolle projection method (CP) is widely used. However, the projection operator has a large computational cost and complex form. By means of the Frank–Wolfe method, this projection operation can be greatly simpliﬁed. We propose two new algorithms, called Chambolle projection based on Frank–Wolfe (CP–FW) and Chambolle projection based on accelerated Frank–Wolfe (CP–AFW). They have a fast convergence rate and low computation cost. Furthermore, we extend the new algorithms to deal with the Poisson noise. The convergence of the new algorithms is discussed, and results of the experiment show their effectiveness and efﬁciency.

For solving the image restoration under the additive Gaussian noise, the Chambolle projection method (CP) is widely used.However, the projection operator has a large computational cost and complex form.By means of the Frank-Wolfe method, this projection operation can be greatly simplified.
In this paper, we propose new accelerated Chambolle projection algorithms for Gaussian denoising, according to the work of Chambolle.We call our new methods Chambolle projection based on Frank-Wolfe (CP-FW) and Chambolle projection based on accelerated Frank-Wolfe (CP-AFW).Meanwhile, we generalize the new algorithms under the Poisson noise condition.The convergence of algorithms is proved.Experiments show the better convergence rate of new algorithms and a low computational cost.
The outline of this paper is organized as follows.In Section 2, we show the related works and declare the notation.In Section 3, new algorithms are proposed under the Gaussian noise condition and a theoretical analysis is given.In Section 4, we extend these algorithms for Poisson noise removal.Section 5 shows the experiments of new algorithms.In Section 6, we summarize this paper.

Related Works
In this section, we first review the Chambolle projection algorithms for minimizing the total variation under the Gaussian noise.Then, the Frank-Wolfe method is briefly discussed.

Chambolle Projection Algorithm
For f = u + n, minimizing total variation [18] under the Gaussian noise has where f is the observed image, u is the clean image, n is the additive white Gaussian noise (AWGN), Ω is the image domain, and λ is a regularization parameter to keep tradeoff between data fidelity term and regularization term.
where ∂ is a sub-differential operator of J(u).Define w ∈ ∂J(u); then, (2) can be written as where J * is the conjugate function of J.Then, We obtain w = f −u λ is the minimizer of Since J * is the indicator function χ K (w), where K is the closure of the set where C 1 c Ω, R 2 represents the space of first-order continuously differentiable functions with compact support in domain Ω, R 2 .The divergence div is defined by div = −∇ * (∇ * is the adjoint of ∇).Thus, (5) becomes The objective function in (6) is a quadratic function, and the constraint set K is convex.We deduce the solution w = π K f λ , where π K is a projection operator, and the solution of ( 1) is shown by Next, we declare the notations in discrete setting.To simplify, we assume that the image f is given by a two-dimensional matrix of size N × N, although it would be easy to adapt all arguments for a more general M × N matrix to use the notations of chambolle, X = R N×N , and Y = X × X .If u ∈ X, the gradient ∇u ∈ Y is given by (∇u) i,j = (∂ x u) i,j , ∂ y u i,j with and and • 2 , it is defined by u 2 = u, u X and u, v X = ∑ i,j u(i, j)v(i, j).In the constraint set K, equivalently, for every p ∈ Y and u ∈ X, −divp, u X = p, ∇u Y .In Y, we use the Euclidean scalar product, defined in the standard way by p, q Y = ∑ 1 i,j N p 1 i,j q 1 i,j + p 2 i,j q 2 i,j for every p = (p 1 , p 2 ), q = (q 1 , q 2 ) ∈ Y.It is easy to check that, for every p = (p 1 , p 2 ) ∈ Y, discrete divergence div is given by Thus, the constraint set K in discrete setting is the closure of Chambolle suggests to solve the nonlinear projection π K of (6) by the following constraint minimization problem: and using the Lagrange multipliers and the fixed-point theory.Especially, let α i,j 0 be the Lagrange multipliers associated with each constraint in (9).We have for each i, j, − (∇(λdivp − f )) i,j + α i,j p i,j = 0, with either α i,j > 0 and p i,j = 1, or p i,j 1 and α i,j = 0. Thus, we have that in any case From this observation, applying the fixed-point algorithm (semi-implicit gradient descent), we compute p n+1 i,j as .
Thus, a sequence of iterations for solving (9) are written as where τ is a time step.The Chambolle projection algorithm is shown as in Algorithm 1.
It is worth noting that the Chambolle projection applying semi-implicit gradient descent method to calculate the nonlinear projection problem (6).It has a convergence rate O 1 k .

Frank-Wolfe Algorithm
We consider the following constraint optimization problem: where F is a smooth convex function and Θ is a non-empty convex and compact constraint set.The Frank-Wolfe method, also known as the conditional gradient method, to solve optimization problem (11), is shown as Algorithm 2.
Algorithm 2 Frank-Wolfe Algorithm (FW) Notice that the core idea is that the target function in Step 2 is a linear approximation of the original target function F(x), which is called the linear minimization oracle (LMO).Assuming that we apply the Frank-Wolfe method to (6), what will happen?By the linear minimization oracle (LMO), we will solve a linear approximation of the target function on the convex set K, which is easy by the special construction of the set K. The details are given in the next section.

Chambolle Projection Algorithms Based on Frank-Wolfe
, then DH = (w − f λ ), where D is the functional derivative operator.By the Algorithm 2, we have It is worth noting that ( 13) is a linear function and the set K in (8) has a special construction; thus, ( 13) is arg min Therefore, Obviously, y k = div p k .To sum up, the solver of ( 12) can be expressed as follows: Compared to (9), simple expressions and a lower computational cost are advantageous.Thus, we have the following Algorithm 3 about the Frank-Wolfe for improving the Chambolle projection to remove the Gaussian noise.It means that CP-FW and CP have the same convergence rate under the worst condition.However, for the same objective function H(w) and the constraint set K, CP-FW has simpler expressions than the expressions of Chambolle projection, and it can be further accelerated by the Nesterov momentum as follows in Algorithm 4, which is named as the Chambolle projection based on accelerated Frank-Wolfe (CP-AFW).

Proposition 3 ([21]
).For the optimization problem (12), the objective function H is convex and has a Lipschitz continuous gradient, and the constraint set K is convex and compact with diameter D. Futher, K is active, that is, DH(w * ) 2 G > 0, where G is constant.The constraint set K is transformed as L p (p 1) norm ball constraint set by (14).Especially, when p = 2, choosing r k = 2 k+3 and θ 0 = 0, AFW guarantees acceleration with the convergence rate where L is the Lipschitz continuous gradient constant, and C and T are the constants depending on L, D, and G. w * is the true solution of objective function H(w).
Proposition 3 indicates the AFW at least has O 1 k convergence rate.The Nesterov momentum is still helpful for accelerating the FW algorithm.
Next, we extend our new algorithms to the Poisson noise removal.
It is noticeable that the total variation along the lines of the ROF model can also be used to the Poisson denoising problem.The total variation model for Poisson noise removal [9] is below: where λ is a regularization parameter and ∇ is a gradient operator.We define J(u) = Ω |∇u| again.
Especially, (18) has the existence and uniqueness of solution.
By the Euler-Lagrange equation for (18), we have Please note that the derivation process of Chambolle projection under the Poisson noise condition is different from the Chambolle projection under the Gaussian noise condition.Equation ( 19) can be written as We let w = f −u λu ; then, we have u = f − λuw.Thus, f λu − w ∈ ∂J * (w) λu .Given by fixed u, we obtain the minimizing of the following optimization problem: For (21), it is a projection process that we attain w = π K f λu , and then we put Thus, the solution of ( 18) can be calculated as (22) is an implicit equation of u; then, we have by the iteration method.
Theorem 1.The sequence u k of (23 Proof.π λK is a firmly non-expansive operator.For Next, we discuss the computation of the projection operator π K .For each iteration, i.e, for fixed k, we adopt (24) to approximate π K ( f /λu k ).
When we denote f /u k as f k , Theorem 2 theoretically describes the convergence of the method.
Proof.The proof is similar to the Proposition 1. See [4].
Under the Poisson noise condition, the Chambolle projection algorithm is a two-level loop algorithm.For improving Algorithm 5, we replace the inner loop with the Frank-Wolfe.
Rethinking the projection process of ( 21), it is and then update u by u = f 1+λw .For solving (25) by Frank-Wolfe, let , then DG(w) = (w − f λu ).We have Thus, we obtained the Chambolle projection based on the Frank-Wolfe for Poisson noise removal in Algorithm 6.

Algorithm 6 Chambolle projection based on Frank-Wolfe for Poisson noise (CP-FW (Poisson))
for s = 1 : n do 3: y s = div p s .
11: end while 12: return u.Theoretical analysis indicates our simplified formulas are convergent.In the next section, numerical experiments show the effectiveness of our new algorithms.

Experiments
For the Chambolle projection and the accelerated Chambolle projection algorithms, we first show the change trend of sequences p k of Algorithm 1 and w k of Algorithms 3 and 4, or sequences p k of Algorithm 5 and w k of Algorithms 6 and 7 under the different noise conditions.Next, we test algorithms performance on the small dataset.

Change Rate of Sequence and Value of Energy Functional
Gaussian denoising: Firstly, we concentrate on the change rate of sequences p k of Algorithm 1 and w k of Algorithms 3 and 4 for Gaussian denoising.The "Papav" image is tested.The change rate of sequence is shown in Figure 1.
From Figure 1, we observe that the change rate of sequences p k of Algorithm 1 and w k of Algorithms 3 and 4: CP-AFW>CP-FW>CP, where horizontal axes represent iterations numbers and vertical axes represent the change rate of sequences (for a sequence of p k of Algorithm 1 and w k of Algorithms 3 and 4, relative error ε k+1 = , and then observe the change trend of sequences).In the experiments, optimal iterations of CP are above 100 iterations, CP-FW is about 60-80 iterations, and CP-AFW is about 20-30 iterations.Therefore, under the theoretical proof of algo-rithm convergence, we verify accelerated Chambolle projection convergent faster than the Chambolle projection algorithm ., and then observe the change trend of sequences).From the Figure 2b, the value of E(u) energy functional decreases with iterations and then remains almost unchanged, where horizontal axes represent outer-loop iterations numbers and vertical axes represent the value of energy functional.The value of the E(u) energy functional of CP-FW is the smallest, followed by CP-AFW, and finally by CP.Note that the tip: the value of energy functional of Equation ( 18) can be calculated by Meanwhile, by the experiments, the optimal inner loop iterations of CP are above 120 iterations, CP-FW is about 50-80 iterations, and CP-AFW is about 20-30 iterations.The outer loop iterations are about 3 iterations by the experiments.
Conclusion: Whether the Gaussian noise condition or Poisson, new accelerated Chambolle projection algorithms have simple expressions, a lower computational cost, and a faster convergence rate.In order to highlight a lower computation cost of the accelerated Chambolle projection, we test one iteration running time on a large image "Airplane", with a size of 512 × 512.The result can be seen in Figure 3a.
One iteration running time of CP-FW is the smallest.The running time of CP-AFW is similar to CP, but the formula of algorithms of CP-FW and CP-AFW is simpler than CP.Meanwhile, CP-AFW converges faster than CP.Thus, CP-FW and CP-AFW are superior to CP in the total computation cost.

Test Algorithms on Dataset
Next, we make a small image dataset for testing all of the algorithms for Gaussian denoising and Poisson denoising in image processing.The dataset includes 10 images, with a size of 256 × 256 pixels.We compare algorithms from the indexes of running time, PSNR and SSIM.
The small image dataset is shown as Figure 4. We apply algorithms for Gaussian denoising and Poisson denoising to the dataset.The experimental data can be seen in Figures 5 and 6.It is worth noting that the average PSNR of the dataset with the Gaussian noise is 22.2906 dB and the average SSIM is 0.3560.The average PSNR of the dataset with the Poisson noise is 26.4548 dB and the average SSIM is 0.5035.From Figures 5 and 6, the relationship of the running time of all algorithms is shown: CP > CP-FW > CP-AFW.That is, CP has the longest running time on the dataset.The relationship of the PSNR of the algorithms is shown: CP > CP-FW > CP-AFW, and the relationship of SSIM of algorithms is also shown: CP > CP-FW > CP-AFW.Meanwhile, we observe that the index of the PSNR value is similar, but SSIM is not.The SSIM value of CP-FW is similar to CP, but CP-AFW is smaller than the other algorithms, which means that the CP-AFW algorithm for solving the minimizing energy functional for denoising may lose image structural information.In conclusion, CP has a high calculation cost and an excellent denoising performance, but CP-AFW has a low calculation cost and slightly poorer denoising performance than CP.CP-FW is a good algorithm to recommend.We take an image "Papav" from the dataset to show the performance of Gaussian denoising with different algorithms in Figure 7 and the image "Peppers" to show the performance of Poisson denoising in Figure 8.The experiments show that the new accelerated Chambolle algorithms have better trade-off between the computational cost and denoising.Author Contributions: Conceptualization, W.W. and X.F.; methodology, W.W. and X.F.; software, .

11 :Algorithm 7
Output: The restoration clean image u.Meanwhile, the Chambolle projection algorithm based on the accelerated Frank-Wolfe for Poisson denoising can further accelerate the algorithm.See in Algorithm 7. Chambolle projection based on Accelerated Frank-Wolfe for Poisson noise (CP-AFW (Poisson))

13 :Theorem 3 .
Output: The restoration clean image u.CP (Poisson), CP-FW (Poisson), and CP-AFW (Poisson) are convergent.Proof.According to Theorems 1 and 2, the Chambolle projection for Poisson denoising converges.According to Theorem 1 and Proposition 2, the Chambolle projection based on the Frank-Wolfe algorithm for Poisson denoising converges.According to Theorem 1 and Proposition 3, the Chambolle projection based on the accelerated Frank-Wolfe for Poisson denoising converges.

Figure 1 .
Figure 1.The change rate of sequences under the Gaussian noise condition.Poisson denoising: Now, we consider the change rate of sequences p k of Algorithm 5 and w k of Algorithms 6 and 7 and the change of value of E(u) energy functional with iterations of three algorithms for Poisson denoising, where the "Papav" image is tested.It can be verified in Figure 2. From the Figure 2a, we can observe the same phenomenon of change trend like the Gaussian noise condition, where the horizontal axes represent inner loop iterations numbers and the vertical axes represent the change rate of sequences (for a sequence of p k of Algorithm 5 and w k of Algorithms 6 and 7, relative error ε k+1 = p k+1 −p k 2 p k 2

Figure 2 .
Figure 2. The change rate of sequences and the change of value of energy functional under the Poisson noise condition.(a) The change rate of sequences; (b) the change of the value of E(u) energy function with iterations.

Figure 3 .
Figure 3. Running time.(a) One iteration running time of three algorithms; (b) airplane.

Figure 5 .
Figure 5. Algorithms of comparation for Gaussian denoising.(a) Running time; (b) the average PSNR of dataset for Gaussian denoising; and (c) the average SSIM of dataset for Gaussian denoising.

Figure 6 .
Figure 6.Algorithms of comparation for Poisson denoising.(a) Running time; (b) the average PSNR of dataset for Poisson denoising; and (c) the average SSIM of dataset for Poisson denoising.