SIAM J. IMAGING SCIENCES Vol. 1, No. 3, pp. 248–272

c 2008 Society for Industrial and Applied Mathematics

A New Alternating Minimization Algorithm for Total Variation Image Reconstruction∗
Yilun Wang†, Junfeng Yang‡, Wotao Yin†, and Yin Zhang†
Abstract. We propose, analyze, and test an alternating minimization algorithm for recovering images from blurry and noisy observations with total variation (TV) regularization. This algorithm arises from a new half-quadratic model applicable to not only the anisotropic but also the isotropic forms of TV discretizations. The per-iteration computational complexity of the algorithm is three fast Fourier transforms. We establish strong convergence properties for the algorithm including finite convergence for some variables and relatively fast exponential (or q-linear in optimization terminology) convergence for the others. Furthermore, we propose a continuation scheme to accelerate the practical convergence of the algorithm. Extensive numerical results show that our algorithm performs favorably in comparison to several state-of-the-art algorithms. In particular, it runs orders of magnitude faster than the lagged diffusivity algorithm for TV-based deblurring. Some extensions of our algorithm are also discussed. Key words. half-quadratic, image deblurring, isotropic total variation, fast Fourier transform AMS subject classifications. 68U10, 65J22, 65K10, 65T50, 90C25 DOI. 10.1137/080724265

1. Introduction. In this paper, we propose a fast algorithm for reconstructing images from blurry and noisy observations. For simplicity, we assume that the underlying images have square domains, but all discussions can be equally applied to rectangle domains. Let 2 2 2 u0 ∈ Rn be an original n×n grayscale image, K ∈ Rn ×n represent a blurring (or convolution) 2 2 operator, ω ∈ Rn be additive noise, and f ∈ Rn be an observation which satisfies the relationship (1.1) f = Ku0 + ω.

Given f and K, the image u0 is recovered from the model n2 (1.2)
∗

min u i=1

Di u +

μ Ku − f 2

2 2,

Received by the editors July 3, 2007; accepted for publication (in revised form) May 16, 2008; published electronically August 20, 2008. http://www.siam.org/journals/siims/1-3/72426.html † Department of Computational and Applied Mathematics, Rice University, 6100 Main, Houston, TX, 77005 (yilun.wang@rice.edu, wotao.yin@rice.edu, yzhang@rice.edu). The work of the first and third authors was supported by NSF Career Grant DMS-0748839 and the latter author’s internal faculty research grant from the Dean of Engineering at Rice University. The fourth author’s work was supported in part by NSF grant DMS-0405831. ‡ Department of Mathematics, Nanjing University, 22 Hankou Road, Nanjing, Jiangsu Province, 210093, People’s Republic of China (jfyang2992@yahoo.com.cn). This author’s work was supported by the Chinese Scholarship Council during the author’s visit to Rice University. 248

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

249

where Di u ∈ R2 denotes the discrete gradient of u at pixel i and the sum Di u is the discrete total variation (TV) of u (the issue of boundary conditions will be discussed later). The origin of (1.2) and related results will be reviewed briefly in section 1.2. Model (1.2) is also widely referred to as the TV/L2 model. In (1.2), the TV is isotropic if the norm · in the summation is the 2-norm, and anisotropic if it is the 1-norm. We emphasize here that, unlike most previous work, our method applies to both isotropic and anisotropic TV. However, we will treat only the isotropic case, · = · 2 , in detail because the treatment for the anisotropic case is completely analogous. Our algorithm is derived from the well-known variable-splitting and penalty techniques in optimization; specifically, at each pixel an auxiliary variable wi ∈ R2 is introduced to transfer Di u out of the nondifferentiable term · 2 , and the difference between wi and Di u is penalized, yielding the following approximation model to (1.2): (1.3) min w,u i

wi

2

+

β 2

wi − Di u i 2 2

+

μ Ku − f 2

2 2,

with a sufficiently large penalty parameter β. This type of quadratic penalty approach can be traced back as early as Courant [13] in 1943. It is well known that the solution of (1.3) converges to that of (1.2) as β → ∞. Clearly, the objective function in (1.3) is convex in (u, w). The motivation for this formulation is that, while either one of the two variables u and w is fixed, minimizing the function with respect to the other has a closed-form formula with low computational complexity and high numerical stability. We will show that, for any fixed β, the algorithm of minimizing u and w alternately has finite convergence for some variables and q-linear convergence rates for the rest. Moreover, the overall convergence of this algorithm can be significantly accelerated by a theoretically well-justified continuation approach on the penalty parameter β. The objective function in (1.3) is quadratic in u and separable in w; therefore, problem (1.3) is half-quadratic according to Geman and Reynolds [17] and Geman and Yang [18]. Following the framework proposed in [17, 18], a number of half-quadratic models have been derived and analyzed. However, model (1.3) cannot be derived from (1.2) by the constructions given in [17, 18]; instead, a new approximate TV function must first be introduced before applying the construction of [18] (see section 2.3 for more details). The approximate TV model (1.3) and the resulting alternating minimization algorithm were first proposed in [42] without a convergence analysis. A similar split method has recently been proposed that uses Bregman iterations [20]. We have implemented the proposed algorithm in MATLAB and compared it with a C++ implementation of the lagged diffusivity algorithm [41], which is one of the fastest algorithms for solving (1.2). Numerical results presented in section 5 show that our algorithm is much faster than the lagged diffusivity algorithm especially when the blurring kernel is relatively large. Compared with a few other deblurring algorithms which solve different models, including ForWaRD (a Fourier and wavelet shrinkage method [32]) and MATLAB Image Processing Toolbox functions “deconvwnr” and “deconvreg,” our algorithm consistently generates higher quality images in comparable running times. In the rest of this section, we will introduce the notation, give a brief review of regularization approaches, and then summarize the contributions and organization of this paper.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

250

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG
2 2

1.1. Notation. Let D(j) ∈ Rn ×n , j = 1, 2, 3, . . . , be finite difference matrices and, in particular, let D(1) and D(2) be the first-order finite difference matrices in the horizontal and vertical directions, respectively. For j ≥ 3, D(j) represents higher-order finite difference matrices. For scalars αi , vectors vi , and matrices Ai of appropriate sizes, i = 1, 2, we let T T a = (α1 ; α2 ) (α1 , α2 )T , v = (v1 ; v2 ) (v1 , v2 )T , and A = (A1 ; A2 ) (AT , AT )T . As is used 1 2 2×n2 is a two-row matrix formed by stacking the ith rows of D (1) and D (2) . in (1.2), Di ∈ R Let ρ(T ) be the spectral radius of matrix T and let P(·) be a projection operator. From here on, the norm · refers to the 2-norm unless otherwise indicated. Additional notation will be introduced as the paper progresses. 1.2. Existing regularization approaches. There are different approaches based on statistics [14, 5], Fourier and/or wavelet transforms [25, 31], or variational analysis [37, 7, 11] for image deblurring. Among them the simplest is the maximum likelihood estimation, which solves the least squares problem minu Ku − f 2 and is also known as the inverse filter. However, the solution of this inverse filter, though it best matches the probabilistic behavior of the data, is often unacceptable because the image deblurring problems are usually severely ill-conditioned and the observed data usually contains noise. To stabilize the restoration process, some a priori knowledge about the unknown image is utilized through the addition of a regularization term, resulting in the model (1.4) min J(u) = Φreg (u) + u μ Ku − f 2

2

,

where Φreg (u) is the regularizer that enforces the a priori knowledge and the parameter μ is used to balance the two terms. Two classes of regularizers are well known. One is the Tikhonov-like class [40] including Φreg (u) = j D(j) u 2 , where D(j) ’s stand for some finite difference operators. In these cases, since the resulting objective functions J(u) are quadratic, it is relatively inexpensive to minimize them by solving linear systems of equations. However, since Tikhonov-like regularizers tend to make images overly smooth, they often fail to adequately preserve important image attributes such as sharp edges. Another well-known class of regularizers are based on TV, which was first proposed for image denoising by Rudin, Osher and Fatemi in [38], and then extended to image deblurring in [37]. In comparison to Tikhonov-like regularizers, TV regularizers can better preserve sharp edges or object boundaries that are usually the most important features to recover. The TV/L2 model (1.2) has been widely studied (see, for example, [7, 8] and the references therein). The superiority of TV over Tikhonov-like regularization was analyzed in [1, 15] for recovering images containing piecewise smooth objects. However, due to the nondifferentiability and nonlinearity of the TV functions, the TV/L2 model (1.2) is computationally more difficult to solve. Despite numerous efforts over the years, algorithms for solving the isotropic TV/L2 model are still much slower than those for solving Tikhonov-like regularization models. In this work, we aim to overcome this difficulty. The most important structure of the TV/L2 model is that both the blurring operator and the finite-difference operators can be treated as discrete convolutions under suitable boundary conditions. Determining how to best exploit this structure is crucial for achieving high computational efficiency.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

