Image Restoration Using a Fixed-Point Method for a TVL2 Regularization Problem

: In this paper, we ﬁrst propose a new TVL2 regularization model for image restoration, and then we propose two iterative methods, which are ﬁxed-point and ﬁxed-point-like methods, using CGLS (Conjugate Gradient Least Squares method) for solving the new proposed TVL2 problem. We also provide convergence analysis for the ﬁxed-point method. Lastly, numerical experiments for several test problems are provided to evaluate the effectiveness of the proposed two iterative methods. Numerical results show that the new proposed TVL2 model is preferred over an existing TVL2 model and the proposed ﬁxed-point-like method is well suited for the new TVL2 model.


Introduction
Image restoration is the fundamental problem in image processing that recovers a true image from a blurry and noisy image. The problem of image restoration usually reduces to find the optimal solution u ∈ R m based on the following model: where A ∈ R m×m is a blurring operator, ε ∈ R m is an unknown white Gaussian noise with variance σ, and f and u denote the observed degraded image and the original image, respectively. Our purpose is to restore the original image u from blurred and noisy image f as well as possible.
Over the past few decades, optimization techniques and various variation models [1][2][3] have been widely studied and applied in many image processing fields. The well-known ROF (Rudin-Osher-Fatemi) total variation model [3] produces the deblurred image given by the following minimization problem: where TV(u) is the total variation (TV) of u and β > 0 is a regularization parameter. Many computational methods for solving problem (2) have been proposed in recent years. For example, the time-marching PDE method, the subgradient descent method, the Newton-like method, the second-order cone programming method, the lagged diffusivity fixed-point method, and the split Bregman method have been proposed by many researchers (see [1,[3][4][5][6][7][8][9][10][11][12][13] for details).
Recently, Chen et al. [14] proposed a fixed-point method for solving the following constrained TVL2 deblurring problem: where α and β are positive constants, and C is a closed convex subset of R n 2 . In this paper, we only consider the case for C = R n 2 . That is, we only consider the unconstrained TVL2 problem (3) with C = R n 2 . This approach motivates us to propose a new TVL2 deblurring model min u∈R n 2 1 2 where A ∈ R n 2 ×n 2 is a blurring matrix, u ∈ R n 2 is an original image, f ∈ R n 2 is a degraded image, α and β are positive regularization parameters, and TV(u) stands for the isotropic TV of u. Note that the new TVL2 model (4) uses a non-smooth term u 2 instead of using a smooth term u 2 2 for the purpose of better preserving the edges and corners in the restored images. The isotropic TV of u is defined by where the discrete gradient operator ∇ : R n 2 → R 2n 2 is defined by Notice that TVL2 problems (3) have a unique solution since its objective function is strictly convex, while the new TVL2 problem (4) may not have a unique solution since its objective function is just convex, not strictly convex.
The purpose of this paper is to propose two iterative methods, which are fixed-point and fixed-point-like methods, using CGLS (Conjugate Gradient Least Squares method [15]) for solving the new proposed TVL2 problem (4). This paper is organized as follows. In Section 2, we introduce some definitions and properties that are used in this paper. In Section 3, we first propose a fixed-point method using CGLS for solving TVL2 problem (4), and then we provide convergence analysis for the fixed-point method. In Section 4, we propose a fixed-point-like method using CGLS for solving TVL2 problem (4). In Section 5, we just provide the split Bregman methods for solving TVL2 problems (3) and (4) to see how efficiently the fixed-point and fixed-point-like methods perform. In Section 6, we describe how to carry out numerical experiments for several test problems in order to evaluate the effectiveness of the proposed two iterative methods. This can be done by comparing their performances for TVL2 problem (4) with those of the fixed-point method proposed in [14] for TVL2 problem (3) and the split Bregman methods for TVL2 problems (3) and (4). In Section 7, we provide numerical results for all test problems. Lastly, some conclusions are drawn.

Preliminaries
In this section, we briefly refer to some definitions and important properties which will be the foundation for the development of algorithms proposed in this paper. Definition 1. Let ϕ : R n → R ∪ {+∞} be a proper, convex, and lower semi-continuous function. The proximity operator of ϕ at x ∈ R n is defined by [16]).
Definition 2. Let ϕ : R n → R ∪ {+∞} be a proper, convex, and lower semi-continuous function. The subdifferential of ϕ at x ∈ R n is defined by Definition 3. A nonlinear operator H : R n → R n is called non-expansive if for any x, y ∈ R n , Definition 4. A nonlinear operator H : R n → R n is called firmly non-expansive if for any x, y ∈ R n , It is easy to show that a firmly non-expansive operator is non-expansive. The following proposition shows a relationship between the proximity operator and the subdifferential of a convex function.
Let ϕ : R 2n 2 → R be a convex function defined by and let B be a 2n 2 × n 2 matrix that represents a discrete gradient operator ∇, which is where I n is the n × n identity matrix, ⊗ denotes the Kronecker product, and D is the first order finite difference matrix of order n Then, the isotropic TV of u ∈ R n 2 can be expressed as We now provide the fixed-point method, called Algorithm 1, for solving TVL2 problem (3), which was proposed by Chen et al. [14].  Notice that line 4 of Algorithm 1 is solved using CGLS instead of using CG (Conjugate Gradient method [19][20][21]) since the linear system in line 4 is equivalent to solving the following least squares problem For all algorithms presented in this paper, maxit denotes the maximum number of iterations, and tol denotes the tolerance value of the stopping criterion.

