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

where $\beta \gg 0$ is a penalty parameter. It is well known that the solution of ( ) converges to that of ( ) as $\beta \to \infty$ . In the following, we concentrate on problem ( ).

Basic algorithm

The benefit of ( ) is that while either one of the two variables $u$ and $\mathbf{w}$ is fixed, minimizing the objective function with respect to the other has a closed-form formula that we willspecify below. First, for a fixed $u$ , the first two terms in ( ) are separable with respect to ${\mathbf{w}}_i$ , and thus the minimization for $\mathbf{w}$ is equivalent to solving

It is easy to verify that the unique solutions of ( ) are

where the convention $0·(0/0)=0$ is followed. On the other hand, for a fixed $\mathbf{w}$ , ( ) is quadratic in $u$ and the minimizer $u$ is given by the normal equations

By noting the relation between $D$ and $D_i$ and a reordering of variables, ( ) can be rewritten as

where

The normal equation ( ) can also be solved easily provided that proper boundary conditions are assumed on $u$ . Since both the finite difference operations and the convolution are notwell-defined on the boundary of $u$ , certain boundary assumptions are needed when solving ( ). Under the periodic boundary conditions for $u$ , i.e. the 2D image $u$ is treated as a periodic function in both horizontal and vertical directions, $D^{\left(1\right)}$ , $D^{\left(2\right)}$ and $K$ are all block circulant matrices with circulant blocks; see e.g. , . Therefore, the Hessianmatrix on the left-hand side of ( ) has a block circulant structure and thus can be diagonalized by the 2D discreteFourier transform $\mathbf{F}$ , see e.g. . Using the convolution theorem of Fourier transforms, the solution of( ) is given by

where the division is implemented by componentwise. Since all quantities but $w$ are constant for given $\beta$ , computing $u$ from ( ) involves merely the finite difference operation on $w$ and two FFTs (including one inverse FFT), once the constant quantities are computed.

Since minimizing the objective function in ( ) with respect to either $\mathbf{w}$ or $u$ is computationally inexpensive, we solve ( ) for a fixed $\beta$ by an alternating minimization scheme given below.

Algorithm :

• Input $f$ , $K$ and $\mu >0$ . Given $\beta >0$ and initialize $u=f$ .
• While“not converged”, Do
  • Compute $\mathbf{w}$ according to ( ) for fixed $u$ .
  • Compute $u$ according to ( ) for fixed $w$ (or equivalently $\mathbf{w}$ ).
• End Do

Optimality conditions and convergence results

To present the convergence results of Algorithm "Basic Algorithm" for a fixed $\beta$ , we make the following weak assumption.

Assumption 1 $\mathcal{N}\left(K\right)\cap \mathcal{N}\left(D\right)=\left\{0\right\}$ , where $\mathcal{N}(·)$ represents the null space of a matrix.

Define

Furthermore, we will make use of the following two index sets:

Under Assumption 1, the proposed algorithm has the following convergence properties.

Theorem 1 For any fixed $\beta >0$ , the sequence $\left\{(w^k,u^k)\right\}$ generated by Algorithm "Basic Algorithm" from any starting point $(w^0,u^0)$ converges to a solution $(w^{*},u^{*})$ of ( ). Furthermore, the sequence satisfies

• $\parallel w_E^{k+1}-w_E^{*}\parallel \le \sqrt{\rho \left({\left(T^2\right)}_{EE}\right)}\parallel w_E^k-w_E^{*}\parallel ;$
• $\parallel u^{k+1}-u^{*}{\parallel }_M\le \sqrt{\rho \left(T_{EE}\right)}{\parallel u^k-u^{*}\parallel }_M;$