Two connections between the SR and HR eigenvalue algorithms

Two Connections Between the SR and HR Eigenvalue Algorithms Peter Benner

and Heike F&bender*

Universitiit Bremen Fachbereich 28357

3-Mathematik

Bremen,

and Informatik

FRG

and David

S. Watkins+

Department

of Pure and Applied

Washington

State University

Pullman,

Washington

Mathematics

99164-3113

Submitted by Volker Mehrmann

ABSTRACT The

SR and

HR

calculating

eigenvalues

connections

between

on

a 2n x 2n

equivalent

algorithms

to an iteration

butterfly

of the i =

of matrices.

of GR This

(1) An iteration

using

shifts

algorithms

paper

makes

for two

of the SR algorithm

kLi, /.L;‘,

on an n X n tridiagonal

i = 1, . . , k. (2) An iteration

Hamiltonian

p;,

matrix

of the HR algorithm

matrix using shifts p, + /I;‘, 2 n X 2 n J-tridagonal matrix using shifts

of the family

subspaces

the SR and HR algorithms:

symplectic

lent to an iteration

are members

and invariant

i = 1,

. . , k, is

sign-symmetric

of the SR algorithm

matrix using shifts p,, - pi, i = 1,

HR algorithm on an n X n tridiagonal . . , k. 0 1998 Elsevier Science Inc.

on a

. . . , k , is equivasign-symmetric

1,.

*E-mail: [email protected], [email protected] 'E-mail: [email protected]&ngaddress:683524thAvenueNE, Seat&WA 981157037.Supported by the National Science Foundation under grant DMS-9403569.

LINEAR ALGEBRA AND ITS APPLICATIONS 0 1998 Elsevier Science Inc. AU rights resewed. 655Avenueoftbe Americas,NewYork,NYlOOlO

272:17-32

(1998) 0024-3795/98/%19.00

PIIsOO24-3795(97)00279-6

P. BENNER,

18

H. FASSBENDER,

AND D. S. WATKINS

INTRODUCTION

The SR algorithm