Fixed-Point Method for TVL2 Problem (4)
In this section, we propose a fixed-point method using CGLS for solving the new TVL2 regularization problem (4). From relation (12), TVL2 problem (4) can be expressed as min u∈R n 2 1 2 Using Proposition 1, we can obtain the following property for a solution of TVL2 problem (13).

Theorem 1.
If ϕ is a real valued convex function on R 2n 2 , B is a 2n 2 × n 2 matrix, A is an n 2 × n 2 matrix, u is a solution to model (13), then, for any γ, λ > 0 there exist vectors a ∈ R n 2 and b ∈ R 2n 2 such that a = (I − prox 1 Conversely, if there exists γ, λ > 0, a ∈ R n 2 , b ∈ R 2n 2 and u ∈ R n 2 satisfying Equations (14)- (16), then u is a solution to model (13).
From Equations (14)- (16) in Theorem 1, we can develop a fixed-point algorithm which converges to a solution to TVL2 problem (4). We now describe how to develop the fixed-point algorithm. Let u be an approximate solution to the ill-conditioned linear system (16) in Theorem 1. Then, u can be expressed as where M is a symmetric positive semi-definite matrix approximating an inverse of the linear system matrix in (16). For example, we can choose M = (A T A) † r , which is a truncated pseudoinverse of A T A using the r largest positive singular values of A T A. Substituting Equation (20) into Equations (14) and (15), one obtains Let us define some operators. For the given convex functions 1 γ · 2 on R n 2 and β λ ϕ on R 2n 2 , we define an operator T : R 3n 2 → R 3n 2 at a vector ( u v ) ∈ R 3n 2 with u ∈ R n 2 and v ∈ R 2n 2 as follows: We also introduce an affine transformation L : R 3n 2 → R 3n 2 defined, for all ( a b ) ∈ R 3n 2 with a ∈ R n 2 and b ∈ R 2n 2 , by and an operator G : R 3n 2 → R 3n 2 defined by Then, Equations (21) and (22) can be expressed as where w = ( a b ).

Proposition 2. The operator G defined by (25) has a fixed point.
Proof. Since a solution of TVL2 problem (13) exists, from Equation (26) and the first part of the proof of Theorem 1, G has a fixed point.

Lemma 1.
If the operator T is defined by Equation (23), then the operator T is non-expansive.
Hence, one obtains T(s) − T(t) 2 ≤ s − t 2 , which means that the operator T is non-expansive. Let Then, Label (24) can be expressed as Proposition 3. If ϕ is a convex function on R 2n 2 , B is a 2n 2 × n 2 matrix, A is an n 2 × n 2 matrix and α, γ, λ are positive constants such that I 3n 2 − PMP T R 2 ≤ 1, then G is non-expansive.
Proof. Since the operator T is non-expansive by Lemma 1, for all w 1 , w 2 ∈ R 3n 2 , we have By the assumption I 3n 2 − PMP T R 2 ≤ 1 and Equation (27), we obtain Hence, G is non-expansive.
Let S : R N → R N be an operator. Then, the Picard iteration of the operator S is defined by for a given vector x 0 ∈ R N . For κ ∈ (0, 1), the κ−averaged operator S κ of S is defined by Proposition 4 (Optial κ-averaged Theorem). Let C be a closed convex set in R 3n 2 and let S : C → C be a non-expansive mapping with at least one fixed point. Then, for any w 0 ∈ C and κ ∈ (0, 1), the Picard iteration of S κ converges to a fixed point of S (see [22]).

Theorem 2.
If ϕ is a convex function, A is an n 2 × n 2 matrix, B is a 2n 2 × n 2 matrix and α, γ, λ > 0 are positive constants such that I 3n 2 − PMP T R 2 ≤ 1, then for any κ ∈ (0, 1) the Picard iteration of G κ converges to a fixed point of G.
Proof. From Propositions 2 and 3, we know that the opertor G has a fixed point and is non-expansive. Hence, for any w 0 ∈ R 3n 2 and κ ∈ (0, 1), the Picard iteration of G κ converges to a fixed point of G by Proposition 4.
For a square matrix K, let ρ(K) denote the spectral radius of K. Then, the following lemma can be obtained. (10) and (11) respectively, and let A be a given n 2 × n 2 blurring matrix. If we choose α, γ and λ such that

Lemma 2. Let ϕ and B be defined by Equations
then I 3n 2 − PMP T R 2 ≤ 1 and thus the operator G is non-expansive.

Theorem 3.
If the assumptions of Lemma 2 hold and κ ∈ (0, 1), then the Picard iteration of G k converges to a fixed point of G.
Proof. The proof follows from Lemma 2 and Theorem 2.
From Theorem 2 and the Picard iteration of the κ-averaged operator G κ = κI + (1 − κ)G, we can obtain a fixed-point method, called Algorithm 2, which converges to a solution to TVL2 problem (4).

Algorithm 2
Fixed-point method for TVL2 problem (4) 1: Given : observed image f , positive parameters α, β, γ, λ and κ ∈ (0, 1) 2: Initialization : a 0 = 0, b 0 = 0 and u 0 = f 3: for k = 0 to maxit do 4:â k+1 = I − prox 1 γ · 2 u k + a k 5: < tol then 10: Stop 11: end if 12: end for The linear system in line 8 of Algorithm 2 is ill-conditioned, so we need to consider how to find an approximate solution to the ill-conditioned linear system. A typical method for finding an approximate solution to the linear system is However, computation of (A T A) † r is very time-consuming when A is large. Thus, we want to propose a different approach for finding an approximate solution to the linear system in line 8 of Algorithm 2. We first split the coefficient matrix A T A into where δ is a positive constant and D = diag(A T A) is a diagonal part of A T A. Then, the ill-conditioned linear system can be solved using the following iterative method:

Inner Solver
Choose y 1 = u k for = 1 to maxl Solve (A T A + δD)y +1 = δDy + A T f − αγa k+1 − λB T b k+1 for y +1 end for u k+1 = y +1 , where u k refers to the restored image computed at the previous kth step, and the optimal parameter δ > 0 is chosen by numerical tries. Semi-convergence analysis for Inner Solver has been studied by Han and Yun [23]. Algorithm 2 is an iterative method that converges to a solution to TVL2 problem (4) as k → ∞, so we do not have to get an accurate solution to the linear system in line 8 of Algorithm 2. For this reason, we have used maxl = 1 for Inner Solver.
Notice that the linear system in Inner Solver is equivalent to solving the following least squares problem: Hence, the linear system in Inner Solver is solved by applying the CGLS to (29).

Fixed-Point-like Method for TVL2 Problem (4)
In this section, we propose a fixed-point-like method using CGLS for solving the new TVL2 problem (4) that can be obtained by modifying Algorithm 2. Notice that Algorithm 2 computesâ k+1 and a k+1 before the solution step of finding u k+1 (see lines 4 and 5 of Algorithm 2). However, the fixed-point-like method to be proposed in this section computesâ k+1 and a k+1 after the solution step of finding u k+1 . Below, we describe how to develop the fixed-point-like method in detail. We first split line 4 of Algorithm 2 intoâ where a k+ 1 2 = a k − prox 1 γ · 2 u k + a k . Replacing the old value u k of Equation (30) with the new value u k+1 , one obtains the following equation:â Then, the solution step (i.e., line 8 of Algorithm 2) is changed to whereâ k+1 is computed using Equation (31) instead of using Equation (30). Substituting Equation (31) into Equation (32), one obtains After finding u k+1 from Equation (33), we computeâ k+1 using Equation (31) and a k+1 = κa k + (1 − κ)â k+1 . By incorporating the above ideas into Algorithm 2, we can obtain a fixed-point-like method, called Algorithm 4, for solving TVL2 problem (4).
In addition, notice that the linear system in line 7 of Algorithm 3 is equivalent to solving the following least squares problem: Hence, the linear system in line 7 of Algorithm 3 is solved using the CGLS instead of using the CG.

Split Bregman Methods for TVL2 Problems (3) and (4)
In this section, we just provide the alternating split Bregman methods for solving the TVL2 problems (3) and (4) in order to evaluate the performance of Algorithms 2 and 3. For more details about the alternating split Bregman method as well as its convergence analysis, we refer to [5,8].
The following Algorithms 4 and 5 are the alternating split Bregman methods corresponding to the TVL2 problems (3) and (4), respectively.
As was done in line 7 of Algorithm 3, the linear systems in line 4 of Algorithms 4 and 5 are solved using the CGLS instead of using the CG. 8:

Numerical Experiments
In this section, we describe how to carry out numerical experiments for several test problems to evaluate the efficiency of two iterative methods, called Algorithms 2 and 3, using CGLS for solving the new proposed TVL2 problem (4). Performance of Algorithms 2 and 3 is evaluated by comparing their numerical results with those of the existing fixed-point method called Algorithm 1 and the split Bregman methods called Algorithms 4 and 5.
All numerical tests have been performed using Matlab R2016a on a personal computer equipped with Intel Core i5-3337 1.8 GHz CPU and 8 GB RAM. For numerical experiments, we have used three types of PSFs (point spread functions) which are Gaussian blur with standard deviation 9 and Average blur and Motion blur of size 9 × 9. PSF arrays P for Gaussian blur with standard deviation 9 and Average blur and Motion blur of size 9 × 9 are generated by the built-in Matlab function f special( Gaussian , [9,9], 9), f special( average , 9) and P = zeros(9); P(4 : 6, :) = f special( motion , 9, 1), respectively. The blurred and noisy image f is generated by where A stands for the blurring matrix that can be generated by the PSF array P according to the reflexive boundary condition, and E is the Gaussian white noise with mean 0 and standard deviation 3 that can generated using Matlab function E = 3 × randn(m, n), where (m, n) denotes the size of true image.
In order to illustrate efficiency of the proposed algorithms, we have used four test images with an intensity range of [0, 255] such as Cameraman, Lena, House, and Boat with pixel size 256 × 256.
To evaluate the quality of the restored images, we have used the peak signal-to-noise ratio (PSNR) between the original image and restored image, which is defined by where · F represents the Frobenius norm,ũ denotes the restored image of the original image u with size m × n, and u ij stands for the value of the original image u at the pixel point (i, j). In general, the larger PSNR stands for the better quality of the restored image. For all numerical experiments, an initial image was set to the blurred and noisy image f , κ = 1 × 10 −6 , maxit = 150, and tol is set to 5 × 10 −4 (for Algorithms 1 and 5), 2 × 10 −4 (for Algorithm 4) or 1 × 10 −3 (for Algorithms 2 and 3). For the CGLS method that is used to solve a linear system every iteration of Algorithms 1-5, the tolerance for stopping criterion is set to 5 × 10 −2 (for Algorithms 1, 4 and 5) or 1 × 10 −3 (for Algorithms 2 and 3), and the maximum number of iterations is set to 60.

Numerical Results
In this section, we provide numerical results for four test images that are listed in Tables 1-4 and Figure 1. In Tables 1-4, "Alg" represents the algorithm number to be used, "P 0 " represents the PSNR values for the blurred and noisy image f , "PSNR" represents the PSNR values for the restored image, "Iter" denotes the number of iterations required for Algorithms 1-5, the values in parentheses under the "Iter" column refer to the average number of iterations for CGLS, and "α, β, γ, λ" and "δ" denote parameters that are chosen by numerical tries. Notice that, according to Theorem 1, the parameters α, β, γ, λ should be chosen appropriately for good performance of the fixed-point methods.
As can be seen in Tables 1-4, Algorithm 3 restores the true image better than Algorithms 1 and 2. This means that the fixed-point-like method for TVL2 problem (4) restores the true image better than the fixed-point methods for TVL2 problems (3) and (4). The linear system in Algorithm 3 that is obtained by computing a k+1 after the solution step of finding u k+1 is well-conditioned, while the linear system in Algorithm 2 is ill-conditioned. This is the reason why Algorithm 3 restores the true image significantly better than Algorithm 2. Since PSNR values of Algorithm 2 are about 0.3 to 1.0 smaller than those of Algorithm 3 for all test images, numerical results of Algorithm 2 are not provided for House and Boat images in Tables 3 and 4. Figure 1 shows the restored images by the fixed-point method, called Algorithm 1, for TVL2 problem (3) and the fixed-point-like method, called Algorithm 3, for TVL2 problem (4).
The fixed-point method (Algorithm 1) for the TVL2 model (3) performs worse than the corresponding split Bregman method (Algorithm 4), while the fixed-point-like method (Algorithm 3) for the TVL2 model (3) performs almost as well as the corresponding split Bregman method (Algorithm 5). Note that both the split Bregman methods for the TVL2 models (3) and (4) perform the same for almost all cases. When considering the number of iterations required for the CGLS, the total number of iterations for Algorithm 3 is less than that of Algorithm 1. Each iteration of Algorithms 1-3 requires one linear system solver CGLS and two matrix-times-vector operations that are the main time-consuming kernels, and each iteration of CGLS requires two matrix-times-vector operations. In addition, there are some additional vector-update operations that can be negligible as compared with matrix-times-vector operation. This means that the execution time of Algorithm 3 for the new TVL2 problem (4) is less than that of Algorithm 1 for TVL2 problem (3). For example, for the Cameraman image with Gaussian blur, the CPU times for Algorithms 1 and 3 are about 68 and 42 seconds, respectively.

Conclusions
In this paper, we first proposed a new TVL2 regularization model (4) for image restoration, and then we proposed the fixed-point method (Algorithm 2) and the fixed-point-like method (Algorithm 3) for solving TVL2 problem (4). According to numerical experiments, the fixed-point-like method for the new TVL2 problem (4) restores true image better than the fixed-point method for TVL2 problem (4). The reason for this is that the linear system in line 7 of Algorithm 3 is well-conditioned and a k+1 is computed after the solution step of finding u k+1 , i.e., a k+1 is computed using the new value of u k+1 instead of using the old value of u k .
Both of the split Bregman methods for TVL2 problems (3) and (4) perform the same for almost all cases, while the fixed-point-like method for TVL2 problem (4) performs better than the fixed-point method for TVL2 problem (3). It can be also seen that the execution time of the fixed-point-like method for TVL2 problem (4) is less than that of the fixed-point method for TVL2 problem (3). Hence, it can be concluded that the new proposed TVL2 model (4) for image restoration is preferred over the TVL2 model (3), and the proposed fixed-point-like method (Algorithm 3) is well suited for the new TVL2 model (4).
The fixed-point-like method and TVL2 model (4) proposed in this paper can be applied to the image inpainting problem or image restoration problem with Poisson noise. Future work will study these kinds of problems.