Stability of Structured Hamiltonian Eigensolvers.
Various applications give rise to eigenvalue problems for which the matrices are Hamiltonian or skew-Hamiltonian and also symmetric or skew-symmetric.  We define structured backward errors that are useful for testing the stability of numerical methods for the solution of these four classes of structured eigenproblems. We introduce the symplectic quasi-QR factorization and show that for three of the classes it enables the structured backward error to be efficiently computed.  We also give a detailed rounding error analysis of some recently developed Jacobi-like algorithms of Fassbender, Mackey, and Mackey [Linear Algebra Appl., to appear] for these eigenproblems.  Based on the direct solution of 4  4, and in one case 8  8, structured subproblems these algorithms produce a complete basis of symplectic orthogonal eigenvectors for the two symmetric cases and a symplectic orthogonal basis for all the real invariant subspaces for the two skew-symmetric cases. We prove that, when the rotations are implemented using suitable formulae, the algorithms are strongly backward stable and we show that the QR algorithm does not have this desirable property.
Introduction
. This work concerns real structured Hamiltonian and skew-
Hamiltonian eigenvalue problems where the matrices are either symmetric or skew-
symmetric. We are interested in algorithms that are strongly backward stable for these
problems. In general, a numerical algorithm is called backward stable if the computed
solution is the true solution for slightly perturbed initial data. If, in addition, this
perturbed initial problem has the same structure as the given problem, then the
algorithm is said to be strongly backward stable.
There are three reasons for our interest in strongly backward stable algorithms.
First, such algorithms preserve the algebraic structure of the problem and hence
force the eigenvalues to lie in a certain region of the complex plane or to occur in
particular kinds of pairings. Because of rounding errors, algorithms that do not
respect the structure of the problem can cause eigenvalues to leave the required region
[26]. Second, by taking advantage of the structure, storage and computation can be
lowered. Finally, structure-preserving algorithms may compute eigenpairs that are
more accurate than the ones provided by a general algorithm.
Structured Hamiltonian eigenvalue problems appear in many scientific and engineering
applications. For instance, symmetric skew-Hamiltonian eigenproblems arise
in quantum mechanical problems with time reversal symmetry [9], [23]. In response
theory, the study of closed shell Hartree-Fock wave functions yields a linear response
eigenvalue equation with a symmetric Hamiltonian [21]. Also, total least squares
problems with symmetric constraints lead to the solution of a symmetric Hamiltonian
problem [17].
# Received by the editors February 23, 2000; accepted for publication (in revised form) by V.
Mehrmann November 24, 2000; published electronically May 3, 2001.
http://www.siam.org/journals/simax/23-1/36800.html
Department of Mathematics, University of Manchester, Manchester, M13 9PL, England
(ftisseur@ma.man.ac.uk, http://www.ma.man.ac.uk/-ftisseur). This work was supported by Engineering
and Physical Sciences Research Council grant GR/L76532.
The motivation for this work comes from recently developed Jacobi algorithms for
structured Hamiltonian eigenproblems [10]. These algorithms are structure-preserving,
inherently parallelizable, and hence attractive for solving large-scale eigenvalue prob-
lems. Our first contribution is to define and show how to compute structured backward
errors for structured Hamiltonian eigenproblems. These backward errors are useful
for testing the stability of numerical algorithms. Our second contribution concerns
the stability of these new Jacobi-like algorithms. We give a unified description of
the algorithms for the four classes of structured Hamiltonian eigenproblems. This
provides a framework for a detailed rounding error analysis and enables us to show
that the algorithms are strongly backward stable when the rotations are implemented
using suitable formulae.
The organization of the paper is as follows. In section 2 we recap the necessary
background concerning structured Hamiltonians. In section 3 we derive computable
structured backward errors for structured Hamiltonian eigenproblems. In section 4,
we describe the structure-preserving QR-like algorithms proposed in [5] for structured
Hamiltonian eigenproblems. We give a unified description of the new Jacobi-like
algorithms and detail the Jacobi-like update for each of the four classes of structured
Hamiltonian. In section 5 we give the rounding error analysis and in section 6 we
use our computable backward errors to confirm empirically the strong stability of the
algorithms.
2. Preliminaries. A matrix P # R 2n-2n is symplectic if P T
-In
In
I n is the n - n identity matrix.
A matrix H # R 2n-2n is Hamiltonian if Hamiltonian
matrices have the form
where E, F, G # R n-n and F G. We denote the set of real Hamiltonian
matrices by H 2n .
A matrix S # R 2n-2n is skew-Hamiltonian if
Skew-Hamiltonian matrices have the form
where E, F, G # R n-n and F are skew-symmetric. We denote the
set of real skew-Hamiltonian matrices by SH 2n .
Note that if H # H 2n , then P
SH 2n , where P is an arbitrary symplectic matrix. Thus symplectic similarities
preserve Hamiltonian and skew-Hamiltonian structure. Also, symmetric and skew-symmetric
structures are preserved by orthogonal similarity transformations. Therefore
structure-preserving algorithms for symmetric or skew-symmetric Hamiltonian
or skew-Hamiltonian eigenproblems have to use real symplectic orthogonal trans-
formations, that is, matrices U # R 2n-2n satisfying U T As
in [10], we denote by SpO(2n) the group of real symplectic orthogonal matrices.
Any U # SpO(2n) can be written as
I and
In

Tables

2.1 and 2.2, we summarize the structure of Hamiltonian and skew-
Hamiltonian matrices that are either symmetric or skew-symmetric, their eigenvalue
STABILITY OF STRUCTURED HAMILTONIAN EIGENSOLVERS 105

Table
Properties of structured Hamiltonian matrices H # H 2n .
Symmetric
real,
Skew-symmetric
pure imaginary,
pairs #, -

Table
Properties of structured skew-Hamiltonian matrices S # SH 2n .
Symmetric
real,
double # D 0
Skew-symmetric
pure imaginary,
double,
properties, and their symplectic orthogonal canonical form. We use D # R n-n to
denote a diagonal matrix and B # R n-n to denote a block-diagonal matrix that is the
direct sum of 1 - 1 zero blocks and 2 - 2 blocks of the form [ 0
-b
0 ]. These canonical
forms are consequences of results in [19].
Next, we show that the eigenvectors of skew-symmetric Hamiltonian matrices can
be chosen to have structure. This property is important when defining and deriving
structured backward errors.
Lemma 2.1. The eigenvectors of a skew-symmetric Hamiltonian matrix H can
be chosen to have the form [ z
-iz ] with z # C n .
Proof. Let
HU be the canonical form of H with
symplectic orthogonal. The matrix
-iI
I
iI ] is unitary and diagonalizes the
canonical form of H:
# .
Hence
is an eigenvector basis for H and this shows that the eigenvectors can be taken to
have the form [ z
-iz ] with z # C n .
Note that an eigenvector of a skew-symmetric Hamiltonian matrix does not necessarily
have the form [ z
-iz ]. For instance, consider
106 FRANC-OISE TISSEUR

