A Fast Algorithm for Solving Special Tridiagonal Systems

9 Springer-Verlag 1994 Printed in Austria

Short Communication/Kurze Mitteilung

9 A Fast Algorithm for Solving Special Tridiagonal Systems

W . - M . Y a n a n d K.-L. C h u n g , T a i p e i

Received June 18, 1993; revised November 11, 1993

Abstract - - Zusammenfassung

A Fast Algorithm for Solving Special Tridiagonal Systems. In this paper, a fast algorithm for solving the special tridiagonal system is presented. This special tridiagonal system is a symmetric diagonally dominant and Toeplitz system of linear equations. The error analysis is also given. Our algorithm is quite competitive with the Gaussian elimination, cyclic reduction, special LU factorization, reversed triangular factorization, and Toeplitz factorization methods. In addition, our result can be applied to solve the near-Toeplitz tridiagonal system. Some examples demonstrate the good efficiency and stability of our algorithm.

AMS Subject Classification: 65F

Key words: Diagonally dominant matrices, error analysis, linear recurrences, Toeplitz matrices, tri- diagonal matrices.

Ein schneller Algorithmus zur Li~sung spezieller tridiagonaler Systeme. In dieser Arbeit wird ein schneller Algorithmus zur L6sung symmetrischer, diagonaldominanter tridiagonaler T6pflitz-Systeme vorgestellt. Auch eine Fehleranalyse liegt vor. Der Algorithmus ist den folgenden Verfahren mindestens gleichwertig: Gauss-Elimination, zyklische Reduktion, spezielle LU-Faktorisierung, umgekehrte Faktorisierung, T6plitz-Faktorisierung. AuBerdem kann unser Vorgehen zur L6sung in tridiagonalen fast-TSplitz- Systemen verwendet werden. Einige Beispiele zeigen die Effizienz und Stabilit/it unseres Algorithmus.

1. I n t r o d u c t i o n C o n s i d e r a n n x n s y s t e m of l i n e a r e q u a t i o n s w h e r e

ax = b,

(1)

A = 7

a n d Ifl] > 2171. T h a t is, t h e m a t r i x A is strictly d i a g o n a l l y d o m i n a n t . S p e c i a l tri- d i a g o n a l s y s t e m s of e q u a t i o n s h a v i n g s y m m e t r i c c i r c u l a n t coefficient m a t r i c e s

appear in many applications ([8], [-9], [-11], E14], [131). Many methods have been

proposed for solving such systems. These methods are Gaussian elimination [41,

cyclic reduction [-9], special

L U

factorization (El4], 1-10]), reversed triangular

factorization ([5], [61, [7]), and Toeplitz factorization with Sherman-Morrison

formula ([81, [12], ~3]). The interested readers are suggested to consult the survey

paper by Boisvert [2].

In this paper, a fast algorithm for solving the special tridiagonal system is presented.

Our algorithm consists of three phases. The first phase is a Toeplitz factorization

of a slightly perturbed system, which is similar to the method of Fischer et al. E8].

The second phase uses a forward and backward substitution procedure to solve the

perturbed problem. In the third phase the solution to the original problem is

recovered from the solution to the perturbed problem; this is called the update

procedure. The major contributions of our algorithm are twofold: (1) based on the

diagonally dominant property, a new but simple update procedure is proposed and

(2) a detailed error analysis is given. Our efficient and stable algorithm is quite

competitive with previous methods [2]. In addition, our result can be applied to

solve near-Toeplitz tridiagonal systems. It also can be applied to solve the quadratic

B-spline curve fitting problem [11] and a parabolic PDE [131.

We begin by presenting our algorithm for solving the symmetric Toeplitz tridiagonal

system. We then give an extension to the symmetric circulant tridiagonal case.

2. Symmetric Toeplitz Tridiagonal Systems

Throughout this paper, matrices are represented by uppercase letters, vectors by

bold lowercase letters, and scalars by plain lowercase letters. The superscripts T

correspond to the transpose operation.

We first consider to solve the symmetric Toeplitz tridiagonal system Bx = b, where

and

I/~l > 2171.

Let d =

-~(Idl > 2)

and 7

d 1

B' = = L U , 1 d 1

(2)

w h e r e L =

1

]

