Document [original]

Technische Universit¨at Berlin

Institut f¨ur Mathematik

Skew-Hamiltonian and Hamiltonian

Eigenvalue Problems: Theory, Algorithms

and Applications

Peter Benner, Daniel Kressner, and

Volker Mehrmann

Preprint 44-2003

Preprint-Reihe des Instituts f¨ur Mathematik

Technische Universit¨at Berlin

http://www.math.tu-berlin.de/preprints

Report 44-2003 2003

SKEW-HAMILTONIAN AND HAMILTONIAN

EIGENVALUE PROBLEMS: THEORY, ALGORITHMS

AND APPLICATIONS∗

Peter Benner

Technische Universit¨at Chemnitz

Fakult¨at f¨ur Mathematik

b[email protected]

Daniel Kressner

Technische Universit¨at Berlin

Institut f¨ur Mathematik

[email protected]erlin.de

Volker Mehrmann

Technische Universit¨at Berlin

Institut f¨ur Mathematik

[email protected]erlin.de

Abstract Skew-Hamiltonian and Hamiltonian eigenvalue problems arise from a

number of applications, particularly in systems and control theory. The

preservation of the underlying matrix structures often plays an impor-

tant role in these applications and may lead to more accurate and more

efficient computational methods. We will discuss the relation of struc-

tured and unstructured condition numbers for these problems as well

as algorithms exploiting the given matrix structures. Applications of

Hamiltonian and skew-Hamiltonian eigenproblems are briefly described.

Keywords: Hamiltonian matrix, skew-Hamiltonian matrix, structured condition num-

bers, structure-preserving algorithms.

∗Supported by the DFG Research Center “Mathematics for key technologies” (FZT 86) in

Berlin.

Contents

Skew-Hamiltonian and Hamiltonian Eigenvalue Problems: Theory, Algo-

rithms and Applications 5

Peter Benner,Daniel Kressner,Volker Mehrmann

1 Preliminaries 10

2 The Skew-Hamiltonian Eigenvalue Problem 12

2.1 Structured Decompositions 12

2.2 Structured Condition Numbers 14

2.3 Algorithms 21

3 The Hamiltonian Eigenvalue Problem 23

3.1 Structured Decompositions 23

3.2 Structured Condition Numbers 24

3.3 Algorithms 28

4 Applications 36

4.1 Stability Radius Computation 36

4.2 H∞Norm Computation 37

4.3 Algebraic Riccati Equations 37

4.4 Quadratic Eigenvalue Problems 38

4.5 Other Applications 39

5 Concluding Remarks 39

References 41

Introduction

Computing eigenvalues and invariant subspaces of matrices with struc-

ture has been an active field of research during the last two decades.

In many instances it has been shown that the exploitation of matrix

structures may give rise to more accurate and more efficient numerical

methods. In this paper we will discuss this issue for two classes of matri-

ces, skew-Hamiltonian and Hamiltonian matrices. A skew-Hamiltonian

matrix has the form

W=·A G

Q AT¸, G =−GT, Q =−QT,(1)

while a Hamiltonian matrix reads as

H=·A G

Q−AT¸, G =GT, Q =QT,(2)

where A, G and Qare real n×nmatrices. A number of applications from

control theory and related areas lead to eigenvalue problems involving

such matrices, with a stronger emphasis on Hamiltonian matrices, see

Section 4.

One of the first questions one should always ask when dealing with

structured eigenvalue problems is what kind of advantages can princi-

pally be expected from exploiting structures. With respect to accuracy

of computed eigenvalues and invariant subspaces this question leads to

the notion of structured condition numbers and their relationship to un-

structured ones. It is interesting to note that the two matrix structures

under consideration differ significantly in this aspect. While it is abso-

lutely necessary to use a structure-preserving algorithm for computing

invariant subspaces of skew-Hamiltonian matrices, the merits of struc-

ture preservation for Hamiltonian matrices are of a more subtle nature

and not always relevant in applications. If one is interested in efficiency

then there is not so much that can be expected. Both matrix classes de-

pend on 2n2+O(n) parameters compared to 4n2parameters of a general

2n×2nmatrix. Hence, a structure-preserving algorithm can be expected

to be at best a decent factor faster than a general-purpose method; for

the matrix classes considered here, this factor is usually in the range of

2–3, see [Benner et al., 2000; Benner and Kressner, 2004; Benner et al.,

1998].

Another important question is whether it is actually possible to de-

sign an algorithm capable to achieve the possible advantages mentioned

above. An ideal method tailored to the matrix structure would

be strongly backward stable in the sense of Bunch described in

[Bunch, 1987], i.e., the computed solution is the exact solution

corresponding to a nearby matrix with the same structure;

be reliable, i.e., capable to solve all eigenvalue problems in the

considered matrix class; and

require O(n3) floating point operations (flops), preferably less than

a competitive general-purpose method.

While for skew-Hamiltonian matrices such a method is known [Van Loan,

1984b], it has been a long-standing open problem to develop an ideal

method for the Hamiltonian eigenvalue problem. So far there is no

method known that meets all three requirements satisfactorily.

The main purpose of this paper is to survey theory and algorithms for

(skew-)Hamiltonian eigenvalue problems. With respect to algorithms,

the account will necessarily be rather incomplete, simply because of the

vast number of algorithms that have been developed. Instead, our focus

will be on methods that are based on orthogonal transformations and

suitable for dense, small to medium-sized matrices. Nevertheless, they

will be related to other existing methods. Another goal in this work

is to describe applications of (skew-)Hamiltonian eigenvalue problems

and identify the extent to which a structure-preserving algorithm may

help to address these applications in a more accurate or more efficient

manner.

The structure of this survey is as follows. After having introduced

some notation and preliminary material in the first section we devote

the second section to the skew-Hamiltonian eigenvalue problem. We re-

view structured Hessenberg-like, Schur-like and block diagonal decom-

positions. This is followed by some recent and new results on struc-

tured condition numbers for the eigenvalues and invariant subspaces.

The section is concluded by a description of the ideal method for skew-

Hamiltonian matrices that was mentioned above. Section 3 contains

similar results for the Hamiltonian eigenvalue problem, with a more ex-

tensive treatment of structure-preserving algorithms. In particular, we

present an explicit version of the Hamiltonian QR algorithm, describe

an alternative derivation for the method given in [Benner et al., 1998],

via an embedding in skew-Hamiltonian matrices, and give an example of

an iterative refinement algorithm. Some applications related to systems

and control theory and how they may benefit from the use of structure-

preserving algorithms are the subject of Section 4.

This paper is accompanied by a Matlab software library for solving

skew-Hamiltonian and Hamiltonian eigenvalue problems. The library is

based on recently developed Fortran 77 routines [Benner and Kressner,

2004] and described in [Kressner, 2003b], which also contains numerical

examples illustrating some of the aspects in this survey.

1. Preliminaries

An ubiquitous matrix in this work is the skew-symmetric matrix

J2n=·0In

−In0¸,(3)

where Indenotes the n×nidentity matrix. In the following we will drop

the subscripts nand 2nwhenever the dimension of the corresponding

matrix is clear from its context. By straightforward algebraic manipula-

tion one can show that a Hamiltonian matrix His equivalently defined by

the property HJ = (HJ)T. Likewise, a matrix Wis skew-Hamiltonian

if and only if WJ =−(WJ)T. Any matrix S∈R2n×2nsatisfying

STJS =SJST=Jis called symplectic, and since

(S−1HS)J=S−1HJS−T=S−1JTHTS−T= [(S−1HS)J]T

we see that symplectic equivalence transformations preserve Hamiltonian

structures. There are cases, however, where both Hand S−1HS are

Hamiltonian but Sis not a symplectic matrix [Freiling et al., 2002]. In

a similar fashion the same can be shown for skew-Hamiltonian matrices.

From a numerical point of view it is desirable that a symplectic ma-

trix U∈R2n×2nis also orthogonal. Such a matrix is called orthogonal

symplectic; the two relations UTJU =Jand UTU=Iimply JUJT=U

which effectively means that every orthogonal symplectic matrix Uhas

the block structure

U=·U1U2

−U2U1¸, U1, U2∈Rn×n.

Two types of elementary orthogonal matrices have this form. These are

2n×2nGivens rotation matrices of the type

Gj(θ) = 





Ij−1

cos θsin θ

In−1

−sin θcos θ

In−j







,1≤j≤n,

for some angle θ∈[−π/2, π/2) and the direct sum of two identical n×n

Householder matrices

(Hj⊕Hj)(v, β) = ·In−βvvT

In−βvvT¸,

where vis a vector of length nwith its first j−1 elements equal to zero

and βa scalar satisfying β(βvTv−2) = 0. Here, ‘⊕’ denotes the direct

sum of matrices.

A simple combination of these transformations can be used to map

an arbitrary vector x∈R2ninto the linear space

Ej= span{e1,...,ej, en+1,...,en+j−1},

where eiis the ith unit vector of length 2n. Such mappings form the

backbone of virtually all structure-preserving algorithms based on or-

thogonal symplectic transformations. They can be constructed using the

following algorithm, where it should be noted that elements 1, . . . , j −1

and n+ 1,...,n+j−1 of the vector xremain unaffected.

Algorithm 1

Input: A vector x∈R2nand an index j≤n.

Output: Vectors v, w ∈Rnand β, γ, θ ∈Rso that

[(Hj⊕Hj)(v, β)·Gj(θ)·(Hj⊕Hj)(w, γ)]Tx∈ Ej.

1 Determine v∈Rnand β∈Rsuch that the last n−jelements of

x←(Hj⊕Hj)(v, β)xare zero, see [Golub and Van Loan, 1996,

p.209].

2 Determine θ∈[−π/2, π/2) such that the (n+j)th element of

x←Gj(θ)xis zero, see [Golub and Van Loan, 1996, p.215].

3 Determine w∈Rnand γ∈Rsuch that the (j+ 1)th to the nth

elements of x←(Hj⊕Hj)(w, γ)xare zero.

The three steps of this algorithm are illustrated in Figure 1. Orthogonal

symplectic matrices of the form

Ej(x)≡Ej(v, w, β, γ, θ) := (Hj⊕Hj)(v, β)·Gj(θ)·(Hj⊕Hj)(w, γ),(4)

as computed by Algorithm 1, will be called elementary.

Let F=h0

0i, then we obtain the following variant of elementary

orthogonal symplectic matrices:

[F·Ej(Fx)·F]Tx∈span{e1,...,ej−1, en+1, . . . , en+j}.

For the sake of brevity we set En+j(x) := F·Ej(Fx)·F, whenever

1≤j≤n.

Figure 1. The three steps of Algorithm 1 for n= 4 and j= 2.

2. The Skew-Hamiltonian Eigenvalue Problem

Imposing skew-Hamiltonian structure on a matrix Whas a number

of consequences for the eigenvalues and eigenvectors of W; one is that

