An approximate inverse preconditioner for Toeplitz systems with multiple right-hand sides

Applied Mathematics and Computation 218 (2012) 11370–11379

Contents lists available at SciVerse ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

An approximate inverse preconditioner for Toeplitz systems with multiple right-hand sides Jie Huang ⇑, Ting-Zhu Huang School of Mathematical Sciences, University of Electronic Science and Technology of China, Chengdu, Sichuan 611731, PR China

a r t i c l e

i n f o

Keywords: Preconditioner Toeplitz matrix Multiple right-hand sides Inverse formula

a b s t r a c t We are interested in preconditioning techniques for symmetric positive definite, but ill-conditioned, Toeplitz systems with multiple right-hand sides TX ¼ B. Here the righthand sides are available simultaneously. Our main idea is to use an approximate inverse of T as a preconditioner for T. The approximate inverse of T is generated by the inverse formula given by Lv and Huang [X.-G. Lv, T.-Z. Huang, A note on inversion of Toeplitz matrices, Appl. Math. Lett., 20 (2007) 1189–1193]. The construction of the preconditioner is of low computational cost, requires only the entries of T and does not require explicit knowledge of the generating function f of T. We show that if given a proper parameter, then the approximate inverse preconditioner approximates T 1 accurately, and thus the spectra of these preconditioned matrices are clustered around 1. It follows that the conjugate gradient method converges very quickly when applied to solve the preconditioned systems. Numerical results are given to demonstrate the effectiveness of our preconditioner. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction An n-by-n matrix T n with the entries t ij is said to be a Toeplitz matrix if it is constant along its diagonals; i.e. tij ¼ tij . In a variety of applications in mathematics and engineering [1], such as signal and image processing, we need to solve Toeplitz systems of the same coefficient matrix and different right-hand sides. When all the right hand sides are available simultaneously, the problem we are concerned with can be expressed as

T n X ¼ B; where T n is an n  n Toeplitz matrix, B ¼ ½b1 ; b2 ; . . . ; bs  and X ¼ ½x1 ; x2 ; . . . ; xs  are rectangular matrices of size n  s with s  n. Here the matrix T n is assumed to be symmetric positive definite. In particular, we are concerned with the development of preconditioners for conjugate gradient (CG) method. A circulant matrix is a special form of a Toeplitz matrix where each row of the matrix is a circular shift of its preceding row. The idea of using the preconditioned conjugate gradient (PCG) method with circulant preconditioners for solving Toeplitz systems was proposed by Strang [2] and Olkin [3] independently. Then circulant preconditioners have been intensively studied by a number of researchers, see for instance Ref. [1] and the references therein. Particularly, Chan et al. [4] constructed the Jackson kernel-type preconditions by convolving the generating function f of T n with the generalized Jackson kernels. This class of preconditioners works quite well for f having finitely many zeros and requires only the entries of the given

⇑ Corresponding author. E-mail addresses: [email protected] (J. Huang), [email protected] (T.-Z. Huang). 0096-3003/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.05.017

J. Huang, T.-Z. Huang / Applied Mathematics and Computation 218 (2012) 11370–11379

11371

