On Behaviors, Rational Interpolation, and the Lanczos Algorithm

On Behaviors, Rational Interpolation, and the Lanczos Algorithm

3a-OI 5 Copyright © 1996IFAC 13th Triennial World Congres~. San Francisco. USA On Behaviors, Rational interpolation, and the Lanczos algorithm" A.C...

357KB Sizes 6 Downloads 46 Views

3a-OI 5

Copyright © 1996IFAC 13th Triennial World Congres~. San Francisco. USA

On Behaviors, Rational interpolation, and the Lanczos algorithm" A.C. Antoulas Department of Electrical and Computer Engineering, Rice Uuiversity Houston, TX 77251-1892, USA E.J. Grimme Coordinated Seieuce Laboratory, Universit.y of Illinois at. Urbana-Champaign Urbana, IL 61801-2307, USA D.C. Sorensen Department. of Computational and Applied Mat.hematics, Rice University HOllston, TX 77251-1892, USA

Abstract. This paper explores the interconnections between two methodl'l which can be used to obtain ra.tional interpolants. The first method, the bchavioral approach, construc:t~ a generating system in the frequency domain which explains a given data set composed of trajectories. The second method, the mh()nal Lanczo.~ algorithm, can be used to c-onstruct a rational interpolant for the transfer function of a linear system defined by (potentially very high-order) statC'-space equations. This paper works to merge the t.heoretical attributes of the behavioral approach with the t.heoretical and computational properties of rational Lanczos. As a result, it lays the foundation for the computation of reduced-order, stabili7.ing controllers through ratiollal inlerpolation. Keywords: linear systems~ hchaviors, rational interpolation, Lanczos algorithm, Krylov sllbspaces.

1. Introduction This paper seeks to implement the mct.hods of [1] and [2] in a numerically robust and efficient way by providing a connection with the Lanczos algorithm. The ability to c.onstruct accnrate rational interpolants and associated controllers regardless of the problem's original size is of great importancp. It is for large systems that *This work was conducted while the s('c,.md author was visiting the Depart.ment or Computational and Applied .\1athematics at Rice University, January-March 199:",.

both rational interpolation and reduced-order models become mo~t competitive compared la more traditional approaches.

The paper reviE-~ws two approaGhe~ for computing rational interpolants in Section 2. The first results from the behavioral approach [1, 2], and constructs a unique (up to equivalence) autonomous system of minimal complexity called the most powe1ful or generating system, ~(e*), which explains the given data set D composed of trajectories. This data set D is equivalent to the rational interpolation data Het (the latter being composed of alTays of complex numbers). The main tool obtained 1 is the square non-singular polynomial matrix 8'(s). The second approach, rational Lanezos [4 , 5], utilize:-; a biorthog:onal projector IT composed of so-called rational Krylov subspaccs [8]. Rational Lanczos is computationally efficient since it never explicitly computes the moments to be matched Connections between the two method are developed in Section 3. It is shown namely, tha.t the results obtained by the rational Lanczos algorithm can be used to construct the polynomial matrix 8* (8). By merging the two approaches, one gains the theoretical properties of the behavioral approarn and the computational attributes of rational Lanczos. In the la.<:,t. Section wc illustrate the interconnections between the behavioral and rational Lanczos approacheH by means of an example.


2. Computational Approaches to Rational Interpolation

is composed of time series taking on the special form [2J


Given an order-n single-input 1 singi('-output linear system described by means of the geupralized state space equations


= Ax(t) + butt),


= eT x(t)

Z(t,CT;,];) = eXPu,lt) LZj(t,CTi,]i) where



= Ai(t) + butt),


= el'i(t) + duet)


Ji = 0, 1, ... Ji -

mj. (cr;)


1, and £ = 1,2, ... ,L The t.erms

