Levinson and Fast Choleski Algorithms for Toeplitz and Almost Toeplitz Matrices RLE Technical Report No. 538 December 1988 Bruce R. Musicus Research Laboratory of Electronics Massachusetts Institute of Technology Cambridge, MA 02139 USA 0· A *IC A- Levinson and Fast Choleski Algorithms for Toeplitz and Almost Toeplitz Matrices Bruce R. lusicus Research Laboratory of Electronics Massachusetts Institute of Technology Cambridge, Mass. 02139 ABSTRACT In this paper, we review Levinson and fast Choleski algorithms for solving sets of linear equations involving Toeplitz or almost Toeplitz matrices. The Levinson-Trench-Zohar algorithm is first presented for solving problems involving exactly Toeplitz matrices. A fast Choleski algorithm is derived by a simple linear transforma- tion. The almost Toeplitz problem is then considered and a Levinson-style algorithm is proposed for solving it. A set of linear transformations converts the algorithm into a fast Choleski method. Symmetric and band diagonal applications are con- sidered. Formulas for the inverse of an almost Toeplitz matrix are derived. The relationship between the fast Choleski algorithms and a Euclidian algorithm is exploited in order to derive accelerated "doubling" algorithms for inverting the matrix. Finally, strategies for removing the strongly nonsingular con- straint of Levinson recursion are considered. We conclude by applying the techniques to several applications, including covari- ance methods of linear prediction, rational Toeplitz matrices, and optimal finite interval ARMA smoothing filters. November 3, 1981 Levinson and Fast Choleski Algorithms for Toeplitz and Almost Toeplitz Matrices Bruce R. Musicus Research Laboratory of Electronics Massachusetts Institute of Technology Cambridge, Mass. 02139 1. Introduction One of the most common problems in numerical calculation is to solve a set of linear equations: RN.N = (1.1) for the (N+1) long vector _N where RN is an (N+1)x(N+1) matrix. Standard methods for solving this problem, such as gaussian elimination or Choleski decomposition, generally require O(N 3) operations. When RN has additional structure, however, the computation can often be significantly reduced. In par- ticular, if RN is Toeplitz with (i,j)tl element Ri j=r(i-j), then the Levinson- Trench-Zohar recursive algorithms1 ,2 ,3' 4 can solve the linear equations using only O(N2) operations and O(N) storage. Similar fast algorithms proposed origi- nally by Friedlander, Morf, Kailath, and Ljung5 and others apply when RN is almost Toeplitz in the sense of having "low displacement rank". These Levinson-style algorithms can be viewed as fast procedures for decomposing the inverse matrix RN1 into a prc luct of Upper triangular, Diagonal, and Lower triangular (UDL) matrices. Applying a simple linear transformation to these algorithms yields a set of "fast Choleski" algorithms which compute an LDU decomposition of RN itself using only O(N2) operations. These algorithms were first discussed by Bareiss6 Morf,7 and Rissanen. 8 In general, the fast Choleski algorithms will either require much more storage or somewhat more computa- -2- tion than the Levinson-style algorithms in order to compute ZN. When the matrix RN is also band diagonal, however, variants of the fast Choleski algo- rithm can be derived which are significantly superior to the Levinson-style band diagonal algorithms suggested by Trench 9 and Dickinson,10 and asymptotically superior to the O(NlogN) matrix splitting and imbedding approaches of Jain11 and Morf, Kailath.12 The fast Choleski algorithms bear a remarkable resemblance to a Euclidian polynomial algorithm. Since a "divide and conquer" strategy combined with Fast Fourier Transforms (FFT's) can be used to accelerate Euclidian algo- rithms,13 it can also be used to accelerate our fast Choleski algorithms. The result is an O(Nlog2N) doubling algorithm for computing RN1, which is similar to those of Gustavson and Yun,14 Bitmead and Anderson,15 and Morf16 Unfor- tunately, the algorithm is relatively complex, so that even for exactly Toeplitz matrices it is only advantageous for matrices of size N>2000. One difficulty with all these algorithms, except that of Gustavson and Yun, is that they require that all the upper left principal minors of RN must be non- singular. In the closing sections of this report, we show how this constraint can be removed. This paper is intended as a coherent summary of Toeplitz matrix algo- rithms, as well as a presentation of several new results. The approach used for deriving the fast Choleski algorithm appears to be new. The displacement rank formalism of Kailath, Kung, Morf17 is used in deriving the Levinson-style almost-Toeplitz algorithms, instead of that of Friedlander et al.18 This simplifies the derivation and the inversion formulas. The band diagonal fast Choleski algorithms using forward and backward recursions to minimize storage appear to be new. Section 8 contains a new partial LDU decomposition formula for RN, which suggests an interesting result for Schur complements. The derivation of -3- the doubling algorithms is completely new, and is considerably more concise and powerful than previous algorithms. The method for dealing with degen- eracy is also new. Finally, in addition to several examples previously considered in the literature, we also present several new applications of the methods, including a very fast finite interval Wiener-Hopf smoothing filter for ARMA models, and a problem in which RN is a rational matrix. This last example has also been considered by Dickinson;19 our approach is much faster, and appears to be numerically stable. 2. The Levinson-Trench-Zohar Algorithm for Exactly Toeplitz Matrices The Levinson-Trench-Zohar (LTZ) algorithm4 is a recursive method for solving the simultaneous linear equations RN.N=/N when the matrix RN is exactly Toeplitz. For simplicity we will consider the case when the entries of R, are scalar, although the case of block Toeplitz matrices can be handled in much the same way.820.21 Let the entries of vectors N and RN be iN and Yi,N for i=0,...,N. Let R be the (n+l)x(n+1) upper left principal minor of RN, and let n n be the vector containing the first (n+1) components of N. Note that because of the Toeplitz structure of RN, we can partition each minor R, so as to show its relationship to the lower order minor Rn_ : 1 Rn,-~ r(-n) r(O) ... r(-(f) gn = : : (2.1) r(n) r(O) (n) R (2.1) This structure suggests an approach for solving the linear equations (1.1) in which we recursively calculate the solution to the following problems for n=O, .. ,.N: Rna = · Rn,, = ; Rn = (2.2) -4- where On is a vector of n zeroes, e will be defined later, and where: n (2.3) = ( 1 1,T bi 1)T bn = ( bnn sn = ( ,n .. n ) The solution for n =0 can be found by inspection: _o =o = (1) ; 0 -r(o) o =I Yo (2.4) The solutions at stage n>O can now be calculated recursively using the solu- tions at stage n-1. To do this, assume that we know the values of x_-l. n -1,b n-l E -1, and Then from (2.1) it is easy to show that: En -1 n1-n -1 Rn~ 0 _ O -1 and Rn l- (2.5) -Cn E -i En -1 where 1 n-I (n = - Z r(n-j)a, n-l (2.6) :n-1 j=0 1 2 V =- - T(-)bnjn n -I j=i Values of a, and b which satisfy (2.2) can thus be computed as appropriate linear combinations of a-1 and b4-: -Z -r TIn1 (2.7) an = + n Direct substitution shows that: (2.8) En = E - ( -n Vn) Finally, from (2.1) we see that: (2.9) n-o where X, =Yn - r(n-j)zjn 2 =o -5- Thus the solution z can be computed as: = + bn (2.10) To summarize, we start with the initial solutions a , bo, x in (2.4), then use the 0 0 recursions in (2.7) and (2.10) to calculate a., b , x for n=O,...N, at which A point we will have found the desired solution EN. Total computation will be about 3N2 operations (1 operation t 1 add + 1 multiply). Because the computa- tion can be done in place, total storage required will be about 5N locations for a, b, L, and r (-N), . . .,r(N). (The solution for z, can be stored in the same location used for y,,.) The algorithm will work correctly provided that En X0 at each stage n. We will see later that this condition is equivalent to requiring that all the principal minors R, of RN must all be non-singular, a condition known as "strong non-singularity". In the terminology of linear prediction, the vectors a and b are known as the forward and backward predictors, e is the n prediction error, and tn and v, are the forward and backward reflection coefficients. Note that if RN is symmetric, with r(n)=r(-n), then the vector a will be identical to b except with the elements in reverse order, a =bj,n. The for- ward and backward reflection coefficients will also be identical, Cn=k,. These relationships can be used to cut the computation required to only 2N2 opera- tions, and cut the storage required to about 3N locations. If the matrix RN is also positive definite, then by Sylvester's criterion all the principal minors R, will be non-singular and thus en will never be zero. An interesting interpretation of the vectors and b can be gained by forming the (n+l)x(n+l) matrices A and B whose jth rows contain the coefficients of the vectors aj and bi. -6- 1 0 1 0 bil 1 B, = (2.11) An = a * b,n " b l,n 1 an,,n " al, 1 From (2.2), and noting that R is a principal minor of RN, it is easy to show that n A R and R BT are upper and lower triangular matrices respectively: n n *.. ,, 4It Eo o0 to a.nd RB = . . (2.12) AnRn = tt ~t e I, This in turn implies that 01 A, = AnRB2 = (2.13) 0 en Rearranging gives: R' = BT-LAn (2.14) The various predictor coefficients generated by the Levinson-Trench-Zohar thus form an Upper triangular, Diagonal, Lower triangular (UDL) decomposition of the inverse matrix R- 1. This interpretation suggests several interesting results. For example, we could calculate the vector xn by exploiting the fact that: (2.15) Wceo ultdh us compute recursively by: We could thus compute xn recursively by: I + AnL (2.16) n where XA = E anj,nYi n J=0 As noted by Zohar, however, this formula for X appears to have no obvious n advantages over the formula in (2.9). The UDL decomposition in (2.14) also implies that: -7- n det(Rn) = 1I ej (2.17) j=0 This proves that the restriction that E O for all n is equivalent to requiring n ~ that all principal minors R of RN must be non-singular. n Several interesting formulas for the inverse matrix R- 1 were suggested by Gohberg and Semencul.22 If R and Rn, are both invertible, then they can be n 1 expressed as sums of products of upper times lower triangular Toeplitz matrices: 1 b ,n "bn ,n 1 0 a ,n R,, = 1 b1 n , (2.18) en 0I 1 ann " al,n 1 0 nn " a l,n 0 0 bb ,n n~ ' ,n 0 b ,,,, " bn,n0 1 n n-l,n ' 1,n' -1 _ 1 * b ,n. :(2.19) 0 1 an-l,n " an 1 an,n "*al,n bn l0 an,n b 1,n bn,n We will derive similar formulas for the more general case of almost Toeplitz matrices in section 8. The important point is that the inverse matrix RP1 can be completely specified by the vectors , bAN and EN Moreover, we can com- pute ZN=R -1N solely from knowledge of aN, bN, EN and do not need to actually compute or store the elements of R 1. In fact, forming the product R-lyN only rvLr~ rrIr~ L F ~1CCLCLI~IC·J VI1\N CC~ VIII~Ll:Y UC ~N LF LL