Toeplitz matrix. Other preconditioners for Toeplitz systems include fast trigonometric transform-based preconditioners (see [5–7] for instance), Hartley transform-based preconditioners (see [8], for instance), and banded preconditioners (see [9–12], for instance). Recently, sparse approximate inverse preconditioners have been studied by a number of researchers; see Kolotilina and Yeremin [13], Benzi et al. [14], Benzi and Tuma [15], Tang [16], Chow [17] and Nikishin and Yeremin [18]. Based on the celebrated Gohberg–Semencul formula [19], Ng et al. [20] proposed to using recursive-based PCG method for solving Toeplitz systems. Afterwards, Ching et al. [21] proposed block diagonal and Schur complement preconditioners for block-Toeplitz systems. Moreover, Lin et al. [22] proposed factorized banded inverse preconditioners for matrices with Toeplitz structure. Furthermore, Lin and Ching [23] presented inverse Toeplitz preconditioners for Hermitian Toeplitz systems. In this paper, we construct an approximate inverse preconditioner, based on the inverse formula given by Lv and Huang [24], for Toeplitz systems with multiple right-hand sides (RHSs). In fact, the inverse formula given by Lv and Huang is a new circulant version of the Gohberg–Semencul formula. It is worth noting that the first circulant version of the Gohberg–Semencul formula was proposed and analyzed by Ammar and Gader [25]. However, unlike the construction of inverse preconditioners recursively by the Gohberg–Semencul formula as in Ng et al. [20], we generate our approximate inverse preconditioner with the two solutions of two fundamental equations directly, according to the inverse formula given by Lv and Huang. We show that if given a proper parameter, then the approximate inverse preconditioner approximates T 1 n accurately, and thus the spectra of these preconditioned matrices are clustered around 1. Therefore the PCG method when applied to solve these systems converges very quickly. It should be noted that to construct our approximate inverse preconditioner, we preliminarily solve two fundamental equations by the PCG method with Jackson kernel preconditioner of order 2 (K m;4 , denoted in the same way as in [4]). It implies that the constructing of our preconditioner needs some additional computational cost. However we will show numerically that the cost for the construction of the inverse preconditioner is acceptable, even preferable for Toeplitz systems with the same coefficient matrix and multiple RHSs. Generally, in our tests, for s > 3, our preconditioner is always superior than the Jackson kernel preconditioner of order 2, especially for very ill-conditioned cases. In addition, to further compare the proposed preconditioner with K m;4 , we use a preconditioned version of the global least squares (Gl-LSQR) method by Toutounian and Karimi [26] to solve T n X ¼ B simultaneously. Again, numerical examples demonstrate the effectiveness of the proposed preconditioner. The outline of the paper is as follows. In Section 2, we introduce the idea of the construction of the approximate inverse preconditioner, and then briefly discuss the spectrum of the preconditioned matrix. Also, the outline of the Gl-LSQR method is given. In Section 3, we give numerical examples to demonstrate the effectiveness of the proposed preconditioner. Finally, concluding remarks are given in Section 4.

2. Approximate inverse preconditioner In this section, we introduce the general idea of the construction of our inverse preconditioner for Toeplitz systems with multiple RHSs. Indeed, for a Toeplitz matrix, a number of reconstruction techniques of inverse have been proposed in the literature (see [19,27–31], for instance). Particularly, in [19], a nice matrix representation of the inverse well-known as ‘‘Gohberg–Semencul formula’’ was presented. The idea of using the Gohberg–Semencul formula recursively to construct approximate inverse preconditioners has been developed by Ng et al. [20]. In this paper, we aim to use the Jackson kernel preconditioner of order 2 for two fundamental equations, and then employ the two obtained numerical solutions to generate an approximate inverse directly according to the inverse formula given by Lv and Huang. 2.1. Approximate inverse preconditioner Lv and Huang [24] proved that the invertibility of a Toeplitz matrix could be determined through the solvability of two standard equations (the so-called fundamental equations):

T n u ¼ f0 and

T n v ¼ e1 ; where

f0 ¼ ð0; t 1n  t 1 ; . . . ; t2  t n2 ; t 1  tn1 ÞT ; e1 ¼ ð1; 0; . . . ; 0ÞT : Moreover, once we obtain the vector

u ¼ ðu1 ; u2 ; . . . ; un ÞT

11372

J. Huang, T.-Z. Huang / Applied Mathematics and Computation 218 (2012) 11370–11379

and the vector

v ¼ ðv 1 ; v 2 ; . . . ; v n ÞT ; then the inverse of T n can be constructed via

T 1 n ¼ C 1 U1 þ C2 U2 ;

ð1Þ

where

0

v1 vn

B Bv B 2 C1 ¼ B B .. @ .

v1 ..

.

... .. . .. .

v2 1

.. C . C C C; A vn C

vn

...

u1

un

B Bu B 2 C2 ¼ B B .. @ .

u1 .. .

1 . . . u2 .. . C . .. C C C; .. C . un A

un

...

u2

0

v2 v1

0 B B B U1 ¼ B B @

1 un 1

0 B B B and U 2 ¼ B B @

1 . . . u2 .. .. C . . C C C; .. C . un A 1

0 v n 0

u1

1 . . . v 2 C .. ... C . C C: .. C . v n A 0

