Image deblurring via adaptive proximal conjugate gradient method

Abstract

It is not easy to reconstruct the geometrical characteristics of the distorted images captured by the devices. One of the most popular optimization methods is fast iterative shrinkage/ thresholding algorithm. In this paper, to deal with its approximation error and the turbulence of the decrease process, an adaptive proximal conjugate gradient (APCG) framework is proposed. It contains three stages. At first stage, a series of adaptive penalty matrices are generated iterate-to-iterate. Second, to trade off the reconstruction accuracy and the computational complexity of the resulting sub-problem, a practical solution is presented, which is characterized by solving the variable ellipsoidal-norm based sub-problem through exploiting the structure of the problem. Third, a correction step is introduced to improve the estimated accuracy. The numerical experiments of the proposed algorithm, in comparison to the favorable state-of-the-art methods, demonstrate the advantages of the proposed method and its potential.

1. Introduction

Image restoration has been studied for a long time [1][2], which is motivated by the spatial (blur) and point degradations of the original scene. This problem is ill-posed problem, varies in different situations. Each elements of the imaging system contributes to the degradations of the system. In the past years, this subject received an increasing amount of interest [42][43] [44][45][46]. The main goal of image restoration is to restore the original image from the degraded scene. A common imaging model consists of three components, i.e. the observed image s , the original scene v and the presented noise ε. This model can be summarized as follows

where K denotes the convolution operator on clean image v . And the noise level ε varies in different conditions, which includes Gaussian, salt-and-pepper noise, speckle noise [3] and other distortion situations [4][5][6][7]. For recovering more edges and textures or improving the stability of the restoration procedure, some regularization terms [10] are developed to obtain a favorable solution. This procedure requires prior acknowledge necessary about the sampling process or the image data.

A popular regularization term is total variation (TV) proposed by Rudin-Osher and Fatemi [8] (ROF), which can extract the image gradient and preserve edges. And its formation can be presented as follows

where ║*#183;║F denotes Frobenius norm on the matrix space D∈Rm×n, the parameter λ is a positive coefficient. It can be noted that TV regularization term would introduce some difficulties, such as non-smooth of objective function and enormous storage requirement. In the past decades, numerous researchers presented vast majority of optimization methods to tackle this problem. On the one hand, some improvements on this term were presented, e.g. non-local based total variation model [9]. On the other hand, various optimization techniques were presented for accelerating the reconstruction procedure, such as gradient projection-like methods, projected Newton method (PN), gradient projection-reduced Newton method (GPRN) and projection-conjugate gradient method [10], which can be viewed as a unified framework of steeped descent algorithm.

Despite of existed optimization schemes, recently, some efficient optimization methods have been developed. These methods consist of proximal splitting schemes [14][15][17], iterative shrinkage/thresholding algorithm (IST) and its variants [11][12][23], sparse reconstruction by separable approximation (SpaRSA) [13]. Among these methods, IST-like algorithms were examined to be simple and efficient minimization methods. Proximal splitting scheme is a key procedure to the success of these methods. Its key idea is to solve the resulting sub-problem. Another mathematical programming method for this problem is the framework of augmented Lagrange method combined with the variable splitting techniques [18]. It has been devoted to image restoration, such as compressive magnetic resonance (CMR) image reconstruction [31], a hot air turbulence removal and object detection method [25], sparse representation (SR) based image restoration [26] and compressed sensing [27][28].

In this paper, an effective image deblurring method, inspired by the woks in [19][30], is proposed. This scheme can provide a generalization representation of proximal conjugate algorithm in [19] for image de-blurring problem. The motivation of the proposed method is to handle the practical limitations in [19][30], e.g. the inaccuracy estimate of the sub-problem and the dimension limitation. Accordingly, the framework of adaptive proximal conjugate gradient (APCG) is proposed. The proposed method can divide into three stages mainly. First, an improved procedure for approximating inverse Hessian matrix is presented for the resulting sub-problem. Second, the sub-problem is solved by an adaptive regularization gradient projection method, characterized by an iterate-varying hyper-ellipsoidal based norm. Then, an efficient solution to this sub-problem is offered, which can be viewed as an incremental method on the problem and backtracking adaptively. Finally, a conjugate gradient method is applied to adjust the approximate error of the second stage. The presented method can exploit the special structure of the problem and improve the estimated accuracy. The convergence properties, presented in [19], guaranteed the effectiveness of the proposed method. At last, the simulation experiments and the performance evaluation, compared with the state-of-the-art methods, demonstrate the effectiveness and efficiency of the proposed method.

