## A new alternating minimization algorithm for total variation image reconstruction (2008)

Venue: | SIAM J. IMAGING SCI |

Citations: | 104 - 17 self |

### BibTeX

@ARTICLE{Wang08anew,

author = {Yilun Wang and Junfeng Yang and Wotao Yin and Yin Zhang},

title = { A new alternating minimization algorithm for total variation image reconstruction},

journal = {SIAM J. IMAGING SCI},

year = {2008},

pages = {248--272}

}

### OpenURL

### 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 isotropic forms of total variation discretizations. The per-iteration computational complexity of the algorithm is three Fast Fourier Transforms (FFTs). 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 total-variation-based deblurring. Some extensions of our algorithm are also discussed.

### Citations

1851 | Compressed sensing
- Donoho
(Show Context)
Citation Context ...artial Fourier transform matrix. In this case, Ku = P F(u) where P ∈ R m×n2 a selection matrix consisting of m < n 2 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 u 0 with a sparse gradient field via solving the model (6.3) min u ∑ i ‖Diu‖ + µ 2 ‖P F(u) − f‖2 . It is easy to verify th... |

1448 |
Nonlinear total variation based noise removal algorithms
- Rudin, Osher, et al.
- 1992
(Show Context)
Citation Context ... important image attributes such as sharp edges. Another well-known class of regularizers are based on total variation (TV), which was first proposed for image denoising by Rudin, Osher and Fatemi in =-=[39]-=-, and then extended to image deblurring in [38]. In comparison to Tikhonov-like regularizers, TV regularizers can better preserve sharp edges or object boundaries that are usually the most important f... |

868 |
Solutions of ill-posed problems
- Tikhonov, Arsenin
- 1977
(Show Context)
Citation Context ... where Φreg(u) is the regularizer that enforces the 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 =-=[41]-=- 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 ... |

793 |
Stable signal recovery from incomplete and inaccurate measurements
- Candès, Romberg, et al.
- 2006
(Show Context)
Citation Context ...artial Fourier transform matrix. In this case, Ku = P F(u) where P ∈ R m×n2 a selection matrix consisting of m < n 2 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 u 0 with a sparse gradient field via solving the model (6.3) min u ∑ i ‖Diu‖ + µ 2 ‖P F(u) − f‖2 . It is easy to verify th... |

365 | An algorithm for total variation minimization and applications
- Chambolle
- 2004
(Show Context)
Citation Context ...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 3 FFTs, that of [23] is dominated by the cost of solving a TV denoising problem in addition to 2 FFTs. According to the... |

275 |
Constrained restoration and the recovery of discontinuity
- Geman, Reynolds
- 1992
(Show Context)
Citation Context ...ified 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 halfquadratic 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 const... |

258 |
Image Recovery via Total Variation Minimization and Related Problems
- Chambolle, Lions
(Show Context)
Citation Context ...ced as the paper progresses. 21.2. Existing Regularization Approaches. There are different approaches based on statistics [14, 5], Fourier and/or wavelet transforms [25, 31], or variational analysis =-=[38, 7, 11]-=- for image deblurring. Among them the simplest is the Maximum Likelihood estimation, which solves the least square problem minu ‖Ku− f‖ 2 and is also known as the inverse filter. However, the solution... |

242 |
Sparse MRI: The application of compressed sensing for rapid
- Lustig, Donoho, et al.
(Show Context)
Citation Context ...ms 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 ∑ i ‖Diu‖ + α‖Φu‖1 + µ 2 ‖P F(u) − f‖2 , where Φ is an orthonormal matrix, can be approximated by ( ∑ (6.5) min ‖wi‖ + u,v,w β 2 ‖wi − Diu‖ 2 ) ( + α ‖v‖1 + γ 2... |

241 | Iteration Methods for Total Variation Denoising
- Vogel, Oman
(Show Context)
Citation Context ...hod has recently been proposed that uses Bregman iterations [20]. We have implemented the proposed algorithm in MATLAB TM and compared it with a C++ implementation of the Lagged Diffusivity algorithm =-=[42]-=-, 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 bl... |

151 |
Image reconstruction and restoration: overview of common estimation structures and problems
- Demoment
- 1989
(Show Context)
Citation Context ...to the 2-norm unless otherwise indicated. Additional notation will be introduced as the paper progresses. 21.2. Existing Regularization Approaches. There are different approaches based on statistics =-=[14, 5]-=-, Fourier and/or wavelet transforms [25, 31], or variational analysis [38, 7, 11] for image deblurring. Among them the simplest is the Maximum Likelihood estimation, which solves the least square prob... |

139 | Nonlinear image recovery with half-quadratic regularization and FFT’s - Geman, Yang - 1999 |

135 |
Total variation based image restoration with free local constraints
- Osher, Rudin
- 1994
(Show Context)
Citation Context ...ced as the paper progresses. 21.2. Existing Regularization Approaches. There are different approaches based on statistics [14, 5], Fourier and/or wavelet transforms [25, 31], or variational analysis =-=[38, 7, 11]-=- for image deblurring. Among them the simplest is the Maximum Likelihood estimation, which solves the least square problem minu ‖Ku− f‖ 2 and is also known as the inverse filter. However, the solution... |

74 | Recovery of blocky images from noisy and blurred data
- Dobson, Santosa
- 1996
(Show Context)
Citation Context ...mportant features to recover. The TV/L 2 model (1.2) has been widely studied (see, for example, [7, 8] and references thereof). The superiority of TV over Tikhonov-like regularization was analyzed in =-=[1, 15]-=- for recovering images containing piecewise smooth objects. However, due to the non-differentiability and non-linearity of the TV functions, the TV/L 2 model (1.2) is computationally more difficult to... |

71 | A fast algorithm for deblurring models with Neumann boundary conditions
- Ng, Chan, et al.
- 2000
(Show Context)
Citation Context ...) + (D (2) ) T D (2) + µ β KT K ) u = (D (1) ) T w1 + (D (2) ) T w2 + µ β KT 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 two-dimensional discrete Fourier transform F. Using the convolution theorem of Fourier transforms, w... |