Based on the above results, we partition our preconditioning technique into two steps: ^ and v ^ of the two standard equations, respectively. Step 1: Obtain the numerical solutions u ^ of Tn X ¼ B with the proposed approximate inverse preconditioner T^ 1 constructed by Step 2: Find the finest solution X n (1). ^ ^ It is worth noting that the approximate inverse preconditioner T^ 1 n need not be computed explicitly, and obtaining u and v ^ ¼ ð^ requires Oðn log nÞ operations by the PCG method with the preconditioner K m;4 . Then, to get X x1 ; ^ x2 ; . . . ; ^ xs Þ with our preconditioner T^ 1 n in Step 2, we solve each T n xi ¼ bi by the PCG method, or solve T n X ¼ B simultaneously by the preconditioned global least squares (PGl-LSQR) method [32]. We will show numerically that, in both cases, the cost of Step 1 is always ^ Hence, numerically, to facilitate convergence of CG or acceptable, even preferable, to obtain the approximate solution X. Gl-LSQR iterations for Toeplitz systems with multiple RHSs, our approximate inverse preconditioner is superior than circulant preconditioners. To solve T n X ¼ B by the above two methods with the proposed preconditioner, we need to consider the spectrum of the preconditioned matrix T^ 1 n T n . In fact, we have the following theorem, which is given in [24, Theorem 2]. Theorem 1. Let T n be a nonsingular Toeplitz matrix and be well-conditioned. Then the formula

T 1 n ¼ C 1 U1 þ C2 U2 ^ and v ^ are perturbed by the normwise relative errors bounded by e: is forward stable. More precisely, if the computed vectors u

^ k2 6 kuk2 ð1 þ eÞ; kv^ k2 6 kv k2 ð1 þ eÞ: ku Then with machine precision ~e, we have

   ^ 1 1      T n  T n  pffiffiffi 1    2 6 nð2e þ n~eÞ 1 þ 2 T n  kf0 k2 þ ~e n:  1  2 T n  2

Thus a simple, but very important, result now follows readily in the following theorem. Theorem 2. Let T n be a nonsingular Toeplitz matrix and T^ 1 n is constructed by means of the inverse formula (1). Then if we choose the value of e properly, the spectrum of the preconditioned matrix T^ 1 n T n is clustered around 1. Proof. By virtue of Theorem 1, we have

   ^ 1 1      T n  T n  pffiffiffi 1    2 6 nð2e þ n~eÞ 1 þ 2 T n  kf0 k2 þ ~e n;  1  2 T n  2

thus,

       ^ 1     ^ 1 1 1  T n T n  In  ¼ T^ 1 n T n  T n T n  6 T n  T n   kT n k2 2

2

2

     pffiffiffi   ~ 6 KðT n Þ  nð2e þ n~eÞ 1 þ 2T 1 n  kf0 k2 þ e n : 2

J. Huang, T.-Z. Huang / Applied Mathematics and Computation 218 (2012) 11370–11379

11373

Here KðT n Þ is the condition number of T n . Therefore, clearly, if T n is well-conditioned, then the spectrum of T^ 1 n T n is clustered around 1 for smaller e. Similar result holds also for ill-conditioned cases when a proper (often very small) value of e is chosen. h Generally speaking, the closer T^ 1 n is to the exact inverse of T n , the higher rate of convergence of iterative methods will be. That is, the smaller value of e is, the more clustered the spectrum will be, and also the more operations will be required to obtain ^ and v ^ in Step 1. Moreover, at each PCG step for the finest solution ^ u xi (i ¼ 1; . . . ; s), computing ðC 1 U 1 þ C 2 U 2 Þr requires six 2nlength FFTs and six n-length FFTs. Thus the total cost of the approximate inverse preconditioner for a Toeplitz system with multiple RHSs is Oðsn log nÞ operations. It is known that the cost of a 2n-length FFT is roughly the twice the cost of an n-length FFT, hence the cost per iteration in solution phase is about 8/3 times that required by circulant preconditioners. While, it is worth stressing that the effectiveness of a preconditioning strategy is strongly problem- and architecturedependent. For instance, a preconditioner which is expensive to compute may become viable if it is to be reused many times, since in this case the initial cost of forming the preconditioner can be amortized over several linear systems. Thus, when a Toeplitz system is very ill-conditioned and has several right-hand sides, setting e small enough to ensure the clustered ^ ^ spectrum of the preconditioned matrix T^ 1 n T n may become acceptable, although it may cost more to obtain u and v of higher accuracy. Numerical tests show that e ¼ 105 is appropriate, whether the matrix T n is well-conditioned or mildly illconditioned. For very ill-conditioned Toeplitz systems, e ¼ 106 provides a better choice. In order to confirm the above statements by several simple numerical tests in Section 3, we need the following definitions. A 2p-periodic continuous real-valued function f is called the generating function of T n if

