Total Variation Image Restoration Method Based on Subspace Optimization

The alternating directionmethod iswidely applied in total variation image restoration.However, the search directions of themethod are not accurate enough. In this paper, one method based on the subspace optimization is proposed to improve its optimization performance. This method corrects the search directions of primal alternating direction method by using the energy function and a linear combination of the previous search directions. In addition, the convergence of the primal alternating direction method is proven under some weaker conditions. Thus the convergence of the corrected method could be easily obtained since it has same convergence with the primal alternating direction method. Numerical examples are given to show the performance of proposed method finally.


Introduction
Digital image restoration has a wide application in various areas including Navigation, Aerospace, and Biomedicine (see [1][2][3] and the references therein).In general, the relationship between the original image f ∈   and the observed image  ∈   is where  ∈   is the additive noise and the spatial-invariant matrix  ∈  × represents the degradation system caused by problems such as motion blur, distortion radiation, and distortion wavelets in seismic imaging,  =  1  2 with  1 and  2 being the number of rows and columns, respectively, when images are expressed as a matrix.In particular, matrix  represents Block Toeplitz-plus-Hankel matrix (with Toeplitzplus-Hankel blocks) or Block Toeplitz matrix (with Toeplitz blocks) when Neumann boundary condition and the zero boundary condition are used, respectively [4].The objective of the image restoration is to estimate the original image f according to some a priori knowledge about the degradation system , the additive noise , and the observed image .However, it often tends to be very ill-conditioned when the reverse process of model ( 1) is only used to get an ideal estimation  * .Thus one of the effective ways to solve these problems is to combine some a priori information of the original image and define the regularization solution; that is,  * is a minimizer of the following cost (energy) function: where  :   →  is the measure of the difference between  and .The regularization term Φ embodies a priori information and a regularization parameter  > 0 is used to control the tradeoff between the terms  and Φ.
In order to overcome the drawback, in [16], the TV regularization is modified as where 0 <  ≪ 1. Obviously, the regularization is differentiable and many classical optimization methods can be applied [24,25].Unfortunately, experimental results show that  must be small enough to keep the restoration quality, but the efficiency of algorithms will be reduced with smaller  [16,20].By using the following modified TV regularization, a primal-dual active set method is proposed in [17], but it has no rotational invariant [8].
Recently, the alternating direction method [26] is used to solve the  2 -TV and  1 -TV image restoration models, and better performance has been obtained [18][19][20][21][22][23]27].It divides original problem into some simple subproblems and solves them alternatively.The downside is that the search directions are not accurate enough which could reduce the performance of algorithms (see the analysis in Section 3).
Based on the considerations above, we will take the algorithms in [18,20] as examples and present a method to improve the performance of the alternating direction method used in  2 -TV and  1 -TV models.By means of the subspace optimization [28], it corrects the search direction of the primal alternating direction method by using the energy function and a linear combination of the previous search directions.In addition, the convergence of the primal alternating direction method is proven under some weaker conditions.Thus the convergence of the corrected method could be easily obtained by the equivalence between them.Numerical examples are given to shown the performance of proposed method.
The outline of this paper is as follows.Some preliminaries are stated in Section 2. A  1 -TV image restoration algorithm based on subspace optimization is proposed in Section 3. The convergence analysis and numerical examples are given in Sections 4 and 5, respectively.Finally, some concluding remarks are given in Section 6.

Preliminaries
In this section, we first propose some basic definitions and properties; then the alternating direction method to solve  1 -TV model in [20] will be introduced.The alternating direction method in [18] ( 2 -TV model) will be omitted here since it is more simpler.But the comparison experimental results of the method and its corrected method still will be listed in Section 5.2.In the following, let  ℎ and  V be the one-side difference matrix on the horizontal direction and vertical direction, respectively,  = ( ℎ ,  V ) ⊤ , Null() is the null space of matric ,   is the th iteration, and  () is the th element of vector .
Definition 1 (see [20]).An operator P is called nonexpansive if for any  1 ,  2 ∈  ⊂   , we have Specially, P is called -averaged nonexpansive if there exists some nonexpansive operator A and  ∈ (0, 1) such that P = (1 − )I + A, where I is the identity operator.
Lemma 5 (see [30]).Let  :   →  be a closed, proper, and coercive function.Then the set of minima of  over   is nonempty and compact.
In [20], the where the parameters  1 > 0 and  2 > 0 are used to ensure the closeness of  and ,  and , respectively.When applied to (13), the alternating direction method is The first step in ( 14) is equivalent to solving the nonsingular linear system: It can be solved directly; that is, To let computational cost be low, many classical optimization and numerical methods such as CG method and preconditioned iteration methods can be used to solve it [24,25,31].As in [4,18], we suppose that the Neumann boundary conditions are used and the blurring function is symmetric in this paper; then  is a (block) Toeplitz-plus-Hankel matrix and it could be diagonalized by a discrete cosine transform matrix.Thus we could solve the nonsingular linear system utilizing lower cost, since the inverse of blurring matrix  can be computed by using fast cosine transforms; please see [4,18,27] for more details. Since the second step in ( 14) is equivalent to solving the minimizers of  functions in the form of ]() = | − | 2 +  1 ||.We can know from Lemma 2 that the th element of its solution  +1 is The third step in ( 14) is a  2 -TV model with the blurring matrix  is the identity matrix, that is, image denoising problem.Many methods can solve it effectively [13][14][15], and Chambolle' projection algorithm [13] is used in [20].

