Cosine Transform Based Preconditioners for Total ... - Semantic Scholar

Cosine Transform Based Preconditioners for Total Variation Deblurring Raymond H. Chan, Tony F. Chan, Chiu-Kwong Wong

Abstract | Image reconstruction is a mathematically illposed problem and regularization methods are often used to obtain a reasonable solution. Recently, the total variation (TV) regularization, proposed by Rudin, Osher and Fatemi (1992), has become very popular for this purpose. In a typical iterative solution of the nonlinear regularization problem, such as the xed point iteration of Vogel or Newton's method, one has to invert linear operators consisting of the sum of two distinct parts. One part corresponds to the blurring operator and is often a convolution; the other part corresponds to the TV regularization and resembles an elliptic operator with highly varying coecients. In this paper, we present a preconditioner for operators of this kind which can be used in conjunction with the conjugate gradient method. It is derived from combining fast transform (e.g. cosine-transform based) preconditioners which the authors had earlier proposed for Toeplitz matrices and for elliptic operators separately. Some numerical results will be presented. In particular, we will compare our preconditioner with a variant of the product preconditioner proposed by Vogel and Oman [28].

remedy this ill-conditioning, several approaches of regularization are used. The Total Variation (TV) regularization method of Rudin, Osher and Fatemi [25], [24] is extremely eective approach for restoring edges of the original image. They consider solving the following constrained minimization problem: min u

Z



jruj dx dy subject to kHu ? z kL2( ) =  (2)

where j
j denotes the R Euclidean norm and  is the noise level. The quantity jruj dx dy is called the total variational norm of u. Instead of solving the constrained problem, Vogel considered the following closely-related regularization problem: 1 kHu ? z k2 +  Z jruj dx dy; (3) min f ( u ) = min L2 ( ) u u 2

I. Introduction see [1], [27]. Here  is a positive parameter which measures In this paper, we apply conjugate gradient precondition- the trade o between a good t and an oscillatory solution. ers to the iterative solution of some large-scale image pro- At a stationary point of (3), the gradient of f vanishes, cessing problems. The quality of the recorded image is giving: usually degraded by blurring and noise. The recorded im  age z and the original image u are often related by the g(u)  H (Hu ? z ) ? r
ru = 0; (x; y) 2 ; (4) equation, jruj

z (x; y) = Hu(x; y) + (x; y)



Z



h(x ? s; y ? t)u(s; t) dt ds + (x; y);

(1) see [12]. Here H denotes the blurring operator for the blurring function h, and  denotes the noise function. The image restoration problem is to obtain a reasonable approximation of the original image. Note that the problem Hu = z is ill-posed and the discretization matrix of H is usually ill-condition and hence the problem is extremely sensitive to noise. Thus, we cannot neglect the eect of noise and simply solve Hu = z . To

@u = 0; (x; y) 2 @ : @n

The Rsecond term in g is obtained by taking the gradient of  jruj dx dy and then applying integration by parts from which Neumann boundary condition results. We remark that the Euler-Lagrange equation for (2) also has a form similar to (4). Due to the term 1=jruj, (4) is a degenerate nonlinear second order diusion equation. The degeneracy can be removed by adding a positive parameter to jruj; see [27]. More precisely, if we let  (u) = p 1 2 ; Lu v = ?r
( (u)rv) (5) jruj + and Au v  (H H + Lu )v; then (4) becomes the following non-degenerate system Au u = H z; (x; y) 2 ; (6) @u with @n = 0; (x; y) 2 @ :

