Stably Computing the Kronecker Structure and Reducing Subspaces of Singular Pencils A-λ for Uncertain Data

:= qr([Bo,Bz])

r "I

0 B ( f + l ) := Rf

3.2 Apply Pi to A2

Next we apply the RZSTR-kernel to A('+l)-XB('+l) and so on until nlfl=O or n,+l#O but r,+l=O, as in theorem 3.3. In steps 1.1 and 2.1 above we use epsu (a bound on the uncertainty in A and B that is given by the user) as a tolerance for deleting small singular values when computing nullities. All singular values smaller than epsu are deleted, but we always insist on a gap between the singular values we interpret as zeros and nonzeros, respectively. In section 5 we further discuss the influence of finite arithmetic to computed results. The second decomposition LISTR computes the left minimal structure and the Jordan structure of the infinite eigenvalues similarly to RZSTR b;r working on B - pA. The main difference is that column compressions are replaced by row compressions and we start working from the lower rightmost (south-east) corner of the pencil. As we said earlier it is also possible to apply RZSTR to B* - pA* to get the left minimal and infinity structures (see section 3.1). However RZSTR working on B*-pA* will produce a block lower triangular form of A-XB, so in order to keep GUPTFU block upper triangular we make use of LISTR. In summary the GUPTFU decomposition is computed by the following 6 reductions: 1: Apply RZSTR to A-XB. Gives the right minimal structure and the Jordan structure of the zero eigenvdue in A,,- XB,,:

ror-ABor* I

J. Dernmel and B. Khgstrom

3 00

0

P1*(A-hB)Q1 =

A1-XB1

The block structure of Ao,-hBo, is exposed by the structure indices nf and rf (column nullities, see theorems 3.1 and 3.3). 2: Separate the right and zero structures by applying RZSTR to Bo,-pAo,. Gives the dght minimal structure in A,-XB,:

YBr* I 0

P z * ( A o ~ - X B O ~ )=Q ~

A ~ - X B ~

Insist on the same right minimal indices as in reduction 1. 3: Apply RZSTR to A2-XB2, which contains the zero eigenvalues of A-XB, giving P ~ * ( A z - X B ~ )= Q ~Ao-hBo Note that reduction 2 destroys the block structure of the zero eigenvalues, but this reduction gives it back by insisting on the same Jordan structure as in reduction 1. Reductions 1 to 3 in summary:

The pencil Al-AB1 comes from reduction 1 and contains the possible left minimal structure and remaining regular part of A -XB. 4: Apply LISTR to A1-XB1. Since A1-XB1 cannot have any zero eigenvalues the reduction produces the block structure associated with &e left indices inAl-ABl:

5 : Apply LISTR to B 2 - g ~ .Gives the Jordan structure of the infinity eigenvalues h Al-XBi: =

P5*(A2-XB2)Q5