tk ¼

1 2p

Z p

f ðhÞeikh dh;

k ¼ 0; 1; 2; . . . ; ðn  1Þ:

p

Clearly t k ¼ tk for all k. It follows that T n is a symmetric matrix. Moreover, if f is an even function with the essential infimum fmin > 0, then the matrix T n is symmetric positive definite. We recall that h0 is a zero of f of order q if f ðh0 Þ ¼ 0 and q is the smallest positive integer such that f ðqÞ ðh0 Þ – 0 and f ðqþ1Þ ðhÞ is continuous in a neighborhood of h0 [1]. We remark that the condition number of T n generated by such a function f grows like Oðnq Þ; see [31]. Because of the ill-conditioning, the PCG method will converge slowly and the number of iterations required for convergence grows like Oðnq=2 Þ. We note that the construction of our inverse preconditioner requires only the entries of T n and does not require explicit knowledge of the generating function f of T n . 2.2. Global least squares method The LSQR algorithm [32] is an iterative method for solving real linear systems of the form

Ax ¼ b; where A is a nonsymmetric matrix of order n and x; b 2 Rn . Recently, Toutounian and Karimi [26] provide a global least squares (Gl-LSQR) algorithm for solving linear systems with multiple RHSs

AX ¼ B; where B and X are n  s matrices. It generalizes the classical LSQR method to handle the linear systems with the same coefficient matrix and different right-hand sides. Chung et al. [33] then applied the Gl-LSQR algorithm for solving the image restoration problem. The procedure of the Gl-LSQR method is written as follows. The Gl-LSQR algorithm 1. 2. 3. 4.

Set X 0 ¼ 0ns .     b1 ¼ kBkF ; U 1 ¼ B=b1 , a1 ¼ AT U 1  ; V 1 ¼ AT U 1 =a1 . F ^ ^ 1 ¼ a1 . Set W 1 ¼ V 1 ; /1 ¼ b1 ; q For k ¼ 1; 2; . . . until convergence Do: i. xk ¼ AV k  ak U k ; bkþ1 ¼ kxk kF ; U kþ1 ¼ xk =bkþ1 . ii. sk ¼ AT U kþ1  bkþ1 V k ; akþ1 ¼ ksk kF ; V kþ1 ¼ sk =akþ1 . ^ 2k þ b2kþ1 Þ1=2 ; ck ¼ q ^ k =qk ; sk ¼ bkþ1 =qk . iii. qk ¼ ðq ^ kþ1 ¼ ck akþ1 . iv. hkþ1 ¼ sk akþ1 ; q ^ k; / ^k. ^ kþ1 ¼ sk / v. /k ¼ ck / vi. X k ¼ X k1 þ ð/k =qk ÞW k . vii. W kþ1 ¼ V kþ1  ðhkþ1 =qk ÞW k . ^ kþ1 j is small enough, then stop. viii. If j/

EndDo.

11374

J. Huang, T.-Z. Huang / Applied Mathematics and Computation 218 (2012) 11370–11379

Here, we use an approximate inverse T^ 1 n to precondition Toeplitz systems with multiple RHSs T n X ¼ B on the right. That ^ ^ 1 is to solve the preconditioned systems T n T^ 1 n Y ¼ B by the Gl-LSQR method, and then compute X ¼ T n Y. Notice that the approximate inverse preconditioner T^ 1 need not be computed explicitly as well, and T^ 1 multiply a vector requires six n n 2n-length FFTs and six n-length FFTs. Hence, in each iteration, compared with circulant preconditioners, our approximate inverse preconditioner will cost additional 30s FFTs. 3. Numerical examples In this section, we apply the two methods in Section 2 to solve symmetric positive definite Toeplitz systems with multiple RHSs T n X ¼ B for different generating functions. The experiments are run in MATLAB (version 7.7.0) with a machine precision 2:2204  1016 . The machine used is a Pentium (R) Dual-Core 2.60 Ghz personal computer with 6G memory. 3.1. Preconditioned CG method 3.1.1. Single RHS We first apply the PCG method to solve Toeplitz systems with single right-hand side

T n x ¼ b; for n ¼ 8192. Clearly, s ¼ 1. The coefficient matrices are generated by f ðxÞ ¼ x4 þ 1; jxj3 þ 0:01; x2 , x2 ðp4  x4 Þ; x4 ðp2  x2 Þ and T x4 , respectively. The right-hand side b 2 Rn is calculated such that x ¼ T 1 n b ¼ ½1; . . . ; 1 . All runs are initiated from the initial vector xð0Þ ¼ 0 and terminated if the current iteration satisfies

