Symplectic Balancing of Hamiltonian Matrices.
We discuss the balancing of Hamiltonian matrices by structure preserving similarity transformations. The method is closely related to balancing nonsymmetric matrices for eigenvalue computations as proposed by Osborne [J. ACM, 7 (1960), pp. 338--345]and Parlett and Reinsch [Numer. Math., 13 (1969), pp. 296--304] and implemented in most linear algebra software packages. It is shown that isolated eigenvalues can be deflated using similarity transformations with symplectic permutation matrices. Balancing is then based on equilibrating row and column norms of the Hamiltonian matrix using symplectic scaling matrices. Due to the given structure, it is sufficient to deal with the leading half rows and columns of the matrix. Numerical examples show that the method improves eigenvalue calculations of Hamiltonian matrices as well as numerical methods for solving continuous-time algebraic Riccati equations.
Introduction
. The eigenvalue problem for Hamiltonian matrices
A G
where A; G; Q 2 R n\Thetan and G, Q are symmetric, plays a fundamental role in many
algorithms of control theory and other areas of applied mathematics as well as computational
physics and chemistry. Computing the eigenvalues of Hamiltonian matrices is
required, e.g., when computing the H1 -norm of transfer matrices (see, e.g., [9, 10]),
calculating the stability radius of a matrix ([13, 29]), computing response functions
[22], and many more. Hamiltonian matrices are also closely related to continuous-time
algebraic Riccati equations (CARE) of the form
with A; G; Q as in (1.1) and X 2 R n\Thetan is a symmetric solution matrix. Many numerical
methods for solving (1.2) are based on computing certain invariant subspaces
of the related Hamiltonian matrices; see, e.g., [19, 21, 26, 28]. For a detailed discussion
of the relations of Hamiltonian matrices and continuous-time algebraic Riccati
equations (1.2) we refer to [18].
In eigenvalue computations, matrices and matrix pencils are often preprocessed
using a balancing procedure as described in [23, 25] for a general matrix A 2 R n\Thetan .
First, A is permuted via similarity transformations in order to isolate eigenvalues, i.e.,
a permutation matrix P 2 R n\Thetan is computed such that
are upper triangular matrices. Then, a diagonal
matrix
Universitat Bremen, Fachbereich 3 - Mathematik und Informatik, Zentrum fur Technomathe-
matik, 28334 Bremen, Germany. E-mail: benner@math.uni-bremen.de
is computed such that rows and columns of D \Gamma1
Z ZDZ are as close in norm as possible.
That is, balancing consists of a permutation step and a scaling step. In the scaling
step, the rows and columns of a matrix are scaled, which usually leads to a decrease
of the matrix norm. This preprocessing step often improves the accuracy of computed
eigenvalues significantly; isolated eigenvalues (i.e., those contained in T 1 and T 2 ) are
even computed without roundoff error.
Unfortunately, applying this balancing strategy to a Hamiltonian matrix H as
given in (1.1) will in general destroy the Hamiltonian structure. This is no problem if
the subsequent eigenvalue algorithm does not preserve or use the Hamiltonian struc-
ture. But during the past fifteen years, several structure preserving methods for the
Hamiltonian eigenproblem have been suggested. In particular, the square-reduced
method [31], the Hamiltonian QR algorithm (if in (1.1), rank
[12], the recently proposed algorithm based on a symplectic URV-like decomposition
[7], or the implicitly restarted symplectic Lanczos method of [5] for large sparse Hamiltonian
eigenproblems are appropriate choices for developing subroutines for library
usage and raise the need for a symplectic balancing routine. Similarity transformations
by symplectic matrices preserve the Hamiltonian structure. Thus, in order to
balance a Hamiltonian matrix and to preserve its structure, the required permutation
matrix and the diagonal scaling matrix should be symplectic.
In Section 2 we will give some necessary background. Isolating eigenvalues of
Hamiltonian matrices without destroying the structure can be achieved using symplectic
permutation matrices. This will be the topic of Section 3. How to equilibrate
rows and norms of Hamiltonian matrices in a similar way as proposed in [25] using
symplectic diagonal scaling matrices will be presented in Section 4. When invariant
subspaces, eigenvectors, or solutions of algebraic Riccati equations are the target of
the computations, some post-processing steps are required. These and some other
applications of the proposed symplectic balancing method are discussed in Section 5.
Some numerical examples on the use of the proposed balancing strategy for eigenvalue
computation and numerical solution of algebraic Riccati equations are given in
Section 6.
2. Preliminaries. The following classes of matrices will be employed in the
sequel.
Definition 2.1. Let
where I n is the n \Theta n identity matrix.
a) A matrix H 2 R 2n\Theta2n is Hamiltonian if (HJ) . The Lie Algebra of
Hamiltonian matrices in R 2n\Theta2n is denoted by H 2n .
b) A matrix H 2 R 2n\Theta2n is skew-Hamiltonian if (HJ) . The Jordan
algebra of skew-Hamiltonian matrices in R 2n\Theta2n is denoted by SH 2n .
c) A matrix S 2 R 2n\Theta2n is symplectic if SJS
The Lie group of symplectic matrices in R 2n\Theta2n is denoted by S 2n .
d) A matrix U 2 R 2n\Theta2n is unitary symplectic if U 2 S 2n and UU . The
compact Lie group of unitary symplectic matrices in R 2n\Theta2n is denoted by US 2n .
Observe that every H 2 H 2n must have the block representation given in (1.1).
In [11], an important relation between symplectic and Hamiltonian matrices is
proved.
Proposition 2.2. Let S 2 R 2n\Theta2n be nonsingular. Then S \Gamma1 HS is Hamiltonian
for all H 2 H 2n if and only if S T
This result shows that in general, similarity transformations that preserve the
Hamiltonian structure have to be symplectic up to scaling with a real scalar.
The following proposition shows that the structure of 2n \Theta 2n orthogonal symplectic
matrices permits them to be represented as a pair of n \Theta n matrices. Hence,
the arithmetic cost and storage for accumulating orthogonal symplectic matrices can
be halved.
Proposition 2.3. [24] An orthogonal matrix U 2 R 2n\Theta2n is symplectic if and
only if it takes the form
2.
We have the following well-known property of the spectra of Hamiltonian matrices
(see, e.g., [18, 21], and the references given therein).
Proposition 2.4. The spectrum of a real Hamiltonian matrix, denoted by oe (H)
is symmetric with respect to the imaginary axis, i.e., if  2 oe (H), then \Gamma 2 oe (H).
The spectrum of Hamiltonian matrices can therefore be partitioned as
oe
When solving Hamiltonian eigenproblems one would like to compute a Schur
form for Hamiltonian matrices analogous to the real Schur form for non-symmetric
matrices. This should be done in a structure-preserving way.
Definition 2.5. a) Let "
H has the form
G
A T
A 2 R n\Thetan is in real Schur form (quasi-upper triangular) and "
is real Hamiltonian quasi-triangular.
b) If H 2 H 2n and there exists U 2 US 2n such that "
real Hamiltonian
quasi-triangular, then "
H is in real Hamiltonian Schur form and U T HU is called
a Hamiltonian Schur decomposition.
If a Hamiltonian Schur decomposition exists such that "
H is as is (2.2), then U
can be chosen such that oe ( "
Most of the structure-preserving methods for the Hamiltonian eigenproblem, i.e.,
those using symplectic (similarity) transformations, rely on the following result. For
Hamiltonian matrices with no purely imaginary eigenvalues this result was first stated
in [24] while in its full generality as given below it has been proved in [20].
Theorem 2.6. Let H 2 H 2n and let iff its pairwise distinct nonzero
purely imaginary eigenvalues. Furthermore, let the associated H-invariant subspaces
be spanned by the columns of U k , Then the following are equivalent.
There exists S 2 S 2n such that S \Gamma1 HS is real Hamiltonian quasi-triangular.
ii) There exists U 2 US 2n such that U T HU is in real Hamiltonian Schur form.
is congruent to J for all is always of the
appropriate dimension.
Note that from Theorem 2.6 it follows that purely imaginary eigenvalues of H 2
must have even algebraic multiplicity in order for the Hamiltonian Schur
form of H to exist.
3. Isolating Eigenvalues by Symplectic Permutations. Let P denote any
n \Theta n permutation matrix. It is easy to see that symplectic permutation matrices
have the form
With matrices of type (3.1) it is possible to transform a Hamiltonian matrix to the
~
22 0
-z -z -z -z -z -z
where A 11 , A 33 are upper triangular and either Q The existence of
such a P s will be proved in a constructive way later by Algorithm 3.4 which transforms
a Hamiltonian matrix to the form given in (3.2). From a Hamiltonian matrix having
the form (3.2) a total of 2(p of H can be read off directly as seen by
the following result.
Lemma 3.1. Let H 2 H 2n and is of the
form (3.2) and either G there exists a permutation matrix
2n\Theta2n such that
are upper triangular with
and
A 22 G 22
is a 2r \Theta 2r Hamiltonian submatrix of H.
Proof. Let H be as in (3.2) and
I
Then
13 A 11 A 12 G 12 G 11 A 13
22 \GammaA T
diag (\Pi q ; I
I
I r 0
I r
I q 07 7 7 7 7 7 5
Thus,
has the desired form. The eigenvalue relation (3.4) follows
from
\GammaI q 0
\GammaI q 0
13 A 11
Lemma 3.1 is merely of theoretical interest and demonstrates that in order to solve the
Hamiltonian eigenvalue problem, we can proceed by working only with H 22 . But the
transformations we have used in the proof are in general non-symplectic. If we want to
compute invariant subspaces, eigenvectors, and/or the Hamiltonian Schur form given
in Theorem 2.6, we can transform the Hamiltonian matrix in (3.2) such that it has
Hamiltonian Schur form in rows and columns
But this can not be accomplished using only symplectic permutation matrices of the
form (3.1). Therefore we need another class of transformation matrices.
Definition 3.2. A matrix P J 2 R 2n\Theta2n is called a J-permutation matrix if
a) it is symplectic, i.e., P T
c) each row and column have exactly one nonzero entry.
As P J 2 US 2n , it is clear that a similarity transformation by a J-permutation
matrix preserves the Hamiltonian structure. In analogy to standard permutations,
similarity transformations with P J can be performed without floating point opera-
tions. Moreover, they can be represented by a signed integer vector IP of length n,
rows and columns k; j are to be
interchanged while the sign indicates if the corresponding entry in P J is +1 or \Gamma1.
The entries of P J in rows to 2n can be deduced from IP and Proposition 2.3.
Furthermore, symplectic permutation matrices as given in (3.1) are J-permutation
matrices.
Lemma 3.3. For any H 2 H 2n having the form (3.2), there exists a J-permutation
matrix P J such that
A 11
A 12
G 11
A 22
G T"
G 22
A T
A T
A T3
A 11 is upper triangular and with the notation in (3.2),
Proof. Let a Hamiltonian matrix H be given as in (3.2). We need a J-permutation
matrix only in the first step. Let
I
Obviously, P 1 is a J-permutation matrix and
A 11 A 12 \GammaG 13 G 11 G 12 A 13
22
Now assume G
~
I q
Then,
13 A T0 A 11 A 12 A 13 G 11 G 12
We thus obtain the form (3.9) by another similarity transformation with
where \Pi q is defined in (3.8). For the other case, i.e.,
~
I q
and
In both cases, "
J-permutation matrix and "
P is a
Hamiltonian matrix having the desired form (3.9).
In order to isolate eigenvalues, it is sufficient to restrict ourselves to symplectic
permutations. But having computed the form (3.2), it is possible that there are still
isolated eigenvalues in H 22 . Applying the same procedure used to isolate eigenvalues
in H to H 22 , we can transform H 22 to the form (3.2). This process can then be
repeated until no more isolated eigenvalues are found. Accumulating all permutations
in a symplectic permutation matrix P s of the form (3.1), this results in a similarity
transformation
~
. A
. 0
. 0
gps
gps
Here, A j;j t, are upper triangular and for
either and the Hamiltonian submatrix
A t;t G t;t
has no isolated eigenvalues. If we now define p :=
then we
can partition ~
H in (3.10) as in (3.2). Then the first step in the proof of Lemma 3.1
can be performed to obtain the form (3.7). Just the block-structure of the upper left
and lower right diagonal blocks in (3.7) are more complicated. But it is still possible
to bring them to upper triangular form using repeatedly the same sequence of permutations
used in the proof of Lemma 3.1. This shows that 2(p + q) eigenvalues of the
Hamiltonian matrix can be read off directly from (3.10) and that ~
H is permutationally
similar to4
~
~
Y ~
Z
Further, we have H t;t
s
s
If only eigenvalues are required, we can continue working only with H t;t . If also eigen-vectors
and/or invariant subspaces are required, the similarity transformations used
to solve the reduced-order eigenproblem for H t;t have to be applied to the whole matrix
~
H . In that case, ~
H should be transformed to the form given in (3.9). Partitioning
~
H from (3.10) as in (3.2), we can perform the first step of the proof of Lemma 3.3
with the J-permutation matrix P 1 . The subsequent steps to achieve upper triangular
form in the first p + q rows and columns have then to be performed for each of the
first distinguishing the cases Q
A procedure to transform a Hamiltonian matrix H to the form in (3.10) is given in
the following algorithm. Note that in the given algorithm,
t.
Algorithm 3.4.
Input: Matrices A; G; Q 2 R n\Thetan , where defining a Hamiltonian
Output: A symplectic permutation matrix P s ; matrices A; G; Q with
defining a Hamiltonian matrix
having the form
~
END IF
END WHILE
~
END IF
END WHILE
END WHILE
END
In each execution of the outer WHILE-loop, we first search a row isolating an
eigenvalue. If such a row is found, we look for a column isolating an eigenvalue. In
this fashion it can be guaranteed that at the end, there are no more isolated eigenvalues
although we always only touch the first n rows and columns of the Hamiltonian matrix.
In an actual implementation one would of course never form the permutation
matrices explicitly but store the relevant information in an integer vector. Multiplications
by permutation matrices are realized by swapping the data contained in the
rows or columns to be permuted; for details, see, e.g., [3].
It is rather difficult to give a complete account of the cost of Algorithm 3.4. If
there are no isolated eigenvalues, the algorithm requires 4n floating point additions
and 2n comparisons as opposed to 8n additions and 4n comparisons for the
unstructured permutation procedure from [25] as implemented in the LAPACK subroutine
xGEBAL [3] when applied to H 2 R 2n\Theta2n . The worst case for Algorithm 3.4
would be that in each execution of the outer WHILE-loop, an isolated eigenvalue is
found in the last execution of the second inner WHILE-loop. In that case, the cost consists
of 4n 3 =3 floating point additions, moving
floating point numbers. But in this worst-case analysis, all eigenvalues are
isolated such that after permuting, there is nothing left to do, and the Hamiltonian
matrix is in Hamiltonian Schur form. A worst-case study for xGEBAL shows that
the permutation part requires 8n 3
and moving floating point numbers. We can therefore conclude that Algorithm
3.4 is about half as expensive as the procedure proposed in [25] applied to a
Hamiltonian matrix.
4. Symplectic Scaling. Suppose now that we have transformed the Hamiltonian
matrix to the form (3.10). Since all subsequent transformations are determined
from H t;t , the scaling parameters to balance H t;t have now to be chosen such that
the rows and columns of H t;t (instead of ~
are as close in norm as possible.
In order to simplify notation we will in the sequel call the Hamiltonian matrix
again H . Let H off be the off-diagonal part of H , i.e.,
We may without loss of generality assume that none of the rows and columns of H off
vanishes identically. Otherwise, we could isolate another pair of eigenvalues.
Now we want to scale H such that the norms of its rows and columns are close
in norm. As noted before, employing the technique of Parlett and Reinsch [25] destroys
the Hamiltonian structure. Diagonal scaling has thus to be performed using a
symplectic diagonal matrix D s . Such a matrix must have the form,
where D 2 R n\Thetan is a nonsingular diagonal matrix.
Let us at first note an obvious result for Hamiltonian matrices. Here and in the
sequel we will use the colon notation (see, e.g., [15]) H(:; k), H(j; :) to indicate the
kth column and jth row, respectively, of a matrix H .
Lemma 4.1. Let H 2 R 2n\Theta2n be a Hamiltonian matrix. Then for all p  1 and
for all
i.e., the p-norms of the ith column equals the norm of the (n + i)th row and the norm
of the ith row equals the norm of the (n
Proof. The result is obvious by noting kxk
for x 2 R 2n and
observing that from the structure of Hamiltonian matrices, we have
and furthermore, Equation (4.3) follows analogously by
noting
We can now conclude that it is sufficient to equilibrate the norms of the first
rows and columns of a 2n \Theta 2n Hamiltonian matrix by using a consequence of
Lemma 4.1.
Corollary 4.2. Let H 2 R 2n\Theta2n be a Hamiltonian matrix. Then for all p  1
and for all
Since a similarity transformation with any diagonal matrix does not affect the
diagonal elements of the transformed matrix, it is in the following sufficient to consider
H off . We will employ the notation
ith column of H off ;
transpose of the ith row of H off :
In the sequel, we will for convenience use 1. The results also hold for any
other p-norm. From a computational point of view it is also reasonable to use the
1-norm, since its computation does not involve any floating point multiplications and
furthermore, reducing the norm of H in one norm usually implies also a reduction in
the other norms.
Equilibrating can now be achieved in a similar way as in the
Parlett/Reinsch method. If fi denotes the base of the floating point arithmetic and
oe i is any signed integer, then they compute fi oe i closest to the real scalar
Thus, with D
H is in general no longer Hamiltonian. Unfortunately, using the
symplectic diagonal matrix
where
i and computing
~
(D (i)
s
we obtain
and thus in general, k ~
Nevertheless, equilibrating the 1-norms of h i and h i can be achieved by requiring
solving the resulting quartic equation
jg ii j:
It remains to show that equation (4.6) has a positive solution.
Theorem 4.3. Let H 2 R 2n\Theta2n be a Hamiltonian matrix and denote its off-diagonal
part by H off . Assume that none of the rows and columns of H off vanishes
identically. Then there exists a unique real number such that for ~
H as in (4.4)
we have
Proof. Solutions of Equation (4.6) are non-zero solutions of
k=0 a k t k . Recalling that g ii = the coefficients
of the polynomial p satisfy
a
a
Since there is at most one change of sign in the coefficients of the polynomial p,
Descartes' rule of signs shows that there is at most one positive zero of p. So if
there exists a positive solution of (4.6), it is unique. By assumption, kh
Therefore, either a ii is part of h i ) and either a 3 ? 0
or a 4 ? 0 (as q ii is part of h i ). Thus, we know that p is a polynomial of degree at
least 3 with positive leading coefficient and hence, lim t!1
On the other hand, if using the mean value theorem it
follows that there exists a positive zero of p. If g ii = 0, then is a zero of p and
The third order polynomial has a positive zero because
of the mean value theorem and
by has again at least one positive zero and hence equation (4.6) has
at least one positive real solution regardless of the value of g ii . On the other hand,
it was already observed that there is at most one such solution and we can conclude
that there exists a unique solving equation (4.6) whence k ~ h
The other equalities follow immediately from Corollary 4.2.
Computing the exact value equilibrating the ith and (n+i)th rows and columns
would require the solution of the fourth-order equation (4.6). Since the diagonal
similarity transformations are to be chosen from the set of machine numbers, it is
sufficient to find the machine number fi oe i closest to ffi i . This can be done similarly
to the computation in the general case as proposed in [25] and implemented in the
Fortran 77 subroutines BALANC from EISPACK [14] or its successor xGEBAL from
LAPACK [3]. That is, starting with the quantities in (4.5) are evaluated and
compared. If k ~ h
then this is repeated for
then we use
. This is achieved by the following algorithm.
Algorithm 4.4.
Input: Hamiltonian matrix H 2 H 2n having no isolated eigenvalues, base
of floating point arithmetic.
Output: Diagonal matrix D s 2 S 2n , H is overwritten by D \Gamma1
row and column norms equilibrated as far as possible.
r
END WHILE
r
END WHILE
END IF
diag
diag
diag
END
One execution of the outer FOR-loop of Algorithm 4.4 can be considered as a
sweep. The algorithm is terminated if for a whole sweep, all D . Usually, the
row and column norms are approximately equal after very few sweeps. Afterwards, the
iteration makes only very limited progress. Therefore, Parlett and Reinsch propose
in [25] a modification, which, translated to our problem, becomes:
Let ffi i be determined by the two inner WHILE-loops of Algorithm 4.4 and
compute
jg ii j:
If (where fl is a given positive constant),
then compute D i as in Algorithm 4.4. Otherwise, set D
For the behavior is essentially the same as for Algorithm 4.4 (in a few
cases, Algorithm 4.4 increases kh which can not happen if 1). For
slightly smaller than one, a step is skipped if it would produce an insubstantial
reduction of kh
In an actual computation, the similarity transformations with the D i 's can be
applied directly to the blocks A, G, and Q of the Hamiltonian matrix without forming
the Hamiltonian matrix itself. Thus, each similarity transformation can be performed
using only 4n \Gamma 4 multiplications. When the standard (not structure preserving)
scaling procedure from [25] is applied to H , each similarity transformation requires
multiplications. (Recall that in Algorithm 4.4, two rows and columns are
equilibrated at a time while only one row and column is treated in each step of the
inner FOR-loop of the standard procedure.)
The number of sweeps required to converge is similar to those for the general case
since the theory derived in [25] only requires the assumption of similarity transformations
with diagonal matrices and that in step i of each sweep, the ith rows and
columns are equilibrated as far as possible with . But this is accomplished by
Algorithm 4.4. Moreover, if ffi i is taken as the exact solution of (4.6), the convergence
of the sequence of similarity transformations to a stationary point can be proved as
in [23, 16]. That is, if
i is the solution of (4.6) in sweep k, then lim k!1 ffi (k)
for all hence in the limit, H is a balanced Hamiltonian matrix.
Note that here, each sweep has length n while in the standard balancing algorithm,
one has to go through each row/column pair of the matrix and thus, each sweep has
length 2n. Thus, the computational cost for scaling a 2n \Theta 2n Hamiltonian matrix
by Algorithm 4.4, assuming k 1 sweeps are required, is 4n n) as opposed
to 8n 2 n) for the standard scaling procedure as given in [25] with assumed
sweeps required for convergence. In general, k 1  k 2 such that the structure-preserving
scaling strategy is about half as expensive as the standard procedure.
These flop counts are based on the assumption that the cost for determining the
can be considered as small (O(1)) compared to the similarity transformations.
Remark 4.5. In [17] it is proposed to solve the matrix balancing problem using
a convex programming approach. To compare the complexity of this approach to
that of Algorithm 4.4, suppose that Algorithm 4.4 terminates after k 1 sweeps with
. For the matrix H to be balanced, let
ng \Theta ng
Eg. Assume that
and  ! e . (Here, Theorem 5 in [17] states that the complexity
of computing a diagonal matrix Y with positive diagonal entries such that the
rows and columns of Y \Gamma1 HY are balanced with the same accuracy as achieved by
Algorithm 4.4 is O
ne
hmin
. From numerical experience, it can be assumed
that with respect to n. Hence, Algorithm 4.4 can be considered
to be of complexity O(n 2 ). This complexity is clearly superior to that of the convex
programming approach which is still the case if k
Algorithm 4.4 requires a careful implementation to guard against over- and underflow
due to a very large/small ffi i . Here, we can use the bounds discussed in [25] and
implemented in LAPACK subroutine xGEBAL [3]; we just have to take into account
that in each step we scale by fi \Sigma2 rather than fi as in xGEBAL.
5. Backtransformation, Ordering of Eigenvalues, and Applications. So
far we have only considered the problem of computing the eigenvalues of a Hamiltonian
matrix. In order to compute eigenvectors, invariant subspaces, and the solutions of
algebraic Riccati equations, we have to transform the Hamiltonian matrix to real Schur
form. As we are considering structure-preserving methods, the goal is to transform
the Hamiltonian matrix to real Hamiltonian Schur form as given in Theorem 2.6 a)
- if it exists.
Assume that we have applied Algorithm 3.4 to the Hamiltonian matrix and obtained
a symplectic permutation matrix P s such that P T
s HP s has the form given in
(3.10). Then, we have applied a J-permutation P J to the permuted Hamiltonian
matrix such that its rows and columns numbered
are in Hamiltonian Schur form, i.e., P T
has the form given in (3.9). (From
Lemma 3.3 we know that such a P J exists.) Next, we have applied Algorithm 4.4 to
the Hamiltonian submatrix H t;t 2 H 2r from (3.11) and obtained a diagonal matrix
diag
I
Then
A 22
G T"
G 22
A T
A T
A T3
A 11 2 R (p+q)\Theta(p+q) is upper triangular and the Hamiltonian submatrix "
H 22 :=
A22
Q22
G22
A
has no isolated eigenvalues and its rows and columns are equilibrated by
Algorithm 4.4. Now assume the Hamiltonian Schur form of "
H 22 exists and we have
computed U
U22
V22
V22
U22
US 2r that transforms "
H 22 into real Hamiltonian Schur
form. Set
I
22
real Hamiltonian quasi-triangular and
. The first n columns of S span a Lagrangian H-invariant subspace. In most
applications, the c-stable H-invariant subspace is desired. Let us assume the method
used to transform "
H 22 to Hamiltonian Schur form chooses U 22 such that the first r
columns of U 22 , i.e., the columns of
U22
V22
, span the "
H 22 -invariant subspace of choice.
But there is no guarantee that the isolated eigenvalues in "
A 11 are the desired ones.
In that case, we have to reorder the Hamiltonian Schur form in order to move the
undesired eigenvalues to the lower right block of "
H and the desired ones to the upper
left block. Assume that we want to compute the Lagrangian H-invariant subspace
corresponding to a set  ae oe (H) which is closed under complex
conjugation. (Note that this is a necessary condition in order to obtain a Lagrangian
invariant subspace [2]). Using the standard reordering algorithm for the real Schur
form of an n \Theta n unsymmetric matrix as given in [15, 30], we can find an orthogonal
matrix ~
U   such that with the orthogonal symplectic matrix U  diag
U   ; ~
U
, we
have that
~
A 11
~
A 12
~
G 11
~
A 22
~
G T~
G 22
A T
A T
A T3
where A 11 ; A 22 are quasi-upper triangular, and
A T
22
A T
Therefore, we have to swap the eigenvalues in ~
A 22 and \Gamma ~
A T
22 . Note that the eigenvalues
to be re-ordered are among the isolated eigenvalues and hence are real. This implies
that ~
A 22 is upper triangular. The re-ordering can be achieved analogously to the
re-ordering of eigenvalues in the real Schur form as given in [15, 30]. The following
procedure uses this standard re-ordering in order to swap eigenvalues within ~
A 22
A T
22 ) and requires rotations working exclusively in rows and columns n and
2n in order to exchange eigenvalues from ~
A 22 with eigenvalues from \Gamma ~
A T
22 . Assume
A 22
A T 22
sn
\Gammas n
cn
be a
Givens rotation matrix that annihilates the second component of
gnn
where
~
Then U n is a symplectic Givens rotation matrix acting in planes n and 2n and
~
A 11
~
A 12  a 1;n
~
G 11
~
A 22
. ~
G T~
G 22
gnn
A T
A T
A T
22 0
Here, the bar indicates elements changed by the similarity transformation.
The next step is now to move n up in the upper diagonal block using again the
standard ordering subroutine such that we obtain again the form given in (5.1), just
A 11 2 R (n\Gammak+1)\Theta(n\Gammak+1) and again, the relations (5.2) hold. This procedure has now
to be repeated until
A 11 ).
Remark 5.1. If the Hamiltonian matrix has the form
h A
which
corresponds to a linear system
stabilizable
and (C; detectable, then each isolated eigenvalue in (5.1) given by the diagonal
elements of ~
A 11 has negative real part. Otherwise, these eigenvalues are unstable or
undetectable and can not be stabilized/detected. Therefore, if we have not mixed
up blocks by the J-permutation matrix P J (i.e., in Algorithm 3.4, i n) and the
c-stable H-invariant subspace is required, no re-ordering is necessary.
Remark 5.2. When solving algebraic Riccati equations using any approach based
on the Hamiltonian eigenproblem, the symplectic balancing strategy proposed here is
often not enough to minimize errors caused by ill-scaling. This is due to the effect that
for a balanced Hamiltonian matrix
h A
G
we still may have kQk AE kGk which
may cause large errors when computing invariant subspaces [27]. Therefore, another
symplectic scaling using a similarity transformation with diag
ae I n
ae 2 R, should be applied to H in order to achieve kAk  kGk  kQk as far as possible;
see [4] for details and a discussion of several heuristic strategies to achieve this.
Remark 5.3. Everything derived so far for Hamiltonian matrices can be applied
in the same way to skew-Hamiltonian matrices. If N 2 SH 2n , then
h A
G
A T
with . The skew-Hamiltonian structure is again preserved
under symplectic similarity transformations. Hence, isolating eigenvalues, re-ordering,
etc., can be achieved in the same way as for Hamiltonian matrices as all considered
transformations do not depend on the signs in the matrix blocks A, G, Q, but only on
the distinction zero/non-zero when isolating eigenvalues and on the absolute values
of the entries when equilibrating rows and norms. Note that Algorithm 4.4 even
simplifies quite a lot for real skew-Hamiltonian matrices. As q
can be computed as in the general balancing algorithm
for non-symmetric matrices because in (4.5) we obtain k ~
Eigenvalues of skew-Hamiltonian matrices as well as a skew-Hamiltonian Schur
form can be computed in a numerically strong backward stable way by Van Loan's
method [31]. It is advisable to balance skew-Hamiltonian matrices using the proposed
strategies prior to applying this algorithm.
Remark 5.4. We have considered so far only real Hamiltonian and skew-
Hamiltonian matrices. Isolating eigenvalues and equilibrating rows and columns
for complex (skew-)Hamiltonian matrices can be achieved in exactly the same way.
A structure-preserving, numerically backward stable (and hence numerically strong
backward stable) method for solving the complex (skew-)Hamiltonian eigenproblem
has recently been proposed [8]. The proposed symplectic balancing method can (and
should) also be used prior to applying this algorithm.
6. Numerical Examples. We have tested the symplectic balancing strategy
for eigenvalue computations. The computations were done in Matlab 1 Version 5.2
with machine precision "  2:2204 \Theta 10 \Gamma16 . Algorithms 3.4 and 4.4 were implemented
as Matlab functions. We used the modified algorithm as suggested by (4.8) where
we set suggested in [25] and implemented in the LAPACK subroutine
xGEBAL [3]. The eigenvalues of the balanced and the unbalanced Hamiltonian matrix
1 Matlab is a trademark of The MathWorks, Inc.
were computed by the square-reduced method using a Matlab function sqred which
implements the explicit version of the square-reduced method (see [31]).
We also tested the effects of symplectic balancing for the numerically backward
stable, structure-preserving method for the Hamiltonian eigenvalue problem presented
in [7]. Like the square-reduced method, this algorithm uses the square of the Hamiltonian
matrix. But it avoids forming the square explicitly using a symplectic URV-type
decomposition of the Hamiltonian matrix.
As reference values we used the eigenvalues computed by the unsymmetric QR
algorithm with Parlett/Reinsch balancing as implemented in the LAPACK expert
driver routine DGEEVX [3], applied to the Hamiltonian matrix and using quadruple
precision.
Moreover, we tested the effects of balancing when solving algebraic Riccati equations
with the structure-preserving multishift method presented in [1] for the examples
from the benchmark collection [6]. We only present some of the most intriguing results

Example 6.1. [6, Example 6] The system data come from an optimal control
problem for a J-100 jet engine as a special case of a multivariable servomechanism
problem. The resulting Hamiltonian matrix H 2 R 60\Theta60 has 8 isolated eigenvalues:
triple eigenvalues at \Sigma20:0 and simple eigenvalues at \Sigma33:3.
Algorithm 3.4 returns and for the permuted Hamiltonian
matrix we have
Hn+1:n+i l \Gamma1;n+1:n+i l
Next, the Hamiltonian submatrix
is scaled using Algorithm 4.4. After six sweeps, we obtain the balanced Hamiltonian
We have
decreased the 2-norm of the matrix used in the subsequent eigenvalue computation by
more than five orders of magnitude. If the eigenvalues are computed by the square-
reduced method applied to the unbalanced Hamiltonian matrix, the triplet of isolated
eigenvalues is returned as a pair of conjugate complex eigenvalues with relative errors
\Gamma11 and a simple eigenvalue with relative error  3:96 \Theta 10 \Gamma11 . For
the simple eigenvalue at 33:3, the relative error is  7:7 \Theta 10 \Gamma15 . For the balanced
version, these eigenvalues are returned with full accuracy since they are not affected
by roundoff errors. The relative errors for the other (not isolated) eigenvalues are
given in Figure 6.1 where we use the relative distance of the computed eigenvalues to
those computed by DGEEVX as an estimate of the real relative error.

Figure

6.1 only contains the relative errors for the eigenvalues with positive real
parts as sqred returns the eigenvalues as exact plus-minus pairs. The '+' for the 26th
eigenvalue is missing as the computed relative error for the balanced version is zero
with respect to machine precision. The eigenvalues are ordered by increasing absolute
values. From Figure 6.1, the increasing accuracy for decreasing ratio kHk 2 =jj is
obvious - with or without balancing. All computed eigenvalues of the balanced
matrix are more accurate than for the unbalanced one. The increase in accuracy is
more significant for the eigenvalues of smaller magnitude. This reflects the decrease
of the ratios kHk 2 =jj which more or less determines the accuracy of the computed
eigenvalues; see [31]. The decrease factor for kHk 2 is about 5 \Theta 10 \Gamma6 . The accuracy
for the eigenvalues of smaller magnitude increases by almost the same factor.
From

Figure

6.2 we see that symplectic balancing also improves the eigenvalues
computed by the method proposed in [7]. As the method does not suffer from the
perturbation, the accuracy for all computed eigenvalues is similar. Also note
that in the unbalanced version, the isolated eigenvalues are computed with a relative
accuracy ranging from 7:0 \Theta 10 \Gamma14 to 1:2 \Theta 10 \Gamma15 .
eigenvalue number
relative
errors
'+' - with symplectic balancing, 'o' - without balancing
Fig. 6.1. square-reduced method.
eigenvalue number
relative
errors
'+' - with symplectic balancing, 'o' - without balancing
Fig. 6.2. symplectic URV
Using the balanced matrix in order to solve algebraic Riccati equations by the
multishift method as described in [1], we obtain the following results: if the multishift
method is applied to the unbalanced data, the computed solution yields a residual
of size 1:5 \Theta 10 \Gamma6 while using the balanced Hamiltonian matrix we get r
This shows that numerical methods for solving algebraic Riccati equations can
substantially be improved employing balancing.
Example 6.2. [6, Example 13] The Hamiltonian matrix is defined as in (1.1)
with
diag (1; 0;
denotes the fourth unit vector. After four sweeps of Algorithm
4.4, kHk 2 is reduced from 10 12 to 1:5 \Theta 10 6 . The accuracy of the computed
eigenvalues did not improve significantly, but for the stabilizing solution of the algebraic
Riccati equation, the Frobenius norm of the residual as defined in (6.1) dropped
from r
7. Concluding Remarks. We have seen that isolated eigenvalues of a real
Hamiltonian matrix can be deflated using similarity transformations with symplectic
permutation matrices, the deflated problem can be scaled in order to reduce the norm
of the deflated Hamiltonian matrix and to equilibrate its row and column norms, and
the remaining (not isolated) eigenvalues can then be determined by computing the
eigenvalues of the deflated, balanced Hamiltonian submatrix. If invariant subspaces
are required, then we can use J-permutation matrices and a symplectic re-ordering
strategy in order to obtain the desired invariant subspaces. The same method can
be applied in order to balance skew-Hamiltonian and complex (skew-)Hamiltonian
matrices.
Numerical examples demonstrate that symplectic balancing can significantly improve
the accuracy of eigenvalues of Hamiltonian matrices as well as the accuracy of solutions
of the associated algebraic Riccati equations computed by structure-preserving
methods.
Final Remark and Acknowledgments. The work presented in this article
continues preliminary results derived in [4]. The author would like to thank Ralph
Byers, Heike Fabender, and Volker Mehrmann for helpful suggestions.



--R

A multishift algorithm for the numerical solution of algebraic Riccati equations


Contributions to the Numerical Solution of Algebraic Riccati Equations and Related Eigenvalue Problems

A collection of benchmark examples for the numerical solution of algebraic Riccati equations I: Continuous-time case
structure preserving method for computing the eigenvalues of real Hamiltonian or symplectic pencils

A bisection method for computing the H1 norm of a transfer matrix and related problems
A fast algorithm to compute the H1-norm of a transfer function matrix
Matrix factorization for symplectic QR-like methods


Matrix Eigensystem Routines- EISPACK Guide Extension
Matrix Computations
Matrix balancing
On the complexity of matrix balancing
The Algebraic Riccati Equation
Invariant subspace methods for the numerical solution of Riccati equations
Canonical forms for Hamiltonian and symplectic matrices and pencils
The Autonomous Linear Quadratic Control Problem
Solution of large matrix equations which occur in response theory

A Schur decomposition for Hamiltonian matrices
Balancing a matrix for calculation of eigenvalues and eigenvec- tors
Computational Methods for Linear Control Systems
Solving continuous-time matrix algebraic Riccati equations with condition and accuracy estimates
Algorithms for Linear-Quadratic Optimization
A fast algorithm to compute the real structured stability radius
Algorithm 506-HQR3 and EXCHNG: Fortran subroutines for calculating and ordering the eigenvalues of a real upper Hessenberg matrix
A symplectic method for approximating all the eigenvalues of a Hamiltonian matrix
--TR

--CTR
Peter Benner , Daniel Kressner , Volker Mehrmann, Structure preservation: a challenge in computational control, Future Generation Computer Systems, v.19 n.7, p.1243-1252, October
Pierluigi Amodio, On the computation of few eigenvalues of positive definite Hamiltonian matrices, Future Generation Computer Systems, v.22 n.4, p.403-411, March 2006
Peter Benner , Daniel Kressner, Algorithm 854: Fortran 77 subroutines for computing the eigenvalues of Hamiltonian matrices II, ACM Transactions on Mathematical Software (TOMS), v.32 n.2, p.352-373, June 2006
hamiltonian matrix;symplectic method;balancing;eigenvalues