every eigenvalue has even algebraic multiplicity and hence appears at

least twice. An easy way to access all these spectral properties is to ob-

serve that for any skew-Hamiltonian matrix Wthere exists a symplectic

matrix Sso that

S−1WS =·W11 0

0WT

11 ¸.(5)

This decomposition – among others – will be described in the following

section.

2.1 Structured Decompositions

As a first application of elementary matrices we obtain Algorithm 2

below which constructs the following structured Hessenberg-like form:

given a skew-Hamiltonian matrix W∈R2n×2nthere is always an orthog-

onal symplectic matrix Uso that UTWU has Paige/Van Loan (PVL)

form, i.e.,

UTWU =·W11 W12

0WT

11 ¸=



@

,(6)

where W11 ∈Rn×nis an upper Hessenberg matrix [Van Loan, 1984b].

Algorithm 2 (PVL decomposition [Van Loan, 1984b])

Input: A skew-Hamiltonian matrix W∈R2n×2n.

Output: An orthogonal symplectic matrix U∈R2n×2n;Wis overwritten

with UTW U having the form (6).

E2E2

E3E3

Figure 2. Illustration of two loops of Algorithm 2 for n= 4.

U←I2n.

for j←1,2,...,n−1

Set x←Wej.

Apply Algorithm 1 to compute Ej+1(x).

Update W←Ej+1(x)TW Ej+1(x), U←UEj+1(x).

end for

A proper implementation of this algorithm requires 40

3n3+O(n2) flops

for reducing Wand additionally 16

3n3+O(n2) flops for computing U.

Figure 2 illustrates Algorithm 2 for n= 4.

An immediate consequence of the PVL form (6) is that each eigenvalue

of Whas even algebraic multiplicity. The same is true for the geometric

multiplicities. To see this we need to eliminate the skew-symmetric

off-diagonal block W12, for which we can use solutions of the following

singular Sylvester equation.

Proposition 3 The matrix equation

W11R−RW T

11 =−W12 (7)

is solvable for all skew-symmetric matrices W12 ∈Rn×nif and only if

W11 ∈Rn×nis nonderogatory, i.e., each eigenvalue of W11 has geo-

metric multiplicity one. In this case, any solution Rof (7) is real and

symmetric.

Proof. This result can be found in [Gantmacher, 1960; Faßbender

et al., 1999]. Actually, the second part is not explicitely stated there but

follows easily from the proof of Proposition 5 in [Faßbender et al., 1999].

We now use this proposition to block-diagonalize a skew-Hamiltonian

matrix in PVL form (6) assuming that W11 is nonderogatory. For this

purpose let Rbe a solution of (7), then the symmetry of Rimplies that

Iiis symplectic. Applying the corresponding symplectic equivalence

transformation yields the transformed matrix

·I R

0I¸−1·W11 W12

0WT

11 ¸· I R

0I¸=·W11 0

0WT

11 ¸.(8)

Note that there is a lot of freedom in the choice of Ras equation (7)

admits infinitely many solutions. From a numerical point of view the

matrix Rshould be chosen so that its norm is as small as possible. The

same question arises in the context of structured condition numbers and

will be discussed in the next section.

It should be stressed that assuming W11 to be nonderogatory is not

necessary and thus, the even geometric multiplicity of eigenvalues also

holds in the general case. In fact, in [Faßbender et al., 1999] it is shown

that any skew-Hamiltonian matrix can be reduced to block diagonal

form (8) using symplectic equivalence transformations. The proof, how-

ever, is much more involved than the simple derivation given above.

Another way to go from a skew-Hamiltonian matrix Win PVL form (6)

to a more condensed form is to reduce W11 further to real Schur form.

This can be achieved by constructing an orthogonal matrix Q1so that

T=QT

1W11Q1is in real Schur form [Golub and Van Loan, 1996,

Thm.7.4.1]:

T=





T11 T12 ··· T1m

0T22 ....

.......Tm−1,m

0··· 0Tmm







,(9)

where all diagonal blocks Tjj of Tare of order one or two. Each scalar

diagonal block contains a real eigenvalue and each two-by-two diagonal

block contains a pair of conjugate complex eigenvalues of W11. Setting

U=U(Q1⊕Q1), we obtain a skew-Hamiltonian Schur decomposition of

UTW˜

U=·T˜

0TT¸,(10)

where ˜

G=QT

1W12Q1is skew-symmetric.

2.2 Structured Condition Numbers

In this section we investigate the change of eigenvalues and certain

invariant subspaces of a skew-Hamiltonian matrix Wunder a sufficiently

small, skew-Hamiltonian perturbation E. Requiring the perturbation to

be structured as well may have a strong positive impact on the sensitivity

of the skew-Hamiltonian eigenvalue problem; this is demonstrated by the

following example.

Example 4 Consider the parameter-dependent matrix

W(ε1, ε2) = 





1 0 0 0

0 2 0 0

ε1ε21 0

−ε20 0 2







The vector e1= [1,0,0,0]Tis an eigenvector of W(0,0) associated with

the eigenvalue λ= 1. No matter how small ε1>0is, any eigen-

vector of W(ε1,0) associated with λhas the completely different form

[0,0,α,0]Tfor some α6= 0. On the other hand, the skew-Hamiltonian

matrix W(0, ε2)has an eigenvector [1,0,0, ε2]Trather close to e1.

Before we deal with structured perturbations, we briefly review stan-

dard perturbation results that apply to general matrices and perturba-

tions, for further details, see e.g. [Stewart and Sun, 1990; Sun, 1998]. Let

A∈Rn×nand let X ⊂ Rnbe a k-dimensional (right) invariant subspace

of A, i.e., AX ⊆ X. If the columns of Xand X⊥span orthonormal bases

for Xand X⊥, respectively, then we obtain a block Schur decomposition

£X X⊥¤TA£X X⊥¤=·A11 A12

0A22 ¸,(11)

where A11 ∈Rk×kand A22 ∈R(n−k)×(n−k). The block A11 satisfies the

relation AX =XA11, which implies that A11 is the representation of A

with respect to X. An important operator associated with the decom-

position (11) is the linear matrix operator T:R(n−k)×k7→ R(n−k)×k

with

T:R7→ A22R−RA11.(12)

One can show that this Sylvester operator Tis invertible if and only if

A11 and A22 have no eigenvalue in common [Golub and Van Loan, 1996,

pp.366–369]. If this condition holds then Xis called a simple invari-

ant subspace. We are now ready to formulate a perturbation expansion

theorem for invariant subspaces and their representations as it can be

found, e.g., in [Sun, 1998, Sec.2.1.2].

In the following we denote by k·k2the Euclidean norm and the asso-

ciated spectral norm for matrices, and by k·kFthe Frobenius norm.

Theorem 5 Let A∈Rn×nhave a block Schur decomposition of the

form (11) and assume that the invariant subspace Xspanned by the

columns of Xis simple. Let E∈Cn×nbe a perturbation of sufficiently

small norm. Then, there is an invariant subspace ˆ

X= span ˆ

Xof A+E

with representation ˆ

A11 satisfying the expansions

A11 =A11 + (YHX)−1YHEX +O(kEk2

2),(13)

X=X−X⊥T−1XH

⊥EX +O(kEk2

2),(14)

where Tis as in (12), the columns of Yform an orthonormal basis for

the left invariant subspace of Abelonging to the eigenvalues of A11 and

XH(ˆ

X−X) = 0.

Bounding the effects of Ein the expansions (13) and (14) can be used

to derive condition numbers for eigenvalues and invariant subspaces. For

example, let λbe a simple eigenvalue of Awith right and left eigenvectors

xand y, respectively. Then Theorem 5 with A11 = [λ] and ˆ

A11 = [ˆ

λ]

yields

|ˆ

λ−λ| ≤ kxk2·kyk2

|yHx|kEk2+O(kEk2

2).

This inequality is attained up to first order by E=εxyHfor any ε > 0,

which shows that the absolute condition number of a simple eigenvalue

λcan be written as

c(λ) := lim

ε→0sup

kEk2≤ε

|ˆ

λ−λ|

ε=kxk2·kyk2

|yHx|.(15)

Note that c(λ) is independent of the choice of xand y.

For a simple invariant subspace Xspanned by the columns of Xwe

obtain

kˆ

X−XkF≤ kT−1k·kEkF+O(kEk2

F),(16)

where kT−1kis the norm induced by the Frobenius norm:

kT−1k:= sup

S6=0

kT−1(S)kF

kSkF

=µinf

R6=0 kA22R−RA11kF

kRkF¶−1

Again, inequality (16) can be attained up to first order by choosing

E=εX⊥V XHwith ε > 0, and a matrix V∈R(n−k)×kwith kVkF= 1

satisfying kT−1(V)kF=kT−1k. Turning (16) into a condition number

for an invariant subspace further requires relating kˆ

X−XkFto quantities

that are independent of the choice of bases for Xand ˆ

X. The matrix

Θ(X,ˆ

X) = diag(θ1(X,ˆ

X), θ2(X,ˆ

X),...,θk(X,ˆ

X)),

where θi(X,ˆ

X) are the canonical angles between Xand ˆ

X[Stewart and

Sun, 1990, p.43], is such a quantity. One can show that XH(ˆ

X−X) = 0

implies

kΘ(X,ˆ

X)kF=kˆ

X−XkF+O(kˆ

X−Xk3

F).

Hence, we obtain the following condition number for an invariant sub-

space X:

c(X) := lim

ε→0sup

kEkF≤ε

kΘ(X,ˆ

X)kF

ε=kT−1k.

Note that kT−1kis invariant under an orthonormal change of basis for

X. A direct (albeit expensive) way to compute this quantity is to express

Tin terms of Kronecker products:

vec(T(R)) = KT·vec(R), KT=Ik⊗A22 −AT

11 ⊗In−k,

with the Kronecker product ‘⊗’ of two matrices [Golub and Van Loan,

1996, Sec. 4.5.5] and the operator vec which stacks the columns of a

matrix into one long vector in their natural order. Then kT−1k−1is

the minimal singular value of the k(n−k)×k(n−k) matrix KT. In

practice, one estimates T−1by solving a few Sylvester equations with

particularly chosen right hand sides, see e.g. [Higham, 1996].

Structured condition numbers for eigenvalues We now turn to

the perturbation theory for an eigenvalue λof a matrix Wunder a

perturbation E, where both Wand Eare skew-Hamiltonian. As λ

is necessarily a multiple eigenvalue we cannot apply Theorem 5 to λ

alone but must consider the eigenvalue cluster containing all copies of λ.

Assuming that λhas algebraic multiplicity two, there are two linearly

independent eigenvectors x1and x2corresponding to λ. Let [x1, x2] =

XR be a QR decomposition with X∈C2n×2and R∈C2×2, then

WX =W[x1, x2]R−1= [x1, x2]A11R−1= [x1, x2]R−1A11 =XA11,