kb  T n xðkÞ k2 6 e: kb  T n xð0Þ k2 ^ and v ^. Here and in the after, we choose the finest level e ¼ 107 for x , and the coarser level e ¼ s0 :¼ 105 for u We note that the first two test functions: f ðxÞ ¼ x4 þ 1 and jxj3 þ 0:01 in Table 1 are positive functions and therefore correspond to well-conditioned systems. The two test functions x2 and x2 ðp4  x4 Þ are nonnegative functions with single or multiple zeros of order 2 on ½p; p. Thus the condition numbers of the Toeplitz matrices are growing like Oðn2 Þ, and hence the number of CG iterations required for convergence without using any preconditioner is increasing like OðnÞ. When the order of the zero is 4, like the two test functions x4 ðp2  x2 Þ and x4 in Table 1, the condition numbers of the resulting Toeplitz matrices will increase like Oðn4 Þ, and the matrices will be very ill-conditioned even for moderate n. In Table 1, we list the number of iterations (denoted by ‘‘IT p ’’), the elapsed CPU time in seconds (denoted by ‘‘CPU p ’’), and the corresponding relative accuracy:

k^x  x k2 kx k2 (denoted by ‘‘ReERR’’) required for convergence by using T^ 1 n ðs0 Þ as the preconditioner for solving n-by-n Toeplitz systems for different generating functions; see the columns under T^ 1 n ðs0 Þ. In addition, the number of iterations and the CPU timings re^ and v ^ in Step 1 are given; see the column under ‘‘IT c ’’ and ‘‘CPU c ’’, respectively. quired to get the two approximate solutions u For comparison purposes, we also give the number of iterations, the CPU timings and the accuracy of the PCG method with no preconditioner (I) and the Jackson kernel preconditioner of order 2 (K m;4 ); see the columns under ‘‘IT’’, ‘‘CPU’’ and ‘‘ReERR’’, respectively. The double asterisk ⁄⁄ signifies that more than 1000 iterations are required. In addition, we denote the CPU time and the accuracy by ‘‘–’’ if the solution obtained in 1000 PCG iterations is inaccurate. Clearly, from Table 1, as predicted, compared with K m;4 , it takes many additional iterations and much CPU time to construct the inverse preconditioner T^ 1 n ðs0 Þ; see the columns under ‘‘IT c ’’ and ‘‘CPU c ’’, respectively. Once the construction is finished, the number of iterations and the CPU time to get the approximate solutions with T^ 1 n ðs0 Þ is less than those with K m;4

Table 1 Behavior of PCG algorithm for various generating functions when

s0 ¼ 105 ; n ¼ 8192 and s = 1.

f ðxÞ

I IT

CPU

ReERR

IT

CPU

ReERR

IT p

IT c

CPU p

CPU c

ReERR

x4 þ 1

75

0.47

2.13e08

5

0.34

5.21e-11

2

12

0.06

0.91

3.57e-14

jxj3 þ 0:01

404

2.64

6.44e-08

5

0.47

8.48e-09

2

15

0.06

1.22

1.34e-12

x2

⁄⁄

–

–

7

1.03

5.51e-08

2

15

0.06

1.56

5.96e-07

⁄⁄

–

–

7

0.70

8.17e-08

2

22

0.06

2.55

1.40e-08

x ðp  x Þ

⁄⁄

–

–

21

2.67

0.7859

3

150

0.09

18.75

0.0172

x4

⁄⁄

–

–

26

3.11

0.7921

4

186

0.14

21.97

0.1141

x2 ðp4  x4 Þ 4

2

2

K m;4

T^ 1 n ðs0 Þ

11375

J. Huang, T.-Z. Huang / Applied Mathematics and Computation 218 (2012) 11370–11379 Table 2 Behavior of PCG algorithm for various generating functions when

s0 ¼ 105 ; n ¼ 8192 and s ¼ 60.

f ðxÞ

I IT

CPU

ERR

IT

CPU

ERR

IT p

IT c

CPU p

CPU c

x4 þ 1

4860

12.54

1.66e-08

412

39.70

2.26e-08

60

12

0.97

1.12

1.57e-07

27648

84.99

1.69e-08

498

58.94

5.37e-08

60

14

1.03

1.65

7.00e-07

