Abstract

This paper proposes a new method, bound alternative direction method (BADM), to address the minimization problems in image deblurring. The approach is to first obtain a bound unconstrained problem through bounding the regularizer by a novel majorizer and then, based on a variable splitting, to reformulate the bound unconstrained problem into a constrained one, which is then addressed via an augmented Lagrangian method. The proposed algorithm actually combines the reweighted minimization method and the alternating direction method of multiples (ADMM) such that it succeeds in extending the application of ADMM to minimization problems. The conducted experimental studies demonstrate the superiority of the proposed algorithm for the synthesis minimization over the state-of-the-art algorithms for the synthesis minimization on image deblurring.

1. Introduction

The mission of image deblurring is to restore an original image from noisy blurred observation modeled as where (stacking an -image into an -vector in lexicographic order), is the matrix representation of a convolution operator, and is Gaussian white noise. This imaging inverse problem has recently inspired a considerable amount of research ([111] and further references therein) in which the optimization problems fall into two varieties: synthesis formulation and analysis formulation [12] which are detailed below.

1.1. Formulations and Algorithms

In the synthesis formulation, the unknown image is represented as , where is a wavelet frame or a redundant dictionary and are the sparse coefficients estimated usually via one sparsity-promoting regularizer, yielding where is the positive regularization parameter and is the regularizer, typically the norm. The deblurred image is then obtained by .

In the analysis formulation, as opposed to (2), it minimizes the cost function with respect to : where is a sparsifying transform (such as wavelet or finite differences) and analyzes the image itself against the coefficients that in (2) works on. If is a wavelet frame, then usually where ; if is a matrix representing finite differences at horizontal and vertical directions, then is the discrete anisotropic or isotropic total variation [8, 9].

In the past several decades, a lot of algorithms have been proposed to solve (2) and (3). Among these, the iterative shrinkage/thresholding (IST, [2]) algorithm can be considered as the standard one. However, IST tends to be slow in particular when in (1) is poorly conditioned. To overcome this problem, several fast variants of IST have been proposed. They are TwIST [13], FISTA [14], and SpaRSA [15]. TwIST is actually a two-step version of IST in which each iterate depends on the two previous iterates, rather than only on the previous one (as in IST); FISTA is a nonsmooth variant of Nesterov’s optimal gradient-based algorithm for smooth convex problems [16, 17]; and SpaRSA adopted more aggressive choices of step size in each iteration. These three algorithms, which only use the gradient information of the data-fitting term, have been shown to clearly outperform IST. Some other efficient algorithms, using the second-order information of the data-fitting term, have also been proposed. A noble representative is the so-called split augmented Lagrangian shrinkage algorithm (SALSA) [10, 11], which was found to be faster than the previous FISTA, TwIST, and SpaRSA when the matrix inversion (where is a positive parameter) can be efficiently computed (e.g., if is a circulant matrix, then the matrix inversion can be efficiently calculated by FFT). Actually, SALSA is an instance of alternating direction method of multipliers (ADMM) [18, 19] which has a close relationship [20] with the Bregman iterations [2124], amongst which, the split Bregman method (SBM) [22] has been recently frequently applied to handle imaging inverse problems [2527].

1.2. Minimization

Convergences of above algorithms are guaranteed, benefiting from the fact that the regularizer in (2) and (3) is usually convex. In this paper, a nonconvex (also nonsmooth) regularizer, that is, the () norm, is adopted since using the norm is able to find sparser solution than using the norm which was demonstrated in many studies [2832]. Thus, there comes the minimization problem of the synthesis formulation: and of the analysis formulation: where .