[" 0"

B3

3-

*

Af-Wf

1

6: Apply the QZalgorithm [27]to As- XB3 and a reordering process to the resulting upper triangular pencil to give the desired pencil Af-XBf: Pg* (A3-AB3)Qa = A / - XBf

*I

Let af and be the diagonal elements of Af and B f , respectively. Then the finite eigenvalues of A-XB are given by a f / p f . Reductions 4 to 6 in summary:

[

[

P6* 0 Ps* 0 0 I] 0 I]

[

Q5

P4*(A1-XBl)Q4

0

0 I]

[

Q6

0

0 I] =

[

Af-XBf

0

*

Af-XBf * 0 Ai-XBi

Note that the transformation matrices Pi and Q, in the 6 reductions above are all unitary and the composite transformations P and Q in GUPTRI are given by

30 1

Stably Coinpiring the Kronecker Structure

[A 4[A 13][ [ [ o 0 0 I]

Q4

Q = Q,

0 0 I]

Q5

Q6

(3.15a) 0 I]

(3.15b)

where the identity matrices in (3.15a-b) are of appropriate sizes. 3.4 Some history and other approaches During the last years w e have seen an intensive study of computational aspects of spectral characteristics of general matrix pencils A-XB. Already in 1977 Vera Kublanovskaya presented her first paper (in Russian) on the AB-algorithm for handling spectral problems of linear matrix pencils. For more recent publications see [16, 171. The AB-algorithm computes two sequences of matrices {Ak} and {Bk} satisfying AkBk+I = BkAk+I , k = 0 , 1 , . * * ;Ao=A , Bo=B where At+ and Bk+ are blocks (one of them upper triangular) of the nullspace of the augmented matrix c k = [Ak Bk] in the following way: N(C&)=

rBk+lI *

By applying the AB-algorithm to a singular pencil A-XB we get the right minimal indices and the Jordan structure of the zero eigenvalues. Different ways to compute N ( C k ) give rise to different algorithms. Kublanovskaya presents the AB-algorithm in terms of the QR (or LR) decomposition. In [18] a modified AB-SVD algorithm is suggested that more stably computes the rangelnullspace separations in terms of the singular value decomposition. In 1979 both Van Dooren and Wilkinson presents papers on computing the KCF. The paper [46] has already been mentioned in sections 2.1-2, where in algorithmic terms he derives the KCF starting from a singular system of differential equations. In [39] Van Dooren presents an algorithm for computing the Kronecker structure which is a straight forward generalization of Kublanovskaya's algorithm for determining the Jordan structure of A-XI, as further developed by Golub and Wilkinson [15] and Ki3gstrt)m and Ruhe [21, 221. His reduction is obtained under unitary equivalence transformations and is similar in form to the one obtained from the RGQZD algorithm. The main difference is that we in RZSTR and LISTR compute n, and nl-rf (di~n(N(A(~-~))) and dirn(N(A('-l)) nN(B('-I))), respectively) from two column compressions, while Van Dooren computes n, and r, from one column compression and one row com ression. Thus Van Dooren never computes a basis for tLe common nullspace of A(f- and B(I-I). The operation counts for one deflation step in Van Dooren's algorithm and in RZSTR are of the same order. For exam le if m = n the first deflation step counts 2n3+O(n2) for Van Dooren and n3+O(n ) for RZSTR. Note that both algorithms use singular value decompositions for the critical rank/nullity decisions. If one consider less stable and robust procedures for determining a rangelnullspace separation (like the QR or LR decompositions) it is possible to formulate faster variants of our algorithms. However this is not recommended except possibly in cases where we for physical or modeling reasons know that the underlying pencil must have a certain Kronecker structure.

!?

f

J. Demmel and B. Kdgstrom

302

4. Perturbation Theory for Pencils 4.1 Review of the Regular Case The perturbation theory of regular pencils shares many of the features of the theory for the standard eigenproblem A-XI. In particular, for sufficiently small smooth perturbations, the eigenvalues and eigenvectors (or deflating subspaces) of a regular pencil also perturb smoothly. One can prove generalizations of both the Bauer-Fike ([6, 10, 11, 361) and Gerschgorin theorems [33] to regular pencils. In order to deal with infinite eigenvalues, one generally uses the chordal metric x(X,X’) =

IX

- X‘I

(1 + X2)” * (1 + X’2)’/2 to measure the distance between eigenvalues, since this extends smoothly to the case A =

00.

In our context we are interested in bounds on the perturbations of eigenvalues and deflating subspaces when the size of the perturbation E is a (not necessarily small) parameter, possible supplied by the user as an estimate of measurement or other error. Thus, we would like to compute a decomposition of the form A-XB = P(S-XT)Q-’ , S- AT block diagonal as in (1.1),which supplies as much information as possible about all matrices in a ball P(E) { A S E - X(B+F) : II(E,F)II.
To illustrate and motivate our approach to regular pencils, we indicate how we would decompose P(E) for various values of E, where A-XB is given by

[:1:]

1100 0 1 0 [o 0 ll where -q is a small number. It is easy to see that IT, the spectrum of A - h B , is a={l/q,O,q}. The three deflating subspaces corresponding to these eigenvalues are spanned by the three columns of the 3 by 3 identity matrix. For E sufficiently small, the spectrum of any pencil in P(E) will contain 3 points, one each inside disjoint sets centered at Vq, 0 and q. In fact, we can draw 3 disjoint closed curves surrounding 3 disjoint regions, one around each X € u such that each pencil in P(E) has exactly one eigenvalue in the region surrounded by each closed m e . Similarly, the three deflating subspaces corresponding to each eigenvalue remain close to orthogonal. Thus, for E sufficiently small, we partition u into three sets, ul={l/q}, u2={0} and A-XB=

u3={d*

0 0 0 -A

v

As E increases to q / f i , it becomes im ible to draw three such curves because there is a pencil almost within distance q/ 2 of A-XB with a double eigenvalue at q/2:

E p 1;H] q;J-A

7

where t; is an arbitrarily small nonzero quantity. Furthermore, there are no longer three independent deflating subspaces, because the 5 causes the two deflating subspaces originally belonging to 0 and q to merge into a single two-dimensional deflating subspace. We can, however, draw two disjoint closed curves, one around

Stably Computing the Kronecker Structure

303

Uq and the other around 0 and q, such that every pencil in P(E) has one eigenvalue inside the curve around Vq and two eigenvalues inside the other curve. In this case {Vq} = u1u 0 2 . we partition u={O,q} As E increases to 1, it no longer becomes possible to draw two disjoint closed curves any more, but merely one around all three eigenvalues, since it is possible to find a pencil inside P(r) with q and Uq having merged into a single eigenvalue near 1, as well as another pencil inside P(E) where 0 and q have merged into a single eigenvalue near q/2, In this case we cannot partition u into any smaller sets. This example motivates the definition of a stable decomposition o f a regular pencil: the decomposition in (1.1) is stable if the entries of P , Q , S and T are continuous and bounded functions of the entries of A and B as A-XB varies inside P(E). In particular, we insist the dimensions nI of the SS-XTil remain constant for A-XB in P(E). This b corresponds to partitioning u=U u, into disjoint pieces which remain disjoint for

u

i=1

A - X B in P(E). We illustrated this disjointness in the example by surround each ul by its own disjoint closed curve. For numerical reasons we will also insist that the matrices P and Q in (1.1) have their condition numbers bounded by some (user specified) threshold TOL for all pencils in P(E). This is equivalent to insisting that the

deflating subspaces belonging to different ui not contain vectors pointing in nearly parallel directions. In the above example, as E grows to q / f i , there are pencils in P(E) where the deflating subspaces belonging to q and 0 become nearly parallel:

E

0 11/2.P q;J - A

[;;;]. T O O

The two right deflating subspaces in question are spanned by [0,1,OIT and [0,-V5,llT, respectively, which become nearly parallel as 6 approaches 0. The numerical reason for constraining the condition numbers of P and Q is that they indicate approximately how much accuracy we expect to lose in computing the decomposition (1.1) [8]. Therefore the user might wish to specify a maximum condition number TOL he is willing to tolerate in a stable decomposition as well as specifying the uncertainty E in his data. With this introduction, we can state our perturbation theorem for regular pencils: it is a criterion for deciding whether a decomposition (1.1) is stable or not. Theorem 4.1: Let A - X B ,

E

and TOL be given. Let u

u into disjoint sets. Define xI for l S l b as

=

b

U uI be some partitioning of

i= 1

wherep,, q,, Dif. and Difr will be explained below. The corresponding decomposition (1.1) is stable if the following two criteria are satisfied: max xl < 1

lslsb

and

304

J. Demmel uml B. Kugstrom

For a proof and further discussion see [ll]. The criterion (4.1) is due essentially to Stewart [32]. If we have no constraint on the condition numbers (i.e. TOL=m), then there is a stronger test for stability (see [ll] for details). The quantities pi,q l , Dif,(u,,u-cr,), and Difl(ul,u-u,) play a central role in the analysis of singular pencils, and will be discussed further in the next section. For now let us say Dif,(u,,u-ui) and Difl(ui,u-mi) measure the "distance" between the eigenvalues in u, and the remainder in u- u,,and that p i and q, measure the "velocity" with which the eigenvalues move under perturbations. Thus the factor multiplying E in (4.1) is velocity over distance, or the reciprocal of how large a perturbation is needed to make an eigenvalue of ui coalesce with one from u--0,. This bound is fairly tight, and in fact the factor b.max (Pi,q,) in (4.3) is essentially the least possible value for ,I the maximum of K ( P ) and K(Q) where P and Q block diagonalize A-AB, the center of P( ). The quantities Dif,, Difl, p , , and q1 may be straightforwardly computed using standard software packages. Also, nearly best conditioned block diagonalizing P and Q in (1.1) can also be computed. Therefore, it is possible to computationally verify conditions (4.1)to (4.3)and so to determine whether or not a decomposition is stable as defined above, as well as to compute the decomposition. 4.2 Why is the Singular Case Harder than the Regular Case? Our analysis of the regular case in the last section depended on the fact that eigenvalues and deflating subspaces of regular pencils generally change continuously under continuous changes of the pencil. This situation is completely different for singular pencils. As mentioned in section 1.2, the features of a nongeneric singular pencil may change abruptly under arbitrarily small changes in the pencil entries. In this section we illustrate this with some examples: As stated in section 1.2, almost any arbitrarily small perturbation of a square singular pencil will make it regular. The singular pencils form a lower dimensional surface called a proper variety in the space of all pencils. A variety is the solution set of a set of polynomial equations. In this case the polynomial equations are obtained by equating the coefficients of all the different powers of A in det(A - AB) to zero [44]. In fact, given a square singular pencil one can find arbitrariIy small perturbations such that the resulting regular pencil has its eigenvalues at arbitrary preassigned points in the extended complex plane [46]. For example, the pencil 0 1 0 -A 1 0 0 0 0 - A O O l = 0 0 - A I 0 [ o 0 11 is singular, and it is easy to see that the perturbed pencil

J [:I A]:-

+

[:a

:]

-"b -c-A 0 has determinant dA3+cA2+bA+a. Clearly we can choose a , b , c and d arbitrarily

0

small so that this polynomial has any three roots desired.

Stably Computiiig the Kroriecker Structure

305

Similarly, almost all nonsquare pencils have the same KCF, which will consists entirely of Lk blocks if there are more columns than rows and LT blocks if there are more row than columns. For example, the following pencil is nongeneric because it has a regular part with an eigenvalue at 1:

[i

0 1 0 [0 0 11 -

=

-[: i

but the following perturbation

-[: i

hO1]

+

a\

I!X]

01

has a generic KC.F consisting of a single L2 block for any nonzero a . In control theory the fact that a generic nonsquare pencil has no regular part implies that a generic control system is completely controllable and observable (SW section 2.3). 4.3 Perturbation Results for Pairs of Reducing Subspaces How do we do perturbation theory for nongeneric pencils in the face of the results of the last section? Clearly, we need to make some restrictions on the nature of the perturbations we &ow, so that the structure we wish to preserve changes smoothly. In particular, we will assume that our unperturbed pencil A-XB has a pair of left and right reducing subspaces P and Q, and that the perturbed pencil ( A + E ) - X ( B + F ) has reducing subspaces PEF and QEFof the same dimensions as P and Q, respectively. Under these assumptions we will ask how much PEF and QEF can differ from the unperturbed P and Q as a function of 11 (E,F)(( E: Theorem 4.2: Let P and Q be the left and right reducing subspaces of A-XB corresponding to the subset ul of the set of eigenvalues u of A-XB. Let E = 11 (E,F)II E . Then if ( A + E ) - X ( B + F ) has reducing subspaces PEFand QEF of the same dimensions as P and Q, respectively, where (where Difu(ul,u-ul), Difl(ul,u-ul), p and q will be defined below) then one of the following two cases must hold: Case 1: emax(P,PEF)Iarctan(x.(p + (p2- I)*~)) and a c t d x * ( q +(q2- 1IU2>). If x
5

ardan(2.x 1 4 ) . In other words, both angles are small, bounded above by a multiple of the norm of the perturbation E . Case 2: Either ~~JQ,QEF) 5

306

J. Demmel and B. K6gstrom

or

In other words, at least one of the angles between perturbed and unperturbed reducing subspaces is bounded away from 0. Note that the definition (4.4) of x is identical to the definition (4.1) for regular pencils. We will see when we sketch the proof of this theorem below how this result for singular pencils is a natural generalization of the result for regular pencils: Without loss of generality we assume A-AB is in the GUPTRI form defined in section 3:

: :I

A,- XB, A-XB = 0 Arcx-XBrcx 1

0

where A,-XB, has only L, blocks in its KCF Al-XB, has only L,' blocks in subscript I-), A,,,- AB,, is upper triangular and regular, pencil. Now repartition (4.5) as follows: All Al2 A - A B = 0 A,,]

[

AI-XBI

(4.5)

(i.e. all right minimal indices, hence the its KCF (i.e. all left minimal indices), and each * is an arbitrary conforming

-

'[

B11 Bl2 0 B22]

where All-ABll includes all of A,-XB,-and possibly part of Areg-XBrc , and Az2-ABzz contains all of Al-XBl and the remainder of A,,-XB,,,. In partidar, we assume All-XBll and Az2-XBzz contain disjoint parts (al and u 2 = u-ul respectively) of the spectrum u of A-XB. We denote the numbers of rows and columns of All-XBl, by m1 and n,;it is easy to see m+nl (with equality if and only if A,-XB, is null) and rn22n2 (with equality if and only if AI-XB, is null). In this coordinate system our unperturbed reducing subspaces P and Q are spanned by the column of [Zm1IOIT and [Zn1IOIT, respectively. Our first step in computing PEF and QEFis to blockdiagonalize the pencil in (4.6), the upper left block containing crl and the lower right block u2. Thus, we seek P and Q of the following special forms such that P-'(A-XB)Q =

(4.7a) (4.7b)

301

Stably Computirrg the Krorzecker Structure

This is a set of 2rn1n2 linear equations in nln2+mlmz unknowns, the entries of L and R . Since rnl
(A11A22;B11J 2 2 )

umin(Zu)

with the trivial consequence:

11 (LsR)llE

11 (A1Z&12)11

E

Dfa(A11A22$11p22)

*

Dif, is specified merely by choosing u1 and uz=cr-ul, permitting us to write Difu(al,a2) if A-XB is known from context or even just Dif, if u1 is known as well. If A-XB is singular and (4.8) has nonunique solutions, we will choose the (L,R) of least norm since it leads to P and Q as well conditioned as possible. Call this minimum norm solution (Lo,Ro) and denote (1+11 Loll 2)m by p and (I+ 11 Roll 2)uz by q . Just as Dif, is specified only by ul,so are 11 Loll , 11 Roll ,p and q . Now consider the perturbed pencil (A+E)-X(B+F). Premultiplying by P - I and postmultiplying by Q yields the pencil p-l

. ((A+E)-h(B+F))

'

Q

Eal

=

]

A 2El2 2 + ~ 2]2- A ~ l lF2l + F 1 lBz2+F22 Fl2 (4.9)

[A1lCEI1

We now seek PFk and Q E F of the forms

and QEF

=

[

- R ~

.

In,

such that premultiplying (4.9) by Pi; and postmultiplying it by Q E F blockdiagonalizes it. It is easy to see that if we can find such PEF and Q E F the perturbed pair of reducing subspaces P E F and QEFwill be spanned by the columns of

so that if Lz and R2 are small, the-angles between the unperturbed and perturbed subspaces will be small. Performing these multiplications and rearranging the result yields

[

(All+~ll-L1EZl) (I+RlR2)

~ 2 ~ ~ 1 1 + ~ l l ~ - ~ ~+ 2 E21 2 +~ L2El2R2 2 2 ~ ~ 2

(All+Ell)Rl-Ll(A22+E22) + El2 - LlEZlRl (A,,+ E22+LZE12) (I+R2Rl)

I

(4.10)

308

J. Demmel and B. Kbgstrom @11+

F11- W 2 1 )

(I+ R P 2 )

( ~ 1 1 + ~ 1 1 ) ~ 1 - ~ 1 ( ~ 2 2 ++ ~ 2F12 2)

- LlF2lRl

@22+ F22+ J52F12) (Z+R2R1)

1

Setting the upper right and lower left corners of the pencil two zeroes yields two sets of equations: (All+Ell)R1-Ll(A22+E22) = -El2+L1E21R1 (4.11a) (4.llb) (4.12a) J52@11+F11)-

@22+F22)R2

= -F21+L2F12R2

.

(4.12b)

(4.11ab) is a set of 2m1n2 nonlinear equations in n l n 2 + m l r n 2 unknowns, and (4.12ab) is a set of 2m2nl nonlinear equations in the same number of unknowns. If A-XB is regular ( m f = n f ) ,both these systems have the same number of equations as unknowns; this was the case considered by Stewart [32]. He shows that as long as Ef, and F!,are small enough, both (4.11ab) and (4.12ab) have unique solutions in a neighborhood of the origin. If A-XB is singular, (4.1lab) will be underdetermined and (4.12ab) will be overdetermined. Overdetermined systems have no solution in general; that this one does is guaranteed by our assumption that the perturbed pencil (A+E)-X(B+F) have reducing subspaces PEF and QEFof the same dimensions as P and Q. To solve (4.11ab) we make the identifications x = [L1hl], Tx =

(4.13a) (4.13b)

and (4.13~)

so that (4.11ab) becomes

Tx = g + $(x) . Expressing T in t e r m of JSronecker products yields

so that [32]

(4.13d)

309

Stably Computing the Kronecker Structure

(4.14a) (4.14b) and (4.14~)

so that (4.12ab) becomes Tx = g

+ +(x)

(4.14d)

Expressing T in terms of Kronecker products yields

I

(All+

EdT@Zrnz

-zn1@

(B 11+ F1d*@zrn2 -1.

I@

( A 2 2 +E22)

(B22 +F Z Z )

Swapping the first mlm2 columns with the last nlnz columns and negating the whole matrix, none of which changes its singular values, yields

[

Znl@(A22

+E22) - (All+ E I Y @ Z r n z

Zn,@P22+F22)

J

-( ~ 1 1 + ~ 1 d * @ Z q

which we recognize as the matrix whose smallest singular value is defined as ~22~11+~11;B22+~22,~11+~11)

2 ~f,(A22,A11~22,B,l) -

v5 II ( ~ l l , E 2 2 , F l l 9 F 2 2 ) 1 1 E *

The quantity Difu(A22,A11;82281J does not generally equal DifU(A11P~2$11,B2~) (unless the A,, and B,, are symmetric). In the interests of retaining our coordinate free formulation of our bounds, we therefore define D i f d U l P Z ) = D i f u ('422 sill $22 9Bld where the fact that Dif, depends only on u1 and u2 follows just as for Dif,. Also, one may prove that Dif,(u1,a2)>0 if and only if u1 and w2 are disjoint, as we have assumed. Thus UrninV) 2

~f,(ul,uz)

* II *

.

( ~ 1 1 , ~ 2 2 ~ 1 1 , ~ 2 E2 ) 1 1

Using these bounds on urn&), we can find bounds on the solution L, and R, of (4.11ab) and (4.12ab): Lemma 4.3: [ll] Consider the equations

J. Demmel and B. Khgstrom

3 10

Tx = g

and

+ +(x)

x = T+(g

+

(4.15a)

+(.I>

(4.15b)

where T + is the pseudoinverse of T. Assume that amin(T),11 $11 , and K

IIgll 11 $11 *

u&(T) <

1

11 gll

satisfy

4 .

Then equation (4.15b) has a unique solution x inside the ball (4.1%)

This solution x of (4.15b) has the following relationship to the solution 2 of (4.15a):

Case 1: If m=n, 2=x‘ is the unique solution of (4.15a) (in the ball). Case 2: If mn, and if (4.15a) has a solution 2 in the ball, then 2=f. Thus, (4.15b) may have a solution whereas (4.15a) may not. Furthermore, in cases 1 and 3, if (4.15a) has a solution which does not lie inside the ball in (4.15c), it must lie outside the ball: (4.15d) Substituting in values for umin(T), 11 gll , and 11 +I[ from equations (4.13abc) and (4.14abc) yields Theorem 4.2. 4.4 Perturbation Results for the Spectrum of the Regular Part If mZ=n2, then Azz-ABz2 is square and hence regular (otherwise it would have both Lk and LJ blocks in its KCF). Then after blockdiagonalizing the perturbed pencil (A+E)-A(B+F) as in (4.10) we see the perturbed pencil also has a square part [(A22+Ezz+LzE12) - w322+F22+L2F12)1 * (Z+RZRl) * If this square part is regular, we can use our bounds on IIL,II and IIR,II from the last section in conjunction with some of the perturbation theory for regular pencils in section 4.1 to derive bounds on the difference between its eigenvalues and the eigenvalues of the regular part of the unperturbed pencil A22-AB22: Theorem 4.4: Suppose that Case 1 of Theorem 4.2 holds. Suppose further that the block A22-AB22 is regular (i.e. rn2=n2>0). Then the spectrum of the perturbed pencil (A +E) -A (B +F) includes the spectrum of (A22+E’22)-Wzz+F’22) where

v5

II (E‘22J7’22)IlE 5 * 4 * II (E,F)llE . Similarly, if we instead assumeA1l-AB1l is regular (ml=nl>O), the perturbed pencil (A +E) -A (B +F) includes the spectrum of (A11+~’11)-~(~11+~’11) where l l ~ ~ ’ l l J ’ l l ~5 llv E5 . P

- II (EJ)II.

*

then the spectrum of

Stablv Computing the Kronecker Structure

311

See [ll] for a detailed proof. 4.5 Some Applications to Linear Systems Theory

Now let us interpret the results of the last two sections in terms of linear system theory: perturbation bounds for controllable subspaces and uncontrollable modes. Let C ( A , C ) denote the controllable subspace of the pair of matrices ( A , C ) as in section 2.3. As discussed in that section, this sFbsp$ce is simply the left reducing subspace corresponding to crl = 0 of the pencil A-XB = [ C k - X I ] . It is easy to s~ that this pencil can have neither infinite eigenvalues nor L l blocks in its KCF since B = [OF] has full rank; hence it can only have finite eigenvalues and L, blocks in its KCF. Also, one can see that the number of L, blocks is a constant equal to the number of columns of C. Thus, the assumption in Theorem 4.2 about the perturbed pencil having reducing subspaces of the same size is implied by the assumption that the perturbed system (A+EA,C+Ec) has a controllable subspace of the same dimension as C ( A , C ) . Also, the algorithms one uses
Then either Case 1: @,,(C(A,C),C(A+ EA,C+Ec))

5

~ctan(~*(p+(p~-l)~~))

with @,,(c(A,c),c(A+EA,c+Ec)) 5 Xctarl(2.X / p )

if x < l / 2 , or Case 2:

For a detailed proof, see [ll]. It is also easy to see that this result applies immediately to observable subspaces as well by duality [48]. Another feature of a control system ( A , C ) for which we can derive perturbation bounds using this approach is the spectrum of the regular part, also called the input decoupling zeroes or uncontrollable modes of the control system. Theorem 4.4 of the last section yields: Corollary 4.6 Assume we are in Case 1 of Corollary 4.5. Then the uncontrollable modes of the perturbed control system (A+EA,C+EC) are the eigenvalues of the pencil (Azz+E122)-~B22

where

II E’zzll E

5

fi

*

4 . II (EA,&)I1

E

*

312

J. Demmel and B. Ktigshom

In [40],P and Q are chosen in (4.1) so that P-IBQ = [rlO], implying Bz2=Z. Thus the problem of finding the eigenvalues of the perturbed pencil (A22+E122)-AB22 above reduces to perturbation theory for the standard eigenproblem. To illustrate these points consider the pencil

A - A i = [Cb-XI]

=

E5

A 2:

8].

0 0 3-A This pencil is in the canonical form (4.2) with rnl=l, rnz=2, n l = 2 and 4 = 2 . In other words All-ABll = 1-A] has KCF L1 and A22- ABzz = diag(2-A,3-A) is

a regular diagonal pencil with eigenvalues 2 and 3. The left reducing subspace P which is spanned by [1,O,OIT is the controllable subspace C(A,C) of the system (AS). The quantities of interest in Corollary 4.5 are Dif,(0,{2,3}) = Difi(0,{2,3}) .I12 , and p =q = 1. Thus Corollary 4.5 implies that if E A and Ec are such that 1 = dim(C(A,C)) = dim(C(A+EA,C+EC)) and 11 (EA,Ec)ll = .028.x where x < l then either Case 1: B,,(C(A,C) , C(A+EA,C+Ec)) 5 ~ c t a n ( x ) or Case 2: B,,(C(A,C) , C(A+EA,C+Ec)) 2 arctan(l/ = arctan(S77) = 5 2 Thus, for example, if we have II(EA,Ec)llE < .001, then any one dimensional controllable subspace of (A+EA,C+Ec) will either be within .036 radians of C ( A , C ) or at least .52 radians away from C ( A , C ) . A simple example of the latter situation is

e)

,o,

[ ,"

0

[c-kE,-w+E~-Az] = lo-'

2-A

0

3J: in which case the new controllable subspace is spanned by [0,1,OITand so is orthogonal to C(A,C). In case the perturbed controllable subspace falls into Case 1, then we can use Corollary 4.6 to bound the uncontrollable modes of the perturbed system: the perturbed modes are eigenvalues of the matrix

[;: :] +

E'

where 11 E ' l f ~ 5 ~ * I (EA,Ec)ll I E . Gershgorin's theorem supplies the simple bound th the perturbed uncontrollable modes lie in disks centered at 2 and 3 with radii 4 l J ( E , , E c ) l J s 5 xm.04. 4.6 Other Approaches

Stably

C'rirnpii turg

313

the Kroriecker Stntcture

Sun uses a somewhat different approach for deriving perturbation theorems for the eigenvalues of singular pencils [38]. He generalizes his results for regular pencils [37] by relating the perturbation of the eigenvalues of a diagonalizable A - XB and the variation of the orthogonal projection onto the column space R([A,B]*) to each other. His assumption that A-XB is diagonalizable means that A-XB may have neither Jordan blocks of dimension greater than 1 nor singular blocks other than Lo and LOT. Sun must also make restrictions on the perturbations in A and B . His assumptions are that the space W=R(Ap) contains no vectors orthogonal to Z=R(A+Ep+F) and that W'=R(A*(B*)* contains no vectors orthogonal to Z'=R[(A+E)*I(B+F)*]*. By using the chordal metric (see section 4.1) to measure the perturbations of eigenvalues he obtains a Wielandt-Hoffman and a Bauer-Fike type of theorem for singular diagonalizable pencils, showing that the eigenvalues of A-XB are insensitive to small perturbations in A and B. 5. Analyzing the error of standard algorithms for computing the Kronecker structure Using the perturbation analysis from section 4 this section analyzes the error of standard algorithms for computing the Kronecker structure from section 3. Here we focus on the GUPTRI form presented in section 3.3, but the analysis hold for any backward stable algorithm. 5.1 Distance from A-XB to the nearest pencil with the computed Kronecker structure Let (Y be a description of the KCF of A-XB (a list of left minimal indices, right minimal indices, and eigenvalues along with the number and sizes of Jordan blocks for each one). Then we define Sing,(A,B) = min{ll (E,F)IIE : (A+E)-X(B+F) has KCF a} (5.1) i.e. Sing,(A,B) is the smallest perturbation of A-XB (measured as II(E,F)IIE) that makes the perturbed pencil (A+E)-X(B+F) have the Kronecker structure a. One estimate of Sing,(A,B) comes from running one of the algorithms in section 3 for computing the KCF. We focus on the GUPTRI form and its RZSTR kernel discussed in section 3.3. The main contribution to errors caused by the finite arithmetic in RZSTR comes from the computation of nl and r,. For example the computation n,:=column-nullity(A(') , epsu) means that all singular values (Tk of A(') less than epsu are interpreted as zeros. In practice we make column-nullity(.;) more robust as follows: if ak
L

J

i.e the right hand side of (5.2) is an exact reduction of a perturbed pencil A(,)+,!?(')- X(B(')+F(')). Here IIE(')IIi is the sum of the squares of the n, singular values interpreted as zeros in the column compression of A(') (see step 1.1of RZSTR

314

J. Dernniel and B. Kagstrom

in section 3.3). Similarly 11 F(’)IIi is the sum of the squares of the n i - q singular values of B1 interpreted as zeros in the column compression of the n, columns of B(l)Vf) (see step 2.1 of RZSTR). We get similar contributions from all reductions. Since we only make use of unitary equivalence transformations when computing the GUPTRI form the resulting decomposition can be expressed in finite arithmetic as

P* ((A +E ) -A (B + F ) ) Q

=

(5.3)

where S2 is the sum of the squares of all singular values interpreted as zeros during the six reductions (see section 3.3). The equations (5.3-4) express that the GUPTRI form is backward stable and it computes a Kronecker structure (Y for a pencil C-AD within distance 6 from A-AB. This computed 6 is an upper bound on Sing,(A,B) for the value of the Kronecker structure (Y the algorithm reports. How good is this upper bound? We cannot expect 6 to be smaller than O(epsu max(llAIIE , llBllE)), since epsu is used as the tolerance for deleting small singular values. However if A-AB has a well-defined Kronecker structure (all a,,!are O(rnacheps max((lAIIE, 11 BIIE)) where macheps is the machine are O(l), where machepssepsu), then 6 will be of the roundoff error, and all size O(macheps max(l1 All , 11 Bll E ) ) too. The described situation shows a case where we have overestimated the uncertainties in A and B. On the other hand examples (see [39]) show that delta may be as bad as a root of the true distance (e.g. cU2 versus E), but some function of it may still provide a crude lower bound on &&(A ,B), much as the output of K%gstrt)m’s and Ruhe’s algorithm for the Jordan canonical form [21] provides a lower bound on the distance from the input matrix to one with reported Jordan form. This problem is under investigation. From the backward stability of the GUPTRI form we see that the computed bases for a pair of reducing subspaces (P,Q) correspond exactly to a nearby pencil C-AD within distance 6 . How to compute error bounds for these bases is the topic of the coming section. 5.2 Perturbation bounds for computed pairs of reducing subspaces In this section we apply the perturbation theory of chapter 4 to analyze the results of the GUPTRI algorithm (or any similar backwards stable algorithm) of chapter 3. Let A-AB be a pencil supplied by the user; it may or may not be singular. The algorithm of chapter 3 had the property that it found an exactly singular pencil A’-AB’ within a small distance 6 of A-XB. Theorem 4.2 showed that if we perturb A’-XB’ such that the perturbed pencil (A’+E)--X(B’+F) has reducing subspaces of the same dimension as A’-XB’, and if the norm of the perturbation 11 (E,F)II is small enough, then we can compute bounds on the maximum angle between the unperturbed and perturbed reducing subspaces. Let A denote the upper bound on 11 (E,F)II of Theorem 4.2:

Then clearly if any singular pencil A”-XB” with reducing subspaces of the same dimensions as A’-XB’ is sufficiently close to A-XB so that it is also within A of

315

Stably Cotnputirig the Kronecker Structure

A‘-XB’, the perturbation bounds of Theorem 4.2 will apply to it. For A”-XB” to be sufficiently close to A’-XB’ it clearly suffices that 11 (A“-A,B”-B)IIE < A-6, if A-6>0. This leads to the following algorithm: Algorithm 5.1: Perturbation Bounds for Computed Reducing Subspaces 1) Reduce A-XB to A‘-XB’ of the form GUPTRI using an algorithm such as the one in section 3.3 based on RGQZD. Let 6 bound the perturbation l l ( A f - A , B ’ - B ) \ l E . Let P’ and Q’ be reducing subspaces of A’-XB’ corresponding to the subset u1of the spectrum of A ’ - X B ’ . 2) Compute the bound A given above in equation (5.5). Let q = A-6. 3) I q>O, then if A”-”” is any singular pencil within distance d = 11 ( A ” - A 7 Z 3 ” - B ) ~ ~ E< q of A-XB with reducing subspaces Q” and P” of the same dimensions as Q’ and P‘, one of the two cases of Theorem 4.2 must

hold: Case 1:

and Omax(Q’,Q”)

5

arctan(-

6+d

A

(q+(q2-1)”))

,

with better bounds if A’ ’ -XB’ ‘ is within q/2 of A - XB Case 2: Either

or

4) If 1110,no perturbation bound can be made because the perturbation 6 made by the algorithm is outside the range of our estimate. One other interpretation of chapters 3 and 4 can be made. If the rank decisions made by the algorithm of chapter 3 are sufficiently well defined, i.e. at each reduction step there is a large gap between the largest singular value (Jk+ less than the threshold epsu and the smallest singular value uk greater than epsu: Uk

>>epsu >>uk + 1

then the algorithm will make the same rank decisions for all pencils in a disk around A - X B , and hence will compute the same Kronecker structure for the left and right indices and zero and infinite eigenvalues. In particular, the computed minimal and maximal reducing subspaces will have the same dimensions for all inputs in the disk around A - X B . Combined with Theorem 4.2, this means that if A-XB is intended to be singular, but fails to be because of roundoff or measurement errors, then the minimal and maximal reducing subspaces will be computed accurately anyway. Applied to linear systems theory (see sections 2.3 and 4 . 3 , this means that the algorithms of chapter 3 compute accurate controllable subspaces and uncontrollable modes of uncontrollable systems.

3 16

J. Demmel and B. Khgstrom

6. Numerical examples A Fortran 77 program implementing the GUPTRI form with accuracy-estimates of computed results (see sections 3.3 and 5 ) is under development. The inputs consist of a matrix pair (A,B), a bound epsu on the size of the uncertainty in ( A $ ) , and flags indicating whether a complete or just partial information (e.g. the minimal indices or the regular part of A - X B ) is desired. The output is a GUPTRI form (3.14) which is true for almost all pencils within distance delta of A-XB (the generic case), or else a nongeneric GUPTRI form lying within distance 6 from the input pencil A - X B . The bound 6 is expressed in terms of small singular values that are deleted ( C e p s u ) during the deflation process to compute the Kronecker structure cx (see section 5.1). Here we report results from a Matlab [26] implementation in double precision arithmetic (the relative machine accuracy macheps = 1.4E-17) of the basic reduction RZSTR used to transform A-XB to a GUPTRI form. The inputs are as for the Fortran program while the output of RZSTR is a block upper triangular form

[

*I

Aor - XBor 0 A1-XB1

where Aor-XBo, contains the Jordan structure of the zero eigenvalue and the right minimal indices. The pencil A1-XB1 contains the possible left minimal structure and the remaining regular part of A - X B . This decomposition is the first of six reductions used for computing the GUPTRI form and corresponds to what we get when replacing GSVD by GQZD in the reduction theorem 3.3. 6.1 A generic pencil Let A and B be 7 = m by 8 = n random matrices (generated in Matlab by the RAND-function) and ~ P S U= m a c h e ~ sm a ( ]I A I] E [I B ]I E ) RZSTR computed the following decomposition corresponding to one L7 block (see section 1.2): 9

Since we mainly are interested of the computed block structure we only display numbers in 4 decimals accuracy. However the zeros displayed above are all accurate to the rounding level. RZSTR does not know in advance that A - X B is generic, so it

317

Stably Compzttiiig the Kronecker Structure

has to go through eight deflation steps in order to realize that A - XB is equivdent to one L7 block. (The generic structure of a k by k+ 1pencil is one Lk block.) Below we display ni,r, and the smallest singular value of [A(')*,B(')*]*in each of the seven deflation steps. i n, ri umin([A(')*$(')*]*) 1 1 1 0.38 2 1 1 0.34 3 1 1 0.33 4 1 1 0.21 5 1 1 0.41 6 1 1 0.36 7 1 1 0.67 8 1 0 Locally at each step i the smallest singular value displayed is the 2-norm of the smallest perturbation of [A(')*$(')*]* that makes A(') and B(') have a common nullspace (ni-ri>O). So at least one of these singular values must be interpreted as zero (
-

A =

0.2055 0.6329 0.1249 0. 0. 0. 0.

0.3407 0.7626 0 5552 0. 0. 0. 0.

0.3004 0.9134 0.1967 0. 0. 0. 0.

0.2294 0.6493 0.2089 0. 0. 0. 0. B=

0.5271 1.0435 0.6309 0. 0. 0. 0.

0.8712 0.9253 0.9119 0. 0. 0. 0.

0.4016 0.5897 0.4095 0. 0. 0. 0.

["d' iz]

rl

0.5034 0.7699 0.5951 0. 0. 0. 0.

=

0.3013 0.7512 0.3976 0. 0. 0. 0.

1.2644 1.8991 1.3820 0.9856 0.3355 0.7987 0.5069

1S207 2.3651 1.7303 1.3820 0.6163 1.2348 0.8677

0.7374 1.1086 0.8401 0.6800 0.2280 0.5545 0.4382

0.5971 0.8535 0.6428 0.0489 0.1748 0.0001 0.0764

0.9958 1.3437 1.1080 0.1379 0.4933 0.0001 0.2155

0.4886 0.6911 0.5103 0.0454 0.1622 0 .oooo 0.0709

=

0.7325 0.9922 0.7249 0. 0. 0. 0.

J. Demmel and B. Khgstrom

3 18

are displayed only to four decimals accuracy. The pencils A1l-XB1l and A22-XB22 have the KCF's diag(Lo,L1,J2(0))and diag(LT,2.J1(w)), respectively. Since A - X B is known to full machine precision RZSTR is called with an epsu similar to the one in the preceding example. In order to obtain the complete Kronecker structure RZSTR is called twice, first with A - XB as input and then with the deflated pencil from the first reduction as input (see the discussion in section 3.3 for more details). The following indices determining a nongeneric Kronecker structure a1are computed:

number of L I - blocks

2 1

1 1

number of Lf-l blocks

3 0

0 1

Table 6.1 During the reduction process all singular values that should be zeros are zeros to the rounding level. What happens if we perturb the nongeneric pencil A-XB? Consider where E and F are full matrices and IIEllE is O(machepsV2 IIA116), IIBllE). In section 4.2 we pointed out that almost all (arbitrarily small) perturbations of a nongeneric pencil will make it generic, so ( A + E ) - X(B +F ) ought to be generic, and it is! When we apply RZSTR with the same epsu as for the unperturbed pencil it computes the generic structure as in section 6.1 (which we denote by a2). However if we consider (A+E)--X(B+F) as a pencil with uncertainties in the data of size 11 Ell and 11 FII , respectively, RZSTR computes the nongeneric Kronecker structure a1 corresponding to the unperturbed pencil A-XB. Below we display the computed reduction P*(A+E)Q = 8.4580e-10 1.5835e-9 3.3293e-9 .1.5798e-9 9.4764e-10 4.0688e-9 6.2048e-9

0.0 2.0000e-15 2.4458e-11 5.4961e-9 -6.8016e-10 2.3264e-10 3.0473~-10 3.9444~-9 -2.9634e-10 2.6389e-9 5.7650e-10 -2.3186e-9 0.0 0.0

["" 0

- P*E’Q

=

(6.la)

-7.0593e-1 1.7643 -4.4274 -4.4983e-1 1.4142e-1 3.7705e-1 -7.0718e-2 -5.2545e-1 -4.9481e-2 -6.7643e-3 -6.4429e-8 2.6144e-8 -6.1760e-1 -9.4730e-2 1.5776e-2 -2.5683e-8 1.9761e-8 1.7517 0.0 0.0 -2.2751e-8 1.0899e-8 1.4153 9.8271e-2 0.0 -4.1396e-8 2.0990e-8 3.6086e-1 1.5134e-2 -6.0993e-2 2.8080e-8 -2.7099e-8 1.5303 -4.0147e-2 -3.2022~-2

(6.lb) 0.0 0.0 -2.6917e-9 -1.0018e-9 -4.6944e-9 7.1359e-10 6.7613e-9

6.6146e-1 0.0 0.0 0.0 0.0 0.0 0.0

-1.4960e-2 8.5456e-2 0.0 0.0 0.0 0.0 0.0

8.3731e-1 4.5289e-2 8.0000e-15 -8.7672e-9 3.7934e-9 -8.5628e-9 1.5586e-8

2.5573 4.2848e-1 5.1278e-1 0.0 0.0 0.0 0.0

-2.4824 3.9301e-2 -3.8112e-1 4.8408e-2 -4.3052e-1 1.7979e-2 9.422e-11 6.8512e-10 -5.4277e-11 -3.6715e-10 2.3142e-9 8.3491e-9 6.0093e-1 -1.3563e-1

-4.2721e-2 2.4430e-2 -7.8591e-3 -9.2871e-10 4.3950e-10 6.6088e-9 -3.9087e-2

rows 1 to 4 and columns 1 to 5

of

Stably Compirtiizg t h e Kronecker Structure

319

P*((A+E)-X(B+F))Q and its block upper triangular form exposes the Jordan structure of the zero eigenvalue and the right minimal structure of A+E--X(B+F). The pencil Ail-XBll is identical to rows 5 to 7 and columns 6 to 8 of P*((A+E)-X(B+F))Q and corresponds to the Jordan structure of the infinite eigenvalues and the left minimal indices. Note that since we applied RZSTR to a transposed pencil in the second reduction All-XBll is block lower triangular. The perturbation matrices E’ and F’ in (6.la-b) correspond to the small but nonzero singular values interpreted as zeros during the deflation process to determine the Kronecker structure. However our Matlab implementation of RZSTR does not delete these small singular values in the process of determining the q s in Table (6.1), so the decomposition (6.la-b) is to the rounding level an exact similarity transformation of the pencil A+E--X(B+F) specified as input. How do we interpret the resulting decomposition? One way is to say that (6.la-b) is exactly nongeneric to the accuracy max (IIE‘II. , llF’IIE)-7e-8, with the Kronecker structure of A-XB (the unperturbed pencil). In the context of section 5.1

[

[

Bor BlZ 0 Bill from (6.la-b) is nongeneric to the rounding level and is an exact equivalence transformation of the perturbed pencil (A+E+E’)-X(B+F+F’). 6.3 Assessment of the (nongeneric) computed Kronecker structures In the following figure we try to visualize our nongeneric example: Aor A12 0 Ail] -

63 =

O(E”2)

A’-XB’ 61 = O(E) /”A-XB

-

Figure 6.1 - A surface of equivalent nongeneric pencils (same KCF) Here epsu=O(E) means that w e claim to know the input pencil to full machine accuracy and epsu=O(Ev2) means that we claim to know the input pencil to around half the machine accuracy. Let tii denote the distance between two matrix pencils (See Fig. 6.1) and corresponds to deleted singular values (see equations 5.3 and 5.4). First the unperturbed nongeneric pencil A-XB and epsu=O(E) were specified as inputs to RZSTR, and a nongeneric GUPTRI form with Kronecker structure a1 was computed. This form corresponds exactly to the pencil A‘-XB’ which is very close to A-XB; thus Sing,,(A,B) 5 E1=O(c). Then RZSTR was called with the perturbed (generic) pencil (A+E)-x(B+F) and epsu=O(e) as inputs, and as a result the generic GUPTRI form corresponding to (A+E)’--X(B+F)’ was produced. In this case Sing,,(A+E,B+F) Ia2=O(e)). Finally (A+E)-X(B+F) and epsu=O(Em) were specified as inputs to RZSTR which computed a nongeneric GUPTRI form with Kronecker structure a1 again, and corresponds to (A+E+E’)-X(B+F+F’). Thus Sing,,(A+E,B+F) 5 63=O(em). To summarize, A’-XB’ and

3 20

J. Dernmel and R. Kagstroni

( A+ E + E ’ )- X(B + F + F ’ )

have the same KCF and belong to the proper variety of nongeneric pencils. Let P’ and Q’ denote the computed pair of reducing subspaces corresponding to Aor-XBor in (6.la-b), and P and Q the corresponding pair of A ’ - X B ’ . What can be said about the maximal angle between these two pairs of reducing subspaces associated with the perturbed pencil (A + E) -X(B + F ) and the unperturbed A -XB , respectively? By applying the error bounds as described in section 5.2 the following results were obtained: max(p,q) = 3.2847 Dif, = 0.0315 Difl = 0.0147 A = 0.0013 Case 1: O,,(P,P’) 5 5.3E-5 radians (approximately 0.003 degrees) O,,(Q,Q’) I1.3E-4 radians (approximately 0.007 degrees) or Case 2: Om,(P,P’) 2 0.19 radians (approximately 10.8 degrees) O,,(Q,Q‘) 2 0.08 radians (approximately 4.6 degrees) Since we computed a GUPTRI form of both ( A + E ) - X ( B + F ) and A - X B we can show the actual angles between the computed pairs of reducing subspaces corresponding to the perturbed and unperturbed pencils: O,,(P,P’) = 1E-7 radians O,,(Q,Q’) = 5E-7 radians 7. Conclusions and future work

With applications of matrix pencils to descriptor systems, singular systems of differential equations, and state-space realizations of linear systems in mind (section 2) we have shown how to stably and accurately compute features of the KCF (for example minimal indices, pairs of reducing subspaces, regular part, and Jordan structures) which are mathematically ill-posed and potentially numerically unstable. The stability is guaranteed by a finite number of unitary equivalence transformations of A - X B , giving a generalized upper triangular (GUPTRI)form (sections 3 and 5.1). We make the mathematical problem well-posed by restricting the allowed perturbations E-XF of A - X B such that A - X B and ( A + E ) - X ( B + F ) have minimal and maximal pairs of reducing subspaces of the same dimension (section 4). This regularization leads to perturbation bounds for pairs of reducing subspaces and for the eigenvalues of the regular part (section 4) which can be used for analyzing the error of standard algorithms (section 5). The accuracy of the computed Kronecker features are confirmed via error bounds (sections 5 and 6 ) . Our future work will be concentrated on the completion of the our Fortran 77 program that implements the GUPTRI form and the error bounds. Then we want to experiment with the software on a large and diverse set of applications and encourage other people to make use of it. Hopefully this will lead to feedback permitting us to build in more application-based knowledge and focus the software on different and specialized applications (see sections 2.3-4). The first version of the program will

Stablji Computing the Kronecker Structure

32 1

treat all pencils as a general A - XB problem.

8. References

[l] D.J. Bender, "Descriptor Systems and Geometric Control Theory", PhD Dissertation, Dept of Electrical Engineering and Computer Science, Univ. of California Santa Barbara, September 1985 [2] D.L. Boley, "Computing the Controllability/Observability Decomposition of a Linear Time-Invariant Dynamic System, A Numerical Approach, PhD Dissertation, Dept of Computer Science, Stanford University, June 1981 [3] K.E. Brenan, "Stability and Convergence of Difference Approximations for Higher Index Differential-Algebraic Systems with Applications in Trajectory Control", PhD Dissertation, Dept of Mathematics, Univ. of California Los Angeles, 1983 [4] S.L. Campbell, Singular systems of differential equations, Vol. 1 and 2, Pitman, Marshfield, MA, 1980 and 1982 [5] S.L. Campbell, "The numerical solution of higher index linear time varying singular systems of differential equations", SIAM J. Sci. Stat. Computing, Vol. 6, NO. 2, 1985, pp 334-348 [6] K.E. Chu, "Exclusion Theorems and the Perturbation Analysis of the Generalized Eigenvalue Problem", NA-report 11/85, Dept of Mathematics, Univ of Reading, UK, 1985 [7] D.J. Cobb, "Controllability, Observability and Duality in Singular Systems", IEEE Trans Aut. Contr., Vol. AC-29, No. 12, 1984, pp 1076-1082 [8] J. Demmel, "The Condition Number of Equivalence Transformations that Block Diagonalize Matrix Pencils", SIAM J. Num. Anal., Vol. 20, 1983 pp 599-610 [9] J. Demmel, "Computing Stable Eigendecompositions of Matrices", to appear in Linear Alg. Appl., 1985 [lo] J. Demmel and B. K&strbm, "Stable Eigendecompositions of Matrix Pencils A-XB", Report UMINF-118.84, Inst of Information Processing, Univ of Umea, Sweden, June 1984 [ll]J. Demmel and B. K%gstrt)m,"Computing Stable Eigendecompositions of Matrix Pencils", Technical report #164, Courant Institute of Mathematical Sciences, 251 Mercer Str., New York, NY 10012, May 1985 [12] F.R. Gantmacher, The theory of matrices, Vol. 1 and 2, (Transl.) Chelsea, Ney York, 1959 [13] C.W. Gear and L.R. Petzold, "Differentid Algebraic Systems and Matrix Pencils", in [23] pp 75-89 [14] C.W. Gear and L.R. Petzold, "ODE Methods for the Solution of DifferentiaYNgebraic Systems", SIAM J. Numer. Anal., Vol. 21, 1984, pp 716-728 [15] G.H. Golub and J.H. Wilkinson, "Ill-conditioned Eigensystems and the Computation of the Jordan Canonical Form", SIAM Review, Vol. 18, No. 4, 1976, pp 578-619 [16] V.N. Kublanovskaya, "An Approach to Solving the Spectral Problem of A - X B " , in [23] pp 17-29 [17] V.N. Kublanovskaya, "AB-algorithm and Its Modifications for the Spectral Problem of Linear Pencils of Matrices", Numer. Math., Vol. 43, 1984, pp 329342 [18] B. K%gstrt)m, "On Computing the Kronecker Canonical Form of Regular

322

J. Dernniel and B. Kagstrorn

(A-XB)-pencils", in [23], pp 30-57 [19] B. mgstrdrn, "The Generalized Singular Value Decomposition and the General A - XB Problem", BIT 24, 1984, pp 568-583 [20] B. K%gstrt)rn, "RGSVD - An Algorithm for Computing the Kronecker Structure and Reducing Subspaces of Singular Matrix Pencils A - XB ," (Report UMINF112.83), to appear in SIAM J. Sci. Stat. Comp., Vol. 7, No. 1,1986 [21] B. K%gstrt)mand A. Ruhe, "An Algorithm for Numerical Computation of the Jordan Normal Form of a Complex Matrix", ACM Trans. Math. Software, Vol. 6, NO. 3, 1980, pp 398-419 [22] B. K8gstrt)m and A. Ruhe, "ALGORITHM 560, JNF, An Algorithm for Numerical Computation of the Jordan Normal Form of a Complex Matrix [El" ACM , Trans. Math. Software, Vol. 6, No. 3, 1980, pp 437-443 [23] B. K%gstrBm and A. Ruhe (eds), Matrix Pencils, Proceedings, Pite Havsbad 1982, Lecture Notes in Mathematics Vol. 973, Springer-Verlag, Berlin, 1983 [24] F.L. Lewis, "Fundamental, Reachability, and Observability Matrices for Descriptor Systems", IEEE Trans. Aut. Contr., Vol. AC-30, No. 5 , May 1985, pp 502-505 [25] D.G. Luenberger, "Dynamic Equations in Descriptor Form", IEEE Trans Aut. Contr., Vol. AC-22, No. 3, 1977, pp 312-321 [26] C. Moler, "MATLAE - An Interactive Matrix Laboratory", Dept of Computer Science, Univ of New Mexico, Albuquerque, New Mexico [27] C. Moler and G.W. Stewart, "An Algorithm for the Generalized Eigenvalue Problem", SIAM J. Num. Anal., Vol. 15, 1973, pp 752-764 [28] C.C. Paige, "Properties of Numerical Algorithms Related to Computing Controllability", IEEE Trans. Aut. Contr., Vol. AG26, No. 1, 1981, pp 130138 [29] C.C. Paige, "A Note on a Result of Sun Ji-Guang: Sensitivity of the CS and GSV Decompositions", SIAM J. Num. Anal., Vol. 21, 1984, pp 186-191 [30] C.C. Paige and M.A. Saunders, "Towards a Generalized Singular Value Decomposition," SIAM J. Num. Anal., Vol. 18, 1981, pp 241-256 [31] R.F. Sincovec, B. Dembart, M.A. Epton, A.M. Erisman, J.W. Manke, E.L. Yip, "Solvability of Large Scale Descriptor Systems", Final Report DOE Contract ET-78-GO1-2876,Boeing Computer Services Co.,Seattle Washington, 1979 [32] G.W. Stewart, "Error and Perturbation Bounds for Subspaces Associated with Certain Eigenvalue Problems", SIAh4 Review, Vol. 15,1973, pp 752-764 [33] G.W. Stewart, "Gershgorin Theory for the Generalized Eigenproblem Ax=zBx", Math. of Comp., Vol. 29, No. 130, April 1975, pp 600-606 [34] G.W. Stewart, "A Method for Computing the Generalized Singular Value Decomposition," in [23], pp 207-220 [35] G.W. Stewart, "Computing the CS Decomposition of a Partitioned Orthonormal Matrix," Num. Math., 1982, pp 297-306 [36] J-G. Sun, "Perturbation Analysis for the Generalized Eigenvalue and the Generalized Singular Value Problem", in [23], pp 221-244 [37] J-G. Sun, "Perturbation Analysis for the Generalized Singular Value Problem", SIAM J. NUm. Anal., Vol. 20, 1983, pp 611-625 [38] J-G. Sun, "Orthogonal Projections and the Perturbation of the Eigenvalues of Singular Pencils", Journal of Computational Mathematics, China, Vol. 1, No. 1, 1983, pp 63-74 [39] P. Van Dooren, "The Computation of Kronecker's Canonical Form of a Singular Pencil", Lin. Alg. Appl., Vol. 27, 1979, pp 103-141 [40] P. Van Dooren, "The Generalized Eigenstructure Problem in Linear System

Stably Computing the Kronecker Structure

3 23

Theory", IEEE Trans. Aut. Contr., Vol. AC-26, No. 1, 1981, pp 111-129 [41] P. Van Dooren, "Reducing Subpaces: Definitions Properties and Algorithms", In [23] pp 58-73 [42] C. Van Loan, "Generalizing the Singular Value Decomposition", SIAM J. Numer. Anal., Vol. 12, 1975, pp 76-83 [43] C. Van Loan, "Computing the CS and the Generalized Singular Value Decompositions", Numer. Math., Vol. 46, No. 4, 1985, pp 479-491 [44] W. Waterhouse, "The Codimension of Singular Matrix Pairs", Lin. Alg. Appl., Vol. 57, 1984, pp 227-245 [45] G.C. Verghese, B.C. Levy, T. Kailath, "A Generalized State-Space for Singular Systems", IEEE Trans. Aut. Contr., Vol. AC-26, No. 4, 1981, pp 811-831 [46] J.H. Wilkinson, "Linear Differential Equations and Kronecker's Canonical Form", In Recent Advances in Numerical Analysis, Ed. C. de Boor, G. Golub, Academic Press, 1978, pp 231-265 [47] J.H. Wilkinson, "Kronecker's Canonical Form and the QZ Algorithm", Lin. Alg. Appl., Vol. 28, 1979, pp 285-303 [48] W.M. Wonham, Linear Multivariable Control: A Geometric Approach, Second ed.,New York, Springer-Verlag, 1979 Acknowledgements The authors would like t o thank several sources for their support of this work. Gene Golub provided both financial support and a pleasant working environment to both authors for several weeks during summer 1985. James Demmel has received support from the National Science Foundation (grant number 8501708). Bo K8gstrBm has received support from the Swedish Natural Science Research Council (contract IWRS-FU1840-101). Both authors have been supported by the Swedish National Board for Technical Development (contract STU-84-5481). Finally, both authors wish to thank Jane Cullurn, Ralph Willoughby, and IBM for their support in helping the authors present their results at the Large Eigenvalue Problems conference of the IBM Europe Institute at Oberlech, Austria, in July 1985.