3

jxj þ 0:01

K m;4

T^ 1 n ðs0 Þ ERR

x2

⁄⁄

–

–

744

90.84

3.69e-08

102

15

1.67

2.06

1.38e-06

x2 ðp4  x4 Þ x4 ðp2  x2 Þ

⁄⁄ ⁄⁄

– –

– –

1001 12804

131.81 1808.68

4.43e-08 6.51e-03

120 180

22 157

1.72 2.39

3.00 24.59

6.19e-08 4.59e-03

x4

⁄⁄

–

–

22660

2890.65

0.0521

360

219

34.40

30.10

0.0474

Fig. 1. Curves of IT by PCG iterations versus the number of the RHSs for various generating functions when n ¼ 8192 and

s0 ¼ 105 .

for all problems to a great extent; see the columns under ‘‘IT p ’’, ‘‘IT’’, ‘‘CPU p ’’ and ‘‘CPU’’, respectively. Hence it is naturally hoped that the additional computation cost is acceptable to handle the multiple RHSs cases.

11376

J. Huang, T.-Z. Huang / Applied Mathematics and Computation 218 (2012) 11370–11379

Fig. 2. Curves of CPU by PCG iterations versus the number of the RHSs for various generating functions when n ¼ 8192 and

s0 ¼ 105 .

Fig. 3. Curves of the relative residual norm versus the iteration number for PGl-LSQR method for f ðxÞ ¼ x4 þ 1 (left) and jxj3 þ 0:01 (right) when n ¼ 8192 and s0 ¼ 105 .

3.1.2. Multiple RHSs Based on the previous claim, we now consider Toeplitz systems with multiple right-hand sides

T n X ¼ B;

11377

J. Huang, T.-Z. Huang / Applied Mathematics and Computation 218 (2012) 11370–11379

where X ¼ ðx1 ; . . . ; xs Þ; B ¼ ðb1 ; . . . ; bs Þ; n ¼ 8192 and s ¼ 60. The coefficient matrices are also generated by f ðxÞ ¼ x4 þ 1; jxj3 þ 0:01; x2 ; x2 ðp4  x4 Þ, x4 ðp2  x2 Þ and x4 , respectively. Here we solve each T n xi ¼ bi by the PCG method. The initial guess is the zero vector and each right-hand side bi ¼ randðn; 1Þ (in MATLAB notation). In Table 2, we list the number of iterations, the CPU time in seconds and the corresponding relative error

^  Bk kT n X 2 kBk2 of the three preconditioners: I; K m;4 and T^ 1 n ðs0 Þ, respectively. Notice that we adopt the above relative error instead of that in ^ see the column under ‘‘ERR’’. Other symbols such as ⁄⁄, ‘‘–’’, ‘‘IT’’, Table 1 to estimate the accuracy of the computed solution X; ‘‘CPU’’ hold the same meaning as those in Table 1. From Table 2, we see that the PCG method with no preconditioner does not work for f ðxÞ ¼ x2 ; x2 ðp4  x4 Þ; x4 ðp2  x2 Þ and x4 . Even though the CG method works, the iterations converge very slowly; see the cases of f ðxÞ ¼ x4 þ 1 and jxj3 þ 0:01.

Fig. 4. Curves of the relative residual norm versus the iteration number for PGl-LSQR method for f ðxÞ ¼ x2 when n ¼ 8192 and

s0 ¼ 105 .

Fig. 5. Curves of the relative residual norm versus the iteration number for PGl-LSQR method for f ðxÞ ¼ x2 ðp4  x4 Þ when n ¼ 8192 and

s0 ¼ 105 .

Fig. 6. Curves of the relative residual norm versus the iteration number for PGl-LSQR method for f ðxÞ ¼ x4 ðp2  x2 Þ (left) and x4 (right) when n ¼ 8192 and s1 ¼ 106 .

11378

J. Huang, T.-Z. Huang / Applied Mathematics and Computation 218 (2012) 11370–11379