- b 1 a n d - b 1 - b 1 U =

i

a , a 1 a

which implies t h a t a - b - - - d a n d - a b = 1. This, in t u r n implies t h a t a = (d +_ xfld 2 - 4)/2. Since we wish the m a t r i x U to be d i a g o n a l l y d o m i n a n t , we will select the sign so t h a t the a b s o l u t e v a l u e of a is g r e a t e r t h a n 1. T h a t is, w h e n d > 2, we chose a = (d + x / d 2 - 4)/2 ( > d/2 > 1), f r o m which

b = a - d = ( - d + @ - 4)/2.

W h e n d < - 2, we chose a = (d - @ - 4)/2 ( < d/2 < - 1), t h e n b = a - d = ( - d - x / d 2 - 4)/2.

Since

labl

= 1 a n d lal > 1, we h a v e t h a t LbL < 1, t h a t is, the m a t r i x L is d i a g o n a l l y d o m i n a n t . T h e c o m p u t a t i o n of a a n d b p r o v i d e s the T o e p l i t z f a c t o r i z a t i o n of the m a t r i x B', which c a n be d o n e in 0 ( 1 ) time. 13 By (2), d = 2 , a n d a - b = d, it is clear t h a t 7 B = 7B' +

i_ b

(3)

T o solve Bx = b, we first solve B ' x ' = l b ( = b') using a f o r w a r d a n d b a c k w a r d s u b s t i t u t i o n p r o c e d u r e (second p h a s e in o u r a l g o r i t h m ) . I t is n o t h a r d to verify t h a t the n u m b e r of f l o a t i n g - p o i n t o p e r a t i o n s r e q u i r e d in this p r o c e d u r e is a b o u t 5n. By (3) a n d b = 7b', we o b t a i n B x ' = 7 B ' x ' -- 7bx'le 1 = b - ? b x ' l e l , (4) w h e r e e l = (1,0 . . . . ,0)% N o t e t h a t (4) c a n b e d e r i v e d f r o m S h e r m a n - M o r r i s o n f o r m u l a [1] a n d the e q u a l i t y B ' - l e l = 7(1 -- b e ~ B ' - l e l ) B - l e ~ w h i c h c a n be d e r i v e d with s o m e effort. Therefore, the s o l u t i o n x to be d e t e r m i n e d is x' + 7 b x ' ~ B - l e ~ . If we can find a v e c t o r p such t h a t Bp is e q u a l to a m u l t i p l e of e l , t h e n we c a n solve B x = b a p p r o x i m a t e l y b y a d d i n g a m u l t i p l e of p to x'.

T o find p, we first solve the r e c u r r e n c e relation: 7Pi-1 + flPi -4- ypi+~ = 0 for 2 _< i n - i. F r o m a - b = d, - ab = 1, a n d d = -fl, it follows b 2 q- db + 1 = 0 a n d 7 b2 +

7

fib + 7 = 0. Hence, b is a z e r o o f the characteristic p o l y n o m i a l of the a b o v e recur- rence relation. N a t u r a l l y , if we t r y p = (b, b 2 . . . b t, 0 . . . . ,0) r, t h e n n - - t 1 B p = (db + b 2,b + db 2 + b a , . . . , b t-1 + db t,b t,O . . . O) r 7 t v" J u - - ' , r - ~ ) t n - - t = (db

+ b 2, 0,..., O,

b t-1 -t- dbt, b t, 0 . . . O) r k I k ) t rl--t If = - - e 1 -- bt+let + btet+l = - e l + bt(e~+l - bet).

(5)