71 | Edge-preserving and scale dependent properties of total variation regularization
- Chan
(Show Context)
Citation Context ...se ratios (will be defined latter) 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 =-=[40]-=-. The issue of how to optimally select µ is important, but outside of the scope of this work. 16As is usually done, we measure the quality of restoration by the signal-to-noise ratio (SNR), defined b... |

66 | Explicit algorithms for a new time dependent model based on level set motion for nonlinear deblurring and noise removal
- Marquina, Osher
(Show Context)
Citation Context ...prehensive 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 presented in =-=[9, 3, 30, 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 re-weighted... |

51 |
Analysis of total variation penalty methods
- Acar, Vogel
- 1994
(Show Context)
Citation Context ...mportant features to recover. The TV/L 2 model (1.2) has been widely studied (see, for example, [7, 8] and references thereof). The superiority of TV over Tikhonov-like regularization was analyzed in =-=[1, 15]-=- for recovering images containing piecewise smooth objects. However, due to the non-differentiability and non-linearity of the TV functions, the TV/L 2 model (1.2) is computationally more difficult to... |

50 | Recent developments in total variation image restoration
- Chan, Esedoglu, et al.
- 2005
(Show Context)
Citation Context ...rizers, TV regularizers can better preserve sharp edges or object boundaries that are usually the most important features to recover. The TV/L 2 model (1.2) has been widely studied (see, for example, =-=[7, 8]-=- and references thereof). The superiority of TV over Tikhonov-like regularization was analyzed in [1, 15] for recovering images containing piecewise smooth objects. However, due to the non-differentia... |

49 | X.C.: Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time
- Lysaker, Lundervold, et al.
- 2003
(Show Context)
Citation Context ...ixed 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 P ∈ R m×n2 a selection matrix consisting of m < n 2 rows of the identity... |

48 | Second-order cone programming methods for total variation-based image restoration
- Goldfarb, Yin
- 2005
(Show Context)
Citation Context ...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 14Relative error Observed convergence rate 10 0 10 −1 10 −2 with continuation without contin... |

45 | On the convergence of the lagged diffusivity fixed point method in total variation image restoration
- Chan, Mulet
- 1999
(Show Context)
Citation Context ...rithm would not be practically effective without continuation. 5. Numerical Experiments. In this section, we present detailed numerical results comparing FTVd to the Lagged Diffusivity (LD) algorithm =-=[42, 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 popul... |

39 |
Fourth-order partial differential equations for noise removal
- You, Kaveh
(Show Context)
Citation Context ...ixed 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 P ∈ R m×n2 a selection matrix consisting of m < n 2 rows of the identity... |

38 | Convex half-quadratic criteria and interacting auxiliary variables for image restoration
- Idier
(Show Context)
Citation Context ...tudied: multiplicative form [17] in which Q(t, s) = 1 2 t2 s and additive form [18] in which Q(t, s) = 1 2 (t − s)2 . For comparisons of these two forms of methods both in theory and in practice, see =-=[24, 2, 34, 12]-=- and 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... |

36 |
The split Bregman algorithm for L1 regularized problems
- Goldstein, Osher
(Show Context)
Citation Context ...l (1.3) and the resulting alternating minimization algorithm were first proposed in [43] 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 TM and compared it with a C++ implementation of the Lagged Diffusivity algorithm [42], which is one of the fastest algorithms for solving (1.2). ... |

26 |
A fast total variation minimization method for image restoration
- Huang, Ng, et al.
(Show Context)
Citation Context ...ered 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) 19=-=[23]-=- which solves the following approximation to the TV/L 2 model (after necessary notational changes for consistency) (5.1) min u,v ∑ i ‖Div‖ + α‖v − u‖ 2 + µ 2 ‖Ku − f‖2 , by alternating minimization. F... |

24 |
MR image reconstruction by using the iterative refinement method and nonlinear inverse scale space methods, UCLA
- He, Chang, et al.
- 2006
(Show Context)
Citation Context ...ms 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 ∑ i ‖Diu‖ + α‖Φu‖1 + µ 2 ‖P F(u) − f‖2 , where Φ is an orthonormal matrix, can be approximated by ( ∑ (6.5) min ‖wi‖ + u,v,w β 2 ‖wi − Diu‖ 2 ) ( + α ‖v‖1 + γ 2... |

23 | Wavelet-based deconvolution for ill-conditioned systems
- Neelamani, Choi, et al.
- 1999
(Show Context)
Citation Context ...dditional notation will be introduced as the paper progresses. 21.2. Existing Regularization Approaches. There are different approaches based on statistics [14, 5], Fourier and/or wavelet transforms =-=[25, 31]-=-, or variational analysis [38, 7, 11] for image deblurring. Among them the simplest is the Maximum Likelihood estimation, which solves the least square problem minu ‖Ku− f‖ 2 and is also known as the ... |

22 |
Analysis of half-quadratic minimization methods for signal and image recovery
- Nikolova, Ng
- 2005
(Show Context)
Citation Context ...tudied: multiplicative form [17] in which Q(t, s) = 1 2 t2 s and additive form [18] in which Q(t, s) = 1 2 (t − s)2 . For comparisons of these two forms of methods both in theory and in practice, see =-=[24, 2, 34, 12]-=- and 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... |

21 |
and nonlinear image deblurring: A documented study
- Carasso, “Linear
- 1999
(Show Context)
Citation Context ...to the 2-norm unless otherwise indicated. Additional notation will be introduced as the paper progresses. 21.2. Existing Regularization Approaches. There are different approaches based on statistics =-=[14, 5]-=-, Fourier and/or wavelet transforms [25, 31], or variational analysis [38, 7, 11] for image deblurring. Among them the simplest is the Maximum Likelihood estimation, which solves the least square prob... |

21 | A fast algorithm for image deblurring with total variation regularization,” Rice University CAAM
- Wang, Yin, et al.
- 2007
(Show Context)
Citation Context ...irst 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 =-=[43]-=- 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 TM and compared it with a C++... |

20 | Iterative image restoration combining total variation minimization and a second order functional
- Lysaker, Tai
- 2006
(Show Context)
Citation Context ... this approach to 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 ∑ i ‖Diu‖ + α‖Φu‖1 + µ 2 ‖P F(u) − f‖2 , where Φ is an orthonormal matrix, can be approximated by ( ∑ (6.5) min ... |

19 | Total variation image restoration: Numerical methods and extensions
- Blomgren, Chan, et al.
- 1997
(Show Context)
Citation Context ...prehensive 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 presented in =-=[9, 3, 30, 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 re-weighted... |

16 |
On global and local convergence of half-quadratic algorithms
- Allain, Idier, et al.
(Show Context)
Citation Context ...tudied: multiplicative form [17] in which Q(t, s) = 1 2 t2 s and additive form [18] in which Q(t, s) = 1 2 (t − s)2 . For comparisons of these two forms of methods both in theory and in practice, see =-=[24, 2, 34, 12]-=- and 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... |

10 |
Variational methods for the solution of problems with equilibrium and vibration
- Courant
- 1943
(Show Context)
Citation Context ...): (1.3) min w,u ∑ i ‖wi‖2 + β 2 ∑ i ‖wi − Diu‖ 2 2 + µ 2 ‖Ku − f‖2 2, with a sufficiently large penalty parameter β. This type of quadratic penalty approach can be traced back to 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 whi... |

9 | Theory and computation of variational image deblurring
- Chan, Shen
(Show Context)
Citation Context ...ced as the paper progresses. 21.2. Existing Regularization Approaches. There are different approaches based on statistics [14, 5], Fourier and/or wavelet transforms [25, 31], or variational analysis =-=[38, 7, 11]-=- for image deblurring. Among them the simplest is the Maximum Likelihood estimation, which solves the least square problem minu ‖Ku− f‖ 2 and is also known as the inverse filter. However, the solution... |

7 |
Deterministic edge preserving regularization in computed imaging
- Charbonnier, Blanc-Féraud, et al.
- 1997
(Show Context)
Citation Context |

3 |
Fourier-wavelet regularized deconvolution for illconditioned systems
- ForWaRd
- 2005
(Show Context)
Citation Context ...orithm 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]-=-), 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 i... |

2 |
Compressed MR imaging based on wavelets and total variation
- Ma, Yin, et al.
- 2008
(Show Context)
Citation Context ...ms 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 ∑ i ‖Diu‖ + α‖Φu‖1 + µ 2 ‖P F(u) − f‖2 , where Φ is an orthonormal matrix, can be approximated by ( ∑ (6.5) min ‖wi‖ + u,v,w β 2 ‖wi − Diu‖ 2 ) ( + α ‖v‖1 + γ 2... |

1 |
Numerical methods for inverse problems and adaptive decomposition (NUMIPAD). http://numipad.sourseforge.net
- Nocedal, Wright, et al.
(Show Context)
Citation Context ...les 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 is enough for attaining a given ac... |