251

1.3. Contributions. This paper has three major contributions: the derivation of a new half-quadratic model and the accompanying algorithm, the analysis of the algorithm’s convergence properties, and the evaluation of its practical performance. 1. To the best of our knowledge, (1.3) is the first half-quadratic method for the isotropic TV/L2 model whose leading computational requirement is three fast Fourier transforms (FFTs) per iteration. Previously, this level of complexity was achievable only for the anisotropic case. Our extension not only enables better image quality, but, more importantly, can be generalized to color (or multichannel) image reconstruction with TV regularization. 2. We established convergence results for our algorithm that are stronger than those proved for existing half-quadratic methods. We established finite convergence for some auxiliary variables and strong q-linear convergence rates for the other variables. The q-linear rates are stronger than usual because they depend on the spectral radii of submatrices rather than of the whole matrix as in usual cases. 3. Based on our convergence results, we propose a continuation scheme on the penalty parameter with a warm start. Our extensive numerical results have confirmed that the proposed algorithm can be orders of magnitude faster than the lagged diffusivity method—one of the state-of-the-art methods for solving the TV/L2 model (1.2). In addition, we can extend our method to other problems involving TV regularization such as compressive sensing. These extensions will be discussed in section 6. 1.4. Organization. The rest of this paper is organized as follows. In section 2, we propose our alternating minimization algorithm, study optimality conditions, and describe its relations with previous works. Convergence properties of the proposed algorithm are analyzed in section 3. In section 4, we propose a continuation scheme, demonstrate its efficiency, and describe the implementation of our complete algorithm. Numerical results in comparison with the lagged diffusivity (LD) and some other typical image restoration methods are presented in section 5. In section 6, we discuss how our method can be extended to more general inverse problems, including some compressed sensing formulations. Finally, concluding remarks are given in section 7. 2. Basic algorithm, optimality, and related works. We first introduce two auxiliary 2 vectors w1 , w2 ∈ Rn to approximate D(1) u and D(2) u, respectively, where we recall that 2 2 D(1) , D(2) ∈ Rn ×n represent the two first-order forward finite difference operators with ap2 propriate boundary conditions. We denote w = (w1 ; w2 ) ∈ R2n and D = (D(1) ; D(2) ) ∈ 2 ×n2 , where the latter is the total finite difference operator. For i = 1, . . . , n2 , we let R2n wi = ((w1 )i ; (w2 )i ) ∈ R2 , which is an approximation of Di u = (D(1) u)i ; (D(2) u)i ∈ R2 . To keep wi close to Di u, we apply a quadratic penalty term, which is the second term in (1.3). 2.1. An alternating minimization algorithm. It is easy to minimize the objective function in (1.3) with respect to either u or w. For a fixed u, the first two terms in (1.3) are separable with respect to wi , so minimizing (1.3) for w is equivalent to solving, for i = 1, 2, . . . , n2 , (2.1) min wi wi +

β wi − Di u 2

2

,

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

252

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

for which the unique minimizer is given by the following two-dimensional (2D) shrinkage formula: Di u 1 , i = 1, . . . , n2 , (2.2) wi = max Di u − , 0 β Di u where the convention 0 · (0/0) = 0 is followed. For the anisotropic case · = · 1 , each component of wi is given by the simpler one-dimensional (1D) shrinkage formula: wi = max{|Di u| − 1/β, 0} sgn(Di u), where all operations are done componentwise. On the other hand, for a fixed w, (1.3) is quadratic in u, and the minimizer u is given by the normal equations
T Di Di + i

μ T K K β

u= i T Di wi +

μ T K f β

or, equivalently, (2.3) (D(1) )T D(1) + (D(2) )T D(2) + μ T μ K K u = (D(1) )T w1 + (D(2) )T w2 + K T f. β β

Under the periodic boundary condition for u, (D(1) )T D(1) , (D(2) )T D(2) , and K T K are all block circulant (see [33, 21], for example). Therefore, the Hessian matrix on the left-hand side of (2.3) can be diagonalized by 2D discrete Fourier transform F. Using the convolution theorem of Fourier transforms, we can write (2.4) u = F −1 F(D(1) )∗ ◦ F(w1 ) + F(D(2) )∗ ◦ F(w2 ) + (μ/β)F(K)∗ ◦ F(f ) F(D(1) )∗ ◦ F(D(1) ) + F(D(2) )∗ ◦ F(D(2) ) + (μ/β)F(K)∗ ◦ F(K) ,

where “∗” denotes complex conjugacy, “◦” denotes componentwise multiplication, and the division is componentwise as well. With a slight abuse of notation, we have used F(K) for the Fourier transform of the function represented by K in the convolution Ku (and similarly for D(1) and D(2) ). Since all quantities but w1 and w2 are constant, computing u from (2.4) involves two FFTs and one inverse FFT, once the constant quantities are computed. Alternatively, under the Neumann boundary condition and assuming that the blurring kernel K is symmetric, the left-hand side of (2.3) is a block Toeplitz-plus-Hankel matrix (see [33]), so (2.3) can be diagonalized by 2D discrete cosine transforms (DCTs), and solving (2.3) requires three DCT calls. In our numerical experiments, we used the periodic boundary condition and FFTs. Since minimizing the objective function in (1.3) with respect to either w or u is computationally inexpensive, we solve (1.3) for a fixed β by an alternating minimization scheme given below. Algorithm 1. Input f , K, μ > 0, and β > 0. Initialize u = f . While “not converged,” Do (1) Compute w according to (2.2) for fixed u. (2) Compute u according to (2.4) for fixed w. End Do More details about this algorithm will be presented later. Next, we derive optimality conditions for (1.3), based on which we will specify the stopping criterion (2.10) for Algorithm 1.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

253

2.2. Optimality conditions. To study the optimality conditions of (1.3) for a fixed β, we need the following proposition. Ax is Proposition 2.1. For any A ∈ Rp×n , the subdifferential of f (x) (2.5) ∂f (x) = {AT Ax/ Ax } AT h : h ≤ 1, h ∈ Rp if Ax = 0, otherwise.

Proof. By the definition of subdifferential for a convex function, (2.6) ∂f (x) {g ∈ Rn : Ay − Ax ≥ g T (y − x) ∀ y ∈ Rn }.

If Ax = 0, then f is differentiable at x, and clearly (2.5) holds. Now assume Ax = 0. In this case, after setting y = x + d, we have ∂f (x) = {g ∈ Rn : Ad ≥ g T d ∀ d ∈ Rn }. We will show that ∂f (x) ≡ S AT h : h ≤ 1, h ∈ Rp . First, for any AT h ∈ S, (AT h)T d ≤ Ad by the Cauchy–Schwarz inequality, implying S ⊂ ∂f (x). Next, we show ∂f (x) ⊂ S by contradiction. Suppose that g ∈ ∂f (x), but g ∈ S. Since S is closed and convex, by the well/ known separation theorem, there exist d ∈ Rn and τ ∈ R such that the hyperplane dT x = τ separates g and S so that dT g > τ > dT p for all p ∈ S. It follows that dT g > τ ≥ sup{dT AT h : h ≤ 1, h ∈ Rp } = Ad , contradicting the assumption that g ∈ ∂f (x). Therefore, we have ∂f (x) = S, and the result is proved. Since the objective function in (1.3) is convex, a pair (w, u) solves (1.3) if and only if the subdifferential of the objective function at (w, u) contains the origin, which gives the following optimality conditions in light of Proposition 2.1 with A = I: (2.7) (2.8) wi / wi + β(wi − Di u) = 0, β Di u ≤ 1, i ∈ I1 i ∈ I2 {i : wi = 0}, {i : wi = 0},

βDT (Du − w) + μK T (Ku − f ) = 0.

Our stopping criterion for Algorithm 1 is based on the optimality conditions (2.7) and (2.8). Let ⎧ i ∈ I1 , ⎨ r1 (i) (wi / wi )/β + wi − Di u, Di u − 1/β, i ∈ I2 , r2 (i) ⎩ T (Du − w) + μK T (Ku − f ), r3 βD