where A11 = diag(λ, λ). An analogous relation holds for the two eigen-

vectors ˆx1,ˆx2belonging to the eigenvalue ˆ

λof the perturbed matrix

W+E. As the spectral norm of ˆ

A11 −A11 is given by |ˆ

λ−λ|, Theorem 5

implies that

|ˆ

λ−λ|=k(¯

XHJX)−1¯

XHJEXk2+O(kEk2

≤ k(¯

XHJX)−1k2+O(kEk2

2),(17)

where we also used the fact that the columns of J¯

Xspan the two-

dimensional left invariant subspace belonging to λ. Note that ¯

Xdenotes

the complex conjugate of X. For real λwe may assume X∈R2n×2and

use the skew-Hamiltonian matrix E=εJT

2nXJ2XTto show that inequal-

ity (17) can be attained up to first order by a skew-Hamiltonian pertur-

bation. This implies that the structured eigenvalue condition number for

an eigenvalue λ∈Rof a skew-Hamiltonian matrix satisfies

cW(λ) := lim

ε→0sup

kEk2≤ε

Eskew-Hamiltonian

|ˆ

λ−λ|

ε=k(XHJX)−1k2.

Likewise for complex λ, we can use perturbations of the form E=

εJT

2n¯

XJ2XH. Note that Esatisfies (EJ)T=−(EJ), i.e., Emay be

regarded as a complex skew-Hamiltonian matrix. It is an open prob-

lem whether one can construct a real skew-Hamiltonian perturbation to

show cW(λ) = k(XHJX)−1k2for complex eigenvalues.

By straightforward computation one can obtain a simple expression

for cW(or an upper bound thereof if λ∈C) in terms of the eigenvectors

x1, x2belonging to λ:

k(XHJX)−1k2=1

|¯xH

1Jx2|qkx1k2

2·kx2k2

2−|xH

1x2|2.

Note that this happens to be the unstructured condition number of the

mean of the eigenvalue cluster containing both copies of λ, see [Bai et al.,

1993] and the references therein.

Structured condition numbers for invariant subspaces The in-

variant subspaces of a skew-Hamiltonian matrix Wthat are usually of

interest in applications are those which are isotropic.

Definition 6 A subspace X ⊆ R2nis called isotropic if X ⊥ J2nX. A

maximal isotropic subspace is called Lagrangian.

Obviously, any eigenvector of Wspans an isotropic invariant subspace

but also the first k≤ncolumns of the matrix ˜

Uin a skew-Hamiltonian

Schur decomposition (10) share this property. Roughly speaking, an

invariant subspace Xof Wis isotropic if Xcorresponds to at most one

copy of each eigenvalue of W. Necessarily, Xis not simple, which makes

the application of Theorem 5 impossible.

Instead, it was shown in [Kressner, 2003d] how one can adapt a tech-

nique developed by Stewart in [Stewart, 1971; Stewart, 1973] in order

to obtain perturbation bounds for X. Here, we describe this approach

only for the important special case when Xhas maximal dimension n,

i.e., Xis Lagrangian. Let the columns of Xform an orthonormal basis

for X. Then [X, JX] is orthogonal and we have the skew-Hamiltonian

block Schur decomposition

£X JX ¤TW£X JX ¤=·W11 W12

0WT

11 ¸.(18)

If Wis perturbed by a skew-Hamiltonian matrix Ethen

£X JX ¤T(W+E)£X JX ¤=·˜

W11 ˜

W12

W21 ˜

11 ¸(19)

where ˜

W12,˜

W21 are skew-symmetric and k˜

W21kF≤ kEkF. Any matrix

X∈R2n×nwith orthonormal columns and ˆ

XTX6= 0 can be written as

X= (X+JXR)(I+RTR)−1/2(20)

for some matrix R∈Rk×k[Stewart, 1973]. The columns of ˆ

Xspan

an invariant subspace ˆ

Xof W+Eif and only if Rsatisfies the matrix

equation

T˜

W(R) + Φ(R) = ˜

W21 (21)

with the Sylvester operator T˜

W:R7→ R˜

W11 −˜

11Rand the quadratic

matrix operator Φ : R7→ R˜

W12R. Moreover, Xis isotropic if and

only if Ris symmetric. The solution of (21) is complicated by the fact

that the dominating linear operator T˜

Wis singular. However, if ˜

W11

is nonderogatory, then Proposition 3 shows that the restricted operator

T˜

W: symm(n)7→ skew(n), where symm(n){skew(n)}denotes the set

of {skew-}symmetric n×nmatrices, is onto. This allows us to define an

operator T†

W: skew(n)7→ sym(n), which maps a skew-symmetric matrix

Q∈Rn×nto the minimal Frobenius-norm solution of R˜

W11−˜

11R=Q,

which must be symmetric according to Proposition 3. The norm of T†

induced by the Frobenius-norm is given by

kT†

Wk:= sup

Q6=0

Q∈skew(n)

kT†

W(Q)kF

kQkF

.(22)

This can be used to estimate the norm of a solution of the nonlinear

equation (21).

Theorem 7 Let the matrices ˜

Wij be defined by (19) and assume that

W11 is nonderogatory. If 4kT†

Wk2·k ˜

W12kF·k ˜

W21kF<1with kT†

Wkas

in (22), then there is a symmetric solution Rof (21) satisfying

kRkF≤2kT†

Wk·k ˜

W12kF.

Proof. This result can be proven along the lines of the proof of The-

orem 2.11 in [Stewart and Sun, 1990] by constructing a sequence

R0= 0, Ri+1 =T†

W(˜

W21 −Φ(Ri))

and applying the contraction mapping theorem [Ortega and Rheinboldt,

1970] to this sequence.

Using the fact that the tangents of the canonical angles between X

and ˆ

Xare the singular values of R[Stewart and Sun, 1990, p.232], The-

orem 7 implies the following perturbation bound for isotropic invariant

subspaces.

Corollary 8 Let W, E ∈R2n×2nbe skew-Hamiltonian matrices, and let

the columns of X∈R2n×nform an orthonormal basis for an isotropic

invariant subspace Xof W. Assume that ˜

W11 =XT(W+E)Xis non-

derogatory and that 4kT†

Wk2·kW+EkF·kEkF<1, with kT†

Wkdefined

as in (22). Then there is an isotropic invariant subspace ˆ

Xof W+E

satisfying

ktan[Θ(X,ˆ

X)]kF≤αkT†

Wk·kEkF,(23)

where α≤2.

It should be remarked that the factor αin (23) can be made arbitrar-

ily close to one if we let kEkF→0. Furthermore, (23) still holds in

an approximate sense if the operator T˜

Wis replaced by TW:R7→

RW11 −WT

11Rcorresponding to the unperturbed block Schur decom-

position (18). This shows that the structured condition number for an

isotropic invariant subspace of a skew-Hamiltonian matrix satisfies

cW(X) := lim

ε→0sup

kEkF≤ε

Eskew-Hamiltonian

kΘ(X,ˆ

X)kF

ε≤ kT†

Wk.

It can be shown that actually cW(X) = kT†

Wkholds [Kressner, 2003d].

An extension of this condition number to lower-dimensional isotropic

invariant subspaces and a discussion on the computation of kT†

Wkcan

also be found in [Kressner, 2003d].

Example 9 For the matrix W(0,0) from Example 4, the structured con-

dition number of the Langrangian invariant subspace Xspanned by the

columns of [I2,0]Tis small, since

cW(X) = k(I2⊗diag(1,2) −diag(1,2) ⊗I2)†k2= 1.

This implies that a strongly backward stable method is guaranteed to

compute an excellent approximation of X. On the other hand, the un-

structured condition number c(X)must be considered as infinite due to

the fact that Xis not simple. Hence, a method which is not strongly

backward stable may return arbitrarily bad results.

Similar remarks apply to the eigenvectors of W.

2.3 Algorithms

In Section 2.1 we used a constructive approach to prove the skew-

Hamiltonian Schur decomposition

UTWU =·T˜

0TT¸,(24)

where Uis orthogonal symplectic and Thas real Schur form. The fol-

lowing algorithm summarizes this construction.

Algorithm 10 (Skew-Hamiltonian Schur decomposition)

Input: A skew-Hamiltonian matrix W∈R2n×2n.

Output: An orthogonal symplectic matrix U∈R2n×2n;Wis overwritten

with UTW U having skew-Hamiltonian Schur form (24).

Apply Algorithm 2 to compute an orthogonal symplectic matrix

Uso that W←UTWU has PVL form (6).

Apply the QR algorithm to the (1,1) block W11 of W to compute

an orthogonal matrix Qso that QTW11Qhas real Schur form.

Update W←(Q⊕Q)TW(Q⊕Q), U←U(Q⊕Q).

This algorithm requires around 20n3flops if only the eigenvalues are

desired, and 44n3flops if the skew-Hamiltonian Schur form and the

orthogonal symplectic factor Uare computed, where we used the flop

estimates for the QR algorithm listed in [Golub and Van Loan, 1996,

p.359]. This compares favorably with the QR algorithm applied to the

whole matrix W, which takes 80n3and 200n3flops, respectively. The

finite precision properties of this algorithm are as follows. Similarly as

for the QR algorithm [Wilkinson, 1965] one can show that there is an

orthogonal symplectic matrix Vwhich transforms the computed skew-

Hamiltonian Schur form ˆ

W=hˆ

TTito a skew-Hamiltonian matrix

near to W, i.e., Vˆ

WV T=W+E, where Eis skew-Hamiltonian, kEk2=

O(u)kWk2and udenotes the unit roundoff. Moreover, the computed

factor ˆ

Uis almost orthogonal in the sense that kˆ

UTˆ

U−Ik2=O(u),

and it has the block representation ˆ

U=hˆ

−ˆ

U1i. This implies that

Uis close to an orthogonal symplectic matrix, see e.g. [Kressner, 2004,

Lemma 4.6].

Once a skew-Hamiltonian Schur decomposition has been computed,

the eigenvalues can be easily obtained from the diagonal blocks of T.

Furthermore, if the (k+ 1, k) entry of Tis zero, then the first kcolumns

of Uspan an isotropic invariant subspace of W. Other isotropic invariant

subspaces can be obtained by swapping the diagonal blocks of Tas

described, e.g., in [Bai and Demmel, 1993; Bai et al., 1993].

Symplectic QR decomposition The following algorithm is only in-

directly related to the skew-Hamiltonian eigenvalue problem. It can be

used, for example, to compute orthonormal bases for isotropic subspaces.

Let A∈R2m×nwith m≥n, then there exists an orthogonal symplectic

matrix Q∈R2m×2mso that A=QR and

R=·R11

R21 ¸, R11 ="@

0#, R21 ="

...

0#,(25)

that is, the matrix R11 ∈Rm×nis upper triangular and R21 ∈Rm×nis

strictly upper triangular. A decomposition of this form is called sym-

plectic QR decomposition [Bunse-Gerstner, 1986] and can be computed

by the following algorithm.

Algorithm 11 (Symplectic QR decomposition)

Input: A general matrix A∈R2m×nwith m≥n.

Output: An orthogonal symplectic matrix Q∈R2m×2m;Ais overwritten

with R=QTAhaving the form (25).

Q←I2m.

for j←1,...,n

Set x←Aej.

Apply Algorithm 1 to compute Ej(x).

Update A←Ej(x)TA,Q←QEj(x).

end for

If properly implemented this algorithm requires 8(mn2−n3/3) + O(n2)

flops for computing the matrix R, and additionally 16

3n3+ 16m2n−

16mn2+O(n2) flops for accumulating the orthogonal factor Qin reversed

order.

Other algorithms Similarly as the Hessenberg form of a general ma-

trix can be computed by Gauss transformations [Golub and Van Loan,

1996, Sec. 7.4.7], it is shown in [Stefanovski and Trenˇcevski, 1998] how

non-orthogonal symplectic transformations can be used to compute the

PVL form of a skew-Hamiltonian matrix. A modification of the Arnoldi

method, suitable for computing eigenvalues and isotropic invariant sub-

spaces of large and sparse skew-Hamiltonian matrices, has been proposed

by [Mehrmann and Watkins, 2000].

Balancing a matrix by a simple and accurate similarity transformation

may have a positive impact on the performance of numerical methods

for computing eigenvalues. A structure-preserving balancing procedure

based on symplectic similarity transformations is presented in [Benner,

2000].

The LAPACK [Anderson et al., 1999] subroutines for computing stan-

dard orthogonal decompositions, such as the QR or Hessenberg decom-

position, attain high efficiency by (implicitly) employing W Y representa-

tions of the involved orthogonal transformations [Bischof and Van Loan,

1987; Dongarra et al., 1989; Schreiber and Van Loan, 1989]. A variant

of this representation can be used to derive efficient block algorithms for

computing orthogonal symplectic decompositions, such as the symplectic

QR and URV decompositions [Kressner, 2003a].

3. The Hamiltonian Eigenvalue Problem

One of the most remarkable properties of a Hamiltonian matrix H=

−ATi∈R2n×2nis that its eigenvalues always occur in pairs {λ, −λ},

if λ∈R,λ∈ıR, or in quadruples {λ, −λ, ¯

λ, −¯

λ}, if λ∈C\(R∪ıR).

The preservation of these pairings in finite precision arithmetic is a ma-

jor benefit of using a structure-preserving algorithm for computing the

eigenvalues of H.

Generally, we will only briefly touch the difficulties that arise when

Hhas eigenvalues on the imaginary axis. Although this case is well-

analyzed with respect to structured decompositions, see [Lin and Ho,

1990; Lin et al., 1999; Freiling et al., 2002] and the references given

therein, it is still an open research problem to define appropriate struc-

tured condition numbers and design satisfactory algorithms for this case.

3.1 Structured Decompositions

A major difficulty in developing computational methods for the Hamil-

tonian eigenvalue problem is that there is so far no O(n3) method for

computing a useful structured Hessenberg-like form known. Although a

slight modification of Algorithm 2 can be used to construct an orthogo-

nal symplectic matrix Uso that

UTHU =·W11 W12

W21 WT

11 ¸=



@ @

@

,

i.e., W11 has upper Hessenberg form and W21 is a diagonal matrix, this

form is of limited use. The Hamiltonian QR algorithm, see Section 3.3

below, only preserves this form if the (2,1) block can be written as W21 =

γeneT

nfor some γ∈R. In this case, UTHU is called a Hamiltonian

Hessenberg form. Byers derived in [Byers, 1983] a simple method for

reducing Hto such a form under the assumption that one of the off-

diagonal blocks Gor Qin Hhas tiny rank, i.e., rank 1, 2 or at most

The general case, however, remains elusive. That it might be diffi-

cult to find a simple method is indicated by a result in [Ammar and

Mehrmann, 1991], which shows that the first column xof an orthogonal

symplectic matrix Uthat reduces Hto Hamiltonian Hessenberg form

has to satisfy the set of nonlinear equations

xTJH2i−1x= 0, i = 1,...,n.

This result can even be extended to non-orthogonal symplectic transfor-

mations [Raines and Watkins, 1994].

A Schur-like form for Hamiltonian matrices is given by the following

theorem [Paige and Van Loan, 1981; Lin et al., 1999].

Theorem 12 Let Hbe a Hamiltonian matrix and assume that all eigen-

values of Hthat are on the imaginary axis have even algebraic multiplic-

ity. Then, there is an orthogonal symplectic matrix Uso that UTHU is

in Hamiltonian Schur form, i.e.,

UTHU =·T˜

0−TT¸,(26)

where T∈Rn×nhas real Schur form (9).

If Hhas no eigenvalues on the imaginary axis, then the invariant

subspace Xbelonging to the n(counting multiplicities) eigenvalues in

the open left half plane is called the stable invariant subspace of H. By a

suitable reordering of the Hamiltonian Schur form, see also Section 3.3,

one can see that Xis isotropic. If the columns of Xform an orthonormal

basis for X, then [X, JX] is orthogonal and we have the Hamiltonian

block-Schur decomposition

£X JX ¤TH£X JX ¤=·A11 G11

0−AT

11 ¸.

3.2 Structured Condition Numbers

An extensive perturbation analysis of (block) Hamiltonian Schur forms

for the case that Hhas no purely imaginary eigenvalues has been pre-

sented in [Konstantinov et al., 2001]. The analysis used therein is based

on the technique of splitting operators and Lyapunov majorants. The

approach used in this section is somewhat simpler; it is based on the

perturbation expansions given in Theorem 5.

Structured condition numbers for eigenvalues Let λbe a simple

eigenvalue of a Hamiltonian matrix Hwith right and left eigenvectors x

and y, respectively. The perturbation expansion (13) implies that for a

sufficiently small perturbation E, there is an eigenvalue ˆ

λof W+Eso

that

|ˆ

λ−λ|=|yHEx|

|yHx|+O(kEk2

2)≤kxk2·kyk2

|yHx|kEk2+O(kEk2

2).(27)

If λis real then we may assume that xand yare real and normalized so

that kxk2=kyk2= 1. For the Hamiltonian perturbation E=ε[y, Jx]·

[x, JTy]Hwe have |yHEx|=ε(1 + |yHJx|2) and

kEk2=εk[x, Jy]k2

2=ε(1 + |yHJx|).

The minimum of (1+|yHJx|2)/(1+|yHJx|) is β= 2√2−2. This implies

that for ε→0 both sides in (27) differ at most by a factor 1/β. Hence,

the structured eigenvalue condition number for a simple eigenvalue of a

Hamiltonian matrix,

cH(λ) := lim

ε→0sup

kEk2≤ε

Eis Hamiltonian

|ˆ

λ−λ|

ε,

satisfies the inequalities

(2√2−2)c(λ)≤cH(λ)≤c(λ),

if λ∈R. This inequality still holds for complex λif one allows complex

Hamiltonian perturbations E, i.e., (EJ)H=EJ. A tight lower bound

for the structured condition number of a complex eigenvalue under real

perturbations is an open problem.

Structured backward errors and condition numbers for eigenvalues

of Hamiltonian matrices with additional structures can be found in

[Tisseur, 2003].

Structured condition numbers for invariant subspaces Let the

columns of X∈R2n×kspan a simple, isotropic invariant subspace Xof

H. By the symplectic QR decomposition (25) there is always a matrix

Y∈R2n×kso that U= [X, Y, JX, JY ] is an orthogonal symplectic

matrix. Moreover, we have the block Hamiltonian Schur form

UTHU =





A11 A12 G11 G12

0A22 GT

12 G22

0 0 −AT

11 0

0Q22 −AT

12 −AT







Assuming that the perturbation Eis sufficiently small, the perturbation

expansion (14) implies that there is a matrix ˆ

Xso that ˆ

X= span ˆ

Xis

an invariant subspace of H+Esatisfying

X=X−X⊥T−1

HXT

⊥EX +O(kEk2

F),

and ˆ

XT(ˆ

X−X) = 0. The involved Sylvester operator THis given by

TH: (R1, R2, R3)7→ 



A22 GT

12 G22

0−AT

11 0

Q22 −AT

12 −AT

22 





R3

−



R3

A11.

If the perturbation Eis Hamiltonian, then XT

⊥EX takes the form

S= [AT

21, QT

11, QT

21]T,

where Q11 ∈symm(k) and A21, Q21 are general (n−k)×kmatrices.

Hence, if we let

kT−1

Hk:= sup

S6=0 (kT−1

H(S)kF

kSkF|S∈R(n−k)×k×symm(k)×R(n−k)×k),

then the structured condition number for an isotropic invariant subspace

of a Hamiltonian matrix satisfies

cH(X) = lim

ε→0sup

kEkF≤ε

EHamiltonian

kΘ(X,ˆ

X)kF

ε=kT−1

Hk.

Obviously, this quantity coincides with the unstructured condition num-

ber if Xis one-dimensional, i.e., Xis spanned by a real eigenvector. A

less trivial observation is that the same holds if Xis the stable invariant

subspace, i.e., the n-dimensional subspace belonging to all eigenvalues

in the left half plane. To show this, first note that in this case

kT−1

Hk= sup

S6=0

S∈symm(n)

kT−1

H(S)kF

kSkF

= inf

S6=0

S∈symm(n)

kA11S+SAT

11kF

kSkF

Using a result in [Byers and Nash, 1987], we have

inf

S6=0

S∈symm(n)

kA11S+SAT

11kF

kSkF

= inf

S6=0 kA11S+SAT

11kF

kSkF

which indeed shows that the structured and unstructured condition num-

bers for the maximal stable invariant subspace coincide.

However, there is a severe loss if we do not require Eto be Hamil-

tonian; the subspace ˆ

Xmight not be isotropic. To obtain a nearby

isotropic subspace one can apply the symplectic QR decomposition to

an orthonormal basis ˆ

Xof ˆ

X. This yields the orthonormal basis Zof

an isotropic subspace Z= span Zso that

kZ−XkF≤2kˆ

X−XkF≤2cH(X)kEkF+O(kEk2

F).

Note that for the original subspace ˆ

Xwe have the desirable property

kˆ

⊥Hˆ

XkF=kEkF, where the columns of ˆ

X⊥form an orthonormal

basis for ˆ

X⊥. For the isotropic subspace Z, however, we can only guar-

antee

k(JZ)THZkF≤4cH(X)·kHkF·kEkF+O(kEk2

F),

which signals a severe loss of backward stability. The following numerical

example demonstrates the undesirable appearance of the factor cH(X)

in k(JZ)THZkF.

Example 13 Let

H=





−10−5−1 1 0

1 0 0 1

0 0 10−5−1

0 0 1 0







and consider the stable invariant subspace spanned by the columns of

X= [I2,0]T, which has condition number 105. If we add a random (non-

Hamiltonian) perturbation Ewith kEkF= 10−10 to H, and compute

(using Matlab) an orthonormal basis ˆ

Xfor the invariant subspace ˆ

of H+Ebelonging to the eigenvalues in the open left half plane, we

observe that

kˆ

⊥Hˆ

XkF≈4.0×10−11.

By computing a symplectic QR decomposition of ˆ

Xwe constructed an

orthonormal basis Zsatisfying ZT(JZ) = 0 and observed

k(J˜

Z)TH˜

ZkF≈4.7×10−6.

3.3 Algorithms

An explicit Hamiltonian QR algorithm The Hamiltonian QR al-

gorithm [Byers, 1983] is a strongly backward stable method for com-

puting the Hamiltonian Schur form of a Hamiltonian matrix Hwith no

purely imaginary eigenvalues. Its only obstacle is that there is no im-

plementation of complexity less than O(n4) known, except for the case

when a Hamiltonian Hessenberg form exists [Byers, 1983; Byers, 1986].

One iteration of the Hamiltonian QR algorithm computes the sym-

plectic QR decomposition of the first ncolumns of the symplectic matrix

M= [(H−σ1I)(H−σ2I)][(H+σ1I)(H+σ2I)]−1,(28)

where {σ1, σ2}is a pair of real or complex conjugate shifts. This yields

an orthogonal symplectic matrix Uso that

UTM=



@

.(29)

The next iterate is obtained by updating H←UTHU. Let us partition

Has follows

H=





A11 A12 G11 G12

A21 A22 GT

12 G22

Q11 Q12 −AT

11 −AT

12 Q22 −AT

12 −AT







,(30)

with A11 ∈R2×2and A22 ∈Rn−2×n−2. Under rather mild assumptions

and a fortunate choice of shifts, it can be shown that the submatrices

A21,Q11 and Q12 converge to zero, i.e., Hconverges to a Hamiltonian

block-Schur form [Watkins and Elsner, 1991]. Choosing the shifts s1, s2

as those eigenvalues of the submatrix hA11

Q11

G11

−AT

11 ithat have positive

real part results in quadratic convergence. If this submatrix has two

imaginary eigenvalues, then we suggest to choose the one eigenvalue

with positive real part twice, and if there are four purely imaginary

eigenvalues, then our suggestion is to choose random shifts.

If the norms of the blocks A21,Q11 and Q12 become less than u·

kHkF, then we may safely regard them as zero and apply the iteration

to the submatrix hA22

Q22

G22

−AT

22 i. This will finally yield a Hamiltonian Schur

form of H. Note that the Hamiltonian QR algorithm will generally not

converge if Hhas eigenvalues on the imaginary axis. In our numerical

experiments, however, we often observed convergence to a Hamiltonian

block-Schur form, where the unreduced block hA22

Q22

G22

−AT

22 icontains all

eigenvalues on the imaginary axis.

Remark 14 One can avoid the explicit computation of the potentially

ill-conditioned matrix Min (28) by the following product QR decompo-

sition approach. First, an orthogonal matrix Qris computed so that

(H+σ1I)(H+σ2I)QT

rhas the block triangular structure displayed

in (29). This can be achieved by a minor modification of the standard

RQ decomposition [Benner et al., 1998]. Secondly, the orthogonal sym-

plectic matrix Uis computed from the symplectic QR decomposition of

the first ncolumns of (H−σ1I)(H−σ2I)QT

Reordering a Hamiltonian Schur decomposition If the Hamil-

tonian QR algorithm has successfully computed a Hamiltonian Schur

decomposition,

UTHU =·T˜

0−TT¸(31)

then the first ncolumns of the orthogonal symplectic matrix Uspan an

isotropic subspace belonging to the eigenvalues of T. Many applications

require the stable invariant subspace, for this purpose the Schur decom-

position (31) must be reordered so that Tcontains all eigenvalues with

negative real part.

One way to achieve this is as follows. If there is a block in Twhich

contains a real eigenvalue or a pair of complex conjugate eigenvalues

with positive real part, then this block is swapped to the bottom right

diagonal block T22 of Tusing the algorithms described in [Bai and Dem-

mel, 1993; Bai et al., 1993]. Now, let G22 denote the corresponding block

in ˜

G; it remains to find an orthogonal symplectic matrix U22 so that

22 ·T22 G22

0−TT

22 ¸U22 =·˜

T22 ˜

G22

0−˜

22 ¸(32)

and the eigenvalues of ˜

T22 have negative real part. If Xis the solution

of the Lyapunov equation T22X−XT T

22 =G22, then Xis symmetric

and the columns of [−X, I]Tspan an isotropic subspace. Thus, there is

a symplectic QR decomposition

·−X

I¸=U22 ·R

0¸

By direct computation, it can be shown that U22 is an orthogonal sym-

plectic matrix which produces a reordering of the form (32). Bai and

Demmel, 1993, show that in some pathological cases, the norm of the

(2,1) block in the reordered matrix may be larger than O(u)kHkF. In

this case, which may only happen if the eigenvalues of T22 are close to

the imaginary axis, the swap must be rejected in order to guarantee the

strong backward stability of the algorithm. A different kind of reorder-

ing algorithm, which is based on Hamiltonian QR iterations with perfect

shifts, can be found in [Byers, 1983].

Conclusively, we have a method for computing eigenvalues and se-

lected invariant subspaces of Hamiltonian matrices. This method is

strongly backward stable and reliable, as long as there are no eigen-

values on the imaginary axis. However, as mentioned in the beginning,

in general it requires O(n4) flops, making it unattractive for decently

large problems.

Algorithms based on H2One of the first O(n3) structure-preserving

methods for the Hamiltonian eigenvalue problem was developed in [Van

Loan, 1984b]. It is based on the fact that H2is a skew-Hamiltonian

matrix, because

(H2J)T= (HJ)THT=HJHT=−H(HJ)T=−H2J.

Thus, one can apply Algorithm 10 to H2and take the positive and neg-

ative square roots of the computed eigenvalues, which gives the eigen-

values of H. An implicit version of this algorithm has been implemented

in [Benner et al., 2000]. The main advantage of this approach is that

the eigenvalue symmetries of Hare fully recovered in finite precision

arithmetic. Also, the computational cost is low when compared to the

QR algorithm. The disadvantage of Van Loan’s method is that a loss

of accuracy up to half the number of significant digits of the computed

eigenvalues of His possible. An error analysis in [Van Loan, 1984b]

shows that for an eigenvalue λof Hthe computed ˆ

λsatisfies

|ˆ

λ−λ|/c(λ)·min{ukHk2

2/|λ|,√ukHk2}.

This indicates that particularly eigenvalues with |λ| ¿ kHk2are affected

by the √u-effect. Note that a similar effect occurs when one attempts to

compute the singular values of a general matrix Afrom the eigenvalues

of ATA, see e.g. [Stewart, 2001, Sec. 3.3.2].

An algorithm that is based on the same idea but achieves numeri-

cal backward stability by completely avoiding the squaring of Hwas

developed in [Benner et al., 1998]. In the following, we show how this

algorithm can be directly derived from Algorithm 10. In lieu of H2we

make use of the skew-Hamiltonian matrix

W=





0A0G

−A0−G0

0Q0−AT

−Q0AT0





∈R4n×4n,(33)

for given H=hA

−ATi. As Wis permutationally equivalent to h0

−H

0i,

we see that ±λis an eigenvalue of Hif and only if ±√−λ2is an eigen-

value of W. Note that the matrix Whas a lot of extra structure besides

being skew-Hamiltonian, which is not exploited if we apply Algorithm 10

directly to W.

Instead, we consider the shuffled matrix ˜

W= (P⊕P)TW(P⊕P),

where

P=£e1e3··· e2n−1e2e4··· e2n¤.

This matrix has the form

W=·˜

WA˜

WQ˜

A¸,

where each of the matrices ˜

WA,˜

WGand ˜

WQis a block matrix composed

of two-by-two blocks having the form

WX=µ· 0xij

−xij 0¸¶n

i,j=1

If an orthogonal symplectic matrix ˜

Qhas the form

Q= (P⊕P)T





V10V20

0U10U2

−V20V10

0−U20U1







(P⊕P),(34)

then ˜

QT˜

W˜

Qis skew-Hamiltonian and has the same zero pattern as ˜

Lemma 15 The orthogonal symplectic factor of the PVL decomposition

computed by Algorithm 2 has the form (34).

Proof. Assume that after (j−1) loops of Algorithm 2 the matrix ˜

has been overwritten by a matrix with the same zero pattern as ˜

W. Let

˜xdenote the jth column of ˜

W. If jis odd then ˜xcan be written as

˜x=x⊗e2and if jis even then ˜x=x⊗e1, where xis a vector of

length 2nand e1, e2are unit vectors of length two. This implies that

Algorithm 1 produces an elementary orthogonal matrix Ej(x) having

the same zero pattern as the matrix ˜

Qin (34), see [Kressner, 2003c].

This shows that the jth loop of Algorithm 2 preserves the zero pattern

of ˜

W. The proof is concluded by using the fact that the set of matrices

having the form (34) is closed under multiplication.

This also shows that the PVL form returned by Algorithm 2 must

take the form

QT˜

W˜

Q= (P⊕P)T





0R22 0R12

−R11 0−RT

12 0

000−RT

0 0 RT

22 0







(P⊕P),(35)

where R11 is an upper triangular matrix and R22 is an upper Hessenberg

matrix. Rewriting (35) in terms of the block entries of ˜

Wand ˜

Qyields

UTHV =·R11 R12

0−RT

22 ¸=



@

(36)

with the orthogonal symplectic matrices U=hU1

−U2

U1iand V=hV1

−V2

V1i,

formed from the entries of ˜

Qin (34). This is a so called symplectic URV

decomposition [Benner et al., 1998].

As Algorithm 2 exclusively operates on the nonzero entries of ˜

Wit

should be possible to reformulate it purely in terms of these entries. This

amounts to the following algorithm [Benner et al., 1998, Alg. 4.4].

Algorithm 16 (Symplectic URV decomposition)

Input: A matrix H∈R2n×2n.

Output: Orthogonal symplectic matrices U, V ∈R2n×2n;His overwritten

with UTHV having the form (36).

U←I2n, V ←I2n.

for j←1,2,...,n

Set x←Hej.

Apply Algorithm 1 to compute Ej(x).

Update H←Ej(x)TH, U ←UEj(x).

if j < n then

Set y←HTen+j.

Apply Algorithm 1 to compute En+j+1(y).

Update H←HEn+j+1(y), V ←V En+j+1(y).

end if

end for

If properly implemented, Algorithm 16 requires 80

3n3+O(n2) floating

point operations (flops) to reduce Hand additionally 16

3n3+O(n2) flops

to compute each of the orthogonal symplectic factors Uand V. Note

that Algorithm 16 does not require Hto be a Hamiltonian matrix, but

even if His Hamiltonian, this structure will be destroyed.

E7E7

E6E6

Figure 3. Illustration of two loops of Algorithm 16 for n= 4.

In the second step of Algorithm 10, the QR algorithm is applied to

the upper left 2n×2nblock of the PVL form (35). In [Kressner, 2003c]

it is shown that this is equivalent to applying the periodic QR algorithm

[Bojanczyk et al., 1992; Hench and Laub, 1994; Van Loan, 1975] to the

matrix product −R22 ·R11, which constructs orthogonal matrices Q1

and Q2so that QT

1R22Q2is reduced to real Schur form while QT

2R11Q1

stays upper triangular. The periodic QR algorithm is a backward stable

method for computing the eigenvalues of R22 ·R11. The positive and

negative square roots of these eigenvalues are the eigenvalues of H.

The procedure, as described above, is a numerically backward stable

method for computing the eigenvalues of a Hamiltonian matrix H. It

preserves the eigenvalue symmetries of Hin finite precision arithmetic

and its complexity is O(n3). As the periodic QR algorithm inherits the

reliability of the standard QR algorithm, this method can be regarded as

highly reliable. Its only drawback is that it does not take full advantage

of the structure of H. It is not clear whether the method is strongly

backward stable or not.

Computation of invariant subspaces based on H2Having com-

puted an invariant subspace for the skew-Hamiltonian matrix H2it

is possible to extract invariant subspaces for Hfrom it [Xu and Lu,

1995; Hwang et al., 2003]. However, we have already observed that the

explicit computation of H2can lead to numerical instabilities and should

be avoided. The above idea of embedding Hin a skew-Hamiltonian ma-

trix Wof double dimension can be extended for computing invariant

subspaces, see [Benner et al., 1997]. However, it should be noted that

this approach might encounter numerical difficulties if Hhas eigenvalues

on or close to the imaginary axis.

Refinement of stable invariant subspaces With all the difficulties

in deriving a strongly backward stable method it might be preferable to

use some kind of iterative refinement algorithms to improve the quanti-

ties computed by a less stable method. This idea is used, for example,

in the multishift algorithm [Ammar et al., 1993] and hybrid methods for

solving algebraic Riccati equations [Benner and Faßbender, 2001].

In the following we describe a method for improving an isotropic sub-

space ˆ

Xthat approximates the stable invariant subspace Xof a Hamil-

tonian matrix H. Let the columns of ˆ

Xspan an orthonormal basis for

Xand consider

£ˆ

X J ˆ

X¤TH£ˆ

X J ˆ

X¤=·˜

A˜

Q−˜

AT¸.

If ˆ

Xhas been computed by a strongly backward stable method then

k˜

Qkis of order u·kHkand it is not possible to refine ˆ

Xmuch further.

However, as we have seen before, if a less stable method has been used

then k˜

Qkmight be much larger. In this case we can apply the following

algorithm to improve the accuracy of ˆ

Algorithm 17

Input: A Hamiltonian matrix H∈R2n×2n, a matrix ˆ

X∈R2n×nso

that [ˆ

X, J ˆ

X]is orthogonal, and a tolerance tol >0.

Output: The matrix ˆ

Xis updated until k(Jˆ

X)THˆ

XkF≤tol · kHkF.

while k(Jˆ

X)THˆ

XkF>tol ·kHkF

Set ˜

A←ˆ

XTHˆ

Xand ˜

Q←(Jˆ

X)THˆ

Solve the Lyapunov equation R˜

A+˜

ATR=−˜

Compute Y∈R2n×nso that [Y, JY ] is orthogonal and

span Y= span ³h I

−Ri´using a symplectic QR decomposition.

Update [ ˆ

X, J ˆ

X]←[ˆ

X, J ˆ

X]·[Y, JY ].

end while

As this algorithm is a special instance of a Newton method for refining

invariant subspaces [Stewart, 1973; Chatelin, 1984; Demmel, 1987; Ben-

ner, 1997] or a block Jacobi-like algorithm [H¨uper and Van Dooren,

2003] it converges locally quadratic. On the other hand, Algorithm 17

can be seen as a particular implementation of a Newton method for solv-

ing algebraic Riccati equation [Kleinman, 1968; Lancaster and Rodman,

1995; Mehrmann, 1991]. By a more general result in [Guo and Lancaster,

1998], this implies under some mild conditions global convergence if H

has no eigenvalues on the imaginary axis and if the iteration is initialized

with a matrix ˆ

Xso that all eigenvalues of ˜

A=ˆ

XTHˆ

Xare in the open

left half plane C−.

In finite precision arithmetic, the minimal attainable tolerance is tol ≈

n2·uunder the assumption that a forward stable method such as the

Bartels-Stewart method [Bartels and Stewart, 1972] is used to solve the

Lyapunov equations R˜

A+˜

ATR=−˜

Q[Higham, 1996; Tisseur, 2001].

Other Algorithms As mentioned in the introduction there is a vast

number of algorithms for the Hamiltonian eigenvalue problem avail-

able. Other algorithms based on orthogonal transformations are the

Hamiltonian Jacobi algorithm [Byers, 1990; Bunse-Gerstner and Faßben-

der, 1997], its variants for Hamiltonian matrices that have additional

structure [Faßbender et al., 2001] and the multishift algorithm [Ammar

et al., 1993]. Algorithms based on symplectic but non-orthogonal trans-

formations include the SR algorithm [Bunse-Gerstner and Mehrmann,

1986; Bunse-Gerstner, 1986; Mehrmann, 1991] as well as related methods

[Bunse-Gerstner et al., 1989; Raines and Watkins, 1994]. A completely

different class of algorithms is based on the matrix sign function, see,

e.g., [Benner, 1999; Mehrmann, 1991; Sima, 1996] and the references

therein. Other Newton-like methods directed towards the computation

of invariant subspaces for Hamiltonian matrices can be found in [Absil

and Van Dooren, 2002; Guo and Lancaster, 1998].

A structure-preserving Arnoldi method based on the H2approach was

developed in [Mehrmann and Watkins, 2000; Kressner, 2004]. There are

also a number of symplectic Lanczos methods available, see [Benner and

Faßbender, 1997; Ferng et al., 1997; Watkins, 2002].

The remarks on balancing and block algorithms at the end of Sec-

tion 2.3, carry over to Hamiltonian matrices. We only note that in

[Benner and Kressner, 2003], a balancing algorithm is described which

is particularly suited for large and sparse Hamiltonian matrices.

4. Applications

Most applications of skew-Hamiltonian and Hamiltonian eigenvalue

problems are in the area of systems and control theory. In the following,

we consider a linear continuous-time system with constant coefficients,

which can be described by a set of matrix differential and algebraic

equations

˙x(t) = Ax(t) + Bu(t), x(0) = x0,

y(t) = Cx(t) + Du(t),(37)

where x(t)∈Rnis the vector of states,u(t)∈Rmthe vector of inputs

(or controls) and y(t)∈Rrthe vector of outputs at time t∈[0,∞). The

system is described by the state matrix A∈Rn×n, the input (control)

matrix B∈Rn×m, the output matrix C∈Rr×nand the feedthrough

matrix D∈Rr×m. It is much beyond the scope of this work to give

an introduction to such systems; for this purpose the reading might

be complemented with any modern, state-space oriented monograph in

this area, see e.g. [Green and Limebeer, 1995; Petkov et al., 1991; Van

Dooren, 2003; Zhou et al., 1996].

4.1 Stability Radius Computation

The system (37) is called (asymptotically) stable if all eigenvalues

λ(A) of the state matrix Alie in C−. It is often important to know how

near the system is to an unstable one, i.e., what is the smallest norm

perturbation E∈Cn×nso that λ(A+E)6⊂ C−. This corresponds to

the computation of the stability radius of A, which is defined as

γ(A) := min{kEk2:λ(A+E)∩ıR6=∅}.

A bisection method for measuring γ(A) can be based on the following

observation [Byers, 1988]: if α≥0, then the Hamiltonian matrix

H(α) = ·A−αIn

αIn−AT¸

has an eigenvalue on the imaginary axis if and only if α≥γ(A). This

suggests a simple bisection algorithm. Start with a lower bound β≥0

and an upper bound δ > γ(A) (an easy-to-compute upper bound is

kA+ATkF/2 [Van Loan, 1984a]). Then in each step, set α:= (β+δ)/2

and compute λ(H(α)). If there is an eigenvalue on the imaginary axis,

choose δ=α, otherwise, set β=α.

The correct decision whether H(α) has eigenvalues on the imaginary

axis is crucial for the success of the bisection method. [Byers, 1988] shows

that if the eigenvalues of H(α) are computed by a strongly backward

stable method, then the computed γ(A) will be within an O(u)·kAk2-

distance of the exact stability radius.

4.2 H∞Norm Computation

A similar problem is the computation of the H∞norm of a stable

system. Consider the transfer function G(s) of a stable system of the

form (37),

G(s) = C(sI −A)−1B+D,

then

kGkH∞= esssup{kG(ıω)k2:ω∈R}.

is the H∞norm of G, see e.g. [Green and Limebeer, 1995; Zhou et al.,

1996].

Let σmax(D) denote the largest singular value of Dand let α∈R

be such that α > σmax(D). Then consider the parameter-dependent

Hamiltonian matrix

H(α) = ·H11(α)H12(α)

H21(α)−H11(α)T¸,

where for R(α) = α2I−DTD,

H11(α) = A+BR(α)−1DTC,

H12(α) = 1

α2BR(α)−1BT,

H21(α) = −CT(I+DR(α)−1DT)C.

The following result can be used to approximate kGkH∞, see e.g. [Zhou

et al., 1996]:

kGkH∞< α ⇔σmax(D)< α and λ(H(α)) ∩ıR=∅.

Using this fact, a bisection algorithm analogous to the stability radius

computation can be formulated, starting with lower bound β=σmax(D)

and upper bound δ > kGkH∞, see [Boyd et al., 1989] for details. Again,

the bisection algorithm benefits if the decisions are based on eigenvalues

computed by a method preserving the eigenvalue symmetries of H(α).

Faster convergent versions of this algorithm, which may also involve

the eigenvectors of H(α), can be found in [Boyd and Balakrishnan, 1990;

Bruinsma and Steinbuch, 1990; Genin et al., 1998].

4.3 Algebraic Riccati Equations

Given a Hamiltonian matrix Has in (2), there is always a correspond-

ing algebraic Riccati equation (ARE)

0 = Q+ATX+XA −XGX. (38)

AREs have played a fundamental role in systems and control theory

since the early 1960’s as they are the major tool to compute feedback

controllers using LQR/LQG (H2) or H∞approaches. The correspon-

dence between feedback controllers and AREs can be found in liter-

ally any modern textbook on control, see, e.g., [Anderson and Moore,

1990; Green and Limebeer, 1995; Zhou et al., 1996] and many others. In

these applications, usually a particular solution of (38) is required which

is stabilizing in the sense that λ(A−GX) is contained in the open left

half plane. This solution is unique if it exists and is related to the Hami-

iltonian eigenproblem as follows. Suppose Xis a symmetric solution of

(38), then it is easy to see that

H·In0

−X In¸=·In0

−X In¸· A−GX G

0−(A−GX)T¸.

Hence, the columns of [In, X]Tspan an H-invariant subspace corre-

sponding to λ(H)∩λ(A−GX). This implies that we can solve AREs by

computing H-invariant subspaces. In particular, if we want the stabiliz-

ing solution, we need the maximal stable H-invariant subspace. Suppose

that a basis of this subspace is given by the columns of [XT

1, XT

2]Twith

X1, X2∈Rn×nthen, under mild assumptions, X1is invertible and X=

−X2X−1

1is the stabilizing solution of (38). Therefore, any algorithm to

compute invariant subspaces of Hamiltonian matrices may be used to

solve AREs. For discussions of this topic see [Benner, 1999; Mehrmann,

1991; Sima, 1996]. It should be noted, though, that often the ARE is

a detour. In feedback control, the solution of the ARE can usually be

avoided by working only with the H-invariant subspaces, see [Benner

et al., 2004; Mehrmann, 1991].

A correspondence between the skew-Hamiltonian eigenproblem (1)

and the anti-symmetric ARE

0 = Q−ATX+XA −XGX, Q =−QT, G =−GT,

is discussed in [Stefanovski and Trenˇcevski, 1998].

4.4 Quadratic Eigenvalue Problems

The quadratic eigenvalue problem (QEP) is to find scalars λand

nonzero vectors xsatisfying

(λ2M+λG +K)x= 0,(39)

where M, G, K ∈Rn×n. It arises, for example, from linear systems

that are governed by second order differential equations, see [Tisseur

and Meerbergen, 2001]. Gyroscopic systems yield QEPs with symmetric

positive definite M, skew-symmetric Gand symmetric K. In this case,

the eigenvalues of (39) have the same symmetries as in the Hamiltonian

eigenvalue problem, i.e., if λis an eigenvalue then −λ, ¯

λand −¯

λare also

eigenvalues.

By [Mehrmann and Watkins, 2000], a linearization of (39) reflect-

ing this property is the skew-Hamiltonian/Hamiltonian matrix pencil

λWLWR−H, where

WL=·I1

0M¸, WR=·M1

0I¸, H =·0−K

M0¸.

From this, it is easy to see that W−1

LHW−1

Ris Hamiltonian and has the

same eigenvalues as (39). Hence, a structure-preserving algorithm ap-

plied to W−1

LHW−1

Rwill preserve the eigenvalue pairings of (39). This

is particularly important for testing the stability of the underlying gy-

roscopic system, which amounts to checking whether all eigenvalues

of (39) are on the imaginary axis, see e.g. [Tisseur and Meerbergen,

2001, Sec.5.3].

However, it should be noted that such an approach is only advis-

able as long as the matrix Mis sufficiently well conditioned. Other-

wise, structure-preserving algorithms that work directly on the pencil

λWLWR−Hshould be preferred [Benner et al., 1998].

Linearizations that lead to skew-Hamiltonian eigenvalue problems are

described in [Mehrmann and Watkins, 2000], and have been used for

computing corner singularities in anisotropic elastic structures [Apel

et al., 2002].

4.5 Other Applications

Other applications for Hamiltonian eigenvalue problems include pas-

sivity preserving model reduction [Antoulas and Sorensen, 2001; Sorensen,

2002], the computation of pseudospectra [Burke et al., 2003b] and the

distance to uncontrollability [Gu, 2000; Burke et al., 2003a].

5. Concluding Remarks

We have presented structured decompositions, condition numbers, al-

gorithms and applications for skew-Hamiltonian and Hamiltonian eigen-

value problems. It is our hope that the reader is now convinced that the

exploitation of such structures is an interesting area of research, not only

from a theoretical point of view but also with respect to applications.

Many problems remain open. In particular, Hamiltonian matrices with

eigenvalues on the imaginary axis require further investigation.

Most of the presented material is already available in the cited liter-

ature. In this survey, the novel pieces of the (skew-)Hamiltonian puzzle

are:

explicit formulas/bounds for the structured eigenvalue condition

numbers;

a relation between the structured and unstructured condition num-

bers for stable invariant subspaces of Hamiltonian matrices;

a new reordering algorithm for the Hamiltonian Schur form based

on symplectic QR decompositions; and

the derivation of the symplectic URV decomposition from the PVL

decomposition;

a structure-preserving iterative refinement algorithm for stable in-

variant subspaces of Hamiltonian matrices.

Acknowledgments

The authors gratefully thank Ralph Byers and Vasile Sima for helpful

discussions.

References

Absil, P.-A. and Van Dooren, P. (2002). Two-sided Grassmann Rayleigh quotient

iteration.

Ammar, G., Benner, P., and Mehrmann, V. (1993). A multishift algorithm for the

numerical solution of algebraic Riccati equations. Electr. Trans. Num. Anal., 1:33–

48.

Ammar, G. and Mehrmann, V. (1991). On Hamiltonian and symplectic Hessenberg

forms. Linear Algebra Appl., 149:55–72.

Anderson, B. and Moore, J. (1990). Optimal Control – Linear Quadratic Methods.

Prentice-Hall, Englewood Cliffs, NJ.

Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J. J., Du Croz,

J., Greenbaum, A., Hammarling, S., McKenney, A., and Sorensen, D. (1999). LA-

PACK Users’ Guide. Society for Industrial and Applied Mathematics, Philadel-

phia, PA, third edition.

Antoulas, A. C. and Sorensen, D. C. (2001). Approximation of large-scale dynamical

systems: an overview. Int. J. Appl. Math. Comput. Sci., 11(5):1093–1121.

Apel, T., Mehrmann, V., and Watkins, D. S. (2002). Structured eigenvalue methods

for the computation of corner singularities in 3D anisotropic elastic structures.

Comput. Methods Appl. Mech. Engrg, 191:4459–4473.

Bai, Z., Demmel, J., and McKenney, A. (1993). On computing condition numbers for

the nonsymmetric eigenproblem. ACM Trans. Math. Software, 19(2):202–223.

Bai, Z. and Demmel, J. W. (1993). On swapping diagonal blocks in real Schur form.

Linear Algebra Appl., 186:73–95.

Bartels, R. H. and Stewart, G. W. (1972). Algorithm 432: The solution of the matrix

equation AX −BX =C.Communications of the ACM, 8:820–826.

Benner, P. (1997). Contributions to the Numerical Solution of Algebraic Riccati Equa-

tions and Related Eigenvalue Problems. Logos–Verlag, Berlin, Germany.

Benner, P. (1999). Computational methods for linear-quadratic optimization. Supple-

mento ai Rendiconti del Circolo Matematico di Palermo, Serie II, No. 58:21–56.

Benner, P. (2000). Symplectic balancing of Hamiltonian matrices. SIAM J. Sci. Com-

put., 22(5):1885–1904.

Benner, P., Byers, R., and Barth, E. (2000). Algorithm 800: Fortran 77 subroutines for

computing the eigenvalues of Hamiltonian matrices I: The square-reduced method.

ACM Trans. Math. Software, 26:49–77.

Benner, P., Byers, R., Mehrmann, V., and Xu, H. (2004). Robust numerical methods

for robust control. Technical Report 06-2004, Institut f¨ur Mathematik, TU Berlin.

Benner, P. and Faßbender, H. (1997). An implicitly restarted symplectic Lanczos

method for the Hamiltonian eigenvalue problem. Linear Algebra Appl., 263:75–

111.

Benner, P. and Faßbender, H. (2001). A hybrid method for the numerical solution

of discrete-time algebraic Riccati equations. Contemporary Mathematics, 280:255–

269.

Benner, P. and Kressner, D. (2003). Balancing sparse Hamiltonian eigenproblems. To

appear in Linear Algebra Appl.

Benner, P. and Kressner, D. (2004). Fortran 77 subroutines for computing the eigen-

values of Hamiltonian matrices II. In preparation. See also http://www.math.

tu-berlin.de/~kressner/hapack/.

Benner, P., Mehrmann, V., and Xu, H. (1997). A new method for computing the

stable invariant subspace of a real Hamiltonian matrix. J. Comput. Appl. Math.,

86:17–43.

Benner, P., Mehrmann, V., and Xu, H. (1998). A numerically stable, structure pre-

serving method for computing the eigenvalues of real Hamiltonian or symplectic

pencils. Numerische Mathematik, 78(3):329–358.

Bischof, C. and Van Loan, C. F. (1987). The W Y representation for products of

Householder matrices. SIAM J. Sci. Statist. Comput., 8(1):S2–S13. Parallel pro-

cessing for scientific computing (Norfolk, Va., 1985).

Bojanczyk, A., Golub, G. H., and Dooren, P. V. (1992). The periodic Schur decompo-

sition; algorithm and applications. In Proc. SPIE Conference, volume 1770, pages

31–42.

Boyd, S. and Balakrishnan, V. (1990). A regularity result for the singular values

of a transfer matrix and a quadratically convergent algorithm for computing its

L∞-norm. Systems Control Lett., 15(1):1–7.

Boyd, S., Balakrishnan, V., and Kabamba, P. (1989). A bisection method for com-

puting the H∞norm of a transfer matrix and related problems. Math. Control,

Signals, Sys., 2:207–219.

Bruinsma, N. A. and Steinbuch, M. (1990). A fast algorithm to compute the H∞-norm

of a transfer function matrix. Sys. Control Lett., 14(4):287–293.

Bunch, J. R. (1987). The weak and strong stability of algorithms in numerical linear

algebra. Linear Algebra Appl., 88/89:49–66.

Bunse-Gerstner, A. (1986). Matrix factorizations for symplectic QR-like methods.

Linear Algebra Appl., 83:49–77.

Bunse-Gerstner, A. and Faßbender, H. (1997). A Jacobi-like method for solving al-

gebraic Riccati equations on parallel computers. IEEE Trans. Automat. Control,

42(8):1071–1084.

Bunse-Gerstner, A. and Mehrmann, V. (1986). A symplectic QR like algorithm for

the solution of the real algebraic Riccati equation. IEEE Trans. Automat. Control,

31(12):1104–1113.

Bunse-Gerstner, A., Mehrmann, V., and Watkins, D. S. (1989). An SR algorithm

for Hamiltonian matrices based on Gaussian elimination. In XII Symposium on

Operations Research (Passau, 1987), volume 58 of Methods Oper. Res., pages 339–

357. Athen¨aum/Hain/Hanstein, K¨onigstein.

Burke, J. V., Lewis, A. S., and Overton, M. L. (2003a). Pseudospectral components

and the distance to uncontrollability. Submitted to SIAM J. Matrix Anal. Appl.

REFERENCES 43

Burke, J. V., Lewis, A. S., and Overton, M. L. (2003b). Robust stability and a criss-

cross algorithm for pseudospectra. IMA J. Numer. Anal., 23(3):359–375.

Byers, R. (1983). Hamiltonian and Symplectic Algorithms for the Algebraic Riccati

Equation. PhD thesis, Cornell University, Dept. Comp. Sci., Ithaca, NY.

Byers, R. (1986). A Hamiltonian QR algorithm. SIAM J. Sci. Statist. Comput.,

7(1):212–229.

Byers, R. (1988). A bisection method for measuring the distance of a stable to unstable

matrices. SIAM J. Sci. Statist. Comput., 9:875–881.

Byers, R. (1990). A Hamiltonian-Jacobi algorithm. IEEE Trans. Automat. Control,

35:566–570.

Byers, R. and Nash, S. (1987). On the singular “vectors” of the Lyapunov operator.

SIAM J. Algebraic Discrete Methods, 8(1):59–66.

Chatelin, F. (1984). Simultaneous Newton’s iteration for the eigenproblem. In Defect

correction methods (Oberwolfach, 1983), volume 5 of Comput. Suppl., pages 67–74.

Springer, Vienna.

Demmel, J. W. (1987). Three methods for refining estimates of invariant subspaces.

Computing, 38:43–57.

Dongarra, J. J., Sorensen, D. C., and Hammarling, S. J. (1989). Block reduction of

matrices to condensed forms for eigenvalue computations. J. Comput. Appl. Math.,

27(1-2):215–227. Reprinted in Parallel algorithms for numerical linear algebra, 215–

227, North-Holland, Amsterdam, 1990.

Faßbender, H., Mackey, D. S., and Mackey, N. (2001). Hamilton and Jacobi come full

circle: Jacobi algorithms for structured Hamiltonian eigenproblems. Linear Algebra

Appl., 332/334:37–80.

Faßbender, H., Mackey, D. S., Mackey, N., and Xu, H. (1999). Hamiltonian square

roots of skew-Hamiltonian matrices. Linear Algebra Appl., 287(1-3):125–159.

Ferng, W., Lin, W.-W., and Wang, C.-S. (1997). The shift-inverted J-Lanczos al-

gorithm for the numerical solutions of large sparse algebraic Riccati equations.

Comput. Math. Appl., 33(10):23–40.

Freiling, G., Mehrmann, V., and Xu, H. (2002). Existence, uniqueness, and parametriza-

tion of Lagrangian invariant subspaces. SIAM J. Matrix Anal. Appl., 23(4):1045–

1069.

Gantmacher, F. (1960). The Theory of Matrices. Chelsea, New York.

Genin, Y., Van Dooren, P., and Vermaut, V. (1998). Convergence of the calculation

of H∞norms and related questions. In Beghi, A., Finesso, L., and Picci, G., edi-

tors, Proceedings of the Conference on the Mathematical Theory of Networks and

Systems, MTNS ’98, pages 429–432.

Golub, G. H. and Van Loan, C. F. (1996). Matrix Computations. Johns Hopkins

University Press, Baltimore, MD, third edition.

Green, M. and Limebeer, D. (1995). Linear Robust Control. Prentice-Hall, Englewood

Cliffs, NJ.

Gu, M. (2000). New methods for estimating the distance to uncontrollability. SIAM

J. Matrix Anal. Appl., 21(3):989–1003.

Guo, C. and Lancaster, P. (1998). Analysis and modification of Newton’s method for

algebraic Riccati equations. Math. Comp., 67:1089–1105.

Hench, J. J. and Laub, A. J. (1994). Numerical solution of the discrete-time periodic

Riccati equation. IEEE Trans. Automat. Control, 39(6):1197–1210.

Higham, N. J. (1996). Accuracy and stability of numerical algorithms. Society for

Industrial and Applied Mathematics (SIAM), Philadelphia, PA.

H¨uper, K. and Van Dooren, P. (2003). New algorithms for the iterative refinement of

estimates of invariant subspaces. Journal Future Generation Computer Systems,

19:1231–1242.

Hwang, T.-M., Lin, W.-W., and Mehrmann, V. (2003). Numerical solution of quadratic

eigenvalue problems with structure-preserving methods. SIAM J. Sci. Comput.,

24(4):1283–1302.

Kleinman, D. (1968). On an iterative technique for Riccati equation computations.

IEEE Trans. Automat. Control, AC-13:114–115.

Konstantinov, M., Mehrmann, V., and Petkov, P. (2001). Perturbation analysis of

Hamiltonian Schur and block-Schur forms. SIAM J. Matrix Anal. Appl., 23(2):387–

424.

Kressner, D. (2003a). Block algorithms for orthogonal symplectic factorizations. BIT,

43(4):775–790.

Kressner, D. (2003b). A Matlab toolbox for solving skew-Hamiltonian and Hamilto-

nian eigenvalue problems. Online available from http:/www.math.tu-berlin.de/

~kressner/hapack/matlab/.

Kressner, D. (2003c). The periodic QR algorithm is a disguised QR algorithm. To

appear in Linear Algebra Appl.

Kressner, D. (2003d). Perturbation bounds for isotropic invariant subspaces of skew-

Hamiltonian matrices. To appear in SIAM J. Matrix Anal. Appl..

Kressner, D. (2004). Numerical Methods and Software for General and Structured

Eigenvalue Problems. PhD thesis, TU Berlin, Institut f¨ur Mathematik, Berlin,

Germany.

Lancaster, P. and Rodman, L. (1995). The Algebraic Riccati Equation. Oxford Uni-

versity Press, Oxford.

Lin, W.-W. and Ho, T.-C. (1990). On Schur type decompositions for Hamiltonian and

symplectic pencils. Technical report, Institute of Applied Mathematics, National

Tsing Hua University, Taiwan.

Lin, W.-W., Mehrmann, V., and Xu, H. (1999). Canonical forms for Hamiltonian and

symplectic matrices and pencils. Linear Algebra Appl., 302/303:469–533.

Mehrmann, V. (1991). The Autonomous Linear Quadratic Control Problem, Theory

and Numerical Solution. Number 163 in Lecture Notes in Control and Information

Sciences. Springer-Verlag, Heidelberg.

Mehrmann, V. and Watkins, D. S. (2000). Structure-preserving methods for comput-

ing eigenpairs of large sparse skew-Hamiltonian/Hamiltonian pencils. SIAM J. Sci.

Comput., 22(6):1905–1925.

Ortega, J. M. and Rheinboldt, W. C. (1970). Iterative solution of nonlinear equations

in several variables. Academic Press, New York.

Paige, C. and Van Loan, C. F. (1981). A Schur decomposition for Hamiltonian ma-

trices. Linear Algebra Appl., 41:11–32.

Petkov, P. H., Christov, N. D., and Konstantinov, M. M. (1991). Computational Meth-

ods for Linear Control Systems. Prentice-Hall, Hertfordshire, UK.

Raines, A. C. and Watkins, D. S. (1994). A class of Hamiltonian–symplectic methods

for solving the algebraic Riccati equation. Linear Algebra Appl., 205/206:1045–

1060.

REFERENCES 45

Schreiber, R. and Van Loan, C. F. (1989). A storage-efficient W Y representation for

products of Householder transformations. SIAM J. Sci. Statist. Comput., 10(1):53–

57.

Sima, V. (1996). Algorithms for Linear-Quadratic Optimization, volume 200 of Pure

and Applied Mathematics. Marcel Dekker, Inc., New York, NY.

Sorensen, D. C. (2002). Passivity preserving model reduction via interpolation of

spectral zeros. Technical report TR02-15, ECE-CAAM Depts, Rice University.

Stefanovski, J. and Trenˇcevski, K. (1998). Antisymmetric Riccati matrix equation. In

1st Congress of the Mathematicians and Computer Scientists of Macedonia (Ohrid,

1996), pages 83–92. Sojuz. Mat. Inform. Maked., Skopje.

Stewart, G. W. (1971). Error bounds for approximate invariant subspaces of closed

linear operators. SIAM J. Numer. Anal., 8:796–808.

Stewart, G. W. (1973). Error and perturbation bounds for subspaces associated with

certain eigenvalue problems. SIAM Rev., 15:727–764.

Stewart, G. W. (2001). Matrix algorithms. Vol. II. Society for Industrial and Applied

Mathematics (SIAM), Philadelphia, PA. Eigensystems.

Stewart, G. W. and Sun, J.-G. (1990). Matrix Perturbation Theory. Academic Press,

New York.

Sun, J.-G. (1998). Stability and accuracy: Perturbation analysis of algebraic eigen-

problems. Technical report UMINF 98-07, Department of Computing Science, Uni-

versity of Ume˚a, Ume˚a, Sweden.

Tisseur, F. (2001). Newton’s method in floating point arithmetic and iterative refine-

ment of generalized eigenvalue problems. SIAM J. Matrix Anal. Appl., 22(4):1038–

1057.

Tisseur, F. (2003). A chart of backward errors for singly and doubly structured eigen-

value problems. SIAM J. Matrix Anal. Appl., 24(3):877–897.

Tisseur, F. and Meerbergen, K. (2001). The quadratic eigenvalue problem. SIAM

Rev., 43(2):235–286.

Van Dooren, P. (2003). Numerical Linear Algebra for Signal Systems and Control.

Draft notes prepared for the Graduate School in Systems and Control.

Van Loan, C. F. (1975). A general matrix eigenvalue algorithm. SIAM J. Numer.

Anal., 12(6):819–834.

Van Loan, C. F. (1984a). How near is a matrix to an unstable matrix? Lin. Alg. and

its Role in Systems Theory, 47:465–479.

Van Loan, C. F. (1984b). A symplectic method for approximating all the eigenvalues

of a Hamiltonian matrix. Linear Algebra Appl., 61:233–251.

Watkins, D. S. (2002). On Hamiltonian and symplectic Lanczos processes. To appear

in Linear Algebra Appl.

Watkins, D. S. and Elsner, L. (1991). Convergence of algorithms of decomposition

type for the eigenvalue problem. Linear Algebra Appl., 143:19–47.

Wilkinson, J. H. (1965). The algebraic eigenvalue problem. Clarendon Press, Oxford.

Xu, H. and Lu, L. Z. (1995). Properties of a quadratic matrix equation and the

solution of the continuous-time algebraic Riccati equation. Linear Algebra Appl.,

222:127–145.

Zhou, K., Doyle, J. C., and Glover, K. (1996). Robust and Optimal Control. Prentice-

Hall, Upper Saddle River, NJ.