Next, a brief survey on the literature of minimization is given. Recently, remarkable attention has been paid to the minimization problem: where . Some sufficient conditions of (6) have been established in [29, 3335]. Many efficient reweighted minimization methods address (6) through iteratively solving where is the estimated signal at the th iteration; is the data-fitting term which preserves the consistency; for instance, for the Gaussian noise; is the regularization term at the iteration, which is aiming to encourage the desirable properties of the target such as sparsity. Most of reweighted minimization methods (see [2, 28, 30, 31, 34, 36, 37] and many references therein) focus on the reweighted (IRL1) and the reweighted (IRL2) minimization, and their corresponding regularizers are and , respectively, where is the component of and is the weight of at the iteration and it is a function of the previous estimated component . In IRL1, usually, , while in IRL2, where . Notice that is always set to 1 (e.g., IRL1: [28] and IRL2: [2]) or (e.g., IRL1: [34, 36, 38] and IRL2: [30, 31, 38]), even [32]. It is worth noting that above weights are separable, even though there also exist some reweighted minimization methods [39, 40] with inseparable weights.

1.3. Proposed Approach

In this paper, the minimization problem will be used for image deblurring, but the resulting problem cannot be efficiently solved by above minimization methods stated in Section 1.2, since they usually adopt slow iterations (such as the generic IST algorithm), and in image deblurring, the matrix-vector products are always computationally expensive. Considering this, the ADMM is used in this paper to tackle the resulting minimization problem because as stated in Section 1.1, the image deblurring problems can be efficiently handled by the ADMM. The proposed approach is to first bound it via a variant of the majorizer proposed in [4, 7], obtaining a bound unconstrained one; then reformulate it into a constrained problem via the technique of variable splitting [41, 42]; and lastly attack this resulting constrained problem using an augmented Lagrangian method (ALM) [43]. If the majorizer that is proposed in [4, 7] is used to bound the original minimization unconstrained problem, then the components of the estimate in the iterations would be stuck at zero forever if they became zero, which very likely prevents convergence to a minimizer. To overcome this problem, in this paper, a variant of this majorizer is proposed by adding a small positive parameter that shrinks in each iteration. The obtained bound unconstrained problem is reformulated into a constrained one through variable splitting which splits the original variable into a pair of variables, serving as the arguments of the data-fitting term and the regularizer, respectively, under the constraint that these two variables have to be equal. This resulting constrained problem is then attacked by an ALM, obtaining a Lagrangian function with two variables resulting from variable splitting. Next, the Lagrangian function is alternatively solved with respect to the two variables, leading to the method called bound alternative direction method (BADM), which is equivalent to the combination of the reweighted minimization method and the ADMM and is able to extend the application of ADMM to the minimization problems.

The proposed BADM has only cost in each iteration in solving the synthesis formulation with a normalized Parseval frame. Experiments on a set of benchmark problems show that the BADM for (4) is favorably competitive with the state-of-the-art algorithms FISTA [14], SALSA [11], and SBM [22] for (2).

1.4. Terminology and Notation

In this section, some useful elements of convex analysis will be given. Let be a real Hilbert space equipped with the inner product and associated norm . Let be a function and let be the class of all lower semicontinuous convex functions that are not equal to everywhere and are never equal to . The proximity operator of is defined as where is positive. If , then is unique [44]; if , the indicator function of a nonempty closed convex set , then becomes the projection of onto , and in this sense, (8) is therefore a generalization of projection operator; and if is the norm, then (8) becomes a well-known soft thresholding: where is the componentwise multiplication between two vectors of the same dimension and sign  is the sign function.

2. Bound Alternative Direction Optimization

Consider a regularized unconstrained model where and is a smooth function with -Lipschitz-continuous gradient and is bounded below. It is unsuitable to directly apply the existing proximal splitting algorithms (such as IST, FISTA, TwIST, and SpaRSA discussed in the section of Introduction) to solve (10), since, for , the nonconvex nature of blocks the use of the proximity operator (see (8)), which is well defined only on the functions which belong to . To overcome this problem, a bound optimization approach will be considered in this paper.

2.1. Bound Optimization
2.1.1. Bound (x) via the Majorizer Proposed in [7]

Since , it is reasonable to bound through where with where . It is easy to see that with equality for , since with equality for . One benefit from the above bound is that the following lemma holds.

Lemma 1. Given , belongs to with respect to .

Proof. Given for , if , then is an affine function of and thus , since ; if , then equals if and otherwise, leading to . Therefore, .

