On the stability of the Bareiss and related Toeplitz factorization algorithms
Adam W. Bojanczyk, Richard P. Brent, Frank R. de Hoog, Douglas R. Sweet
aa r X i v : . [ m a t h . NA ] A p r On the stability of the Bareiss and relatedToeplitz factorization algorithms ∗ A. W. BojanczykSchool of Electrical EngineeringCornell UniversityIthaca, NY 14853-5401, USA R. P. BrentComputer Sciences LaboratoryAustralian National UniversityCanberra, ACT 0200, AustraliaF. R. de HoogDivision of Maths. and Stats.CSIROCanberra, ACT 2601, Australia D. R. SweetElectronics Research LaboratoryDSTOSalisbury, SA 5108, AustraliaReport TR-CS-93-14November 1993
Abstract
This paper contains a numerical stability analysis of factorization algorithms for com-puting the Cholesky decomposition of symmetric positive definite matrices of displacementrank 2. The algorithms in the class can be expressed as sequences of elementary downdating steps. The stability of the factorization algorithms follows directly from the numerical prop-erties of algorithms for realizing elementary downdating operations. It is shown that theBareiss algorithm for factorizing a symmetric positive definite Toeplitz matrix is in the classand hence the Bareiss algorithm is stable. Some numerical experiments that compare behav-ior of the Bareiss algorithm and the Levinson algorithm are presented. These experimentsindicate that in general (when the reflection coefficients are not all positive) the Levinsonalgorithm can give much larger residuals than the Bareiss algorithm.
We consider the numerical stability of algorithms for solving a linear system
T x = b, (1.1)where T is an n × n positive definite Toeplitz matrix and b is an n × ǫ by first computing theCholesky factor of T . Hence the emphasis of the paper is on factorization algorithms for thematrix T .Roundoff error analyses of Toeplitz systems solvers have been given by Cybenko [10] andSweet [22]. Cybenko showed that the Levinson-Durbin algorithm produces a residual which,under the condition that all reflection coefficients are positive, is of comparable size to thatproduced by the well behaved Cholesky method. He hypothesised that the same is true evenif the reflection coefficients are not all positive. If correct, this would indicate that numericalquality of the Levinson-Durbin algorithm is comparable to that of the Cholesky method. ∗ Copyright c (cid:13) A TEX n his PhD thesis [22], Sweet presented a roundoff error analysis of a variant of the Bareissalgorithm [2], and concluded that the algorithm is numerically stable (in the sense specified inSection 7). In this paper we strengthen and generalize these early results on the stability ofthe Bareiss algorithm. In particular, our approach via elementary downdating greatly simplifiesroundoff error analysis and makes it applicable to a larger-than-Toeplitz class of matrices.After introducing the notation and the concept of elementary downdating in Sections 2 and 3,in Section 4 we derive matrix factorization algorithms as a sequence of elementary downdatingoperations (see also [4]). In Section 5 we present a first order analysis by bounding the firstterm in an asymptotic expansion for the error in powers of ǫ . By analyzing the propagation offirst order error in the sequence of downdatings that define the algorithms, we obtain boundson the perturbations of the factors in the decompositions. We show that the computed uppertriangular factor ˜ U of a positive definite Toeplitz matrix T satisfies T = ˜ U T ˜ U + ∆ T , || ∆ T || ≤ c ( n ) ǫ || T || , where c ( n ) is a low order polynomial in n and is independent of the condition number of T .Many of the results of Sections 2–5 were first reported in [5], which also contains some resultson the stability of Levinson’s algorithm.In Section 6 we discuss the connection with the Bareiss algorithm and conclude that theBareiss algorithm is stable for the class of symmetric positive definite matrices. Finally, inSection 7 we report some interesting numerical examples that contrast the behaviour of theBareiss algorithm with that of the Levinson algorithm. We show numerically that, in caseswhere the reflection coefficients are not all positive, the Levinson algorithm can give muchlarger residuals than the Bareiss or Cholesky algorithms. Unless it is clear from the context, all vectors are real and of dimension n . Likewise, all matricesare real and their default dimension is n × n . If a ∈ ℜ n , k a k denotes the usual Euclidean norm,and if T ∈ ℜ n × n , k T k denotes the induced matrix norm: k T k = max k a k =1 k T a k . Our primary interest is in a symmetric positive definite Toeplitz matrix T whose i, j th entryis t ij = t | i − j | . We denote by e k , k = 1 , . . . , n , the unit vector whose k th element is 1 and whose otherelements are 0. We use the following special matrices: Z ≡ n − X i =1 e i +1 e Ti = · · · · · ·
01 0 · · · · · · ,J ≡ n X i =1 e n − i +1 e Ti = · · · · · · · · · · ...0 1 · ...1 0 · · · · · · . Z is known as a shift-down matrix. We also make use of powers of the matrix Z ,for which we introduce the following notation: Z k = ( I if k = 0, Z k if k > J is called a reversal matrix, because the effect of applying J to avector is to reverse the order of components of the vector: J x x ... x n = x n x n − ... x . The hyperbolic rotation matrix H ( θ ) ∈ ℜ × is defined by H ( θ ) = 1cos θ " − sin θ − sin θ . (2.1)The matrix H ( θ ) satisfies the relation H ( θ ) " − H ( θ ) = " − , and it has eigenvalues λ ( θ ), λ ( θ ) given by λ ( θ ) = λ − ( θ ) = sec θ − tan θ. (2.2)For a given pair of real numbers a and b with | a | > | b | , there exists a hyperbolic rotation matrix H ( θ ) such that H ( θ ) " ab = " √ a − b . (2.3)The angle of rotation θ is determined bysin θ = ba , cos θ = √ a − b a . (2.4) In this section we introduce the concept of elementary downdating. The elementary downdatingproblem is a special case of a more general downdating problem that arises in Cholesky factor-ization of a positive definite difference of two outer product matrices [1, 6, 7, 12]. In Section 4,factorization algorithms are derived in terms of a sequence of downdating steps. The numericalproperties of the algorithms are then related to the properties of the sequence of elementarydowndating steps.Let u k , v k ∈ ℜ n have the following form: k ↓ u Tk = [0 . . . × × × . . . × ] , v Tk = [0 . . . × × . . . × ] , ↑ k + 13hat is: e Tj u k = 0 , j < k , and e Tj v k = 0 , j ≤ k . Applying the shift-down matrix Z to u k , we have k + 1 ↓ u Tk Z T = [0 . . . × × . . . × ] , v Tk = [0 . . . × × . . . × ] . ↑ k + 1Suppose that we wish to find u k +1 , v k +1 ∈ ℜ n to satisfy u k +1 u Tk +1 − v k +1 v Tk +1 = Z u k u Tk Z T − v k v Tk , (3.1)where k + 1 ↓ u Tk +1 = [0 . . . × × . . . × ] , v Tk +1 = [0 . . . × . . . × ] , ↑ k + 2that is e Tj u k +1 = 0 , j < k + 1 , and e Tj v k +1 = 0 , j ≤ k + 1 . (3.2)We refer to the problem of finding u k +1 and v k +1 to satisfy (3.1), given u k and v k , as the elementary downdating problem. It can be rewritten as follows:[ u k +1 v k +1 ] " − u Tk +1 v Tk +1 = [ Z u k v k ] " − u Tk Z T v Tk . From (2.1), (2.3) and (2.4), it is clear that the vectors u k +1 and v k +1 can be found by using ahyperbolic rotation H ( θ k ) defined by the following relations:sin θ k = e Tk +1 v k / e Tk u k , (3.3a)cos θ k = q − sin θ k , (3.3b)and " u Tk +1 v Tk +1 = H ( θ k ) " u Tk Z T v Tk . (3.4)The elementary downdating problem has a unique solution (up to obvious sign changes) if | e Tk u k | > | e Tk +1 v k | . The calculation of u k +1 , v k +1 via (3.4) can be performed in the obvious manner. Followingcommon usage, algorithms which perform downdating in this manner will be referred to as hyperbolic downdating algorithms.Some computational advantages may be obtained by rewriting (3.1) as follows:[ u k +1 v k ] " u Tk +1 v Tk = [ Z u k v k +1 ] " u Tk Z T v Tk +1 . G ( θ k ), G ( θ k ) = " cos θ k sin θ k − sin θ k cos θ k , where cos θ k and sin θ k are defined by (3.3b) and (3.3a), respectively. Then it is easy to checkthat G ( θ k ) " u Tk +1 v Tk = " u Tk Z T v Tk +1 , (3.5)or, equivalently, " u Tk +1 v Tk = G ( θ k ) T " u Tk Z T v Tk +1 . (3.6)Thus, we may rewrite (3.6) as v k +1 = ( v k − sin θ k Z u k ) / cos θ k , (3.7a) u k +1 = − sin θ k v k +1 + cos θ k Z u k . (3.7b)Note that equation (3.7a) is the same as the second component of (3.4). However, (3.7b) differsfrom the first component of (3.4) as it uses v k +1 in place of v k to define u k +1 . It is possible toconstruct an alternative algorithm by using the first component of (3.5) to define u k +1 . Thisleads to the following formulas: u k +1 = ( Z u k − sin θ k v k ) / cos θ k , (3.8a) v k +1 = − sin θ k u k +1 + cos θ k v k . (3.8b)We call algorithms based on (3.7a)–(3.7b) or (3.8a)–(3.8b) mixed elementary downdating al-gorithms. The reason for considering mixed algorithms is that they have superior stabilityproperties to hyperbolic algorithms in the following sense.Let ˜ u k , ˜ v k be the values of u k , v k that are computed in floating point arithmetic with relativemachine precision ǫ . The computed values ˜ u k , ˜ v k satisfy a perturbed version of (3.1), that is,˜ u k +1 ˜ u Tk +1 − ˜ v k +1 ˜ v Tk +1 = Z ˜ u k ˜ u Tk Z T − ˜ v k ˜ v Tk + ǫG k + O ( ǫ ) , (3.9)where the second order term O ( ǫ ) should be understood as a matrix whose elements are boundedby a constant multiple of ǫ || G k || . The norm of the perturbation G k depends on the precisespecification of the algorithm used. It can be shown [6] that the term G k satisfies k G k k ≤ c m (cid:16) k Z u k k + k v k k + k u k +1 k + k v k +1 k (cid:17) (3.10)when a mixed downdating strategy is used (here c m is a positive constant). When hyperbolicdowndating is used the term G k satisfies k G k k ≤ c h k H ( θ k ) k ( k Z u k k + k v k k ) ( k u k +1 k + k v k +1 k ) , (3.11)where c h is a positive constant [6]. (The constants c m and c h are dependent on implementationdetails, but are of order unity and independent of n .) Note the presence of the multiplier k H ( θ k ) k in the bound (3.11) but not in (3.10). In view of (2.2), k H ( θ k ) k could be large. The significanceof the multiplier k H ( θ k ) k depends on the context in which the downdating arises. We considerthe implications of the bounds (3.10) and (3.11) in Section 5 after we make a connection betweendowndating and the factorization of Toeplitz matrices.5t is easily seen that a single step of the hyperbolic or mixed downdating algorithm requires4( n − k )+ O (1) multiplications. A substantial increase in efficiency can be achieved by consideringthe following modified downdating problem. Given α k , β k ∈ ℜ and w k , x k ∈ ℜ n that satisfy e Tj w k = 0 , j < k and e Tj x k = 0 , j ≤ k , find α k +1 , β k +1 and w k +1 , x k +1 ∈ ℜ n that satisfy α k +1 w k +1 w Tk +1 − β k +1 x k +1 x Tk +1 = α k Z w k w Tk Z T − β k x k x Tk , with e Tj w k = 0 , j < k and e Tj x k = 0 , j ≤ k . If we make the identification u k = α k w k and v k = β k x k , then we find that the modified elementary downdating problem is equivalent to the elementarydowndating problem. However, the extra parameters can be chosen judiciously to eliminatesome multiplications. For example, if we take α k = β k , α k +1 = β k +1 , then from (3.3a), (3.3b)and (3.4), sin θ k = e Tk +1 x k / e Tk w k , (3.12a) α k +1 = α k / cos θ k , (3.12b)and w k +1 = Z w k − sin θ k x k , (3.13a) x k +1 = − sin θ k Z w k + x k . (3.13b)Equations (3.12a)–(3.13b) form a basis for a scaled hyperbolic elementary downdating algorithmwhich requires 2( n − k ) + O (1) multiplications. This is about half the number required by theunscaled algorithm based on (3.4). (The price is an increased likelihood of underflow or overflow,but this can be avoided if suitable precautions are taken in the code.)Similarly, from (3.7a) and (3.7b) we can obtain a scaled mixed elementary downdating algo-rithm via sin θ k = β k e Tk +1 x k /α k e Tk w k ,α k +1 = α k cos θ k ,β k +1 = β k / cos θ k , and x k +1 = x k − sin θ k α k β k Z w k , w k +1 = − sin θ k β k +1 α k +1 x k +1 + Z w k . The stability properties of scaled mixed algorithms are similar to those of the correspondingunscaled algorithms [12]. 6
Symmetric Factorization
We adopt the following definition from [18].
Definition 4.1: An n × n symmetric matrix T has displacement rank 2 iff there exist vectors u , v ∈ ℜ n such that T − ZT Z T = uu T − vv T . (4.1) ✷ The vectors u and v are called the generators of T and determine the matrix T uniquely.Whenever we want to stress the dependence of T on u and v we write T = T ( u , v ).In the sequel we will be concerned with a subset T of all matrices satisfying (4.1). Thesubset is defined as follows. Definition 4.2:
A matrix T is in T iff (a) T is positive definite, (b) T satisfies (4.1) with generators u and v , (c) v T e = 0, i.e., the first component of v is zero. ✷ It is well known that positive definite n × n Toeplitz matrices form a subset of T . Indeed, if T = ( t | i − j | ) n − i,j =0 , then T − ZT Z T = uu T − vv T , where u T = ( t , t , . . . , t n − ) / √ t , v T = (0 , t , . . . , t n − ) / √ t . The set T also contains matrices which are not Toeplitz, as the following example shows. Example:
Let T =
25 20 1520 32 2915 29 40 , u = and v = . It is easy to check that T is positive definite. Moreover, T − ZT Z T =
25 20 1520 7 915 9 8 =
25 20 1520 16 1215 12 9 − = uu T − vv T . Hence T = T ( u , v ) ∈ T , but T is not Toeplitz. ✷ We now establish a connection between the elementary downdating problem and symmetricfactorizations of a matrix from the set T .Let T = T ( u , v ) ∈ T . Set u = u , v = v and, for k = 1 , . . . , n − , solve the elementary downdating problem defined by (3.1), u k +1 u Tk +1 − v k +1 v Tk +1 = Z u k u Tk Z T − v k v Tk , k . On summing over k = 1 , . . . , n − n − X k =1 u k +1 u Tk +1 − n − X k =1 v k +1 v Tk +1 = n − X k =1 Z u k u Tk Z T − n − X k =1 v k v Tk . If we now observe that, from (3.2), Z u n = v n = 0 , we arrive at the following relation: n X k =1 u k u Tk − Z n X k =1 u k u Tk ! Z T = u u T − v v T , (4.2)which implies that P nk =1 u k u Tk ∈ T . Moreover, as matrices having the same generators areidentical, we obtain T = n X k =1 u k u Tk = U T U , where U = n X k =1 e k u Tk is upper triangular, and hence is the Cholesky factor of T . We have derived, albeit in a ratherindirect manner, the basis of an algorithm for calculating the Cholesky decomposition of a matrixfrom the set T .We now return to the question of existence of a solution to the elementary downdatingproblem for each k = 1 , . . . , n −
1. It is easy to verify that, if T ∈ T , then | e T u | > | e T v | .Using (4.2) and (3.1), it can be shown by induction on k that | e Tk u k | > | e Tk +1 v k | , k = 2 , . . . , n − . Consequently, | sin θ k | < k = 1 , . . . , n − T = T ( u , v ) ∈ T . Algorithm
FACTOR( T ) :Set u = u , v = v .For k = 1 , . . . , n − u k +1 , v k +1 such that u k +1 u Tk +1 − v k +1 v Tk +1 = Z u k u Tk Z T − v k v Tk , e Tk +1 v k +1 = 0 . Then T = U T U, where U = P nk =1 e k u Tk . ✷ In fact we have not one algorithm but a class of factorization algorithms, where each al-gorithm corresponds to a particular way of realizing the elementary downdating steps. Forexample, the connection with the scaled elementary downdating problem is straightforward. Onmaking the identification u k = α k w k and v k = β k x k , (4.3)we obtain T = W T D W , W = n X k =1 e k w Tk ,D = n X k =1 α k e k e Tk . It is clear from Section 3 that Algorithm
FACTOR( T ) requires 2 n + O ( n ) multiplicationswhen the unscaled version of elementary downdating is used, and n + O ( n ) multiplicationswhen the scaled version of elementary downdating is used. However, in the sequel we do notdwell on the precise details of algorithms. Using (4.3), we can relate algorithms based on thescaled elementary downdating problem to those based on the unscaled elementary downdatingproblem. Thus, for simplicity, we consider only the unscaled elementary downdating algorithms. In this section we present a numerical stability analysis of the factorization of T ∈ T viaAlgorithm FACTOR( T ) . The result of the analysis is applied to the case when the matrix T isToeplitz.Let ˜ u k , ˜ v k be the values of u k , v k that are computed in floating point arithmetic with relativemachine relative precision ǫ . The computed quantities ˜ u k and ˜ v k satisfy the relations˜ u k = u k + O ( ǫ ) , ˜ v k = v k + O ( ǫ ) , (5.1)and the aim of this section is to provide a first order analysis of the error. By a first order analysiswe mean that the error can be bounded by a function which has an asymptotic expansion inpowers of ǫ , but we only consider the first term of this asymptotic expansion. One should thinkof ǫ →
0+ while the problem remains fixed [19]. Thus, in this section (except for Corollary 5.1)we omit functions of n from the “ O ” terms in relations such as (5.1) and (5.2).The computed vectors ˜ u k , ˜ v k satisfy a perturbed version (3.9) of (3.1). On summing (3.9)over k = 1 , . . . , n − T − Z ˜ T Z T = ˜ u ˜ u T − ˜ v ˜ v T − ( Z ˜ u n ˜ u Tn Z T − ˜ v n ˜ v Tn ) + ǫ n − X k =1 G k + O ( ǫ ) , where ˜ T = ˜ U T ˜ U , ˜ U = n X k =1 e k ˜ u Tk . Since Z ˜ u n = O ( ǫ ) , ˜ v n = O ( ǫ ) , we find that ˜ T − Z ˜ T Z T = ˜ u ˜ u T − ˜ v ˜ v T + ǫ n − X k =1 G k + O ( ǫ ) . (5.2)Now define ˜ E = ˜ T − T. (5.3)9hen, using (4.1), (5.2) and (5.3),˜ E − Z ˜ EZ T = ˜ u ˜ u T − uu T + ˜ v ˜ v T − vv T + ǫ n − X k =1 G k + O ( ǫ ) . In a similar manner we obtain expressions for Z j ˜ EZ Tj − Z j +1 ˜ EZ Tj +1 , j = 0 , . . . , n −
1. Summingover j gives˜ E = n − X j =0 Z j (cid:16) (˜ u ˜ u T − u u T ) + (˜ v ˜ v T − v v T ) (cid:17) Z Tj + ǫ n − X j =0 n − X k =1 Z j G k Z Tj + O ( ǫ ) . (5.4)We see from (5.4) that the error consists of two parts – the first associated with initial errorsand the second associated with the fact that (5.2) contains an inhomogeneous term. Now k ˜ u ˜ u T − uu T k ≤ k u k k ˜ u − u k + O ( ǫ ) , k ˜ v ˜ v T − vv T k ≤ k v k k ˜ v − v k + O ( ǫ ) . Furthermore, from (4.1),
T r ( T ) − T r ( ZT Z T ) = k u k − k v k > , and hence (cid:13)(cid:13)(cid:13) n − X j =0 Z j (˜ u ˜ u T − uu T + ˜ v ˜ v T − vv T ) Z Tj (cid:13)(cid:13)(cid:13) ≤ n k u k (cid:16) k ˜ u − u k + k ˜ v − v k (cid:17) + O ( ǫ ) . (5.5)This demonstrates that initial errors do not propagate unduly. To investigate the double sumin (5.4) we require a preliminary result. Lemma 5.1
For k = 1 , , . . . , n − j = 0 , , , . . . , k Z j v k k ≤ k Z j +1 u k k . ✷ Proof
Let T k = T − k X l =1 u l u Tl = n X l = k +1 u l u Tl . It is easy to verify that T k − ZT k Z T = Z u k u Tk Z T − v k v Tk and, since T k is positive semi-definite, T r (cid:16) Z j T k Z Tj − Z j +1 T k Z Tj +1 (cid:17) = k Z j +1 u k k − k Z j v k k ≥ . ✷ We now demonstrate stability when the mixed version of elementary downdating is used inAlgorithm
FACTOR( T ) . In this case the inhomogeneous term G k satisfies a shifted versionof (3.10), that is k Z j G k Z Tj k ≤ c m (cid:16) k Z j +1 u k k + k Z j v k k + k Z j u k +1 k + k Z j v k +1 k (cid:17) , (5.6)where c m is a positive constant. 10 heorem 5.1 Assume that (3.9) and (5.6) hold. Then k T − ˜ U T ˜ U k ≤ n k u k (cid:16) k ˜ u − u k + k ˜ v − v k (cid:17) + 4 ǫc m n − X j =0 T r ( Z j T Z Tj ) + O ( ǫ ) . ✷ Proof
Using Lemma 5.1, k Z j G k Z Tj k ≤ c m (cid:16) k Z j +1 u k k + k Z j u k +1 k (cid:17) . Furthermore, since
T r ( Z j T Z Tj ) = n X k =1 k Z j u k k , it follows that (cid:13)(cid:13)(cid:13) n − X j =0 n X k =1 Z j G k Z Tj (cid:13)(cid:13)(cid:13) ≤ c m n − X j =0 T r ( Z j T Z Tj ) . (5.7)The result now follows from (5.4), (5.5) and (5.7). ✷ For the hyperbolic version of the elementary downdating algorithms a shifted version of theweaker bound (3.11) on G k holds (see [6]), namely k Z j G k Z Tj k ≤ c h k H ( θ k ) k ( k Z j +1 u k k + k Z j v k k )( k Z j u k +1 k + k Z j v k +1 k ) . (5.8)By Lemma 5.1, this simplifies to k Z j G k Z Tj k ≤ c h k H ( θ k ) k k Z j +1 u k k k Z j u k +1 k . (5.9)The essential difference between (3.10) and (3.11) is the occurence of the multiplier k H ( θ k ) k which can be quite large. This term explains numerical difficulties in applications such as thedowndating of a Cholesky decomposition [6]. However, because of the special structure of thematrix T , it is of lesser importance here, in view of the following result. Lemma 5.2
For k = 1 , , . . . , n − j = 0 , , . . . , n − k , k H ( θ k ) k k Z j u k +1 k ≤ n − k − j ) k Z j +1 u k k . ✷ Proof
It is easy to verify from (3.4) that1 ∓ sin θ k cos θ k ( u k +1 ∓ v k +1 ) = Z u k ∓ v k , and from (2.1) that k H ( θ k ) k = 1 + | sin θ | cos θ . Thus, k H ( θ k ) k k Z j u k +1 k ≤ k H ( θ k ) k k Z j v k +1 k + k Z j +1 u k k + k Z j v k k≤ k H ( θ k ) k k Z j +1 u k +1 k + 2 k Z j +1 u k k , where the last inequality was obtained using Lemma 5.1. Thus k H ( θ k ) k k Z j u k +1 k ≤ n − k X l = j +1 k Z l u k k , and the result follows. ✷ emark Lemma 5.2 does not hold for the computed quantities unless we introduce an O ( ǫ )term. However in a first order analysis we only need it to hold for the exact quantities. Theorem 5.2
Assume that (3.9) and (5.8) hold. Then k T − ˜ U T ˜ U k ≤ n k u k ( k ˜ u − u k + k ˜ v − v k ) + 8 ǫc h n − X j =1 ( n − j ) T r ( Z j T Z Tj ) + O ( ǫ ) . ✷ Proof
Applying Lemma 5.2 to (5.9) gives k Z j G k Z Tj k ≤ c h ( n − j − k Z j +1 u k k , and hence (cid:13)(cid:13)(cid:13) n − X j =0 n − X k =1 Z j G k Z Tj (cid:13)(cid:13)(cid:13) ≤ c h n − X j =1 n − X k =1 ( n − j ) k Z j u k k ≤ c h n − X j =1 ( n − j ) T r ( Z j T Z Tj ) . (5.10)The result now follows from (5.4), (5.5) and (5.10). ✷ Note that, when T is Toeplitz, T r ( Z j T Z Tj ) = ( n − j ) t . Hence, from Theorems 5.1 and 5.2, we obtain our main result on the stability of the factorizationalgorithms based on Algorithm
FACTOR( T ) for a symmetric positive definite Toeplitz matrix: Corollary 5.1
The factorization algorithm
FACTOR( T ) applied to a symmetric positivedefinite Toeplitz matrix T produces an upper triangular matrix ˜ U such that T = ˜ U T ˜ U + ∆ T , where k ∆ T k = O ( ǫt n ) when mixed downdating is used, and k ∆ T k = O ( ǫt n ) when hyperbolicdowndating is used. ✷ In his 1969 paper [2], Bareiss proposed an O ( n ) algorithm for solving Toeplitz linear systems.For a symmetric Toeplitz matrix T , the algorithm, called a symmetric Bareiss algorithm in [22],can be expressed as follows. Start with a matrix A (0) := T and partition it in two ways: A (0) = U (0) T (1) ! , A (0) = T ( − L (0) ! , where U (0) is the first row of T and L (0) is the last row of T . Now, starting from A (0) , computesuccessively two matrix sequences { A ( i ) } and { A ( − i ) } , i = 1 , . . . , n − , according to the relations A ( i ) = A ( i − − α i − Z i A ( − i +1) , A ( − i ) = A ( − i +1) − α − i +1 Z Ti A ( i − . (6.1)12or 1 ≤ i ≤ n − , partition A ( i ) and A ( − i ) as follows: A ( i ) = U ( i ) T ( i +1) ! , A ( − i ) = T ( − i − L ( i ) ! , where U ( i ) denotes the first i + 1 rows of A ( i ) , and L ( i ) denotes the last i + 1 rows of A ( − i ) . Itis shown in [2] that (a) T ( i +1) and T ( − i − are Toeplitz, (b) for a proper choice of α i − and α − i +1 , the matrices L ( i ) and U ( i ) are lower and uppertrapezoidal, respectively, and (c) with the choice of α i − and α − i +1 as in (b), the Toeplitz matrix T ( − i − has zero elementsin positions 2 , . . . , i + 1 of its first row, while the Toeplitz matrix T ( i +1) has zero elementsin positions n − , . . . , n − i of its last row.Pictorially, A ( i ) = U ( i ) T ( i +1) ! = × · · · · · · · · · · · · · · · · · · × × × ... . . . × × ...0 · · · × · · · · · · · · · ×× · · · × · · · · · · × ... . . . . . . . . . ...... . . . . . . . . . ... × · · · · · · × · · · × A ( − i ) = T ( − i − L ( i ) ! = × · · · × · · · · · · × ... . . . . . . . . . ...... . . . . . . . . . ... × · · · · · · × · · · ×× · · · · · · × × · · · × × · · · · · · · · · · · · · · · × After n − A ( n − and A ( − n +1) are lower and upper triangular, respec-tively. At step i only rows i + 1 , . . . , n of A ( i ) and rows 1 , , . . . , n − i of A ( − i ) are modified; theremaining rows stay unchanged. Moreover, Bareiss [2] noticed that, because of the symmetryof T , T ( i +1) = J i +1 T ( − i − J n and α i − = α − i +1 , (6.2)Here J i +1 and J n are the reversal matrices of dimension ( i + 1) × ( i + 1) and n × n respectively.Now, taking into account (6.2), it can be seen that the essential part of a step of the Bareissalgorithm (6.1) can be written as follows: t ( i +1) i +2 t ( i +1) i +3 . . . t ( i +1) n t ( − i − i +3 . . . t ( − i − n ! = − α i − − α i − ! t ( i ) i +2 t ( i ) i +3 . . . t ( i ) n t ( − i ) i +2 t ( − i ) i +3 . . . t ( − i ) n ! , (6.3)13here ( t ( − i ) i +2 , t ( − i ) i +2 , . . . , t ( − i ) n ) are the last ( n − i −
1) components of the first row of T ( − i ) , and( t ( i ) i +2 , t ( i ) i +3 , . . . , t ( i ) n ) are the last ( n − i −
1) components of the first row of T ( i ) .Note that (6.3) has the same form as (3.13a)–(3.13b), and hence a connection between theBareiss algorithm and algorithm F ACT OR ( T ) is evident. That such a connection exists wasobserved by Sweet [22], and later by Delosme and Ipsen [11]. Sweet [22] related a step of theBareiss algorithm (6.3) to a step of Bennett’s downdating procedure [3]. Next, he derived the LU factorization of a Toeplitz matrix as a sequence of Bennett’s downdating steps. Finally, heestimated the forward error in the decomposition using Fletcher and Powell’s methodology [12].This paper generalizes and presents new derivations of the results obtained in [22]. We adopt from [17] the following definitions of forward and backward stability.
Definition 7.1:
An algorithm for solving the equation (1.1) is forward stable if the computedsolution ˜ x satisfies || x − ˜ x || ≤ c ( n ) ǫ cond( T ) || ˜ x || , where cond( T ) = || T || || T − || is the condition number of T , and c ( n ) may grow at most as fastas a polynomial in n , the dimension of the system. Definition 7.2:
An algorithm for solving the equation (1.1) is backward stable if the com-puted solution ˜ x satisfies || T ˜ x − b || ≤ c ( n ) ǫ || T || || ˜ x || , where c ( n ) may grow at most as fast as a polynomial in n , the dimension of the system.It is known that an algorithm (for solving a system of linear equations) is backward stableiff there exists a matrix ∆ T such that( T + ∆ T )˜ x = b , || ∆ T || ≤ c ( n ) ǫ || T || , where c ( n ) may grow at most as fast as a polynomial in n .Note that our definitions do not require the perturbation ∆ T to be Toeplitz, even if thematrix T is Toeplitz. The case that ∆ T is Toeplitz is discussed in [13, 24]. The reader isreferred to [9, 14, 19] for a detailed treatment of roundoff analysis for general matrix algorithms.It is easy to see that backward stability implies forward stability, but not vice versa. This ismanifested by the size of the residual vector.Cybenko [10] showed that the L norm of the inverse of a n × n symmetric positive definiteToeplitz matrix T n is bounded bymax n Q n − i =1 cos θ i , Q n − i =1 (1 + sin θ i ) o ≤ k T − n k ≤ n − Y i =1 | sin θ i | − | sin θ i | , where {− sin θ i } n − i =1 are quantities called reflection coefficients . It is not difficult to pick thereflection coefficients in such a way that the corresponding Toeplitz matrix T n satisfiescond( T n ) ≈ /ǫ . One possible way of constructing a Toeplitz matrix with given reflection coefficients {− sin θ i } n − i =1 is by tracing the elementary downdating steps backwards.An example of a symmetric positive definite Toeplitz matrix that can be made poorly con-ditioned by suitable choice of parameters is the Prolate matrix [21, 23], defined by t k = ( ω if k = 0, sin(2 πωk ) πk otherwise , ≤ ω ≤ . For small ω the eigenvalues of the Prolate matrix cluster around 0 and 1.We performed numerical tests in which we solved systems of Toeplitz linear equations usingvariants of the Bareiss and Levinson algorithms, and (for comparison) the standard Choleskymethod. The relative machine precision was ǫ = 2 − ≈ − . We varied the dimension of thesystem from 10 to 100, the condition number of the matrix from 1 to ǫ − , the signs of reflectioncoefficients, and the right hand side so the magnitude of the norm of the solution vector variedfrom 1 to ǫ − . In each test we monitored the errors in the decomposition, in the solution vector,and the size of the residual vector.Let x B and x L denote the solutions computed by the Bareiss and Levinson algorithms. Also,let r B = T x B − b and r L = T x L − b . Then for the Bareiss algorithms we always observed thatthe scaled residual s B ≡ k r B k ǫ k x B kk T k was of order unity, as small as would be expected for a backward stable method. However, wewere not able to find an example which would demonstrate the superiority of the Bareiss algo-rithm based on mixed downdating over the Bareiss algorithm based on hyperbolic downdating.In fact, the Bareiss algorithm based on hyperbolic downdating often gave slightly smaller errorsthan the Bareiss algorithm based on mixed downdating. In our experiments with Bareiss algo-rithms, neither the norm of the error matrix in the decomposition of T nor the residual error inthe solution seemed to depend in any clear way on n , although a quadratic or cubic dependencewould be expected from the worst-case error bounds of Theorems 5.1–5.2 and Corollary 5.1.For well conditioned systems the Bareiss and Levinson algorithms behaved similarly, andgave results comparable to results produced by a general stable method (the Cholesky method).Differences between the Bareiss and Levinson algorithms were noticeable only for very ill-conditioned systems and special right-hand side vectors.For the Levinson algorithm, when the matrix was very ill-conditioned and the norm of thesolution vector was of order unity (that is, when the norm of the solution vector did not reflectthe ill-conditioning of the matrix), we often observed that the scaled residual s L ≡ k r L k ǫ k x L kk T k , was as large as 10 . Varah [23] was the first to observe this behavior of the Levinson algorithm onthe Prolate matrix. Higham and Pickering [16] used a search method proposed in [15] to generateToeplitz matrices for which the Levinson algorithm yields large residual errors. However, thesearch never produced s L larger than 5 · . It plausible that s L is a slowly increasing functionof n and 1 /ǫ .Tables 7.1–7.3 show typical behavior of the Cholesky, Bareiss and Levinson algorithms forill-conditioned Toeplitz systems of linear equations when the norm of the solution vectors is oforder unity. The decomposition error was measured for the Cholesky and Bareiss algorithms bythe quotient || T − L · L T || / ( ǫ · || T || ), where L was the computed factor of T . The solution errorwas measured by the quotient || x comp − x || / || x || , where x comp was the computed solution vector.Finally, the residual error was measured by the quotient || T · x comp − b || / ( || T || · || x comp || · ǫ ).decomp. error soln. error resid. errorCholesky 5 . · − . · − . · Bareiss(hyp) 3 . · . · − . · − Bareiss(mixed) 2 . · . · . · Levinson 5 . · . · n = 21, ω = 0 . cond = 3 . · decomp. error soln. error resid. errorCholesky 1 . · − . · − . · − Bareiss(hyp) 2 . · . · − . · − Bareiss(mixed) 3 . · . · − . · − Levinson 5 . · − . · Table 7.2: Reflection coefficients | sin θ i | of the same magnitude | K | butalternating signs, | K | = 0 . n = 41, cond = 8 . · decomp. error soln. error resid. errorCholesky 8 . · − . · − . · − Bareiss(hyp) 8 . · . · − . · − Bareiss(mixed) 6 . · . · − . · − Levinson 2 . · − . · Table 7.3: Reflection coefficients | sin θ i | of the same magnitude butalternating signs, | K | = 0 . n = 92, cond = 2 . · This paper generalizes and presents new derivations of results obtained earlier by Sweet [22].The bound in Corollary 5.1 for the case of mixed downdating is stronger than that given in [22].The applicability of the Bareiss algorithms based on elementary downdating steps is extended toa class of matrices, satisfying Definition 4.2, which includes symmetric positive definite Toeplitzmatrices. The approach via elementary downdating greatly simplifies roundoff error analysis.Lemmas 5.1 and 5.2 appear to be new. The stability of the Bareiss algorithms follows directlyfrom these Lemmas and the results on the roundoff error analysis for elementary downdatingsteps given in [6].The approach via downdating can be extended to the symmetric factorization of positivedefinite matrices of displacement rank k ≥ k the factorization algo-rithm uses elementary rank- k downdating via hyperbolic Householder or mixed Householderreflections [8, 20].We conclude by noting that the Bariess algorithms guarantee small residual errors in thesense of Definition 7.2, but the Levinson algorithm can yield residuals at least five orders ofmagnitude larger than those expected for a backward stable method. This result suggests that,if the Levinson algorithm is used in applications where the reflection coefficients are not knownin advance to be positive, the residuals should be computed to see if they are acceptably small.This can be done in O ( n log n ) arithmetic operations (using the FFT).It is an interesting open question whether the Levinson algorithm can give scaled residualerrors which are arbitrarily large (for matrices which are numerically nonsingular). A relatedquestion is whether the Levinson algorithm, for positive definite Toeplitz matrices T without arestriction on the reflection coefficients, is stable in the sense of Definitions 7.1 or 7.2.16 eferences [1] S.T. Alexander, C-T. Pan and R.J. Plemmons, “Analysis of a Recursive least SquaresHyperbolic Rotation Algorithm for Signal Processing”, Linear Algebra and Its Applications ,vol 98, pp 3-40, 1988.[2] E.H. Bareiss, “Numerical Solution of Linear Equations with Toeplitz and Vector ToeplitzMatrices”,
Numerische Mathematik , vol 13, pp 404-424, 1969.[3] J.M. Bennett, “Triangular factorization of modified matrices”,
Numerische Mathematik, vol 7, pp 217-221, 1965.[4] A.W. Bojanczyk, R.P. Brent and F.R. de Hoog, “QR Factorization of Toeplitz Matrices”,
Numerische Mathematik , vol 49, pp 81-94, 1986.[5] A.W. Bojanczyk, R.P. Brent and F.R. de Hoog, “Stability Analysis of Fast Toeplitz LinearSystem Solvers”, Report CMA-MR17-91, Centre for Mathematical Analysis, The AustralianNational University, August 1991.[6] A.W. Bojanczyk, R.P. Brent, P. Van Dooren and F.R. de Hoog, “A Note on Downdatingthe Cholesky Factorization”,
SIAM J. Sci. Statist. Comput. , vol 8, pp 210-220, 1987.[7] A.W. Bojanczyk and A. Steinhardt, “Matrix Downdating Techniques for Signal Process-ing”,
Proceedings of the SPIE Conference on Advanced Algorithms and Architectures forSignal Processing III , vol 975, pp 68-75, 1988.[8] A.W. Bojanczyk and A.O. Steinhardt, “Stabilized Hyperbolic Householder Transforma-tions”,
IEEE Trans. Acoustics, Speech and Signal Processing , vol ASSP-37, 1989, pp 1286-1288.[9] J.R. Bunch, “The Weak and Strong Stability of Algorithms in Numerical Linear Algebra”,
Linear Algebra and Its Applications , vol 88/89, pp 49-66, 1987.[10] G. Cybenko, “The Numerical Stability of the Levinson-Durbin Algorithm for Toeplitz Sys-tems of Equations”,
SIAM J. Sci. Statist. Comput. , vol 1, pp 303-319, 1980.[11] J-M. Delosme and I.C.F. Ipsen, “From Bareiss’s algorithm to the stable computation ofpartial correlations”,
Journal of Computational and Applied Mathematics, vol 27, pp 53-91,1989.[12] R. Fletcher and M.J.D. Powell, “On the Modification of
LDL T Factorizations”,
Mathemat-ics of Computation , vol 28, pp 1067-87, 1974.[13] I. Gohberg, I. Koltracht and D. Xiao, “On the solution of the Yule-Walker equation”,
Proceedings of the SPIE Conference on Advanced Algorithms and Architectures for SignalProcessing IV , vol 1566, July 1991.[14] G.H. Golub and C. Van Loan,
Matrix Computations , second edition, Johns Hopkins Press,Baltimore, Maryland, 1989.[15] N.J. Higham, “Optimization by Direct Search in Matrix Computations”, Numerical Analy-sis Report No. 197, University of Manchester, England, 1991; to appear in
SIAM J. MatrixAnal. Appl. [16] N.J. Higham and R.L. Pickering, private communication.1717] M. Jankowski and H. Wozniakowski, “Iterative Refinement Implies Numerical Stability”,
BIT , vol 17, pp 303-311, 1977.[18] T. Kailath, S.Y. Kung and M. Morf, “Displacement Ranks of Matrices and Linear Equa-tions”,
J. Math. Anal. Appl. , vol 68, pp 395-407, 1979.[19] W. Miller and C. Wrathall,
Software for Roundoff Analysis of Matrix Algorithms , AcademicPress, 1980.[20] C.M. Rader and A.O. Steinhardt, “Hyperbolic Householder Transformations”,
IEEE Trans-action on Acoustics, Speech and Signal Processing, vol ASSP-34, 1986, pp 1584-1602.[21] D. Slepian, “Prolate spheroidal wave functions, Fourier analysis, and uncertainty V: thediscrete case”,
Bell Systems Tech. J. , vol 57, 1978, pp 1371-1430.[22] D. Sweet, “Numerical Methods for Toeplitz Matrices”, PhD thesis, University of Adelaide,1982.[23] J.M. Varah, “The Prolate Matrix: A Good Toeplitz Test Example”,
SIAM Conference onControl, Signals and Systems , San Francisco, 1990. Also