(2.9)

where I1 and I2 are defined as in (2.7). Algorithm 1 is terminated once (2.10) Res max max{ r1 (i) }, max{r2 (i)}, r3 i∈I1 i∈I2 ∞

≤

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

254

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

is met, where Res measures the total residual and > 0 is a prescribed tolerance. In (2.9), condition (2.7) is scaled by a factor of 1/β, but (2.8) is not, because in practice the latter can be met quickly even without this scaling. Combining (2.7) and (2.8) to eliminate the w-variable, we can derive (2.11) i∈I1 T Di

Di u + Di u

T Di hi + μK T (Ku − f ) = 0, i∈I2

∗ where hi βDi u satisfies hi ≤ 1. Let u∗ be any solution of (1.2). Define I1 = {i, Di u∗ = 0} ∗ = {1, . . . , n2 } \ I ∗ . Then, from Proposition 2.1, there exist {h∗ ∈ R2 : h∗ ≤ 1, i ∈ I ∗ } and I2 1 2 i i such that

(2.12)
∗ i∈I1

T Di

Di u∗ + Di u∗

T Di h∗ + μK T (Ku∗ − f ) = 0. i
∗ i∈I2

Equation (2.11) differs from (2.12) only in the index sets over which summations are taken. ∗ As β increases, I1 should approach I1 . In section 3, we will study the convergence properties of Algorithm 1. 2.3. Related works. While (1.3) is an instance of the classical quadratic penalty method, it is also a half-quadratic model according to Geman and Reyonlds [17] and Geman and Yang [18]. However, this model cannot be directly derived from the original frameworks in [17, 18] without introducing a new approximation to the isotropic TV function. Consider the following general framework of recovering an image u from its corrupted measurements f : (2.13)
2

min u i

T φ(gi u) +

μ Ku − f 2

2

,

T where gi ∈ Rn is a local finite difference operator, φ(gi ·) is convex and edge-preserving, and K is a convolution operator. Instead of solving (2.13) directly, the authors of [17] and [18] propose solving the equivalent problem

(2.14)

min u,b i

T Q(gi u, bi ) + ψ(bi ) +

μ Ku − f 2

2

,

where Q(t, s) and ψ(s) are chosen such that Q(t, s) is quadratic in t and (2.15) φ(t) = min(Q(t, s) + ψ(s)) ∀t ∈ R. s∈R The name “half-quadratic” comes from the fact that the objective function in (2.14) is quadratic in u but not in b. As such, under certain conditions such as convexity (2.14) can be solved by minimizing with respect to u and b alternately. Computationally, it is important that the constructed Q and ψ allow fast solutions for u and b, respectively. For this purpose, two forms of half-quadratic formulations have been widely studied: the multiplicative form [17] in which Q(t, s) = 1 t2 s and the additive form [18] in which Q(t, s) = 1 (t − s)2 . For 2 2

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

255

comparisons of these two forms of methods in both theory and practice, see [24, 2, 34, 12] and the references therein. Unfortunately, the constructions given in [17] and [18] cannot be directly applied to either the anisotropic or the isotropic TV because some technical conditions fail to hold. Over the years, numerous half-quadratic models have been proposed to approximate the TV functions. Some of these models are of the additive form [24, 34], while others are of the multiplicative form [7, 24]. However, in the half-quadratic models of the multiplicative form, the partial Hessians with respect to u depend on b and thus vary from one iteration to another. In addition, such partial Hessians do not have a block circulant structure. As a result, halfquadratics of the multiplicative form cannot be efficiently minimized by FFTs. Our half-quadratic model in (1.3) approximates the isotropic TV function and is of the additive form. The Hessian of the quadratic is constant and has a block circulant structure (under appropriate boundary conditions). The only other half-quadratic model for TV approximation that shares these attributes is the one by Nikolova and Ng [34], but it is for the anisotropic rather than the isotropic TV function. In [34], the authors use the Huber function
1 2 2 t

(2.16)

φ(t) =

|t| −

2

if |t| ≤ , otherwise,

where 0 < 1, to approximate φ(t) = |t| in (2.13) and apply the additive form halfquadratic construction. The resulting problem (2.14) is quadratic in u and separable in b, and the Hessian of the quadratic is constant and has a block circulant structure, enabling the use of fast transforms. However, [34] does not include a convergence analysis or extensive numerical experiments on the half-quadratic model for this anisotropic TV approximation. In [34], two approaches are proposed to approximate the isotropic TV function. One is based on the multiplicative construction; the other employs the additive construction but does not yield a half-quadratic model. Neither approach leads to a fast alternating minimization algorithm. Thus far all of the edge-preserving functions φ used in half-quadratic constructions that we are aware of take only scalar variables. In order to approximate the isotropic TV function, one should allow φ to take vector variables; then our model (1.3) can be derived from a similar additive construction, which we describe below. First, (1.2) is approximated by (2.17) min u i

φ(Di u) +

μ Ku − f 2

2

,

where φ : R2 → R is a 2D Huber-like function defined as φ(t) = β 2

t

2 1 2β

t −

if t ≤ 1/β, otherwise

for β 0. Then a similar additive construction in two dimensions, instead of in one dimension, gives Q(t, s) = β t − s 2 and ψ(s) = s , where s, t ∈ R2 , which satisfy the 2D version of 2 (2.15). Clearly, s and β t − s 2 are nothing but the first and second terms in our model 2 (1.3) for s = wi and t = Di u. Therefore, (1.3) turns out to be a long overdue extension to Nikolova and Ng’s half-quadratic model [34].

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

256

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

Since the work of Geman and Reynolds [17] and Geman and Yang [18], the convergence of alternating minimization algorithms for various half-quadratic models has been extensively studied, including that for multiplicative formulations in [24, 2, 34, 12] and for additive formulations in [24, 2, 34]. However, these analyses were done under some strict convexity assumptions on objective functions; see [24, 34, 2, 12]. Moreover, the convergence rate results obtained were r-linear at best; see [34, 2], for example. Unlike these previous results, our convergence results, presented in the next section, require neither nonsingularity of the operator K nor strict convexity of the objective function as a whole (notice that φ = · 2 is not strictly convex). Moreover, we establish a finite convergence property for some variables and strong q-linear convergence for the rest. 3. Convergence analysis. The convergence of the quadratic penalty method as the penalty parameter goes to infinity is well known (see Theorem 17.1 in [35], for example). That is, as β → ∞, the solution of (1.3) converges to that of (1.2). However, it is generally difficult to determine theoretically how large a β value must be to attain a given accuracy. We will later choose β based on numerical experiments. In this section, we analyze the convergence properties of Algorithm 1 for a fixed β > 0. First, we prove that the sequence {(wk , uk )} generated by Algorithm 1 from any initial point converges to a solution of (1.3). Then, using a finite convergence property for {wk }, we establish q-linear convergence rates for Algorithm 1. In what follows, we omit β for notational simplicity. For a ∈ R2 , the 2D shrinkage operator s : R2 → R2 is defined as (3.1) s(a) max a − 1 ,0 β a , a

where the convention 0 · (0/0) = 0 is followed. It is easy to see that s(a) = a − P(a), where P(·) PB (·) : R2 → R2 is the projection onto the closed disc B For vectors u, v ∈ RN , N ≥ 1, we define S(u; v) : R2N → R2N by S(u; v) (s(a1 ); . . . ; s(aN )) , where ai = ui vi ; {a ∈ R2 | a ≤ 1/β}.

i.e., S applies 2D shrinkage to each pair (ui , vi ) ∈ R2 , for i = 1, 2, . . . , N . The first convergence result is Theorem 3.4, and we take a few steps to prove it. The first step is to prove the nonexpansiveness of the shrinkage operator s. The following proposition is known to hold for projections P onto any closed convex set. We include a proof for the sake of completeness. Proposition 3.1. For any a, b ∈ R2 , it holds that s(a) − s(b)
2

≤ a−b

2

− P(a) − P(b) 2 .

Furthermore, if s(a) − s(b) = a − b , then s(a) − s(b) = a − b.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

257

Proof. Since B is convex, P(·) satisfies (a − P(a))T (P(a) − P(b)) ≥ 0 ∀ a, b ∈ R2 .

Exchanging a and b gives us −(b − P(b))T (P(a) − P(b)) ≥ 0. Adding the two yields (3.2) (a − b)T (P(a) − P(b)) ≥ P(a) − P(b) 2 .