[7, 61 and the HR algorithm

family of GR algorithms

[25] for calculating

spaces

oldest

of matrices.

algorithm

[M-20],

The

member

structures.

matrices,

The

SR algorithm

and the HR algorithm

sub-

of the family is Rutishauser’s

LR

are useful because preserves

preserves

ces. In this paper we prove two interesting

of the

and invariant

and the most widely used is the QR algorithm

24, 261. The SR and HR algorithms special

[4, 5] are members

eigenvalues

they preserve

Hamiltonian

sign-symmetric connections

19, 11, 23, certain

and symplectic tridiagonal

between

matri-

SR and HR

algorithms: (1) An iteration

of the SR algorithm

matrix using shifts Z.L~,&‘, HR

algorithm

pi + p,Il,

on a 2n X 2n symplectic

i = 1,. . . , k, is equivalent

on an n X n tridiagonal

butterfly

to an iteration

sign-symmetric

of the

matrix using shifts

i = 1, . . . , k.

(2) An iteration

of the SR algorithm

nian matrix using shifts pi, -pi, the HR algorithm

on a 2n X 2n J-tridiagonal

i = 1, . . . , k, is equivalent

on an n X n tridagonal

sign-symmetric

Hamilto-

to an iteration

of

matrix using shifts

CL&,i = 1,. . . , k.

1.

SR BASICS

The

SR

algorithm

2n X 2n. Throughout

is applicable

to real

matrices

of even

the paper we express such matrices

dimensions

as block matrices

A=

in which the blocks

Aij are always n X n. Let

where I denotes the n X n identity matrix. A matrix S E (W2nx27’is symplectic if STJS = J ( or equivalently, SJS’ = J). If S is symplectic, then S is nonsingular, and S - ’ = JS’J T. The symplectic matrices form a group. The eigenvalues of symplectic matrices occur in reciprocal pairs: If A is an eigenvalue

of S with right eigenvector

x, then A-’ is an eigenvalue

of S with

SR AND

HR EIGENVALUE

left eigenvector

19

ALGORITHMS

(IX>‘. Symplectic

eigenvalue

problems

arise in discrete-time

control, filtering, and estimation problems (see, e.g., [12, 14, 16, 211 and the references given therein) and the computation of discrete stability radii [lo]. A matrix

Rll

is said to be J-triangular

R,,

of a J-triangular

The product of ]-triangular gular matrices For

R,,

Ril are all upper triangular and

if the submatrices

is strictly upper triangular.

and columns

1

Rl2

R= [ R,,

If one performs

a perfect

shuffle of the rows

matrix, one gets an upper triangular

matrices

is J-triangular.

The nonsingular

matrix. J-trian-

form a group.

the purposes

of this paper,

trivia2 if it is both symplectic

a matrix

and J-triangular.

M E [W2nx2n will be called M is trivial if and only if it has

the form

McC [0 where C and F are diagonal matrices, Almost

every matrix

exists, then

other

C nonsingular.

A E RznX2”

A = SR, where S is symplectic

sition

F

c-’ 1 ’

can be decomposed

and R is J-triangular

SR decompositions

of

trivial matrix, S = SM, and fi = Mof A. If A is nonsingular,

SR decompositions.

’R,

then

A = SE is another

In other words, the SR decomposition

sition at each iteration.

SR decom-

then this is the only way to create

trivial factors. The SR algorithm is an iterative algorithm

tion function

A can be built from it by S and R. That is, if M is a

passing trivial factors back and forth between position

into a product

[S]. If this SR decovnpo-

that performs

is unique

other up to

an SR decompo-

If B is the current iterate, then a spectral transformu-

q is chosen

and the SR decomposition

of q(B)

is formed,

if

possible:

q(B) Then the symplectic

= SR.

factor S is used to per-forma

similarity transformation

on

B to yield the next iterate, which we will call B:

6 = S-'BS.

(1)

20 We

P. BENNER, H. FASSBENDER,

shall assume throughout this paper that q(B)

AND D. S. WATKINS

is nonsingular.

Nothing bad

happens in the singular case [3, 251; we are avoiding it here solely to simplify the discussion. Since q(B) is nonsingular, S is determined up to a trivial factor, so B is determined

up to similarity transformation

by a trivial matrix.

HR BASICS

2.

Now consider matrices in R “’ “. A signature matrix is a diagonal matrix D = diag{d,, . . . , d,} such that each di is either 1 or - 1. Given a signature matrix D, we say that a matrix A E [wn Xn is D-symmetric if ( DA)T = DA. A tridiagonal matrix T is D-symmetric for some D if and only if Iti + 1,i 1= 1ti, i + 11 for i = I,..., n - 1. Every irreducible tridiagonal matrix is similar to a D-symmetric matrix (for some D> by a diagonal similarity with positive main-diagonal

entries.

the unsymmetric

D-symmetric

Lanczos

process

tridiagonal

matrices

can be generated

by

[I3], for example.

Almost every A E (w”‘” has an HR decomposition A = HU, in which JJ triangular and H satisfies the hyperbolic property H TDH = D, D is another signature matrix [5]. For nonsingular A the HR decom-

is uppe; where position

is unique

up to a signature

insisting that the upper triangular The

HR algorithm

matrix.

factor

[4, 51 is an iterative

decomposition.

Choose

is nonsingular,

and form the HR decomposition

use H to perform

a similarity

can make

it unique

by

function

on the HR p for which p(A)

of p(A),

if possible:

process

a spectral transformation

p(A) Then

We

U satisfy uii > 0, i = 1, . . . , n. based

= HU. transformation

on

A to get the next

iterate: fi = H-‘AH. The HR algorithm has the fo$wing structuze preservation property: If A is D-symmetri: and H TDH = D, then A is D-symmetric. If A is tridiagonal, then so is A.

3.

THE

SYMPLECTIC

We return our focus symplectic matrices form structure. symplectic.

CASE

to symplectic a group, the

matrices

in Rznx2”. Because the preserves symplectic

SR algorithm

That is, if the initial matrix is symplectic,

then all iterates

will be

SR AND

HR EIGENVALUE

A symplectic

21

ALGORITHMS

matrix

is called a bzttierfly

matrix if B,, and B,, are diagonal, and B,; and B,, are Banse and Bunse-Gerstner [l, 21 showed that for every symplectic

tridiagonal.

matrix M, there exist numerous is a symplectic

butterfly

form: If B is a symplectic An unre&ced submatrix

B,,

symplectic

matrix. The butterfly

symplectic

is irreducible

S such that B = S ’MS preserves

the butterfly

matrix, then so is 6 in (1) [I, 21.

butterfly

matrix is one for which the tridiagonal

[3]. Using the definition

one easily verifies that if B is unreduced, nonsingular.

matrices

SR algorithm

of a symplectic

matrix.

then the diagonal submatrix

This allows a decomposition

of B into two simpler

is

B,,

svmplectic

matrices:

where

is an unreduced,

T = B,‘B,,

decomposition

of an unreduced

matrix times a butterfly

3.1.

A Canonical

Form for Symplectic

The outcome

it is determined and J-triangular)

of an SR iteration It is therefore

matrix

0 I I

Butterjly

butterfly

matrix. into

This

a trivial

-I is unique T I

[3].

Matrices

form is preserved

by the SR

is not quite uniquely determined;

up to a similarity transformation matrix.

tridiagonal

butterfly

matrix of the special form

We have noted that the symplectic algorithm.

symmetric,

symplectic

by a trivial (i.e. symplectic

of interest

to develop

a canonical

form for butterfly matrices under similarity transformations by trivial matrices. We restrict our attention to unreduced symplectic butterfly matrices. since every butterfly unreduced

matrix can be decomposed

into two or more

smaller

ones.

THEOREM 1. Let B be an unreduced symplectic butterfly matrix. Then there exists a symplectic J-triangular matrix X such that B = X-‘l?X has the canonical

form

22

P. BENNER,

where D is a signature nal matrix.

matrix,

D is uniquely

transformation

H. FASSBENDER,

and T is a D-symmetric,

determined,

by a signature

AND

matrix,

T is determined

D. S. WATKINS

war-educed tridiagoup to a similarity

and X is unique up to multiplication

by

a signature

matrix of the form diag{C, C}. The eigenvalues of T are hi + Ai’, n, where hi, A,:‘, i = 1, . . . , n, are the eigenvalues of B. >.*.>

i=l

We are motivated by the decomposition

Proof. gular matrix

transformation

B,,

is used as a pivot to eliminate

(2>, in which the nonsin-

B,,. We now seek a similarity

that achieves a similar end. Let

[‘(r’ ;“I

Xc

(3)

be a trivial matrix. We shall determine conditions on Y and F under which the desired canonical form is realized. Focusing on the first block column of the similarity transformation

B = X-‘gX,

we have

We see that B,, = 0 if and only if YI?,, + FZ?,, = 0, which implies F = -Yl?,, k2<‘. Thus F is uniquely determined, once Y has been chosen. We have B,, inertia.

= Y-‘B -21Y-‘, Thus

the best

which shows that B,, we can do is take

B,,

and B,,

must have the same

= D = sign(&),

which

is

achieved by choosing Y = I I?,, ll”. In summary, we should take X as in (3), where

Y =

and 1~2111’ 2

F = -Ye,,&‘.

The resulting B has the desired form. The only aspect of the computation that is not completely straightforward is showing that B,, = -D. However, this becomes easy when one applies the following fact: If B is a symplectic The D-symmetry of T is also an matrix with B,, = 0, then B,, = -B&l. easy consequence of the symplectic structure of B. The uniqueness statements are easily verified.

1

is an eigenvector x [Y Ty = (A + A-‘)y. If

of

B with eigenvalue

A, then

y # 0, and ??

SR AND

HR EIGENVALUE

23

ALGORITHMS

REMARKS 1.

The

canonical

subdiagonal 2.

form

could

entries be positive:

The decomposition

be

made

unique

by insisting

t,, I, i > 0, i = 1,. . . , n -

(2) of the canonical

form

that

T’s

1.

B is

B=[::-:I =[‘o’ :I[:’21. 3. Theorem 1 is a theoretical result. From the standpoint of numerical stability, it might not be advisable to transform a symplectic butterfly matrix into

canonical

condensed

form.

In the

process,

the

spectral

information

into T as h + A-‘. The original information

A, A-’

can be recovered

is via

the inverse transformation

v

V

v+-+

However,

eigenvalues

not Lipschitz

near

continuous

2-

\ii

2

2

1

- 1.

f 1 will be resolved

poorly because

at v = * 2. The behavior

this map is

is similar to that of Van

Loan’s method for the Hamiltonian eigenvalue problem [22] (see also the remarks in Section 4.1). One may lose up to half of the significant digits as compared eigenvalues

to the standard

QR algorithm.

of the symplectic

1+6

S = G’

where

G is a randomly

For instance,

try to compute

the

matrix

I

generated

0

o

1 1+6

Givens

1 G,

rotation

and

6 is less than the

square root of the machine precision, once by applying the QR algorithm to S and once to S + SP1 followed by the inverse transformation given above. 4.

Since

= diag{TT, T}. Th us, forming T is equivalent to adding S + S + S-’ was used in similar fashion in 171 to compute the eigenvalues of a symplectic pencil.

we have B + B-’

BP’ to B. The transformation [IS,

“013

z(,_-/

+ 7-q - (,_a

+

a) = (,-fQ

- z)W -

Ed

[/&] put? [khliY-] g 30 slo$aama&a

pauF$qo

SNIBLVM

uay4 ‘(0 z h) h(r-~ + u) = 4~ 31 :,I 30 asoy aqJ lai\omJ 04 aIqpsod OS@ s;r 41 ‘8 30 asoy 111013 aye ~wp u~oys amy aM 3oold aq3 UI 0s

uw J 30 slopama&a

3 'a aNv

‘wraN38ssvd

'H ‘~3~~88

-d

PZ

SR AND

HR EIGENVALUE

THEOREM 2.

he an unreduced iteration of degree an HR iteration D-symmetric Pmof.

ALGORITHMS

Let

symplectic of degree

matrix in canonical form.

Then an SR

k with

shifis

pi + p,:‘.

i = 1,. . . , k,

to

on the

matrix T. The SR iteration

q is the Laurent

Notice that q(B)

has the form

B^=

= SR

S-‘BS,

polynomial

q(A) =

Since

butterfly

2k with shi$s pi, pi- ‘, i = 1, . . . , k, on B is equivalent

q(B)

where

25

= p(B

fI[P+

A--‘)

- (/-+ + /J-l)].

i=l

+ B-l>,

w h ere p is the ordinary

polynomial

B + B -’ = diag{T ‘, T}, we have

q(B) = An HR iteration

p(T’) P(T)

on T with shifts pi + I,-‘,

p(T)

= HU,

I/ is upper triangular, and H satisfies n . matrix, and T is D-symmetric.

= DHUD

(4)

i = 1, . . . , k, has the form

f = H-‘TH. H’DH

Now let us relate this to the SR iteration we have p(T’)

1.

= fi, where on B. Since

= H-“fiUD.

fi is a signature

Dp(T)

= p(T”)D,

26

P. BENNER,

H. FASSBENDER,

AND D. S. WATKINS

Thus

q(B) =

dTT) 1’(T)] =[H-T H][kJD

This is the SR decomposition

is symplectic,

u].

i of q(B),

for

and

is J-triangular.

Using this SR decomposition

to perform

the SR iteration,

we

obtain

z.

’

B^=S-‘BS= [ Thus the HR iteration on T is equivalent In principle

we can compute

the

i?

1

-D . 9

to the SR iteration on B. spectrum

of a symplectic

?? butterfly

matrix by putting it into canonical form, calculating the eigenvalues of T, then inverting the transformation A + A + A-l. Conversely, we can calculate the eigenvalues of a D-symmetric tridiagonal matrix T by embedding T and D in a symplectic matrix B, calculating the eigenvalues of B, and applying the transformation

A -+ A + A ‘.

These transformations are not necessarily advisable from the standpoint of numerical stability. The first will resolve eigenvalues near f 1 poorly because as we already mentioned, the inverse transformation is not Lipschitz-continuous. The second transformation is perhaps less objectionable. However, any eigenvalues of T that are near zero will have poor relative accuracy, because cancellation will occur in the transformation A + A + A-‘.

SR AND

HR EIGENVALUE

4.

HAMILTONIAN

THE

ALGORITHMS CASE

A matrix

A=

is called

A,,

A,,

A,,

A,,

I

1

if it satisfies (JA)“

Hamiltonian

Hamiltonian

if and only if A,,

Eigenvalues

of real, Hamiltonian

E

[wL,IX2,,

= ]A. One easily checks

= -AT,,

and

matrices

A,,

and

A,,

that

A is

are symmetric.

appear in plus-minus

pairs: If A is

an eigenvalue

of A with right eigenvector X, then -A is an eigenvalue of A with left eigenvector (IX) I’. If A is Hamiltonian and S is symplectic, then SiAS is Hamiltonian. Thus the Hamiltonian form is preserved urlder iterations

of the SR algorithm.

Hamiltonian

eigenvalue

problems

variety of continuous-time control, filtering, and estimation e.g., [12, 14, 16, 211 and the references given therein. A Hamiltonian

matrix is in J-tridiagonal

diagonal, and A,, such that S’AS

is tridiagonal. is J-tridiagonal

tonian J-tridiagonal An unreduced and A,,

4.1.

There

form if A,, , A,,,

exist numerous

preserves

are

matrices

S

the Hamil-

form. J&diagonal

is unreduced,

A Canonical

matrix is one for which

that is, its subdiagonal

entries

A,,

is nonsingular

are all nonzero.

Form for Hamiltonian J-ttidiagonal Matrice.r case, we now introduce a canonical

Just as we did in the symplectic for unreduced

see,

and A,,

symplectic

[6]. The SR algorithm

arise in a

problems:

J-tridiagonal

matrices

under

similarity

form

transformations

b>.

trivial matrices. THEOREM 3.

Let A be an unreduced

Then there exists a symplectic

J-triangular

Hamiltonian

J-tridiagonal

matr_+x.

matrix X such that A = X -‘A X

has the canonical form

where D is a signature matrix, and V is a symmetric, irreducible tridiagonal V is determined up to a similarity matrix. D is uniquely determined, transformation

by a signature

matrix,

and X is unique up to multiplication

by

a signature matrix of the form diag{C, Cl. Let T denote the D-symmetric matrix DV. The eigenvalues of T are hf , i = 1, . . . , n, where h,, - Ai, i = l,..., n, are the eigenvalues of A.

28

P. BENNER,

Just

Proof.

transforming

as in the

H. FASSBENDER,

proof

of Theorem

AND

1, we can

D. S. WATKINS

show

that

the

matrix

with

lA2y2

Y =

and

F = -YA,,&i’.

results in an A whose first block column is of the desired form. The fact that the other block column also has the desired form follows from the fact that A is Hamiltonian

and other elementary

As in the symplectic If

x [Y

considerations.

case, the uniqueness

1

is an eigenvector

of

statements

A with eigenvalue

are easily verified.

A, then

y # 0, and

Ty = DVy = h”y.

??

REMARKS. 1.

The canonical

or V’s subdiagonal 2.

form could be made unique by insisting that either T’s

entries be positive.

From the standpoint

to transform process,

a Hamiltonian

the spectral

eigenvalues extremely

of numerical

stability, it might not be advisable

J&diagonal

matrix into canonical

information

+h

of A are transformed vulnerable

to roundoff

is condensed

form. In the

into T as A2. Any small

to tiny eigenvalues

of T, which are then

errors in any subsequent

computations

on

T is tantamount

to

T. 3. squaring

We note that A. Squaring

A2 = diag{Tr, a Hamiltonian

the basis of Van Loan’s

T}. Thus,

forming

matrix to compute

square-reduced

method

its eigenvalues

is also

[22]. An error estimate

for

retrieving the eigenvalues of a Hamiltonian matrix from their squares computed by Van Loan’s method is given in [22] and indicates that one may lose up to half the significant method

digits as compared

as the QR algorithm.

to a numerically

The same limitations

transform a Hamiltonian J-tridiagonal pute its eigenvalues via those of T.

backward stable

in accuracy

matrix into canonical

apply if we

form and com-

HR EIGENVALUE

SR AND 4.

As in the symplectic

case, the eigenvectors

from those of T. If Ty = A2y

are eigenvectors

4.2.

29

ALGORITHMS

(y

of A associated

of A can be recovered

z O), then

with eigenvalues

k A.

Equivalence of the HR and Hnmiltonian SR Algorithms an SR iteration on a Hamilton matrix A. Since the eigenvalues

Consider

occur in plus-minus

pairs, it is reasonable

pairs. If we wish to effect

to choose

an SR iteration

the shifts in plus-minus

of degree

2.k with shifts

+ F!,

>..., k, we use the polynomial

i=l

q(A)

Again we restrict insist that complex

= h(A i=l

-

ourselves

to the nonsingular

p,l)(A

+ p,l)

= AA2 i=

shifts be present

A=

in conjugate

v

(1

[D

0

-

/$I.

I

case for simplicity.

LYc also

pairs, so that q( A) is ~~11.

1

he an unreduced Hamiltonian]-tridiagonal matrix in canonical form. Then at1 SR iteration of degree 2k with shifts + pi, i = 1, . . . , k, OILA is eyuicalrnt to an HR iteration of degree k with shifts &, i = 1, , . , k. on the D-.ynr~~tric~ matrix T = DV. Proof.

The SR iteration

q(A) where

has the form = SR,

i

= SF’AS,

q is the polynomial

q(A) =

I?(^ - /-%)(A+ I-%).

i=l

30

P. BENNER,

Notice that q(A)

Since

AND

D. S. WATKINS

p(T)},

in analogy with

= p( A’), where

A2 = diag{Tr,

(4). Proceeding

H. FASSBENDER,

T}, we have q(A)

= diag{ p(Tr>,

as in the proof of Theorem

on T with shifts pCLf,i = 1,.

2, we recall that an HR iteration

. . , k, has the form

rE’= H-ITH.

p(T) = HU,

U is upper triangular, H satisfies HTDH = 6, where 6 is a signature 1 1 matrix, and T is D-symmetric. Just as in the proof of Theorem 2, we have

p(TT)

= DHUD = H-TLkJD,

40) = which

dTT) p(Tj] [“’

I

z

is an SR decomposition

perform

the SR iteration,

*

A

LI

so

of q(A).

Using

H][fim u], this

SR decomposition

to

we obtain

where T = DV. Thus the HR iteration on A.

on T is equivalent

to the SR iteration ??

In principle we can compute the spectrum of a Hamiltonian, J-tridiagonal matrix by putting it into canonical form, calculating the eigenvalues of

T = DV, then taking square roots. We have already noted the dangers of this approach. Conversely, we can calculate the eigenvalues of a D-symmetric tridiagonal matrix T by embedding V = DT and D in a Hamiltonian J-tridiagonal matrix A in canonical form, calculating the eigenalues of A, and squaring them.

HR

SR AND

5.

EIGENVALUE

31

ALGORITHMS

CONCLUSIONS

We have derived connections J-tridiagonal

matrices.

I-tridiagonal

matrices

for symplectic

Transforming

syrnplectic

into the canonical

forms

it can be shown that the SR iterations special

choice

matrix

of half

eigenvahles

of shifts

are equivalent

the

Using

size.

and eigenvectors

the HR iteration

between

and the SR algorithms

ric matrices

this approach

so obtained

Hamiltonian

matrices

3 and 3, with

a

on a sign-symmetric

it is possible

butterfly

agonal matrices by applying the HR algorithm

and

in Sections

to an HR iteration

of symplectic

and Hamiltonian

butterfly

introduced

for the

for sign-symmet-

butterfly

to obtain

and Hamiltonian

to the associated

thr

J-trid-

sign-symmet-

ric matrix. The suffer

results

from

are

mainly

a possible

tion to carionical

of theoretical

interest,

loss of half the significant

as the digits

resulting

during

methods

the transforma-

form.

KEFERENCES 1

2

3

4

3 6

7 8 9 IO

<:. Banse, Symplektische Eigenwcxrtvrrfahren zur Liisung Zritdiskretrr optimalrr Steucrungsprohleme, Ph.D. Thesis, Univ. of Bremrn, 199rj. G. Bansr and A. Bunse-Gerstnrsr, A condensed form tar the solution of the, ‘~‘hroq s!-mplectic eigenvalue problem, in Sy.sterns and Networks: Mnthonatical nrd Applications (U. Helmke, K. Menniken, and J. Saut~r, Eds.), Akademic> \‘erlag, 1994, pp. 613-616. P. Bcnner and H. F&bender, The syrnplectic eigenvalue problem, the hutterfl!. form, the SR algorithm, and the I,anczos method, Linccrr Algcbm Appl., to qq’e:lr. M. A. Brebner and J. Grad, Eigenvalues of As = ABs for real sylnmrtric matrices A and B computed t)y reduction to pselltlosylrlmetric kxm and tlw HH prowss, Linear Algebra Appl. 43:99-l 18 (1982). A. Blmsr-Gerstner, An analysis of the HR algorithm for colnputing the r+gcanvalIWS of a trlatrix, Linear Algebrcl Appl. 3.5:155-173 (1981). A. Hnnsr Gerstner and V. Mehrrnann, A s\;lllplectic OR-like algorithm for the, solution of the real algebraic Riccati equation. IEEE Trrrrl.~. AILtowwt. (,‘ontro/ AC-:31:1104-1113 (lYs6). J. Della-l)ora, Nrunerical linear algorithms and grollp th(sor?,. Litwtrr .4/p1~rx App/. 10:267-283 (19X5). L. Elsner, On some algebraic prohlems in connection with general eigenvaluc algorithms, Linear Algebra A&. 26:123-138 (1979). J. (:. F. Francis, The QR tranaformatiolt, parts I. II, Cornput. /. -l:%E-272, 3.3~345 (1961). I). Hinrichsen and Iv. K. Son, Stability radii of’linear discwte-time systems anal s~mplectic pencils, Interr~cd. /. ~~ohf~st Xonlincwr (:orltro/ I :79-97 (199 I).

32

P. BENNER,

11 V. N. Kublanovskaya, eigenvalue

On

problem,

12 P. Lancaster 13 C. Lanczos,

some

USSR

An iteration

algorithms

Comput.

and L. Rodman,

linear differential

H. FASSBENDER,

AND

for the

D. S. WATKINS

solution

Math. and Math.

of the

complete

Phys. 3:637-657

(1961).

The AZgebruic Riccati Equation,

method

for the solution

and integral operators,

Oxford U.P., 1995.

of the eigenvalue

I. Res. Nut. Bur.

problem

Stundurds

of

45:255-281

(1950).

14 A. J. Laub,

Invariant

equations, Eds.),

subspace

in The Riccuti

Springer-Verlag,

15 W.-W.

Lin,

discrete-time

methods

Equation

for the

(J. Bittanti,

numerical A. J. Laub,

solution

of Riccati

and J. C. Willems,

pp. 163-196.

A new

method

algebraic

for computing

Riccati

equation,

the

Linear

closed

Algebra

loop

eigenvalues

of a

(1987). Theory and

AppZ. 6:157-180

16 V. Mehrmann, The Autonomous Linear Quadratic Control Problem, NumericuZ Solution, Lecture Notes in Control and Inform. Sci. 163, SpringerVerlag,

Heidelberg,

July 1991.

17 R. V. Patel, On computing the eigenvalues of a symplectic pencil, Linear Algebra AppZ. 188/189:591-611 (1993); see also Proceedings of CDC-31, Tuscan, Ariz., 1992, pp. 1921-1926.

18 H. Rutishauser, Phys. 5:233-251 19 H. Rutishauser, Math.

Der

7, Birkhauser,

20 H. Rutishauser, Nat. Bur.

21 V. Sima, Marcel

Der

2.

Anger.

Math.

Inst.

Angew.

Solution

Algorithms

Mitt.

Quotienten-DifSerenzen-Algorithmus, Base],

Standards

Dekker,

Quotienten-Differenzen-Algorithmus,

(1954). 1957. of eigenvalue

AppZ. Math.

problems

Ser. 49:47-81

for Linear-Qcmdrutic

New York,

with

the

LR-transformation,

(1958).

Optimization,

Pure Appl. Math.

200,

1996.

22 C. F. Van Loan, A symplectic method for approximating all the eigenvalues of a Hamiltonian matrix, Linear Algebra AppZ. 16:233-251 (1984). of the QR algorithm, SIAM Rev. 24:427-440 23 D. S. Watkins, Understanding (1982). 24 D. S. Watkins, Fnndumentuls of Matrix Computations, Wiley, New York, 1991. 25 D. S. Watkins and L. Elsner, Convergence of algorithms of decomposition type for the eigenvalue

26 J. H. Wilkinson,

problem,

Linear

The Algebraic

Algebra

Eigenvulue

AppZ. 143: 19-47 Problem,

Clarendon,

Received 2.5 April 1997;final manuscript accepted 15 ]une 1997

(1991). Oxford,

1965.