The rest of the paper is organized as follows. In section 2, a review on the methods of IST and FISTA is presented. Section 3 demonstrates our framework and a solution to the resulting sub-problem. In Section 4, extensive experiments and discussion in details under several degradation situations are performed. At last, the conclusion on the proposed method is carried out.

 

2. Preliminaries

2.1 An overview of proximal splitting methods

In this section, a short description of the schemes, i.e. the basic definitions of the proximal gradient scheme, the proximal Newton splitting scheme [19] and fast iterative/thresholding algorithm [12], are presented.

Definition 2.1: proximal gradient scheme is expressed as follows

where ρ is a positive scale parameter, viewed as penalty parameter, and sometimes can be configured with a constant parameter, such as Lipschitz constant. And E denotes the matrix space E∈Rm×n. Its theoretical results have been discussed and analyzed in [10].

Definition 2.2: proximal Newton splitting scheme [19], which can be defined as follows

where H is a generalization representation to the parameter ρ and assumed to be symmetric and positive definiteness matrix. And x can be viewed as approaching point. The symbol H can be chosen to be a proper matrix. Solving the problem (4) can be considered as a procedure to weight different structures locally, and may bring in an alternate way toward the optimal solution. A general solution to this problem may be difficult because of TV (x) and the determination of the weighting matrix H . It can be observed that H may have an influence on the storage requirement and computation complexity, which have not been investigated fully in the work [11][12]. A general scheme for minimizing TV based prlobem is presented in Algorithm 1.

Algorithm 1.Proximal splitting scheme for minimizing TV

However, the performance of proximal splitting scheme depends on the implementation, including step size ti, stopping condition, and initialization parameters despite of its apparently simplicity. A simple solution to ℓ1 -based minimization problem by proximal splitting scheme is IST. The key procedure of this scheme can be given by

where ρ is a chosen threshold for discriminating different coefficients. It can be seen that the procedure is piece-wise. However, this scheme suffers from low convergence rate and approximate accuracy.

To accelerate the convergence rate of IST, a fast and simple numerical scheme, named fast iterative shrinkage/thresholding algorithm (FISTA) [12], is proposed. It should be noted that the proximal splitting method and IST-like algorithms involve with the operation on the approximation point vi-1 - ti▽f(xi-1) and the mechanism to obtain the step size ti.

2.2 Conjugate gradient method

The conjugate gradient method is popular for solving large scale optimization problems. In this subsection, this scheme is intruded to remedy the local increment in FISTA [12]. Then, the backward step is employed for the proposed algorithm. Many of the special effects have been devoted to the conjugate gradient method [32]. A detailed discussion of the theoretical and performance picture of conjugate gradient method have been presented. The convergence results of the proposed method are similar to the work in [7]. So, the proofs of these results are ignored. A general procedure of conjugate gradient method, presented in algorithm 2, is introduced to the proposed algorithm.

Algorithm 2.General congujate gradient method for non-smooth function

The symbol C(·) denotes the procedure for obtaining the conjugate gradient update parameter that belongs to a real number.

was developed by Fletcher and Reeves [35]. proposed by Polak and Ribiere and Polyak [36][37], behaves a better way in a practical way. Although there are some classical methods for choosing this parameter, Eq.(7) is mainly employed in this paper. We extended the algorithm by combining classical general descent algorithm and line search by choosing adaptive step size. The details of these steps are presented in algorithm 3.

 

3. Proposed method

3.2 Adaptive proximal-conjugate gradient method