Table
t: Number of parameters defining H.
Hamiltonian Skew-Hamiltonian
is an eigenvector of H, corresponding to the eigenvalue -id, that
is not of the form [ z
3. Structured backward error. We begin by developing structured backward
errors that can be used to test the strong stability of algorithms for our classes of
Hamiltonian eigenproblems.
3.1. Definition. For notational convenience, the symbol H denotes from now
on both Hamiltonian and skew-Hamiltonian matrices. Let (#x,
#) be an approximate
eigenpair for the structured Hamiltonian eigenvalue problem
R 2n-2n . A natural definition of the normwise backward error of an approximate
eigenpair is
where we measure the perturbation in a relative sense and # denotes any vector
norm and the corresponding subordinate matrix norm. Deif [8] derived the explicit
expression for the 2-norm
#x -H#x is the residual. This shows that the normwise relative backward
error is a scaled residual. The componentwise backward error is a more stringent
measure of the backward error in which the components of the perturbation #H are
measured individually:
#x, |#H| #|H| # .
Here inequalities between matrices hold componentwise. Geurts [12] showed that
1#i#2n
The componentwise backward error provides a more meaningful measure of the stability
than the normwise version when the elements in H vary widely in magnitude.
However, this measure is not entirely appropriate for our problems as it does not
respect any structure (other than sparsity) in H. Bunch [2] and Van Dooren [25] have
also discussed other situations when it is desirable to preserve structure in definitions
of backward errors.
The four classes of structured Hamiltonian matrices we are dealing with are defined
by t # real parameters that make up E and F (see Table 3.1). We write
this dependence as extend the
notion of componentwise backward error to allow dependence of the perturbations on
STABILITY OF STRUCTURED HAMILTONIAN EIGENSOLVERS 107
a set of parameters and they define structured componentwise backward errors. Following
their idea and notation we define the structured relative normwise backward
error by
-(#x,
implies that #H has the same structure as H. The
structured relative componentwise backward error #x,
#) is defined as in (3.1) but
with the constraint #H#F #H#F replaced by |#H| #|H|.
In our case, the dependence of the data on the t parameters is linear. We naturally
require (#x,
#) to have any properties forced upon the exact eigenpairs, otherwise
the backward error will be infinite. In the next subsections, we give algorithms for
computing these backward errors. We start by describing a general approach that
was used in [13] in the context of structured linear systems and extend it to the case
where the approximate solution lies in the complex plane.
3.2. A general approach for the computation of -(# x,
#). Let
and
-+i#. By equating real and imaginary parts, the constraint
#x
in (3.1) becomes
# u
or equivalently #H [
# u
. Applying the vec operator (which stacks the
columns of a matrix into one long vector), we obtain
where# denotes the Kronecker product. We refer to Lancaster and Tismenetsky [18,
Chap. 12] for properties of the vec operator and the Kronecker product. By linearity
we have
-t of full rank and where #p is the t-vector of parameters defining #H.
There exists a diagonal matrix D 1 depending on the structure of H (symmetric/skew-
symmetric Hamiltonian/skew-Hamiltonian) such that
# u
. Using (3.4) we can rewrite (3.3)
as Y BD using (3.5),
-(#x,
y
This shows that the structured backward error is given in terms of the minimal 2-norm
solution to an underdetermined system. If the underdetermined system is consistent,
then the minimal 2-norm solution is given in terms of the pseudo-inverse by
In this case
-(#x,
When H is a symmetric structured Hamiltonian, we can assume that
# and
# x are
real. Therefore
and from (3.2) we have [
Applying the vec operation gives
# u
I 2n # vec(#-I
As
-I - H is also a symmetric structured Hamiltonian, we have by linearity that
vec(#-I
#- is the t-vector of parameters defining
#- lies in the range of Y BD Therefore, the underdetermined
system in (3.6) is consistent for symmetric Hamiltonians and for symmetric skew-
Hamiltonians. For a skew-symmetric Hamiltonian, we can again prove consistency
for pure imaginary approximate eigenvalues and approximate eigenvectors of the form
in Lemma 2.1. We have not been able to prove that the underdetermined system is
consistent for the skew-symmetric skew-Hamiltonian case.
As the dependence on the parameters is linear, in the definition of the structured
relative componentwise backward error #x,
#), we have the equivalence
|#H| #|H| #p| #|p|.
q. Then the smallest # satisfying |#p| #|p| is
#q# . The minimal #-norm solution of Y BD 2 can be approximated by
minimizing in the 2-norm. We have
#).
By looking at each problem individually, it is possible to reduce the size of the
underdetermined system. Nevertheless, solution of the system by standard techniques
still takes O(n 3 ) operations. In the next section, we show that by using a symplectic
quasi-QR factorization of the approximate eigenvector and residual (or some appropriate
parts) we can derive expressions for -(#x,
#) that are cheaper to compute for all
the structured Hamiltonians of interest except for skew-symmetric skew-Hamiltonians.
First, we define a symplectic quasi-QR factorization.
3.3. Symplectic quasi-QR factorization. We define the symplectic quasi-QR
factorization of an 2n -m matrix A by
where Q is real symplectic orthogonal, T 1 # R n-m is upper trapezoidal, and T 2 #
R n-m is strictly upper trapezoidal. Such a symplectic quasi-QR factorization has
also been discussed by Bunse-Gerstner [3, Cor. 4.5(ii)].Before giving an algorithm to
STABILITY OF STRUCTURED HAMILTONIAN EIGENSOLVERS 109
compute this symplectic quasi-QR factorization, we need to describe two types of elementary
orthogonal symplectic matrices that can be used to zero selected components
of a vector.
A symplectic Householder matrix H # R 2n-2n is a direct sum of n-n Householder
matrices:
where
diag # I k-1 , I n-k+1 -v T v
I n otherwise,
and v is determined such that for a given x # R n , P (k,
A symplectic Givens rotation G(k, # R 2n-2n is a Givens rotation where the
rotation is performed in the plane (k, k G(k, #) has the form
where # is chosen such that for a given x # R 2n , G(k,
We use a combination of these orthogonal transformations to compute our symplectic
quasi-QR factorization: symplectic Householder matrices are used to zero large
portions of a vector and symplectic Givens are used to zero single entries.
Algorithm 3.1 (symplectic quasi-QR factorization). Given a matrix
with A 1 , A 2 # R n-m , this algorithm computes the symplectic quasi-QR factorization
(3.8).
For
End
Determine G
End
We illustrate the procedure for a generic 6 - 4 matrix:
-#
-#
G3
-#
3.4. Symmetric Hamiltonian eigenproblems. Let
be the residual vector and be the symplectic quasi-QR factorization (3.8)
with Q symplectic orthogonal and
. 0
e n+1,2
We have Q
which is equivalent to
e 110
e 22e n+1,2#
still a symmetric Hamiltonian matrix. Equation (3.10)
defines the first column of #
H. As |e 11
e 22 /e 11 , # - h
E) T
STABILITY OF STRUCTURED HAMILTONIAN EIGENSOLVERS 111
and #
be such that
e 11
e 22 -0 -
e 11
where the -'s are arbitrary real coe#cients. Then, any symmetric Hamiltonian of the
F
F -#
#x. The Frobenius norm of #H is minimized by setting the
-'s to zero in the definition of #
F . We obtain the following lemma.
Lemma 3.2. The backward error of an approximate eigenpair of a symmetric
Hamiltonian eigenproblem is given by
-(#x,
|e 11 |
is the quasi-triangular factor in the symplectic quasi-QR
factorization of [#x r] with
#I -H)#x. We also have
-(#x,
where e 2 is the second unit vector.
3.5. Skew-symmetric Hamiltonian eigenproblems. For skew-symmetric
Hamiltonian eigenproblems the technique developed in section 3.4 needs to be modified
as in this case r,
# x are complex vectors and we want to define a real skew-symmetric
Hamiltonian perturbation
so that (H
#x.
In the definition of the structured backward error (3.1), we now assume that
# is
pure imaginary and that
# x has the form [
(see Lemma 2.1). Taking the
plus sign in
# x, the equation (H
#x can be written as
E)#z.
Multiplying (3.12) by -i gives (3.11). Hence, we carry out the analysis with (3.11)
only. Setting
in (3.11) and equating real and imaginary
parts yields
which is equivalent to
Using
we show that w and s are orthogonal:
For the other choice of sign with
, the equation
#x
is equivalent to
and
we can show that w T
We can now carry on the analysis as in section 3.4. Let be the
symplectic quasi-QR factorization of [w s]. As w T we have that e
obtain #H by solving the underdetermined system
e 110
#e 22e n+1,2#
Lemma 3.3. The backward error of an approximate eigenpair (#x,
#) of a skew-symmetric
Hamiltonian eigenproblem with
pure imaginary and
x of the form
is given by
-(#x,
|e 11 |
is the quasi-triangular factor in the symplectic quasi-QR factorization
of [w s] with
if
We also have
-(#x,
where e 2 is the second unit vector.
3.6. Symmetric skew-Hamiltonian eigenproblems. The analysis for symmetric
skew-Hamiltonian eigenproblems is similar to that in section 3.4. The only
di#erence comes from noting that
F
STABILITY OF STRUCTURED HAMILTONIAN EIGENSOLVERS 113
using v T
#I-E)#x 1 . Instead of computing a symplectic
quasi-QR factorization of [#x r], we compute a symplectic quasi-QR factorization of
# x r] in order to introduce one more zero in the triangular factor R. We summarize
the result in the next lemma.
Lemma 3.4. The backward error of an approximate eigenpair of a symmetric
skew-Hamiltonian eigenproblem is given by
-(#x,
#) =|e 11 |
x r] is the quasi-triangular factor in the symplectic quasi-QR factorization
of [J
x r] with
#I -H)#x. We also have
-(#x,
#H#F .
3.7. Comments. Lemmas 3.2-3.4 provide an explicit formula for the backward
error that can be computed in O(n 2 ) operations.
For skew-symmetric skew-Hamiltonian matrices H, the eigenvectors are complex
with no particular structure. The constraint (H
#x in (3.1) can be written
in the form #H[#x), is the residual.
We were unable to explicitly construct matrices #H satisfying this constraint via a
symplectic QR factorization of [#x), #x), #(r), #(r)]. Thus, in this case, we have to
use the approach described in section 3.2 to compute -(#x,
#), which has the drawback
that it requires O(n 3 ) operations.
4. Algorithms for Hamiltonian eigenproblems. A simple but ine#cient
approach to solve structured Hamiltonian eigenproblems is to use the (symmetric or
unsymmetric as appropriate) QR algorithm on the 2n - 2n structured Hamiltonian
matrix. This approach is computationally expensive and uses 4n 2 storage locations.
Moreover, the QR algorithm does not use symplectic orthogonal transformations and
is therefore not structure-preserving.
Benner, Merhmann, and Xu's method [1] for computing the eigenvalues and invariant
subspaces of a real Hamiltonian matrix uses the relationship between the
eigenvalues and invariant subspaces of H and an extended 4n - 4n Hamiltonian ma-
trix. Their algorithm is structure-preserving for the extended Hamiltonian matrix
but is not structure-preserving for H. Therefore, it is not strongly backward stable
in the sense of this paper.
4.1. QR-like algorithms. Bunse-Gerstner, Byers, and Mehrmann [5] provide
a chart of numerical methods for structured eigenvalue problems, most of them based
on QR-like algorithms. In this section, we describe their recommended algorithms for
our structured Hamiltonian eigenproblems. In the limited case where rank(F
Byer's Hamiltonian QR algorithm [6] based on symplectic orthogonal transformations
yields a strongly backward stable algorithm.
For symmetric Hamiltonian eigenproblems, the quaternion QR algorithm [4] is
suggested. The quaternion QR algorithm is an extension of the Francis QR algorithm
for complex or real matrices to quaternion matrices. This algorithm uses exclusively
quaternion unitary similarity transformations so that it is backward stable. Compared
with the standard QR algorithm for symmetric matrices, this algorithm cuts the
storage and work requirements approximately in half. However, its implementation
requires quaternion arithmetic and it is not clear whether it is strongly backward
stable.
A skew-symmetric Hamiltonian H is first reduced via symplectic orthogonal transformations
to block antidiagonal form [ 0
-T
the blocks are symmetric tridi-
agonal. The complete solution is obtained via the symmetric QR algorithm applied
to T . The whole algorithm is strongly backward stable as it uses only real symplectic
orthogonal transformations that are known to be backward stable.
For symmetric skew-Hamiltonian problems, the use of the "X-trick" is suggested:
# with X =# 2
I I
-iI iI
# .
The eigenvalues of H are computed from the eigenvalue of the Hermitian matrices
using the Hermitian QR algorithm for instance. One drawback of
this approach is that it uses complex arithmetic and does not provide a real symplectic
orthogonal eigenvector basis. Hence the algorithm does not preserve the "realness"
of the original matrix.
Finally, for the skew-symmetric skew-Hamiltonian case, H is reduced to block-diagonal
form via a finite sequence of symplectic orthogonal transformations. The
blocks are themselves tridiagonal and skew-symmetric. Then Paardekooper's Jacobi
algorithm [22] or the algorithm in [11] for skew-symmetric tridiagonal matrices can
be used to obtain the complete solution. The whole algorithm is strongly backward
stable.
4.2. Jacobi-like algorithms. Byers [7] adapted the nonsymmetric Jacobi algorithm
[24] to the special structure of Hamiltonian matrices. The Hamiltonian Jacobi
algorithm based on symplectic Givens rotations and symplectic double Jacobi rotations
of the form
J# I 2n , where J is a 2 - 2 Jacobi rotation, preserves the Hamiltonian
structure. This Jacobi algorithm, when it converges, builds a Hamiltonian
Schur decomposition [7, Thm. 1]. For symmetric H, this Jacobi algorithm converges
to the canonical form [ D0
-D ] and is strongly backward stable. For skew-symmetric
Hamiltonian H, this Jacobi algorithm does not converge as the symplectic orthogonal
canonical form for H is not Hamiltonian triangular.
Recently, Fa-bender, Mackey, and Mackey [10] developed Jacobi algorithms for
structured Hamiltonian eigenproblems that preserve the structure and produce a complete
basis of symplectic orthogonal eigenvectors for the two symmetric cases and a
symplectic orthogonal basis for all the real invariant subspaces for the two skew-symmetric
cases. These Jacobi algorithms are based on the direct solution of 4 - 4,
and in one case 8 - 8, subproblems using appropriate transformations. The algorithms
work entirely in real arithmetic. Note that "realness" of the initial matrix can
be viewed as additional structure that these Jacobi algorithms preserve. We give a
unified description of these Jacobi-like algorithms for the four classes of structured
Hamiltonian eigenproblems under consideration.
Let H # R 2n-2n be a structured Hamiltonian matrix (see Table 2.1 and 2.2).
These Jacobi methods attempt to reduce the quantity (o#-diagonal norm)
STABILITY OF STRUCTURED HAMILTONIAN EIGENSOLVERS 115
where S is a set of indices depending on the structure of the problem using a sequence
of symplectic orthogonal transformations H # SHS T with S # R 2n-2n . The aim
is that H converges to its canonical form. In the following, we note A i,j,i+n,j+n the
restriction to the (i, n) plane of A.
Algorithm 4.1. Given a structured Hamiltonian matrix H # R 2n-2n and a
tolerance tol > 0, this algorithm overwrites H with its approximate canonical form
orthogonal and o#(PHP T
while
Choose (i,
Compute a symplectic orthogonal S
such that (SHS T ) i,j,i+n,j+n is in canonical form.
preserving structure
preserving structure
Note that the pair (i, uniquely determines a 4 - 4 principal submatrix
that also inherits the Hamiltonian or skew-Hamiltonian structure together with the
symmetry or skew-symmetry property. There are many ways of choosing the indices
(i, j) but this choice does not a#ect the rest of the analysis. We refer to n(n - 1)/2
updates as a sweep. Each sweep must be complete, that is, every part of the matrix
must be reached. We see immediately that any complete sweep of the (1, 1) block of
H consisting of 2-2 principal submatrices generates a corresponding complete sweep
of H.
For each 4 - 4 target submatrix, a symplectic orthogonal matrix that directly
computes the corresponding canonical form is constructed and embedded into the
in the same way that the 4 - 4 target has been extracted.
For skew-symmetric skew-Hamiltonians, the 4-4 based Jacobi algorithm does not
converge. The aim of these Jacobi algorithms is to move the weight to the diagonal
of either the diagonal blocks or o#-diagonal blocks. That cannot be done for a skew-symmetric
skew-Hamiltonian because these diagonals are zero. There is no safe place
where the norm of the target submatrix can be kept. However, if an 8 - 8 skew-symmetric
skew-Hamiltonian problem is solved instead, the 2 - 2 diagonal blocks of
H become a safe place for the norm of target submatrices and the resulting 8 - 8
based Jacobi algorithm is expected to converge. The complete sweep is defined by
partitioning blocks along the rightmost and
lower edges when n is odd. Hence, in this case we must also be able to directly solve
subproblems.
Immediately, we see that the di#cult part in deriving these algorithms is to define
the appropriate symplectic orthogonal transformation S that computes the canonical
form of the restriction to the (i, n) plane of H. Fa-bender, Mackey, and
Mackey [10] show that by using a quaternion representation of the 4 - 4 symplectic
orthogonal group, as well as 4 - 4 Hamiltonian and skew-Hamiltonian matrices in the
tensor square of the quaternion algebra, we can define and construct 4 - 4 symplectic
orthogonal matrices R that do the job. These transformations are based on rotations
of the subspace of pure quaternions.
We need to give all the required transformations in a form suitable for rounding
error analysis and also to facilitate the description of the structure preserving Jacobi
algorithms. We start by defining two types of quaternion rotations. This enables us
to encode the formulas in [10] into one. Let e s #= e 1 be a standard basis vector of R 4
and p # R 4 such that p #= 0, e T
(p is a pure quaternion), and p/#p# 2 #= e s . Let
s
# .
We define the left quaternion rotation by
# .
QL is symplectic orthogonal and not di#cult to compute. We have x
and the other components of x are just permutations of the coordinates of p.
We define the right quaternion rotation by
# .
The matrix QR is orthogonal. It is symplectic when s #= 3 and x
R 4 be nonzero. Following [10], we define the 4 - 4
symplectic orthogonal Givens rotation associated with p by
# .
We now have all the tools needed to define the symplectic orthogonal transformations
that directly compute the canonical form for each of the 4 - 4 structured
Hamiltonian eigenproblems of interest. We refer to [10] for more details about how
these transformations have been derived.
4.2.1. Symmetric Hamiltonian. Let H # R 4-4 be a symmetric Hamiltonian
matrix. The canonical form of H is obtained in two steps: first H is reduced to 2 - 2
block diagonal form and then the complete diagonalization is obtained by using a
double Jacobi rotation.
For the first step we consider the singular value decomposition of the 3-3 matrix
Let u 1 and v 1 be the left and right singular vectors corresponding to the largest
singular value # 1 and let
v1 ]. We have A T
STABILITY OF STRUCTURED HAMILTONIAN EIGENSOLVERS 117
so that e T
v, the vector x in (4.3) is such that
which implies that the right quaternion rotation QR (v, 2) is symplectic
orthogonal. As shown in [10], the product diagonalizes
H, that is, QHQ
E). Complete diagonalization is obtained by using
a double Jacobi rotation
is chosen such that
sin #
sin #
diagonalizes
In summary, the symplectic orthogonal transformation S used in Algorithm 4.1 is
equal to the identity matrix except in the (i, plane, where the (i, j, n
j)-restriction matrix is given by
4.2.2. Skew-symmetric Hamiltonian. Let H # R 4-4 be a skew-symmetric
Hamiltonian matrix and let p # R 4 be defined from the elements of H by
It is easy to verify that for
4.2.3. Symmetric skew-Hamiltonian. Let H # R 4-4 be a symmetric skew-
Hamiltonian matrix and let p # R 4 be defined from the elements of H by
diagonalizes H and
4.2.4. Skew-symmetric skew-Hamiltonian. For the convergence of the Jacobi
algorithm to be possible we need to solve an 8 - 8 subproblem. The matrix
H # R 8-8 is block diagonalized with three 4 - 4 symplectic Givens rotations of the
form (4.6) and one symplectic Givens rotation of the form (3.9). Let G be the product
of these rotations. We have
where
tridiagonal and skew-symmetric. The complete 2 - 2 block-
diagonalization is obtained by directly transforming
its real Schur form as
follows. In [20], Mackey showed that the transformation
directly computes the real Schur form of
E, that is,
(Q# I 2 )G is the symplectic orthogonal
transformation that computes the real Schur form of the 8 - 8 skew-symmetric skew-
Hamiltonian H:
# .
When n is odd, we have to solve a 6-6 subproblem for each complete sweep of the
Jacobi algorithm. As for the 8-8 case, the 6-6 skew-symmetric skew-Hamiltonian H
is first reduced to the form (4.7), where
tridiagonal and skew-symmetric.
This is done by using just one 4 - 4 symplectic Givens rotation followed by one
symplectic Givens rotation. Let
and
. Then computes directly the
real Schur form of
Moreover, we have e T
Q) and
# with
5. Error analysis of the Jacobi algorithms. In floating point arithmetic,
Algorithm 4.1 computes an approximate canonical form
T such that
where P is symplectic orthogonal, and an approximate basis of symplectic orthogonal
eigenvectors
P . We want to derive bounds for #H#
- I#, and
5.1. Preliminaries. We use the standard model for floating point arithmetic
[16]
where u is the unit roundo#. We assume that (5.1) holds also for the square roots
operation. To keep track of the higher terms in u we make use of the following result
[16, Lem. 3.1].
Lemma 5.1. If |# i | # u and #
We define
where p denotes a small integer constant whose value is unimportant. In the following,
computed quantities will be denoted by hats.
First, we consider the construction of a 4 - 4 Givens rotation and left and right
quaternion rotations.
Lemma 5.2. Let a 4 - 4 Givens rotation constructed according to
(4.6) with p # R 4 . Then the computed
G satisfies | #
G-G| # 5 |G|.
Proof. This result is a straightforward extension of Lemma 18.6 in [16] concerning
Givens rotations.
The rounding error properties of right and left quaternion rotations require more
attention. When p s < 0, the computation of #p# therefore the computation
of QL (p, s) or QR (p, s) is a#ected by cancellation. This problem can be overcome by
using another formula as shown in the next lemma.
Lemma 5.3. Let 4 - 4 left and right quaternion rotations
constructed according to
where
and
with p # R 4 given. Then the computed
QL and
| #
QL -QL | #
Proof. It is straightforward to verify that the expressions for QL (p, s) and QR (p, s)
in (5.2) and (5.3) agree with the definitions in (4.4) and (4.5).
We have f
As p s # 0, there exists # 5 such that f
using the same argument we have
We also have
f l
##
##
Using [16, Lem. 3.3] we have
and
Hence, we certainly have
In the following we use the term elementary symplectic orthogonal matrix to
describe any double Givens rotation, 4-4 Givens rotation, or left or right quaternion
rotation that is embedded as a principal submatrix of the identity matrix I # R 2n-2n .
We have proved that any computed elementary symplectic orthogonal matrix
used by the Jacobi algorithm satisfies a bound of the form
| #
Lemma 5.4. Let x # R 2n-2n and consider the computation of
Px, where
P is
a computed elementary symplectic orthogonal matrix satisfying (5.5). The computed
y satisfies
where P is the exact elementary symplectic orthogonal matrix.
Proof. The vector
y di#ers from x only in elements We have
We obtain similar results for
y n+i , and
y n+j . Hence,
As
Finally, we define
and note that #x#
Now, we consider the pre- and postmultiplication of a matrix H by an approximate
elementary symplectic orthogonal matrix
Lemma 5.5. Let H # R 2n-2n and P # R 2n-2n be any elementary symplectic
orthogonal matrix such that f l(P ) satisfies (5.5). Then,
f
STABILITY OF STRUCTURED HAMILTONIAN EIGENSOLVERS 121
Proof. Let h i be the ith column of H. By Lemma 5.4 we have
The same result holds for h j , h n+i , and h n+j and the other columns of H are un-
changed. Hence, f
Similarly,
B)P T with #
B#F . Then, with
with #
As a consequence of Lemma 5.5, if H k+1 is the matrix obtained after one Jacobi
update with S k (which is the product up to six elementary symplectic orthogonal
matrices), we have
is the exact transformation for
Up to now, we made no assumption on H. If H is a structured Hamiltonian
matrix, the (i, j)-restriction of RHR T is in canonical form. For instance,
if H is a skew-symmetric Hamiltonian matrix, in a computer implementation the
diagonal elements of H are not computed but are set to zero. Also, h ij , h i,j+n and by
skew-symmetry h ji , h j+n,i are set to zero. But by forcing these elements to be zero,
we are making the error smaller so the bounds still hold.
Because of the structure of the problem, both storage and the flop count can be
reduced by a factor of four. Any structured Hamiltonian matrix needs less than n 2 +n
storage locations. If only the t parameters defining H are computed, the structure in
the error is preserved and #H has the same structure as H. It is easy to see that
the bounds in Lemma 5.6 are still valid with the property that #H has the same
structure as H.
Theorem 5.6. Algorithm 4.1 for structured Hamiltonians H compute a canonical
T such that
where #H has the same structure as H and #H#F # k #H#F , where k is the
number of symplectic orthogonal transformations S i applied for each Jacobi update.
The computed basis of symplectic orthogonal eigenvectors
satisfies
Proof. From (5.6), one Jacobi update of H satisfies
For the second update we have
122 FRANC-OISE TISSEUR
Continuing
in this fashion, we find that, after k updates,
1 . S T
k with #H k #F #
In a similar way, using the first part of Lemma 5.5 we have
After k updates,
readily.
Theorem 5.6 shows that the computed eigenvalues are the exact eigenvalues of a
nearby structured Hamiltonian matrix and that the computed basis of eigenvectors is
orthogonal and symplectic up to machine precision. This proves the strong backward
stability of the Jacobi algorithms.
6. Numerical experiments. To illustrate our results we present some numerical
examples. All computations were carried out in MATLAB, which has unit roundo#
For symmetric Hamiltonians, symmetric skew-Hamiltonians, and skew-symmetric
Hamiltonians with approximate eigenvector
# x of the form [ z
-iz ], computing -(#x,
#) in
involves a symplectic quasi-QR factorization of a 2n - 2 matrix, which can be
done in order n 2 flops, a cost negligible compared with the O(n 3 ) cost of the whole
eigendecomposition.
For skew-symmetric Hamiltonians with approximate eigenvector
x not of the form
-iz ], and for skew-symmetric skew-Hamiltonians, the computation of -(#x,
#) requires
O(n 3 ) flops as we have to find the minimal 2-norm solution of a large underdetermined
system in (3.6). Thus, in this case, -(#x,
#) is not a quantity we would compute
routinely in the course of solving a problem.
Note that in our implementation of the Jacobi-like algorithm for skew-symmetric
Hamiltonians we choose the approximate eigenvectors to be the columns of P [ I
-iI
I
where P is the accumulation of the symplectic orthogonal transformations used by
the algorithm to build the canonical form. In this case, the approximate eigenvectors
x are guaranteed to be of the form [ z
To test the strong stability of numerical algorithms for solving structured Hamiltonian
eigenproblems, we applied the direct search maximization routine mdsmax of
the MATLAB Test Matrix Toolbox [15] to the function
1#i#2n
are the computed eigenpairs. In this way we carried out a search for
problems on which the algorithms performs unstably.
As expected from the theory, we could not generate examples for which the structured
backward error for the Jacobi-like algorithms is large: -(#x,
#) < nu#H#F in all
our tests.
The symmetric QR algorithm does not use symplectic orthogonal transformations
and is therefore not structure-preserving. To our surprise, we could not generate examples
of symmetric Hamiltonian and symmetric skew-Hamiltonian matrices for which
STABILITY OF STRUCTURED HAMILTONIAN EIGENSOLVERS 123

Table
Backward error of the eigenpair for of the 4-4 skew-symmetric Hamiltonian defined
by (6.1).
#) -max (#x,
#)
Jacobi-like algorithm

Table
Backward errors of the approximation of the eigenvalue 0 for a 30-30 random skew-symmetric
skew-Hamiltonian matrix.
|
#)
any of the eigenpairs computed by the symmetric QR algorithm has a large backward
error. However, the QR algorithm does not compute a symplectic orthogonal basis
of eigenvectors and also, it is easy to generate examples for which the -# structure
for symmetric Hamiltonians and eigenvalue multiplicity 2 structure for symmetric
skew-Hamiltonians is not preserved. If we generalize the definition of the structured
backward error of a single eigenpair to a set of k eigenpairs, the symmetric QR algorithm
is likely to produce sets of eigenpairs with an infinite structured backward
error. The QR-like algorithm for symmetric skew-Hamiltonians is likely to provide
eigenvectors that are complex instead of real, yielding an infinite structured backward
error in (3.14).
The good backward stability of individual eigenpairs computed by the QR algorithm
does not hold for the skew-symmetric Hamiltonian case. For instance, we
considered the skew-symmetric Hamiltonian eigenproblem
# , with
whose eigenvalues are distinct: In Table
6.1, we give the normwise, componentwise, and structured normwise backward error
of the eigenpair for computed by the unsymmetric QR algorithm and
the skew-symmetric Jacobi algorithm. The QR algorithm does not use symplectic
orthogonal transformations and the computed eigenvectors do not have the structure
-iz ]. Therefore, for the computation of - max (#x,
#), we use the general formula (3.7).
In the skew-symmetric skew-Hamiltonian case, when n is odd, 0 is an eigenvalue
of multiplicity two and is not always well approximated with the unsymmetric QR
algorithm. We generated a random 15 - 15 E and F . We give in Table 6.2 the
backward errors associated with the approximation of the eigenvalue 0 for both the
QR algorithm and Jacobi algorithm.
7. Conclusion. The first contribution of this work is to extend existing definitions
of backward errors in a way appropriate to structured Hamiltonian eigen-
problems. We provided computable formulae that are inexpensive to evaluate except
for skew-symmetric skew-Hamiltonians. Our numerical experiments showed that for
symmetric structured Hamiltonian eigenproblems, the symmetric QR algorithm computes
eigenpairs with a small structured backward error but the algebraic properties
of the problem are not preserved.
Our second contribution is a detailed rounding error analysis of the new Jacobi
algorithms of Fa-bender, Mackey, and Mackey [10] for structured Hamiltonian eigen-
problems. These algorithms are structure-preserving, inherently parallelizable, and
hence attractive for solving large-scale eigenvalue problems. We proved their strong
stability when the left and right quaternion rotations are implemented using our formulae
(5.2), (5.3). Jacobi algorithms are easy to implement and o#er a good alternative
to QR algorithms, namely, the unsymmetric QR algorithm, which we showed to be
not strongly backward stable for skew-symmetric Hamiltonian and skew-Hamiltonian
eigenproblems, and the algorithm for symmetric skew-Hamiltonians based on applying
the QR algorithm to (4.1), which does not respect the "realness" of the problem.

Acknowledgments

. I thank Nil Mackey for pointing out the open question
concerning the strong stability of the Jacobi algorithms for structured Hamiltonian
eigenproblems and for her suggestion in fixing the cancellation problem when computing
the quaternion rotations. I also thank Steve Mackey for his helpful comments
on an earlier manuscript.



--R

A new method for computing the stable invariant subspace of a real Hamiltonian matrix
The weak and strong stability of algorithms in numerical linear algebra
Matrix factorizations for symplectic QR-like methods
A quaternion QR algorithm
A chart of numerical methods for structured eigenvalue problems
A Hamiltonian QR algorithm
IEEE Trans.
A relative backward perturbation theorem for the eigenvalue problem

Hamilton and Jacobi come full circle: Jacobi algorithms for structured Hamiltonian eigenproblems
Accurately counting singular values of bidiagonal matrices and eigenvalues of skew-symmetric tridiagonal matrices
A contribution to the theory of condition
Backward error and condition of structured linear systems
Structured backward error and condition of generalized eigenvalue problems
The Test Matrix Toolbox for Matlab (version 3.0)
Accuracy and Stability of Numerical Algorithms
Oxford University Press
The Theory of Matrices
Canonical forms for Hamiltonian and symplectic matrices and pencils
Hamilton and Jacobi meet again: Quaternions and the eigenvalue problem
Solution of the large matrix equations which occur in response theory
An eigenvalue algorithm for skew-symmetric matrices

A Jacobi-like algorithm for computing the Schur decomposition of a non-Hermitian matrix
Structured linear algebra problems in digital signal processing
A symplectic method for approximating all the eigenvalues of a Hamiltonian matrix
--TR
quaternion rotation;structure-preserving;symplectic;jacobi algorithm;backward error;skew-Hamiltonian;hamiltonian;symmetric;rounding error;skew-symmetric