9.1 A class of fast algorithms for total variation image restoration  (Page 3/6)

for $i=1,...,n^2$ (with certain boundary conditions assumed for $i>n^2-n$ ). Then the discrete form of TV defined in ( ) is given by

We will refer to

with discretized TV regularization ( ) as TV/L ${}^2$ . For impulsive noise, we replace the ${\ell }_2$ fidelity by ${\ell }_1$ fidelity and refer to the resulted problem as TV/L ${}^1$ .

Now we introduce several more notation. For simplicity, we let ${\sum }_i$ be the summation taken over all pixels. The two first-order global finite difference operators in horizontal and verticaldirections are, respectively, denoted by $D^{\left(1\right)}$ and $D^{\left(2\right)}$ which are $n^2$ -by- $n^2$ matrices (boundary conditions are the same as those assumed on $D_i$ ). As such, it is worth noting that the two-row matrix $D_i$ is formed by stacking the $i$ th row of $D^{\left(1\right)}$ on that of $D^{\left(2\right)}$ . For vectors $v_1$ and $v_2$ , we let $v=(v_1;v_2)\triangleq {(v_1^{\top },v_2^{\top })}^{\top }$ , i.e. $v$ is the vector formed by stacking $v_1$ on the top of $v_2$ . Similarly, we let $D=(D^{\left(1\right)};D^{\left(2\right)})={\left({\left(D^{\left(1\right)}\right)}^{\top },,,{\left(D^{\left(2\right)}\right)}^{\top }\right)}^{\top }$ . Given a matrix $T$ , we let $\mathrm{diag}\left(T\right)$ be the vector containing the elements on the diagonal of $T$ , and $\mathcal{F}\left(T\right)=\mathbf{F}T{\mathbf{F}}^{-1}$ , where $\mathbf{F}\in n^2\times n^2$ is the 2D discrete Fourier transform matrix.

Existing methods

Since TV is nonsmooth, quite a few algorithms are based on smoothing the TV term and solving an approximation problem. The TV of $u$ is usually replaced by

where $ϵ>0$ is a small constant. Then the resulted approximate TV/L ${}^2$ problem is smooth and many optimization methods are available. Among others, the simplest method is the gradientdescent method as was used in . However, this method suffers slow convergence especially when the iterate point is closeto the solution. Another important method is the linearized gradient method proposed in for denoising and in for deblurring. Both the gradient descent and the linearized gradient methods are globally and at best linearlyconvergent. To obtain super linear convergence, a primal-dual based Newton method was proposed in . Both the linearized gradient method and this primal-dual method need to solvea large system of linear equations at each iteration. When $ϵ$ is small and/or $K$ becomes more ill-conditioned, the linear system becomes more and more difficult to solve. Anotherclass of well-known methods for TV/L ${}^2$ are the iterative shrinkage/thresholding (IST) based methods . For IST-based methods, a TV denoising problem needs to be solved ateach iteration. Also, in the authors transformed the TV/L ${}^2$ problem into a second order cone program and solved it by interior point method.

A new alternating minimization algorithm

In this section, we derive a new algorithm for the TV/L ${}^2$ problem

In ( ), the fidelity term is quadratic with respect to $u$ . Moreover, $K$ is a convolution matrix and thus can be easily diagonalized by fast transforms (with proper boundary conditionsassumed on $u$ ). Therefore, the main difficulty in solving ( ) is caused by the nondifferentiability and the universal coupling of variables of the TV term. Our algorithm is derived fromthe well-known variable-splitting and penalty techniques in optimization. First, we introduce an auxiliary variable ${\mathbf{w}}_i\in {\mathbb{R}}^2$ at pixel $i$ to transfer $D_{i}u$ out of the nondifferentiable term $\parallel ·\parallel$ . Then we penalize the difference between ${\mathbf{w}}_i$ and $D_{i}u$ quadratically. As such, the auxiliary variables ${\mathbf{w}}_i$ 's are separable with respect to one another. For convenience, in the following we let $\mathbf{w}\triangleq [{\mathbf{w}}_1,...,{\mathbf{w}}_{n^2}]$ . The approximation model to ( ) is given by