The principle of the proposed framework is to solve the problem Eq.(2) through striking a balance between the effort expended at each iteration to solve the corresponding sub-problem and the number of iterations required for convergence. The main challenges of the proposed framework have three aspects. On the one hand, when the dimensions growing, especially, the calculation of the approximate inverse Hessian matrix is required to be a suitable formation. Unfavorably, FISTA has not been explored this problem. The accuracy of the proposed framework is completely determined by the approximate accuracy of the sub-problem. Furthermore, the relative cost of the steps of many optimization methods changes when the problem becomes large. Then, certain difficulties arise.

The proposed method for total variation minimization problem, including three stages, is presented in algorithm 3. First stage in the step 3 is provided for overcoming the limitations in [19][30]. It should be noted that this matrix is assumed to be symmetric and positive matrix. Second stage in the step 4 can convert Eq.(3) into a general formation Eq.(4). Or more precisely, the resulting problem can be viewed as hyper-ellipsoidal norm based sub-problem defined in Eq.(8). It should be noted that this generalized representation can product ellipsoidal-norm changes locally. The third stage, including the steps from 5 to 8, is employed for adjusting the approximate error through introducing conjugate gradient method.

Algorithm 3.The framework of APCG

The matrix is not presented to approximate the toeplitz matrix simply. Its goal is to provide a more accurate and positive definiteness penalty parameter matrix to the sub-problem in step 4. The practical importance of can be presented in many ways. First, based the presentation in the literature [39], the matrix generalizes the representation of penalty parameter and dominates the progress toward the feasible solution. Second, empirical experiments, presented in FISTA and alternative direction multiplier method (ADMM) [40], indicate that the selection of is very important. Some theoretical analyses were presented in [19][39][41]. can be calculated in Algorithm 4. It can be noted that this matrix is updated iterate-to-iterate. This algorithm can provide various formations of approximate inverse Hessian matrix.

Remark 1: An approximate solution of the step 4 is presented in algorithm 4. This solution of the sub-problem is improved from the work in [12], which can be expressed as follows

We provide this algorithm without proofs. The advantage of this algorithm can exploit the structure of matrix directly, indicates the computation of the true Hessian matrix is bypassed. All the key aspects of algorithm 4 is the main difference to the work in [19][30].

Algorithm 4.A general approximation method to

Remark 2: The proposed method is improved from the work in [12]. It is extended to a more feasible solution by saving the storage requirement. This weighting scheme can improve the performance and reliability of the algorithm through searching the descent direction adaptively.

Remark 3: The symbol O(·) stands for the operation on which can produce some special structure matrices [20], such as tri-diagonal matrix, up-arrowhead and down-arrowhead. In this paper, the tri-diagonal matrix is taken mainly.

Remark 4: The proposed framework is different from the works [19][31] in several ways. First, the proposed framework extends the original scheme [31] by introducing the conjugate gradient method, which can reduce the approximate error in some degree. Second, the proposed method deals with the resulting problem by taking account into the special pattern of the matrix. A further benefit of this modification is that lower storage requirement and computational cost is required. Third, the computation of exact Hessian matrix is avoided, which can reduce the computational demand significantly when the dimension growing.

Remark 5: For ADMM [40] and FISTA [12] based methods, there are many differences in many ways. First, the main difference lies in the update and the form of the penalty parameter [12]. In contrast, the proposed algorithm generalizes the penalty parameter into a full matrix penalty formation. Second, the solution to ℓ1 -based sub-problem [40] is different to the proposed method. The application of IST in [40] may smooth local structures or texture components in some extent. The geometric illustration, offered in the reference [39], indicates the choice and the form of the penalty parameter has an important influence on the convergence process. In summary, the essential features of the proposed method are the bounded and hyper-ellipsoidal adaption algorithm. The convergence analysis of the proposed framework confirmed that it should be a novel algorithm for this problem.

3.4 Some theoretical results of APCG

Some notations: An iterative sequence vi converges to v*. Based on the most effective technique by comparing the improvement at each step to the improvement at the previous step. An algorithm converges q -quadratically is denoted as follows:

and the and the function f (v) assumed to be locally convex and ▽2f to be locally Lipschitz continuous. Similarly, q -superlinearly is defined as follows

The symbols Q and Qs are related to the rates of convergence of the presented method. And a method would be practicable if iterative process is the closeness to the required point with some rapidity.

Convergence analysis of APCG: some assumptions and theoretical results against the iterative sequence generated by the APCG are presented.

Lemma 1 [19]. Suppose f (v) is twice continuously differentiable and ▽2f is locally Lipschitz continuous with constant Φ. IF satisfies Dennis-More criterion and their eigenvalues are bounded, then the unit step length satisfies the sufficient descent condition after sufficiently many iterations.

Theorem 2. Suppose (1) f (v) is locally strongly convex with constant m , (2) ▽2f is locally Lipschitz continuous with constant Φ , and (3) then the hyper-ellipsoidal adaption algorithm converges q -superlinearly to v* .

Proof. The assumptions of Lemma 1 are satisfied; then the stage of the conjugate gradient satisfy the sufficient descent condition after many iterations:

where ti assumed to be equal to 1. Based on the convergence results (c.f. Theorem 3.3 in [19]), i.e. converges q -quadratically

where denotes the proximal hyper-ellipsoids direction. It should be noted that the eigenvalues of hyper-ellipsoids matrix provided by Algorithm 3 is box-constrained, then the last term can be reformatted as follows

The details of Eq. (13) referred the Proposition 3.5 in [19]. It is notable that the ▽2f and di+1, then it can be observed that

Due to the results (Lemma 2 in [38]) and q -quadratically convergence (Theorem 3.3 [19]), then we can get

putting these expressions into Eq.(13), then we can get

considering the Eq.(16) and Eq.(12), then we can obtain

at last, it can be seen that vk converges to v* q -superlinearly.

3.5 An approximate solution to the problem (2)

The fundamental issues of the framework of APCG are the presence of and the corresponding sub-problem. Moreover, the computational expenses of this framework are referred to two estimated procedures mainly, i.e. the solution to the sub-problem Eq.(6) and the conjugate gradient method. Based on provided by algorithm 4, an improved scheme from the work [12], is proposed in algorithm 5.

Algorithm 5.A hyper-ellipsoids regularization gradient projection method

Remark 6: The operation P stands for the orthogonal projection operator on the convex set C . In this paper, the set C is defined as a bound condition. For the symbol D , it denotes similar convex projection on the dual space. The real definition of the L and the adjoint operator LT can be found in [12].

 

4. Numerical experiments and discussion

4.1 The experiment settings

In this subsection, the experimental settings for evaluating the referred methods are presented. The IST [23], the FISTA [12], FTVD based method [40] (i.e. ADMM) and the proposed method are implemented in Matlab. The assessment procedures are performed on a work station shipped with an Intel processor Core2 Quad 2.66 GHz and 3G RAM. Some synthetic examples are applied to the proposed method, i.e. Gaussian (Gau) blur, uniform (Uni) blur and disk based blur. The Gaussian noise with zero mean and a variance σ2 = 10-3 is added to the blurred image. In this paper, two criteria are used to measure the performance numerically, i.e. signal-to-noise (SNR) and structural similarity (SSIM) index [34].

The settings of the experiments are described in Table 1. All experiments are designed and performed for some goals, such as numerical behavior of the sub-problem, sensitivity analyses of penalty parameter and noise-levels, a wide range of datasets and blurring operations. In order to have a fair comparison between FTVD based method and our method, the regularization parameter λ is chosen to be = 5 × 10-4 from the row 3 to 5.

Table 1.The general configuration of the experiments

All the required methods are performed with default parameters. Meanwhile, all the pixel intensity values are rearranged to the box bound from 0 to 1. And the inner iterative numbers for the sub-problem are chosen to 10. The image dataset can be downloaded from http://www.cipr.rpi.edu/resource/stills/index.html. CIPR datasets, contain about 101 images, including aerial (7), brodatz (53), canon (18) and Kodak (23). The code of FTVD is obtained from the web site: http://www.caam.rice.edu/~optimization/L1/ftvd.

