IEEE TRANSACTIONS O N CIRCUITS AND SYSTEMS, VOL. c~s-26, NO. 3, M A R C H 1 9 7 9 1 8 0 1 2 1 131 141 PI 161 171 1 8 1 191
M. E. Van Valkenburg, Network Ana&sis.
Prentice-Hall, 1955. Englewood Cliffs, NJ:
C. A. Desoer a n d E. S. Kuh, Basic Circuit Theov. New York: McGraw-Hill, 1969.
G. C. Temes .and S. K. Mitra, Modern Filter Theory and Design.
New York: W ilev. 1973.
C. F. Kurth, “Steady-state analysis of sinusoidal time-variant networks applied to equivalent circuits for transmission networks,”
IEEE Trans. Circuits Syst., vol. CAS-24, pp. 610-624, Nov. 1977. S. K. Mitra, Ana&sis and Sjvzthesis of &ear Act& Networks.
New YOI’ --. 2~: W iley, 1969.
G. S. MC jschytz, Linear Integrated Networks: Fundamentals. New
York: Van Nostrand-Reinhold, 1974.
F4.kA Jury, Sampled-Data Control Systems. New York: W iley, B. J. Hosticka, R. W. Broderson, a n d P. R. Gray, “MOS Sampled- Data Recursive Filters Using State Variable Techniques,” in Proc. Int. Symp. Circuits Syst., Phoenix, AZ, pp. 525-529, Apr. 1977.
[lo] I. A. Young, 0. A. Hodges, a n d P. R. Gray, “Analog NMOS sampled-data recursive filters,” 1977 IEEE Solid State Circuits Conf., pp. 156-157, Feb. 1977.
[11] G. S. Moschytz, Linear Integrated Networks: Design. New York: Van Nostrand-Reinhold, 1975.
+
Carl F. Kurtb (SMW-F’79), for a photograph a n d biography please see p a g e 1 0 4 of the February issue of this TRANSACTIONS.
+
George S. Moschytz (M’65-SM’77-F78), for a photograph a n d biogra- phy ,please see p a g e 1 0 4 of the February issue of this TRANSACTIONS.
O n the Determ ination of the Smith-Macmillan
F o rm’ of a Rational M a trix F rom
Its Laurent Expansion
PAUL M . VAN DOOREN, STUDENT MEMBER, IEEE, PATRICK DEWILDE, MEMBER, IEEE, AND
JOOS VANDEWALLE, MEMBER, IEEE
Abstmct-A novel method is presented to determine the Smitb- Macmillaa form of a ratfonai m X n matrix R(p) from Laurent expansions in its poles a n d zeros. Based o n that method, a numerically stable alge rithmisdedueed,whichusesonlyaminimalnumberoftermsofthe
Lmrent expansion, bence providing a shortcot with respect to cumbersome a n d unstable procedo~ based o n elementary transformations witb oni- modalar matrices.
The method can b e viewed as a generalization of Kublauovkaya’s algorithm for the complete solution of the eigemWucWe problem for
AI-A. From II system’s point of view it provides a bandy a n d numerically stable way to determhe the degree of a zero of a transfer function a n d unifies a number of results from multivariable realization a n d invertibiity theory. The paper presents a systematic treatment of the relation between the efgen-information of a traosfer function a n d the information contained in partial fraction or Laurent expansions. Although a number of results are kaown, they are presented in a systematic way which considerably sim- plifies the total picture a n d fhoduces in a natural way a number of novel tedlniqu~.
I. INTR~OUCTI~N
T
HE PROBLEM of efficient determination of the Smith-Macmillan form of a rational m X n matrixManuscript received June 19, 1978; revised July 17, 1978.
P. M. Van Dooren is with the Department of Electrical Engineering Systems, the University of Southern California, Los Angeles, CA.
P. Dewilde is with Afdeling Elektrotechniek, T. H. Delft, The Nether- lands.
J. Vandewalle is with the Department of Electrical Engineering, the University of California at Berkeley, Berkeley, CA.
R(p) does not s e e m to have received a great deal of attention in the past, although its importance as a key element in systems analysis a n d design can hardly b e denied. T h e classical method of using unimodular matrix m a n ipulations is cumbersome a n d not suited for numeri- cal computations, because it results in a n extraordinarily large number of polynomial m a n ipulations. In all methods based hereon, numerical stability is lost because pivoting is not based o n the coefficients of p but o n its power [I].
O n the other hand, a number of papers are devoted to the realization problem for system transfer functions a n d a host of algorithms have b e e n devised [2]-[7]. Another set of algorithms were proposed for system inversion both in the case of systems over a finite field [8]-[ lo] a n d in the case of systems over @ or Iw [ 111, a n d criteria for system invertibility where derived [12], [13]. Most of these methods require the handling of large size matrices a n d n o n e devote any attention to the numerical stability prob- lem.
An answer to the invertibility problem is needed, e.g., in coding .theory where o n e is interested in deciding whether the transfer function has a unique zero at infinity (of large degree) a n d if so, in determining the degree of that zero a n d the inverse of the matrix. Also, in invertibility theory o n e wishes to know whether there is actually a zero at infinity in which case the system cannot b e inverted. In 0098-4094/79/0300-0180$00.75 0 1 9 7 9 IEEE
both cases, the structure of the zero at a given point a n d
its numerical determination is of crucial importance. S(i) t Problems of system factorization from the synthesis
[14], [15] a n d from the analysis [16] point of view have also b e e n studied for coprime factorization [17] a n d spectral factorization [ 181, [19]. Most of these techniques are based in s o m e way or the other o n properties of the Toeplitz matrices resulting from a Laurent expansion of the transfer function ,3(p).
In this paper, we show that the Laurent expansion of R(p) in o n e of its poles/zeros (in this case they m a y coincide) provides in a very simple way all the informa- tion n e e d e d to determine the Smith-Macmillan form of R(p). T h e method as deduced is local a n d recursive, treats each pole/zero individually, a n d uses only the exact amount of information n e e d e d from a smallest possible partial expansion. As soon as all the relevant information is collected, the algorithm proceeds to a next point.
Fig. 1.
Moreover we present a n implementable version of the algorithm which is fast (because of the work savings) a n d numerically stable (only unitary transformations are per- formed o n the data). This fast version is also shown to boil down to Silverman’s structure algorithm [l l] when used o n state-space descriptions a n d also to the eigen- structure algorithm of Kublanovskaya (for a n excellent discussion of this algorithm, see [20]), when used o n U-A in o n e of the eigenvalues cx .of A. T h e present theory gives additional insights in both specific applica- tions. T h e introduced concept of rank function is also closely related to F o m e y ’s use of valuations [25]. T h e algorithm presented is in fact a handy way for computing these valuations.
W ith respect to a given R(p) a n d o n e of its poles or zeros (Y, we define a function S (i) which characterizes the Smith-Macmillan form completely as far as the oc- currence of the pole/zero (Y in its entries is concerned. Let A(p) b e factorized as A(p)=A,(p).G(p) with
A,=
I
% --r,r 1 0 m --r,n--r 1
a n d G=diag { gl,+**,g,,} (3) where the gi have n o pole or zero at cr (put gi = 1 for r<i<n). Let, for O<i<r:
S (i) = ui for i integer II. hELIMINARIES
Suppose R(p) is a general rational m X n matrix, then it has a unique Smith-Macmillan form A(p) [21], given by
R(P) = MP)A(P)N(P) (1)
where M a n d N are m x m , respectively, n x it unimodular matrices a n d
S (i) = a,, for i noninteger a n d where i + is the
upwards rounded version of i. J4)
T h e picture of S(i) is .thus a step function as shown in F ig. 1. Because of the divisibility properties of e, a n d -f,
S(i) is a nondecreasing staircase. This function also gives the order w,(w,) a n d the degree S,(S,) of the pole (zero) (Y ]71,]231:
T h e e, a n d A are m a n ic polynomials, whereby e, divides ei+i,J divides&,, a n d e, is mutually prime withA [21]. If a finite point (Y EC is a zero of any e,, then it is called a zero of R(p). If it is a zero of any&, then it is called a pole of R(p). (These definitions are. not exactly standard, but c o m m o n a n d logical; for a n extensive discussion see [22].) It should b e stressed at this point that a single point (Y can b e at the s a m e tim e both a pole a n d a zero of R(p), each with specific order a n d degree.
u p = -u, (if (I, < 0) wz = 0, (if a, > 0)
c 6, = - 2 ai
I-
O ,
<Q
8, = x ui. (5) a,>0T h e degrees can also b e defined as areas delimited by the graph of S(i) (shaded areas in F ig. 1).
ExampIe I: Let A,(p)=diag {(p - a)-*, 1, (p - CX)~, (p - CX)~}. T h e n S (i) is given as shown in F ig. 2. n
Suppose
R(p)= 5 Ri(p-a)i (6)
iE-1
182 IEEE TRANSACTIONS O N CIRCUITS A N D SYSTEMS, VOL. CAS-26, NO. 3, M A R C H 1979
I S(i)
Fig. 2. w P = 2I I
6p= 2 wz=3 6z= 6q(R, CX) the Toeplitz matrix (i > - r)
In the sequel we show how to obtain S(i) at any point (Y by working on T(R,a). Our main theorem retrieves this information from a smallest possible number of terms in the Laurent expansion of R(p) at the point (Y.
III. SMITH-MACMILLAN INFORMATION
We demonstrate that the rank information of the Toep- litz matrices T(R,a) completely determines A,(p). We introduce first some terminology.
Definition 3. I
Starting with (6) and using the notation defined by (7) we define the rank index p,(R,a) as (we drop R and (Y when superfluous because of the context)
pi P rank q-rank T-,. (8)
Hereby, we suppose the nonexisting q (i < - I) to have rank zero. It is easy to check then that
rank T. = C pi.
j= -, (9)
Definition 3.2
Let R(p) and Q(p) be two mXn rational matrices related by
R(P) = %)-Q(P)-~;(P) (10)
where E and F are m x m , respectively, n x n rational matrices. We then say that R and Q are similar at a if E and F are regular at (Y (E(a) and F(a) invertible).
Proposition 3.3
If R(p) and Q(p) are similar at (Y then they have the same rank indices at (Y.
Proof: Since E and F are regular at (Y they both have a Taylor expansion at a:
where E,= E(a) and FO= F(a).
Expanding R(p) and Q(p) in Laurent series at (Y, it is easy to see from (10) that they have the same order I for the pole (Y and that their coefficients satisfy
(12)
for i > - 1. Since E, and F, are invertible, the Toeplitz matrices built on E(p) and F(p) are also. From (12) it follows then that
rank q(R,cw)=rank q(Q,cy) (13)
for i > - I. n
Corollas 3.4
A rational matrix R has the same rank indices as its Smith-Macmillan form in any finite point.
Proof: Unimodular matrices have no finite poles or zeros and are thus, regular in any finite point (Y EC. A rational matrix is thus similar to its Smith-Macmillan form at any finite point (Y. The proof is then completed by
Proposition 3.3. H
Corollary 3.5
A rational matrix R has the same rank indices at (Y as A,(p) (where A, is given by (3)).
Proof R and A have the same rank indices at any finite point, and A and A, are obviously similar at (Y (see (3))
since G is regular at CX. n
Since R(p) and A,(p) have the same rank indices at (Y we can deduce their properties from the Toeplitz matrices q(Aa,o). Those have special properties because of the specific form of A,(p): 1) all rows of Ti(A,,cr) are either zero or have only one nonzero entry, and 2) the nonzero rows of Ti(Aa,o) are linearly independent. Because of this last property the following holds:
p,=rank c(A,,cr)-rank q-,(A,,a)
=rank [A-,;. . ,Ai] (14)
where A, is thejth coefficient in the Laurent expansion of A,(p) at CY. From the higher mentioned properties of TJA,,cw) it follows then that this rank equals the number of l’s in [A-,;** ,Ai] or also the number of powers 5 smaller than i in A,(p). This can immediately be read off from S(i) of R(p) at CY.
Example 2: Resuming example 1: A,(p) =diag {(p - CI-~, 1, (p - (Y)~, (p - a)‘} (see Fig. 3). H
If we associate a rank function C%(i) to the rank indices pi as follows (staircase function):
S(i) I I 3-p3 3 p3 1, IT P3 = 4 2- p2 2 p2 w P2 = 2 1 p1 p1 z l- Pl = 2 P O i c p-1 12 3 4’ PO=2 -1 - p- wP Pel= 1 -2 4
!
t Pm2’ 1 P3 = 4 w P2 = 2 z PI = 2 i P O = 2 wP Pel= 1 Pm2' 1 Pm3= 0 Fig. 3. R(i) L rank r=4 1 , 6P' , , i -4 -3 -2 -1 01 2 3 4 Fig. 4.a(i) = pi for i integer
S(i) = pi- for i noninteger where i - is the downwards rounded version of i
then obviously the % staircase is the reflection of the s staircase with respect to the bisectrice (see examples 1 a n d 2) except that 9l is also defined for i < -up a n d ,i >o,. Remark that 9 % is also nondecreasing.
Example 3: Using example 2 we obtain a(i) 1s shown
in F ig. 4. n
From this similarity between 9, a n d S the following results immediately follow trivially.
Corollary 3.6
W P - -min {ilp,#O} 1 (if>O) -
(if > 0) (where r is the rank of R(p), a n d can b e <m,n). Corollary 3.7
/
a,= 2 pi i= --a P
These easy results were derived earlier (e.g., [23]) but are put here in a form convenient for the sequel.
A rank search of the Toeplitz matrices q( R, a) gives all information about the occurrence of the pole/zero (Y in the Smith -Macmillan form of R(p). As soon as a rank index pk equals the normal rank r of R(p) the search can b e terminated. By then, a m inimal number of coefficients of the Laurent expansion at (Y has b e e n e m p loyed in the computation. In the next section we develop a fast recur- sive algorithm that performs this rank search in a numeri- cally stable way.
IV. TOEPLITZ RANK SEARCH
In this section we describe a recursive algorithm to transform a Toeplitz matrix ZJ to a matrix whose nonzero rows are linearly independent. W e use, therefore, row transformations to reduce a n arbitrary matrix A to the form
A- >P u.A= -
1 1
0whereby A - has p linearly independent rows. In the sequel we call such a transformation a row compression of the matrix A.
Since each Toeplitz matrix q(R, a) is a n embedding of any lower order Toeplitz matrix TJR, a) (for j <i), it can b e expected that the compression of these smaller Toeplitz matrices can b e used to transform T. This idea of making full use of the Toeplitz structure of q is applied in the following recursive Toeplitz rank search.
T h e algorithm, acting o n the sequence R-,, * * * ,&, de- fines a set of indices {p _ I,. . . , pi} which will b e proven to b e the rank indices of q.
Algorithm 4.1
1) j= -I; l&=R, fork=-l;..,i.
2) Construct-a row transformation q that computes the rank pj of Ri a n d puts it in compressed form
q.xj
=
3) M u ltiply & (for k =j+ 1; . . ,i) left with q a n d partition analogously to 2)
Rk- > fj qk= -
I 1
Rk+ b ’j’ 4) Update the & aslTk= Rky 1
I 1
Rk+’ for k=j+ 1;. * ,i. 5) Ifj<ithenj=j+l; g o to 2 else stop. nRemark that the symbols ik, Rk-, a n d Rk+ denote in each step a different object. In a more rigorous notation we should use a n additional index referring to the current step j, giving iko), RkwQ, Rk?. For sake of conciseness, we only a d d this index when confusion is possible.
184 IEEE TRANSACTIONS O N CIRCUITS A N D SYSTEMS, VOL. CAS-26, NO. 3, MARCH 1979
For better understanding we draw the actions of step j on the matrix
FA [ $A 1 lq, 1 . . . 1 R”i’il, ( R”,o’) 1. (15) Steps 2 and 3 then multiply f left with C$ such that
(16)
with Rj-Q having full rank 4. The new f matrix is then (by steps 4 and 5) one block shorter and looks like
(17)
which is (16) with the 9 bottom rows shifted to the- left and the last block then chopped off. After the definition of the ip+i) blocks this is also equal to
f= [ ⅈ” 1 &;I’ 1 . . . 1 R”i’jf,‘) 1 R”io’+‘) ]. (18)
This is repeated until f is vanishing (at j= i). In the algorithm the old version of ? is deleted and replaced by its new version (18) which also justifies the dropping of the step indices.
Theorem 4.2
Algorithm 4.1 produces the rank indices of the Toeplitz matrix q defined on the sequence R-,, * * * , Ri.
Proof: We prove this inductively. In the first step G= - I) the required rank index p _, = rank (R _ ,) is indeed computed by compressing the rows of R -,. We then multiply left each block row of Ti with the transformation
U-,(let UPdiag[U-,;.*,U-,I):
v-q=
and (18): -R--,(-l) ~RQ$--;)( . . . . (R;(-1) p*u*q= o 0 0p Iq-+
)pi.
(20) 1 I 0 ) 0 1 b-1By premultiplying K with the transformation P-U we have compressed the Toeplitz matrix T-, to the full rank matrix V-, and separated the remaining rows in the form of a new Toeplitz matrix (without any additional com- putations). We show by induction that this is maintained in each step.
Induction step j: Suppose that at the beginning of step j, ‘T,. is already transformed to j-1 >$-I= kzm,fk j-1 ‘kg,vk
R”(j)
witi
e)=
1 1
0’ R”,o”
\ do,
C21)
J RI,
I I
0 Ri- Ri+ RiI, RiT1 Rfo1 Rfh 1 RI, 0 >P-l 3 a 3 3 >v-1 (19)After permuting block rows of equal block length (the + and with q:.- i of full rank I?$-. , which is the rank of Tj- ,. blocks are permuted with the - blocks as indicated by Step j acts on c(J) analogously to (19) and (20) and the arrows), we obtain, using the indexed notation of (17) transforms it to:
0 i [ fy+ 1) ‘I * (22) I 0 0 1 9
Embedding this in (21), we again obtain (21) for j incre- mented by 1. Thereby, kJ. has the form (for s o m e X)
l&j-I > @ '
Since both Rim0 a n d I$, have full rank, 3 also has full rank equal to pi + %- ,. T h e row operations applied o n T i d o not affect the original column ordering of q nor of any of its e m b e d d e d Toeplitz submatrices (e.g., q). This m e a n s that q a n d vj have the s a m e row space a n d hence the s a m e rank. Therefore, b is indeed the required rank index a n d 5 is a row-compression of T j. n
Corollary 4.3
Algorithm 4.1 implicitly computes a row transformation U which transforms q to the following compressed form where the RjmQ blocks have full rank pi.
. . . 0 &:(-I) . . . R,O’) .
1
R ; (0
If-l Ipi . IPi i ’ k?-lvk (23) n It is useful at this point to discuss several aspects related to the algorithm.1) For efficient rank determination o n e is forced to use the singular value decomposition (SVD) [20]. T h e C$ matrices are then derived from the SVD of $ (* denoting Hermitian conjugate):
O n putting the singular values smaller than E, equal to zero o n e gets
(25) where the rows of Rj- are orthogonal to each other a n d
have norms greater than l . By using the row transforma-
tions q ? the overall transformation matrix U in (23) will also b e unitary, which guarantees the numerical stability of the Toeplitz rank search.
2) An SVD of T would require 0[ m%(f + i + 1)3] opera- tions while the fast version using the Toeplitz structure requires only 0[ m % ( I + i + 1)2] operations (moreover it yields the ranks of all lower order Toeplitz matrices).
3) Stop criterion:
a) W h e n only the rank indices p -,, . . . , pi are derived, then Algorithm 4.1 will compute them, when acting o n the sequence R-,; . . , Ri. This is for instance the case when o n e is only interested in the polar structure at (Y [2].
b) If the total Smith-Macmillan information is re- quired there are different possibilities.
i) T h e order d of the zero (Y is known. T h e n Algo- rithm 4.1 acting o n R -,, . * * , Rd utill give all required information.
ii) T h e order d of the zero (Y is not known but R(p) has a finite expansion at (Y. This is, e.g., always the case with polynomial matrices (see Section V). T h e n the algo- rithm will b e performed o n this finite sequence a n d can b e continued a n arbitrary number of steps, if o n e cares to shift in zeros at the end. It should b e remarked that the Toeplitz rank search does not stop when reaching the e n d of the sequence. Example 4 (Section V) wiIl show a n instance where Toeplitz ranks keep increasing even after the length of the sequence has b e e n reached. Yet, the algorithm will only process the significant blocks while the transformations o n zero blocks remain trivially zero.
iii) T h e order d is unknown a n d R(p) has a n infinite expansion. T h e n we start with a n arbitrary sequence R-,; . * ,& (i can also b e equal to - /) a n d perform Algorithm 4.1. If at the e n d of this, pi is not equal to the normal rank r, the sequence should b e enlarged to, e.g., R 3. * * ,Ri,Ri+l* Algorithm 4.1 can be adapted for such gr&ing sections but if o n e wants to avoid double work, s o m e intermediary results a n d transformations must b e kept in memory. This can b e d o n e efficiently with a m inimal memory content by following algorithm (using indexed notation).
Algorithm 4.4 1) j,= - I; 2) R!-‘, = R:
3) if>>-Z&enfori=-Iuntilj-1 do: a) m u ltiply ij@ with Vi a n d partition rows
186 IEEE TRANSACTIONS O N CIRCUITS A N D SYSTEMS, VOL. ~~-26, NO. 3, M A R C H 1979
4) construct a row transformation q that computes
the rank Q of &j(j) and puts it in compressed form: R(i)
i S(i)
5) if pj = r then stop
elsej=j+ I; go to 2.
The step indices used in the algorithm refer to earlier Fig. 5.
notations and help the reader to check the procedure. While Algorithm 4.1 computes the form (23) row by row, the algorithm above does it column by column. Therefore,
u _ A at x = 2 is
the procedure does not need to know a priori how many 1 -1 1
blocks to process. 2 -2 2
4) Extension for p = CCJ 1 -1 1
The Smith-Macmillan form of a rational matrix does
I
+I*@--2).
not reveal the pole/zero structure of the point p = 00 since Since the only eigenvalue of A is X=2, we can construct the transformations performed are not regular at infinity. the Smith form by looking only at the rank indices of Our approach with Laurent expansions however can still X = 2. As remarked previously, we can work on T = [B I I] go through. If R(p) has an expansion at infinity of the only, keeping in m ind that null blocks are following (refer
type to formulas (15) and following for notation).
R(p)=Ds’+D,-,p’-‘+...D,+D-,p-‘+... (26) Step 0:
1 0 0
then the rank search of Toeplitz matrices of the type u,F=
I I
-2 1 0 *F
-1 0 1
D/
- D-i
q=
I I
\A
I(27)
will give the pole/zero structure at p = ao. It is easy to
prove that this corresponds to the Smith-Macmillan redefine structure of the point A= 0 of the transformed matrix
R(X - ‘) (this way of coping with the point at infinity is not standard, but seems logical in present context).
V. EIGENSTRUCTLJRIZ ALGORITHM
and pa=1
(zeros are shifted in).
Suppose (Y is an eigenvalue of the r X r matrix A and we want to compute the Jordan structure of hl- A at (Y. It is known [I] that this is the same as the behavior at the point (Y of the Smith-form A(A) of that polynomial matrix. The Laurent expansion at (Y of that matrix is
Step UI f=
/I
1 0 0 0 1 -1 -1 0-f 1I
((III-A)+(&cw)ZA B+(h-a)Z. (28)By computing the rank indices of the sequence B,Z,O; . . ,O we thus obtain the required information. The rank function at(i) will give us the function S (i) of XI- A at (Y which also determines the Jordan structure of A at (Y since each elementary divisor (h-a)’ corresponds to a Jordan block of size 9 [l].
Example 4: In order to illustrate the basic idea, let us compute the eigenstructure for
A=[ I; ; I;]
which has a unique eigenvalue at A=2. Also we use elementary rather than orthogonal transformations for illustrative simplicity. In numerical practice, orthogonal transformations would be used. The Laurent expansion of
=[ -d -i i 1-p. i !l andp,=2
redefine
?= [ li -d % I i i g] (zeros are shifted in).
Step 2: F has leading block of full rank: p2=3. The rank function and the function S(i) are thus given as in Fig. 5. This means that the Smith form is A(A)=diag { 1, (h - 2), (h - 2)2} and A has a Jordan canonical form as (Jordan sizes 1 and 2):
2 0 0
A=M-‘.
H-1
0 2 I *M. IW e now show that the eigenstructure algorithm of Kublanovskaya [20] is similar to our rank search (except for s o m e additional saving of work because of the special type of expansion B + (A i a)Z).
This algorithm can be written as follows, keeping the same numbering as in Algorithm 4.1.
Algorithm 5.1
1) j=O; B”=B,
2) compute a SVD of B” a n d its rank 5 : B” = UZ V*, 3) compress rows of B” to full rank by a similarity
transformation,
5
nj
4) update &B,, ,
5) if nj>O thenj=j+ 1; g o to 2
else stop. n
T h e nj are the null ranks at each step a n d are equal to our 5. T h e 5 are not the rank indices (since the size of the block B” is decreasing) but are related to them:
In order to show the similarity between both algorithms
we remark that multiplication with a constant invertible
column transformation does not affect the rank search. This m e a n s that each block of f m a y b e right-multiplied by a n invertible matrix. Let us put Algorithm 4.1 in the form (15)-(17) using the matrix ?=[jIZ] a n d allow slight m o d ifications in order to obtain Algorithm 5.1.
Steps 2 a n d 3: T h e SVD of Z ? gives 5 a n d q.
a n d allowing a right transformation Uj
Steps 4 a n d 5: Redefine
An additional row transformation can zero out B,,*, a n d shows that we can work equally well o n ?= [B, llZ,] = [B 1 Z,] (this size reduction preserves null ranks, not rank indices).
Since pi = r - 5 = r - nj the stop criterion is indeed ni = 0
or pi = r, the normal rank of AZ, - A.
T h e use of similarity transformations in Alg_oiithm 5.1
maintains the unit matrix as second block in T, whereby o n e can work only o n B”. Also, when proceeding to a next
eigenvalue, an easy updating of B” is possible without
restoring the original size r [20]. These modifications
cannot b e performed in the general case treated in this paper, because at each new step novel information m ight
be shifted in, which is not the case here. Therefore, we
m a y assert that Algorithm 4.1 is a straightforward gener-
alization of Algorithm 5.1.
T h e connection with Kublanovskaya’s algorithm allows us to m a k e s o m e considerations about numerical proper- ties of the generalized eigenvalue problem. It is known that for the eigenvalue problem AZ- A a n expansion (cxZ- A) + (A - a)Z is n e e d e d in order to retrieve the Jordan information of A [20]. Indeed, when such a n expansion is not known, every attempt to retrieve a Jordan chain causes the m u ltiple eigenvalue to split u p in a cluster of not necessarily close eigenvalues [20]. Analo- gously in the case of the Smith-Macmillan information, we would e n d u p with all first-order poles a n d zeros when n o expansions are used as starting point.
W h e n using Laurent expansion of the matrix R(p), the problem can b e converted to a rank-definition problem for which the singular value .decomposition can b e successfully used. Just as in the eigenvalue problem [20] we see how the SVD allows us, in a certain sense, to clear off the errors that would cause the poles or zeros to split up. Therefore, the knowledge of the eigenvalue a n d a n expansion in that point are required.
VI. APPLICATIONS
In this section we discuss briefly a number of applica- tions of the theory previously treated. In each case we will refer to the relevant papers in the specific application. T h e numerical algorithm presented is to b e introduced wherever a Toeplitz reduction is to b e performed.
During revision of this paper other papers were pub- lished by Emre a n d Silverman [26]’ a n d by Van Dooren a n d Dewilde [2], covering part of the results of Section IV. However, the connection with the Smith-Macmillan in- formation a n d the generalized eigenvalue problem are not shown. It is exactly this relation that generates the wide variety of applications s u m m e d u p hereafter.
A) Coprime Factorization
W h e n R(p) is given in a partial fraction expansion P W
R(p)= i: qP)+%+R,(P)
i=l (30)
then we have the polar parts of the Laurent expansion in each pole explicitly given. T h e coprime factorization
R(P)= D -‘(P)+(P) (31)
can b e viewed [ 1 7 1 as the construction of a m inimal operator D(p) that displaces the finite poles of R(p) at infinity in the product
188 IEEE TRANSACTIONS O N CIRCUITS A N D SYSTEMS, VOL. CAS-26, NO. 3, M A R C H 1979
The m inimality of D(p) ensures the coprimeness between D(P) and N(P) 1171.
One proves [17], [2] that this can be done with the present structure algorithm when performed on each polar section R,(p) consecutively. The structure of the poles cui will indeed be reflected in D(p) since the poles cui must cancel in the product (32). The rank search of the Toep- litzes T- ,(R,,(Y~) will give all necessary information to build up D(p) and will also guarantee its m inimality [17]. B) Realization
Using the PFE (30) the problem of realizing D(p) boils down to determining its polar structure in each of its finite poles. In [2] we prove that the present ideas lead to a fast and stable algorithm for realizing ,each finite pole of R(p). A total realization is then merely a block arrangement of these partial realizations. The superiority of this algorithm above other realization algorithms, both in fastness and stability, is also demonstrated in that paper.
C) Inversion of Systems
The problem here [ 111, [24] is to find an inverse (left or right) for the system
R(p)= C(pZ-A)-‘B+ D (33)
where D is constant. If D is not invertible (left or right) transformations have to be performed on a growing Toeplitz matrix
q=[DeqJIJ (34)
until a modified D is found that can be inverted. It is easy to show that D + CBp i ’ + CABp -’ + . . . is an expansion of R(p) at infinity and that the Toeplitz search of Silver- man [l l] is nothing but our rank search at the point p = cc. This will thus deliver the zero structure of R(p) at p = co. As soon as D has full rank (corresponding to a rank index which equals the normal rank), the zero struc- ture is completely determined, and the system can be inverted. This alternative way of looking at the problem gives an easy understanding of the inversion problem. Moreover, when one has to cope with a polynomial D in (33) the extension is trivial since the expansion of R(p) is
R(p) = D&+ . . . D,p+D,+CBp-‘+CABp-‘+...
(35)
and the Toeplitz matrix becomes
D) Linearization of Polynomial Matrices
When using present ideas on the polynomial matrix
D(p)=D,+D,p+...Dp’ (37)
we can obtain a m inimal linearization pS- T which has the same zero structure as D(p). In order to insure m inimality, a rank search of
T-,= D t\t
[ 1
I
(38)
is required [27]. Linearizations of polynomial matrices have the advantage of bringing the probl.em in a suitable form for computation of its zero structure; numerically stable algorithms have been developed [28], [29] for com- putation of generalized eigenvalue problems of the type AS-T.
E) Eigenstructure of Polynomial and Rational Matrices When using jointly the realization algorithm B) and the linearization algorithm D) on an arbitrary rational matrix we obtain a linearization of R(p) in the sense that it has the same zero structure as R(p) [27]. When the m x n matrix R(p) has normal rar&r smaller than m and/or n this linearization also yields interesting results about the role of a left and/or right dual basis for the zero structure of the matrix [27], as well as a unifying theory for inver- sion of rational matrices.
F) Factorization of Rational Matrices
When factoring R(p) as R,(p) *R,(p) (eventually spectral factorization) a choice of poles and zeros must be done. Also the spaces which accompany those points [ 141, [ 161, [ 181 have to be determined. A handsome way of doing this is by determining the realization of R(p) with the present algorithm [2], and by computing the zeros with the linearization of R(p). A factorization of R(p) can then be computed (if it exists) by updating its realization so that it splits into the tandem system R,(p)*R,(p) [30].
VII. CONCLUSIONS
The numerically stable algorithm presented in this paper is derived in a particularly simple way from the Smith-Macmillan theory. The major concern of the paper is to provide a good insight in the algebraic structure of this generalized eigenvalue problem. The connection with the Kublanovskaya algorithm, e.g., shows how fragile the Smith-Macmillan information is and how heavily it relies on the knowledge of poles and zeros and expansions in these points. In some applications, where only partial results are requested, such information is indeed available (see applications C) and D)). In other cases one can solve the problem without the knowledge of poles and zeros but with the implications thereof (absence of multiple poles/zeros, eventual unstability, etc.). The algorithm has been implemented on computer and has been used in a
number of applications. Another interesting aspect of the theory is the linkage it provides between a number of seemingly disjoint topics: 1) the Kublanovskaya algorithm in numerical analysis as described in [20], 2) the realiza- tion theory of a transfer function in ABCD form, 3) the inversion theory as used in coding theory a n d the so- called structure algorithm, 4) factorization theory spectral or general, a n d 5) network realization theory through coprime factorization [19]. S o m e features of the theory m a y appear almost trivial, but seen against the prolifera- tion of difficult arguments found in the literature o n this very topic, it seems useful to stress the basic simplicity obtained.
ACKNOWLEDGMENT
W e thank G . Verghese for drawing our attention to L. M . Silverman’s work a n d for the useful remarks h e made.
1 1 1 PI 131 141 1 5 1 PI 171 PI 191 f101 f111 1 1 2 1 v31 P41 USI 1 1 6 1 1 1 7 1 UsI P91 PO1 W I REFERENCES
F. R. Gantmacher, Theoty of Matrices Z & ZZ. New York: Cl * relsea, 1959.
P. Van Dooren a n d P. Dewilde, “State-space realization of a general rational matrix. A numerically stable algorithm,” presented at Midwest Symp. Circuits Syst., Aug. 1977.
W. A. Wolovtch. Linear Multivariable Svstems. New York: Springer Verlag, 1974.
Y. L. Kuo, “O n the irreducible Jordan form realization a n d the d e
cr -17, pp. 322-332, Aug. 1970. e e of a rational matrix,” IEEE Trans. Circuit Theov, vol J. Rissanen, Realization of Matrix Sequences. San Jose, CA.: IBM Res. Lab., Apr. 1972.
B. L. Ho a n d R. E. Kalman, “Effective construction of linear state-variable models from input
technik, vol. 14, Heft 12, pp. 545- 4 48, 1966. output functions,” Regehngs- H. H. Rosenbrock, State-space and Multiuariable Theory.
London: Nelson, 1970.
E. R. Berlekamp, Algebraic Coding Theory. New York: MacGraw-Hill, 1968.
G. D. Fomey, “Convolutional codes I. Algebraic structure,” IEEE Trans. Inform. Theory, vol. H-16, pp. 720-738, Nov. 1970. M. K. Sain a n d J. S. Massey, “Invertibility of linear time-invariant dynamical systems,” IEEE Trans. Automat. Contr., vol. AC-14, pp. 141-149, Apr. 1969.
L. M. Silverman, “Inversion of multivariable linear systems,”
IEEE Trans. Automat. Cant:, vol. AC-14, pp. 270-276, June 1969. S. H. W a n g a n d E. J. Davtson, “A new invertibility criterion for linear multivariable systems,” IEEE Trans. Automat. Contr., vol. AC-18, p
A. S. W J
. 338-339, Oct. 1973.
sky, “O n the invertibility of linear systems,” IEEE Trans. Automat. Contr., vol. AC-19, p
J. Vandewalle, “Synthesis a n B . 272-274, June 1974. analysis of multivariable systems through factorixatton in frequency domain,” Ph.D. dissertation, Kath. Univ. Leuven, Louvain, Mar. 1976.
F. M. Callier a n d C. D. Nahum, “Necess
“il a n d sufficient condi- tions for the corn
in series using tl! e coprime factonxation of a rational matrix,” lete controllability a n d o servability of systems
IEEE Trans. Circuit $vst., vol. CAS-22, pp. 90-95, Feb. 1975. P. Dewilde a n d J. Vandewalle, “O n the factorization of a non-sin- gular rational matrix,” ZEEE Trans. Circuit Syst., vol CAS-22, pp. 637-645, Aug. 1975.
P. Van Dooren a n d P. Dewilde, ‘Coprime factorization of rational matrices using elementary factors,” Intern. Rep. Div. Appl. Math. &
Camp. SC., Kath. Univ. Leuven, 1975.
J. Vandewalle a n d P. Dewilde, “O n the minimal spectral factorixa- tion of nonsingular positive rational matrices,” IEEE Trans. Zn-
form. Th., vo;; IT-21, pp. 612-618, Nov. 1975.
zzf;tils A novel algorithm for spectral factorization,” Proc.
1975. ’ ummer Sch. Alg. Syst. New York: Sprmger, June G. H. Golub a n d J. H. W ilkinson, “Ill-conditioned eigensystems a n d the computation of the Jordan canonical form,” Siam Review,
vol. 18, n o 4; Oct. 1976.
B. Macmillan, “Introduction to formal realizability theory,” Bell. Syst. Tech. J., vol. 31, 1952.
1 2 2 1 1 2 3 1 1 2 4 1 1 2 5 1 P61 [271 Psi v91 1 3 0 1
A. G. Mac Farlane a n d N. Karcanias, “Poles a n d zeros of linear multivariable systems: A survey of the algebraic, geometric a n d complex-variable theory,” Znt. J. Control, vol. 24, pp. 33-74, 1976. J. Vandewalle a n d P. Dewilde, “O n the determination of the order a n d the degree of a zero of a rational matrix,” IEEE Trans. Automat. Contr., vol. AC-19, pp. 608-609, Oct. 1974.
P. J. Moylan, “Stable inversion of Linear Systems,” IEEE Trans. Automat. Contr., vol. AC-22, pp. 74-78, Feb. 1977.
G. D. Fomey, Jr., “Minimal bases of rational vector spaces with applications to multivariable linear systems,” Siam J. Control, vol. 13, pp. 493-520, May 1975.
E. Emre a n d L. M. Silverman, “New criteria a n d system theoretic interpretations for relatively prime polynomial matrices,” IEEE Trans. Automat. Contr., vol. AG22, pp. 239-242, A r. 1977. P. Van Dooren, ‘The generalized eigenvalue pro t lem. Applica- tions in Systems Theory,” Intern. Re
f;” Dept. El. Eng-Systems, Univ. Southern CA., Los Angeles, 1 9 7 .
C. B. Moler a n d G. W. Stewart, “An algorithm for the generalixed matrix eigenvalue problem,” Siam J. Numer. Anal., vol. 10, pp. 241-256, Apr. 1973.
R. C. Ward, “The combination shift Q Z algorithm,” Siam J. Numer. Anal., vol. 12, n o 6, pp. 835-853, Dec. 1975.
P. Dewilde, “M$imal multivariable cascade system synthesis: A ;tr$ure theory, presented at Midwest Symp. Circuit Cyst., Aug.
+
Paul M. Van Dooren (S’74) was born in Tienen, Belgium, o n November 5, 1950. He received the engineer degree in computer science from Louvain University, Louvain, Belgium, in 1974.
From 1 9 7 4 to 1 9 7 7 h e was a Teaching a n d Research Assistant at the Division of Applied Mathematics a n d Computer Science, University of Louvain. He is currently a Research Associate in the Department of Electrical En- gineering-systems, University of Southern Cali- fornia, Los Angeles, California. His main inter-
_. . . . .
ests lie in the areas of lmear system theory a n d numerical analysts. Mr Van Dooren is a student member of the Association for Comput- ing Machinery.
+
Patrick Dewllde (S’66-M’68) was born in Korbeek-Lo, Belgium, o n January 17, 1943. He received the engineer degree from the University of Louvain, Louvain, Belgium, in 1966, a n d the Ph.D. degree from Stanford University, Stan- ford, CA in 1970.
He has held research a n d teaching positions with the University of Louvain, the University of California, Berkeley, a n d the University of Lagos, Nigeria. He is currently Professor of Network Theory at the T.H. Delft, The Nether- lands. His main interests are in applied algebra (systems a n d network synthesis, factorization properties for several types of matrix functions, scattering matrix theory, a n d numerical analysis) a n d in teaching applied algebra.
+
Joos VandewaBe (S’71-M’77) was born in Kortrijk, Belgium, o n August 31, 1948. He re- ceived the engineering degree a n d a doctorate in applied sciences, both from the Katholieke Uni- versiteit, Leuven, Belgium, in 1 9 7 1 a n d 1976, respectively.
From 1 9 7 2 to 1 9 7 6 h e was Assistant at the E.S.A.T. Laboratory of the Katholieke Uni- versiteit Leuven, Belgium. From 1 9 7 6 to 1 9 7 8 h e was Research Associate at the University of California, Berkeley, where h e is currently Visit- ing Assistant Professor. His research interests are mainly in decomposi- tion a n d interconnection problems in system theory a n d general proper- ties of nonlinear networks.
