A GENERALIZED FOR THE OF A STRUCTURED

T. B o a o s (0,

SCHUR-TYPE ALGORITHM JOINT FACTORIZATION MATRIX AND ITS INVERSE

A . H . SAYED(2),

H. LEV-ARI (3),

T.

K A I L A T H (4)

- A Schur-type algorithm is presented for the simultaneous triangular factorization of a given (non-degenerate) structured matrix and its inverse. The algorithm takes the displacement generator of a Hermitian, strongly regular matrix / / as an input, and computes the displacement generator of the inverse matrix R -1 as an output. From these generators we can directly deduce the L D - 1 L * (lower-diagonal-upper) decomposition of R, and the U D - 1 U * (upper-diagonallower) decomposition of R -1. The computational complexity of the algorithm is O ( r n 2) operations, where n and r denote the size and the displacement rank of R, respectively. Moreover, this method is especially suited for parallel (systolic array) implementations: using n processors the algorithm can be carried out in (2(n) steps.

ABSTRACT

1.

Introduction

T h e present paper is concerned with the factorization of Hermitian, strongly regular (s) matrices R E C ~• that obey a displacement equation of the form

R-FRF*

(1.1)

=GJG*,

where F 6 C '~• is a general lower triangular array, G 6 C T M generator matrix,

J=Ip69-Iq

(1.2)

p > 0,

is a full-rank

q> 0

is a signature matrix with p positive and q negative entries on the main diagonal (p + q = r << n), and the superscript '*' denotes the Hermitian transpose of a

-

-

(1) (2) (3) (4) (g)

Received: november 16, 1996 Stanford University University of California Los Angeles Northeastern University Stanford University A matrix is strongly regular if all of its leading minors are non-zero.

132

BOROS,SAYED, LEV-ARI, KAILATH: A Generalized Schur-Type Algorithm

matrix. It is assumed that the diagonal entries/o,/1, ..-, fn of F satisfy (1.3)

1--fir; #0,

i, j E { O , 1 , . . . , n - - 1 } ,

so that the Lyapunov equation (1.1) has a unique solution for R in terms of {Y, G, J}. The triplet {F, G, J } is known as a displacement generator of R (see

[151, [91, [17], [5], [20], [21], [161, [1]). By now, it is widely known that for the problem of factoring a strongly regular matrix R that satisfies the displacement equation (1.1) the generalized Schur algorithm gives a nice parallelizable solution (see, eg, [19], [20], [22] and [16]). On the other hand, similar results for computing the triangular factors of the inverse matrix R -I have not yet been obtained for the general case. For Toeplitz matrices, the first and best known algorithm for factoring the inverse is the celebrated Levinson algorithm. Early attempts at extending this algorithm beyond the Toeplitz case were made in IS], [6], [7], [17], but the formulas were rather complicated. After the (re)discovery of the Schur algorithm for directly factoring Toeplitz matrices rather than their inverses, it was realized by several authors (see [14] and [11]) that the Levinson algorithm could be replaced by a two-step procedure: use the Schur algorithm to compute the so-called reflection coefficients, and then use a simple recursion to compute the triangular factors of the inverse. This two-step procedure requires slightly more computation than the Levinson algorithm on a serial machine, but it is significantly less expensive on a parallel machine. The reason is that the classical Levinson algorithm obtains the reflection coefficients via certain inner products, which cannot be parallelized in an efficient manner. This hybrid algorithm was partially exiended in [4] and [13] to matrices obeying displacement equations of the form R - F R F * = GJG*, where F has some strictly lower triangular structure, such as F = Zk (9 Zn_~. However, these extended algorithms generally require a large displacement generator. In the present paper all these restrictions are removed. However, for simplicity, we shall assume that the following nondegeneracy condition holds: (1.4)

i~ j

==~

fi r fj ,

i, j E {O, 1 , . . . , n - 1 }

.

This condition of course fails, eg, when F is strictly lower triangular; this case will be discussed in [3]. Thus, the main problem that we solve here can be described as follows: PROBLEM 1.1. [PROBLEM STATEMENT] Let R be a Hermitian strongly regular matrix so that and are satisfied Given the displacement gen erator {F, G, J } o / R derive a computationally ej~cient parallelizable algorithm that simultaneously computes an Z D - 1 Z * decomposition for R and a U D - 1 U * decomposition for R -1 The following lemma is well known, see eg, [12]; we repeat it here because

for structured matrices

133

its proof will guide us to the desired joint factorization algorithm (see Section 4). LEMMA 1.1. If R is an invertible matrix that satisfies (1.1), then there must exist a full-rank matrix H E @r• so that (1.5)

PROOF.

R -1 - F* R -1 F = H* J H .

The block matrix

decompositions: F"

1--

[RF 1

F.R-1

F*

I

admits the following block-triangular

R -~

0

R-I-F'R-IF

F * R -1

I

Now Sylvester's law of inertia [10, p.416] implies that (1.7)

Inertia{R -~ - F * R - 1 F } = I n e r t i a { R - F R F ' } .

It follows from (1.1) ghat there must exist a full-rank matrix H E @r• (1.5) is valid.

SO that 9

The objective of the main algorithm presented in Section 4 is to map the displacement generator ( F , G, J } of R to the displacement generator { F ' , H*, J } of R -1. This objective is essentially achieved in two steps. First, the displacement generator {F, G, J } is reduced to a set of ( D -1 ~ J ) - u n i t a r y rotations 12~ by using the generalized Schur algorithm. Inversion in this new domain can be performed simply by changing signs and taking hyperbolic conjugates, see formula (5.3) below. Subsequently, the inverted rotations are mapped back to a displacement generator {F*, H*, J } . 2.

Notation

The triangular decompositions of R and R -I will be sought in the somewhat non-traditional forms (2.1)

R = L D -1 L* ,

g/-1 _ U D -1 U* ,

where D----diag{ do, dl, . . . , d~-i }. The reason for requiring D -1 in both factorizations in (2.1) will appear later (see (5.2a-b) in Section 5). The strong regularity assumption on R means that the diagonal entries of D will all be nonzero. The lower triangular factor L is normalized in such a way that do, dl, ..., d,_l appear on its main diagonal, while the diagonal entries of U are normalized to be 1. It is then easy to see that L and U must satisfy the identity (2.2)

L D -1 U* -- I n .

134

BOP~OS,SAYED, L~.v-ARI, KAILATI{: A GeneralizedSchur-Type Algorithm

The nonqzero part of the consecutive columns of L and U will be denoted by li and ui:

(2.3)

L=

..

In-,'

],

(2.4) 9"

',

0

n--i--1

Let P i E C i• be the leading principal submatrix of R, and let R / E C (=-i)x(n-i) be the Schur complement of P i in R: (6)

(2.5)

R = [Pi Qi Q~] s~ '

l~=Si_Qip71Q:,

i=0,1,...,n.

The well-known block inversion formula (2.6)

R-l=

P~ Qi Qi Si

p.~l + pr~Q~R;~Q~p;1 _p;1Q~R.-(1 -R.~QiPT 1 R['

then shows that, after inversion, R~-1 will be the trailing principal submatrix of R -1. Moreover, the Schur complement of R~-1 in R -1 turns out to be p~-1 :

A~ = (p:l + p-(~Q;R~IQp:I) _ (p.-(~O;R.~l) pq (R-(1Q, p-(1)--_ p:(~. An important property is that the matrices P i and tl4 along with their inverses all inherit the displacement structure of R. LEMMA 2.1 (PRINCIPAL BLOCKS AND THEIR SCHUR COMPLEMENTS) Let F i E C ixi and F i E (~(~-i)• denote the leading and trailing principal submatrices of F , respectively9 If (1.1) holds then Pi, Ri -1, Ri and p~-i must satisfy the displacement equations (2.7a)

Pi - Fi Pi P; = Zi J X; ,

(2.75)

_P~ - F ~ R ~ F ;

(2.7c)

p : l _ ~,*ip-(l~,i = H;

(2.7d)

=

G~JG;, J Hi ,

R~-t - _F* R~-~_F~ -- Y; J Y , ,

for some matrices Xi, H* E @i•

and Gi, Y~ e G(~-i)•

(i = 0, 1 , . . . , n).

PROOF. The results follow easily from the fact that/:/-1 can be written in the form (2.6). Indeed, let R, F and G be partitioned as

(6) The boundary conditions are: P0 = Rn = [ ] (void matrix), Pn ---- Ro = R. 9The strong regularity assumption of R implies that Pi and R4 axe non-singular for all i.

135

for structured matrices

where the ' . ' - s denote "don't-care" entries. Substituting (2.8) into (1.1) and evaluating the leading principal blocks on both sides of the equation yields (2.7a), which shows that P i has displacement structure. Now Lemma 1.1 states that p~-I must inherit the displacement structure of Pi- Therefore there must exist a matrix H i so that (2.7c) is satisfied. In a similar manner, let R -1 and H be partitioned as 9

R[ t

,

F=

9

_F i

,

H=

*

Yi

9

Substituting (2.9) into (1.5) and evaluating the trailing principal blocks on both sides of the equations yields (2.7d), which shows that the trailing principal submatrix R~-1 of R -I has displacement structure. Now L e m m a 1.i implies that P~ must inherit the displacement structure of R~-t. Therefore there must exist a matrix Gi so that (2.7b) is satisfied. 9 Note that X i and Y i can be readily obtained from G and H as shown in (2.8) and (2.9). The main algorithm in Section 4 will describe how to construct suitable matrices Gi and H i so as to satisfy (2.7b) and (2.7c). 3.

A General Factorization Procedure

The inverse of any strongly regular, Hermitian matrix can be computed recursively by deflating R and inflating R -1 as shown below in A l g o r i t h m 3.1. The following scheme provides an L D - 1 L * decomposition for R as well as a U D - 1 U * decomposition for R -1 (see [1, p.ll]): ALGORITHM 3.1 (GAUSS/SCHUR RECURSIONS) Start the procedure with Ro = R and A o = [ ], and repeat Steps 1 - 4 below for i = 0 , 1 , . . . , n - 1; S t e p 1: Let li = "first column of Ri" and di ="upper-left-corner element of Ri". S t e p 2: Perform the Gauss/Schur reduction step: (3.1a)

[0o] 0 Ri+l

=R/

-

li di

li 9

S t e p 3: Given {/0, 11,..., l/} and {do, d l , . . . , di}, compute ui by solving a triangular system of equations (see equation (2.2) above). S t e p 4: Perform the Gauss/Schur construction step: (3.1b)

Ai+l__

[

Ai 0

0 0

]

+uld i u i .

At the end of the procedure, A ~ = R -1.

136

BOI~OS,SAYED, LEV-ARI, KAILATH: A Generalized Schur-Type Algorithm

PROOF. Writing out the factorization of R as R = L D - 1 L * = v.,~-I 9 z-.,i----0 ~,t-1~* ~z~i ~ shows that the Schur complements P~ and P~+I are related via (3.1a). Similarly, R - I = U D - 1 U . = ~ =.-1 ~ - i ~ui. implies that the Schur complements A i and 0 u~di z~i+l are related via

(3.2)

[ A~ 0 ] ~- /~/-~l -- '~ di-1 "~~ 0

0

"

Equation (3.1b) is obtained by rearranging the terms in (3.2).

9

Observe that Step 2 is a reduction step ie, a rank-1 matrix is subtracted from P~. In contrast, Step 4 is a construction step in which a rank-1 matrix is added to Ai (after bordering the matrix with zeros). Thus the procedure recursively eliminates the rows and columns of R and successively constructs the rows and columns of R -1. In general, Algorithm 3.i requires O(n 3) operations. However, when R has displacement structure, the computational burden can be significantly reduced by exploiting the fact that the Schur complements P~ and Ai inherit the displacement structure of R as indicated by Lemma 2.1 above. Indeed, the solvability condition (1.3) ensures that the displacement generators (_Fi, Gi, J } and {$'7, H i , J } uniquely determine the Schur complements P~ and p~-l. Therefore, instead of explicitly computing the Schur complement sequences R0 -+ R1 -+ ... -+ P~ and A o --+ A 1 --+ ... --+ A , , it will be sufficient to compute the generator sequences Go --+ G1 -+ ... --+ G,,, and H~ -+ H E --+ ... --> H~. 4.

The Fast Joint Factorization Algorithm

The development of the fast algorithm that we present in this paper was motivated by the state-space approach of [20]; the factorization procedure for a strictly positive definite matrix R was explicitly described there, while the factorization of R -1 was only alluded to in the concluding remarks. Here we present a fast factorization procedure for general (non-definite) matrices that simultaneously computes the triangular factors of R and R -I (see also [1, p.13]). ALGORITHM 4.1 (MAIN ALGORITHM) Start the recursion with G o = G , H o = [ ], _ F o = F , ~'o=[ ], and perform the following steps for i = 0 , . . . , n - 1: S t e p 1 Let gi be the first row of Gi. Compute li by solving the linear system

(4.1)

(

z ) l, =

In particular, obtain the first element of li as d~ = g i J g * / ( 1 - fi f~) 9 Step 2

Given f~ and g~ find matrices h i E C r•

and k i E r r•

so as to satisfy

for structured matrices

137

S t e p 3 Compute G~+I from Gi as

G~+I

= _-Vil~ h~ + G~ J k~.

S t e p 4 Let ~'i+1 and ui be partitioned as qo~ fi

'

ui =

1

"

Compute ~2i by solving the linear system (v) (4.5)

( F ; - Z Ii ) ~2i-- H ; J g; - ~ ; .

S t e p 5 Compute H~+ 1 from H~ as (4.6)

H*+l=uih*+[H~lJk* 0

The next section contains the detailed derivation of the joint factorization algorithm. We shall use a square-root-array-type argument similar to the one used in [21], [23], and [22]. We remark that several other derivation techniques have been used in the past to obtain generalized Schur-type algorithms for the factorization of R, including a generating function (transform-domain) approach [18, 19] and a state-space cascade factorization approach [20]; however, they do not appear to be so directly applicable to the joint factorization problem. 5.

Algorithm Development

The approach used in [23] to derive the Generalized Schur Algorithm for the factorization of R alone was based on combining the expressions R = FR~'* qG J G * and R = L D - 1 L * to obtain the equality

[o,o] o

1

[o o] o

1"

from which it follows that there must exist a (D -t ~ J)-unitary matrix f~ such that

For the joint factorization of R and R -t we follow a similar route but now starting with the block matrix

(5.1)

[.Fa. n_lFj

(T} Note that ~20---[ ].

138

BOROS, SAYED, LEV-ARI, KAILATtt: A Generalized S c h u r - T y p e A l g o r i t h m

t h a t appeared in the proof of L e m m a 1.1. First use the equations R = L D - I L *, R -I - F*R-1F=H*JH, and R - F R F * = G J G * , R-I=UD-1U* t o obtain the block-triangular factorizations

=

F*U

H*

0

J

F*U

H*

T h e center m a t r i x in both (5.2a) and (5.2b) is the same diagonal m a t r i x ( D - I ~ J ) ; it was to achieve this that we started with the non-traditional factorizations R=LD-1L * and R - I = U D - 1 U *. The next step is to observe t h a t the above relations imply (see A p p e n d i x A in [2]) that there must exist a ( D - 1 5 J ) - u n i t a r y matrix f~ such t h a t U

0

F*U

H*

"

Formula (5.3) suggests that if we can find any matrix ~ t h a t lower-triangularizes the left-hand array in (5.3), then we can read out L, H * , and F * U from the result. At first, it is difficult to see how (5.3) can be used to compute L, D , U and H when F and G are given, since unknown quantities appear on b o t h sides of the equation. This apparent difficulty is resolved by proceeding recursively. Essentially ~2 is a ( D -1 ~ J ) - u n i t a r y transformation, t h a t maps a secondary block upper triangular (s) p r e a r r a y to a block lower triangular p o s t a r r a y by eliminating all entries in the upper right corner of the matrix. In order to obtain a recursive algorithm, one can implement ~ as a sequence of n elementary transformations denoted by ~20, f21, ..., Y~-I, so that each transformation eliminates exactly one row of G. This procedure is illustrated below (see [2] for the details). A n e x a m p l e : Let n = 3 following zero-patterns:

.o oOli. **[0i 9

I

U

and r = 2 .

0

*

* * 0 * O0

In this case the pre- and post-arrays have the

*

* * *00

.o Oloooo **.100 *

'

o

F*U

i,o**oooo).** [;:ooooo

H"

=

*

0

* * 0 * O0

* *

* *

* *

*

*

*

[i~ ~176

The (1,2) block entry of the pre-array can be eliminated recursively as follows:

9 9

* *

0 *

0 0

0 0

0 0

* *

* *

0 *

* 0

* 0

o oooo

* 0

* *

* *

0 *

* *

* *

0 0

(8) Upper triangular with respect to the secondary diagonal.

* *

* ] * * *

0 *

*

"

for structured matrices

139

The light gray boxes mark the entries to be eliminated at each step of the recursion by using elementary (D -1 $J)-unitary transformations (scaled column permutations, Givens rotations, and Householder projections). The dark gray boxes mark the position of the pivot elements. The ultimate result of the recursion is that the (1,2) block of the pre-array is eliminated row-by-row ("reduction procedure"), while the (2,2) block is filled up with non-zero elements ( "construction procedure"). 9 The results of the recursive procedure can be guessed from the relation (5.2). Let L, F , D and U be partitioned as

:], .=rv.._..o] o=[ o, o], where Li, $'i, Di, Ui E C ~• ingly as

_<,':':]

The pre- and post-arrays can be partitioned accord-

I~,-<,,+_~,~ ~,-~

o.1_ s.,,

o,

~,

o

_ ~.

' [~%~J-F;o o,

After the i-th step, the first i columns of the post-array are already computed, while the last (i + 1)-th ..., n-th columns of the pre-array remain unchanged. Therefore, the intermediate array obtained after the i-th step has the following form:

(5.4)

[ F 1L 2 G0 ~] f~~ - I U"

Li

0

Li

_FiDi

0 G~

~)'i

H~

0

U~

o

where we have deliberately used Gi and H* to denote the non-zero elements that appear in the (1, 2) and (2, 2) blocks. The reason is that, as is easily shown, Gi is indeed a displacement generator for the Schur complement Ra, while H ~ is a displacement generator for the Schur-complement Ai = p~-i (see Theorem 1.5.2 in [1]). In general the (D -1 @ J ) - u n i t a r y rotations 12i can be implemented as

(5.5)

I2i =

I~ 0 0 0

0 ai 0 c~

0 0

0 b~

In-i-1 0

0 si

where ai E (U, b i- ~C "• , c~EC l• (5.6)

[ ][ ai bi c~ si

d~-1 0 0 J

'

and si E (7"• must satisfy

][ ][ ai bi eq si

=

d~-1 0 0 J

]

"

140

BOROS, SAYED, LEV-ARI, KAILATH: A Generalized Schur-Type Algorithm

Thus the i-th step of the elimination procedure can be succinctly described as _-Fill

(5.7)

[

Gi

04 bi

H;

ca s,

o

=

li

Gi+l

F;+lui

H;+I

-

where the irrelevant columns and rows of the pre- and post-arrays (ie, those rows and columns t h a t remain unchanged) were omitted. Equation (5.7) can be further decomposed into its block components, viz, (5.8a)

_Fi I~ 04 + Gi ca =

(5.8b)

_Fili bi -~ Gi si

0

si

=

=

li,

[o] Gi+l

'

H~+ t 9

In order to proceed, one must first choose suitable rotation parameters 04, bi, ca and si so as to satisfy (5.4) and (5.6). The first equation (5.8a) leads to (I,~_i - F i 04)li = G i ca. However, it can be seen from the displacement equation P~ - FiR~F~ = G i J G ~ that the first column li of P~ also obeys (5.9)

(I,~-i - _Fi f~) li = Gi J g ; .

Thus a possible choice for 04 and ca is given by (a) (5.10)

04 = Z ,

ca = J g ; .

Let us now introduce the new parameters (5.11)

hi=b;,

ki=s*J .

Equation (5.8b) along with (5.10) leads to Gi+t

= -Filih~ + G i J k

i 9

It follows from (5.6) that s g~, hi and k~ must satisfy (4.2). If ~'i+1 and ui are partitioned as shown in (4.4) then equation (5.8c) can be written in the form

Oj[lj=[ 0 JZg~" (9) If i __
.for structured matrices

141

Therefore, (~'~ - ]i" Ii)~2~ = H~- J g~ - ~p~. This equation has a solution for ,2i whenever fi • {f0, f l , . . . , s } (ie, whenever the non-degeneracy condition (1.4) is valid). Finally, equation (5.8d) along with (4.4) and (5.11) lead to (5.i3)

H:+l=u~h*+[H~]jk *.

This concludes the derivation of Algorithm 4.i. Equations (5.12) and (5.13) above are called the generator recursions, since they describe how to update the generator matrices Gi and H r. Specifically, (5.12) is referred to as the forward generator recursion (see [17], [20], [21], [22], [1]), while (5.13) will be referred to as the backward generator recursion (see [20], [1]). The forward generator recursion speeds up the Gauss/Schur reduction procedure by exploiting the displacement structure of R, while the backward generator recursion speeds up the Gauss/Schur construction procedure by exploiting the displacement structure of R -I. 6.

P r o p e r F o r m of t h e G e n e r a t o r R e c u r s i o n s

It turns out (see Theorem. 6.1 below) that the solution {h*, k~} of the (di~J)unitary matrix completion problem (4.2) is not unique. Thus Algorithm 4.1 gives several different specialized algorithms depending on the choice of the parameters hi and ki. The objective of this section is to describe an important so-called proper form of the generator recursions. Let us begin by characterizing the solution set of the (di ~ J)-unitary matrix completion problem. THEOREM 6.1 (SOLUTION OF THE (di 0 J)-UNITARY M C P ) The pair {hi, ki} satisfies (4.2) if, and only if (6.1 a)

hi

(6.1b)

k{ --- 0 i [ J - 1 - / i+E , / ; g:d:(tgi ] ,

=

" *

O i g i d i -1, 9

for some J-unitary parameter Oi that satisfies OiJO~ = J. PROOF. Theorem 6.1 was first published in [20] for the special case when R is positive definite and IAi(~')] < 1 for all i (see Theorem 1 in [20]). A straightforward generalization extends it to non-definite strongly regular matrices, see [21],

[16]. By using Theorem 6.1, the forward and backward generator recursions (4.3) and (4.6) can be reduced to the following form:

gi g" +

]

142

BOP~OS, SAYED,LEV-ARI, KAILATtt: A Generalized Schur-Type Algorithm

where (6.3a)

"hi

=

(I=_i - .f~*_Fi)-:(_Fi

(6.3b)

"I'i

=

(~;

(6.3c)

r

=

- ( ~ ; - f~9 Ii) _, qoi, ,

- f,* Z,)-'(X,

-

- k I,,-i) ,

k~;),

and d~ = giJg*/(1 - A f~). The proper form of the generator recursions can be reached by imposing further restrictions on the J-unitary parameter Oi. Namely, Oi must be chosen so as to transform the vector gi into proper form (see [17], [4], [21], [16]). This means that O~ will eliminate all elements of g~ with the exception of a single pivot element. In particular, whenever giJg~ > 0 holds, there must exist a J-unitary rotation Oi so that (1~

(6.4)

= [

o

o ],

=

With this particular choice the quantities h~ and k~ in (6.1) become

9 " "

'

I~-i

0

'

where ai = 5i/d~, and the generator recursions (6.2a) and (6.2b) can be further written as (6.6a)

[0

(6.6b)

H~+I-- 0 e ' o

Gi+l

[0o o]

+ `hi Gi Oi

[: o]

I,._: +

o OJLO j

o

/r-1

0 0

'

+

'

o~ o].

Equation (6.6a) admits the following simple interpretation: S t e p 1: Transform Gi into proper form by using a J-unitary rotation Oi. S t e p 2: Multiply the first column by `hi and keep all other columns unchanged. Note that the first step eliminates all entries in the first row of G~ except the pivot entry in the upper left corner position (see Figure 1 below). The second step is needed to eliminate the pivot entry (see Figure I below). On the other hand, equation (6.6b) yields the following slightly more complicated backward generator recursion (see [1, p.25]):

Step 1: Multiply H* by O~. Step 2: Multiply the first column by ff2i and keep all other columns unchanged. (lo) If g~Jg~ < 0 was valid for some i, then there would exist a J-unitary matrix Oi so that g i O i = [ 0 ... 0 5i ], 5i = ~ . Thus the subsequent derivation could be repeated with minor modifications.

143

for structured matrices

st~p>x

,00 * 9 ,

00

0

* ]= I ~

st,ep>2

G i ~-~

:

Gi+l

Figure 1: Illustrating the proper form of the forward generator recursion Step 3: Attach a zero row to the bottom of the array. Step 4: Add the correction term ~7i[ r

1 ]* to the first column.

Note that initially H~ is in proper form. Multiplying the array by O~ will destroy this properness. After attaching a zero row to the bottom of the matrix and adding a correction term to the first column, the resulting matrix H,*+I will emerge in proper form again (see Figure 2 below).

Steps ~-2

Step)3

step)4

= H;+ I

Figure 2: Illustrating the proper form of the backward generator recursion

7. Concluding R e m a r k s A fast algorithm was developed for the triangular factorization of the inverse of a given non-degenerate matrix R. The algorithm maps the displacement generator of R to a displacement generator of R -I while simultaneously computing an L D - X L * decomposition for R and a U D - X U * decomposition for R -1. The main algorithm was obtained by combining displacement structure arguments with the Gauss/Schur reduction and construction procedures for a given matrix and its inverse. Specifically, the key formula (5.3) proves that the Gauss/Schur reduction of the displacement generator {F, G, J } and the Gauss-Schur construction of the displacement generator { F ' , H*, J } can be performed simultaneously by using a sequence of (D -1 @ J)-unitary rotations. For brevity and simplicity, the main results were derived for Hermitian, strongly regular matrices; with a little effort, the presented ideas can also be extended to the non-Hermitian case. Moreover, to keep the discussion simple, we imposed the non-degeneracy condition (1.4) on the displacement generator of R. This requirement can be eliminated at the expense of a slight increase in computational complexity, but a greater increase in conceptual complexity, as will be shown in [3]).

144

BOROS, SAYED, LEV-ARI, KAILATtt: A Generalized Sehur-Type Algorithm

Acknowledgment A major part of this paper was prepared while the first and third authors were visiting ~he Electrical and Communication Engineering Department of the Indian Institute of Science, Bangalore. The hospitality of Prof. V. U. Reddy and of the Jawaharlal Nehru Center is gratefully acknowledged.

REFERENCES [1] T. BOLOS, Studies in Displacement Structure Theory. PhD thesis, Stanford University, Stanford, CA, June 1996. [2] T. BOROS, A. H. SAYED, T. KAILATH, AND H. LEv-Aru, A generalized Schurtype algorithm for the joint factorization of a structured matrix and its inverse, Part I: Non-degenerate case, In preparation. [3] T. BOROS, A. H. SAYED, T. KAILATlt, AND H. LEV-ARI, A generalized Schurtype algorithm for the joint factorization of a structured matrix and its inverse, Part II: Degenerate case, In preparation. [4] J. CI-IUN, Fast Array Algorithms for Structured Matrices. PhD thesis, Stanford University, Stanford, CA, 1989. [5] J. CHUN AND T. KAILATH,Displacement structure for Hankel, Vanderrnonde and related (derived) matrices. Lin. Alg. and Its Appl. 151, (1991) 199-227. [6] J. M. DELOSME, Algorithms and Implementations for Linear Least-Squares Estimation. PhD thesis, Stanford University, Stanford, CA, 1982. [7] P. DELSAKTE, Y. GENIN, AND Y. KAMP, A polynomial approach to the generalized Levinson algorithm. IEEE Trans. Inform. Theory 29(2), (March 1983) 268-278. [8] B. FRIEDLANDER, T. KAILATI{, M. MORF, AND L. LJUNG, Extended Levinson and Chandrasekhar equations for general discrete-time linear estimation problems. IEEE Trans. Automatic Control 23(4), (August 1978) 653-659. [9] Y. GENIN, P. VAN DOOREN, T . KAILATH, J. DELOSME, AND M. MORE, On E-lossless transfer functions and related questions. Lin. Alg. and Its Appl. 50, (1983) 251-275. [10] G. H. GOLUB AND C. F. VAN LOAN, Matrix Computation. The Johns Hopkins University Press, 2-nd edition, 1989. [11] T. KAILATH, Signal processing in the VLSI era. In S. Y. Kung, H. Whitehouse, and T. Kailath, editors, VLSI and Modern Sign. Proc., pp. 5-24. Prentice Hall, 1985. [12] T. KAILATH, Signal processing applications of some moment problems. In Proc. of Symp. in Appl. Math. vol. 37, pp. 71-109. American Math. Soc., 1987. [13] T. KAILATH AND J. CHUN, Generalized displacement structure for block-Toeplitz, Toeplitz-bIock, and Toeplitz-derived matrices. SIAM J. Matrix Anal. Appl. 15(1), (January 1994) 114-128.

for structured matrices

145

[14] S. Y. KUNG AND Y. H. Hu, A highly concurrent algorithm and pipelined architecture for solving Toeplitz systems. IEEE Trans. on Acoust., Speech, and Sign. Proc. 31, (1983) 66-76. [15] T. KAILATH, S. Y. KUNG, AND M. MORF, Displacement ranks of a matrix. Bulletin of the American Mathematical Society 1(5), (September 1979) 769773. [16] T. KAILATIt AND A. H. SAYED, Displacement structure: Theory and applications. SIAM Review 37(3), (1995) 297-386. [17] H. LEv-AaI, Nonstationary Lattice-Filter Modeling. PhD thesis, Stanford University, Stanford, CA, December 1983. [18] H. LEv-ARI AND T. KAILATH, Lattice filter parametrization and modeling of nonstationary processes. IEEE Trans. Inform. Theory 30, (1984) 2-16. [19] H. LEv-Aa! AND T. KAILATH, Triangular fuctorization of structured Hermitian matrices. Operator Theory: Advances and Applications 18, (1986) 301-324. [20] H. LEv-ARI AND T. KAILA'rH, State-space approach to factorization of lossless transfer functions and structured matrices. Lin. Alg. and Its Appl. 162-164, (1992) 273-295. [21] A. H. SAYED, Displacement Structure in Signal Processing and Mathematics. PhD thesis, Stanford University, Stanford, CA, August 1992. [22] A. H. SAYED, T. CONSTANTINESCU, AND T. KAILATH, Square-root algorithms for structured matrices, interpolation, and completion problems. In Lin. Alg. for Sign. Proc., vol. 69 of IMA Volumes in Math. and its Appl., pp. 153184. Springer, 1995. [23] A. H. SAYED, H. LEv-ARI, AND T. KAILATH, Time-varying displacement structure and triangular arrays. IEEE Trans. Signal Proc. 42, (May 1994) 10521062.