Banded Matrix Fraction Representation of Triangular Input Normal Pairs

Banded Matrix F raction Represen tation of T riangular Input Normal P airs Andrew P . Mullhaupt ∗ and Kurt S. Riedel † ∗ September 25, 201 8 Abstract An input pair ( A, B ) is triangular input normal if and only if A is triangular and AA ∗ + B B ∗ = I n , where I n is the i dentit y matrix. Input normal pairs generate an orthonormal basis for the impulse resp onse. Ev ery inpu t pair ma y b e tr ansformed to a triangular in put n ormal p air. A new system representati on is giv en: ( A, B ) is triangular normal and A is a matrix fraction, A = M − 1 N , where M and N are triangular matrices of lo w b andwidth. F or single inp ut pairs, M and N are bidiagonal and an explicit parameterization is giv en in terms of th e eigen v al ues of A . This b and fraction structur e allo ws f or fast up d ates of state sp ace systems and fast system iden tification. When A has only real eigen v alues, one state adv ance requires 3 n multi plications f or the single inpu t case. EDICS Num b ers: 2-SREP , 2-SYSM Key w ords. System represen tations, state space, balanced systems , orthonormal represen tations, system iden tification 1 INTR OD UCTION. An increasingly p opular tec hnique in sys tem identification is the use of orthonor- mal forms with fixed denominators [12, 19, 2, 13, 17, 15, 14]. In this approach, the ∗ S.A.C. Capital Manag emen t, LLC, 540 Madison Av e, New Y o r k, NY 1002 2, † Present Address: Millennium Partners, 66 6 Fifth Ave, New Y ork 10103 1 unkno wn impulse resp onse, G ( q ; { c k } ), is represen ted as a linear sum of orthonormal basis functions: G ( q , { c k } ) = n X k =1 c k B k ( q ) , (1.1) where q is the forw ard shift op erator . Here B k ( q ) are prescrib ed rational transfer functions with p o les { λ j , j = 1 . . . n } . The unkno wn co efficien ts are then iden tified using least squares identification or a robust ana lo g [6]. The utility of this approach is no w widely do cumen ted [14, 15]. In system iden tification, the use of orthonormal pa ir s impro v es the condition of the estimate [14, 1 1]. W e sp ecialize to the case where B k ( q ) corresp onds to the k th row in the semi-infinite impulse resp onse matrix, Ω( A, B ): Ω( A, B ) ≡ ( B , AB , A 2 B . . . ) . (1.2) Here A is the n × n state transition mat rix and B is a n × d matrix with n ≥ d . Let C b e t he p × n matrix whose k th column is c k . The j t h term or lead of the impulse resp onse in (1.1) is C A j − 1 B . In this pap er, we dev elop a computationa lly con v enien t r epresen t ation of the state transition matrices that arise fr om orthonormal bases. W e then sho w ho w these represen tations allow the rapid iden tification of or thonormal single input mo dels with prescribed p oles. O ur system represen tation results are general and we ll-suited for filter design as well. In the single input case, o ur filters are a pa r t icular realization of the orthogonal filters of [2, 13]. In system represen tation theory , one seeks families of system triples, ( A, B , C ), that parameterize the set o f all transfer functions sub ject to the similarit y equiv - alence. (Tw o sys tems are equiv alent if there exists an in v ertible matrix, T , suc h that ( A ′ , B ′ , C ′ ) = ( T AT − 1 , T B , C T − 1 )). W e examine r epresen tations of the input pair, ( A, B ). The input pair is in triangular input normal (TIN) fo rm if and only if AA ∗ + B B ∗ = I n , with A lower triangular (L T) and I n is the identit y matr ix. The input normal conditio n is equiv a lent to the row orthogonality of the semi-infinite impulse respo nse matrix, Ω. The state transition ma t rix A is represen ted as a matrix fraction: A = M − 1 N . Our results are not related to the extensiv e literature on represen tations of the transfer function G ( q ) as a matrix fraction: G ( q ) = N ( q ) /D ( q ). W e show tha t generic linear state space systems with forcing ǫ t ha v e an equiv alent represen tation of the form: M z t +1 = N z t + ˆ B ǫ t , (1.3) where M and N are lo we r triangular matrices with ba ndwidth d . In addition, M and N are chos en suc h that the c ovarianc e matrix of the state ve ctor z t tends to the identity ma trix f or white noise forcing. 2 W e b eliev e that our w ork is the first systematic study of sys tem represen tations, where A is a matrix fraction of banded matrices, A = M − 1 N . The ba nd fraction, A = M − 1 N , is neve r computed in practice. T o a dv ance t he state v ector, w e first compute r t ≡ N z t + ˆ B ǫ t and then solv e t he ba nded matrix equation M z t +1 = r t . Constructing r t requires n ( d + 1) m ultiplications and solving for z t +1 requires nd m ultiplications since M i,i = 1. Thus the t otal m ultiplication count for a state v ector up date is (2 d + 1) n . F or SISO sys tems with real eigen v alues, the state tra nsition matr ix is identical to the cascade systems prop osed in [2]. Ho w ev er, the parameterization/realizations of the ( A, B ) pair differ. In Section 2, the existence and iden tifiabilit y of TIN pairs is examined. G eneri- cally , TIN pairs are unique up to reordering the eigen v alues and phase ro tations of the co ordinates: B j,k ← B j,k exp i ( θ j ), A i,j ← A i,j exp i ( θ i − θ j ). In Sections 3-4, w e pro v e that generic TIN pairs ha v e a band fraction represe nta- tion: A = M − 1 N , B = M − 1 ˜ B , where M and N are triangular matrices of bandwidth d . W e deriv e the eigenv ectors of A in the d = 1 case and pro v e the n umerical stability of the represen tat io n when the eigen v alues of A are order in increasing magnitude. In Section 5, we discuss the numerical details of the band f raction represen tation and other realizations o f IN filters. W e relate our band fr action represen tation in the d = 1 case, to the orthonormal basis functions of [13]. In Section 6, we sho w the adv a n tage of our represen tations in system iden tification. Notation: In the remainder of the article, we suppres s the word ‘lo w er”. (‘T rian- gular” means lo w er triangula r .) Here A is a n × n matrix with eigen v alues { λ i } By A i : j,k : m , w e denote the ( j − i + 1) × ( m − k + 1) subblo ck of A fro m row i to row j and from column k to column m . W e abbreviate A i : j, 1: n b y A i : j, : . The matrix A has upp er bandwidth d if A i,j = 0 when j > i + d . W e denote b y ( B | A ) the n × ( n + d ) matrix formed b y the concatenation o f B and A . ‘Stable” means | λ | < 1. The n × n identit y matrix is I n and e k is the unit v ector in the k th co ordinate. W e define the equiv alence ( A, B ) ≈ ( A 2 , B 2 ) when ( A 2 ≡ E AE − 1 , B 2 ≡ E B ), where E n is a diagonal unitar y matrix: E j,k = exp( iθ j ) δ j − k . 2 TRIANGULAR INPUT NORMAL P AIRS In this section, w e presen t the fundamental represen tation results for the r eduction to triangular normal form. Definition 2.1 A n input p a i r, ( A, B ) , i s input no rm al (IN) if and only if I n − AA ∗ = B B ∗ . (2.1) The inp ut p a i r is a TIN p a ir if it is input normal and A is lower triang ular. 3 IN pairs are not required to b e stable o r con trolla ble. ( F rom (2.1), A must b e at least marginally stable.) In [8], ‘input balanced” has a more restrictiv e definition of (2.1) and the additional requiremen t that the observ abilit y Grammian b e diagonal. W e do not imp ose an y suc h conditio n on the observ a bilit y Grammian. W e choose this language so that ‘normal” denotes restrictions o n only one Grammian while ‘balanced” denotes sim ultaneous restrictions on b oth G rammians. Th us ( A, B ) is an IN pair if a nd only if the concatenated n × ( n + d ) matrix ( A | B ) is ro w orthogo nal. When ( A, B ) is input normal, then the identit y matrix solve s Stein’s equation (also kno wn as the discrete Ly apuno v equation): P − AP A ∗ = B B ∗ . (2.2) F o r stable A , P = P ∞ j =0 A j B B ∗ A j ∗ , So ( A, B ) is input norma l when Ω( A, B ) ∗ Ω( A, B ) = I n or equiv alen tly when the basis, {B k ( q ) } , constitutes a n orthonorma l set. Lemma 2.2 L et ( A, B ) and ( A ′ , B ′ ) = ( T AT − 1 , T B ) b e e quivalent IN p airs, with ( A, B ) s table and c ontr ol lable. Then T is unitary. Pr o of: Both I n and T − 1 T −∗ solv e the Stein equation. Since ( A, B ) is stable and con trollable, the solution of Stein’s equation is unique. Theorem 2.3 Every stable, c o ntr ol labl e input p air ( A, B ) , is similar to a lower tri- angular input norm al p air ( ˜ A ≡ T AT − 1 , ˜ B ≡ T B ) with k ˜ B k 2 ≤ 1 . The or der of the eigenvalues of ˜ A may b e sp e cifie d arbitr arily. If ( A, B ) is r e al and A ha s r e a l eigenvalues, then ( ˜ A, ˜ B ) a nd T may b e chosen to b e r e al. Pr o of: Let L b e the unique Cholesky low er triangular factor of P with p ositiv e diagonal en tries: P = LL ∗ . Here L is in v ertible. W e set ˆ A = L − 1 AL and ˆ B = L − 1 B . Note ( ˆ A, ˆ B ) satisfies (2.1 ) . By Sc h ur’s unitary triangularizat io n theorem (Horn and Johnson 2 .3.1), there exists a unitary matrix, Q suc h that Q ˆ AQ ∗ is low er triangular. The pro of of the real Sc hur fo r m is in Sec. 7.4.1 of [1] and Theorem 2.3.4 of [4]. The eigen v alues of ˜ A ma y b e placed in an y order [4]. Clearly , ˜ A = Q ˆ AQ ∗ and ˜ B = Q ˆ B satisfy (2.1) and therefore k ˜ B k 2 ≤ 1. W e now study the n umber o f equiv alent TIN pairs: Theorem 2.4 L et A an d ˆ A b e n × n L T matric es with A ii = ˆ A ii and le t A and ˆ A b e unitarily similar: ˆ A U = U A with U unitary. L et m b e the numb er of distinct eigen- values. Partition A , ˆ A an d U into m blo cks c orr esp onding to the r ep e ate d ei g envalue blo cks. L et n i b e the multiplicity of the i th eigenvalue, 1 ≤ i ≤ m . Then U has blo ck diagonal form: U = U 1 ⊕ U 2 ⊕ . . . ⊕ U m , wher e U i is a unitary m atrix of size n i × n i . 4 Pr o of: F rom ˆ A m,m U m, 1 = U m, 1 A 1 , 1 . If m > 1, t hen ˆ A m,m and A 1 , 1 ha v e no common eigen v alues. By Lemma 7.1.5 of [1], U m, 1 ≡ 0. Rep eating t his argumen t sho ws U m − k, 1 ≡ 0 f o r k = 0 , 1 . . . < m − 1. By orthogonality , U 1 ,j = 0 for 1 < j < m . W e con tin ue this c hain sho wing that U i, 2 = 0 for i 6 = 2, etc. Pro of by finite induction. Corollary 2.5 L et A b e an n × n matrix with distinct eigenva l ues. Then A is unitarily similar to triangular matrix ˆ A with or der e d eige nvalues and ˆ A is unique up to di a g onal unitary simi l a rities: ˆ A ← E ˆ AE ∗ , wher e E i,j = exp ( iθ j ) δ i − j . If A is inv ertible and ( A, B ) is TIN, then A is the Cholesky factor of I n − B B ∗ times a diagonal unitary matrix: Theorem 2.6 L et ( A, B ) b e a TIN p air with r ank ( I n − B B ∗ ) = n , then A = Cholesky ( I n − B B ∗ ) E n , (2.3) wher e E n is a diagon al unitary matrix. I n the r e al c ase, E n is a signatur e matrix. Pr o of: Let L b e t he Cholesky facto r of ( I n − B B ∗ ). Clearly A = LQ some unitary Q with L ii ≥ 0. By the uniquene ss o f the LQ factorization, L and Q are unique. Note Q = L − 1 A is lo w er triangula r and therefore diagonal. 3 BAND FRA C TION REPRESENT A TIONS OF SINGLE INPUT TIN P AIRS. F o r the single input case ( d = 1 ) , TIN pairs hav e an explicit band f r a ction rep- resen tation A = M − 1 N , B = ρ 1 M − 1 e 1 , where M and N are bidiagonal. The repre- sen tation is parameterized b y the eigenv alues of A . F or eac h eigen v alue, λ k , of A , we define ρ k = p 1 − | λ k | 2 , µ k = ρ k +1 /ρ k and γ k = λ ∗ k µ k . W e define the ma t r ices: M ≡ bidiag 1 1 1 · · · γ 1 γ 2 γ 3 · · · ! , (3.1) N ≡ bidiag λ 1 λ 2 λ 3 · · · µ 1 µ 2 µ 3 · · · ! , (3.2) where the top row contains the diagonal elemen ts o f M or N a nd the b ottom row con tains the ( n − 1) elemen ts of first sub diagonal of the bidiagonal matrices, M and N . Theorem 3.1 L et ˆ A ≡ M − 1 N an d ˆ B = M − 1 ρ 1 e 1 , wher e M and N ar e given by (3.1- 3.2 ) with | λ k | < 1 . Then ( ˆ A, ˆ B ) is TIN with eigenvalues { λ k } . 5 Pr o of: Explicit ev aluation sho ws M M ∗ − N N ∗ = ρ 2 1 e 1 e ∗ 1 . The Hautus criterio n states that if a stable matrix, A , is nonderogatory (There is only one Jordan blo c k for eac h eigenv alue of A in the eigendecompo sition of A .), then there is a v ector B suc h that ( A, B ) is con trolla ble. T hus the Hautus criterion together with Theorem 2.3 imply that Corollary 3.2 L et A b e a stable, nonder o gatory matrix. Then ther e exists a similar- ity tr ansfo rmation, T , such that ˆ A = T AT − 1 has the bidiagonal fr action r ep r esen ta- tion of The or em 3.1. The bidiagonal fraction can b e casted in a more aestetic for m A = M − 1 N = M − 1 0 N 0 , where M 0 and N 0 are defined as M 0 ≡ bidiag c 1 c 2 c 3 · · · s ∗ 1 s ∗ 2 s ∗ 3 · · · ! , (3.3) N 0 ≡ bidiag s 1 s 2 s 3 · · · c 1 c 2 c 3 · · · ! , (3.4) with c k = 1  q 1 − | λ k | 2 and s k ≡ λ k / p 1 − | λ k | 2 . In a series of excellen t articles [1 9, 2 , 13, 17, 15, 1 4], a nu mber of researc hers ha v e constructed IN filters using cascade realizations. W e now consider the case single input case where A has real eigen v alues. In this case, the cascade realization of the IN filter is triangular. By Corollary 2.5, there is a unique triangular state adv ance matrix with the given eigen v alue ordering ( up to diago na l similarit y transformations with ± 1 elemen ts). Thus a ll state space constructions of IN filters are simply different realizations of the same ( A, B ) pair (up to sign flips of the co ordinates). In Section 5, w e discuss the n umerics of the v ar io us realizatio ns of IN filters. W e now give a more constructiv e pro of of equiv alence of the band fraction repre- sen tation with the orthogonal basis f unctions o f [13]. Giv en a set of deca y rat es/p oles of the frequency resp onses, { λ n } , Ninness and Gustafsson deriv e a set of o r t honormal basis functions in the frequency domain: ˆ z n ( w ) = ∞ X t =0 z n ( t ) w − t = w α q 1 − | λ n | 2 w − λ n n − 1 Y k =1  1 − λ ∗ k w w − λ k  , (3.5) where α is 0 o r 1. F rom this, w e deriv e the relation: ( w − λ k ) ˆ z k ( w ) / q 1 − | λ k | 2 = (1 − λ ∗ k − 1 w ) ˆ z k − 1 ( w ) / q 1 − | λ k − 1 | 2 . (3.6) 6 for k ≥ 1 with λ 0 ≡ 0 and ˆ z 0 ( w ) ≡ z α . Equating lik e p ow ers of w yields z k ( t + 1) + λ ∗ k − 1 µ k − 1 z k − 1 ( t + 1) = λ k z k ( t ) + µ k − 1 z k − 1 ( t ) , (3.7) with z 0 ( t ) ≡ δ 0 , − α . This sho ws that the bidia g onal matrix fraction represen tation generates the same basis functions as in [13]. The k -th comp o nen t of z t is to b e iden tified with the t -th lead of the k -th orthonormal transfer function, B k ( q ). W e sho w that if the eigen v alues are sort ed in order of increasing magnitude, then M − 1 do es not ha ve la rge elemen ts. W e no w give tw o results concerning the bidiag onal matrix fraction represen tation ( d = 1 ) . The first result sho ws that if the eigenv alues are sorted in order of increasing ma g nitude, then the bidiagonal factor ization is w ell- conditioned. W e then giv e the eigen v ectors of A . A. Exp li cit M inversion. W e can in terpret the bidiagonal ma t rix fraction represen tation as an LR factor- ization of the augmen ted sys tem: ( B | A ) = M − 1 ( e 1 | N ). W e now explicitly in v ert M in (3.1) and show that the n umerical conditioning of the matrix inv ersion is go o d without piv oting. W e define the x − condition n um b er of an inv ertible matrix, G , as κ x ( G ) ≡ k G k x k G − 1 k x , where x = 1 , 2 , or ∞ [1, 3]. Here k G k 1 is the maximu m column sum norm and k G k ∞ is the maxim um ro w sum norm. Theorem 3.3 L et ( A, B ) b e a stable TIN p air ( d = 1 ). If the eigenvalue mag n itudes ar e in asc ending or der, | λ k +1 | ≥ | λ k | , then | ( M − 1 ) i,j | < 1 for i > j and κ 2 ( M − 1 ) ≤ 2 n . Pr o of: The LR factorizat io n is ( B | A ) = ( M − 1 ) ( e 1 | N ). Now M is unit lo w er t r ia n- gular with ( M ) k +1 ,k = λ ∗ k . So ( M − 1 ) i,j = ( − 1) i − j Y j ≤ k j . Note k M − 1 k x ≤ n and k M k x ≤ 2 for x = 1 and x = ∞ . Thus κ 2 ( M − 1 ) ≤ κ 1 ( M − 1 ) κ ∞ ( M − 1 ) ≤ 2 n . This implies that the LR factorizat io n of ( B | A ) is n umerically stable without piv oting. B. Eigenve c tors W e now ev a luate the eigen v ectors of the single input TIN system using the pa- rameterization (3 .1)-(3.2) for the case of distinct eigen v alues. It is w ell kno wn that the eigen v ectors of triangular ma t r ices satisfy a recursion f orm ula. Let V ij b e the i -th comp onen t of the j -th eigen v ector with V ij = 0 for i < j . W rite V = [ V ij ] and Λ ≡ diag ( λ 1 , . . . λ n ). The eigen v ector equation AV = V Λ b ecomes N V = M V Λ or elemen t b y elemen t: λ k  λ ∗ j − 1 µ j − 1 V j − 1 ,k + V j k  = ( µ j − 1 V j − 1 ,k + λ j V j k ) , (3.9) 7 W e set V k k = 1, V j k = 0 for j < k and solv e the recursion for V j k when j > k : V j k =  1 − | λ j | 2  1 / 2  1 − | λ k | 2  1 / 2 λ k − λ j " Y k k . Pr o of: By Theorem 3.2.1 of [1], ( B | A ) = L − 1 R has a unique ( B | A 1: n, 1:( n − d ) ) = L − 1 R has a unique represen tation where R is UT and L is L T with L i,i = 1. Let ˜ R b e the submatrix of R con taining columns ( d + 1) through ( n + d ). Since A is L T, ˜ R = LA is L T of bandwidth d . Note R ( I d ⊕ D ) R ∗ = LD L ∗ . By Theorem 4.3.1 of [1], L has bandwidth d . W e set M = L and N = ˜ R F o r d = 1, ev ery con trollable TIN pair has no nv anishing principal minors. The condition that ( B | A ) ha ve nonv anishing principal minors for k < n is generically true. If ( B | A ) has an LR decomp osition and a v anishing principal minor for k < n , then the induced TIN pair, ( A, B ), do es not ha v e a unique LR decomp osition of ( B | A ) with M i,i = 1 . F r o m M M ∗ = N N ∗ + ˆ B ˆ B ∗ , the singular v alues of M are singular v alues of ( ˆ B | N ). When the condition that M i,i = 1 is relaxed, other band fra ction represen tations may b e generated by ( ˆ B , N ) ← D ′ ( ˆ B , N ), M ← D ′ M , where D ′ is a nonsingular diagonal matrix. 8 The band fractio n structure implies tha t the n × ( n + d ) matrix Y ≡ ( ˆ B | N ) is upp er triangular with upp er ba ndwidth d . W e now examine parameterizations of the band fraction represen t a tion. W e define Definition 4.2 L et Y denote the set of al l upp er triangular n × ( n + d ) matric es with upp er b andwidth d . L et Y 1 b e the set of Y ∈ Y such that Y i,i = 1 for i < n and k Y n, : k = 1 with a p ositive first nonz er o element in the n th r ow. Note Y has n ( d + 1) no nzero co ordinat es while Y 1 is an nd dimensional manifold. W e no w show that w e can construct a TIN pair from any Y in Y using the LQ decomp osition. Theorem 4.3 L et Y ∈ Y h ave r ank ( n ) and define M and ˆ Q as the LQ d e c omp o sition of Y . D efine B = ˆ Q : , 1: d and A = ˆ Q : , ( d +1):( n + d ) . Then M has lower b and width d and is invertible and ( A, B ) is a TIN p a i r. Pr o of: Note rank( M ) =rank( Y ) = n . Let ( ˆ B | N ) = Y . Note M M ∗ = Y Y ∗ = N N ∗ + ˆ B ˆ B ∗ . By Theorem 4.3.1 of [1], M has bandwidth d . F rom A = M − 1 N , A is triangular and AA ∗ + B B ∗ = I n implies ( A, B ) is a TIN pair. W e no w examine the map f ( Y ) = ( A = ˆ Q : , ( d +1):( n + d ) , B = ˆ Q : , 1: d ), where M ii ( Y ) > 0. F rom Theorem 9.1 of [3], ( B | A ) = M − 1 Y has nonzero principal minors f or k < n if and only if Y do es. W e can r escale Y and M by a nonsingular diagonal matrix D : Y ← D Y a nd M ← D M , and preserv e the induced TIN pair , ( A, B ). This f r eedom and complex phase equiv alence allows us to restrict our consideration to Y ∈ Y 1 : Theorem 4.4 F or e a ch TIN p air, ( A, B ) , with nonvanishi n g princ i p al minors in ( B | A ) for k < n , ther e e x i s ts a unique Y ∈ Y 1 and a unique diagonal unitary matrix E such that ( E AE ∗ , E B ) is ge n er ate d by Y : f ( Y ) = ( E AE ∗ , E B ) . Pr o of: Let ( ˜ B , ˜ N , ˜ M ) be the unique band fra ction decomp o sition of ( B | A ) with ˜ M ii = 1. F or an arbitrary diagona l unitary mat rix E , the set of band fraction repre- sen tations of ( E AE ∗ , E B ) is { ( D E ˜ B , D E ˜ N E ∗ , D E ˜ M E ∗ ) |} , where D is nonsingular and diagonal. Supp ose b ot h Y 1 and Y 2 ∈ Y 1 generate an equiv alen t v ersion of ( A, B ), i.e. f ( Y i ) = ( E i AE ∗ i , E i B ). Th us Y i = ( D i E i ˜ B | D i E i ˜ N E ∗ i ) for some D i and E i where D i is nonsingular and diagonal and E i is diagonal and unita ry . Since M i = D i E i ˜ M E ∗ i , the condition M i,i > 0 sho ws tha t D 1 and D 2 are p ositiv e matrices. The requireme nt that Y i ∈ Y 1 forces | D 1 | = | D 2 | and E 1 = E 2 . W e can parameterize input normal systems using Y ∈ Y 1 for d > 1. There are t w o difficulties with t his par ameterization. First, the represen tations ( M − 1 N , M − 1 ˆ B ) f or 9 ( ˆ B | N ) in Y 1 con tain man y equiv alen t TIN represen tations corresp onding to differen t orderings of the eigenv alues of A . Second, in the Y 1 parameterization, the eigen v alues of A are only know n after one solve s M b y computing the LQ decomp osition of Y . 5 NUMERICS W e now discuss the computational sp eed o f v arious system represe ntations. The bidiagonal band fraction represen tat ion, A = M − 1 N , requires only 3 n m ul- tiplications for a state v ector up date in a single input system. First, w e compute r t ≡ N z t + e 1 ǫ t in 2 n m ultiplications. W e then solv e the bidiagonal matrix equation M z t +1 = r t in n mu ltiplications since M i,i = 1. m ultiplications. The matrix op era- tions in the band fraction adv ance may be implemen ted using the dtbmv and dtbtrs routines from the optimized LAP AC K soft w are [5]. W e caution that if A has nonreal eigen v alues, then M and N are complex. Thus the computational adv antage of the bidiagonal matrix fraction r epresen t ation is primarily limited to the case when the eigen v alues of A are real. In [2], Heub erger et al. prop ose a cascade realizatio n of IN filters. This cascade represen tation r equires 4 n m ultiplications for a state v ector up date in a single input system. A differen t cascade realization o f SISO IN pairs is giv en in Figure 1 of [13]. The Ninness and Gustafson represen tation uses a lossless cascade follo w ed b y a set of AR ( 1) mo dels . This realization requires at least 5 n m ultiplications for a state up date. The direct f orm represen tations [12] hav e ev en faster state adv ances, but these represen tations a re not IN a nd can hav e v ery ill-conditio ned Grammians [11, 14]. W e advise aga inst using high order mo dels that are not in input normal form. The tridiagonal form [7] is another fast represe ntation that is not input normal. In [12, 18], em b edded lossless represen ta tions are constructed. F or d = 1, these em b edded lossless filters 8 n m ultiplications p er adv ance [12]. These embedded filters include the cost of ev aluating C · z t while the IN filter represen tation has an additional cost of n multiplications to compute C · z t , so we should credit the em b edded filters with n multiplic ations, yielding 7 n m ultiplications. Unfortunately , embedded lossles s systems are not input normal. In [10], we giv e a differe nt represen tation of Hes sen b erg and tria ngular input pairs b y treating ( B , A ) as a pro jection of a pro duct o f nd Giv ens rotatio ns. The m ultiplication coun t for these Giv ens pro duct represen tations is 4 dn , whic h for large d is t wice as larg e as t he band fraction represen tations. If Householder transformations are used , the multiplication coun t is again 2 nd asymptotically . The G iv ens pro duct 10 represen tations are a MIMO generalization of the cascade architec ture. This approac h ma y b e more natural for multiv ariate case. In the d = 1 case, the band fraction represen tation has the adv an tage that the eigenv alues of A are para meters of the system. T his eigen v a lue parameterization is con v enien t when the eigen v alues are presp ecified or adaptiv ely estimated. The Giv ens pro duct represen tations in [1 0 ] ma y b e less sensitiv e to r o undoff error since they are o rthogonal matrices [12]. The computational adv an tage of the band fraction represen ta tion is that fast n umerical algorithms for band ma t r ix inv ersion and m ultiplication ar e readily av ailable. 6 F AST SYSTEM IDENTIFICA TION The band fraction represen tations may b e used for the rapid iden tification of impulse resp onses. W e remain in the framew ork where a fixed basis of ort ho gonal basis functions are given as describ ed in [1 2 , 19, 13, 17, 15, 14] and summarized in the in tro duction. W e then use the band fraction repres entation of Sections 3- 4, so ( A, B ) ≡ M − 1 ( N , ˆ B ) is giv en. The state v ector ev o lv es according to (1.3). W e effectiv ely compute the sec ond momen t matrices, ˆ P δ t ≡ P t i =1 δ t − i z i z ∗ i and ˆ d δ t = P t i = i δ t − 1 z i y ∗ i , where δ is the for g etting fa ctor. The unkno wn co efficien ts, ˆ C , using recursiv e least squares (RLS). A t eac h time step, w e reestimate ˆ C by solving ˆ P δ t ˆ C t = ˆ d δ t . This is the normal equations for the least squares estimate of ˆ C t . The k th comp onen t of z t represen ts the k th orthonormal transfer function, B k ( q ), applied to the sequence { u 1 , u 2 , . . . u t } . The resulting estimate of the the j th lead of the impulse resp onse is ˆ C A j − 1 B . T o solve for ˆ C t ., w e recomme nd using a QR up date o f the square ro ot of normal equations [16]. Ev en faster metho ds are a v a ilable that use the displacemen t ra nk structure of ˆ P δ t [9, 16]. The in terested reader can consult [6, 14] for a comprehen siv e description of adaptiv e estimation. The input normal filter represen tations ar e adv an tageous fo r man y reasons: First, ˆ P δ t t →∞ − → constant × I n , when the ‘true” state space mo del is used and the forcing noise is white. Th us the regression for ˆ C t is w ell-conditioned. Similarly , IN filter structures are resistan t to r o undoff error [1 2]. Th us IN represen tations will t ime asymptotically satisfy the ansatz need b y least mean squares (LMS) iden tification alg orithms . This leads to a second adv an tage of IN filters: Gradient algorit hms suc h as the least mean squares (LMS) a lg orithm o f ten p erfo r m w ell enough in certain applications to obviate the need for more complicated and computationally in tensiv e RLS algorithms. Third, an O ( n ) up date o f ˆ C t is p ossible [9]. Finally , when the a dv ance matrix is triangular, 11 the IN pairs form nested families: ( A 1: k , 1: k , B 1: k , : ) is TIN for k ≤ n when ( A, B ) is TIN. This nesting ma y b e used in adaptiv e order selection. These a r e adv antages of all TIN filter represen tations. F o r single input pairs, The adv antages of the band fraction represen tation o v er the cascade represen tation a r e primarily n umerical, the faster computational sp eed (3 n m ultiplications v ersus 4 n or more) and the a v ailability of f ast linear algebra soft w are for banded matrices. 7 SUMMAR Y W e ha v e sho wn tha t TIN pairs generically hav e a matrix f raction represen tation, A = M − 1 N , where M and N ha v e lo w er bandwidth d . This structure allo ws fast state space upda t es, M z t +1 ← N z t + ˆ B ǫ t , and fast solution of matr ix equations: Ax = f . The t o tal op erations coun t is (2 d + 1) n multiplic atio ns. W e b eliev e that our w ork is the first study of system represen tations, where A is a matrix fraction of banded matrices, A = M − 1 N . F o r single input filter design, our represe ntations are parameterized explicitly in terms of the eigen v alues, { λ k } . When M is bidiagonal, solving M z t +1 = f t requires only O ( n ) m ultiplications and may b e implemen ted as a systolic arr a y with or der indep enden t lat ency . As sho w ed in Theorem 3.3, if the eigen v a lues are ordered in ascending magnitude, then the matrix in v ersion is w ell-conditioned and suitable for fixed p oint op erations. F or syste m iden tification, the co v ariance of z t tends to the iden tit y matr ix when the forcing noise is white. Th us least squares r egression is b oth w ell-conditioned and computationally fast. Since A is t riangular, one can compute the ev olution of all em b edded mo dels of dimension n ′ < n simply b y pro jecting t he mo del of dimension n onto the first n ′ co ordinates. Th us the triang ula r structure simplifies the ev aluation o f this nested family of mo dels for mo del order selection. References [1] G. H. Go lub and C. F. V an Loan, Matrix C o m putations , third edition, Baltimore: John Hopkins Univ ersit y Press, 199 6 . [2] P .S.C. Heub erger, P .M.J. V an den Hof, and O.H. Bosgra, ‘A generalized or- thonormal basis for linear dynamical systems,” IEEE T r ans. Aut. Cont. , v ol. 40, pp. 451-465, 1995. [3] N. J. Higham, A c cur acy an d Stability of Numeric al Algorithms , Philadelphia: SIAM Press , 1996. 12 [4] R. A. Horn and C. R. Johnson, Matrix analysis , Cam bridge: Cam bridge Univ er- sit y Press, 1985. [5] E. Ande rson, et a l. LAP ACK User’s Guide Philadelphia: SIAM Press, 1999 . [6] L. Ljung and T. S¨ oderstr¨ om, The o ry and pr actic e of r e cursive identific ation , Cam bridge, MA: The MIT Press, 1983. [7] T. McKelv ey and A. Helmersson, ‘State-space parameterizations of multiv ari- ate linear systems using tridiagonal matrix f orms,” Pr o c. 3 5th IEEE Conf. on De cision and Co ntr ol , pp. 1666- 1 671, IEEE Press, 199 6. [8] B. Mo ore, ‘Principal comp onen ts analysis in linear systems: controllabilit y , ob- serv ability , mo del reduction,” IEEE T r ans. Aut. Cont. , v ol. 26, pp. 17-32 . [9] A. Mullhaupt and K. S. Riedel, ‘F ast iden tification of innov ations filters,” IEEE T r ans . Signal Pr o c essing , v ol. 45, pp. 2616-26 1 9, 1997. [10] A. Mullhaupt a nd K. S. Riedel, ‘Hessen b erg and Sc h ur output normal pair rep- resen tations,” submitted. [11] A. Mullhaupt and K. S. Riedel, ‘Bounds on the condition n umber of solutions of the Stein equation,” submitted. [12] R. A. Roberts and C. T. Mullis, Digital Signal Pr o c essing , Addison W esley , Reading, MA, 1 987. [13] B. Ninness and F. G ustafsson, ‘A unifying construction of orthonormal bases for system iden tification,” IEEE T r an s . A ut. Cont. vol. 42, pp. 51 5 -521, 1997. [14] B. Ninness, ‘The utility of orthonormal bases in system iden tification,” T ec hnical Rep ort 9802, Dept. of EEC E, Univ ersity of New castle, Australia, 19 98. [15] B. Ninness, H. Hjalmarsson and F. G ustafsson, ‘The fundamen tal role of or- thonormal bases in system iden tification,” IEEE T r ans. A ut. C ont. vol. 44, pp. 1384-14 07, 1999. [16] A. H. Say ed and T. Kailath, ‘A state-space approach to adaptiv e RLS filtering,” IEEE Signa l Pr o c essing Magazine , v ol. 11, pp. 1 8-60, 1994. [17] P .M.J. V an den Hof, P .S.C. Heub erger, and J. Bok or, ‘System identific atio n with generalized o rthonormal basis functions,” A utomatic a , vol. 31, pp. 1821-1 831, 1995. 13 [18] A.J. v an der V een and M. Vib erg, ‘Minimal con tin uous state-space par ametriza- tions,” In Pr o c. Eusip c o , T rieste (Italy), pp.523-526 , 1996. [19] B. W ahlb erg, ‘System identification using Laguerre models,” IEEE T r ans. A ut. Cont. , v ol. 36, pp. 551-562, 1991. 14