Clearly, the PCG with K m;4 and T^ 1 n ðs0 Þ converge very quickly with comparable accuracy to the exact solution for all cases. For the convergent iterations, the overall iterations (IT p þ IT c ) of T^ 1 n ðs0 Þ are often much smaller than iterations (IT) of K m;4 , especially for the very ill-conditioned cases f ðxÞ ¼ x4 ðp2  x2 Þ and x4 . The total CPU timings (CPU p þ CPU c ) by T^ 1 n ðs0 Þ are also less than those of K m;4 to a great extent. As expected, the cost for constructing the inverse preconditioner is acceptable, even preferable for multiple RHSs. Hence, it is numerically evident that the proposed approximate inverse preconditioner is much more robust, stable and efficient than the Jackson kernel preconditioner of order 2 when they are used to precondition Toeplitz systems with multiple RHSs. 3.1.3. Variant number of RHSs The above observations are further confirmed by the curves of the number of iterations versus the number of the RHSs, and the curves of the CPU timings versus the number of the RHSs plotted in Figs. 1 and 2 for various generating functions, respectively. When n ¼ 8192; s0 ¼ 105 , and s ranging from 1 to 60, the two figures clearly show that both the number of iterations and the CPU timings (including the cost to construct the inverse preconditioner) by T^ 1 n ðs0 Þ are much less than those by K m;4 , whether the coefficient matrix T n is well-conditioned, mildly ill-conditioned, or very ill-conditioned. 3.2. Preconditioned Gl-LSQR method The numerical results in this subsection concern the preconditioned Gl-LSQR (PGl-LSQR) method for solving systems with multiple RHSs simultaneously, with no preconditioner (I), the Jackson kernel preconditioner of order 2 (K m;4 ), and the approximate inverse preconditioner (T^ 1 n ðsÞ). We plot the curves of the relative residual norm (in 10 logarithm):

log10

^  Bk kT n X 2 kBk2

versus the iteration number for PGl-LSQR method for various generating functions in Figs. 3–6. Evidently, the PGl-LSQR method with K m;4 converges very slowly for all problems, similarly to the case without preconditioning technique. However, ^ for all problems. We PGl-LSQR with the inverse preconditioner converges very fast and accurately to the exact solution X ^ ^ note that the iteration number for T^ 1 ð s Þ includes the iterations for u and v . Hence, numerically, the inverse preconditioner n is much more efficient than the Jackson kernel preconditioner for the PGl-LSQR method. 4. Conclusions In this paper, we have presented an approximate inverse preconditioner for symmetric positive definite, but illconditioned, Toeplitz systems with multiple right-hand sides. The construction of our preconditioner is based on the inverse formula given by Lv and Huang [24]. We have numerically tested several generating functions, no matter the resulting coefficient matrices are well-conditioned, mildly ill-conditioned, or very ill-conditioned. The numerical experiments have shown that the approximate inverse preconditioner outperforms the Jackson kernel preconditioner of order 2 on aspects of solution accuracy, computation time and number of iterations when they are used to accelerate the CG method and Gl-LSQR for Toeplitz systems with multiple right-hand sides. Acknowledgements The authors would like to express their thanks to the reviewers and editor-in-chief Dr. Melvin Scott for their much constructive, detailed and helpful advice regarding revising this paper. This research is supported by NSFC (60973015, 61170311), Chinese Universities Specialized Research Fund for the Doctoral Program (20110185110020), Sichuan Province Sci. & Tech. Research Project (12ZC1802). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

R.H. Chan, M.K. Ng, Conjugate gradient methods for Toeplitz systems, SIAM Rev. 38 (1996) 427–482. G. Strang, A proposal for Toeplitz matrix calculations, Stud. Appl. Math. 74 (1986) 171–176. J. Olkin, Linear and nonlinear deconvolution problems, Ph.D. Thesis, Rice University, Houston, TX, 1986. R.H. Chan, A.M. Yip, M.K. Ng, The best circulant preconditioners for Hermitian Toeplitz systems, SIAM J. Numer. Anal. 38 (2001) 876–896. E. Boman, I. Koltracht, Fast transform based preconditioners for Toeplitz equations, SIAM J. Matrix Anal. Appl. 16 (1995) 628–645. T. Kailath, V. Olshevsky, Displacement structure approach to discrete-trigonometric-transform based preconditioners of G. Strang type and of T. Chan type, Calcolo 33 (1996) 191–208. S. Serra, A Korovkin-type theory for finite Toeplitz operators via matrix algebras, Numer. Math. 82 (1999) 117–142. D. Bini, P. Favati, On a matrix algebra related to the discrete Hartley transform, SIAM J. Matrix Anal. Appl. 14 (1993) 500–507. R. Chan, Toeplitz preconditioners for Toeplitz systems with nonnegative generating functions, IMA J. Numer. Anal. 11 (1991) 333–345. R. Chan, P. Tang, Constrained minimax approximation and optimal preconditioners for Toeplitz systems, Numer. Algorithms 5 (1993) 353–364. D. Noutsos, P. Vassalos, New band Toeplitz preconditioners for ill-conditioned symmetric positive definite Toeplitz systems, SIAM J. Matrix Anal. Appl. 23 (2002) 728–743. S. Serra, Optimal quasi-optimal and superlinear band-Toeplitz preconditioners for asymptotically ill-conditioned positive definite Toeplitz systems, Math. Comput. 66 (1997) 651–665.