_eT {(A - cr;E)--lE


·r"_·_ )I' ,J1 J 1.

and Zj = 0 for j 2: li' Each of the time series in D is associated with a specific interpolation point ai_ The

(A - IT;E)-lb


arc respectively the jI moments of the original and reduced-order systems at the frequency ai E C. The number of moments matched, E!= I jil is at least 2k. If the rational flUlctions g(s) and yes) arc the transfer functions corresponding to (1) and (2) respectively, then m;.(o-;) = g(;;)(o-;) andmj.(IT;) = 911 ')(<7;) wheregU)(s) is the j'h derivative of g(s). Th1ls fI(") interpolates g(s). In thi8 sectioIl 1 two approachps for computing a ratioIlal interpolant are reviewed. \Vc will only concentrate on those results of importa.nce in the sequel. The reader is referred to the original papers for greater detail. 1. Behavioral approach

In the behavioral approach [1, 21, dynamical systems are determined which explain, i.e., could have generated, a collection of time series Z(t, i) contained in a certain data set D. For our problem, Z(t, i) is equivalent to [,,;(t) y;(t)V where ,,;(tl, y;(t) is a specific input-output pair satisfying (1). A fiystem is said to be an llnfal.~ified model of the data ~et D, if D C B where the linear subspace B is the behavior of the system, that is the collection of all traje(:tories which satisfy the system equations (1). The system defined by the smallest linear suus pace B'" covering D is said to be the must powerful unfalsifted model (MPUM) of D. Rational interpolation turns out to be a special case of the general modeling problem. The data set

= {Z(t,"'l,JtI, ... ,7,(t,IT,,],)}




g(j)(,,-il j!

+ 1) =








quantity), is the number of pieces of data (moments) to be matched at the it" interpolation point. Note that the terms Zj making up Z(t, CTil.ii) depend on the number of data pieces being matched. One can easily verify the recursions

which satisfies


Zj(t,CTi,li)= (

tJi -1

1») (); -



the rational interpolation problem can be phrased as finding an order-A; < n model descrihed by E:i:(t)



and Zj(t,<7;,j) = 0 for j > };. It is also straightforward to verify that the inverse reclIrsion


= (:t

- "';)Zj(t,IT;,j

+ 1)


holds for any nonnegative j. Corresponding to any behavior B generated as the span of polynom'ial - €:J.:punential trajectories as above, one can atlsociate a differential operator with constant coefficients, e(d/dt), such that H = ker(e), i.e., e(d/dt)Z = 0 for all Z E l3. The system whose behavior is B can be described by the differential equations corresponding to e(d/dt) and is denoted E(e). Corresponding to B* is the so-called generating matrix e'(d/dt) which, for all D c H, satisfies

Dc B' = ker(e') C ker(EJ) = H Some of the salient features associated with 8'" are summarized below. See [2] for details. Theorem 1 Given the data. set D and an associated generating matrix 8"', the following hold true: (a) Any system which model .• D takes the form E(fEJ') whe1·e r is an arbitrary polynomial matrix. If r is unimodular' then E(rEJ') = E(W).

(b) The sy,.tern E(EJ') has minimal complexity among all systems modeling D. This minimal complexity is deg det W(s).



(e) The roots of det 0'(8) form the ,
(d) The generating system of the data ,.et D U D+ can be written as 0+(s)0'(8) where 0+(,,) is an appropriately comput~d polynomial matrix. Theorem l(a) justifies the name of generating matrix for e.... Theorem led) .suggcst.~ a recursive approach to generating a minimal complexity model for D. First, construct f:li,] such that 0i,] (d/dt)Z(t, <71,1) = O. An update 8 2 is then eomput~~d so that 8 21 := 8 28i" 1 annihilates Z(t, al, 2) as well; in other ~ords 0 21 (d/dt)Z(t,<71,2) = 0, One then proceeds to repeatedly update 0;_],1 through left IIlultiplication by 0; until a e;k I is found which annihilat.es all the trajecto-des of the 'data in (3). An explicit construction of and examples are provided in [1,2], By simple, repetitive updates, one can construct the 2 x 2 polynomial matrix e~.{"l (8) so that E(e"') solves the rational interpolation prohlern. One ~hooses this generating matrix so that 0;.,1 (d/dtJ[lli(t) Yi(t)]T = 0 for appropriat.e Ui(t) and Yi(t) corresponding to D. Then by Theorem I (a),



0]0; •. 1(8) =:



yields a strictly proper transfer function ~ which is a <> rational interpolant of g( s), If 2k+ 1 pieces of data are to be matched, d oF and hence the product [1 0]0;k+1.1 leads to only a proper rational intcl'polant. .

where tJ! is a matrix and x is a column vector. Ji turns out to be the number of moments matched by the kth order rational interpolant about the interpolation frequency O"i. If Vk and Wk satisfy the above conditions, it can be shown (proof by construction; see [4] and Algorithm 1 below) that the relations

As is well knowll, llumerical difficulties arise by explicitly computing the moments in t he rational interpolation problem. The Lancl,Os algorithm avoids these problems by comput.ing t.he intcrpolant.s without computing the moments [7], The rational Lanczos method [4. 5] generates a rational intel'polant by applying a specific biorthogonal projector rh = VdV[ to the original system defined by (1). Specifically, Fk E Rnxk and ll'k E Rnxk are chosen to satisfy (i) the biorthogonality conditioIl, W[Vk = I, and (ii) the rational Krylov sllhsp"rce condition [8]: colsp(V.) =

UKit ((A -

a;E)-l E, (A - aiE)-lb)


colsp(Wk )

= UKit+1 (ET(A-<7iE)-T,C)



hold, where the matrices Hk,k and Kk,k are upperHessenbcl'g (and in fact, essentially tridiagonal [5]) matrices of size k x k. The structure and features of the rational Lanczos algorithm are explored in detail in [5]. At every iteration of the algorithm, another column is added to both the vectors F and VV. The interpolation point used in the kth iteration, it" is chosen from a set of I (user specified) points, A pseudo-optimal strategy is developed in [5J for choosing it; so as to most improve the ktn order model over the (k - 1)", This stalegy has the additional advantage that it tends to drive the algorithm away from breakdovms, i.e., divisions by zero. The relation (7) and the matrices and K.,k play a key role in the reduced-order model of rational Lanc7,os. The reduced-order mod(~1 in (2) is chosen as


A = K.,. + ai;H.,k, E = H.,h, cl:= d = b = W[C4 - ai;Et 1b, c = cTV.K.,k


2. Rational Lanczos

= EV'+1Kk+1,k

Hu = W[(A - <7i;E)-lVk K.,k

(A - ai; )V'+1 Hk+1,k

° (8)

The choices made in (8) were first motivated in [4], For our purposes, it is most important to note that this reduced-order model does indeed solve the rational interpolation problem, For a proof see [5], Theorem 2 Given the reduced-order model (8) generated by the mtional Lanczos algorithm,


= mj. (Ui), for i = 1, ... , O'~ and ji =

1, ... ,:li.

The value of Ji is equal to twice the number of times appears in the finite sequence {a- i" , U i" , ... , a- i* }.




The rat.ional Lanczos algorithm provides an alternative and numerically stable approach for computing a rational interpolant. It should be noted however that the rational Lanc7.os approach can only match an even number of moment.~ about any interpolation point due to to the fact that Ji must be even. In general, rational interpolation (and also t.he bchavioral approach) allows Ji to be an odd number. This difference turns out to Le of little fOIlsequence in general.


3, Relating the Rational Lanczos and Behavioral Approaches

where colsp denotes the column span. A Krylov subspace is defined as


Since both the behavioral and rational Lanczos approaches lead to a rational interpolant , it would seem that connections can he drawn betwet".Il them . Wc lay out these connect.ions in thi.~ section.

where for the tridiagonal matrix one typically uses the nota.tion {Jk, Ok and 'Yk in place of (Jh+ i ,k, /3k, k 1 /3k ,k+ 1· Denoting the characteristic polynomial of H". as p.(8) , it is known (see [6)) that thes" polynomials are generatccl from the three term recurrence

1. Single point case

Theorem 3 Let "' - '(') and ",(I'»~ be th e transfer fun eqJ.-J $ q" ill tions corresponding to the reduced-order models generated by k - 1 and k iteration,.; re.<;pectivdYI of the. rational Lanczos algorithm about a 8ingle interpolation point (J. Th e determinant of the matri:r \)!

() k.1

s =


From this expression it is easily seen that the reverse polynomials (10) satisfy the reverse three term recur-


In matrix form, this recursio n becomes


(s _ ".)'P. _ l (S )

is equal (up to a constant) to (s - aj2k.


Proof: We begin by recalling (8), the defini tion of the reduced-order model arising from k rational Lancws it-

H", i(t) = (J + uHu)x(t)

+ bu t t),

:,,(t) = (ri;(I) (9)

1- (A -"')"'j (

= >,~(I'

(A - 0" )'

so t hat

since f(k,k = I. It is well-knmvu that the poles of this reduced-order system are the generalized eigenvalues of the pencil (A Hk ,k - (l + aH.,d)· The characteristic polynom ial of the matrix H k,k is

p..(8) = det(D! - H,.,.)


Now, define Po 0 and P_,{A) (A - a)-'. These choices are consistent with the initialization of the seQuence of Lauczos polynomiab. l\1oreover) define

erations. For a single interpolat.ion point u 1 this model simplifies (,0

P utting ()

(9 - "'HI )1,.(11) - {3n,p'-1 (9)

PHI (8) =

First we examin e, for simplicit.y, lhe case where tbere is only onc interpolating point (1. The following result holds true.

det(>jij(A) = (Jj-Jlj- l(A - 0")2




(}]>jij(,\)) = (}] (Jj- I'YJ-I) (A_a)"

Finally) inLrocluce the 2 x 2 matrix poly nomial recursion

gives the n:vel'se polynomial


:= -(A _ a)'Pk

(A - a )kdet (Hk,k -

'MA) ) = (}11Ji k j (A) ) 1Ji 0 , 1

(_1_) = _I-I) =


A- a


where \)!0.1 is initialized wit.h t/Jo = Po = 0 and 0 = ,) = deg(p.) + 1 and

A-a det ( A - rr)Hk,k - I ) =

det (AHk ,k - (aHu

+ I))


T hus the poles of the kth order model (9) are also the roots of t he polynomial pdA). For the single interpolation point case, the matrix recur,ion (6) leading to (9), can he simplified to

(.4 - "E) - I),.

= V.Hu + Jk+ I" T,e[

But this expression is a standard Lanczos recurrence so that H k,' is a partial t ridiagona lization of (A -". E) - 1E. Hence:

where the q2j are the reverse polynomials of the second kind, generated by the three term recurrence corresponding to Hj,j with the initial conditions qo = 1 and

4_ 1 = O. It is a trivial matter to complete the proof. Simply observe that det( \)!.,l ) = - -Yn{3n(A - a)" due to (11) and det(\)!k,J) = 'iJ,Pk - .q" from (12). 0 Thus the closed-loop sys te m with a reduced-order plant t' (' j anu a feedbad< ';,q,, (s) "" _1d·) (_ ) has all of its poles at


(1 .


2. Multiple point case Assume now t hat there are t distinct interpolation points {0"1' ... , ar} . The following generalization of theorem 3 holds. Theorem 4 Let ~l: - ll$~ and ~kl(")} be the transfer funcq.l: _ IS qkS tions corresponding to the T'educed-order mode.1.~ genf.rated by k - 1 and k iterations respectively of the rational Lanezos alg071th11l. If the set of chosen interpolation points {ai-i , ... )O'i;) ~ {ql •... ,ud , then the 2k roots of the determinant of

are all elements of the set {"I , ' . . , ",}. Proof: Wc will show that E(q;k.d

= E(e;•.I)

2k and is minimal (2J. Bnt. deg det('1'k,,) deg(s<7" )2(Vk
Combining the facts that E(8i. I) is the most powerful model of the dat" set D, 'deg dct,('1' •. I) ~ deg deg(8;k,,) and B (0'k ,l) <:;; B('1'.,,) implies that E(>Ji •.1 ) = 1:(8;,,1)' 0 From the proof of Theorem 2, it is clear that tbe polynomial matrix >Ji (s) generated from rational Lanczos is in fact It valid generating matrix 9;", (s) for t.he data set D . Due to its structure, rational Lanczos can only lead to generating matrices \\'hkh correspond to an even number of momeuts. Out unlike the methods of [2], rational Lanczos computes a generating matrix primarily through st.ate-space operations.



4. An Example

9 Zk ,l (s) is the generating matrix corresponding to the IlIo,t powerful llnfalsified model of the data set D = {Z(t,O'I,jd , ... ,Z(t,O'" J;j. By [2, Theorem 6.5], E(Wk.l) = E(0;k,l) implies that '1' •. I(S) and 0,.,1(') share invariant factors. But the proonct of the invariaut factors of 8 21.:,1 (5) is a polynomial whose roots are elements of t he desired set {"I , """' } [2, expression

In the sequel we present a simple example of rational interpolation illustrating th" relationship between the behavioral and Lanczos approaches. We are given the rational function


The first step in showing E(>Ji '·.d = E(9;, .,) is to demonstrate lhat. B(0;.,) <:;; B (wu). Ily definition, 0 2k ,] Z(t, ai, :h) = 0 for i = 1~ ... 1 f; wc must. prove t.hat


[p, qk 1Z(I,Oi,.J;)

= 0

for any Z(t,I1j. ,]d E D . Similarly, the expression [P'_I Q,_dZ(t, <7"j;) is zel'O, for all the time series except Z(t , OiZ ,)i). Bec.ause I-wo lIew moments were matched at Or.- b~t.wecn the (k - l)l
(')( ) _ .!. dkg(s) [J







= 0,

1, ... , 7

at s = -1 , turn out to be

1 1 5 35 7 21 15 16' 8' 32' 256 ' 64 ' 256' 256 Applying the behavioral algorithm (see [2]), we get the following sequenc.c of m at ric(~8 after 2,4,6 and 8 steps:

8 2(.,) =


c: 32

0.(.) = But this relation is in fact sufficient for our purposes because it and t.he inverse rec.ursion of (5) imply that

(s - 1)4

"Ve wish to recover it by means of interpolation at the point ." = -1. The quantitieR


for i = 1, ... , i as well. The value of j,j is twice the nUIllber of times 0i is chos(!fl as the intcrpolatlon point. By Theorem 2, ~ interpolates the moments corresponding to D so thal.


9(S) :=


-81 ) ,(-).(s) =



C+~ 8~OJ •



~ -200)'(~R(')= ("+ ~ 8-

1 2

_..:.. 160


Multiplication of the above updates gives 0 , .I(S) =

8 2(8)

for i = 1,2, ... , i. lIenee (13) does hold , The second item of importanc.e is t he complexity associated with the polynomial matrices e;/;:,1 (s) and '1' •. 1 (s). The compJcxity of E(!-l;k,li is deg det0;,,! =



-8(21 8' - 148 + 5)

,2+ l06 + 4 9 160


8,.! (8) ,5 4


+ lfs'l + 243 +

.,3 _ US2



+ !lS 10


_ .l.

T,-I F'.IT, = A,,!, T,- ' aR,1 = B"I, H",T. = C"l


This equivalence transformation is given by

= 88(S)8 •. 1(S) = 8,(S)<,-),(8) 8 4(8)8,(8) 256(75 3 - L4s 2

+ 12~:· + 708 2 + 2688 + 769


+ lIs



s-l-4 5 3 +6s 2 -48+1


The en~ries of the s~coIld row of the matrices Elk and 6 k ,l , for k = 2,4,6 ,8, are the quantities obtainerl hy the rational Lanczos algorithm . The state-space representation o f t.he a.bove calcu lations i.t; as follows. The rational matrix E>s,l (.<;) has a minimal realization CS,l} Aa,l, B s,l as follows:

0 1 0 0 0 0 0 0



0 0 0 0 1 0 0 1 0 0 I)



0 0


-769 -268 - 70 - 12 1 0

0 0 0 0 0 1 0 0

() ()

0 I)

0 0



0 0 Cl 0 (j 0 1 0 0 1

where lz anu O2 denote the 2 x 2 identity and zero matrices , respectively. The realizations CS ,ll As,], B8 , 1 and Ha,l, Fs,l, Gg,l are equivalent; this means that there exists a nonsingular 8 x 8 matrix T, such tha t

768 -2816 3584 -1792 -1 4 -6

where 91, g2 are the first and second columns of GS,l' In a similar way we can find equivalence transformations Tk hetween the realizatioIlf:! C k,l , Ak,ll Bk,l and Hk ,l, Fk ,l l Gk ,l for k = 4 and k = O.

References [1) A.C. Antoulas, Recursive 1Il0deling of discrete-time series, in P. Van Dooren and B.F. Wyman editors ,

[MA volume on Linear Algebra for Control, pp. 122 , 1993. [2) A.C. Antoulas and J.C. WiIlems, A behavioral approach to linear exact modeling, IEEE 1tans. Au-


tomatic Control, vo!. 38. pp. 1776- 1802,1993. Furthermore only the (1,1), (5,2) elements of B"I , and t.he (1 , 4), (2,8) elements of CB" are equal to 1, the rest

[3) K. Gallivan, E. Grimme and P. Van Doorcn, Asymptotic wavefo rm evaluation via a Lanczos

being 7.ero.

This realization can also be pllt together from the

method , Appl. Math. Lett. , vo!. 7, pp. 7[,-80, 1994.

realizations of the individual updates C'n A.l;, Bk l k = 2, 4, 6, 8. First notice th at C. = ilk = h k = 2,4 ,6,8.

[4) K. Gallivan, E. Grimrnc- and P. Van Dooren, Pad"

From th e formulM above it follows Ihat

A, = - (: 32

A. = Since

~,18 ) .

.4, =


(! -~:O),



:: )


_ (


= - (

8,'8, le6 s,' 1


Approximation of large-scale dynamic systems with

Lanczos meLhous , Proe. 88rd fEEE Conf. on Deci8ion and Control, (Lake Buena Vista, FL) , 1994.

[5) K. Gallivan, E. Grimrroe and P. Van Dooren, A


rational Lan czo,~ algorithm for model reduction, CSRD Report No. 1411, University of Illinois at Urbana-Champaign, IL, 1995.


[6) W.B. Gragg and A. Limlquist , On the partial realit lo llows t hat a real-

i.ation problem, Linear Algebra and Appl., vo!. 50, pp. 277-319, 1983.

ization H a, l, FR,l ) G foI ,l of e8,~(s) eall be obtained as



H ."

[7] B .N. Parlett, Reuuction to tridiagonal form and .4,










0, 0, GB,!















minimal realizations, SIAM J. Matrix Anal. Appl. , vo!. 13 , pp. 567- 593, 19!)2.

[81 A. Ruhe, Rational Krylov sequence methods for eigenvalue computation, Linea,- Algebra and Appl.)

vo!. 58, pp. 391-405, 19H4.

0, )