4.2 The first experiment on Airplane

In this subsection, the visual perception and numerical performance of the experiment on the airplane (512×512 ) are demonstrated. The second row of Table 2 exhibits the maximum scores of the SNR and the SSIM, which are obtained by performing 100 iterates. It can be notable that the proposed method is better than other methods separately, which is more than 1dB in term of the SNR. In another viewpoint on these objective results, the proposed method performs a better behavior. Fig. 1 demonstrates the observed data and the visual results of these methods. The restoration image is more similar to the original image. It can be noted that the proposed algorithm restores more local features, for example the textures, edges and shapes. As an evaluation of the reconstruction process, Fig.2 displays the curves of the SNR and SSIM at every iteration. The quantitative comparison of these methods indicates that our method is better than the other methods.

Table 2.The numerical results of the 1st and 2nd experiments (λ = 1×10-3)

Fig. 1.The reconstruction results of the first experiment (λ = 10-3). (a) The observation image, (b) for IST, (c) for FISTA and (d) for APCG.

Fig. 2.The SNR and the SSIM of the first experiment (λ = 1×10-3). (a) for the SNR, (b) for the SSIM.

4.3 An evaluation on the solutions to the sub-problem (8)

In this subsection, an assessment on the sub-problem was performed. It can be noted that is a identify matrix in FISTA based method. First, the reconstruction results were examined in Lena image. Second, the times of inner iteration for solving was measured and evaluated. These experiments can provide a comprehensive view on the numerical behavior of the proposed method.

The visual outcomes are presented in Fig.3. In comparison to the other methods, the proposed method can restore a sharper and clear latent image. As shown in Table 2, it is notable that the numerical results, based on the indices of the SNR and SSIM, imply that the proposed method is better than the other methods again. The total iterate numbers of APCG based method are presented in five column of Table 2.

Fig. 3.The reconstruction results of the 2nd experiment. (a) The observation image, (b) for IST, (c) for FISTA and (d) for APCG.

Fig. 4.The SNR and SSIM of the 2nd experiment (λ = 1×10-3). (a) for the SNR, (b) for the SSIM.

The influence on the resulting sub-problem based is presented in Fig. 5. The five column of Table 2 presents the total iterate numbers by summarizing all the minor iterates Ii when solving Eq.(6). It is evident that the IST conducts the lowest iterate numbers. However, compared to the FISTA (Itol = 837), the proposed method achieves a more performance improvement and involves with a lower total iterate numbers (Itol = 770). This phenomenon can be observed in Fig. 5 carefully. Although the curve of FISTA decreases from 10 to 6 slowly, it can be noted that the outcome of the APCG is mainly equal to 6 and hold it for many times, which can be observed in Fig. 5 about Ii=50-100.

Fig. 5.The inner iteration numbers of the first and second experiments (λ = 1×10-3).

4.4 Sensitivity analysis of penalty parameters

In this subsection, some experiments are performed for observing the numerical behavior of FISTA and FTVD in various penalty parameters. In another view on these numerical results, it can be seen that the update process may not be effective in any iterate. In Table 3, the symbol * denotes the fail of the algorithm. The symbol ρ defined in FISTA is penalty parameter, and β for FTVD. It is notable that the generalization representation of the proposed framework suggests an alternative approach to this problem. When β in FTVD decreases, the behavior of the SNR history would act a different way.

Table 3.The evaluation of the FISTA(ρ) and FTVD(β) on 'Lena' image(λ = 5×10-4)

4.5 Sensitivity analysis of noise-levels

Fig. 6 demonstrates the influence of the noise-levels on these optimization methods. It also can be noted that the IST generates a more stable curve, which implies that it is not sensitive to the parameter σ2. These appearances may relate to the hard thresholding on the latent image. This operation may keep regional features away. For the APCG, although it behaves quite more sensitive to noise-levels increasingly, but it should be noted that its reconstruction results, evaluated by the indices of SNR and SSIM, are better than other methods roughly.