J. Huang, T.-Z. Huang / Applied Mathematics and Computation 218 (2012) 11370–11379

11379

[13] L.Y. Kolotilina, A.Y. Yeremin, Factorized sparse approximate inverse preconditionings I: Theory, SIAM J. Matrix Anal. Appl. 14 (1993) 45–58. [14] M. Benzi, C.D. Meyer, M. Tuma, A sparse approximate inverse preconditioner for the conjugate gradient method, SIAM J. Sci. Comput. 17 (1996) 1135– 1149. [15] M. Benzi, M. Tuma, A sparse approximate inverse preconditioner for nonsymmetric linear systems, SIAM J. Sci. Comput. 19 (1998) 968–994. [16] W.P. Tang, Toward an effective sparse approximate inverse preconditioner, SIAM J. Matrix Anal. Appl. 20 (1999) 970–986. [17] E. Chow, A priori sparsity patterns for parallel sparse approximate inverse preconditioners, SIAM J. Sci. Comput. 21 (2000) 1804–1822. [18] A. Nikishin, A.Y. Yeremin, Prefiltration technique via aggregation for constructing low density high-quality factorized sparse approximate inverse preconditionings, Numer. Linear Algebra Appl. 10 (2003) 235–246. [19] I. Gohberg, A. Semencul, On the inversion of finite Toeplitz matrices and their continuous analogs, Mat. Issled. 7 (1972) 201–233. [20] M.K. Ng, H.W. Sun, X.Q. Jin, Recursive-based PCG methods for Toeplitz systems with nonnegative generating functions, SIAM J. Sci. Comput. 24 (2003) 1507–1529. [21] W.K. Ching, M.K. Ng, Y.W. Wen, Block diagonal and Schur complement preconditioners for block-Toeplitz systems with small size blocks, SIAM J. Matrix Anal. Appl. 29 (2007) 1101–1119. [22] F.R. Lin, M.K. Ng, W.K. Ching, Factorized banded inverse preconditioners for matrices with Toeplitz structure, SIAM J. Sci. Comput. 26 (2005) 1852– 1870. [23] F.R. Lin, W.K. Ching, Inverse Toeplitz preconditioners for Hermitian Toeplitz systems, Numer. Linear Algebra Appl. 12 (2005) 221–229. [24] X.-G. Lv, T.-Z. Huang, A note on inversion of Toeplitz matrices, Appl. Math. Lett. 20 (2007) 1189–1193. [25] G. Ammar, P. Gader, A variant of the Gohberg–Semencul formula involving circulant matrices, SIAM J. Matrix Anal. Appl. 12 (1991) 534–540. [26] F. Toutounian, S. Karimi, Global least squares method for solving general linear systems with several right-hand sides, Appl. Math. Comput. 178 (2006) 452–460. [27] G. Heinig, K. Rost, Algebraic methods for Toeplitz-like matrices and operators, in: Operator Theory: Advances and Applications, Birkhäuser, Basel, 1984. [28] A. Ben-Artzi, T. Shalom, On inversion of Toeplitz and close to Toeplitz matrices, Linear Algebra Appl. 75 (1986) 173–192. [29] M.K. Ng, K. Rost, Y.W. Wen, On inversion of Toeplitz matrices, Linear Algebra Appl. 348 (2002) 145–151. [30] G. Heinig, On the reconstruction of Toeplitz matrix inverses from columns, Linear Algebra Appl. 350 (2002) 199–212. [31] S. Serra, On the extreme eigenvalues of Hermitian (block) Toeplitz matrices, Linear Algebra Appl. 270 (1998) 109–129. [32] C.C. Paige, M.A. Saunders, LSQR: an algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Softw. 8 (1) (1982) 43–71. [33] S.Y. Chung, S.Y. Oh, S.J. Kwon, Restoration of blurred images by global least squares method, Chungcheong Math. Soc. 22 (2009) 177–186.