Therefore, from s(a) = a − P(a), we get s(a) − s(b)
2

= a − b − (P(a) − P(b)) = a−b ≤ a−b
2 2

2 2

− 2(a − b)T (P(a) − P(b)) + P(a) − P(b) − P(a) − P(b) 2 ,

where the last inequality follows from (3.2). Moreover, if s(a) − s(b) = a − b , then P(a) = P(b) and s(a) − s(b) = a − b − (P(a) − P(b)) = a − b. We will make use of the following assumption in our convergence analysis. It is a mild assumption that has been used in most, if not all, previous analysis of a similar type. Assumption 1. N (K) ∩ N (D) = {0}, where N (·) represents the null space of a matrix. The following two symmetric positive definite matrices will be used in our analysis: (3.3) M = DT D + μ T K K and T = DM −1 DT . β

Assumption 1 ensures the nonsingularity of M ; thus T is well defined. It is worth noting that 2 2 ρ(T ) ≤ 1. We also define a linear operator h : R2n → R2n as h(w) = (h(1) (w); h(2) (w)), where μ h(j) (w) = D(j) M −1 DT w + K T f , j = 1, 2. β Using the definitions of S and h, we can rewrite the two iterative steps in Algorithm 1 as (3.4) (3.5) wk+1 = S(D(1) uk ; D(2) uk ) = S ◦ h(wk ), μ uk+1 = M −1 DT wk+1 + K T f . β

Since the objective function in (1.3) is convex, bounded below, and coercive (i.e., its value goes to infinity as (w, u) → ∞), (1.3) has at least one minimizer pair (w∗ , u∗ ) that cannot be decreased by our alternating minimization scheme and thus must satisfy (3.6) (3.7) w∗ = S(D(1) u∗ ; D(2) u∗ ) = S ◦ h(w∗ ), μ u∗ = M −1 DT w∗ + K T f . β

Particularly, (3.6) means that w∗ is a fixed point of S ◦ h. We need to show that h is nonexpansive. 2 Proposition 3.2. For any w = w in R2n , it holds that ˜ h(w) − h(w) ≤ w − w , ˜ ˜

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

258

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

and the equality holds if and only if h(w) − h(w) = w − w. ˜ ˜ Proof. From the definitions of M and T , it can be shown that ρ(T ) = ρ(DM −1 DT ) ≤ 1 and h(w) − h(w) = T (w − w) ≤ ρ(T ) w − w ≤ w − w . ˜ ˜ ˜ ˜ Let T = QT ΛQ be the eigendecomposition of T , where Q is an orthogonal matrix and Λ is a diagonal matrix with elements 0 ≤ λi ≤ 1, for all i. If h(w) − h(w) = w − w , namely, ˜ ˜ QT ΛQ(w − w) = w − w , then ΛQ(w − w) = Q(w − w) . Since Λ is diagonal and ˜ ˜ ˜ ˜ 0 ≤ λi ≤ 1, ΛQ(w − w) = Q(w − w) holds true. Multiplying both sides by QT completes the ˜ ˜ proof. The following lemma gives a useful property for fixed points of the operator S ◦ h. Lemma 3.3. Let w be any fixed point of S ◦ h. For any w, we have S ◦ h(w) − S ◦ h(w) < ˆ ˆ w − w unless w is a fixed point of S ◦ h. ˆ Proof. From Propositions 3.1 and 3.2 and the definition of S, it holds that S ◦ h(w) − S ◦ h(w) ≤ w − w . If equality holds, then, again from Propositions 3.1 and 3.2, we get ˆ ˆ S ◦ h(w) − S ◦ h(w) = h(w) − h(w) = w − w. ˆ ˆ ˆ Since w = S ◦ h(w), we get w = S ◦ h(w). ˆ ˆ Now we are ready to prove the convergence of Algorithm 1. Theorem 3.4. For any fixed β > 0, the sequence {(wk , uk )} generated by Algorithm 1 from any starting point (w0 , u0 ) converges to a solution (w∗ , u∗ ) of (1.3). Proof. We prove the convergence of {wk } in three steps. First, from (3.4) and the nonexpansiveness of S and h, it is easy to show that the sequence {wk } lies in a compact region and thus has at least one limit point, say w∗ = limj→∞ wkj . Letting w be any fixed point of ˆ S ◦ h, i.e., w = S ◦ h(w), we get ˆ ˆ ˆ ˆ ˆ wk − w = S ◦ h(wk−1 ) − S ◦ h(w) ≤ wk−1 − w . Therefore, the following limit exists: (3.8) k→∞ ˆ ˆ ˆ lim wk − w = lim wkj − w = w∗ − w , j→∞ which means that all limit points of {wk }, if more than one, have an equal distance to w. By ˆ the continuity of S ◦ h, we have S ◦ h(w∗ ) = lim S ◦ h(wkj ) = lim wkj +1 . j→∞ j→∞

Hence, S ◦ h(w∗ ) is also a limit point of {wk } that must have the same distance to w as w∗ ˆ does; i.e., w∗ − w = S ◦ h(w∗ ) − w = S ◦ h(w∗ ) − S ◦ h(w) . ˆ ˆ ˆ It follows from Lemma 3.3 that w∗ = S ◦h(w∗ ). Since w is any fixed point of S ◦h, by replacing ˆ w with w∗ in (3.8), we establish the convergence of {wk }: limk→∞ wk = w∗ . The convergence ˆ of {uk } to some u∗ follows directly from (3.5). Therefore, (w∗ , u∗ ) satisfies (3.6)–(3.7) and is a solution of (1.3).

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

259

We note that Theorem 3.4 does not require the objective function in (1.3) to be strictly convex in (u, w) (even though it is so with respect to each individual variable). As such, (1.3) may have multiple solutions, and the limit (w∗ , u∗ ) of the sequence generated by Algorithm 1 depends on starting point (w0 , u0 ). Next we develop a finite convergence property for the auxiliary variable w. Let hi (w) = (hi (w); hi (w)) ∈ R2 ,
(1) (2)

i = 1, . . . , n2 ;

namely, hi (w) is a vector formed by stacking the ith components of h(1) (w) and h(2) (w). We will make use of the following two index sets: (3.9) L= i : Di u∗ = hi (w∗ ) < 1 β and E = {1, . . . , n2 } \ L.

Theorem 3.5 (finite convergence). The sequence {(wk , uk )} generated by Algorithm 1 from k ∗ any starting point (w0 , u0 ) satisfies wi = wi = 0, for all i ∈ L, for all but a finite number of 0 − w ∗ 2 /ω 2 , where iterations that does not exceed w (3.10) ω min i∈L 1 − hi (w∗ ) β

> 0.

Proof. By the nonexpansiveness of s(·), for all i, (3.11) k+1 ∗ wi − wi = s ◦ hi (wk ) − s ◦ hi (w∗ ) ≤ hi (wk ) − hi (w∗ ) .