Therefore, the proximity operator of given is obtained by where with

Notice that soft for any .

From above, a closed-form solution of the proximity operator of can be obtained after the bound operation. However, in many iteratively minimization algorithms discussed in the section of Introduction, the proximity operator of the regularization term is commonly used. For , a nature way is to set as the previous estimate and obtain the current estimate by computing the proximity operator of ; that is, where is an iteratively temporary variable which differs in each algorithm. According to (13) with (14), if a component of becomes zero forever, then this component of is set to zero and will be stuck at zero forever, which may prevent convergence to a local minimizer, letting alone a global one. To overcome this shortcoming, a new bound method is proposed below.

2.1.2. Proposed Bound Method

A method is presented that bounds through where with where and is a small positive parameter. It is clear that with equality for (where is a vector composed of all zeros) or with . has the following property.

Lemma 2. belongs to with respect to .

Proof. For any ,  ; thus is an affine function of and thus . Therefore, .

Therefore, a closed-form solution of the proximity operator of can be obtained where .

2.1.3. Iterative Bound Optimization (IBO)

Now, it is ready to propose a framework of IBO as in Algorithm 1.

(1)Set , and choose an arbitrary , and .
(2)repeat
(3)
(4)
(5)
(6)until some stopping criterion is satisfied.

A key observation is that the sequence , generated by the IBO, approaches to zero as , such that (10) can be solved by the IBO which iteratively solves a sequence of problems: where is the th component of . Any accumulation point of the sequence generated by the IBO is a first-order stationary point of (19), which is guaranteed by the following theorem.

Theorem 3. Let the sequence be generated by above IBO and suppose that is an accumulation point of ; then is a first-order stationary point of (19).

Proof. IBO is actually a specific case of the reweighted () method proposed in [45], such that this theorem is also a specific case of Theorem  3.1 in [45].

2.2. Bound Alternative Direction Method

Considering the unconstrained optimization problem that corresponds to Step 3 of IBO, Using the technique of variable splitting, (20) becomes a constrained optimization problem: The rationale behind variable splitting is that it may be easier to solve (21) than it is to solve (20). The augmented Lagrangian function of (21) is where is the vector of Lagrangian multipliers and is the penalty parameter. According to the augmented Lagrangian method (ALM) [43], (21) can be solved through repeating the following iterative process until some stopping criterion is satisfied: minimize (22) with respect to and fixing and setting to the previous estimate of and then update , yielding a variant of ALM called bound ALM (BALM) (Algorithm 2).

(1)Set , and choose , , , and an arbitrary .
(2)repeat
(3)
(4) =
(5) =
(6)
(7)until some stopping criterion is satisfied

Notice that the terms added to in the definition of (see (22)) can be reformulated as a single quadratic term, leading to a variant of BALM (Algorithm 3).

(1)Set , and choose , , , and an arbitrary .
(2)repeat
(3)
(4) =
(5) =
(6)
(7)until some stopping criterion is satisfied.

In (Algorithm 3), corresponds to that these two parameters are used in BALM. It is usually difficult to simultaneously obtain and in Step 3. To overcome this difficulty, the technique of nonlinear block Gauss-Seidel (NLBGS) [46] is naturally used, in which the objective function of Step 3 is solved by alternatively minimizing it with respect to and , while keeping the other variable fixed, yielding the following bound alternative direction method (BADM) (Algorithm 4).

(1)Set , and choose , and an arbitrary .
(2)repeat
(3)
(4)
(5)
(6)
(7)
(8)until some stopping criterion is satisfied.

In (Algorithm 4), Step 3 is equivalent to , while Step 4 (see (13) and (9)) is equivalent to Moreover, since the objective function of (10) is nonconvex for , BADM cannot be guaranteed to converge to a global optimum. Nevertheless, as stated in Theorem 3, the proposed algorithm is able to obtain a stationary point, which in practice is always a good quality deblurred image.

3. Image Deblurring Using Synthesis Formulation and BADM