x = x' - b x l p , (6) t h e n b y (6), (4), a n d (5), we o b t a i n B x - b = - b t + l T x i ( e t + 1 - b e t ) . T h e u p d a t e p r o c e d u r e in (6) is the third p h a s e in o u r a l g o r i t h m , It is n o t h a r d to verify t h a t the n u m b e r o f f l o a t i n g - p o i n t o p e r a t i o n s r e q u i r e d in the u p d a t e p r o c e d u r e (6) is a b o u t 2t. T a k i n g s u p - n o r m o n b o t h sides, b e c a u s e o f Ibl < 1, it gives

Hnx - bll = [b'+avx;[ < Ib~+l~'llx'll. (7)

T o estimate Irx'll, we need the following l e m m a .

1 L e m m a 1. I f B ' x ' = b', then

IIx'l/~

2 lal I Itb'll. Proof: B e c a u s e o f B ' x ' = b', s u p p o s e ]1 x' II = I x~l for s o m e i, 2 < i _< n - 1, we h a v e a x ; = h; - x ; - 1 - x ; + l .

U s i n g the triangle inequality, it follows t h a t Ildx'll = [dx;I _< Ib[[ + IXs -< tlb'll + 2llx'H l IIx'll -< - - I I b ' l l 9 I d l - 2

T h e cases Ilx'll =

Ixll

a n d

Irx'll = Ix;I

are c o n s i d e r e d next. S u p p o s e

I]x'll

= Ixil, we h a v e

: t a x i = b l - x z .

U s i n g the triangle i n e q u a l i t y again, it follows t h a t I[ax']l = lax'll -< ]bl[ + Ix;I _< Ilb'll + Ilx'll 1 IIx'll - < - IIb'll.

i a ] - 1

