3aOI 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 772511892, USA E.J. Grimme Coordinated Seieuce Laboratory, Universit.y of Illinois at. UrbanaChampaign Urbana, IL 618012307, USA D.C. Sorensen Department. of Computational and Applied Mat.hematics, Rice University HOllston, TX 772511892, 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 construct a rational interpolant for the transfer function of a linear system defined by (potentially very highorder) 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 reducedorder, 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, JanuaryMarch 199:",.
both rational interpolation and reducedorder 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 nonsingular polynomial matrix 8'(s). The second approach, rational Lanezos [4 , 5], utilize:; a biorthog:onal projector IT composed of socalled 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.
3957
2. Computational Approaches to Rational Interpolation
is composed of time series taking on the special form [2J
=
Given an ordern singleinput 1 singi('output linear system described by means of the geupralized state space equations
gt(t)
= Ax(t) + butt),
yet)
= eT x(t)
Z(t,CT;,];) = eXPu,lt) LZj(t,CTi,]i) where
ZO(t,CT;,];)
(1)
= Ai(t) + butt),
Y(t)
= el'i(t) + duet)
(2)
Ji = 0, 1, ... Ji 
mj. (cr;)
;=
1, and £ = 1,2, ... ,L The t.erms
_eT {(A  cr;E)lE
r
·r"_·_ )I' ,J1 J 1.
l:S:j
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
z],(t,cr;,),
arc respectively the jI moments of the original and reducedorder 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 inputoutput 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,,],)}
I)!
t],j1
)
g(j)(,,il j!
+ 1) =
(
,lJ)~Ud
)
Ji!
h
D
gio;
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
for
Zj(t,CTi,li)= (
tJi 1
1») (); 
=(
0
the rational interpolation problem can be phrased as finding an orderA; < n model descrihed by E:i:(t)
(4)
j=O
and Zj(t,<7;,j) = 0 for j > };. It is also straightforward to verify that the inverse reclIrsion
Zj(t,lJ;,j)
= (:t
 "';)Zj(t,IT;,j
+ 1)
(5)
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 socalled 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).
(3)
3958
(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 trajectodes 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),
Si
[1
0]0; •. 1(8) =:
[)l.(s)
q.(s)]
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)
i=l
colsp(Wk )
= UKit+1 (ET(A<7iE)T,C)
(6)
(7)
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 pseudooptimal 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 reducedorder model of rational Lanc7,os. The reducedorder mod(~1 in (2) is chosen as
H.,.
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 reducedorder model does indeed solve the rational interpolation problem, For a proof see [5], Theorem 2 Given the reducedorder model (8) generated by the mtional Lanczos algorithm,
mj,(ai)
= 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* }.
"
Ui
,
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.
i=l
3, Relating the Rational Lanczos and Behavioral Approaches
where colsp denotes the column span. A Krylov subspace is defined as
3959
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 reducedorder 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
rence
In matrix form, this recursio n becomes
Pk(S)
(s _ ".)'P. _ l (S )
is equal (up to a constant) to (s  aj2k.
=
Proof: We begin by recalling (8), the defini tion of the reducedorder 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 wellknmvu that the poles of this reducedorder 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) = (JjJlj l(A  0")2
det
and
(11)
(}]>jij(,\)) = (}] (Jj I'YJI) (A_a)"
Finally) inLrocluce the 2 x 2 matrix poly nomial recursion
gives the n:vel'se polynomial
peA).
:= (A _ a)'Pk
(A  a )kdet (Hk,k 
'MA) ) = (}11Ji k j (A) ) 1Ji 0 , 1
(_1_) = _II) =
",.(A)
A a
(12)
where \)!0.1 is initialized wit.h t/Jo = Po = 0 and 0 = ,) = deg(p.) + 1 and
Aa det ( A  rr)Hk,k  I ) =
det (AHk ,k  (aHu
+ I))
(10)
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 closedloop sys te m with a reducedorder plant t' (' j anu a feedbad< ';,q,, (s) "" _1d·) (_ ) has all of its poles at
3960
(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'educedorder mode.1.~ genf.rated by k  1 and k iterations respectively of the rational Lanezos alg071th11l. If the set of chosen interpolation points {aii , ... )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.atespace operations.
'.1
where
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
6.7aJ.
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 Iwo lIew moments were matched at Or. b~t.wecn the (k  l)l
(')( ) _ .!. dkg(s) [J
8

k!
d.'3k
'
k
= 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:~
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
(13)
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.
1
9(S) :=
32
81 ) ,().(s) =
s+
2'
C+~ 8~OJ •

160
~ 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;,,! =
3961
33
8(21 8'  148 + 5)
,2+ l06 + 4 9 160
(
8,.! (8) ,5 4
Tt
+ lfs'l + 243 +
.,3 _ US2
,
LO
+ !lS 10
)
_ .l.
T,I F'.IT, = A,,!, T, ' aR,1 = B"I, H",T. = C"l
10
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

3»)
sl4 5 3 +6s 2 48+1
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 statespace 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
=
A,,!
0 0 0 0 1 0 0 1 0 0 I)
()
0
0 0
()
769 268  70  12 1 0
0 0 0 0 0 1 0 0
() ()
0 I)
0 0
I)
I)
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 discretetime 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
4
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, =
160
(! ~:O),
eo,:
=
:: )
~.
_ (
As
=  (
8,'8, le6 s,' 1
_
Approximation of largescale 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
2
rational Lan czo,~ algorithm for model reduction, CSRD Report No. 1411, University of Illinois at UrbanaChampaign, IL, 1995.
~!o
[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. 277319, 1983.
ization H a, l, FR,l ) G foI ,l of e8,~(s) eall be obtained as
follows:
Fs,1
H ."
[7] B .N. Parlett, Reuuction to tridiagonal form and .4,
f,
0,
04
0,
,4,
I,
04
=
=
0, 0, GB,!
0,
02
0,
0,
I,
0,
0,
=
A.
I4
0,
02
A.
I,
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. 391405, 19H4.
0, )
3962