In this section, the synthesis formulation (see (4)) and the BADM are applied to image deblurring. Note that using the analysis formulation can be naturally extended and will be left as future work. Therefore, here , and is the proposed majorizer of   at iteration, which are inserted into the BADM yielding (Algorithm 5), where Step 4 is derived from (18). Assume that represents a (periodic) convolution and is a normalized Parseval frame (i.e., , where is the identity matrix and possibly ). According to the Sherman-Morrison-Woodbury matrix inversion formula, Moreover, under the periodic boundary condition for , is block circulant, such that can be diagonalized by two-dimensional discrete Fourier transform (DFT). Let ; then is equivalent to a filter in the Fourier domain and the cost of products by using FFT algorithms is [10]. Thus, where . As the computational complexity analysis in [10], the cost of computing is . Moreover, computing is the soft thresholding whose cost is and computing also has cost. Therefore, each iteration of the BADM for (4) has cost.

(1)Set , and choose , , , and an arbitrary .
(2)repeat
(3) =
(4) = soft
(5) =
(6) =
(7)
(8)until some stopping criterion is satisfied.

4. Experiments

In this section, the proposed algorithm BADM for (4) is compared with the state-of-the art algorithms: FISTA [14], SALSA [11], and SBM [22] for (2) (from now on, the BADM is only for (4) while FISTA, SALSA, and SBM are for (2) if without specification). Consider the low-frequency images (Cameraman), high-frequency images (Mandril), and both low- and high-frequency images (Lena) (see Figure 1), with size pixels, corrupted by the following three benchmark cases [5, 11]: (I) uniform blur of size and noise variance (termed UNI); (II) Gaussian blur kernel with (termed GAU); (III) the point spread function of the blur operator is for and (termed PSF). All the experiments were performed using MATLAB on a 64-bit Windows 7 PC with an Intel Core i7 3.07 GHz processor and 6.0 GB of RAM. In order to measure the performance of different algorithms, the following four metrics (where is the original image and and are the observed image and the estimated image at the iteration, resp.) are employed: (a) consumed CPU time (Time); (b) number of iterations (Iterations); (c) mean square error (); and (d) improvement in SNR () and the stopping criterion is chosen as , where , and , as well as other necessary parameters (such as , for the proposed BADM algorithm), is hand tuned for each algorithm in each experiment for the best ISNR; that is, for experiment (I); for (II); and for (III) and in (I), (II), and (III).

The obtained results of four metrics for the Cameraman, Mandril, and Lena images are listed in Tables 1, 2, and 3, respectively, and the deblurred images by BADM are shown in Figures 2, 3, and 4, respectively. And the evolutions of objective function and MSE by different algorithms in experiments (I), (II), and (III) (to avoid redundancy, only for the results of Cameraman images) are shown in Figures 5, 6, and 7, respectively. Moreover, the results of time, MSE, and ISNR over are shown in Figure 8. From above results, it is clear that the BADM outperforms the FISTA, SALSA, and SBM in terms of MSE and ISNR.

5. Conclusions

This paper has proposed a new bound alternative direction method (BADM) for the minimization problems in image deblurring. In order to solve the unconstrained optimization problem, the idea of BADM is to first bound the regularizer to obtain a bound unconstrained problem, which is then reformulated into a constrained one by variable splitting. The resulting constrained problem is further addressed by an augmented Lagrangian method, more specifically, the alternating direction method of multipliers (ADMM). Therefore, the BADM is an extension of the ADMM to solve the minimization problems. Experiments on a set of image deblurring problems have shown that the BADM for the synthesis formulation is favorably competitive with the state-of-the-art algorithms for the synthesis formulation.

In future work, the BADM will be applied to the analysis formulation and other applications such as in painting and magnetic resonance imaging.

Conflict of Interests

The author declares that there is no conflict of interests regarding to the publication of this paper.

Acknowledgments

This work was partially supported by the China Scholarship Council (CSC: 2010611017). Xiangrong Zeng would like to thank the anonymous reviewers who have helped to improve the quality of this paper.