Finally, s u p p o s e IIx'll = Ix'.[, we h a v e

t t d x ; = b. - x . - 1 . Similarly, it follows t h a t Ildx'll = Idx',[ -< Ib'.l + Ix',-ll -< Ilb'll + [Ix'll

1

I I x ' l l - < - Nb'll.

I d l - 1

W e c o n c l u d e t h a t l[x'li _< m a x - 2 ' L a l - l ' L d L - 1 Nb'lI. S i n c e a - b - - d a n d Lb[ < 1, we h a v e [dl < lal + Ibl < la[ + 1, a n d h e n c e 0 < Idl - 2 < lal - 1, b e c a u s e of Id[ > 2. H e n c e , we o b t a i n m a x I d l - 2 ' l a l - l ' l d l 1 - [d[ 2" T h e r e f o r e ,

1

we h a v e IIx'll -< [ d l ~ IIb'll. [ ] By (7), L e m m a 1, a n d b = 7b', we h a v e Llnx - b]l < Ib'+ly[llx'll < ib,+~Ti IIb'll I d l - 2 Ibm+a[ Ilbll

I , / I - 2

I m m e d i a t e l y , we h a v e the following t h e o r e m . Jb'+llllbJI T h e o r e m 2. ]]Bx - bll _< I d l - 2

Recall that Ibj < 1, so the b o u n d of IlBx - bll is d e c a y e d e x p o n e n t i a l l y in terms of b. F o r e x a m p l e , if we h o p e Ilnx - b]l < r t h e n

log(IdL - 2) + l o g ~

t(~) _> - 1. (8)

The smallest t(~) of (8) will be small enough for some relative tolerance ~ and sufficiently large diagonal dominance ratio Idl. Using (8), Table 1 illustrates some relations between ]dl's and the smallest t's for some 4. Note that the tolerance cannot be achieved when t(f) > n, and hence the algorithm will break down when the required relative tolerance is sufficiently small and/or the diagonal dominance ratio is sufficiently close to 2.

ldl 2.001 2.01 2.05 2.1 2.5 4 6 8 Table 1 t(10 -2) t(10 -4} t(10 -6) 364 92 34 21 7 2 1 1 509 138 54 36 14 6 4 3

t(10 -8)

655 800 184 230 75 95 51 65 20 27 9 13 7 9 5 8

In summary, based on our three-phase algorithm described in this section, solving the symmetric Toeplitz tridiagonal system Bx = b takes about 5n + 2t. For some relative tolerance and highly diagonal dominance ratio, we have that t is O(1). That is, it takes about 5n to solve Bx = b for this case. The number of floating point operations required in our algorithm is the same as the previous fastest ones such as the special L U factorization and reversed triangular factorization [2 3.

3. Symmetric Circulant Tridiagonal Systems

We now consider the symmetric circulant tridiagonal system Ax = b, where

ii i]

A = Y By (2), d = fl, and a - b = d, we have 7 - T b A = 7B' + 7 (9)

(lo)

1 T o s o l v e A x = b, we first s o l v e B ' x ' = = b ( = b') as d e s c r i b e d in S e c t i o n 2. F r o m (10), we o b t a i n A x ' = 7 B ' x ' + V(x', - bx'~)e~ + ,/x'le, = b + 7(x'~ - b x l ) e ~ + 7 x l e , . (11) F r o m b 2 + db + 1 = 0 ( d e r i v e d in S e c t i o n 2), it f o l l o w s t h a t 1 A p = (db + b z, b + db 2 + b 3 . . . b ' - i + db t, b ~, 0 , . . . , 0, b) r k - M t k ) t nY--l: = (db + bZ, 0 , . . . , 0 , b t - t + dbt, U , O , . . . , O , b ) r k " v I k y " _ _ J t n - t

= - e I - bt+let + btet+l + be,

= ( - e 1 + be,) + b~(et+~ - be,). 0 2 )

W i t h r e s p e c t to the r e v e r s e f o r m of p, n a t u r a l l y , we let q = ( 0 , . . . , O, b t . . . . , b 2, bl) r,

k,.~,~f~) k Y )

r l - - t t

f r o m w h i c h

1

- A q = ( - % + be~) + bt(e,_t - ben-~+l) 7

= (bex - en) + br(e,_t - be._t+~). (13)

1 1

L e t u = (p + bq) a n d v = L ( b p + q). B y (12) a n d (13), we h a v e

7 7

A u = (b 2 - 1)e 1 + bt(et+l - be t + ben_ , - b2e._t+x);

( 1 4 ) A v = (b 2 - 1)e,, + b ' ( b e , + i - b Z e t + e . - t - b e . - t + i ) . L e t ? x = x ' b 2 - 1 [ ( x " - b x l ) a + x l v ] ' b y (15), (11), a n d (14), we o b t a i n _ _ 7b~ , , b ) x l ) e , _ t A x b b 1 [ x " e ' + l - b x . e t + (bx'. + (1 - 2 , (b2x', + (b 3 , - - b ) x l ) e n _ , + l ] . S u p p o s e {t + 1, t} c~ {n - t,n - t + 1} = ~ . Since [b] < 1, we h a v e 7b t IIAx - btl = ~

max(lx'.l, lbx'. +

(1 - b 2 ) x i l ) yb' -< I ' ~ - - i - l ' ( I b l I + l1 - b 2 l ) l l x ' [ [ . (15)

Using the s a m e a r g u m e n t s as L e m m a 1, we also have llx' II -< - - I d l - 2

JIb'll

Id] - 2 fl b' I[. Therefore, 7b ~ [ l A x -

bll ~ ~

(Ibl + l1 - b2[) -

(Ib[

+ I1 - b21)lbtl

= ~ 2 -- ]-I~d~ Z 2) tlbll.

(16)

So, the b o u n d of

[lAx - bll

is decayed exponentially in terms of b. F o r example, if we require ]lAx -

bll

< ~llblJ, then

log([dt - 2) + l o g l b 2 - 11 - log(lb[ + [1 - b2[) + l o g ~

t(() _> - 1. (17)

loglbl

T a b l e 2 illustrates s o m e relations between

Idrs

a n d the smallest t's for s o m e 4. N o t e that the tolerance ~ c a n n o t be achieved when t(~) > n, a n d hence the a l g o r i t h m will b r e a k d o w n when the required relative tolerance is sufficiently small a n d / o r the diagonal d o m i n a n c e ratio is sufficiently close to 2.

Idl Table 2 t(10 -2) t(10 -4) t(10 -6) 2.00l 454 2.0i 111 2.05 40 2.1 25 2.5 9 4 4 6 2 8 2 599 I57 60 40 16 7 5 4 745 203 81 55 22 11 8 6 t(10 -8) 891 249 102 69 29 14 10 9 By (15), we have

7 ((x'~- bxl)u +

x l v ) X = X' b 2 - 1 x ' b 2 - l ((x" - bx~)(p + b q ) +

x'~(bp +

q))

1

= x' b 2 - 1 (x',p +

(bx', +

(1 - bE)x~)q)

= X, b2X"

(

bx'l )

_ l p +

x[

b~_-- 1- q.

It is easy to verify that the n u m b e r of floating-point o p e r a t i o n s required in the a b o v e u p d a t e p r o c e d u r e is a b o u t 4t. In s u m m a r y , based on o u r three-phase algorithm, solving the s y m m e t r i c circulant Toeplitz tridiagonal system Ax = b takes a b o u t 5n + 4t. F o r s o m e relative tolerance a n d highly d i a g o n a l d o m i n a n c e ratio, we h a v e

that t is O(1) too. It follows that the number of floating-point operations needed in our algorithm can complete with the previous fastest ones [2]. Following the derivations of our three-phase algorithm, the special near-Toeplitz tridiagonal

system Cx = b can be solved in a similar way, where

C =

and C is diagonally dominant.

Acknowledgements

The authors are indebted to the referee for several insightful comments. These comments improved the quality and presentation of this paper.

References

[1] Bartlett, M. S,: An inverse matrix adjustment arising in discriminant analysis. Ann. Math. Statist. 22, 107-111 (1951).

[2] Boisvert, R. F.: Algorithms for special tridiagonal systems. SIAM J. Sci. Stat. Comput. 12, 423-442 (1991).

[3] Chen• M. K.: •n the s••uti•n •f •ir•u•ant •inear systems• SIAM J. Numer. Ana•. 24• 668-683 ( •987). [4] Dongarra, J. J., Moler, C. B., Bunch, J. R., Stewart, G. W.: LINPACK user's guide. SIAM Press,

1979.

[5] Evans, D. J., Forrington, C. V. D.: Note on the solution of certain tri-diagonal systems of linear equations. Comput. J. 5, 327-328 (1963).

[6] Evans, D. J.: An algorithm for the solution of certain tri-diagonal systems of linear equations. Comput. J. 15, 356-359 (1972).

[7] Evans, D. J.: On the solution of certain Toeplitz tridiagonal linear systems. SIAM J. Numer. Anal. I7, 675-680 (1980).

[8] Fischer, D., Golub, G., Hald, O., Levia, C., Winlund, O.: On Fourier-ToepIitz methods for separable elliptic problems. Math. Comput. 28, 349-368 (1974).

[9] Hockney, R. W.: A fast direct solution of Poisson~ equation using Fourier analysis. J. ACM. 12, 95-113 (1965).

[10] Malcolm, M. A., Palmer, J.: A fast method for solving a class of tridiagonal linear systems. Comm. ACM. 17, 14-17 (1974).

[11] Pham, B.: Quadratic B-splines for automatic curve and surface fitting. Comput. Graphics 13, 471-475 (1989).

[12] Rojo, O.: A new method for solving symmetric circulant tridiagonal systems of linear equations. Comput. Math. Appl. 20, 61-67 (1990).

[13] Smith, G. D.: Numerical solution of partial differential equations: finite difference methods, 3rd edn. New York: Oxford University Press 1985.

[14] Widlund, O. B.: On the use of fast methods for separable finite difference equations for the solution of general elliptic problems. In: Sparse matrices and their applications (Rose, D. J., Willoughby, R. A., eds.), pp. 121-131. New York: Plenum Press 1972.

Mr. W.-M. Yan

Department of Computer Science and Information Engineering National Taiwan University Taipei, Taiwan 10764, R.O.C. Taiwan

Dr. K.-L. Chung

Department of Information Management National Taiwan Institut of Technology Taipei, Taiwan 10672, R.O.C.