Raymond H. Chan is with Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong Kong (http://www.math.cuhk.edu.hk/~ rchan). Research supported in part by CUHK DAG 2060110 and Summer Research Grant. Tony F. Chan is with the Department of Mathematics, University of California, Los Angeles, 405 Hilgard Avenue, Los Angeles, CA 90095-1555 (http://www.math.ucla.edu/~ chan). Supported by grants ONR-N00017-96-1-0277 and NSF DMS-9626755. Part of this work was performed during a visit to the Department of Mathematics at the CUHK. Chiu-Kwong Wong is with Department of Mathematics, University of California, Los Angeles, 405 Hilgard Avenue, Los Angeles, CA A recent survey of related PDE approach to image analysis 90095-1555 (http://www.math.ucla.edu/~ ckwong). can be found in [2].

In [27], Vogel introduced the \lagged diusivity xed point iteration", which we denote by FP, to solve the system (6). If Au , H and Lu denote respectively the discretization matrices of Au , H and Lu , then the FP iteration will produce a sequence of discrete approximations fuk g to the solution u and can be expressed as Au uk+1  (H  H +Lu )uk+1 = H  z; k = 0; 1; ::: : (7) Note that obtaining uk+1 from the solution of uk requires solving a linear system with coecient matrix H  H +Lu . For the image restoration problem in (1), H corresponds to a discretization of the convolution operator, and often H will be a Toeplitz matrix. Thus, the coecient matrix in (7) corresponds to a sum of a convolution operator and an elliptic operator. We emphasize that it is not easy to devise fast iterative algorithms to solve this linear system. For example, the technique of applying multigrid method to solve such linear system is not yet well developed, see [9]. Vogel and Oman [28] has recently proposed using a \product" preconditioner for (7) which allows the deblurring part H  H and the PDE part Lu to be preconditioned separately. An alternative approach to solving the gradient equation (7) is to directly solve the minimization problem (3) by non-smooth optimization techniques; see for example [21], [22]. In this work, we apply the preconditioned conjugate gradient (PCG) method to solve (7) and we concentrate on nding a good preconditioner for (7). Given a matrix A, there are two criteria for choosing a preconditioner for A; see [18]. First, a preconditioner should be a \good" approximation to A. Secondly, it must be easily invertible. Recall that Au corresponds to a sum of a convolution operator and an elliptic operator. For matrices arising from elliptic boundary value problem, a \good" preconditioner must retain the boundary condition of the given operator [23]. Based on this idea, optimal sine transform preconditioners were constructed [15] for elliptic problems with Dirichlet boundary condition. If the boundary is rectangular, it was proved [15], [29] that the convergence rate of the PCG method with this preconditioner is independent of the grid size. In our present problem, Neumann boundary condition is imposed, see (4). Since the discrete Laplacian on a unit square with Neumann boundary conditions can be diagonalized by the discrete cosine transform matrix, this motivates us to use the optimal cosine transform approximation [11] to Lu as a preconditioner for the elliptic part in (7). In addition, R. Chan, Ng and Wong [14] applied the sine transform approximation to construct preconditioners for Toeplitz systems. It gives rise to fast convergence of the PCG method. It has been shown in [20] that the PCG method with the optimal cosine transform approximation can also produce the same good convergence result. By combining the results mentioned in the previous two paragraphs, we here propose a preconditioner for the system (7) by taking the sum of the cosine transform approximations to the matrices H H and Lu separately. The resulting approximation can still be diagonalized by the k

k

k

k

k

k

k

k

k

k

k

discrete cosine transform matrix and therefore easily invertible. This preconditioner was originally proposed in [10] where preliminary results and numerical experiments were given. In this paper, we will give a detailed discussion and analysis of the preconditioner, including a comparison of our preconditioner with the product preconditioner proposed by Vogel and Oman [28]. The outline of the paper is as follows. In the next section, we will de ne and construct the optimal cosine transform approximation for a general matrix. In Section III, we will use the approximation to construct a preconditioner for the system (7). In Section IV, we will introduce Vogel and Oman's product preconditioner and some of its possible variants. In the nal section, numerical performance of the preconditioner will be presented. II. Optimal Discrete Cosine Transform Preconditioner

The concept of optimal transform approximation was rst introduced by T. Chan [16]. Since preconditioners can be viewed as approximations to the given matrix An , it is reasonable to consider preconditioners which minimize jjBn ? An jj over all Bn belonging to some class of matrices and for some matrix norm jj
jj. T. Chan [16] proposed optimal circulant preconditioner that is the minimizer of the Frobenius norm kBn ? An kF over the class of all circulant matrices Bn . These preconditioners have been proved to be very ecient preconditioners for solving Toeplitz systems with the PCG method; see [5]. Analogously, R. Chan, Ng and Wong [14] de ned the optimal sine transform preconditioner to be the minimizer of kBn ? An kF over all matrices Bn which can be diagonalized by the discrete sine transform. They proved that for a large class of Toeplitz system, the PCG method with the sine transform preconditioner converges at the same rate as the optimal circulant one. Following the same approach, we will construct in Section II-A the optimal cosine transform preconditioner for general matrices. We remark that although the derivation is almost the same as that of the sine transform preconditioner, we present it here for the seek of charity. An alternative derivation using displacement structure can be found in [20]. The preconditioner will be applied to precondition both H H and Lu in (7) separately. For a survey on fast transform type preconditioners, we refer the reader to [8]. k

A. Construction of One-dimensional Preconditioner Let us denote Cn to be the n-by-n discrete cosine transform matrix. If ij is the Kronecker delta, then the (i; j )th entry of Cn is given by

r





2 ? j1 cos (2i ? 1)(j ? 1) ; 1  i; j  n; n 2n

see Sorensen and Burrus [26, p.557]. We note that the

Cn 's are orthogonal, i.e. Cn Cnt = In . Also, for any nvector v, the matrices-vector multiplication Cn v and Cnt v can be computed in O(n log n) real operations; see [30].

An

general Toeplitz banded

cost of constructing c(An ) O(n2 ) O(n) O((bl + bu )n) TABLE I

Cost of Constructing c(An ).

Ann , we denote (Ann )i;j;k;l to be the (i; j )th entry of the (k; l)th block of Ann . Let R be a permutation matrix which simply reorders Ann in another coordinate direction. More precisely, R satis es (Rt Ann R)i;j;k;l = (Ann )k;l;i;j ; 1  i; j  n; 1  k; l  n: Then the Level-2 cosine transform preconditioner c2 (Ann ) for Ann is de ned by

c2 (Ann ) = Rc1 (Rt c1 (Ann )R)Rt : Let Bnn be the vector space containing all matrices that can be diagonalized by Cn . More precisely, It is easy to show that the preconditioner c2 (Ann ) can Bnn = fCn n Cnt j n is an n?by?n be diagonalized by Cn Cn : real diagonal matrixg: c2 (Ann ) = (Cn Cn ) diag((Cn Cn )t Ann (Cn Cn )) (Cn Cn )t : For an n-by-n matrix An , we choose our preconditioner c(An ) to be the minimizer of jjBn ? An jjF in the Frobenius Hence, c2 (Ann ) can be inverted easily. Note that in the norm in the space Bnn . Following the terminology used construction of our preconditioner, we will use the Levelin T. Chan, the approximation is called the optimal cosine 2 approximation instead of the Level-1. It is because the transform preconditioner for An and denoted by c(An ). It Level-1 preconditioner has a relatively expensive initializacan be shown that c(An ) is linear and preserves positive tion cost for block Toeplitz matrices with Toeplitz block, de niteness; see [19]. see [7]. We will show in Appendix A that c(An ) can be obFor elliptic problem, it can be proved that the Level-2 optained optimally in O(n2 ) operations for general matrices. timal cosine transform preconditioner c2 (Ann ) is a \good" The cost can be reduced to O(n) operations when An is a preconditioner. banded matrix or a Toeplitz matrix. This is the same cost Theorem 1: Let Ann be the 5-point centered discretizaas that for constructing the optimal circulant precondition- tion matrix of ers. We recall that in our case, Lu and H are banded matrix and Toeplitz matrix respectively. We summarize the ?(a(x; y)ux)x ? (b(x; y)uy )y + !u = f (x; y) on [0; 1]2 construction cost of c(An ) in Table I. In Table I, we denote with homogeneous Neumann boundary condition. Assume bl and bu the lower and upper band width of An . ! > 0 and the mesh is uniform with size 1=n. Then we B. Construction of Two-dimensional Preconditioner have (c2 (Ann )?1 Ann )  ( ccmax )2 For 2D n  n images, the matrices H  H and Lu in (7) min are block matrices of the following form: where 0 < cmin  a(x; y); b(x; y)  cmax . 0 A A ::: A 1 Proof: We refer the reader to Appendix B. 1;1 1;2 1;n B A2;1 A2;2 : : : A2;n C Optimal cosine transform preconditioner can also be C B Ann = B shown to be good for solving Toeplitz system. For one@ ... . . . . . . ... CA : dimensional cases, i.e. point-Toeplitz matrix systems, the An;1 An;2 : : : An;n convergence proof has been established in [20]. For 2dimensional cases, we can prove the following result. Here Ai;j are square matrices of order n. Theorem 2: Let Ann be a block Toeplitz matrix with In [17], T. Chan and Olkin proposed the Level-1 and Assume the underlying generating seLevel-2 circulant preconditioners for such block matrices. Toeplitz (blocks. Following their approach, we de ne the Level-1 and Level-2 quence akj) of Ann is absolutely summable, i.e., cosine transform preconditioners for Ann . The idea of the 1X 1 X Level-1 and Level-2 preconditioners is to approximate the ja(kj) j  G < 1 matrix Ann in one direction and two directions respectively. j =0 k=0 Speci cally, the Level-1 preconditioner c1 (Ann ) is de ned by and 0 c(A ) c(A ) : : : c(A ) 1 1;1 1;2 1;n a(kj) = (Ann )(j?1)n+k;1 for 1  j; k  n: B C c ( A ) c ( A ) : : : c ( A ) 2 ; 1 2 ; 2 2 ;n c1 (Ann ) = B B C: .. C ... ... If min(Ann )   > 0 where  is independence on n, then @ ... . A for all  > 0, there exists N > 0 such that for all n > N , at c(An;1 ) c(An;2 ) : : : c(An;n ) most O(n) of the eigenvalue values of the preconditioned To de ne the Level-2 cosine transform preconditioner, let matrix c2?1 (Ann )Ann have absolute values larger than . us rst give some notations. For any n2 -by-n2 block matrix Thus the PCG method converges in at most O(n) steps. k

Proof: The proof is similar to that of Corollary 1 in [7] and will be omitted. We remark that the numerical results in [7] for the Level-2 optimal circulant preconditioner is better than the theoretical result and the numbers of iterations do not grow with n. With Theorems 1 and 2, it seems reasonable to use the Level-2 cosine approximation to construct preconditioners for the linear system (7). III. Cosine Transform Preconditioner for TV denoising and deblurring

A straightforward preconditioner for Au is k

vector multiplication
?1=2 v can be computed in O(n2 ) operations because
?1=2 is a diagonal matrix. Since Lu is banded, Lu v can also be done in O(n2 ) operations. For H being a block Toeplitz matrix with Toeplitz blocks, Hv can be calculated in O(n2 log n) operations; see [7]. Therefore, the matrix vector multiplication can be done in O(n2 log n) operations. The system MD y = b can be solved in O(n2 log n) operations by exploiting the Fast Cosine Transform. Therefore, the total cost of each PCG iteration is bounded by of O(n2 log n). k

k

IV. The Product Preconditioner of Vogel and Oman

c2 (Au ) = c2 (H  H + Lu ) = c2 (H  H ) + c2 (Lu ): In this section, we introduce the product preconditioner proposed in Vogel and Oman in [28] and discuss some of However, computing c2 (H  H ) according to the formula in its possible variants. The product preconditioner for the Corollary 1 in Appendix A, requires computing all the en- system (7) is de ned as tries of H  H and is costly. Another way is to approximate c2 (H  H ) by c2 (H ) c2 (H ). More precisely, a preconditioner P~ = 1 (H  H + I )1=2 (L + I )(H  H + I )1=2 : (10) u for Au in (7) can be de ned as

k

k

k

k

k

M = c2 (H ) c2 (H ) + c2 (Lu ): (8) Here, is a parameter and will be chosen to be p as suggested in [28]. Note that computing P~ ?1 v in each PCG One problem with the preconditioner M is that it does iteration requires solving a convolution problem (H  H + not capture the possibly large variation in the coecient I )?1=2 v and an elliptic problem (Lu + I )?1v which are of the elliptic operator in (7) caused by the vanishing of not straightforward. To make the product preconditioner jruj in (4). To cure this problem, we apply the technique P~ more practical, Vogel and Oman assumed that the blurof diagonal scaling. More precisely, if we denote (
) to be ring function h is periodic and hence (H  H + I )?1=2 v k

k

the spectral radius of a matrix and we de ne
 (H  H )I +  diag(Lu )

k

k

then we consider solving the equivalent system A~u u~k+1 = H~  z k

where A~u = H~  H~ + L~ u , H~ = H
?1=2 , L~ u =
?1=2 Lu
?1=2 and u~k =
1=2 uk . In summary, the Level2 cosine preconditioner with diagonal scaling is given by MD = H^  H^ + c2 (L~ u ) (9) k

k

k

k

k

where H^ = c2 (H )c2 (
?1=2 ). We note further that if 1 , 2 and 3 are respectively the eigenvalue matrices of c2 (H ), c2 (
?1=2 ) and c2 (L~ u ), then the preconditioner can be expressed as k

MD = (Cn Cn ) ( 1  2 + 3 ) (Cn Cn )t : 1

2

Hence, the preconditioner can be easily invertible. Finally, we comment on the cost of constructing MD and of each PCG iteration. We note that Lu is a sparse matrix with only ve nonzero bands. Also since H corresponds to a discretization of the convolution operator (1), H often is a block Toeplitz matrix with Toeplitz blocks. By using Table I, the construction cost of c2 (H ) and c2 (Lu ) can be shown to be O(n2 log n) operations; see [7], [17]. The cost of one PCG iteration is bounded by the cost of the matrix vector multiplication A~u v =
?1=2 (H  H + Lu )
?1=2 v and the cost of solving the system MD y = b. The matrix k

k

k

can be computed by the FFT. Moreover, when solving (Lu + I )?1v, it requires an inner PCG iteration in which a multigrid preconditioner is used, see [28]. More precisely, there are three nested iterations in their method: the outermost FP iterations, the inner PCG iterations with preconditioner P~ and the innermost PCG iterations to invert (Lu + I ). Since we do not make any assumptions on h and it would not be fair to compare only the inner PCG iteration numbers and ignore the work required in the innermost PCG iterations (i.e. the work on solving (Lu + I )?1 v which is substantial), this makes the comparison between the two preconditioners dicult. In the experiments below, we will instead use the following preconditioner for comparison: P = 1 (c2 (H ) c2 (H ) + I )1=2
(11) (c2 (Lu ) + I )
(c2 (H ) c2 (H ) + I )1=2 ;

k

k

k

k

p

where is again set to . This preconditioner is obtained by taking cosine transform approximations of the three factors in P~ . The resulting preconditioner P can be diagonalized by the 2-D cosine transform and hence solving P ?1 v at each PCG step requires just about the same cost as our preconditioner M . Moreover, there is no need for the innermost PCG iterations. We note that there is a connection between the product preconditioner P and the preconditioner M given in (8). First, note that the product preconditioner P in (11) can be viewed as an operator-splitted approximation of our preconditioner M in (8). Since the three factors in the right

hand side of (11) are commutative, we have

P ? M =  c2 (H ) c2 (H )c2 (Lu ) + I:

(12)

k

The right hand side matrix pwill be called the splitting error of P . If we again take = , the splitting error becomes p(c (H ) c (H )c (L ) + I ): 2 2 2 u Now, let us investigate the contribution from this splitting operator by using Fourier analysis. Denote H (f ) and L(f ) to be the spectrum of the blurring function h and the dierential operator Lu respectively. Then the splitting error in (12) will have spectrum p(jH (f )j2 L(f ) + 1); (13) where f denotes the frequency variable. In general, the second order dierential operator Lu has L(f ) proportional to f 2 . In fact Lu , being a high pass lter, has large eigenvalues corresponding to the large frequency modes and small eigenvalues corresponding to the low frequency modes (c.f. the case when Lu is the Laplacian operator). On the other hand, the blurring operator H being a low pass lter has the large eigenvalues corresponding to the low frequency modes and the small eigenvalues corresponding to the high frequency modes. Hence, we expect that the large eigenvalues of the dierential operator will be damped by the small eigenvalues of the blurring operator when they are multiplied together as in (13). In particular, we expect c2 (H ) c2 (H )c2 (Lu ) to be a bounded operator (Note that the cosine transform approximation c2 (
) only changes the boundary condition of the operators, and in general will not change the ordering of the eigenvalues). Thus, for small , the splitting error in (12) will be small and hence the performance of P and M should be about the same. However, since L(f ) is proportional to f 2 , when the spectrum of the blurring function H (f ) decays slowly, L(f ) may not be suciently damped. As a result, the splitting error in (12) will be large. Therefore, in this case, we expect M to outperform P for 's that are not very small. In particular, let us consider the limiting case, the denoising problem with H (f )  1 (i.e. h(x; y) is the delta function). In this case, M = c2 (Lu ) + I; p p P = 1 + (c (L ) + I ) = (1 + )( c (L ) + I ); k

k

k

k

k

k

k

and hence

2

2

uk

p

P ? M = (c2 (Lu ) + I );

uk

k

which is not small. We remark that even in this denoising p case, P 6= M . The reason is due to the choice of = . Of course, when = 1, the performance of the preconditioners P and pM are exactly the same (since P = 2M ). However, =  is the best choice over a wide range of problems, see [28]. In conclusion, our preconditioner M is good and robust for solving deblurring problems and we expect it to outperform the product preconditioner P especially when the

splitting error in (12) is large. That is the case when h is not a very signi cant blur or when h tends to the delta function (i.e., H tends to I ). This will be demonstrated in Section V-B when we compare the performance of P and M for various blurring functions h and 's. As in the case for the preconditioner M , we can apply a diagonal scaling to the matrix P to capture the possible large variation in the coecient matrix. One possible way is to de ne it analogous to (9), i.e. P = 1 (H^  H^ + I )1=2 (c (L~ ) + I )(H^  H^ + I )1=2 : D

2

uk

Similarly, we have

^ 2 (L~ u ) + I; PD ? MD =  H^  Hc k

which will be called the splitting error of PD . V. Numerical Results

In this section, we compare the numerical performance of dierent preconditioners in solving the linear system (7) for dierent images. A. Test Problem 1 In this test, the original image is given by

u(x; y) = [1=5;2=5][1=4;3=4] + [3=5;4=5][1=4;3=4] where (x; y) 2 [0; 1]  [0; 1] and [a;b] denotes the characteristics function for the interval [a; b]. The blurring function h in (1) is chosen to be a truncated Gaussian, given by h(x; y) =



e? (x2+y2 ) if jxj; jyj  1=4 : 0 otherwise

(14)

Here  is a parameter which controls the severity of the blurring. More precisely, the smaller  is, the more significant the blurring will be. In this experiment, we choose  = 200 so that h is a moderate blurring function. The noise function  has normal distribution and is scaled such that the noise-to-signal ratio (NSR), kkL2 =kHukL2 , is 0:5. The true image and the observed image are shown in Figure 3. We will perform the FP iterations, starting with u0 = z , until the gradient g in (4) satis es kg(uk )k2 = kg(u0)k2 < 10?3. We will apply the CG method to solve the linear system (7) and the initial guess for the CG method in the kth FP iteration is chosen to be the (k ? 1)th FP iterate. The iteration is stopped when the residual vector rk of the linear system (7) at the kth CG iteration satis es jjrk jj2 =jjr0 jj2 < 10?3. In our numerical experiment, we will focus on the performance of dierent choices of preconditioners for various of parameters n,  and . Here n is the number of pixels in one direction, i.e. the matrix Au is of size n2 -by-n2. Tables II and III show the number of iterations required for convergence of the FP iteration and the CG iteration for dierent choices of preconditioners and parameters. Note k

that the CG iteration numbers shown in Tables II and III are the average number of CG iterations per FP step. The symbol \#" denotes the number of iterations for FP. The notations I ,
denote respectively no preconditioner and the diagonal scaling preconditioner. Some of the data are plotted in Figures 1 and 2. We observe from Figure 1 that the MD and PD preconditioners require signi cantly fewer iterations than other preconditioners for all values of and . Moreover, we can observe that the smaller is, the more ill-conditioned the system is. From Figure 2, we observe that the number of iterations corresponding to I grows like O(n0:9 ), which from standard convergence theory for CG implies that (Au )  O(n1:8 ). If the preconditioner MD or PD is used, the number of iterations grows like O(n0:22 ) which implies that (MD?1 Au )  O(n0:44 ). However, the preconditioners
, M and P reduce the growth of the number of iterations only to O(n0:67 ), O(n0:58 ) and O(n0:54 ) respectively. Therefore, MD and PD as preconditioners are much more eective than other three preconditioners. From Table III, we observe that the performance of the preconditioners M and P (resp. MD and PD ) are almost the same. This is in agreement with our observation in Section IV that for a smooth kernal function h (as in our case, the Gaussian function in (14)), the splitting error P ? M will be small for small  (e.g. < 10?2 ) and hence we expect the performance of the preconditioners M and P to be close. In Figure 4, we show the recovered images for various . The smaller is, the closer the recovered image is to the true image. Figure 5 shows how the recovered images depend on the value of . In this case, the FP method produces the best image when  = 10?3. k

k



I

# 93 70 49 27

10?3 10?2 10?1 100

n = 15 M P MD PD


26 16 11 8

44 27 16 10

43 52 25 32 15 20 9 13 n = 31

24 15 10 7

27 18 12 10

# I
M P MD PD 10?3 112 108 42 72 77 23 21 10?2 101 62 26 41 44 16 14 10?1 20 36 18 24 25 11 10 8 100 52 21 13 13 14 8 n = 63 # I
M P M D PD 10?3 266 197 59 112 114 25 24 10?2 211 113 42 64 65 18 16 10?1 164 64 29 35 36 13 11 100 115 37 21 18 18 10 8 n = 127 # I
M P M D PD 10?3 420 > 400 99 169 169 32 29 10?2 340 213 67 93 94 22 19 10?1 257 118 47 49 50 15 13 8 100 177 69 34 25 25 9 TABLE II  = 10?3

+ = I,

x = Delta,

* = M,

solid = M_D,

o=P

dash = P_D

2.4

2.2

B. Test Problem 2 We will now perform two experiments to illustrate that there are situations where our preconditioner M is better and more robust than the product preconditioner P . In this test problem, we will basically repeat the experiment in Section V-A with two dierent type of blurring functions. In the rst experiment, we will choose the Gaussian blurring function in (14) with several parameters  , see Figure 6. We remark that the larger the  is, the less signi cant is the blurring (see Figure 7), and we expect that our preconditioner to give a better performance than the product preconditioner P . The parameter n is xed to 63, to 0:1 and NSR=0.5. We report the number of PCG iterations required for various preconditioners and 's in Table IV. We observe from Table IV that when  = 200, the performance of M and P are almost identical, as expected. When  is increased to 2000, the number of PCG iterations required by P are much more than that by M for  2 [10?2; 10?4 ]. As we mentioned in Section IV, when Fig. 1. Observation (top gure): the smaller the is, the more the blurring is not very signi cant (i.e., when h is close ill-condition of Au . Here n = 127 and  = 10?3 . Observation (bottom gure): For various , preconditioners MD and PD reto the delta function), the splitting error in (12) becomes quires signi cantly less PCG iterations compared to the others ? 5 large unless  is very small (e.g. 10 when  = 2000). We preconditioners. Here n = 127 and = 0:1. remark that the performance of MD and PD is still close log(no. of cg iterations)

2

1.8

1.6

1.4

1.2

1

0.8 −3

−2.5

+ = I,

−2

x = Delta,

−1.5 log β

* = M,

−1

solid = M_D,

−0.5

o=P

0

dash = P_D

log(no. of cg iterations)

2

1.5

1 −5

−4.5

k

−4

−3.5 log

α

−3

−2.5

−2

−1.5

+ = I,

x = Delta,

* = M,

solid = M_D,

o=P

dash = P_D

log(no. of cg iteration)

2

1 1

1.2

+ = I,

1.4

x = Delta,

1.6 log n

* = M,

1.8

solid = M_D,

2

o=P

2.2



# I
76 52 31 80 36 18 92 48 36 94 183 183



# I
137 98 48 164 64 29 104 64 38 214 115 102



# 211 257 181 136

10?2 10?3 10?4 10?5

dash = P_D

2.2

2

log(no. of cg iteration)

# I
36 25 19 49 16 11 32 37 36 20 123 125

10?2 10?3 10?4 10?5

1.5

1.8

1.6

10?2 10?3 10?4 10?5

1.4

1.2

1 1

1.2

1.4

1.6 log n

1.8

2

2.2

Fig. 2. We plot the no. of PCG iterations versus n in log scale for various preconditioners and . Here  = 10?3 is xed. Top : = 0:01, slopes of I ,
, M , P , MD , PD are 0:88, 0:67, 0:58, 0:54, 0:22, 0:22 respectively. Bottom: = 0:1, slopes of I ,
, M , P , MD , PD are 0:84, 0:68, 0:51, 0:49, 0:22, 0:19 respectively. Observation: Preconditioning by MD or PD require the smallest number of PCG iterations. The lines correspond to MD and PD have the smallest slope and therefore PCG has faster convergence rate. 5

5

10

10

15

15

20

20

25

25

30

n = 15 M P MD PD



10?2 10?3 10?4 10?5

I

189 118 109 107


73 47 36 75

14 20 15 20 21 22 15 17 n = 31

10 10 19 15

20 25 24 25 35 40 78 79 n = 63

13 11 27 74

29 30 35 36 35 36 72 73 n = 127

16 13 21 60

10 12 22 17

M P MD PD 12 10 30 76

M P M D PD 17 11 20 60

M P MD PD 37 49 50 54

TABLE III = 0:1

37 50 50 55

17 15 20 35

22 13 17 33

30 5

10

15

20

25

30

5

10

15

20

25

30

Fig. 3. True image (left) and the blurred noisy image, NSR=0.5 (right). 5

5

10

10

5

5

10

10

15

15

20

20

25

25

30

30 5

10

15

20

25

30

15

15

20

20

5

5

25

10

10

15

15

20

20

25

25

25

30

30 5

10

15

20

25

30

5

5

10

15

20

25

5

10

15

20

25

30

5

10

15

20

25

30

30

5

10

10 30

30 5

15

10

15

20

25

30

15

Fig. 5. Deblurred images for various : 10?5 (top left), 10?4 (top right), 10?3 (bottom left), 10?2 (bottom right). In this experiment, = 0:1, n = 31. Observation: when  is too small (e.g 10?4 ), the recovered images do not get enough regularization and become irregular. When  is large (e.g. 10?2 ), edges in the recovered images seem to be smooth out. In this case, the FP Fig. 4. Deblurred images for various : 0.01 (left top), 0.1 (right iteration produces the best image when  = 10?3 . top), 1 (bottom left), 10 (bottom right). In this experiment, ? 3  = 10 , n = 31. Observation: the smaller the is, the shaper the recovered image is . 20

20

25

25

30

30

5

10

15

20

25

30

5

10

15

20

25

30



10?2 10?3 10?4 10?5

# 137 164 104 214



#

 = 200 M P M D PD 29 35 35 72

30 16 36 13 36 21 73 60  = 2000

17 11 20 60

10?2 10?3 10?4 10?5

M P M D PD

10?2 111 31 42 10?3 61 31 47 10?4 136 46 55 10?5 82 96 92

18 17 29 83  = 20000



18 14 30 80

10?2 10?3 10?4 10?5

 # M P M D PD 10?2 25 23 37 17 16 10?3 19 14 36 9 12 10?4 6 5 22 4 10 10?5 3 2 13 2 6



10?2 10?3 10?4 10?5

TABLE IV

Gaussian blur, n = 63 and = 0:1

 1=( 2) px2 + y2 < 

17 12 23 99

PD 19 12 34 72

PD 19 15 17 13

TABLE V

10

10

10

20

20

20

30

30

30

40

40

40

50

50

60

50

60 10

20

30

40

50

60

60 10

20

30

40

50

60

10

20

30

40

50

60

Fig. 7. Blurred images when  = 200; 2000; 20000 (left to right).

:

elsewhere

0

PD

Out of focus blur, n = 63 and = 0:1

in this case. When we further increase  to 20000, both preconditioners M and MD outperform P and PD for a wide range of  (e.g. [10?2 ; 10?5]). In the second experiment, we choose the out of focus blurring functions. More precisely,

h(x; y) =

 = 0:1 # M P MD 141 28 29 16 146 36 36 13 90 38 38 24 90 118 110 100  = 0:05 # M P MD 126 32 35 18 86 35 37 14 101 56 66 34 44 78 80 69  = 0:01 # M P MD 31 26 34 19 31 18 35 13 10 7 23 5 5 3 14 3



In this case, the smaller  is, the less signi cant is the blurring, see Figures 9 and 10. From Table V, we can make the same observation as for the Gaussian blurring function. Namely, when the blurring is very signi cant (e.g. when  = 0:1), the performance of M and P are almost the same. However when the severity of blur decreases (e.g.  = 0:05 Fig. 8. Blurred and noisy images (observed images) when  = or 0:01), our preconditioners M and MD perform better. 200; 2000; 20000 (left to right) . In conclusion, the cosine transform preconditioners M and MD are more robust than the product preconditioners P and PD over a wide range of blurring functions. 10

10

10

20

20

20

30

30

30

40

40

40

50

50

50

60

60

10

10

10

10

20

20

20

30

30

30

40

40

40

50

60 20

30

40

50

60

40

50

60

60

10

20

30

40

50

60

10

10

20

20

20

30

30

30

40

40

40

50

50

50

60

60 20

30

40

50

60

10

20

30

40

50

60

10

20

30

40

50

60

60 10

20

30

40

50

60

Fig. 9. Out of focus PSF's when  = 0:1; 0:05; 0:01 (left to right).

50

60 10

30

10

10

50

20

60 10

20

30

40

50

60

10

20

30

40

50

60

Fig. 6. Gaussian PSF's when  = 200; 2000; 20000 (left to right).

C. Test Problem 3 We consider a 2D image restoration problem arising in ground-based atmospheric imaging. In this test, we will compare the quality of the recovered images for several

10

10

10

20

20

20

30

30

30

40

40

40

50

50

50

60

60 10

20

30

40

50

60

60 10

20

30

40

50

60

10

20

30

40

50

60

Fig. 10. Blurred images when  = 0:1; 0:05; 0:01 (left to right).

10

10

10

20

20

20

30

30

30

40

40

40

50

50

50

60

60 10

20

30

40

50

60

60 10

20

30

40

50

60

10

20

30

40

50

60

Fig. 11. Observed images when  = 0:1; 0:05; 0:01 (left to right). Fig. 14. Restored images for various FP iterations. Left: 1st FP with

kg(u)k = 1:1  10?2 . Middle: 3rd FP ?with kg(u)k = 3:3  10?3 . numbers of FP iterations. The problem consists of a 2563 . Here  = 10?4 . Right: 5th FP with k g ( u ) k = 1 : 6  10 by-256 image of an ocean reconnaissance satellite observed by a simulated ground-based imaging system together with a 256-by-256 image of a guide star observed under similar circumstances, see Figure 13. The data are provided by the Phillips Air Force Laboratory at Kirkland AFB, NM [4]. The imaging system detects the atmospheric distortions using a natural guide star image, see Figure 12 (right). A wavefront sensor measures the optical distortions which can then be digitized into a blurred image of the guide star pixel. We refer the reader to [13] on how to form the Restored images for various FP iterations. Left: 1st FP with blurring matrix H . In Figures 14-17, we present restored Fig. 15. kg(u)k = 1:4  10?3 . Middle: 3rd FP ?with kg(u)k = 3:4  10?4 . images for various values of the parameter  after one, three Right: 5th FP with kg(u)k = 1:7  10 4 . Here  = 10?5 . and ve FP iteration(s). Here, we x = 0:1 and use the MD preconditioner when solving the linear system (7). We perform the PCG iterations until the relative residual less than 10?2. The value kg(u)k after the FP iterations are presented. In Table VI, we compare the number of CG iterations with and without applying our preconditioner MD for various of 's. We note that the preconditioner MD can signifciantly speed up the convergence rate of the CG method. We observe from Figures 14-17 that after performing only Fig. 16. Restored images for various FP iterations. Left: 1st FP with kg(u)k = 2:7  10?4 . Middle: 3rd FP ?with kg(u)k = 5:4  10?5 . 3 FP iterations, we obtain very good recovered images and 5 . Here  = 10?6 . Right: 5th FP with k g ( u ) k = 2 : 9  10 most of the noise in the observed image are removed. Also, we note that when  is too large (e.g.  = 10?4), the recovered image looks \ at" and lost most of the features in the original image. When  got smaller (e.g.  = 10?5 or 10?6), an antenna appears in the recovered images. We remark that the nonsmoothless of the images appearing in Figure 17 is due to insucient regularization. Fig. 17. Restored images for various FP iterations: Left: 1st FP with kg(u)k = 1:9  10?4 . Middle: 3rd FP ?with kg(u)k = 1:1  10?5 . Right: 5th FP with kg(u)k = 6:8  10 6 . Here  = 10?7 . Fig. 12. Original image (left) and guide star image (right)



10?4 10?5 10?6 10?7

I

176 152 105 154

TABLE VI

Fig. 13. Observed image

MD 19 21 30 40

Average number of PCG iterations in the first five FP steps. The column \I" denotes no preconditioning.

Here ei denotes the ith unit vector in IRn . Note that Tn (
) Hn ((
)) are linear operators. Therefore, an n-byIn this paper, we propose cosine transform precondition- and n matrix to Bnn if and only if there exist ers for solving linear systems arising from total variation w ; : : : ; w B2n IRbelongs such that 1 n image deblurring problem. Our analysis and the numerical n n results indicate that the cosine transform preconditioner is X X an ecient and robust preconditioner over a wide class of wj [Tn (ej ) + Hn ((ej ))] Bn = wj Qj = j =1 j =1 image deblurring problems. VI. Conclusion Remarks

= Tn (

VII. Appendix A

In this appendix, we present how to construct c(An ) eciently, as stated in Table I in Section II-A. The approach is basically the same as that for the sine transform preconditioner, see [14]; we present it here for the sake of clairty. An alternative derivation can be found in [20]. An essential idea is to make use of a sparse and structured basis for Bnn , see [14]. Lemma 1: (Boman and Koltracht [3]) Let Qk , k = 1; : : : ; n, be n-by-n matrices with the (i; j )th entry given by 8 1 if ji ? j j = k ? 1, > < j = 2n ? k + 2, Qk (i; j ) = > 11 ifif ii + + j = k, : 0 otherwise. Then fQk gnk=1 is a basis for Bnn. In other words, every matrix in Bnn can be expressed as a linear combination of the n matrices fQk gnk=1 . We display the basis for the case n = 6.

01 0 Q1 = @ 00 0 0

0 1 0 0 0 0

0 0 1 0 0 0

0 0 0 1 0 0

0 0 0 0 1 0

0 0 0 0 0 1

0 0

1 0 0 1 0 0

1 0 0 0 1 0

0 1 0 0 0 1

0 0 1 0 0 1

0 0 0 1 1 0

1 0

0 0 1 0 0 1

0 1 0 0 0 1

1 0 0 0 1 0

1 0 0 1 0 0

0 1 1 0 0 0

00 1 Q3 = @ 01 00 0 Q5 = @ 10

1 01 A ; Q2 = @ 001 0 0

1 0 1 0 0 0

0 1 0 1 0 0

0 0 1 0 1 0

0 0 0 1 0 1

0 0 0 0 1 1

0

0 1 0 0 1 0

1 0 0 0 0 1

1 0 0 0 0 1

0 1 0 0 1 0

0 0 1 1 0 0

1

0 0 0 1 0 1

0 0 1 0 1 0

0 1 0 1 0 0

1 0 1 0 0 0

1 1 0 0 0 0

1 00 A ; Q4 = @ 110 0

1 00 A ; Q6 = @ 000 1

1 A;

1 A; 1 A:

n X j =1

wj ej ) + Hn ((

= Tn (w) + Hn ((w))

P

n X j =1

wj ej ))

with w = nj=1 wj ej . Now computing the optimal cosine transform approximation can be reformulated as solving the n-dimensional minimization problem, min

w=(w1 ;:::;wn )2IRn

kTn (w) + Hn ((w)) ? An kF :

The minimum can be calculated by setting

@ 2 @wi kTn (w) + Hn ((w)) ? An kF = 0; for i = 1; : : : ; n. For the sake of presentation, let us il-

lustrate the procedure of computing the minimum by considering the simple case n = 6. By the de nition of the Frobenius norm, we express kTn (w) + Hn ((w)) ? An k2F in terms of the matrix entries and then carry out the partial derivative with respect to wi . Then, we see that w satis es the following linear system

06 BB 2 BB 0 BB 2 @0

1

2 0 2 0 2 12 4 0 4 0 C C 4 12 4 0 4 C Cw 0 4 12 4 0 C C 4 0 4 12 4 A 0 2 0 4 0 a4 +12a + a BB a12 + a23 + a34 + a4511 + a5622 + a2133 B a13 + a24 + a35 + a46 + a31 + a42 = B BB a14 + a25 + a36 + a41 + a52 + a63 @ a15 + a26 + a51 + a62 + a14 + a23 a16 + a61 + a15 + a24 + a33 + a42 +a44 + a55 + a66 +a32 + a43 + a54 + a65 + a11 + a66 +a53 + a64 + a12 + a21 + a56 + a65 +a13 + a22 + a31 + a46 + a55 + a64 +a32 + a41 + a36 + a45 + a54 + a63 +a51 + a26 + a35 + a44 + a53 + a62

(15) In general, each Qk has at most 2n non-zero entries. In order to give a precise description of Bnn, we intro1 duce the following notations. t De nition 1: Let w = (w1 ; : : : ; wn ) be an n-vector. DeCC ne the shift of w as (w)  (w2 ; w3 ; : : : ; wn ; 0)t . De ne CC : (16) Tn (w) to be the n-by-n symmetric Toeplitz matrix with w CC as the rst column and Hn (w) to be the n-by-n Hankel A matrix with w as the rst column and (wn ; : : : ; w1 )t as the last column. Lemma 2: Bnn = fTn (w) + Hn ((w)) j w = We observe that the kth entry of the right hand side (w1 ; : : : ; wn )t 2 IRn g. vector in (16) is obtained by adding those aij for which the Proof: From (15), we observe that every Qi (i = (i; j )th position of Qk is nonzero. For general n, let rn be 1;


; 6) is a sum of a Toeplitz matrix and a Hankel matrix. an n-vector with the k-th component given by In fact, by using Lemma 1, it is not dicult to prove that X ai;j : (17) [rn ]k = Qi = Tn (ei ) + Hn ((ei )): (Q ) 6=0 k i;j

An

general Toeplitz banded

if n is even; and

cost of constructing rn O(n2 ) O(n) O((bl + bu )n)

w = 2n1 2 Rnt 2(n ? 1)U

+1 2 2nU n?2 1 ; n+1 2

?

TABLE VII

Cost of Constructing rn .

+ nI

+1 2

n

n

?2nU

+ ne1et1

+1 ; n?1 2 2 2nU n?2 1 + I n?2 1 n

!

Rn rn

if n is odd. Proof: Here we just give the proof for the case n is even. The proof for odd n is similar. To minimize kBn ? An k2F over Bnn, we set

If An has no special structure, then clearly by (17), rn can be computed in O(n2 ) operations because each Qi has only @ 2 O(n) non-zero entries. If An = [ai;j ] is a Toeplitz matrix @wi kc(An ) ? An kF = 0; for i = 1; : : : ; n: (correspond to H in (7)), then the sum in (17) can be computed without explicit addition because summing ai;j for We obtain a linear system that has the same structure as constant value of jj ? ij can be reduced to a scalar multi- that in (16). Permutating the system by Rn yields plication. Similarly, for banded matrix An with lower and  nV 2V U  upper band width bl and bu , the cost of forming rn can be 2 2U 2 V t 2nI 2 Rn w = Rn rn : reduced to O((bl + bu )n). We summarize the construction cost of rn in Table VII. Here V is an n2 -by- n2 matrix given by We now go back to the solution of the linear system V = 2I 2 ? e1 et1 : (16). We rst reorder the unknowns wi of w in such a way that the odd index entries and even index entries appear respectively in the upper half and lower half of the resulting By direct veri cation, vector. For simplicity, this leads to the following de nition. nV 2nU + nI + ne1 et1 + 2V U ?2nU  = 2n2 I ; 2 2 2 2 2 De nition 2: Let Rn be the n-by-n permutation matrix      with the (i; j )th entry given by 2U 2 V t ?2nU 2 + 2nI 2 2(n ? 1)U 2 + nI 2 = 2n2 I 2 ;     8 1 if 1  i  d n e and j = 2i ? 1, nV ?2nU 2 + 2V U 2 2(n ? 1)U 2 + nI 2 = 0; < 2 [Rn ]i;j = : 1 if d n2 e < i  n and j = 2i ? 2d n2 e, and 0 otherwise.     After permutation, (16) becomes a block system, 2U 2 V t 2nU 2 + nI 2 + ne1et1 + 2nI 2 ?2nU 2 = 0: n

n

n

n

n

n

n

n

n

n

06 B 0 B B 0 B B 2 B @2

0 0 2 2 2 12 0 4 4 4 0 12 4 4 4 4 4 12 0 0 4 4 0 12 0 2 4 4 0 0 12

10 w CC BB w13 CC BB w5 CC BB w2 A@ w 4 w6

1 CC CC = R r : CC 6 6 A

n

n

n

2

+ ne1 et1

?2nU 2

n

2

 ?2nU 2 2(n ? 1)U 2 + nI 2 Rn rn n

n

n

n

n

with 2

n

n

n

n

n

n

n

n

 nV 2V U  2  22UnU2 V t+ nI2nI+2 ne et

c(An ) = Tn (w) + Hn ((w))

n

n

n

n

Therefore, we get

The following theorem and corollary prove that in general if rn is known in advance, then the block system can be solved in O(n) operations. Theorem 3: Let An = [ai;j ] be an n-by-n matrix and c(An ) be the minimizer of kBn ? An kF over all Bn 2 Bnn . Denote Ui;j to be the i-by-j matrix with all its entries being one and if i is equal to j , then the matrix is just denoted by Ui . Let e1 to be the rst unit vector of length d n2 e. Then

w =2nU2n12 R+nt nI

n

n

n

n

2

?2nU 2

n

1 1

 ?2nU 2 2 2(n ? 1)U 2 + nI 2 = 2n In : n

n

n

Combining together with the fact that Rn is orthogonal, (18) follows. Before going on, let us rst emphasize the relationship between the rst column of matrices B 2 Bnn and their eigenvalues. For any matrix B 2 Bnn , we have B = Cn Cnt where  is the eigenvalue matrix of B . If D denotes the diagonal matrix whose diagonal is equal to the rst column of Cnt and 1n denotes the n-vector of all ones, then we have Cnt e1 = D1n . Therefore, the relation is given by

D?1 Cnt Be1 = 1n: In particular if the rst column of c(An ) is known, we can obtain the eigenvalue matrix  of c(An ) in O(n log n) operations. Hence, the matrix-vector multiplication c(An )?1 v can be computed as v ? Cnt v, v ? ?1 v, v ? Cn v which costs O(n log n) operations. The following corollary

(18) gives the explicit formula for the entries of the rst column

of c(An ). The proof follows directly from the expressions (18) and therefore we omit it. Corollary 1: Let An be an n-by-n matrix and c(An ) be the minimizer of kBn ? An kF over all Bn 2 Bnn. Denote by q the rst column of c(An ). If so and se are de ned respectively to be the sum of the odd and even index entries of rn , then we have, for n even, [q]1 = 2n1 2 (2n[rn ]1 + n[rn ]2 ? 2se ) [q]i = 2n1 2 (n[rn ]i + n[rn ]i+1 ? 2se ) i = 2; : : : ; n ? 1 [q]n = 2n1 2 (?2nso + (2n ? 2)se + n[rn ]n ) and for n odd, [q]1 = 2n1 2 (2n[rn ]1 + n[rn ]2 ? 2so) [q]i = 2n1 2 (n[rn ]i + n[rn ]i+1 ? 2so ) i = 2; : : : ; n ? 1 [q]n = 2n1 2 (?2nse + (2n ? 2)so + n[rn ]n ): From Corollary 1 and Table VII, we see that c(An ) can be obtained in O(n2 ) operations for general matrix An and O(n) operations for band matrix An and Toeplitz matrix An . VIII. Appendix B: Proof of Theorem 1

and (! ? cmax )I are positive de nite, using (19) and (20) we have  xt(A + cmaxI )x ! ? cmax  xt Ann x  max t (c2 (A) + cmin I )x ; ! ? cmin  xt c2 (Ann )x   xcmax t  max c xxt ((AAL ++ II))xx ; !! ?? ccmax

and

[1] [2] [3]

cminAL  c2 (A)  cmaxAL :

and

xt (c2 (A) + cmax I )x + (! ? cmax )xt x = xt c2 (Ann )x = xt (c2 (A) + cminI )x + (! ? cmin )xt x:

Since the inequality

a b  a b a + b  max ; min ; 

[10]

(19) [11]

(20)



min

 t I )x ; ! ? cmin   min xtx(c ((AA+) +cmin I )x ! ? cmax  c 2xt (A +cmax min L I )x ! ? cmin   min c xt (A + I )x ; ! ? c  max  max ?Lcmin  = min ccmin ; !! ? cmax max c min ?! c as  ?! 0:

Hence, Theorem 1 is proved.

[7] By noting that AL can be diagonalized by Cn Cn , we can [8] prove similar to [15, equation (17)] that [9]

xt (A + cminI )x + (! ? cmin)xt x = xt Ann x = xt (A + cmax I )x + (! ? cmax )xt x

L

max

cminAL  A  cmaxAL :

2

min

min

xt Ann x t x c2 (Ann )x

Let A and AL be the 5-point discretization matrices of ?(a(x; y)ux )x ? (b(x; y)uy )y and the negation of the [4] [5] Laplacian respectively, both with homogeneous Neumann boundary condition. It follows immediately that Ann = [6] A + !I: Similar to [6, (4.2)], we can prove that

! and x 2 IRn , we have Now for any 0 <  < cmax



cmax max ccmax ; !! ? ? cmin min c max ?! c as  ?! 0 =

[12] [13] [14]

[15] c d c+d c d holds for any a; b; c; d > 0 and the matrices A + cminI , [16] A + cmax I , c2 (A) + cmax I , c2 (A) + cminI , (! ? cmin)I

References R. Acar and C. Vogel, Analysis of Bounded Variation Penalty Methods, Inverse Problems, 10 (1994), pp. 1217{1229. L. Alvarez and J. Morel, Formalization and Computational Aspects of Image Analysis, Acta Numerica (1994), pp. 1{59. E. Boman and I. Koltracht, Fast Transform Based Preconditioners for Toeplitz Equations , SIAM J. Matrix Anal. Appl., 16, 628{645, 1995. R. Carreras, Personal Correspondence, 1993. R. Chan, The Spectrum of a Family of Circulant Preconditioned Toeplitz Systems , SIAM J. Numer. Anal., 26 (1989), pp. 503{ 506. R. Chan and T. Chan, Circulant Preconditioners for Elliptic Problems J. Numer. Linear Algebra Appls., 1 (1992), pp. 77{ 101. R. Chan and X. Jin, A Family of Block Preconditioners for Block Systems, SIAM J. Sci. Stat. Comput., 13 (1992), pp. 1218{1235. R. Chan and M. Ng, Conjugate Gradient Methods for Toeplitz Systems, SIAM Review, 38 (1996), pp. 427-482. R. Chan, T. Chan and W. Wan, Multigrid for DierentialConvolution Problems Arising from Image Processing, Proceedings of the Workshop on Scienti c Computing, Hong Kong, 97, Springer-Verlag, Ed: G. Golub, S.H. Lui, F. Luk and R. Plemmons. R. Chan, T.F. Chan, and C.K. Wong, Cosine Transform Based Preconditioners for Total Variation Minimization Problems in Image Processing, Iterative Methods in Linear Algebra, II, V3, IMACS Series in Computational and Applied Mathematics, Proceedings of the Second IMACS International Symposium on Iterative Methods in Linear Algebra, Bulgaria, June, 1995, pp. 311{329, Ed: S. Margenov and P. Vassilevski. R. Chan, W. Ching and C. Wong, Optimal Trigonometric Preconditioners for Elliptic and Queueing Problems, SEA Bull. Math., 20 (1996), pp. 117{124. R. Chan, J. Nagy and R. Plemmons, FFT-Based Preconditioners for Toeplitz-Block Least Squares Problems, SIAM J. Numer. Anal., 30 (1993), pp. 1740{1768. R. Chan, M. Ng and R. Plemmons, Generalization of Strang's preconditioner with applications to iterative deconvolution, Numer. Linear Algebra Appls., 3 (1996), pp. 45{64. R. Chan, M. Ng, and C. Wong, Sine Transform Based Preconditioners for Symmetric Toeplitz Systems, Linear Algebra Appls., 232 (1996), 237-260. R. Chan and C. Wong, Sine Transform Based Preconditioners for Elliptic Problems, accepted for publication by Numerical Linear Algebra and Applications. T. Chan, An Optimal Circulant Preconditioner for Toeplitz Systems , SIAM J. Sci. Stat. Comput., 9 (1988), pp. 766{771.

[17] [18] [19] [20]

[21] [22] [23] [24] [25] [26] [27] [28] [29] [30]

T. Chan and J. Olkin, Circulant Preconditioners for Toeplitzblock Matrices, Numer. Algorithms, 6 (1994) pp. 89{101. G. Golub and C. Van Loan, Matrix Computations , 2nd Ed., The Johns Hopkins University Press, Maryland, 1989. T. Huckle, Circulant/skewcirculant Matrices for Solving Toeplitz Matrix Problems, SIAM J. Matrix Anal. Appl., 13 (1992), pp. 767 { 777. T. Kailath and V. Olshevsky, Displacement structure approach to discrete-trigonometric-transform based preconditioners of G. Strang and T. Chan types, appears in the Italian Journal Calcolo, devoted to the Proceedings of the International Workshop on \Toeplitz matrices, algorithms and applications", Cortona, Italy, September 1996. I. Ito and K. Kunisch, An Active Set Strategy for Image Restoration Based on the Augmented Lagrangian Formulation, Preprint, Center for Research in Scienti c Computation, North Carolina State University. Y. Li and F. Santosa, An Ane Scaling Algorithm for Minimizing Total Variation in Image Enhancement, Technical Report CTC94TR201, Cornell Theory Center, Cornell University. T. Manteuel and S. Parter, Preconditioning and Boundary Conditions, SIAM J. Numer. Anal., 27 (1990), pp. 656{694. L. Rudin, S. Osher, Total Variation Based Image Restoration with Free Local Constraints, IEEE International Conference on Image Processing, Austin, TX, 1994. L. Rudin, S. Osher and E. Fatemi, Nonlinear Total Variation Based Noise Removal Algorithms, Physica D., 60 (1992), pp. 259{268. H. Sorensen and C. Burrus, Fast DFT and Convolution Algorithms , Handbook of Signal Processing edited by S. Mitra and J. Kaiser, New York, Wiley. C. Vogel and M. Oman, Iterative Methods for Total Variation Denoising, SIAM J. Sci. Stat. Comput., 17 (1996) pp. 227{238. C. Vogel and M. Oman, Fast Total Variation-Based Image Reconstruction, in the Proceedings of the 1995 ASME Design Engineering Conferences, V3, part C, pp. 1009-1015. C. Wong, Block Toeplitz Type Preconditioners for Elliptic Problem, M.Phil. thesis, The University of Hong Kong, 1994. P. Yip and K. Rao, Fast Decimation-in-time Algorithms for a Family of Discrete Sine and Cosine Transforms , Circuits, Syst., Signal Process, 3 (1984), pp. 387{408.