Fig. 6.Sensitivity analysis of noise-levels on Lena (λ = 1×10-3), (a) for the SNR, (b) for the SSIM.

4.6 The numerical results of CIPR still images

The performance of the proposed method is verified in many extended experiments. These experiments are performed in three different blur operations. Similar to the previous experiments, the better numerical results of the referred methods are displayed in bold.

The outcomes of all these methods for typical degradation conditions are presented in Table 4 and Table 5, separately about SNR and SSIM. It can be seen that our method fulfills a substantial performance improvement in comparison to the other methods. The last two rows present the mean and variance of the numerical results of these blurring conditions. The outcomes of these experiments about these four datasets indicate that our method can be applied to a wider variety of images. Meanwhile, it can be seen that the means and variances presented in Table 4 and Table 5 demonstrate that the proposed method behaves more stable than other methods. For the outcomes of the variance of SSIM, FTVD based method is better than other methods. Our approach is a natural extension in the context of proximal splitting scheme. In part, this is due to the introduction of APCG framework.

Table 4.The average SNR of CIPR still images(λ = 5×10-4)

Table 5.The average SSIM of CIPR still images(λ = 5×10-4)

4.7 Discussion

In this paper, we have show that the efficiency of APCG, which involves the unified framework and its approximate solution. It is obvious that the presented method acts a better way than other methods in term of the improvement on the SNR and SSIM, the significant decrement in total iterate numbers Itol, sensitivity analysis of penalty parameters and noise-levels. The numerical results and the visual inspection of the restoration pictures demonstrate an involvement of the generalization representation of the proximal term and an efficient strategy to linear search effectively, leading to the success of the proposed framework. In another viewpoint on this condition, the association of the generalized proximal term with the conjugate gradient restores more high frequencies information from the degraded image. Taken together, these results indicate that the combination of these elements is essential to the achievement of the reconstruction performance. In term of the indices of the SNR and the SSIM in CIPR datasets and blurring conditions, the outcomes of the proposed method are available for a variety of datasets. In the meanwhile, an insight on the inner iterate numbers Ii at a single iterate of the APCG, presented in Fig. 5(a) and Fig. 5(b), implies that the approximate error at every iterate can be alleviated in some degree. The summary of total iterations exhibited in Table 2 verifies the effectiveness of the proposed method again. In contrast to other methods, this difference may relate to the method used. It is possible that the composite framework with the hyper-ellipsoidal norm based proximal term and the backtracking line search mechanism lead to the success of the APCG.

However, there are some limitations on our method. First, the FTVD based method, implemented by ADMM, provides a well balance between the complexity and the run-time of the algorithm. However, we argued that the proposed framework can provide a unified representation of penalty parameter. Second, the proposed method is required to calculate the Algorithm 4. Although this may introduce some computational expense, the reduction in the total iterate numbers of the solution to the sub-problem compensates this extra computational cost to some extent. Third, the storage of the matrix may be required to be huge when the size of the problem growing.

In this paper, an alternate approach is proposed to handle the existed problems of the work in [19][30]. Despite of the extra issues of APCG based method, the evaluation and outcomes of the experiments confirm us that the proposed method is superior to the other methods and is suitable for large-scale optimization problems. Thanks to the accurate estimation of the iterate variable, the proposed method behaves a better way in these experiments.

 

5. Conclusion

In this paper, an efficient optimization framework for the image de-blurring problem, named adaptive proximal conjugate gradient method, is presented. This framework can be viewed as a variable hyper-ellipsoidal norm based method, which is formulated by generating weighting matrix adaptively and offering a good linear direction. The proposed framework may be of functional significance to the other large scale problems. To deal with the practical limitations of the proposed framework, an approximation solution to the sub-problem is presented. Observed from the experiments, the proposed algorithm can provide an alternate approach to total variation minimization problem. The experiment results based on the pictorial and numerical results of the referred methods, indicate that the proposed method is superior to the other methods.