The Method Based on Subspace Optimization
Based on the principle of alternating direction method above, it is easy to know that the solution  +1 = [ +1 ;  +1 ;  +1 ] of ( 14) is different from the solution of ( 13), so the performance of solving the minima of ( 13) may be lowered.Thus a method based on subspace optimization [28] will be proposed to improve the optimization performance of ( 14) in this section.Firstly, we introduce subspace optimization method.Let  +1 = [ +1 ;  +1 ;  +1 ] denote the solution corrected by subspace optimization method, and where Supposing  +1 ∈  +1 is selected as (1, 0, . . ., 0) in (20), then the equation  +1 =  +1 is valid.In fact, we can know from (21) that  +1 is a minimum point of the cost function: ( +1 ) = (  +  +1 ); thus ( +1 ) ≤ ( +1 ) must be true.It is helpful to improve the performance of original alternating direction method on every iteration, and the key problem is how to solve (21) now.To reduce the computational cost, let  +1 = [ +1 0 ] in this paper.Next, the correction method of (14) will be proposed according to three cases.

Case 1 (𝜔 (𝑖) ̸
=  () and  1   ̸ = 0 or  2   ̸ = 0 for ∀ ∈ ).Obviously, the function (, , ) is smooth; we have where   denotes the unit column vector with its th element being 1 and the others being zero, It is expensive to compute  2 ‖‖ TV / 2 , so we let  2 ‖‖ TV / 2 = 0 to improve the efficiency of the algorithm because ∇ 2  is independent of variable  now.
On the other hand, since  +1 is an optimum of the convex function , we can know from (21) and Taylor series expansion that Equation (23) implies that Obviously, vector  +1 ̸ = 0; then Thus and  +1 −   =  +1 −   .The two equalities do not hold in image restoration since the dimension of vector  is very huge.Therefore ( +1 ) ⊤ ∇ 2  +1 > 0, and from (24).
In summary, the proposed algorithm is as follows.

Convergence Analysis
As known from ( 20) and ( 26), when lim →∞ ‖  0 ‖ 2 = 0 and one of the sequences {  } and {  } is convergent, then the other sequence is also convergent, where . Thus a simple proof of the convergence of algorithm (14) will be given below without the condition that H  1 (H  (⋅, )) and H TV (H  (, ⋅)) are asymptotically regular which is needed in [20], where Lemmas 6 and 7 below appeared in [18,20], but the proofs given by us are simpler.Lemma 6. Suppose () ∩ () = {0}; then the sets of fixed points of operators T 1 and T 2 are nonempty, respectively.
Proof.It is easy to see that the objection function  in (13) satisfies the conditions of Lemma 5 [18,20].Then  has at least one minimizer (, , ) that cannot be decreased by the alternating scheme (14).Thus  Proof.From Lemmas 2 and 3, we have This completes the proof.
Remark 10.If  ∈ Null(), then  () =  for all  ∈ , where  is a nonzero constant.Since matrix  ≥ 0 in image restoration, so  is a nonzero vector.It follows that the assumption Null() ∩ Null() = {0} holds in general.

Experimental Results
In this section, some numerical results will be provided to show the performance of our method.As in [18,20], the tested images are selected as Cameraman of 256 × 256 and Lena of 512 × 512.The algorithm in [20] and ours will be compared in Section 5.1 ( 1 -TV model).The algorithm in [18] and its corrected algorithm will be compared in Section 5.2 ( 2 -TV model).All the computational tasks are performed using MATLAB 2016a with Core(TM)2CPU with 2.83 GHz and 3.87 GB of RAM.The average value of ten tests would be selected.
The tested blurring function is chosen to be truncated 2D Gaussian function: Here three sets of parameters are chosen: In all runs, CPU time is used to compare the efficiency, signal-to-noise ratio (SNR), and peak signal-to-noise ratio (PSNR): which are used to measure the quality of the restored images.The stopping criterion of the algorithms should satisfies 5.1.Comparison Experiment for the  1 -TV Model.In this subsection, a comparison of the algorithm in [20] with the proposed method is made under different non-Gaussian additive noises.Firstly, for the uniform noise, we let  = 0.008,  1 = 0.05, and  2 = 0.2 as in [20].The uniform random numbers appear in the intervals (0, 0.05), (0, 0.1), and (0, 0.2), respectively.Specially, [20] points out that the noise added to the blur image is independent and identically distributed, so the same set of parameters could be selected for different images when the same kind of noise is considered.This could reduce the computational cost of searching regularization parameters.For the speckle noise, as in [20,23], let  ≡ 1,  1 = 2.5 ×10 4 , the continuation scheme about  2 will be used, easier subsequences with smaller  2 can be solved quickly, and the later subproblems can also be solved relatively quickly with warm starts from previous solutions.The noise variance will be selected as 0, 0.01, and 0.05 respectively.To be fair, the same methods will be used to solve the primal alternating direction method in this paper.Now, we select two cases that have significant improvement of the visual sense to illustrate our method's restoration quality.Under different conditions, more image restoration results with uniform noise and speckle noise will be summarized in Tables 1 and 2, respectively.Figure 1  From Figures 1 and 2, we see that the quality of images restored by our method is improved obviously compared to those restored by the method of [20] in the two cases.For the uniform noise and speckle noise, Tables 1 and 2 show that the SNR and PSNR of the images restored by our method both are higher than the results of [20] in most cases.For the same requirement of relative error, our proposed algorithm is clearly faster than the competing algorithm in [20].Referring to Sections 2 and 3, the proposed algorithm uses subspace optimization method to correct the search direction, which could ensure ( +1 ) ≤ ( +1 ) is valid on the th iteration.
So it is clear that the proposed method is more efficient.The analysis is verified by the experiments results above.

