Bound Alternative Direction Optimization for Image Deblurring

This paper proposes a new method, bound alternative direction method (BADM), to address the lp (p ∈ (0, 1)) minimization problems in image deblurring.The approach is to first obtain a bound unconstrained problem through bounding the lp 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 l1 minimization method and the alternating direction method of multiples (ADMM) such that it succeeds in extending the application of ADMM to lp minimization problems. The conducted experimental studies demonstrate the superiority of the proposed algorithm for the synthesis lp minimization over the state-of-the-art algorithms for the synthesis l1 minimization on image deblurring.


Introduction
The mission of image deblurring is to restore an original image x from noisy blurred observation y ∈ R  modeled as where x ∈ R  (stacking an  × -image into an ( = )-vector in lexicographic order), A ∈ R × is the matrix representation of a convolution operator, and n ∈ R  is Gaussian white noise.This imaging inverse problem has recently inspired a considerable amount of research ( [1][2][3][4][5][6][7][8][9][10][11] 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 x is represented as x = Ws, where W is a wavelet frame or a redundant dictionary and s are the sparse coefficients estimated usually via one sparsitypromoting regularizer, yielding where  is the positive regularization parameter and  is the regularizer, typically the ℓ 1 norm.The deblurred image is then obtained by x * = Ws * .In the analysis formulation, as opposed to (2), it minimizes the cost function with respect to x: where T is a sparsifying transform (such as wavelet or finite differences) and (Tx) analyzes the image x itself against the coefficients s that (s) in (2) works on.If T = W is a wavelet frame, then usually (Tx) = ‖Tx‖ 1 where ‖k‖ 1 = ∑  V  ; if T is a matrix representing finite differences at horizontal and vertical directions, then (Tx) 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 A 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 (A  A + I) −1 (where  is a positive parameter) can be efficiently computed (e.g., if A 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 [21][22][23][24], amongst which, the split Bregman method (SBM) [22] has been recently frequently applied to handle imaging inverse problems [25][26][27].
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 ℓ 1 minimization method and the ADMM and is able to extend the application of ADMM to the ℓ  minimization problems.
The proposed BADM has only O( log ) 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).

Terminology and Notation.
In this section, some useful elements of convex analysis will be given.Let H be a real Hilbert space equipped with the inner product ⟨⋅, ⋅⟩ and associated norm ‖ ⋅ ‖.Let  : H → [−∞, +∞] be a function and let Γ(x) 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  ∈ Γ(x), then Prox  is unique [44]; if (x) =   (x), the indicator function of a nonempty closed convex set , then Prox  (x) becomes the projection of k onto , and in this sense, ( 8) is therefore a generalization of projection operator; and if  is the ℓ 1 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.

Bound Alternative Direction Optimization
Consider a ℓ  regularized unconstrained model where  ∈ (0, 1) and  : R  → R 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  ∈ (0, 1), the nonconvex nature of ‖x‖   blocks the use of the proximity operator (see ( 8)), which is well defined only on the functions which belong to Γ(x).To overcome this problem, a bound optimization approach will be considered in this paper.[7].Since

Bound 𝐹(x) via the Majorizer Proposed in
where where Ω( x) = | x| −1 .It is easy to see that (x) ≤ (x, x) with equality for x = x, since ‖x‖   ≤ Φ(x, x) with equality for x = x.One benefit from the above bound is that the following lemma holds.Therefore, the proximity operator of Φ(x, x) given x is obtained by where Notice that soft(, +∞) = 0 for any .
From above, a closed-form solution of the proximity operator of Φ(x, x) 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 Φ(x, x), a nature way is to set x as the previous estimate x  and obtain the current estimate x +1 by computing the proximity operator of Φ(x, x  ); that is, where ^ is an iteratively temporary variable which differs in each algorithm.According to (13) with (14), if a component of x  becomes zero forever, then this component of x +1 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.

Proposed Bound Method.
A method is presented that bounds (x) through where where Ω  ( x) = (| x| + ) −1 and  is a small positive parameter.It is clear that (x) ≤   (x, x) with equality for x = x = 0 (where 0 is a vector composed of all zeros) or x = x with  = 0. Φ  (x, x) has the following property.(1) Set  = 0, and choose an arbitrary x 0 ,  0 > 0 and  > 1.

Iterative Bound Optimization (IBO). Now, it is ready to propose a framework of IBO as in Algorithm 1.
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 x  .Any accumulation point of the sequence {x  } generated by the IBO is a first-order stationary point of (19), which is guaranteed by the following theorem.
Theorem 3. Let the sequence {x  } be generated by above IBO and suppose that x * is an accumulation point of {x  }; then x * is a first-order stationary point of (19).
Proof.IBO is actually a specific case of the reweighted ℓ  ( = 1 or 2) method proposed in [45], such that this theorem is also a specific case of Theorem 3.1 in [45].

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 k and x fixing  and setting x to the previous estimate of x and then update , yielding a variant of ALM called bound ALM (BALM) (Algorithm 2).Notice that the terms added to   (x, x) in the definition of L(k, x; x, , , ) (see (22)) can be reformulated as a single quadratic term, leading to a variant of BALM (Algorithm 3).
In (Algorithm 3), d  corresponds to   / that these two parameters are used in BALM.It is usually difficult to simultaneously obtain k +1 and x +1 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 k and x, while keeping the other variable fixed,   yielding the following bound alternative direction method (BADM) (Algorithm 4).
In (Algorithm 4), Step 3 is equivalent to Prox / (x  + d  ), while Step 4 (see (13) and ( 9)) is equivalent to Moreover, since the objective function of ( 10) is nonconvex for  ∈ (0, 1), 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.

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   (18).Assume that A represents a (periodic) convolution and W is a normalized Parseval frame (i.e., WW  = I, where I is the identity matrix and possibly W  W ̸ = I).According to the Sherman-Morrison-Woodbury matrix inversion formula, Moreover, under the periodic boundary condition for s, AA  is block circulant, such that (AA  + I) −1 can be diagonalized by two-dimensional discrete Fourier transform (DFT).Let F = A  (AA  + I) −1 A; then F is equivalent to a filter in the Fourier domain and the cost of products by F using FFT algorithms is O( log ) [10].Thus, where u  = (AW)  y + (s  + d  ).As the computational complexity analysis in [10], the cost of computing k +1 is O( log ).Moreover, computing s +1 is the soft thresholding whose cost is O() and computing d +1 also has O() cost.Therefore, each iteration of the BADM for (4) has O( log ) cost.

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 512 × 512 pixels, corrupted by the following three benchmark cases [5,11] )) and the stopping criterion is chosen as ‖x +1 − x  ‖/‖x +1 ‖ ≤ , where  = 0.9, 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,  = 0.0075 for experiment (I); 0.09 for (II); and 0.25 for (III) and  = /10 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.

Conclusions
This paper has proposed a new bound alternative direction method (BADM) for the ℓ  ( ∈ (0, 1)) 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 ℓ  ( ∈ (0, 1)) 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 ℓ 1 formulation.
In future work, the BADM will be applied to the analysis ℓ  formulation and other applications such as in painting and magnetic resonance imaging.

Table 1 :
Results of deblurring the Cameraman images.

Table 2 :
Results of deblurring the Mandril images.

Table 3 :
Results of deblurring the Lena images.  (s, s  ) is the proposed majorizer of ‖s‖   at  iteration, which are inserted into the BADM yielding (Algorithm 5), where Step 4 is derived from