IEEE T R A N S A C T I O N S O N I N F O R M A T I O N T H E O R Y . VOL. 37, N O . I , J A N U A R Y 1991

43

A Class of Least-Squares Filtering and Identification Algorithms with Systolic Array Architectures
Seth Z. Kalson, Member, IEEE, and Kung Yao,
Member, IEEE
Abstract -A unified approach is presented for deriving a large class of new and previously known time and order recursive least-squares algorithms with systolic array architectures, suitable for high throughput rate and VLSl implementations of space-time filtering and system identification problems. The geometrical derivation given here is unique in that no assumption is made concerning the rank of the sample data correlation matrix. Our method utilizes and extends the concept of oblique projections, as used previously in the derivations of the leastsquares lattice algorithms. Both the growing and sliding memory, exponentially weighted least-squares criteria are considered.
Index Terms-Least-squares systolic arrays.

tions of the least-squares estimation problem: 1) the filtering problem is to find the filtered output y , , ( t ) , where n .

Y,!(t)S Cgl'(t)xi(t), i=l 1ItIT;

(1.2)

2) the identification problem is to find the filter weights g ; ( t ) , i = 1;. ., n, for any t I. T

This generalization of the least-squares estimation problem is important whenever practical space-time or multichannel filtering arises, such as in adaptive antenna arrays, I. INTRODUCTION decision feedback and fractionally spaced channel equalizINIMUM mean-square estimation is an old and ma- ers, etc. In the previous formalization of the least-squares ture subject, pervading throughout much of the com- problem, we do not need formal stochastic characterizations munication and signal processing literature [l]. Specifically, of the sequences but deal only with observed deterministic various versions of the classical discrete time Wiener filtering data. Furthermore, we have restricted T to be finite so as to [2] problem have found widespread applications to such avoid dealing with infinite dimensional Hilbert spaces. This topics as channel equalization, adaptive antenna arrays, noise is not really restrictive from a practical engineering point of canceling, spectrum estimation, and system identification, to view since T can be thought of as some arbitrarily large but name just a few. In the classical discrete time Wiener filter- finite integer. ing problem, given two random discrete time sequences The filtering and identification problems are intimately denoted by y ( t ) and x ( t ) , we desire the optimum finite connected. This can be made precise by constructing an impulse linear filter G that minimizes E { l y ( t > - G * x ( t ) 1 2 ) , isomorphism between suitably chosen spaces so that time where E denotes the statistical ensemble averaging operator recursive algorithms can be easily derived for the filter and G . x ( t ) denotes time-domain filtering. In this paper, we weights g f ( t ) [3]-[5]. However, such an approach will not be will be interested in generalizing this problem where y ( t ) taken here, but rather, we shall primarily be concerned with and x , ( t ) , i = 1,. . . ,n are deterministic complex-valued dis- filtering algorithms since these possess a highly concurrent crete time sequences (or realizations of random sequences) architecture. Once such filtering algorithms are obtained, it with time indices t = 1,2;. * , T . For each time index, let will be shown that the filter weights may be found with g ; ( t ) , i = 1 . .,ti denote a set of n time varying filter weights additional processing, thus allowing one to use the same ; that satisfies the exponentially weighted, growing memory systolic array for both filtering and identification. least-squares criterion As old as the least-squares problem is, new algorithms continue to be developed. This is in part due to the pressing P'- Y ( 7 ) gj'(t)Xr(T) needs of modern communication, radar, and sonar systems in 7=1 which real-time signal processing is important. As a result, the signal processing community has been interested in least-squares algorithms that are either fast [31, [41, [61-[131, in the sense that the total number of numerical computations where p is a real-valued weight factor satisfying 0 < p I 1. performed each time update is in some sense minimized, or For convenience, we shall make the following two classifica- that can take advantage of concurrent processing, such as systolic arrays [13]-[21]. Henceforth, algorithms for systolic arrays will be referred Manuscript received September 16, 1986; revised July 17, 1989. This work was supported in part by the NASA/AMES Contract No. NAG-2- to as systolic algorithms. The differences between fast and 304 and the NSF grant NCR-8814407. systolic algorithms are not always clearly defined. For examS . Kalson is with MIT Lincoln Laboratory, ANI, P.O. Box 73, Lexing- ple, the least-squares lattice algorithm (LSLA) of [31 and [91 ton, MA 02173. K. Yao is with the Electrical Engineering Department, 56-125 B, for the single channel case meets the definition of systolic arrays. Also, the fast, least-squares filtering algorithms of Engr. IV, University of California, Los Angeles, CA 90024-1594. [ l l ] and [16] can be implemented as a cascade of systolic IEEE Log Number 9038860.

filtering, least-squares identification,

M

I

c

(1

0018-9448/91/0100-0043$01.00 0 1991 IEEE

~

44

l t k E T K A N S A C T I O N S O N I N F O R M A T I O N T H E O R Y , VOL. 37, NO. I . J A N U A R Y I V Y 1

arrays, and the so-called Circular Lattice algorithm of [13] nearly meets the definition of being systolic. Many of these fast algorithms solve a special case of the least squares problem of (1.1) in which some or all of the x,(t), i = 1;. ., n are related by time-shifts. However, in this paper we are only concerned with filtering and identification algorithms that solve the general least-squares problem of (1.1 ) and that meet the requirements of systolic arrays exactly. We shall use in our derivation the powerful geometrical concept of oblique projection, which has proven so useful for deriving fast least squares algorithms, specifically, the LeastSquares Lattice Algorithm [3], [9]. However, unlike previous geometrical derivations, we have generalized the concept of oblique projections so that the n x n sample-data correlation matrix R ( t ) of the time sequences x,(t), i = I : . . , n , with elements R,,(t) given by
I

1,2;. ., T, a corresponding U E C ' , with tth component denoted as [ u ] ~ and equal to u ( t ) , t = 1,2; . ., T . By convention, we define [ U ] , 0, t < 1, sometimes referred to as the = pre-windowed assumption. Similarly, with the time sequences y(t) and x,(t), i = 1;. ., n , we associate y E C T and x, E c', i = 1 . . , T I , respectively. : Let ( , ),, t I, be a bilinear functional mapping C' X C7 T into the complex field, defined for any U E CT, U E C T as t ( U , U ) ,=
P

pI--T[*]_[Lq:,
8

(2.1)

-I

where p is a real number satisfying 0 < p I The above 1. equation can be written in the time recursive form

Rij(t)= r= I

p L I P r x , ( ~ ) x , ( ~ ) * , j = l ; * . , n , (1.3) i

which will be used later. Define the following functional, denoted as (1 I(,, mapping C7- into the nonnegative reals by
IIUII,

need not be full rank. Here, the asterisk denotes complex conjugation. Depending upon the physical significance of the time sequences x , ( t ) , the sample-data autocorrelation function R j j ( t )may represent the ij element of a sample-data spatial correlation matrix of spatially separated sequences x , ( t ) and x , ( t > , or any combination of space-time correlation. For example, if for some time sequence x ( t ) we have x , ( t ) = x ( t - i + 11, i = 1,2; . ., n, then R j j ( t )is the sample-data time autocorrelation function associated with x ( t ) . It should be clear how any combination of space-time correlation can be realized. We have made no implicit assumption concerning how the sequences x , ( t ) are related to one another, so that the least-squares problem considered here is very general. As expected, the geometrical approach provides considerable intuition and insight into these algorithms, which are somewhat complex due to the interplay of a large number of variables. More importantly, the geometrical derivation leads rigorously and systematically to a set of vector order and time updates, from which a large class of new and previously known least squares systolic arrays can easily be derived in a unified manner. We have conceptually grouped the least-squares systolic algorithms of this paper into a single class since their structure is a direct consequence of the way in which the order updates are constructed. This is the only class considered in this paper. It can be shown that nearly all published leastsquares systolic arrays solving the general least-squares problem of (1.1) belong to this class. That is to say, they can be obtained directly from the order and time updates presented here, or with at most some algebraic manipulations. There does, however, exist another (perhaps lesser known) set of order updates that lead to systolic algorithms solving the general least-squares problem, [ 131, and consequently these algorithms belong to a different class.
11. NOTATION N D PROBLEM A STATEMENT

=

Jm,

(2.3)

T for t I. It is easily seen that 11 IIT satisfies all the properties of a norm except that ~ ( u I I , = 0 does not imply U = 0, where 0 denotes the null vector. That is, ((U V I ( , = 0 does not imply U = U . Consequently, 11 11, is a pseudonorm [22], and we shall refer to ( , ), as a pseudo inner product. Specifically, 11 11, and ( , ), are a growing memory, exponentially weighted pseudonorm and pseudo inner product, respectively. The exponentially weighted, sliding memory case is discussed in the Appendix. To continue with the preliminaries, we define the nested set of subspaces S ' c C T , i = O , l ; . . , n , as

s"= 0,
S i = linspan { x , , x2;

(2.4) for i = 1;

. ., x i }

. . ,n . (2.5)

We have thus defined a family of pseudo inner product spaces { S i , ( , ) t } , i = 0,l; . ., n , in which the pseudo inner product is indexed by t . The least-squares problem can now be posed as follows: Given y E C'r, x, E C', i = 1; . ., 17, consider at time t the set q( t )

= ( U : U E S " , Ily - uII,

=

q E S"

inf lly - 411,).

(2.6)

Let U E C ' be any member of cp(t). Then for each time t = 1,2; . ., T , the filtering problem finds the filtered output y,,(t)= [U],, while the identification problem finds a set of filterweights g,?(t), i = l ; . . , n , where C:=,g,!(t)x,=u. This is equivalent to the least-squares estimation problem of Section I.
111. OBLIQUE PROJECTIONS

Now, consider the definition of an oblique projection. Definition I : For every U E C T and subspace S c CT, the oblique projection of U onto S with respect to ( , ),, denoted as l ' $ , , is the set of complex T-vectors ~, l',\,/U = { U : U E S,(U

-U,X)/ =

0

vX

E S).

It is most convenient to cast the least-squares problem in a geometric setting. To do this, some preliminary definitions and concepts are first introduced. We shall use the somewhat unique geometric formalism proposed in [3]. To this end, let C T denote the space of all complex valued T-vectors. Associate with any complex-valued time sequence u(t), t =

The definition of P\.,u is formally equivalent to the definition of the conventional orthogonal projection of U onto S. To make this statement precise, associate with any U E C' the t-vector U , E C', where [ u , ] = [ u ] ~ .= 1,2; . ., t . Let S, ~ T be the linear span of all [-vectors U , where U E S. With a slight abuse of notation, let ( , ), also denote an inner

~

KALSON AND YAO: CLASS OF LEAST-SQUARES FILTERING A N D ALGORITHMS

45

product on S ,

X

S , defined i? the obvious way so that

( u , , ~ , ) = ( u , ~ ) , Finally, let P,,,,: C‘ 4 C‘ be the orthogo, .

nal projection operator onto S , with respect to ( , ),. Then is a complete inner product space, it is clear since { S , , ( , from Definition 1 and the well-known properties of orthogonal projections that we have th: following three properties. Property 1: Y E Ps, U v, = Ps,,U ,. Property 2: {U: U E S , 11 - V I ( / = infqEsJ1u- qll,} = P,,,,U. 11 Property 3: U , w E P.,./u * IIu - wll, = 0. Property 1 relates conventional orthogonal projections with oblique projections. Property 2 gives the least-squares interpretation of oblique projections. Property 3 simply states that the first t components of all T-vectors belonging to P,.,u are unique, which easily follows from Property I and the uniqueness of orthogonal projections. These properties bring out the salient features of oblique projections. Since P.y,,v is a set of vectors, we shall define, for any U , U E C’, the pseudo inner product ( U ,Ps,tv)rto be ( u , w ) , , where w is any member of Ps,,v. Note that by Property 3, this definition is independent of which w in Ps,,v is chosen. We give similar meaning to the pseudo inner products (P.y,ru,u>r and ( P s , , ~ , P s , y ) ,That is, ( P s . , ~ , v > t( z , v > t . = and ( P . y , r Ps,,v), = (.z,w),, ~, where z and w are any members of Ps,,u and Ps.,v, respectively. Then by Property 1, and the self-adjoint and idempotent property of orthogonal projections, we have for any u , u E C T ,

-

ered in [9] (for the case of a full rank R ( t ) ) is essentially the first t components of ,/ l l ~ ; ~ , ~ / x lIIP;’l,‘x,llr + 0, ll~, any finite scalar, l l P , ~ , , , x ; I=/ 0. l as a self-adjoint and idempotent operator. Equation (3.1) embodies the self-adjoint and idempotent “character” of oblique projections, and will be used freely throughout much of the derivations in this paper. For any U E C T , the oblique projection P.y,,u is a set of T-vectors, but we shall sometimes refer to any member of P s , p as the “oblique projection of v onto S with respect to ( , ),.” Remark I : In spite of the formal similarities between oblique and orthogonal projections, we remark that oblique projections contain more information than the latter. Specifically, in developing time recursions, we seek a solution to a least-squares problem at time t in terms of the least-squares solution at time t - 1. This solution requires quantities of the form [Ps,,-lul, for appropriate S and U. This information is cot readily available from the orthogonal projection P s , r - I ~ , -As a result, the use of oblique projections simpli,. fies the derivation of time updates and leads to a more compact notation. To avoid using subscripts with superscripts, we shall let P/.,U denote the oblique projection of U onto SI. Note that Po,,u = (O}, the set containing the null vector. Thus, from (2.6) and Property 2 we have

Proofi Let U E P;-,,,u+P j l , , , x j . It is clear that U € K S‘. Consider that (u-Y,x;)/=(UP;-,,,U-~;~l,rX,,Xj)r
=

( P ; ? l , t u , x ; > t- K ( P ; ? I , / x ; , x j ) /

- ( u , P ; ? I , , ~ j )- K ( ~ ; , P , ? l , , x ; ) t . t

Clearly for j = l ; . . , i - l , P j I I , , x j = Oand thus ( u - v V , x j ) , = 0. For the case IIP,? ,,,xiJJt 0, it is easily seen that ( U # U , x i ) , = 0 for K = ( U , P j i , , r ~ i ) t / J I P j ~ I , , xWhereas illf. ~~Pjil,,xj~~,=Oimplies[Pji,,,rj],=Ofor 7=1,2;..,t,which in turn implies (U- U , x i ) ‘ = 0, regardless of the choice for U K . Thus, by Definition 1, v E P;,,U. From Theorem 1, the following corollary follows immediately. Corollary I: For every U E C’, we have for i = 1 . ., n , ;

P‘?,,/U Kp,i,,,x;c Pi+, with K given as in Theorem 1. Our goal in this section is to construct a vector order update algorithm that will yield a particular vector in q ( t )= P,,,y. This algorithm will also require a number of auxiliary vectors, to be defined shortly. The “order update” part of the systolic algorithms derived in this paper will then be obtained by considering various components of these vectors. Only after applying the time updates of the next section to various time varying pseudo inner products occurring in this vector algorithm will order and time recursive systolic algorithms be derived for least-squares filtering. The particular vector in P,,,,y given by our vector algorithm depends on how the scalar K is chosen when the pseudonorm of P j ? , . , x i is zero. We choose K = 0 whenever this occurs. The vector P,,,,y so determined will be denoted

d t ) = p,,.lY and the filtered output ,y,,(t)= [U],,U E P,l,ly, is unique by Property 3. Remark 2: For any U E P,l.,y,we have an expansion U = E ~ l ~ l g , ’ Y f )where the g : ’ ( t ) satisfy the well-known set of xr, normal equations. However, for a singular R(r), this expansion is not unique, and P,,.,y may contain more than one vector. This was the need for defining an oblique projection as a set of vectors. For a full rank R(t), our oblique projection is identical to the oblique projection defined in [3]. Furthermore, the oblique projection of a vector U as consid-

40

IEEE T R A N S A C T I O N S O N I N F O R M A T I O N 'THEORY. VOL.. 37, N O . I.J A N U A K Y 1991

as y , , ( t ) . However, it is more convenient to construct algorithms that give the residual T-vectors r , ( t ) , i = 1,2; . . , n , where

y

From the definition of P,,:y we obviously have y , , ( t ) = - r , J I ) . Again, r , ( t ) is a particular vector in P,,:y determined by our choice of K . Order updates for r , ( t ) are given by direct application of Corollary 1. To express this, and other vector updates concisely, we introduce the following auxiliary vectors and scalars: We follow the convention that an ordered equation with indexes j = i + 1 , . . . ,n is not performed for i = n. The vector algorithm has a redundant set of updates since its purpose is to provide a method by which a large class of algorithms can be derived. For example, if (4.7) is used to find k , , ( t ) ,then (4.6) and (4.9) are not required. We remark that the vector algorithm, excluding (4.61, (4.91, and (4.101, is essentially a modified Gram-Schmidt (MGS) orthogonalization procedure. This observation follows from the fact that b , , ( t ) , i = l ; . . , n is an orthogonal basis (with respect to ( , ),I for the space S",and the upper triangular form of the evaluations of the indexes in (4.7) and (4.8) is identical to that in the MGS algorithm [24]. The desired oblique projection r , ( t ) is given by the simple recursion relation of (4.10). Other orthogonalization procedures can be applied to x,, i = 1,. . ., n. For example [13] represents a different procedure for orthogonalization, leading to a different class of systolic algorithms. Another orthogonalization technique, similar to the construction of forward and backward innovations processes, is given in [23], and leads to lattice type algorithms under various assumptions for the x,, i = 1; . . , n.

In particular, b , , ( t ) can be interpreted as the "innovation" part of x, relative to Sf-' to time t . We also introduce a up function 6: C X C + C, where for mnemonic reasons we take the liberty to express its dependence upon C X C with the notation 6 { x / y } , where both x and y belong to C.

With these definitions, and using Corollary 1 with U replaced with y, we have the vector order update for r , ( t )

Note that r & t ) = P & y = y. In this equation, we need b,,(t),i = I; . ., n , which can be found from order updates for b , , ( t ) in the index i. Applying Corollary 1 with U replaced with xi,we obtain

V. TIME UPDATES
Detailed consideration shows that the geometric theorems found in [3] and [9] are essentially identical, except for some subtle differences. Since our vectors are of fixed dimension T , our derivation of the time updates will be in the spirit of [3]. However, the theorem given here is a generalization of that given in [3] and [9] since no assumption concerning the rank of the sample-data correlation matrix is made. The extension of this theorem to the sliding memory pseudo inner product is given in the Appendix. The time update theorem requires the use of an auxiliary oblique projection. To this end, let d t ) denote the unit T-vector given by

Note that b , , ( t )= x,. We shall give the time updates for k , , ( t ) in the next section. However, k , , ( t ) can also be obtained from order updates via an auxiliary variable. Furthermore, these order updates are necessary for deriving normalized algorithms, see [5]. These auxiliary variables are

From (4.21, we have k , , ( t )= e,,(t). Order updates for e , , ( t ) in the index i follow readily by forming the pseudo inner product of both sides of (4.3) with x, to yield

The auxiliary oblique projection needed is simply the oblique projection of d t ) onto Sf with respect to ( , ),. Since any vector belonging to the oblique projection will do, let c , ( t ) be such a vector satisfying c ; ( r >E P , , , p ( t ) .

By collecting the order updates, along with initialization, the previous results can be summarized in the following algorithm. Vector Algorithm: Let y, x,, i = 1; . ., n be given. Set r O ( t ) = Y, b , , ( t ) = x,, e , , ( [ ) = llblf(t)ll~, i = 1;. . , n . Then r , ( f ) , i = 1 , . . ., n are given by performing the following or-

The basic property of c , ( t ) that makes it so useful for developing time updates is that its pseudo inner product with any vector U is

( e,( t ) , U>,=(

U),=

( T(f ) >

.,41=

[~,.,UI,*

9

(5.1) where (3.1) has been used. In this equation, [P,,,u], be is to

~

KALSON AND YAO: CLASS OF LEAST-SQUARES FILTERING AND ALGORITHMS

47

interpreted as the tth component of any T-vector in l , . (5.3) = From its definition, c , ( t ) is real and nonnegative. FurtherCorollary 2 provides all the time updates needed in derivmore, by applying the Schwartz inequality to d t ) and c ; ( t ) , ing a large number of least-squares algorithms. However, it easily is shown that (5.7) requires knowledge of c,(t), which can be found by i = 0,l;. - , n . 0 I ci(t) I 1, using Theorem 1 with U replaced with d t ) to obtain Remark 3: Since c,,(t) is an oblique projection onto S”, there exists the expansion c,,(t)= Cy=,Gj(t)*xj.If we let G,(t), i = 1,2; . ., n be the n components of the n dimenBy considering the tth component of this update equation, be the n sional column vector G ( t ) and x j ( f ) , i = 1,2;..,n we have the order update for c , ( t ) , components of the n dimensional column vector X ( t ) , then it easily follows from the definition of c , ( t ) that G ( t ) = R ( t ) - ’ X ( t )(provided R ( t ) is full rank), the so-called Kalman gain vector. Furthermore, it follows that c , , ( t ) = X ( t ) + R ( t ) -‘ X ( t ) ,the log-likelihood variable.

If c , ( t ) = 1, then (5.9) shows that [P,,+ulr 0 for any U , in = which case (5.7) is equivalent to either (5.5) or (5.6). If c , ( t ) + 1 , then using (5.9) to solve for [P,.:_,uIrand substitutU ing in (5.5) verifies (5.7).

Theorem 2: For the growing memory, exponentially weighted pseudo inner product of (2.1), we have
P,,+-jU

VI.

LEAST-SQUARES

FILTERING ALGORITHMS AND SYSTOLIC ARRAYS

The vector algorithm, the three time update equations of Corollary 2, and the order update equation (5.10) for c , ( t ) Proof: Let U denote a member of the set given by provide all the tools needed to derive a large number of and least-squares filtering algorithms with systolic array architecthe left hand side of (5.4). Since P,,:-,u = U - P j , , - l ~ c , ( t ) E S’, it is easily seen that U is of the form U = U - a , tures. Quite a few algorithms can be derived, depending where a E Si. Thus, U E P;,;U if we can show ( u , , x ) / = 0, for upon the choice of which time update of Corollary 2 is used and whether k , , ( t ) is given by time updates, or from the all x E Si. But this is easily verified, for if x E SI,then order update for e , , ( ? ) as given in the vector algorithm. We (U,X)f = (P;;-Iu>x)f - [ p;,:-lu],(c;(t),x)/ shall derive only two algorithms, but other algorithms can be derived in a similar way [5]. Also, in [5], a normalized =P(p;.:-IU,x)l-l +[P;;-lu]/[Xlf* - [ P ; ; - l ~ l l [ ~ l / * algorithm is presented in which nearly all internal variables = 0, are in magnitude less than one, and thus may be especially where we applied (2.2) to ( P j , ~ _ l ~ ,and l(5.1) to ( c i ( t ) , x ) r . well suited for fixed point arithmetic. x) 0 Our first algorithm is relatively straight forward to derive. Coro/laty 2 For the growing memory, exponentially It is obtained by using the time update (5.7) to calculate k , ( t ) : weighted pseudo inner product of (2.1), and for any T-vec- and k , , ( t ) . Direct application of (5.7) to the defining equators U and U , the following time updates follow from Theo- tion for k , ( t ) in (4.1), yields rem 2: k , ( t ) = P k , ( t - 1) + 6 { r l - 1 ( ~ ) 4 , ( ~ ) * -(c , - l ( t > ) } ? / 1 ( U ,P , . : d / = P ( U , P;$ I l l ) / - I + [ P;; U ] [ P;,$ IU], , (5.5) * where we have made the definitions ( U ,P,.:U>/ P ( U , p- l U > f - I + [ P,:- IU] / [ p ; ; 4 = ; (5.6) r , ( t >= [ T , ( t ) I / ? P;.:.)/ = P ( ~ , f $ , u ) I - l h,,(t) = [ b , , ( t ) ] / . + 6 { [ P ; U / [ 4: U /*/( 1 - ci ( t ) 1} . i . (5.7) Similarly, applying (5.7) to (4.2) gives the time update for Proof: Forming the pseudo inner product of both sides k , , ( t ) , of (5.4) with U yields k , , W = p k , , ( t -I)+{ ~ , , ( ~ ) ~ , , ( f ) * / ( l - ~ , - l ( ~ ) ) } . ~

- [ P;,+-lu]/c;() c

1

1

>), +I;.
.,n,

4,+1),(4 b , , ( t ) =

s{k,,(r)/k,,(t))b,,(t),
] = I

+ I;..,n,

k , ( f )= p k , ( t - I ) + i ( , + l ) , ( f )= & , , ( t )-

r,-l(t)i,,(t)*,
-

W , ( t

l)/kl,(t

-

1 ) ] m > f l , ‘ . ‘ , H ,

c,(t) = ~ , - 1 ( ~ ) + 6 ( 1 ~ , , ( ~ ) 1 2 / ~ , 1 ( ~ ) } ~ r , ( t ) = r,-l(t)
-

J = l

s{k,(t)/k,,(t)}b,,(t>.

k,,W

=Pk,@ -

1) + b , , ( f ) h , l ( t ) * .

In this algorithm, k , , ( t ) and k , ( t ) are initialized as k,,(O) = k,(O)= 0. Our second algorithm IS obtained by using the time update (5.51, and thus the c , ( t ) , I = 1,2; . ., n are not needed. Direct application of (5.5) to (4.1) and (4.2) yields the time updates kl(j) =

j=i,i+I;..,n, b(,+ I ), ( t 1 = b,,( t 1 - 6{k , , ( t 1/ k l l ( t ) )bl,( t 1
?

j = i + l ; . . ,11 rl( 1

t 1 = r, - I( f 1 - 6 { k I ( t 1/ k I , ( 1 h,, t 1.

1

Pk,(t

-

11 + r l - l m l 1 ( t ) * >

k,,(t) =Pk,(t -l)+b,,(Qi,,(Q*?

Again, we initialize k,,(O) = k,(O) = 0. In practice, with regard to application of S ( x / y ) to both algorithms, the selection of 6 ( x / y ) = 0 is taken when 14’1 < E , where E is clearly .

KALSON AND YAO: CLASS OF LEAST-SQUARES FILTERING AND ALGORITHMS

49

each of the processing cells can perform their functions, and consequently the throughput is independent of n. During each clock time, all cells of the array are processing data. The overall latency (the time to process a set of inputs) is 2n clock cycles. This latency is illustrated in Fig. 1, where it is shown that y ( t +2n), xi(r +2n), i = 1,2,.. . , n are ready to enter the array and the residuals r;(t),i = 1,2,. . n have all been computed. The relative order in which these residuals are computed is also illustrated in Fig. 1, where it is shown that r , ( t )was found n - 1 clock times before r,,(f).Note that since the total number of processing cells of the array is O b 2 ) , the actual number of computations performed per each time update to solve the least-squares filtering problem is O(n2). e, VII.

SYSTEM iDENTIFlCATlON

The system identification problem can be generalized to that of finding the n sets of filter weights ( g T ( t ) , i = 1,2; . ., m ) for the various orders m = 1,2; . ., n, such that

I
1 ,

m

m

i=l

For the present, let us assume that k , , ( t )# 0 for all i (this is equivalent to the assumption that the sample data correlation matrix R ( t ) is full rank). The removal of this assumption will be discussed at the end of this section. As discussed in Section IV, b,,(t), i = 1,2; * n can be considered as an orthogonal basis for x,, i = 1,2,*.. , n in which case, under our current assumption, b , , ( t ) / d m , i = 1,2; . n would be an orthonormal basis. With this orthonormal basis, one can use standard techniques (such as "QR" decomposition [24]) for solving (7.1). To this end, let us define the n-vector d t ) with elements a , ( t ) equal to k,( t )/ k J t ), i = 1,2,. n and the n x n unit upper triangular matrix A t ) with elements A , , ( t ) equal to k i , ( f ) / k l l ( f ) , 1I i I j I n. Also, let g " ( t ) denote the m-vector with elements g,"(t), i = 1,2;..,m, the m-vector a m ( t ) be the first m elements of a ( t ) , and A"(t) be the m X m principal submatrix of A ( t ) . Then (7.1) is equivalent to the system of matrix equations a , a ,

chosen to be slightly greater than the smallest nonzero positive machine representable number. These least-squares filtering algorithms, and others, can be implemented by the generic two dimensional systolic array given in Fig. 1. The overall general structure of the communication paths between the processing cells in Fig. 1 are independent of which algorithms are implemented. However, since Algorithm 2 does not require calculation of c , ( t ) , it does not need the communication path with unit delay cells between the circular cells. The array consists of three types of processing cells whose functions depend upon which algorithm is implemented. As an example, the description of these cells for Algorithm 1 are given in Figs. 2, 3, and 4. In these figures, the time dependence is omitted for clarity. The time varying scalars k , ( t ) , k , , ( t ) for j > i, and k , , ( t ) are stored within the cells defined in Figs. 2, 3, and 4, respect ively . As shown in Fig. 1, the input data y ( r ) , x , ( t ) , i = 1,2;. ', n first past through delay cells. These delay cells are needed to skew the incoming data so that the array can work in a fully systolic manner. The throughput of the array is the rate that

a'"(t) = A m ( t ) g M ( t ) ,

m=1,2,...,n.

(7.2)

A systolic array solving (7.2), along with input and output, is illustrated in Fig. 5, where the square processing cells are described in Fig. 6. The circular cell in Fig. 5 is just a simple unit delay cell. Again, the time dependence is omitted in these figures for clarity. In Fig. 5, the first row of output, ( g f ( t ) ,g;(t),. g,"(t)), is available n 1 clock cycles after the input (0,1,1; 1) enters the array. The last output, g ; ( t ) , exits the array 2n clock cycles from when the input enters. Thus, the total latency for solving the system of matrix equations represented by (7.2) is 2n. The structure of the systolic array of Fig. 5 is essentially identical to the systolic array of Fig. 1, except delay cells and some communication paths are not needed. Thus, by proper control and management of memory, the systolic array used for least-squares filtering can also be used for system identification. For example, if the filter vector g " ( t ) is desired, then the least-squares filtering operation of the array in Fig. 1 must be suspended, and the array must be operated as

-

a ,

+

e,

SO

IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 37, NO. I . J A N U A R Y IY91

sliding memory cost criteria. In our derivation, we have generalized the definition of oblique projections so that the case of singular sample data correlation matrix can be considered. The construction of these algorithms follow easily from the general time and order update equations derived in this paper. We have presented only two such algorithms explicitly, but it should be straight forward to derive others. APPENDIX In this Appendix, we derive time updates for the sliding memory pseudo inner product. For notational convenience, we still use ( , >, for the sliding memory pseudo inner product, and since the growing memory case is not considered in this appendix, there should be no cause for confusion. For any U E C', U E C', we define ( u , u > , as f (U,U>f

=
T =f

C
-

1n

+I

Pr--TIUl-TIUl-~>

(A.1)

where p is the exponential weight factor as used in the growing memory inner product and m is the length of the memory. That is, data that is older than m units of time is completely ignored. Just as for the growing memory pseudo inner product, it is necessary to introduce some auxiliary vectors. Let ~ , ~ , ( t ) denote the T-vector defined as
U
Fig. 5. Systolic array for system identification.

Also, we let
Pl.fT/??(f)7

s,(t)

be some T-vector in the oblique projection s,(t> E P;,f%(f),

I1

If

f',,,,,

=

I,,,,

+

where we note that all oblique projections are now with respect to the pseudo inner product of (A.1). As in (5.11, we have for s , ( t ) and any T-vector U
Ail

(111

I,,,,,(

=

Fig. 6 . Cell descriptions for system identification.

( s,( f 1, U )

/

=

P+;(f)Ilf2.

PI.,U

I,*_

I,)

+I

'

(A.2) as ('4.3)

Furthermore, we define the real-valued scalar indicated in Fig. 5. The systolic array can therefore be "switched" between the modes of filtering and identification. So far in this section, we have assumed that k l l ( t )# 0 for all i. If this assumption is removed, then for sufficiently large order m , there will exist some nonvoid subset of {1,2;. . , m } for which k , , ( t )= 0 for i belonging to this subset. Denote this subset by N J t ) and its complement by N,,(t). If we set a , ( t )= G { k , ( t ) / k , , ( t ) } and A / , ( t )= G { k l , ( f ) / k l l ( f ) ) , then it can be shown that the system of matrix equations (7.2) solves the following reduced order least-squares identification problem sl(r> sl(t)

Just as for c , ( t ) , it is easy to verify that s , ( t ) satisfies 0 I s i ( [ )I 1. We shall need the following simple, but important lemma. Lemma 1: For any T-vector U , c , ( t ) = 1 implies [ P ! , + u ]= 0, , and s , ( t ) = 1 implies [PI.';U ] , ~,,, = 0. +

,

Proof: If c , ( t ) = 1, then from (5.2) and (5.3) we have
~ ~ c l = f1 ) and [~ i ( t ) ] = 1. Together, these last two equa( ~ ~ c , = tions imply [ c , ( t ) ] , 0, t - m l I 7 I t - l , and thus for

+

any vector

U

f
~Ll~r[Cl(t)]TIUIT*=[Ul/*~

(

c l ( t ) , u ) f =

T =r

-in

+I

for each m

=

1,2; . . , n .

However, from (5.1) and this equation, we have [ P i , r u ]=r[ U ] , , from which it follows that 0 = [ U ] , - [ P , , , u ] ,= [ P,,:u],. The first part of the lemma is seen to be true and the latter half U can be Droved in the same wav. The vector time update theorem ic now given. Tlzeorem 3: For the sliding memory pseudo inner product of (A.l), and tor any T-vector U , we have
PI,; - I~

VIII. CONCLUDING REMARKS
A unified approach is presented for deriving a large class of new and previously known least-squares filtering and system identification algorithms with systolic array architectures. Time updates are presented for both the growing and

+ As,( t - 1 ) - Bc,( t ) c PI: U ,

(A.4)

KALSON AND YAO: CLASS OF LEAST-SQUARES FILTERING AND ALGORITHMS

SI

where

Using both (A.14) and (A.15) simplifies (A.13) to

A

= p('?J I Wg( -

[ p, ; - l ~ ] , - l , i / ( ~ - ~ ; ( ~ - l ) ) } ~

(A.5)

(SI('

-

w),

=

[ s , ( t - l)l,[xIl*

B = [P;;-IU],+p(')J-I)/*
.6{ [ S i ( f - 1

+ p ( f 7 J + 1 ) / 2 [ ~ ] :l--,s,,,((r - 1 ) ) .
/
~

(A.16)

)I ,[ p,:

~

IU]

,,,/( 1 - S i ( t

- 1)

11.

(A.6)

Finally, using (A.12), (A.161, and (A.10) for ( P,.:-,u,x)t, ( s , ( t - l ) , ~ ) ,and ( c f ( f ) , x ) , , respectively, will verify (A.11). ,
I3

Proof: For the pseudo inner product of (A.1), we have the recursion relation
(U,U>=P(U,U)/-I /

+ [ u I / [ . l , * - PffJ[U1r -'n[~lr*-',l. (A.7)

From Theorem 3, we have the following corollary. Corollary 3: For the sliding memory pseudo inner product of (A. 1),

Consider first the case s J t - 1) = 1. Then from ( A S ) and (A.6), A = 0 and B = [Pj.:- Iu],, in which case (A.4) reduces to
P ; ; - l u - [ P ; , + - , u ] / C ; (c)) : . tP , .

( u , P ' ; u > /= p ( U , P , ; - I U ) / - l - p ( 1 i 1 + 1 ) / 2 A * [ P , ; _ I u ] , - , n

+ B* [ p,;

4,

7

('4.17)

with A and B given in (AS) and (A.6), respectively.

('4.8)

To verify (A.81, we need to show that for all x E S',
0 = r
=

Outline of Proof: For s , ( t - 1) = 1, the vector time update Theorem 3 reduces to (A.8), which is formally equivalent to the vector time update Theorem 2 for the growing memory pseudo inner product. In this case, the derivation of (A.17) is identical to that of (5.5) Corollary 2. For the case of s , ( t - 1 ) < 1, form the pseudo inner product of both sides of (A.4) with U to obtain
(u,Pl,:u)/ = ( U , P , , : - J ,

+A*(U,S,(t

-1)h - B * ( u , c , ( t ) ) r .

[P;.:-,u]/[xl/*.

From 6.11, we see that the second pseudo inner product on the right-hand side of (A.9) is, for x E S', ( . ; ( ~ ) ? ~ ) r = [ ~ , , r x=[XI/*. lf* (A.lO)

This equation can be shown to be equivalent to (A.17) by performing the following algebraic steps: Apply the recursion relation (A.7) to ( U ,P , , ~ - l u ) and ( v , s , ( t - l)),. The terms , ( u , s , ( t - 1))/-1 and [ s , ( t - l)l,-mwill appear in the expression for ( U ,s , ( t - l)),, which can be evaluated by using (A.2) and (A.lS), respectively. Finally, use (5.1) to evaluate
( U ,c , ( t ) ) , .

0

The previous two equations verify (A.9), which in turn verifies (A.8). The case of s,(f - 1) < 1 is not conceptually different. Again, (A.4) is verified provided we can show that for all x E S',
~ = ( P , , : . - ~ U , X )A + ~ j ( t - - l ) , ~ ) / - ( c ; ( t ) , x ) , . ( A . l l ) ,( B

To verify this, we need alternate expressions for the three pseudo inner products on the right-hand side of (A.11). Applying (A.7) to the first pseudo inner product, and again making use of ( Pj,+- I u ,x), - I = 0, for all x E S', one obtains
( P: ; .
-IU

Of course, other similar scalar time updates follow from (A.17) by interchanging the roles of U and U , complex conjugation, and using the self-adjoint property of the oblique projection operator. Finally, we present one more scalar time update that takes a simple form. Corollary 4: For the sliding memory pseudo inner product of (A.11, we have (u,P,;u)f
-p - (U,~,,:-l~)f-l

,

=

[ P;

-

Iu] / [ x I,* - .',I[

P ;

- 1111 t

- , , , [ X I ~ "n-. i (A.12)

+6
-P

{ [ 4;

4 ,[ P,;I ,*I(1U
-

C/( t

1) 1
1- S I ( - 1)) } . (A.18)

Applying (A.7) to the second pseudo inner product in (A.11) gives (s;(f-l),x)f =p(s;(f-1),X),-l+[s;(f
- Pf7'[Si(f -

f n d[ 4 ;

1 4 , ,n [ PI;-

IU]

,*- m / (

-l)]/[xl/*

1)1/ - l , l ~ ~ l : - l , l(A.13) .

Outfine of Proof To verify (A.18), let us first assume that c , ( t ) and s , ( t ) are less than one. Then, taking the tth component of both sides of (A.4) in Theorem 3 yields after some algebra [P,,+I//(lc , ( t ) ) = B.

This equation can be simplified as follows. Using (A.2) with replaced by x and t replaced by t - 1 (remembering that x E SI), we have
U

( s i ( f - 1 ) , x ) , - , =p("1-'1)/2[x]1*-f1)1. (A.14)
Also, from (A.2) and (A.3) it is easy to see that ~j(~-l)=(sj(~-l),sf(r-l)),-l p ( " J - 1 ) " 2 [ ~-(If) ] , l

Substitution of this expression for B and (AS) for A into (A.17) yields (A.18). The cases of one or both of c , ( t ) and s , ( t ) equal to one can be verified in a similar way, noting that I3 Lemma 3 must be used where appropriate.

(A.15)

Remark 4: The scalar time update equations of Corollaries 3 and 4 can be used to derive sliding memory leastsquares filtering algorithms in essentially the same way as the

52

IEEE T R A N S A C T I O N S O N I N F O R M A T I O N T H E O R Y . VOL. 37. NO. I . J A N U A R Y l < J Y I

scalar time update of Corollary 2 was used to derive growing memory least-squares algorithms, [5]. The additional order updates for [ s , ( t - l)], and s , ( t ) are required, which are easily derived from the vector order update for s,.This order update is given below, and is easily obtained from Theorem 1 where U is replaced with r , J t ) .

ACKNOWLtDGMtN I

The authors would like to thank the reviewers for their comments and suggestions, which have helped to make this paper more readable and precise.

REFEKENCES
T. Kailath, “ A view of three decades of linear filtering theory,” IEEE Trans. Inform. Theory, vol. IT-20, no. 2, pp. 146-181, Mar. 1974. N. Wiener, Time Series. Cambridge, MA: The M.I.T. Press, 1977. M. J. Shensa, “Recursive least-squares lattice algorithms: A geometrical approach,” IEEE Trans. Airtomat. Contr., vol. AC-26, pp. 695-702, June 1981. H. Lev-Ari, T. Kailath, and J. Cioffi, “Least-sqkares adaptive lattice and transversal filters: A unified geometric theory,” IEEE Trans. Inform. Theory, vol. IT-30, no. 2, pp. 222-236, Mar. 1984. S. Kalson, “Recursive least-squares filtering algorithms with systolic array architectures: A geometrical approach,” Ph.D. dissertation, Univ. of California, Los Angeles, 1985. M. Morf, T. Kailath, and L. Ljung, “Fast algorithms for recursive identification,” in Proc. 1976 C o n t Decision Conrr., Clearwater Beach, FL, Dec. 1976, pp. 916- 321. M. Morf, A. Vierra, and D. T. Lee. “Ladder forms for identification and speech processing,” i n Proc. 1977 IEEE Cotit Deckion and Contr.. New Orleans, LA, Dec. 1977, pp. 1074- 1078. L. Ljung, M. Morf, and D. D. Falconer, “Fast calculation of gain matrices for recursive estimation schemes,” Int. J. Contr., vol. 27, pp. 1-19, 1978. D. Lee, M. Morf, and B. Friedlander, “Recursive least-squares ladder estimation algorithms,” IEEE Trans. Acoirst. Speech Signal Processing, vol. ASSP-29, no. 3, pp. 627-641, June 1981.

[IO] B. Porat and T. Kailath, “Normalized lattice algorithms for least-squares FIR system identification,” IEEE Trtrns. Acoir.~. S p c w l r Signtrl Procc~ssirrg, ASSP-3 I , no. I , pp. 122- 128, Feb. vol. 1983. [ I I ] F. Ling and J. G. Proakis, “A generalized multichannel leastsquares lattice algorithm based o n sequential processing stages,” IEEE Trans. Acoirst . Speech Signul Processing, vol. ASSP-32, no. 2, pp. 381-38‘). Apr. 1984. [ 121 J. Cioffi and T. Kailath, “Fast, recursive-least-squares transversal filters for adaptive filtering,” IEEE Trans. Acoirst. Sp:.c~!i Signal Processing, vol. ASSP-32, no. 2, pp. 304-337, Apr. 1984. [I31 T. Kawase, H. Sakai, and H. Tokumaru, “Recursive leastsquares circular lattice and estimation algorithms,” IEEE Trans. Acoirst. Speech Signal Processing, vol. ASSP-31, no. I , pp. 228-231, Feb. 1983. [14] W. M. Gentleman and H. T. Kung, “Matrix triangularization by systolic arrays,” Proc. SPIE, Keal-Time Signal P w c . IV. vol. 298, 1981, pp. 19-26. [15] J. G. McWhirter, “Recursive least-squares minimization using a systolic array,” in Proc. SPIE Real Time Signal Proc. VI, vol. 431, 1983, pp. 105-1 12. [ I61 H. Lev-Ari, Wave-front array processing architectures for multichannel ladder algorithms,” Integrated System Inc., Palo Alto, CA, technical memo 5016-07, June 1982. [I71 H. T. Kung and C. E. Leiserson, “Systolic arrays (for VLSI),” Sparse Matrix Proc. 1978, I. S. Duff and G. W. Stewart, Eds., SIAM 1979, pp. 256-282. [18] S. Y. Kung, K. S. Arun, R. J. Gal-Ezer, and D. V. Bhaskar Rao, “Wavefront array processor: Language, architecture, and applications,” IEEE Trans. Comput. (Special Issue on Parallel und Distributed Computers), vol. C-31, no. 11, pp. 1054-1066, Nov. 1982. [19] S. Y. Kung and R. J. Gal-Ezer, “Synchronous vs. asynchronous computation in VLSI array processors,” in Proc. SPlE Conf., Arlington, VA, 1982. [20] F. Ling, D. Manolakis, and J. G. Proakis, “A recursive modified Gram-Schmidt algorithm for least-squares estimation.” IEEE Trans. Acoust. Speech Signal Processirrg, vol. ASSP-34, no. 4, pp. 829-836, Aug. 1986. [21] S. Kalson and K. Yao, “Systolic array processing for order and time recursive generalized least-squares estimation,” in Proc. SPIE Real Time Signal Proc. VIII, vol. 564, 1985, pp. 28-38. [22] H. L. Royden, Real Analysis Second ed. New York: MacmilIan Publishing Co., Inc., 1968. [23] G. Cybenko, “A general ortliogonalization technique with applications to time series analysis and signal processing,” Math. of Comput., vol. 40, no. 161, pp. 323-336, Jan. 1983. [24] L. Lawson and R. J. Hanson, Soking Least-Squares Problems. Englewood Cliffs, NJ: Prentice-Hall, Inc., 1974.
“

c.