Comparison
Experiment for the  2 -TV Model.In the subsection, the algorithm in [18] and ours will be compared.
Being different from [19,20], [18] involves a fitting of the auxiliary variable to  which has superior performance compared to fitting  [18,23,32].The additive noise is selected as Gaussian noise.Here four sets of the standard deviation of noise are chosen: 0.01, 0.05, 0.1, and 0.2.Same methods appearing in [18] will be used to solve the primal alternating direction method.Similarly, we will first select two cases that have significant improvement of the visual sense to illustrate our method's restoration quality.Figures 3(a) and 4(a) show the original  Figure 3(c) shows the restored Camera image by algorithm in [18]; Figure 3(d) shows the restored Camera image by our method, where the standard deviation of noise is 0.1 and the support is equal to 9 × 9. Figure 4(c) shows the restored Lena image by algorithm in [18]; Figure 4(d) shows the restored Lena image by our method, where the standard deviation of noise is 0.2 and the support is equal to 9 × 9. Table 3 shows more image restoration results with the different support and standard deviation of noise.
As in Section 5.1, Figures 3 and 4 show that the quality of restoration images by using our method is better than those restored by the algorithm of [18] in the two cases.We also can know from Table 3 that the SNR and PSNR of the images restored by our method both are higher than the results of [18] in most cases.With the same reasons mentioned in  1 -TV model, all computational time required by our method is significantly less than that required by the algorithm in [18].So the performance of the algorithm in [18] also has been improved through the correction.

Concluding Remarks
Taking the  1 -TV model in image restoration as an example, a modified method is proposed to overcome drawbacks of the primal alternating direction method in this paper.We have illustrated that our method about how to correct the search direction improves the optimization performance.In addition, the convergence of the primal alternating direction method has been proven under some weaker conditions, and thus the convergence of proposed method is easily obtained by the equivalence between them.The experimental results based on two models in [18,20] show that the proposed method could enhance the quality of restored images in most cases and the efficiency of algorithms has been significantly improved.In fact, our method can be applied to many other cases optimized alternatively.

( 31 )Lemma 7 .
Therefore ,  are the fixed points of T 1 and T 2 , respectively.The operators T 1 and T 2 are nonexpansive.

Figure 1 :
Figure 1: The restored Camera images by different methods with the uniform noise.The random numbers are in the interval (0, 0.2), the support being equal to 9 × 9. (a) Original image.(b) Observed image.(c) Image restored by algorithm in [20].(d) Image restored by our method.

Figure 2 : 2 2
Figure 2: The restored Lena images by different methods with the speckle noise.The noise variance is 0.01, the support being equal to 7 × 7. (a) Original image.(b) Observed image.(c) Image restored by algorithm in [20].(d) Image restored by our method.

Figure 3 :
Figure 3: The restored Camera images by different methods with Gaussian noise.The standard deviation of noise is 0.1; the support is 9 × 9. (a) Original image.(b) Observed image.(c) Image restored by algorithm in [18].(d) Image restored by our method.

Figure 4 :
Figure 4: The restored Lena images by different methods with Gaussian noise.The standard deviation of noise is 0.2; the support is 9 × 9. (a) Original image.(b) Observed image.(c) Image restored by algorithm in [18].(d) Image restored by our method.
(a)  shows the original Camera image.Figure1(b)shows the observed Camera image, where the noise is the uniform noise, the random numbers are in the interval (0, 0.2), and the support is equal to 9 × 9. Figure1(c)shows the restored Camera image by the algorithm in[20].Figure1(d)shows the restored Camera image by our method.Figure2(a)shows the original Lena image.Figure2(b)shows the observed Lena image, where the noise is the speckle noise, the noise variance is 0.01, the support is equal to 7 × 7.

Figure 2 (
c) shows the restored Lena image by the algorithm in[20].Figure2(d)shows the restored Lena image by our method.

images. Figures 3
(b) and 4(b) show the observed images.

Table 1 :
Summary results for the uniform noise.

Table 2 :
Summary results for the speckle noise.

Table 3 :
Summary results for Gaussian noise.