# Bound Alternative Direction Optimization for Image Deblurring

Authors: Xiangrong Zeng
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

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


---

## Body

## 1. Introduction

The mission of image deblurring is to restore an original imagex from noisy blurred observation y∈Rm modeled as
(1)y=Ax+n,
where x∈Rn (stacking an M×N-image into an (n=MN)-vector in lexicographic order), A∈Rm×n is the matrix representation of a convolution operator, and n∈Rm is Gaussian white noise. This imaging inverse problem has recently inspired a considerable amount of research ([1–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 imagex 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 sparsity-promoting regularizer, yielding
(2)s*=argmin⁡s12∥AWs-y∥22+λφ(s),
where λ is the positive regularization parameter and φ is the regularizer, typically the l1 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:
(3)x*=argmin⁡x12∥Ax-y∥22+λφ(Tx),
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 ∥v∥1=∑i‍vi; 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 (ATA+μ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 ofalternating direction method of multipliers (ADMM) [18, 19] which has a close relationship [20] with the Bregman iterations [21–24], amongst which, the split Bregman method (SBM) [22] has been recently frequently applied to handle imaging inverse problems [25–27].


### 1.2.lp 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 lp (p∈(0,1)) norm, is adopted since using the lp norm is able to find sparser solution than using the l1 norm which was demonstrated in many studies [28–32]. Thus, there comes the lp minimization problem of the synthesis formulation:
(4)min⁡s12∥AWs-y∥22+λ∥s∥pp
and of the analysis formulation:
(5)min⁡x12∥Ax-y∥22+λ∥Tx∥pp,
where ∥v∥pp=∑i‍|vi|p.Next, a brief survey on the literature oflp minimization is given. Recently, remarkable attention has been paid to the lp minimization problem:
(6)min⁡x∥x∥pp:s.t.Ax=y,
where p∈(0,1). Some sufficient conditions of (6) have been established in [29, 33–35]. Many efficient reweighted minimization methods address (6) through iteratively solving
(7)xk+1=argmin⁡xf(x)+λψk+1(x),
where xk+1 is the estimated signal at the (k+1)th iteration; f is the data-fitting term which preserves the consistency; for instance, f(x)=(1/2)∥Ax-y∥22 for the Gaussian noise; ψk+1(x) is the regularization term at the (k+1)th iteration, which is aiming to encourage the desirable properties of the target x such as sparsity. Most of reweighted minimization methods (see [2, 28, 30, 31, 34, 36, 37] and many references therein) focus on the reweighted l1 (IRL1) and the reweighted l2 (IRL2) minimization, and their corresponding regularizers are ψk+1(x)=∑i‍wik+1|xi| and ∑i‍wik+1|xi|2, respectively, where xi is the ith component of x and wik+1 is the weight of xi at the kth iteration and it is a function of the previous estimated component xik. In IRL1, usually, wik=(|xik|+ϵ)q-1, while wik=(|xik|2+ϵ)q/2-1 in IRL2, where q∈(0,1]. Notice that q is always set to 1 (e.g., IRL1: [28] and IRL2: [2]) or p (e.g., IRL1: [34, 36, 38] and IRL2: [30, 31, 38]), even 1/2 [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, thelp minimization problem will be used for image deblurring, but the resulting problem cannot be efficiently solved by above lp 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 lp 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 lp 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 lp 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 calledbound alternative direction method (BADM), which is equivalent to the combination of the reweighted l1 minimization method and the ADMM and is able to extend the application of ADMM to the lp minimization problems.The proposed BADM has onlyO(nlog⁡⁡n) cost in each iteration in solving the synthesis lp 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. LetH be a real Hilbert space equipped with the inner product 〈·,·〉 and associated norm ∥·∥. Let f: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 f is defined as
(8)Proxτf(v)=argmin⁡x(τf(x)+12∥x-v∥2),
where τ is positive. If f∈Γ(x), then Proxτf is unique [44]; if f(x)=ιC(x), the indicator function of a nonempty closed convex set C, then Proxτf(x) becomes the projection of v onto C, and in this sense, (8) is therefore a generalization of projection operator; and if f is the l1 norm, then (8) becomes a well-known soft thresholding:
(9)soft(v,τ)=sign⁡(v)⊙max⁡⁡{|v|-τ,0},
where ⊙ is the componentwise multiplication between two vectors of the same dimension and sign (·) is the sign function.


## 1.1. Formulations and Algorithms

In the synthesis formulation, the unknown imagex 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 sparsity-promoting regularizer, yielding
(2)s*=argmin⁡s12∥AWs-y∥22+λφ(s),
where λ is the positive regularization parameter and φ is the regularizer, typically the l1 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:
(3)x*=argmin⁡x12∥Ax-y∥22+λφ(Tx),
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 ∥v∥1=∑i‍vi; 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 (ATA+μ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 ofalternating direction method of multipliers (ADMM) [18, 19] which has a close relationship [20] with the Bregman iterations [21–24], amongst which, the split Bregman method (SBM) [22] has been recently frequently applied to handle imaging inverse problems [25–27].


## 1.2.lp 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 lp (p∈(0,1)) norm, is adopted since using the lp norm is able to find sparser solution than using the l1 norm which was demonstrated in many studies [28–32]. Thus, there comes the lp minimization problem of the synthesis formulation:
(4)min⁡s12∥AWs-y∥22+λ∥s∥pp
and of the analysis formulation:
(5)min⁡x12∥Ax-y∥22+λ∥Tx∥pp,
where ∥v∥pp=∑i‍|vi|p.Next, a brief survey on the literature oflp minimization is given. Recently, remarkable attention has been paid to the lp minimization problem:
(6)min⁡x∥x∥pp:s.t.Ax=y,
where p∈(0,1). Some sufficient conditions of (6) have been established in [29, 33–35]. Many efficient reweighted minimization methods address (6) through iteratively solving
(7)xk+1=argmin⁡xf(x)+λψk+1(x),
where xk+1 is the estimated signal at the (k+1)th iteration; f is the data-fitting term which preserves the consistency; for instance, f(x)=(1/2)∥Ax-y∥22 for the Gaussian noise; ψk+1(x) is the regularization term at the (k+1)th iteration, which is aiming to encourage the desirable properties of the target x such as sparsity. Most of reweighted minimization methods (see [2, 28, 30, 31, 34, 36, 37] and many references therein) focus on the reweighted l1 (IRL1) and the reweighted l2 (IRL2) minimization, and their corresponding regularizers are ψk+1(x)=∑i‍wik+1|xi| and ∑i‍wik+1|xi|2, respectively, where xi is the ith component of x and wik+1 is the weight of xi at the kth iteration and it is a function of the previous estimated component xik. In IRL1, usually, wik=(|xik|+ϵ)q-1, while wik=(|xik|2+ϵ)q/2-1 in IRL2, where q∈(0,1]. Notice that q is always set to 1 (e.g., IRL1: [28] and IRL2: [2]) or p (e.g., IRL1: [34, 36, 38] and IRL2: [30, 31, 38]), even 1/2 [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, thelp minimization problem will be used for image deblurring, but the resulting problem cannot be efficiently solved by above lp 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 lp 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 lp 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 lp 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 calledbound alternative direction method (BADM), which is equivalent to the combination of the reweighted l1 minimization method and the ADMM and is able to extend the application of ADMM to the lp minimization problems.The proposed BADM has onlyO(nlog⁡⁡n) cost in each iteration in solving the synthesis lp 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. LetH be a real Hilbert space equipped with the inner product 〈·,·〉 and associated norm ∥·∥. Let f: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 f is defined as
(8)Proxτf(v)=argmin⁡x(τf(x)+12∥x-v∥2),
where τ is positive. If f∈Γ(x), then Proxτf is unique [44]; if f(x)=ιC(x), the indicator function of a nonempty closed convex set C, then Proxτf(x) becomes the projection of v onto C, and in this sense, (8) is therefore a generalization of projection operator; and if f is the l1 norm, then (8) becomes a well-known soft thresholding:
(9)soft(v,τ)=sign⁡(v)⊙max⁡⁡{|v|-τ,0},
where ⊙ is the componentwise multiplication between two vectors of the same dimension and sign (·) is the sign function.


## 2. Bound Alternative Direction Optimization

Consider alp regularized unconstrained model
(10)min⁡xF(x)∶=f(x)+λ∥x∥pp,
where p∈(0,1) and f:Rn→R is a smooth function with L-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 p∈(0,1), the nonconvex nature of ∥x∥pp 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.


### 2.1. Bound Optimization

#### 2.1.1. BoundF(x) via the Majorizer Proposed in [7]

Since∥x∥pp(p∈(0,1))∉Γ(x), it is reasonable to bound F(x) through
(11)G^(x,x~)=f(x)+Φ^(x,x~),
where Φ^(x,x~)=λ∑i‍ϕ^(xi,x~i) with
(12)ϕ^(x,x~)={Ω^(x~)|x|+(1-p)|x~|p,ifx~≠0,+∞,ifx~=0,x≠0,0,ifx~=0,x=0,
where Ω^(x~)=p|x~|p-1. It is easy to see that F(x)≤G(x,x~) with equality for x=x~, since λ∥x∥pp≤Φ^(x,x~) with equality for x=x~. One benefit from the above bound is that the following lemma holds.Lemma 1.
Givenx~, Φ^(x,x~) belongs to Γ(x) with respect to x.Proof.
Givenx~i for i=1,…,n, if x~i≠0, then ϕ^(xi,x~i) is an affine function of |xi| and thus ϕ^(xi,x~i)∈Γ(xi), since |xi|∈Γ(xi); if x~i=0, then ϕ^(xi,x~i) equals 0 if xi=0 and +∞ otherwise, leading to ϕ^(xi,x~i)∈Γ(xi). Therefore, Φ^(x,x~)=λ∑i‍ϕ^(xi,x~i)∈Γ(x).Therefore, the proximity operator ofΦ^(x,x~) given x~ is obtained by
(13)ProxΦ^(x,x~)(v)=argmin⁡x(Φ^(x,x~)+12∥x-v∥2)=soft(v,ω^(x~)),
where ω^(x~)=[w1(x~1),…,wn(x~n)]T with
(14)wi(x~i)={λΩ^(x~i),ifx~i≠0,+∞,ifx~i=0.Notice that soft(x,+∞)=0 for any x.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 xk and obtain the current estimate xk+1 by computing the proximity operator of Φ^(x,xk); that is,
(15)xk+1=ProxΦ^(x,xk)(νk)=soft(νk,ω^(xk)),
where νk is an iteratively temporary variable which differs in each algorithm. According to (13) with (14), if a component of xk becomes zero forever, then this component of xk+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.


#### 2.1.2. Proposed Bound Method

A method is presented that boundsF(x) through
(16)Gϵ(x,x~)=f(x)+Φϵ(x,x~),
where Φϵ(x,x~)=λ∑i‍ϕϵ(xi,x~i) with
(17)ϕϵ(x,x~)=Ωϵ(x~)|x|+(1-p)|x~|p,
where Ωϵ(x~)=p(|x~|+ϵ)p-1 and ϵ is a small positive parameter. It is clear that F(x)≤Gϵ(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.Lemma 2.
Φ
ϵ
(
x
,
x
~
) belongs to Γ(x) with respect to x.Proof.
For anyx~i,i=1,…,n,  Ωϵ(x~i)>0; thus ϕϵ(xi,x~i) is an affine function of |xi| and thus ϕϵ(xi,x~i)∈Γ(xi). Therefore, Φϵ(x,x~)=λ∑i‍ϕϵ(xi,x~i)∈Γ(x).Therefore, a closed-form solution of the proximity operator ofΦϵ(x,x~) can be obtained
(18)ProxΦϵ(x,x~)(v)=argmin⁡x(Φϵ(x,x~)+12∥x-v∥2)=soft(v,ωϵ(x~)),
where ωϵ(x~)=[λΩϵ(x~1),…,λΩϵ(x~n)]T.


#### 2.1.3. Iterative Bound Optimization (IBO)

Now, it is ready to propose a framework of IBO as in Algorithm1.Algorithm 1:IBO.
(1)Set k=0, and choose an arbitrary x0, ϵ0>0 and γ>1.
(2)repeat
(3) xk+1=arg⁡min⁡x⁡Gϵk(x,xk)
(4) ϵk+1=ϵk/γ
(5) k←k+1
(6)until some stopping criterion is satisfied.A key observation is that the sequence{ϵk}, generated by the IBO, approaches to zero as k→+∞, such that (10) can be solved by the IBO which iteratively solves a sequence of problems:
(19)min⁡xFϵ(x)∶=f(x)+λ∑i=1n‍(|xik|+ϵk)p,
where xik is the ith component of xk. Any accumulation point of the sequence {xk} generated by the IBO is a first-order stationary point of (19), which is guaranteed by the following theorem.Theorem 3.
Let the sequence{xk} be generated by above IBO and suppose that x* is an accumulation point of {xk}; then x* is a first-order stationary point of (19).Proof.
IBO is actually a specific case of the reweightedlα (α=1or2) 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,(20)min⁡xGϵ(x,x~)∶=f(x)+Φϵ(x,x~).
Using the technique ofvariable splitting, (20) becomes a constrained optimization problem:
(21)min⁡v,xf(v)+Φϵ(x,x~)s.t.v=x.
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
(22)L(v,x;x~,α,β,ϵ)=f(v)+Φϵ(x,x~)+αT(v-x)+β2∥v-x∥22,
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 v and x fixing α and setting x~ to the previous estimate of x and then update α, yielding a variant of ALM calledbound ALM (BALM) (Algorithm 2).Algorithm 2:BALM.
(1)Set k=0, and choose β>0, γ>1, ϵ0>0, α0 and an arbitrary x0.
(2)repeat
(3) (vk+1,xk+1)∈arg⁡minv,x⁡L(v,x;xk,αk,β,ϵk)
(4) αk+1 = αk-β(vk+1-xk+1)
(5) ϵk+1 = ϵk/γ
(6) k←k+1
(7)until some stopping criterion is satisfiedNotice that the terms added toGϵ(x,x~) in the definition of L(v,x;x~,α,β,ϵ) (see (22)) can be reformulated as a single quadratic term, leading to a variant of BALM (Algorithm 3).Algorithm 3:Variant of BALM.
(1)Set k=0, and choose β>0, γ>1, ϵ0>0, d0 and an arbitrary x0.
(2)repeat
(3) (vk+1,xk+1)∈arg⁡minv,x⁡f⁡(v)+Φϵk(x,xk)+(β/2)∥v-x-dk∥22
(4) dk+1 = dk-(vk+1-xk+1)
(5) ϵk+1 = ϵk/γ
(6) k←k+1
(7)until some stopping criterion is satisfied.In (Algorithm3), dk corresponds to αk/β that these two parameters are used in BALM. It is usually difficult to simultaneously obtain vk+1 and xk+1 in Step 3. To overcome this difficulty, the technique ofnonlinear 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 v and x, while keeping the other variable fixed, yielding the followingbound alternative direction method (BADM) (Algorithm 4).Algorithm 4:BADM.
(1)Set k=0, and choose β>0,γ>1,ϵ0>0, d0 and an arbitrary x0.
(2) repeat
(3) vk+1=arg⁡⁡minv⁡f(v)+(β/2)∥v-xk-dk∥22⁡
(4) xk+1=arg⁡⁡minx⁡Φϵk(x,xk)+(β/2)∥vk+1-x-dk∥22
(5) dk+1=dk-(vk+1-xk+1)
(6) ϵk+1=ϵk/γ
(7) k←k+1
(8)until some stopping criterion is satisfied.In (Algorithm4), Step 3 is equivalent to Proxf/β(xk+dk), while Step 4 (see (13) and (9)) is equivalent to
(23)ProxΦϵk(x,xk)/β(vk+1-dk)=soft(vk+1-dk,ωϵk(xk)β).
Moreover, since the objective function of (10) is nonconvex for p∈(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.


## 2.1. Bound Optimization

### 2.1.1. BoundF(x) via the Majorizer Proposed in [7]

Since∥x∥pp(p∈(0,1))∉Γ(x), it is reasonable to bound F(x) through
(11)G^(x,x~)=f(x)+Φ^(x,x~),
where Φ^(x,x~)=λ∑i‍ϕ^(xi,x~i) with
(12)ϕ^(x,x~)={Ω^(x~)|x|+(1-p)|x~|p,ifx~≠0,+∞,ifx~=0,x≠0,0,ifx~=0,x=0,
where Ω^(x~)=p|x~|p-1. It is easy to see that F(x)≤G(x,x~) with equality for x=x~, since λ∥x∥pp≤Φ^(x,x~) with equality for x=x~. One benefit from the above bound is that the following lemma holds.Lemma 1.
Givenx~, Φ^(x,x~) belongs to Γ(x) with respect to x.Proof.
Givenx~i for i=1,…,n, if x~i≠0, then ϕ^(xi,x~i) is an affine function of |xi| and thus ϕ^(xi,x~i)∈Γ(xi), since |xi|∈Γ(xi); if x~i=0, then ϕ^(xi,x~i) equals 0 if xi=0 and +∞ otherwise, leading to ϕ^(xi,x~i)∈Γ(xi). Therefore, Φ^(x,x~)=λ∑i‍ϕ^(xi,x~i)∈Γ(x).Therefore, the proximity operator ofΦ^(x,x~) given x~ is obtained by
(13)ProxΦ^(x,x~)(v)=argmin⁡x(Φ^(x,x~)+12∥x-v∥2)=soft(v,ω^(x~)),
where ω^(x~)=[w1(x~1),…,wn(x~n)]T with
(14)wi(x~i)={λΩ^(x~i),ifx~i≠0,+∞,ifx~i=0.Notice that soft(x,+∞)=0 for any x.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 xk and obtain the current estimate xk+1 by computing the proximity operator of Φ^(x,xk); that is,
(15)xk+1=ProxΦ^(x,xk)(νk)=soft(νk,ω^(xk)),
where νk is an iteratively temporary variable which differs in each algorithm. According to (13) with (14), if a component of xk becomes zero forever, then this component of xk+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.


### 2.1.2. Proposed Bound Method

A method is presented that boundsF(x) through
(16)Gϵ(x,x~)=f(x)+Φϵ(x,x~),
where Φϵ(x,x~)=λ∑i‍ϕϵ(xi,x~i) with
(17)ϕϵ(x,x~)=Ωϵ(x~)|x|+(1-p)|x~|p,
where Ωϵ(x~)=p(|x~|+ϵ)p-1 and ϵ is a small positive parameter. It is clear that F(x)≤Gϵ(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.Lemma 2.
Φ
ϵ
(
x
,
x
~
) belongs to Γ(x) with respect to x.Proof.
For anyx~i,i=1,…,n,  Ωϵ(x~i)>0; thus ϕϵ(xi,x~i) is an affine function of |xi| and thus ϕϵ(xi,x~i)∈Γ(xi). Therefore, Φϵ(x,x~)=λ∑i‍ϕϵ(xi,x~i)∈Γ(x).Therefore, a closed-form solution of the proximity operator ofΦϵ(x,x~) can be obtained
(18)ProxΦϵ(x,x~)(v)=argmin⁡x(Φϵ(x,x~)+12∥x-v∥2)=soft(v,ωϵ(x~)),
where ωϵ(x~)=[λΩϵ(x~1),…,λΩϵ(x~n)]T.


### 2.1.3. Iterative Bound Optimization (IBO)

Now, it is ready to propose a framework of IBO as in Algorithm1.Algorithm 1:IBO.
(1)Set k=0, and choose an arbitrary x0, ϵ0>0 and γ>1.
(2)repeat
(3) xk+1=arg⁡min⁡x⁡Gϵk(x,xk)
(4) ϵk+1=ϵk/γ
(5) k←k+1
(6)until some stopping criterion is satisfied.A key observation is that the sequence{ϵk}, generated by the IBO, approaches to zero as k→+∞, such that (10) can be solved by the IBO which iteratively solves a sequence of problems:
(19)min⁡xFϵ(x)∶=f(x)+λ∑i=1n‍(|xik|+ϵk)p,
where xik is the ith component of xk. Any accumulation point of the sequence {xk} generated by the IBO is a first-order stationary point of (19), which is guaranteed by the following theorem.Theorem 3.
Let the sequence{xk} be generated by above IBO and suppose that x* is an accumulation point of {xk}; then x* is a first-order stationary point of (19).Proof.
IBO is actually a specific case of the reweightedlα (α=1or2) method proposed in [45], such that this theorem is also a specific case of Theorem  3.1 in [45].


## 2.1.1. BoundF(x) via the Majorizer Proposed in [7]

Since∥x∥pp(p∈(0,1))∉Γ(x), it is reasonable to bound F(x) through
(11)G^(x,x~)=f(x)+Φ^(x,x~),
where Φ^(x,x~)=λ∑i‍ϕ^(xi,x~i) with
(12)ϕ^(x,x~)={Ω^(x~)|x|+(1-p)|x~|p,ifx~≠0,+∞,ifx~=0,x≠0,0,ifx~=0,x=0,
where Ω^(x~)=p|x~|p-1. It is easy to see that F(x)≤G(x,x~) with equality for x=x~, since λ∥x∥pp≤Φ^(x,x~) with equality for x=x~. One benefit from the above bound is that the following lemma holds.Lemma 1.
Givenx~, Φ^(x,x~) belongs to Γ(x) with respect to x.Proof.
Givenx~i for i=1,…,n, if x~i≠0, then ϕ^(xi,x~i) is an affine function of |xi| and thus ϕ^(xi,x~i)∈Γ(xi), since |xi|∈Γ(xi); if x~i=0, then ϕ^(xi,x~i) equals 0 if xi=0 and +∞ otherwise, leading to ϕ^(xi,x~i)∈Γ(xi). Therefore, Φ^(x,x~)=λ∑i‍ϕ^(xi,x~i)∈Γ(x).Therefore, the proximity operator ofΦ^(x,x~) given x~ is obtained by
(13)ProxΦ^(x,x~)(v)=argmin⁡x(Φ^(x,x~)+12∥x-v∥2)=soft(v,ω^(x~)),
where ω^(x~)=[w1(x~1),…,wn(x~n)]T with
(14)wi(x~i)={λΩ^(x~i),ifx~i≠0,+∞,ifx~i=0.Notice that soft(x,+∞)=0 for any x.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 xk and obtain the current estimate xk+1 by computing the proximity operator of Φ^(x,xk); that is,
(15)xk+1=ProxΦ^(x,xk)(νk)=soft(νk,ω^(xk)),
where νk is an iteratively temporary variable which differs in each algorithm. According to (13) with (14), if a component of xk becomes zero forever, then this component of xk+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.


## 2.1.2. Proposed Bound Method

A method is presented that boundsF(x) through
(16)Gϵ(x,x~)=f(x)+Φϵ(x,x~),
where Φϵ(x,x~)=λ∑i‍ϕϵ(xi,x~i) with
(17)ϕϵ(x,x~)=Ωϵ(x~)|x|+(1-p)|x~|p,
where Ωϵ(x~)=p(|x~|+ϵ)p-1 and ϵ is a small positive parameter. It is clear that F(x)≤Gϵ(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.Lemma 2.
Φ
ϵ
(
x
,
x
~
) belongs to Γ(x) with respect to x.Proof.
For anyx~i,i=1,…,n,  Ωϵ(x~i)>0; thus ϕϵ(xi,x~i) is an affine function of |xi| and thus ϕϵ(xi,x~i)∈Γ(xi). Therefore, Φϵ(x,x~)=λ∑i‍ϕϵ(xi,x~i)∈Γ(x).Therefore, a closed-form solution of the proximity operator ofΦϵ(x,x~) can be obtained
(18)ProxΦϵ(x,x~)(v)=argmin⁡x(Φϵ(x,x~)+12∥x-v∥2)=soft(v,ωϵ(x~)),
where ωϵ(x~)=[λΩϵ(x~1),…,λΩϵ(x~n)]T.


## 2.1.3. Iterative Bound Optimization (IBO)

Now, it is ready to propose a framework of IBO as in Algorithm1.Algorithm 1:IBO.
(1)Set k=0, and choose an arbitrary x0, ϵ0>0 and γ>1.
(2)repeat
(3) xk+1=arg⁡min⁡x⁡Gϵk(x,xk)
(4) ϵk+1=ϵk/γ
(5) k←k+1
(6)until some stopping criterion is satisfied.A key observation is that the sequence{ϵk}, generated by the IBO, approaches to zero as k→+∞, such that (10) can be solved by the IBO which iteratively solves a sequence of problems:
(19)min⁡xFϵ(x)∶=f(x)+λ∑i=1n‍(|xik|+ϵk)p,
where xik is the ith component of xk. Any accumulation point of the sequence {xk} generated by the IBO is a first-order stationary point of (19), which is guaranteed by the following theorem.Theorem 3.
Let the sequence{xk} be generated by above IBO and suppose that x* is an accumulation point of {xk}; then x* is a first-order stationary point of (19).Proof.
IBO is actually a specific case of the reweightedlα (α=1or2) 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,(20)min⁡xGϵ(x,x~)∶=f(x)+Φϵ(x,x~).
Using the technique ofvariable splitting, (20) becomes a constrained optimization problem:
(21)min⁡v,xf(v)+Φϵ(x,x~)s.t.v=x.
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
(22)L(v,x;x~,α,β,ϵ)=f(v)+Φϵ(x,x~)+αT(v-x)+β2∥v-x∥22,
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 v and x fixing α and setting x~ to the previous estimate of x and then update α, yielding a variant of ALM calledbound ALM (BALM) (Algorithm 2).Algorithm 2:BALM.
(1)Set k=0, and choose β>0, γ>1, ϵ0>0, α0 and an arbitrary x0.
(2)repeat
(3) (vk+1,xk+1)∈arg⁡minv,x⁡L(v,x;xk,αk,β,ϵk)
(4) αk+1 = αk-β(vk+1-xk+1)
(5) ϵk+1 = ϵk/γ
(6) k←k+1
(7)until some stopping criterion is satisfiedNotice that the terms added toGϵ(x,x~) in the definition of L(v,x;x~,α,β,ϵ) (see (22)) can be reformulated as a single quadratic term, leading to a variant of BALM (Algorithm 3).Algorithm 3:Variant of BALM.
(1)Set k=0, and choose β>0, γ>1, ϵ0>0, d0 and an arbitrary x0.
(2)repeat
(3) (vk+1,xk+1)∈arg⁡minv,x⁡f⁡(v)+Φϵk(x,xk)+(β/2)∥v-x-dk∥22
(4) dk+1 = dk-(vk+1-xk+1)
(5) ϵk+1 = ϵk/γ
(6) k←k+1
(7)until some stopping criterion is satisfied.In (Algorithm3), dk corresponds to αk/β that these two parameters are used in BALM. It is usually difficult to simultaneously obtain vk+1 and xk+1 in Step 3. To overcome this difficulty, the technique ofnonlinear 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 v and x, while keeping the other variable fixed, yielding the followingbound alternative direction method (BADM) (Algorithm 4).Algorithm 4:BADM.
(1)Set k=0, and choose β>0,γ>1,ϵ0>0, d0 and an arbitrary x0.
(2) repeat
(3) vk+1=arg⁡⁡minv⁡f(v)+(β/2)∥v-xk-dk∥22⁡
(4) xk+1=arg⁡⁡minx⁡Φϵk(x,xk)+(β/2)∥vk+1-x-dk∥22
(5) dk+1=dk-(vk+1-xk+1)
(6) ϵk+1=ϵk/γ
(7) k←k+1
(8)until some stopping criterion is satisfied.In (Algorithm4), Step 3 is equivalent to Proxf/β(xk+dk), while Step 4 (see (13) and (9)) is equivalent to
(23)ProxΦϵk(x,xk)/β(vk+1-dk)=soft(vk+1-dk,ωϵk(xk)β).
Moreover, since the objective function of (10) is nonconvex for p∈(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.


## 3. Image Deblurring Using Synthesislp Formulation and BADM

In this section, the synthesislp formulation (see (4)) and the BADM are applied to image deblurring. Note that using the analysis lp formulation can be naturally extended and will be left as future work. Therefore, here f(s)=(1/2)∥AWs-y∥22, and Φϵk(s,sk) is the proposed majorizer of λ∥s∥pp  at k iteration, which are inserted into the BADM yielding (Algorithm 5), where Step 4 is derived from (18). Assume that A represents a (periodic) convolution and W is a normalized Parseval frame (i.e., WWT=I, where I is the identity matrix and possibly WTW≠I). According to the Sherman-Morrison-Woodbury matrix inversion formula,
(24)[(AW)TAW+βI]-1=1β[I-WTAT(AAT+βI)-1AW].
Moreover, under the periodic boundary condition for s, AAT is block circulant, such that (AAT+βI)-1 can be diagonalized by two-dimensional discrete Fourier transform (DFT). Let F=AT(AAT+βI)-1A; then F is equivalent to a filter in the Fourier domain and the cost of products by F using FFT algorithms is O(nlog⁡n) [10]. Thus,
(25)vk+1=1β(I-WTFW)uk,
where uk=(AW)Ty+β(sk+dk). As the computational complexity analysis in [10], the cost of computing vk+1 is O(nlog⁡n). Moreover, computing sk+1 is the soft thresholding whose cost is O(n) and computing dk+1 also has O(n) cost. Therefore, each iteration of the BADM for (4) has O(nlog⁡n) cost.Algorithm 5:BADM for (4).
(1)Set k=0, and choose β>0, γ>1, ϵ0>0, d0 and an arbitrary s0.
(2)repeat
(3) vk+1 = [(AW)TAW+βI]-1[(AW)Ty+β(sk+dk)]
(4) sk+1 = soft(vk+1-dk,ωϵk(sk)/β)
(5) dk+1 = dk-(vk+1-sk+1)
(6) ϵk+1 = ϵk/γ
(7) k←k+1
(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 512×512 pixels, corrupted by the following three benchmark cases [5, 11]: (I) uniform blur of size 9×9 and noise variance σ2=0.562 (termed UNI); (II) Gaussian blur kernel with σ2=8 (termed GAU); (III) the point spread function of the blur operator is hij=1/(1+i2+j2) for i,j=-7,…,7 and σ2=8 (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 x is the original image and yk and xk are the observed image and the estimated image at the k iteration, resp.) are employed: (a) consumed CPU time (Time); (b) number of iterations (Iterations); (c) mean square error (MSE=∥x-xk∥2/n); and (d) improvement in SNR (ISNR=10log⁡10(∑k‍∥x-yk∥2/∑k‍∥x-xk∥2)) and the stopping criterion is chosen as ∥xk+1-xk∥/∥xk+1∥≤η, where p=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).Test images.
(a)
Cameraman
(b)
Mandril
(c)
LenaThe obtained results of four metrics for the Cameraman, Mandril, and Lena images are listed in Tables1, 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 p 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.Table 1
Results of deblurring the Cameraman images.
Algorithm
Time (seconds)
Iterations
MSE
ISNR
UNI
GAU
PSF
UNI
GAU
PSF
UNI
GAU
PSF
UNI
GAU
PSF
FISTA
16.52
11.01
5.02
150
98
47
102.71
217.85
112.32
7.19
2.67
4.39
SALSA
0.80
2.17
1.72
4
8
6
90.40
214.41
109.25
7.80
2.75
4.52
SBM
0.95
2.53
1.50
4
8
6
90.40
214.36
109.25
7.80
2.75
4.52
BADM
1.44
3.28
1.55
4
9
6
85.42
212.29
105.18
8.05
2.78
4.66Table 2
Results of deblurring the Mandril images.
Algorithm
Time (seconds)
Iterations
MSE
ISNR
UNI
GAU
PSF
UNI
GAU
PSF
UNI
GAU
PSF
UNI
GAU
PSF
FISTA
56.82
49.11
17.22
110
94
33
119.23
288.45
142.12
5.50
1.42
2.19
SALSA
3.41
5.41
3.63
4
8
5
119.12
255.30
82.98
5.59
1.95
4.51
SBM
3.06
4.51
2.81
4
8
5
119.63
255.30
82.98
5.59
1.95
4.51
BADM
4.26
9.34
3.60
4
8
5
117.24
246.12
81.15
5.69
2.11
4.60Table 3
Results of deblurring the Lena images.
Algorithm
Time (seconds)
Iterations
MSE
ISNR
UNI
GAU
PSF
UNI
GAU
PSF
UNI
GAU
PSF
UNI
GAU
PSF
FISTA
53.13
40.8
21.41
98
81
42
43.89
82.84
53.90
6.01
2.82
2.17
SALSA
2.91
4.95
4.09
4
6
6
41.15
72.42
34.32
6.15
3.34
4.14
SBM
2.52
4.07
3.42
4
6
6
41.15
72.42
34.32
6.15
3.34
4.14
BADM
4.22
5.21
6.66
4
5
6
38.06
69.10
34.10
6.47
3.59
4.18Figure 2
Deblurred Cameraman images by BADM.Figure 3
Deblurred Mandril images by BADM.Figure 4
Deblurred Lena images by BADM.Evolutions in experiment (I) (Cameraman): (a) objective function and (b) MSE.
(a)
(b)Evolutions in experiment (II) (Cameraman): (a) objective function and (b) MSE.
(a)
(b)Evolutions in experiment (III) (Cameraman): (a) objective function and (b) MSE.
(a)
(b)Performance of BADM overp on deblurring the corrupted Cameraman images: (a) time; (b) MSE; and (c) ISNR.
(a)
(b)
(c)


## 5. Conclusions

This paper has proposed a newbound alternative direction method (BADM) for the lp(p∈(0,1)) minimization problems in image deblurring. In order to solve the unconstrained lp optimization problem, the idea of BADM is to first bound the lp 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, thealternating direction method of multipliers (ADMM). Therefore, the BADM is an extension of the ADMM to solve the lp(p∈(0,1)) minimization problems. Experiments on a set of image deblurring problems have shown that the BADM for the synthesis lp formulation is favorably competitive with the state-of-the-art algorithms for the synthesis l1 formulation.In future work, the BADM will be applied to the analysislp formulation and other applications such as in painting and magnetic resonance imaging.



---

*Source: 206926-2014-08-13.xml*