k+1 Suppose that at iteration k there exists at least one index i ∈ L such that wi = s ◦ hi (wk ) = ∗ ) < 1/β, h (w k ) ≥ 1/β, and w∗ = s ◦ h (w ∗ ) = 0. Therefore, 0. Then, hi (w i i i

(3.12)

k+1 ∗ wi − wi

2

= s ◦ hi (wk )

2

= ( hi (wk ) − 1/β)2
2 2

≤ { hi (wk ) − hi (w∗ ) − (1/β − hi (w∗ ) )}2 ≤ hi (wk ) − hi (w∗ ) ≤ hi (wk ) − hi (w∗ ) − (1/β − hi (w∗ ) )2 − ω2,

where the first inequality is a triangular inequality, the second follows from the fact that hi (wk ) − hi (w∗ ) ≥ 1/β − hi (w∗ ) > 0, and the last uses the definition of ω in (3.10). Combining (3.11), (3.12), and the nonexpansiveness of h(·), we obtain wk+1 − w∗
2

= i k+1 ∗ wi − wi

2

≤ i hi (wk ) − hi (w∗ )
2

2

− ω2

= h(wk ) − h(w∗ )

2

− ω 2 ≤ wk − w∗

− ω2.

k+1 Therefore, the number of iterations k with wi = 0 does not exceed w0 − w∗ 2 /ω 2 . k = w∗ for i ∈ L, we next show the linear convergence of Given the finite convergence of wi i k and wk for i ∈ E. Denote w = ((w ) ; (w ) ), where (w ) and (w ) are the subvectors u 1 E 2 E 1 E 2 E E i of w1 and w2 , respectively, with components (w1 )i , (w2 )i , i ∈ E. Theorem 3.6 (q-linear convergence). Let M and T be defined as in (3.3) and let TEE = [Ti,j ]i,j∈E∪(n2 +E) . Under Assumption 1, the sequence {(wk , uk )} generated by Algorithm 1 satisfies

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

260

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

k+1 ∗ k ∗ 1. wE − wE ≤ ρ((T 2 )EE ) wE − wE ; 2. D(uk+1 − u∗ ) ≤ ρ((T 2 )EE ) D(uk − u∗ ) ; 3. uk+1 − u∗ M ≤ ρ(TEE ) uk − u∗ M for all k sufficiently large. Proof. From (3.4)–(3.7) and the nonexpansiveness of S, we get

(3.13) and (3.14) wk+1 − w∗
2

uk+1 − u∗ = M −1 DT (wk+1 − w∗ )

= S(D(1) uk ; D(2) uk ) − S(D(1) u∗ ; D(2) u∗ )

2

≤ D(uk − u∗ ) 2 .

Combining the recursion (3.13), (3.14) and the definition of T , it holds that wk+1 − w∗
2

≤ T (wk − w∗ ) 2 .

Since we are interested in the asymptotic behavior of Algorithm 1, without loss of generality, k k k ∗ hereafter we assume that wL ((w1 )L ; (w2 )L ) = wL = (0; 0). Therefore, the above inequality becomes k+1 ∗ wE − wE 2 k ∗ k ∗ k ∗ ≤ (wE − wE )T (T 2 )EE (wE − wE ) ≤ ρ((T 2 )EE ) wE − wE 2

.

k ∗ Multiplying D on both sides of (3.13), from wL = wL = 0 and (3.14), we get

D(uk+1 − u∗ )

2

≤ ρ((T 2 )EE ) wk+1 − w∗

2

≤ ρ((T 2 )EE ) D(uk − u∗ ) 2 .

k The above two inequalities imply assertions 1 and 2 of this theorem. From (3.13) and wL = ∗ wL = 0, we have

uk+1 − u∗

2 M

= (wk+1 − w∗ )T T (wk+1 − w∗ ) ≤ ρ(TEE ) wk+1 − w∗ 2 .

Considering (3.14) and the definition of M , the above inequality becomes uk+1 − u∗
M

≤

ρ(TEE ) D(uk − u∗ ) ≤

ρ(TEE ) uk − u∗

M;

namely, the last conclusion holds, and thus uk converges to u∗ q-linearly. Theorem 3.6 states that Algorithm 1 converges q-linearly with the convergence rate depending on the spectral radii of the submatrices TEE and (T 2 )EE rather than on that of the whole matrix T . Since ρ(T ) ≤ 1 and TEE is a minor of T , it holds that ρ(TEE ) ≤ ρ(T ) ≤ 1. Similarly, ρ((T 2 )EE ) ≤ ρ(T 2 ) ≤ 1. In the next section, we will demonstrate that in practice ρ(TEE ) can be much smaller than ρ(T ). 4. A continuation scheme. In this section, we develop a continuation scheme on the penalty parameter β based on our convergence results and demonstrate its effectiveness. It is easy to see from (3.3) that a smaller β generally gives a smaller ρ(T ). As indicated in (3.9), a smaller β should also give a smaller set E and hence a smaller ρ(TEE ) (this fact will be verified by our numerical experiments). These observations strongly suggest the use of a continuation scheme on β in which β is initially small and gradually increases to the final

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

261

value. In this scheme, earlier subproblems corresponding to smaller β values can be solved quickly, and the later subproblems can also be solved relatively quickly by warm starting from the previous solutions. To demonstrate the need for continuation and its effectiveness, we carried out a set of experiments on the 128 × 128 image Checkerboard, a blocky image available in MATLAB. Our experiments were done under MATLAB using its Image Processing Toolbox. A motion blurring kernel of motion distance “len = 5” and angle “theta = 135” was applied to the image with additive Gaussian noise of zero mean and standard deviation 10−3 through MATLAB functions “fspecial,” “imfilter,” and “imnoise.” In our experiments, we set μ = 5 × 104 . In the first experiment, we compared the quantities ρ(T ) and ρ(TEE ) for different β values, both related to q-convergence rates. In addition, we also compared the theoretical upper bound, N0 w0 − w∗ /ω 2 , for finite convergence (see Theorem 3.5) with the observed k ˆ ˆ iteration number, N0 , for finite convergence (i.e., after N0 iterations wi , i ∈ L, remained zero until termination). To do these calculations, we ran Algorithm 1 for three β values: β = 24 , 29 , and 214 . For ∗ each β, an “exact” solution (u∗ , wβ ) of (1.3) was obtained using an extremely tight stopping β tolerance = 10−15 in (2.10). Afterwards, we reran the algorithm for each β and stopped ˆ once uk − u∗ / u∗ < 1.2 × 10−6 . The calculated values of ρ(T ), ρ(TEE ), N0 , and N0 β β for the three β values are given in Table 1.
Table 1 Comparison of theoretical convergence rates for different β. β ρ(T ) ρ(TEE ) N0 ˆ N0 24 0.9999 0.5188 5.6850e+006 5 29 0.9999 0.9955 1.1990e+017 139 214 0.9999 0.9999 1.2103e+018 3658

Data in Table 1 show that on the tested image (i) ρ(TEE ) can be much smaller than ρ(T ) when β is small and (ii) the theoretical worst-case upper bound for finite convergence can be vastly overconservative in practice, especially when β is small. In Figure 1, we plot the histories of the relative error, uk − u∗ / u∗ , and the observed β β convergence factor, uk+1 − u∗ M / uk − u∗ M , generated by Algorithm 1 for the three tested β β β values. The plots clearly illustrate that the convergence speed of Algorithm 1 is much faster for smaller β values. Specifically, when β increased from 24 to 214 , the required iteration number for the given accuracy jumped from less than 10 to over 3000. Similarly, the observed q-convergence factor increased from around 0.2 to close to 1. Our theoretical and numerical results strongly suggest the use of a continuation scheme on β. How to do this continuation most efficiently should be an interesting research issue in its own right, but will not be studied in depth in the current paper. Instead, we implement a very basic version, likely still far from optimal, to test the viability of our proposed framework. We call the resulting algorithm “fast total variation deconvolution” or FTVd, which is given below.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

262

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

Relative error
10
0

Observed convergence rate β=24 1 0.9 0.8 0.7 0.6 0.5 β=24 β=29 β=214 β=2
9

10

−1

β=214

10

−2

10

−3

10

−4

0.4 0.3 0.2

10

−5

10

−6

10

0

10

1

10

2

10

3

10

4

0.1 0 10

10

1

10

2

10

3

10

4

Figure 1. Convergence behavior for β = 24 , 29 , and 214 . The left plot is the relative error ek (β) = uk −u∗ / u∗ , and the right one is the observed q-linear convergence factor qk (β) = uk+1 −u∗ M / uk −u∗ M . β β β β In each plot, the horizontal axis represents the number of iterations.

Algorithm 2 (FTVd). Input f , K, μ > 0, β0 > 0, and βmax > β0 . Initialize u = f , up = 0, β = β0 , and > 0. While β ≤ βmax , Do (1) Run Algorithm 1 until (2.10) is met. (2) β ← 2 ∗ β. End Do One can modify the above framework to make it more flexible, though as it is the above basic implementation already works surprisingly well. For example, one could adaptively increase β and choose from one outer iteration to another (or use another stopping criterion for Algorithm 1). We tested the effectiveness of FTVd for β0 = 1 and βmax = 214 and in Figure 2 present the results in comparison to the results without doing continuation. From this comparison, it is abundantly clear that the algorithm would not be practically effective without continuation. 5. Numerical experiments. In this section, we present detailed numerical results comparing FTVd to the LD algorithm [41, 10], a state-of-the-art method for solving the isotropic TV deblurring model (1.2), as well as to some non-TV deblurring algorithms. 5.1. Overall assessment of several algorithms. We tested several popular algorithms for solving the TV deblurring model (1.2). Based on our experience, we have reached the following overall impression. The artificial time-marching algorithm used by Rudin and Osher in [37] and by Rudin, Osher, and Fatemi in [38] is easy to implement, but usually takes a large number of iterations, each with a small step size governed by the CFL condition, to reach a modest accuracy. On the other hand, the second-order cone programming (SOCP) approach [19] recovers solutions with a high accuracy in a small number of iterations (typically 30–50), but takes much more running time and memory in each iteration even on medium-sized images. A comprehensive comparison between the artificial time-marching algorithm and the SOCP approach is presented in [19]. We did not test the interior-point primal-dual implicit quadratic methods

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

263

Relative error
10
0

Observed convergence rate
1 0.9 0.8 0.7 0.6

with continuation without continuation 10
−1

10

−2

0.5 0.4

10

−3

0.3 0.2 0.1 with continuation without continuation 10
1

10

−4

0

500

1000

1500

2000

2500

3000

0 0 10

10

2

10

3

10

4

Figure 2. Continuation vs. no continuation: u∗ is an “exact” solution corresponding to β = 214 . The left plot is the relative error ek = uk − u∗ / u∗ , and the right one is the observed q-linear convergence factor qk = uk+1 − u∗ M / uk − u∗ M . In each plot, the horizontal axis represents the number of iterations.

presented in [9, 3, 8, 30], but expect their performance to be similar to SOCP because they also require solving large systems of linear equations at each iteration. We also tested the recently proposed iteratively reweighted least squares algorithm [43] and obtained speed and accuracy consistent with that reported by the authors. Specifically, it was slower than the LD method devised by Vogel and Oman [41, 10] for solving the TV/L2 model (1.2). Therefore, among all tested algorithms, LD is the most efficient for solving (1.2). Clearly, there are other deblurring algorithms not included in our discussion, but a more exhaustive comparison is beyond the scope of this work. 5.2. Test images and test platforms. We used two images, Lena and Man, in our experiments; see Figure 3. Image Lena is a good test image because it has a nice mixture of detail, flat regions, shading area, and texture. Image Man also consists of complex components in different scales, with different patterns and under inhomogeneous illuminations. Both images are suited for our experiments. We tested several kinds of blurring kernels including Gaussian, average, and motion, and found that the running times of FTVd and LD remained essentially constant for different blurring kernels as long as other conditions (e.g., kernel size or image size) were the same. Therefore, without loss of generality, we present only results using a “Gaussian” blurring kernel with different kernel sizes. In all tests, the additive noise used was Gaussian noise with zero mean and standard deviation 10−3 . This level of noise is substantial in the context of deblurring. Detailed information of our experiment set-up is summarized in Table 2. FTVd was implemented in MATLAB. All blurring effects were generated using the MATLAB function “imfilter” with periodic boundary conditions. For LD, we used an efficient C++ code NUMIPAD [36]. NUMIPAD uses Dirichlet boundary conditions for TV calculation (i.e., areas outside of the image boundary are treated as 0) instead of periodic boundary conditions. This difference caused the two codes to generate images with minor differences near boundaries, but should not noticeably affect their running times.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

264

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

Figure 3. Test images: Lena (left) and Man (right).

Table 2 Information on test images and blurring setting. Test no. 1. 2. Image Lena Man Size 512 × 512 1024 × 1024 Blurring type Gaussian Gaussian Blurring kernel parameters hsize={3, 5, . . . , 19, 21}, σ = 10 hsize={3, 5, . . . , 19, 21}, σ = 10

All comparisons between FTVd and LD were performed under GNU/Linux Release 2.6.955.0.2 and MATLAB v7.2 (R2006a) running on a Dell Optiplex GX620 with Dual Intel Pentium D CPU 3.20GHz (only one processor was used by MATLAB) and 3 GB of memory. The C++ code of LD was a part of the library package NUMIPAD [36] and was compiled with gcc v3.4.6. Since the code of ForWaRD had compiling problems on our Linux workstation, we chose to compare FTVd with ForWaRD, and also MATLAB functions “deconvwnr” and “deconvreg” under Windows XP and MATLAB v7.5 (R2007b) running on a Lenovo laptop with an Intel Core 2 Duo CPU at 2 GHz and 2 GB of memory. 5.3. Parameter values. In model (1.2), the parameter μ controls the amount of penalty applied to the squared L2 -distance between Ku and f . According to (1.1), an appropriate μ should give a solution u of (1.2) satisfying Ku − f ≈ ω . In our experiments, we used μ = 0.05/σ 2 , where σ is the standard deviation of the additive Gaussian noise ω in (1.1). This formula is based on the observation that μ should be inversely proportional to the noise variance, while the constant 0.05 was determined empirically so that the restored images had reasonable signal-to-noise ratios (to be defined later) and relative errors. More practical techniques exist for choosing μ, often by testing on a small window of image. The reader interested in this topic is referred to [39]. The issue of how to optimally select μ is important, but is outside of the scope of this work. As is usually done, we measure the quality of restoration by the signal-to-noise ratio

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

265

23 22 21 20

SNR (dB)

19 18 17 16 15 14 13 0 2 4 6 8 10 12 14 Lena Man

log β
2

16

18

Figure 4. SNRs of images recovered from (1.3) for different β.

(SNR), defined by SNR 10 ∗ log10 u0 − u ˜ 0−u u
2 2

,

where u0 is the original image and u is the mean intensity value of u0 . Generally, the quality ˜ of the restored image is expected to increase as β increases (hence, (1.3) becomes a closer approximation to (1.2)). In practice, we observed that the SNR value of recovered images from (1.3) stabilized once β reached a reasonably large value. To see this, we plot the SNR values of restored images corresponding to β = 20 , 21 , . . . , 218 in Figure 4. In this experiment, “Gaussian” out-of-focus blurring was applied to Lena with blurring size hsize = 11 and standard deviation σ = 5, and the “motion” blurring was applied to Man with motion distance “len = 21” and “theta = 135.” As can be seen from Figure 4, the SNR values on both images essentially remain constant for β ≥ 27 . This suggests that β need not be excessively large from a practical point of view. In our numerical experiments comparing FTVd with LD, we set β0 = 1 and βmax = 27 in Algorithm 1. For each β, the inner iteration was stopped once (2.10) was satisfied with = 0.05. 5.4. Comparison of FTVd and LD. In this subsection, we compare FTVd with LD in terms of both convergence speed and quality of recovery. As mentioned in subsection 5.2, the two methods use different boundary conditions, which causes minor differences. The blurry images with the circular boundary condition usually have slightly higher SNRs than those blurred with the Neumann boundary condition, as do the restored images. To have a fairer comparison, we measure the restoration quality by the so-called improvement SNR (ISNR) defined as f − u0 2 , ISNR 10 ∗ log10 u − u0 2

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

266

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

9

10

4

8.5 10 8
3

FTVd: Test No.1 LD: Test No.1 FTVd: Test No.2 LD: Test No.2

Running times (s)
FTVd: Test No.1 LD: Test No.1 FTVd: Test No.2 LD: Test No.2 4 6 8 10 12 14 16 18 20

ISNR (dB)

7.5

10

2

7

10

1

6.5

6

10

0

4

6

8

10

12

14

16

18

20

hsize

hsize

Figure 5. ISNR (left) and CPU time (right) comparisons for FTVd and LD methods. In both figures, the horizontal axis represents the size of the blurring kernel.

where u0 , u, and f are the original, restored, and observed images, respectively. This quantity measures the quality of a restored image u relative to the blurry and noisy observation f . The ISNRs of the recovered images computed by the FTVd and LD methods for the two test images are given in the left-hand plot of Figure 5. We observe that FTVd produces slightly lower ISNRs when hsize is small and slightly higher ISNRs when hsize is large. However, since both FTVd and LD solve the same model (boundary conditions aside), the restored images have few visible differences. Therefore, we will not display the restored images here. The right-hand plot in Figure 5 shows the comparison of running time between FTVd and LD. While generating the restored images with similar qualities, our MATLAB code FTVd is much faster than the C++ code for LD. The running time of FTVd essentially remains constant for different values of hsize and increases only moderately when image size is doubled. However, the running time for LD increases dramatically with the increase in both the image size and the blurring kernel size. As hsize becomes relatively large, the CPU times consumed by FTVd and LD become orders of magnitude apart. The reason is clear: Larger hsize makes the matrix K denser and more ill-conditioned; hence, solving a linear system involving K becomes more difficult for the preconditioned conjugate gradient (CG) method used in NUMIPAD [36]. On the other hand, the increase of kernel size does not noticeably increase the cost of doing FFT on the function represented by K. 5.5. Observations on FTVd’s behavior and a note. As has been mentioned, the periteration computation of FTVd is dominated by three FFTs, each at the cost of O(n2 log(n)). The question is how many iterations are needed in general for FTVd to attain a required accuracy for images of different sizes. In our experiments, the number of total inner iterations taken by the basic version of FTVd, with the given “default” parameters, is almost always around 12. These total inner iteration numbers are quite reasonable, given that there were 8 outer iterations corresponding to the β-sequence {20 , 21 , . . . , 27 }. For each β > 1, since FTVd has already obtained a good starting point, on average it takes only one or two inner

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

267

iterations to reach the prescribed accuracy. Given that three FFTs are required at each inner iteration, the total number of FFTs taken by FTVd is around 40 for all tests. Upon examining detailed timing information, we observed that over two-thirds of the CPU time was spent on the FFT and related calculations in (2.4). Finally, FTVd is numerically stable and insensitive to ill-conditioning, apparently because it does not require any matrix operations. After we finished the writing of this paper, we found a newly released paper by Huang, Ng and Wen [23] which solves the following approximation to the TV/L2 model (after necessary notational changes for consistency): (5.1) min u,v i

Di v + α v − u

2

+

μ Ku − f 2

2

,

by alternating minimization. For a fixed v, the minimization with respect to u involves two FFTs, one less than what is required by the corresponding subproblem in our method. However, for a fixed u, the minimization with respect to v is nothing but a TV denoising problem which does not have a closed-form solution with linear complexity as our corresponding subproblem has (see formula (2.2)). In [23], the TV denoising problem is solved iteratively by Chambolle’s projection algorithm [6]. While the per-iteration computational complexity of our method is dominated by three FFTs, that of [23] is dominated by the cost of solving a TV denoising problem in addition to two FFTs. According to the reported numerical results in [23], their algorithm appears to require at least as many outer iterations as ours. 5.6. Comparison of FTVd with other methods. We also conducted tests to show how FTVd compares with some other state-of-the-art, non-TV methods. A hybrid Fourier-wavelet regularized deconvolution (ForWaRD) algorithm by Neelanmani, Choi, and Baraniuk [32] uses shrinkage in both Fourier and wavelet domains to recover images while preserving edges and removing noise. Their code was written in the C programming language with a MATLAB interface.1 We compared the quality of images recovered by FTVd, ForWaRD, and two MATLAB deblurring functions “deconvwnr” (Wiener filter) and “deconvreg” (Laplacian lowpass filter). These methods are all based on distinct models. In this experiment, we used the image Lena, which was blurred with a “Gaussian” kernel of hsize = 21 and standard deviation σ = 11. We ran FTVd twice, once with β = 25 and tolerance = 0.05 in (2.10), and again with β = 27 and = 0.002. The blurry and noisy and the restored images are presented in Figure 6. SNRs of the restored images and running times of different algorithms are also given. As can be seen from Figure 6, FTVd is comparable to ForWaRD in speed and attains a better image quality. The MATLAB functions are faster as they solve simpler models, but they generate images of lower quality. Specifically, the artifacts of multiple rings or ripples are clearly more visible in the restored images by the three methods than in the ones by FTVd. By visually comparing the restored images, we give preference to the results of FTVd. 5.7. Summary of numerical results. First, FTVd can be orders of magnitude faster than the LD method. The performance gap between them will continue to grow as the image and blurring kernel sizes increase. Since LD compared favorably to many other existing algorithms
1

The code ForWaRD can be downloaded at http://www.dsp.rice.edu/software/ward.shtml.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

268

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG
Blurry&Noisy. SNR: 5.19dB ForWaRD: SNR: 12.46dB, t = 4.88s

FTVd: β = 25, SNR: 12.58dB, t = 1.83s

FTVd: β = 2 , SNR: 13.11dB, t = 14.10s

7

deconvwnr: SNR: 11.51dB, t = 0.05s

deconvreg: SNR: 11.20dB, t = 0.34s

Figure 6. The first row contains the blurry and noisy image and the image recovered by ForWaRD; the second row contains the results of FTVd (left: β = 25 , = 5 × 10−2 ; right: β = 27 , = 2 × 10−3 ); and the third row is recovered by deconvwnr (left) and deconvreg (right).

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

269

for solving (1.2) with isotropic TV, the lead of FTVd obviously extends to those algorithms. Second, compared to some other non-TV deblurring algorithms such as ForWaRD and the MATLAB deblurring functions, FTVd generally produces images of higher quality in our experiments on numerous natural images. Under the “default” parameter setting, the speed of FTVd is not as fast as the MATLAB functions “deconvwnr” and “deconvreg.” However, with a less restrictive stopping tolerance and a smaller βmax , the running time of FTVd can be faster than that of “deconvreg” on large images and close to that of “deconvwnr” while still producing images of higher quality. 6. Extensions. There are some obvious extensions in which one can apply the variablesplitting and alternating-minimization methodology to solving other problems. Here we mention a few. Color or other multichannel image reconstruction with TV regularization. This extension is practically important, as straightforward as it may be, because it enables multichannel TV image reconstruction problems to be solved at a speed not seen before. We are currently working on this extension and will report our findings in a forthcoming paper. In addition, we can also extend this methodology to models where the fidelity term Ku − f is not measured in the 2-norm but in the 1-norm in order to handle impulsive noise (such as “salt and pepper”) other than white noise. The use of higher-order finite difference operators. In this case, we may have models in the form of (6.1) min u i

Di u +

(p)

μ Ku − f 2

2

,

for some p > 2, which can be approximated by (6.2) min w,u i

wi +

β (p) wi − Di u 2
(p)

2

+

μ Ku − f 2

2

.

As long as the difference operators Di are linear and shift-invariant, minimization with respect to u is as simple as a few FFTs for a fixed w, while minimization with respect to w can be done by shrinkage in an appropriate dimension (depending on the order of the finite difference involved). Higher-order derivatives of u are used in [27, 44], for example, to reduce staircasing effects. Operator K is a partial Fourier transform matrix. In this case, Ku = P F(u), where 2 P ∈ Rm×n is a selection matrix consisting of m < n2 rows of the identity matrix. This case occurs, for example, in compressive sensing [4, 16] where a subset of noisy Fourier coefficients, f , is used to recover a signal u0 with a sparse gradient field via solving the model (6.3) min u i

Di u +

μ P F(u) − f 2

2

.

It is easy to verify that Algorithm 1 can be applied to the above problem with minor modifi2 cations to formula (2.4), including replacing F(K) by a section vector p ∈ Rn whose entries are 1 corresponding to rows selected by P and 0 otherwise. We have applied this approach to

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

270

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

solving simulation problems in compressive magnetic resonance (MR) imaging with promising results. Multiple regularization terms. Such examples arise in both image processing (e.g., [28, 8]) and compressive MR imaging (e.g., [22, 26, 29]). For example, the model (6.4) min u i

Di u + α Φu

1

+

μ P F(u) − f 2

2

,

where Φ is an orthonormal matrix, can be approximated by (6.5) min wi + i u,v,w

β wi − Di u 2

2

+α

v

1

+

γ v − Φu 2

2

+

μ P F(u) − f 2

2

.

In an alternating minimization setting, w can be solved by (2.2) as before, v can be solved by 1D shrinkage, and u can be solved by four FFTs (two for w1 and w2 , one for ΦT v noting that v − Φu = ΦT v − u , and one inverse FFT). In general, one needs to introduce an auxiliary variable for each regularization term and manages to obtain easily solvable subproblems. Finally, we mention that even in the case where the quadratic part in (1.3) cannot be minimized with respect to u via fast transforms as in formula (2.4), model (1.3) may still lead to an efficient algorithm if the quadratic in u can be effectively decreased by a few iterations of an iterative method such as a preconditioned CG method. This topic is of interest for further investigation. 7. Concluding remarks. We proposed, analyzed, and tested a new algorithm FTVd for solving the TV/L2 model (1.2). Our approximation model (1.3) was derived from the classic quadratic penalty method first proposed by Courant in 1943, but also falls into the more recent half-quadratic class for image processing after our extensions. In fact, it can be regarded as a long missing, isotropic TV model in the half-quadratic class whose quadratic portion has a block circulant Hessian, thus allowing fast transforms to be utilized in quadratic minimization. We established strong convergence results for our algorithm and validated a continuation scheme. Our numerical results convincingly show that FTVd can be orders of magnitude faster than the LD algorithm. In particular, our results demonstrate that the TV/L2 model can be solved just as quickly as some other more widely used models, putting it on an equal footing with others in terms of solution speed. The results of this paper can be extended to a number of models involving TV regularization. More studies are underway to further explore the potential of the proposed approach. Acknowledgments. We are grateful to Paul Rodr´ iguez and Brendt Wohlberg at the Los Alamos National Lab for valuable discussions and for making their algorithm library NUMIPAD [36] available to us.
REFERENCES [1] R. Acar and C. R. Vogel, Analysis of total variation penalty methods, Inverse Problems, 10 (1994), pp. 1217–1229. [2] M. Allain, J. Idier, and Y. Goussard, On global and local convergence of half-quadratic algorithms, IEEE Trans. Image Process., 15 (2006), pp. 1130–1142.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

AN ALTERNATING ALGORITHM FOR TOTAL VARIATION MINIMIZATION

271

[3] P. Blomgren, T. F. Chan, and C. K. Wong, Total variation image restoration: Numerical methods and extensions, in Proceedings of the IEEE International Conference on Image Processing, Vol. 3, 1997, pp. 384–387. ` [4] E. Candes, J. Romberg, and T. Tao, Stable signal recovery from incomplete and inaccurate measurements, Comm. Pure Appl. Math., 59 (2006), pp. 1207–1223. [5] A. S. Carasso, Linear and nonlinear image deblurring: A documented study, SIAM J. Numer. Anal., 36 (1999), pp. 1659–1689. [6] A. Chambolle, An algorithm for total variation minimization and applications, J. Math. Imaging Vision, 20 (2004), pp. 89–97. [7] A. Chambolle and P. L. Lions, Image recovery via total variation minimization and related problems, Numer. Math., 76 (1997), pp. 167–188. [8] T. F. Chan, S. Esedoglu, F. Park, and A. Yip, Total variation image restoration: Overview and recent developments, in Handbook of Mathematical Models in Computer Vision, edited by N. Paragios, Y. Chen, and O. Faugeras, Springer-Verlag, New York, 2006, pp. 17–31. [9] T. F. Chan, G. H. Golub, and P. Mulet, A nonlinear primal-dual method for total variation-based image restoration, SIAM J. Sci. Comput., 20 (1999), pp. 1964–1977. [10] T. F. Chan and P. Mulet, On the convergence of the lagged diffusivity fixed point method in total variation image restoration, SIAM J. Numer. Anal., 36 (1999), pp. 354–367. [11] T. F. Chan and J. Shen, Theory and Computation of Variational Image Deblurring, IMS Lecture Notes, Institute of Mathematical Statistics, Beachwood, OH, 2006. ´ [12] P. Charbonnier, L. Blanc-Feraud, G. Aubert, and M. Barlaud, Deterministic edge preserving regularization in computed imaging, IEEE Trans. Image Process., 6 (1997), pp. 298–311. [13] R. Courant, Variational methods for the solution of problems with equilibrium and vibration, Bull. Amer. Math. Soc., 49 (1943), pp. 1–23. [14] G. Demoment, Image reconstruction and restoration: Overview of common estimation structures and problems, IEEE Trans. Acoust. Speech Signal Process., 37 (1989), pp. 2024–2036. [15] D. C. Dobson and F. Santosa, Recovery of blocky images from noisy and blurred data, SIAM J. Appl. Math., 56 (1996), pp. 1181–1198. [16] D. Donoho, Compressed sensing, IEEE Trans. Inform. Theory, 52 (2006), pp. 1289–1306. [17] D. Geman and G. Reynolds, Constrained restoration and the recovery of discontinuities, IEEE Trans. Pattern Anal. Mach. Intell., 14 (1992), pp. 367–383. [18] D. Geman and C. Yang, Nonlinear image recovery with half-quadratic regularization, IEEE Trans. Image Process., 4 (1995), pp. 932–946. [19] D. Goldforb and W. Yin, Second-order cone programming methods for total variation-based image restoration, SIAM J. Sci. Comput., 27 (2005), pp. 622–645. [20] T. Goldstein and S. Osher, The Split Bregman Algorithm for L1 Regularized Problems, UCLA CAM Report 08-29, UCLA, Los Angeles, 2008. [21] R. Gonzalez and R. Woods, Digital Image Processing, Addison-Wesley, Reading, MA, 1992. [22] L. He, T.-C. Chang, S. Osher, T. Fang, and P. Speier, MR Image Reconstruction by Using the Iterative Refinement Method and Nonlinear Inverse Scale Space Methods, UCLA CAM Report 06-35, UCLA, Los Angeles, 2006. [23] Y. Huang, M. K. Ng, and Y.-W. Wen, A fast total variation minimization method for image restoration, Multiscale Model. Simul., 7 (2008), pp. 774–795. [24] J. Idier, Convex half-quadratic criteria and interacting auxiliary variables for image restoration, IEEE Trans. Image Process., 10 (2001), pp. 1001–1009. [25] A. K. Katsaggelos, ed., Digital Image Restoration, Springer-Verlag, Berlin, 1991. [26] M. Lustig, D. Donoho, and J. Pauly, Sparse MRI: The application of compressed sensing for rapid MR imaging, Magn. Reson. Med., 58 (2007), pp. 1182–1195. [27] M. Lysaker, A. Lundervold, and X. C. Tai, Noise removal using fourth-order partial differential equations with applications to medical magnetic resonance images in space and time, IEEE Trans. Image Process., 12 (2002), pp. 1579–1590. [28] O. M. Lysaker and X.-C. Tai, Iterative image restoration combining total variation minimization and a second-order functional, Internat. J. Comput. Vision, 66 (2006), pp. 5–18.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

272

YILUN WANG, JUNFENG YANG, WOTAO YIN, AND YIN ZHANG

[29] S. Ma, W. Yin, Y. Zhang, and A. Chakraborty, Compressed MR imaging based on wavelets and total variation, in Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’08), 2008. [30] A. Marquina and S. Osher, Explicit algorithms for a new time dependent model based on level set motion for nonlinear deblurring and noise removal, SIAM J. Sci. Comput., 22 (2000), pp. 387–405. [31] R. Neelamani, H. Choi, and R. G. Baraniuk, Wavelet-based deconvolution for ill-conditioned systems, in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing Vol. 6, 1999, pp. 3241–3244. [32] R. Neelamani, H. Choi, and R. G. Baraniuk, ForWaRD: Fourier-wavelet regularized deconvolution for ill-conditioned systems, IEEE Trans. Signal Process., 52 (2004), pp. 418–433. [33] M. K. Ng, R. H. Chan, and W.-C. Tang, A fast algorithm for deblurring models with Neumann boundary conditions, SIAM J. Sci. Comput., 21 (1999), pp. 851–866. [34] M. Nikolova and M. K. Ng, Analysis of half-quadratic minimization methods for signal and image recovery, SIAM J. Sci. Comput., 27 (2005), pp. 937–966. [35] J. Nocedal and S. J. Wright, Numerical Optimization, Springer-Verlag, New York, 1999. [36] P. Rodr´ iguez and B. Wohlberg, Numerical Methods for Inverse Problems and Adaptive Decomposition (NUMIPAD), http://numipad.sourceforge.net/. [37] L. Rudin and S. Osher, Total variation based image restoration with free local constraints, in Proceedings of the 1st IEEE International Conference on Image Processing, Vol. 1, 1994, pp. 31–35. [38] L. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Phys. D, (1992), pp. 259–268. [39] D. Strong and T. F. Chan, Edge-preserving and scale-dependent properties of total variation regularization, Inverse Problems, 19 (2003), pp. 165–187. [40] A. Tikhonov and V. Arsenin, Solution of Ill-Posed Problems, Winston, Washington, DC, 1977. [41] C. R. Vogel and M. E. Oman, Iterative methods for total variation denoising, SIAM J. Sci. Comput., 17 (1996), pp. 227–238. [42] Y. Wang, W. Yin, and Y. Zhang, A Fast Algorithm for Image Deblurring with Total Variation Regularization, CAAM Technical Report 07-10, Rice University, Houston, TX, 2007. [43] B. Wohlberg and P. Rodr´ iguez, An iteratively reweighted norm algorithm for minimization of total variation functionals, IEEE Signal Process. Lett., 14 (2007), pp. 948–951. [44] Y. L. You and M. Kaveh, Fourth-order partial differential equations for noise removal, IEEE Trans. Image Process., 9 (2000), pp. 1723–1730.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.