Fast Adaptive Identification of Stable Innovation Filters

F ast Adaptiv e Iden tifi ca ti o n of Stable Inno v ation Filters A. P . Mullhaupt and K. S. Riedel Coura n t Inst itute of Mathem atica l Sciences New Y ork Univ er sity New Y ork, NY 1 0012 -1185 Abstract The adaptiv e id en tification of th e impu lse resp onse of an inno v ation filter is con- sidered. The impulse resp onse is a finite su m of kno wn basis fun ctions with unknown co efficien ts. These unknown coefficien ts are estimated usin g a pseudolinear regression. This estimate is i m p lemen ted using a square ro ot alg orithm based o n a d isplacemen t rank structur e. When the initial conditions ha ve lo w displacemen t rank, the fi lter up d ate is O ( n ). If the filter arc hitecture is chosen to b e triangular input balanced, the estimation problem is w ell-conditioned and a s imple, low r ank initializa tion i s a v ailable. 1 INTR ODUCTION W e consider inno v ation mo dels f or state space systems with a n unknow n and p ossibly time dep enden t impulse resp onse. Inno v ation mo dels use the prediction fit errors a s the sto c hastic input in to the state space ev olutio n and the feedbac k/gain matrix is estimated empirically . W e use a pseudolinear regression (PLR) [LS, Section 3.7.3] t o id entify the unkno wn co efficien ts of t he state space impulse resp onse. PLR is a least squares estimate in whic h t he unkno wn co efficien ts are recursiv ely up da t ed from past estimates of the r esiduals. W e refer the reader to [LS] for an exce llent exp o sition of the theoretical prop erties o f inno v ation mo dels and PLR. W e show that inno v ation filter systems p ossess a displaceme nt structure [KS]. Using this displace ment structure, w e construct a fast O ( n ) square ro ot filter. Our sq uare ro ot 1 displacemen t filter (SRDF) effects a time-dep enden t change of co ordinates that preserv es the imp ulse resp onse while transforming the empirical co v ariance matrix to the identit y . A canonical represe ntation of the state space system is used to simplify the SRD F initialization and to improv e the condition of the estimation of these parameters. In Section 2, the inno v ations filter represen tation of the impulse r esp o nse is giv en and pseudo-linear regression is used to estimate the unknown co efficien ts. Se ctions 3 and 4 presen t a fast square ro ot v ersion of the PLR estimate that is based on lo w displacemen t rank. Section 5 desc rib es a new matrix c anonical represen tation, triangula r input b a la nced (TIB) form and its applications to filter arc hitecture. When the system adv ance matrix is in TIB form, the system is easily initia lized and is alwa ys w ell-conditioned. Sections 6- 7 discuss and summarize our results. Except whe re explicitly noted, our results apply for the m ultiple input m ultiple output (MIMO) case. 2 INNO V A TIONS FIL TER F OR SYS TEM IDENTI- FICA TION Let y t b e a sequence of d -dimensional measuremen t v ectors. W e consider the syste m iden tification problem of asso ciating to y t an unknow n m ulti-input/m ulti-o utput (MIMO) time inv ariant linear system of the form: z t +1 = Az t + b ǫ t , (2.1) y t = c ∗ z t + ǫ t , (2.2) where ∗ denotes the Hermitian tra nsp ose. Here z t is the n -dimensional stat e ve ctor , A is the n -dimensional system matrix, and b and c a re n × d -dimensional matrices. Equations (2.1) and (2.2) are in i n novations mo del / pr e diction err or mo del form, i.e. the measuremen t noise and the system noise are completely correlated via the d -dimensional inno v ations v ector ǫ t . Ev ery stable linear t ime-inv ariant n -dimensional state space system ma y b e represen ted in this f orm [HD]. Our g o al is to iden tify the underlying impulse resp onse of the system. It mak es sense to assume that ( A , b , c ) is minimal , but instead we will simply assume that it is c ontr ol lab l e . In addition, w e will assume that A is stable and, to reduce the complexit y of exposition, nonsingular . This implies that ( A , b ) is completely realizable. 2 The impulse respo nse is preserv ed b y the c hange of state space co ordinates: A → T A T − 1 , b → Tb , c → T −∗ c , (2.3) where T is an arbitrary nonsingular n × n mat r ix. Since the impulse resp onse, h j ≡ c ∗ A j b , is preserv ed by this c hang e of co ordinates, we ma y c ho ose an observ ationally equiv alent state space system in whic h the system adv ance has desirable prop erties. In particular, by an appro pr ia te c hoice of T , w e can represen t the system in TIB f o rm. ( See Sec. 5.) In practice, w e prescrib e t he eigenv alues of A to corresp o nd to the c har acteristic re- sp onse times of t he syste m. W e then c ho ose ( A , b ) to b e in TIB fo rm. This uniquely sp ecifies ( A , b ) up to the transformation group: ( A , b ← (Θ A Θ ∗ , Θ bU ), where U is a d × d unitary tra nsformation and Θ is a n × n diag onal unitar y transformation. Th us, when the eigen v alues are prescrib ed, only c need b e estimated to identify the impulse resp onse. A t eac h time step, w e observ e y t , a nd estimate the co efficien ts using the PLR o f [LS]. W e denote the estimate of z t and c using data up to a nd including y t b y ˆ z t and ˆ c t . The predictiv e residual error, ˆ ǫ t = y t − ˆ c ∗ t − 1 ˆ z t , is ev a lua ted, where ˆ z t is the estimate of z t giv en b y using the t ime adv ance equation, (2 .1) with the substitutions, z t ← ˆ z t , ǫ t ← ˆ ǫ t and c ← ˆ c t − 1 . Let δ b e a ‘forgetting factor ’ with 0 < δ ≤ 1 and 1 − δ << 1. W e define the (w eigh ted) empiric a l co v ariance ˆ P δ t and the empirical cross-co v ariance d δ t : ˆ P δ t ≡ δ t ˆ P δ 0 + t X k =1 δ t − k ˆ z k ˆ z ∗ k , d δ t ≡ δ t d δ 0 + t X k =1 δ t − k ˆ z k y ∗ k , (2.4) where ˆ P δ 0 is t he initial cov ariance, and the initial cross-cov ariance is d δ 0 = ˆ P δ 0 ˆ c 0 . The PLR estimate of c is ˆ c t = h ˆ P δ t i − 1 d δ t . Note that the v alue of y t is used in the estimate ˆ c t . In summary , the PLR estimate of (2.1)-(2.2) giv en the measuremen ts { y s | 1 ≤ s ≤ t } is ˆ z t +1 = A ˆ z t + b ˆ ǫ t , (2.5) ˆ ǫ t = y t − ˆ c ∗ t ˆ z t , (2.6) ˆ c t = h ˆ P δ t i − 1 d δ t . (2.7) In the “pre-windo w ed” case, the filter is initialized as ˆ z 1 ≡ ˆ z t =1 = ˆ c t =0 = 0. In order to pr e d ict y t +1 , we may use ˆ y t +1 | t = ˆ c ∗ t ˆ z t +1 and ˆ z t +1 = A ˆ z t + b  y t − ˆ c ∗ t − 1 ˆ z t  . The empirical co v ariances satisfy ˆ P δ t = δ ˆ P δ t − 1 + ˆ z t ˆ z ∗ t , d δ t = δ d δ t − 1 + ˆ z t y ∗ t . (2.8) 3 Let ˆ P δ 0 b e inv ertible and define Φ t = h ˆ P δ t i − 1 . By the matrix in v ersion identit y , Φ t = δ − 1 h Φ t − 1 − Φ t − 1 ˆ z t ˆ z ∗ t Φ ∗ t − 1 δ + ˆ z ∗ t Φ t − 1 ˆ z t i , and the PLR up date for ˆ c t is ˆ c t = ˆ c t − 1 − Φ t − 1 ˆ z t  y t − ˆ c ∗ t − 1 ˆ z t  ∗ δ + ˆ z ∗ t Φ t − 1 ˆ z t . (2.9) The PLR estimate determines ˆ c t in terms of the previous estimates of the prediction errors ˆ ǫ t = y t − ˆ c ∗ t − 1 ˆ z t . F or large times , t > 1 / (1 − δ ), the exp ectation of empirical co v ariance tends to E[ ˆ P δ t ] = P ∞ / (1 − δ ), where P ∞ satisfies St ein’s equation [L T ]: P ∞ − A P ∞ A ∗ = σ 2 bb ∗ . (2.10) 3 DISPLA CEMENT STR UCTURE of t h e CO V ARI- ANCE In Section 4, w e describ e a fa st square r o ot algorithm for computing ˆ c t . The fa st algo rithm is based on t he displac ement structur e of ˆ P δ t . This is similar to, but not the same as, the fast up date of Sa y ed and Kailath [SK1]. Here, w e compute a generator and its signature for the displacemen t of the empirical cov ariance, ˆ P δ t − δ − 1 A ˆ P δ t A ∗ , this determines an upp er b ound on the ra nk o f minimal generators for this displacemen t. Any matrix X t suc h that ˆ P δ t − δ − 1 A ˆ P δ t A ∗ = X t SX ∗ t . (3.1) is called a generator of the displacemen t, and S is a corresp onding signature. W e assume that a generator X 0 and signat ure S 0 for the initia l cov ariance, X 0 S 0 X ∗ 0 = ˆ P δ 0 − δ − 1 A ˆ P δ 0 A ∗ , ha v e b een computed f r o m ˆ P δ 0 b y the singular v alue decomposition o r b y more sp ecialized metho ds [SK1, SK2], see also Section 5. By substituting (2.4) in to (3 .1) a nd collapsing the summation using (2.1), w e hav e X t SX ∗ t = − δ − 1 A ˆ z t ˆ z ∗ t A ∗ + δ t − 1 ˆ z 1 ˆ z ∗ 1 + δ t X 0 SX ∗ 0 + t − 1 X j =1 δ t − j − 1  b ˆ ǫ j ˆ z ∗ j A ∗ + A ˆ z j ˆ ǫ ∗ j b ∗ + b ˆ ǫ j ˆ ǫ ∗ j b ∗  = − δ − 1 A ˆ z t ˆ z ∗ t A ∗ + δ t − 1 ˆ z 1 ˆ z ∗ 1 + ˆ g t − 1 ˆ g ∗ t − 1 − ˆ h t − 1 ˆ h ∗ t − 1 + δ t X 0 S 0 X ∗ 0 , (3.2) where w e ha v e define d the n × d -dimensional matrices: ˆ f t ≡ t X j =1 δ t − j  A ˆ z j + 1 2 b ˆ ǫ j  ˆ ǫ ∗ j , ˆ g t ≡ 2 − 1 / 2  ˆ f t + b  , ˆ h t ≡ 2 − 1 / 2  ˆ f t − b  . (3.3) 4 Let I d b e the d - dimensional iden tity matr ix. W e define the signature matrix S = ( − 1) ⊕ 1 ⊕ I d ⊕ ( − I d ) ⊕ S 0 , and define the dis plac ement r ank , α ≡ rank( S ). One particular c hoice of g enerator is X t =  δ − 1 / 2 A ˆ z t | δ ( t − 1) / 2 ˆ z 1 | ˆ g t − 1 | ˆ h t − 1 | δ t/ 2 X 0  . (3.4) This sho ws that X t ma y b e c hosen with rank ( X t ) ≤ α ≡ rank( S ). Note tha t α dep ends on the initialization of t he filter and it is desirable that the initial displacemen t ha v e hav e minimal rank. The case ˆ z 1 = 0 is called the pr ewindowe d case. As described in Section 5, w e recommend c ho osing the initial conditions such that ra nk ( X 0 ) = 1. F rom (2 .8), an y pair of g enerators, X t and X t +1 , satisfy the up date equation: X t +1 SX ∗ t +1 = δ X t SX ∗ t + ˆ z t +1 ˆ z ∗ t +1 − A ˆ z t +1 ˆ z ∗ t +1 A ∗ . (3.5) T o o bt a in a minimal rank generator of the displacemen t, w e p erform a n eigen v ector de- comp osition (EVD) of the up dated displacemen t: U t +1 Λ t +1 U ∗ t +1 = δ X t SX ∗ t + ˆ z t +1 ˆ z ∗ t +1 − A ˆ z t +1 ˆ z ∗ t +1 A ∗ , (3.6) where U t +1 is unitary and Λ t +1 is diagonal. W e choose X t +1 = U t +1 | Λ | 1 / 2 t +1 where the eigen v alues of Λ t +1 are ordered to resp ect the signature matrix, S . W e refer to this choice of generator as the E V D − induced g enerator. Giv en the EVD of X t SX ∗ t , the EVD of (3.6) ma y b e computed in O ( α 2 n ) op erat io ns using the up date/do wndate alg orithm of [GE]. If A is a structured matrix suc h that Az t is computable in O ( n ) op erations, then the displacemen t generator, X t , can b e up da ted in O ( α 2 n ) op erations. The EVD up da te of X t +1 ma y be replaced b y an y alternativ e update of the generator whic h is computable in O ( α 2 n ) flops and whic h guaran tees that rank X t +1 ≤ α . W e recommend the EVD up date of X t +1 to minimize k X t +1 k and preserv e the ortho g onalit y of t he columns of X t +1 . The EVD up date/dow ndate of [GE] is backw ards stable and ma y b e computed to mac hine precision. 4 SQUARE R OOT DISPLA CE ME NT FIL TER W e presen t a square ro ot displaceme nt v ersion of the PLR updat e of c t when A is upp er triangular (UT). The lo w er triangula r case is similar. The standard Cholesky up date of c t using (2.9) fo r a system requires O ( n 2 d ) time and O ( n 2 ) space, while our v ersion has O ( nα 2 ) time up dates. Ho we ver, the straig h tforward initializatio n of the SRDF is 5 computationally in tensiv e, O ( n 3 ), in the general case. In the next section, w e describ e a no v el filter architec ture that can b e initialized in O ( α 2 n ) op erations and that is alw ays w ell-conditioned. W e define R t and W t b y ˆ P δ t = R t R ∗ t , W t = R − 1 t R t − 1 , (4.1) where R t is upp er triangular and R t,j j > 0 for 1 ≤ j ≤ n . W e mak e the time dep endent c hange of co ordina t es: u t = R − 1 t − 1 ˆ z t , ´ c t = R ∗ t − 1 ˆ c t . The transformed state space mo del is u t +1 =  R − 1 t AR t  W t u t + R − 1 t b ˆ ǫ t , (4.2) y t = ´ c ∗ t u t + ˆ ǫ t . (4.3) By h yp othesis, R t and A are UT, R − 1 t AR t is UT and the dia gonal eleme nts of R − 1 t AR t matc h those of A :  R − 1 t AR t  j j = A j j . The co efficien t update is ´ c t = W −∗ t − 1 ´ c t − 1 + u t  y t − ´ c ∗ t − 1 u t  ∗ δ + u ∗ t u t . (4.4) These up dates can b e carried o ut in O ( n ) time and space as follow s: Conjugating equation (2.8) b y R − 1 t − 1 yields W ∗ t W t = δ − 1  I − u t u ∗ t δ + u ∗ t u t  , W − 1 t W −∗ t = δ I + u t u ∗ t , (4.5) where I is the n × n identit y matrix. Since W t is UT with p ositiv e diagonal elemen ts, it is a TIB mat rix in the sense of Sec. 5 and it is uniquely determined by (4.5). W e define β t = R − 1 t b and up dat e it by β t = W t β t − 1 . Conjugating equation (3 .1 ) by R − 1 t and defining Y t = R − 1 t X t yields  R − 1 t AR t   R − 1 t AR t  ∗ = δ ( I − Y t SY ∗ t ) . (4.6) Giv en Y t and the diagonal elemen ts of A (or at least t he complex phases of the diag o nal elemen ts), R − 1 t AR t is uniquely determined by (4.6). In the app endix, w e demonstrate that a represen tatio n of R − 1 t AR t can b e calculated from Y t in O ( α 2 n ). The resulting represen tation o f R − 1 t AR t allo ws the matrix v ector pro duct,  R − 1 t AR t  v , to b e computed in O ( α n ) time, where v is an arbitra ry n -ve ctor. Conjugating (3.5) by R − 1 t yields the displacemen t up date in the tra nsfor med co ordi- nates: Y t +1 SY ∗ t +1 = δ W t +1 Y t SY ∗ t W ∗ t +1 + W t +1 u t +1 u ∗ t +1 W ∗ t +1 − δ − 1 µ t +1 ∗ µ t +1 , (4.7) 6 where µ t +1 ≡ W t +1 ˜ A t u t +1 with ˜ A t ≡ R − 1 t AR t , T o enhance the n umerical stabilit y , w e c ho ose the EVD- induced generator: Y t = V t | D t | 1 / 2 where V t is unitary and the eigen- decomp osition of the displacemen t is Y t SY ∗ t = V t D t V ∗ t . S imilar to (3.6), we compute Y t +1 b y up da ting the EVD of (4.7 ) . W e claim t ha t W t = δ − 1 / 2 ( I − γ t u t u ∗ t ) Q t , where Q t is an unita ry matrix a nd γ t satisfy 2 γ t − γ 2 t k u t k 2 = 1 / ( δ + k u t k 2 ). (Pro of: Since δ − 1 / 2 ( I − γ t u t u ∗ t ) and W t are b o th square ro o ts o f the same matrix defined in (4.5), Q t is an unitary mat r ix.) Th us W t +1 Y t = δ − 1 / 2  I − γ t +1 u t +1 u ∗ t +1  Q t +1 Y t . The EVD of W t +1 Y t SY ∗ t W ∗ t +1 is a rank tw o p erturbation of the kno wn EVD of Q t +1 V t D t V ∗ t Q ∗ t +1 . Since W t +1 is a rank-one TIB matrix, Q t +1 V t ma y be rapidly ev aluated in in O ( α n ) op erations as δ 1 / 2  I − γ t +1 u t +1 u ∗ t +1  − 1 W t +1 U t . Th us the eigendecomp osition of (4.7) is a rank four p erturbation of the kno wn eigen v ectors, Q t +1 V t and ma y b e ev aluated in O ( α 2 n ) op erations (or few er if the fast m ultip ole algorithm is used) using the n umerically stable algorithm of [GE]. The eigendecomp osition-induced up date minimizes t he norm o f Y t and ensures that(4.6) ha s a positive semi-definite solution. Although the time up date of the SRDF is O ( α 2 n ), the straightforw ard initialization requires O ( n 3 ) op erations to compute the initial v alue of Y t . The initialization reduces an arbitrary inno v ations filter to an equiv alen t problem o f displacemen t rank α . Equation (3.2) sho ws tha t the mapping preserv es this lo w displacemen t rank as the filter ev olve s. Initialization algorithms for related displacemen t structures are describ ed in [SK1, SK2]. T o simplify the initialization, we select the initial co v ariance, ˆ P δ t =0 , suc h that ˆ P δ 0 − δ − 1 A ˆ P δ 0 A ∗ = ρ 2 bb ∗ , (4.8) where ρ is a scalar free parameter whic h determines the size of the initial co v ariance. The resulting SRDF initialization is Y 0 = ρ b . Although Y 0 do es not require the solution of (4.8), the solution o f (4.8) is required for β 0 = R − 1 0 b and u 1 = R − 1 0 z 1 . Giv en the cost of the initialization, t here are three cases where the SRDF will b e of v alue: 1) when the num ber of time steps is large relativ e to n 2 ; 2 ) when the initialization ma y b e done o ff - line and real time p erformance is imp ortant; 3) when the stat e space mo del has a sp ecial structure whic h simplifies the initialization. In the next section, we concen trate on filter arc hitectures which simplify the initialization. 7 5 TRIANGULAR INPUT BA L A NC ED F ORM The initialization of the square ro ot displacemen t filter may b e simplified b y c ho osing the state space co ordinates so that the system is in tria ng ular input balanced (TIB) for m. As described in [MR], the system, ( A , b ), is in UTIB form if I − AA ∗ = σ 2 r bb ∗ , (5.1) with A upp er triang ular and σ 2 /r > 0 arbitr a ry sub ject to σ 2 k b k 2 /r ≤ 1. TIB sys tems ha v e a n umber of adv an tag es [MR]: F ir st, the app endix sho ws ho w to construct a sequence of r ank-one TIB systems F k suc h A = F 1 . . . F d in O ( nd ) op erations. Giv en this represen- tation, m ultiplying A b y a v ector or solving a matrix system is O ( nd ). Second, the solution of Stein’s equation, (2.10), is a m ultiple of the iden tity , P ∞ = r I . Third, the solution of (4.8) is ˆ P δ t =0 = ( µ 2 r /σ 2 ) I + O ((1 − δ ) ). This expansion allows a fast initialization via (4.8) if w e drop the O ((1 − δ )) correction. F ourth, the exp ectation of empirical co v ariance tends to E[ ˆ P δ t →∞ ] = P ∞ / (1 − δ ) = r I / (1 − δ ). Th us, the iden tification o f ˆ c t b ecomes w ell-conditioned asymptotically . The or em 1 [MR]: Ev ery exp onentially asymptotically stable, completely realizable system ( A , b ) is similar to a UTIB sys tem ( A ′ , b ′ ), in whic h the diag onal elemen ts of A ′ can b e c hosen in an y order. The impulse respo nse , h j ≡ c ∗ A j b , of the PLR inno v ations filter is a linear comb inat ion of the n ro ws of t he n × ∞ matrix M ≡ [ b , Ab , A 2 b , ... ]. F or TIB systems, the n basis functions, M i, · , of the impulse resp onse are ortho no rmal: MM ∗ = ∞ X j =0 A j bb ∗ A ∗ j = 1 σ 2 P ∞ = r I . (5.2) F or SISO systems, when the TIB form A has only one eigen v alue (and therefore only one Jordan blo c k), the row s of M a r e given b y orthono r mal functions of the form λ t p k ( t ) . The resulting p olynomials, p k , ar e the Meixner p olynomials [NSU], a discrete analo g o f the Laguerre p olynomials. This orthonormalit y impro ve s the conditioning of the estimates of the co efficien ts, ˆ c t . Given the advantages of TIB form, we r e c ommend cho osing the T IB a r chite ctur e when the filter desig ner knows the e i g envalues of the system r esp onse and is no t c onstr aine d to a p articular state matrix. The following theorem pro vides an explic it construction of TIB forms for t he case of a SISO system: 8 The or em 2[MR, GGMS]: Supp o se ( A , b ) are UTIB and A is nonsingular then A j j = λ j and A ij = g i b j for i < j and 1 ≤ i, j ≤ n , where λ j = | λ j | e iθ j , | λ j | 2 = 1 + α j | b j | 2 and g j = e iθ j α j ¯ b | λ j | − 1 . The α j satisfy α j +1 = α j | λ j | − 1 , α j = α 1 − P j − 1 k =1 | g | 2 k and α 1 = − 1. This theorem enables the user to sp ecify the eigenv alues of A and c ho ose the r emainder of the filter arc hitecture suc h that A is TIB. An alternativ e construction o f the TIB f o rm sp ecifies the diagonal en tries of A and then determines a sequence of n hyperb o lic Giv ens rotations [GL] so that A satisfies (5 .1) [MR]. This factorization can b e o bt a ined from a P otap ov cascade embedding as in [LK]. W e anticipate that the TIB arc hitecture will hav e man y additional applications in filter design. In our applications, our goal is to iden tify the impulse resp onse, a nd the eigen v alues of A sim ply determine a basis for mo deling the impulse resp onse. Th us the filter designer need only choose the eigen v alues accurately enough to capture the characteristic time scales of the impulse resp o nse. 6 LMS and CHANDRASEKHAR FIL TERS The least mean square error ( L MS) estimate replaces Φ t b y the iden tity matrix in (2.9 ): ˆ c t = ˆ c t − 1 − µ t z t ǫ ∗ t . The p erformance of the LMS algorithm dep ends strongly on the condition nu mber of the stat e co v ariance. If A is a TIB filter, then E [Φ t ] tends to a multiple of the iden tity matrix plus a random fluctuatio n of relativ e size (1 − δ ) / (1 − δ k A k 2 ). Th us the LMS update sacrifices little accuracy since the second order correction is small. Th us the TIB arc hitecture a llo ws b o th the PLR a nd LMS up dates to p erform comparably with the adv an tages of b oth filters: well-conditioning, (near) second order accuracy , and O ( n ) time a dv ances . Assuming tha t the innov ations w ere zero, the fa st filter of Say ed and Kailath [SK1] can up date the displacemen t of ˆ P δ t +1 − A ˆ P δ t A ∗ . In con trast, our f a st filters are based on up dates to the displacemen t ˆ P δ t − δ − 1 A ˆ P δ t A ∗ . The tw o displacemen ts ar e related through (3.1). If the innov ations ǫ t w ere zero, then our filter w ould satisfy the structural conditions of [SK 1]. The presence of the inno v ation, b ǫ t , in the stat e up dat e introduces an additional lo w rank term. W e ha v e not compared our fast filter with the g eneralized Chandrasekhar filter of [SK1] since w e do not know how to treat the inno v ations ǫ t b in t he generalized Chandrasekhar filter. Chandrasekhar filters a re commonly used in conjunction with finite impulse resp onse (FIR) (time la g) mo dels. FIR mo dels hav e the disadv an tage that the c hara cteristic time 9 of the impulse resp onse is b ounded b y the num b er of time lag s/free par a meters. As the duration of the impulse r esp o nse b ecomes longer, more time lag co efficien ts need to b e iden tified. In contrast, the innov ation filters (2.1) ha v e the impulse resp onse duration con trolled by the eigen v alues of A while the num b er of free pa rameters is c hosen b y dim( A ) . This freedom a llows the filter designer to sp ecify the c hara cteristic time scales of the impulse resp onse and then to optimize the n um b er of terms in the orthonormal basis expansion of the impulse resp o nse. 7 SUMMAR Y W e estimate the co efficien ts of the impulse resp onse of the inno v ations f orm of a state space system using pseudolinear regression. In our case, this is iden tical to using recursiv e least squares when the innov ations a re estimated sequen tially and not up dated. Th e key adv an tage of ada ptiv e estimation using the innov ations filter is that the estimation of c t is a linear problem. The state space system has small displacemen t rank if the initial system do es. The co efficien t estimation is implemen ted with a square ro ot displacemen t filter. The square ro ot displacemen t filter constructs a state space co ordinate c hange at eac h time step to appro ximate TIB form. The square roo t displacemen t filter dep ends on the state matrix, A , through the current v alues of R − 1 t AR t and u t or equiv alen tly Y t and u t . The time adv ance r equires O ( α 2 n ) op erations. The difficulty in a pplying the SRDF is that the starting v alues of Y t and u t are expensiv e to compute for an arbitrary inno v ations filter. When the state system is in tr ia ngular input bala nced fo r m, an approximate initialization is a v ailable in O ( n ). If the filter designer kno ws the eigen v alues of A , Theorem 2 constructs the TIB form for SISO systems. This corr esp o nds to expanding the impulse resp onse in a n orthonormal basis of exp onen tial deca ying p olynomials, thereb y improv ing the conditio ning of the co efficien t estimate, ˆ c t . In practice, the eigen v alues of A need no t b e kno wn exactly , but only well enough that the series represen tation of the impulse respo nse is r easonable. The fast square ro ot displacemen t PLR filter is: ˆ ǫ t = y t − ˆ c ∗ t − 1 u t (7.1) ˆ c t = W −∗ t − 1 ˆ c t − 1 + u t ˆ ǫ t δ + u ∗ t u t (7.2) β t = W t β t − 1 (7.3) 10 u t +1 = ˜ A t W t u t + β t ˆ ǫ t (7.4) W ∗ t +1 W t +1 = δ − 1  I − u t +1 u ∗ t +1 δ + u ∗ t +1 u t +1  , (7.5) W − 1 t +1 W −∗ t +1 = δ I + u t +1 u ∗ t +1 (7.6) V t +1 D V ∗ t +1 = δ W t +1 V t D V ∗ t W ∗ t +1 + W t +1 u t +1 u ∗ t +1 W ∗ t +1 − δ − 1 µ t +1 ∗ µ t +1 (7.7) ˜ A t +1 ˜ A ∗ t +1 = δ  I − Y t +1 SY ∗ t +1  (7.8) where V t is unitary , Y t = V t | D t | 1 / 2 , ˜ A t = R − 1 t AR t , µ t +1 ≡ W t +1 ˜ A t u t +1 . The up date of the eigendecompo sition is p erformed using W t +1 V t = δ − 1 / 2  I − γ t +1 u t +1 u ∗ t +1  ˜ V t where ˜ V t = δ 1 / 2 ( I − γ t u t u ∗ t ) − 1 W t +1 is unitary . This algorithm may b e mo dified b y replacing the prediction error, ˆ ǫ t , with a p osteriori errors (residuals), ǫ t = y t − ˆ c ∗ t u t , in the up date of ˆ c t [LS]. In this case, (7.2) is replaced by ˆ c t =  I − u t u ∗ t δ + u ∗ t u t  − 1  W −∗ t − 1 ˆ c t − 1 − u t y ∗ t δ + u ∗ t u t  . (7.9) 8 APPENDIX: F AS T F A CTORIZA TION OF R − 1 t AR t T o achiev e a n O ( α 2 n ) factorization of R − 1 t AR t , we rewrite (4.6)  R − 1 t AR t   R − 1 t AR t  ∗ = δ I − α X k =1 s k y k ,t y ∗ k ,t ! , (8.1) where y k ,t is the k th column of Y t . W e define a sequence of α UT factor izat io n of the partial sums: ˜ A ( k ) t ˜ A ( k ) ∗ t = ˜ A ( k − 1) t ˜ A ( k − 1) ∗ t − δ s k y k ,t y ∗ k ,t , (8.2) with ˜ A (0) t = δ 1 / 2 I . W e rewrite ˜ A ( k ) ˜ A ( k ) ∗ t in the pro duct form ˜ A ( k ) t = ˜ A ( k − 1) t F ( k ) t = F (1) t F (2) t . . . F ( k ) t . ( 8 .3) Here F ( k ) t is the UT fa cto r of F ( k ) t F ( k ) ∗ t = I − δ k ξ k ,t ξ ∗ k ,t , with ξ k ,t defined by ˜ A ( k − 1) t ξ k ,t ≡ y k ,t . Eac h F ( k ) t is TIB and its in v erse is computable in O (3 n ) op erations giv en ξ k ,t . Giv en F (1) t . . . F ( k − 1) t , then O ( k n ) op erations are required to compute ξ k ,t . Th us O ( α 2 n ) op erations are required to complete the factorization. G iv en the factorizatio n, matrix v ector pro ducts are computable in O ( αn ) op erations. A t eac h step in this factorizatio n of ˜ A ( k ) t , w e gua r a n tee that F ( k ) t F ( k ) ∗ t is nonsingular b y reordering the ch o ice of columns in (8.2) suc h that all of the negativ e s k are pro cessed prior to the pro cessing of the p ositiv e s k . 11 Both ˜ A ( α ) t and R − 1 t AR t are triangular and satisfy (8.1). Since R t is also t riangular, the similarity transformatio n preserv es the order of the eigenv alues of A . Th us w e reco ver R − 1 t AR t b y righ t m ultiplying by the appropriate diagonal unitary matrix D to preserv e the complex phases of the diagonal of A . A CKN O WLEDGEME N T The authors thank Prof. A. Say ed a nd t he referees for their helpful commen ts. The w ork of KSR w as funded b y U.S. Dept. of Energy Grant D E-F G 0 2-86ER532 23. References [GGMS] P . E. Gill, G. H. Golub, W. Murray and M. A. Sa unders, “Metho ds for mo difying matrix factorizations,” Math. Comp . , vol. 28, pp. 505-535 , 1974. [GL] G. H. Golub and C. F. v an Loan, Matrix Computations , second edition, Section 12.6.4, Baltimore: Johns Hopkins Univ erssit y Press, 1 9 89. [GE] M. Gu and S. C. Eisenstat, “A Stable and Efficien t Algorithm for the rank-one mo dification of the symmetric eigen v alue problem,” SIAM J. Matrix Anal. Appl. , v ol. 15, pp. 1266-127 6, 1 994. [HD] E. J. Hannan and M. D eistler, The Statistic al The ory of Line ar Systems , p. 16, New Y ork: John Wiley and Sons, 19 88. [KS] T. Kailath and A. H. Say ed, “Displacemen t structure: theory and applications,” SIAM R eview , v ol. 37, pp. 297-3 86 1995 . [L T] P . La ncaster and M. Tismenetsky , T h e ory of Matric es , second edition, Boston : Academic Press, 19 85. [LK] H. Lev-Ari and T. Kailath, “Sta t e space approach t o factorization of lossless transfer functions and structured matrices,” Line ar Algebr a and its Applic ations , v ol. 162,pp. 273-29 5, 1992. [LS] L. Lj ung a nd T. S¨ oderstr¨ om, The ory and Pr actic e of R e cursive Iden tific ation , Cam bridge, MA: MIT Press, 1983. 12 [MR] A. P . Mullhaupt and K. S. Riedel, “Band matrix represen tation of triangular input balanced ,” submitted. [NSU] A. F. Nikiforo v, S. K. Suslo v and V. B. Uv aro v, Classic al O rtho gonal Polynomials of a Discr ete V ariable , pp. 30, 44, New Y ork: Springer V erlag, 1 9 91. [SK1] A. H. Say ed a nd T. Kailath, “ Extended Chandrasekhar recursions,” IEEE T r a n s. A utomatic Contr ol , v ol. 39, pp. 619-6 23, 1994 . [SK2] A. H. Sa yed and T. Kailath, “A state-space approach t o adaptiv e RLS filtering,” IEEE Signa l Pr o c essing Magazine , v ol. 11 , pp. 18- 60, 1994. 13