Balancing Skew-Hamiltonian/Hamiltonian Pencils
With Applications in Control Engineering
Vasile Sima
National Institute for Research & Development in Informatics, 8–10 Bd. Mares¸al Averescu, Bucharest, Romania
Keywords:
Deflating Subspaces, Generalized Eigenvalues, Numerical Methods, Skew-Hamiltonian/Hamiltonian Matrix
Pencil, Software, Structure-preservation.
Abstract:
Badly-scaled matrix pencils could reduce the reliability and accuracy of computed results for many numeri-
cal problems, including computation of eigenvalues and deflating subspaces, which are needed in many key
procedures for optimal and H
control, model reduction, spectral factorization, and so on. Standard balancing
techniques can improve the results in many cases, but there are situations when the solution of the scaled
problem is much worse than that for the unscaled problem. This paper presents a new structure-preserving
balancing technique for skew-Hamiltonian/Hamiltonian matrix pencils, and illustrates its good performance
in solving eigenvalue problems and algebraic Riccati equations for large sets of examples from well-known
benchmark collections with difficult examples.
1 INTRODUCTION
Computing eigenvalues and bases of certain associ-
ated invariant or deflating subspaces of matrices or
matrix pencils is central to various numerical tech-
niques in control engineering and other domains.
When the corresponding eigenproblems have special
structure, implying structural properties of their spec-
tra, it is important to use structure-preserving and/or
structure-exploiting algorithms. General numerical
algorithms cannot ensure that the theoretical prop-
erties are preserved during computations (Paige and
Van Loan, 1981; Van Loan, 1984; Kressner, 2005).
Common special structures are Hamiltonian and sym-
plectic matrices or matrix pencils, which are encoun-
tered in optimal or H
control (e.g., solution of alge-
braic Riccati equations, or evaluation of the H
- and
L
-norms of linear dynamic systems), spectral factor-
ization, model reduction, etc., see, for instance, (Bru-
insma and Steinbuch, 1990; Laub, 1979; Mehrmann,
1991; Sima, 1996). While structured matrices are en-
countered for standard linear dynamic systems, gen-
eralized and descriptor systems may involve struc-
tured matrix pencils. Often such pencils can be re-
cast as skew-Hamiltonian/Hamiltonian pencils, for
which dedicated algorithms have been developed, see
e.g., (Benner et al., 2007; Benner et al., 2002; Ben-
ner et al., 2005; Benner et al., 2012a; Benner et al.,
2012b; Benner et al., 2013a; Benner et al., 2013b;
Jiang and Voigt, 2013), and the references therein.
Quite often, the pencil matrices have large norms
and elements with highly different magnitude. Such
pencils imply potential numerical difficulties for
software implementations of eigensolvers, with
negative consequences on the reliability and accuracy
of the results, see, e.g., (Sima and Benner, 2015a).
Balancing procedures can be used to improve the
numerical behavior. Ward (1981) proposed a bal-
ancing technique for general matrix pencils, which
has been incorporated in state-of-the-art software
packages, such as LAPACK (Anderson et al., 1999).
(This will be referred below as standard balancing.)
The data matrices are preprocessed by equivalence
transformations, in two optional stages: the first
stage uses permutations to find isolated eigenvalues
(which are available by inspection, with no rounding
errors), and the second stage uses diagonal scaling
transformations to make the row and corresponding
column 1-norms as close as possible. Balancing
may reduce the 1-norm of the scaled matrices, but
there is no guarantee in general. This is the reason
why full balancing is either avoided or provided as
an option in the LAPACK subroutines; some expert
driver routines, such as DGEES, DGEESX, DGGES,
DGGESX, and DGGEV use permutations only, while
other drivers, e.g., DGEEVX and DGGEVX, have an
input argument allowing either permutations, scaling,
both permutations and scaling, or no balancing at all.
Sima, V.
Balancing Skew-Hamiltonian/Hamiltonian Pencils - With Applications in Control Engineering.
DOI: 10.5220/0005981201770184
In Proceedings of the 13th International Conference on Informatics in Control, Automation and Robotics (ICINCO 2016) - Volume 1, pages 177-184
ISBN: 978-989-758-198-4
Copyright
c
2016 by SCITEPRESS Science and Technology Publications, Lda. All rights reserved
177
A structure-preserving balancing technique for
Hamiltonian matrices has been developed in (Benner,
2001), but no special procedure has been available
for skew-Hamiltonian/Hamiltonian matrix pencils.
This paper presents a new, structure-preserving
balancing procedure for skew-Hamiltonian/Hamil-
tonian pencils, and illustrates its performance for
some control engineering applications, in comparison
with the general approach, and with the case when
no balancing is used. This procedure includes several
optional enhancements which allow to get meaning-
ful results when standard balancing (possibly even a
structured variant) fails, as shown for some numerical
examples.
2 BALANCING
SKEW-HAMILTONIAN/
HAMILTONIAN PENCILS
Let αS βH be a complex skew-Hamiltonian/
Hamiltonian pencil, with α, β C, S C
2n×2n
a
skew-Hamiltonian matrix and H C
2n×2n
a Hamil-
tonian matrix, i.e.,
(SJ )
H
= SJ , (H J )
H
= H J , J :=
0 I
n
I
n
0
,
(1)
where M
H
denotes the conjugate transpose of a ma-
trix M, M
H
=
¯
M
T
, the overbar denotes the complex
conjugate, M
T
denotes the transpose of M, and I
n
is
the identity matrix of order n. These definitions imply
the following structure for S and H :
S =
A D
E A
H
, H =
C V
W C
H
, (2)
where D and E are skew-Hermitian, i.e., D
H
= D,
E
H
= E, and V and W are Hermitian matrices.
When S and H are real matrices, then D and E are
skew-symmetric, V and W are symmetric, the H su-
perscript is replaced by T, and some derivations be-
low simplify.
When the matrices S and H are badly scaled, e.g.,
when their elements have moduli with highly differ-
ent magnitude, the accuracy of computed solutions of
related eigenproblems can be very poor. The accuracy
may be improved by using a balancing procedure, as
described in (Ward, 1981) for general matrix pencils,
and implemented in the LAPACK package (Anderson
et al., 1999). Similarly to the techniques in (Anderson
et al., 1999; Benner, 2001; Ward, 1981), a structure-
exploiting balancing procedure for the structured pen-
cil matrices in (2) would involve, in a first step,
permuting αS βH by a symplectic equivalence
transformation to isolate eigenvalues in the elements
1 : 1 and n+1 : n+1 on the diagonals of S and
H , and in a second step, applying a diagonal equiva-
lence transformation to rows and columns : n and
n + : 2n, to make the rows and columns as close
in 1-norm as possible. (A MATLAB-style notation
for array indexing (MATLAB, 2016) is used.) Both
steps above are optional. Balancing may reduce the
1-norms of the matrices S and H .
Note that it is enough to equilibrate the 1-norms
of the rows and columns 1:n of the matrices S and H ,
since kS
i+n,:
k
p
= kS
:,i
k
p
and kS
:,i+n
k
p
= kS
i,:
k
p
, for
any p-norm, for p 1, and similarly for H (see (Ben-
ner, 2001)). The balancing procedure performs the
following transformation:
e
S = LS R ,
e
H = LH R , (3)
where
L = L
n
···L
P
T
1
···P
T
1
, R = P
1
···P
1
R
···R
n
,
(4)
where P
k
, k = 1 : 1, are symplectic permutation or
J -permutation matrices (Benner, 2001), and L
k
and
R
k
are left and right diagonal scaling matrices (with
L
k
(k,k), L
k
(k+n,k+n), R
k
(k,k), and R
k
(k+n,k+n)
the only diagonal elements different from 1), chosen
to equilibrate the 1-norms of the k-th row and col-
umn of S and H , for k = ℓ, · · · ,n. A permutation
matrix has exactly one nonzero entry, which is 1, in
each row and column. A J -permutation matrix has a
similar structure, but the nonzero entries may also be
1. When possible, the permutations are chosen in
the form
P =
P 0
0 P
, (5)
where P is an n×n permutation matrix (hence, P
T
P =
PP
T
= I
n
). Clearly, P is symplectic, since P J P
T
=
J . As shown in (Benner, 2001), it is generally not
possible to preserve the (skew-)Hamiltonian structure
using symplectic permutations only, but symplectic
J -permutations are also needed.
Applying a permutation (5) in (3) is easy. Specifi-
cally, the transformedmatrices become P
T
AP, P
T
DP,
P
T
EP, P
T
CP, P
T
VP, and P
T
WP.
Consider applying a J -permutation, P, which has
the following structure, defined by an index i, i n,
P =
I
i1
0 0 0 0 0
0 0 0 0 1 0
0 0 I
ni
0 0 0
0 0 0 I
i1
0 0
0 1 0 0 0 0
0 0 0 0 0 I
ni
, (6)
to the skew-Hamiltonian matrix S, and Hamiltonian
ICINCO 2016 - 13th International Conference on Informatics in Control, Automation and Robotics
178
matrix H . Using the partition of H in (2), the trans-
formed matrix P
T
H P is
P
T
H P =
v
1:i1,i
w
i,1:i1
¯c
ii
w
H
i+1:n,i
v
H
i,i+1:n
c
H
i,1:i1
c
i,1:i1
v
ii
c
i,i+1:n
c
H
i,i+1:n
c
1:i1,i
c
H
1:i1,i
w
ii
c
H
i+1:n,i
c
i+1:n,i
w
H
i,1:i1
v
H
1:i1,i
c
ii
v
i,i+1:n
w
i+1:n,i
, (7)
where denotes submatrices which do not change.
Matrix S is transformed similarly.
The permutation procedure first scans the columns
1 : n and then the rows 1 : n of S and H , and uses
row and column permutations with matrices of the
form (5) and (6) so that
e
S and
e
H finally have the
first 1 columns upper triangular, i.e., with zeros
below the diagonals. This way, the first 1 eigen-
values of the pencil are defined by the diagonal ele-
ments 1 :1 in
e
S and
e
H . The upper triangular struc-
ture is built column by column, starting from the first.
Therefore, the permutationsabove can take advantage
of the zero subdiagonal part in the already processed
columns. For instance, when a J -permutation (6) is
formed and applied, the vectors c
i,1:i1
, and w
i,1:i1
in (7) (and similarly for a
i,1:i1
, e
i,1:i1
in S), as well
as all vectors under the diagonals in the first i 1
columns are already zero.
The balancing procedure for (skew-)Hamiltonian
matrices must use symplectic scaling matrices as well,
in order to preserve the structure (Benner, 2001). For-
tunately, this is not needed for skew-Hamiltonian/Ha-
miltonian pencils, if diagonal scaling is used, as men-
tioned above. Indeed, without loss of generality, con-
sider that permutations are not wanted, and only scal-
ing should be applied. Let L and R be the real diago-
nal matrices which scale the first n rows and the first
n columns, respectively, of S and H . Since R and L
will be the scaling matrices for the last n rows and
columns, respectively, the global scaling transforma-
tion is
e
S =
L 0
0 R
A D
E A
H
R 0
0 L
,
e
H =
L 0
0 R
C V
W C
H
R 0
0 L
. (8)
Then,
e
SJ =
LAR LDL
RER RA
H
L
0 I
n
I
n
0
=
LDL LAR
RA
H
L RER
= (
e
SJ )
H
, (9)
since D = D
H
, E = E
H
. Similarly,
e
H J =
LCR LVL
RWR RC
H
L
0 I
n
I
n
0
=
LVL LCR
RC
H
L RWR
= (
e
H J )
H
, (10)
since V = V
H
and W = W
H
. Therefore, the struc-
ture of S and H is preserved for the transformation
in (8). Note that the left and right transformations
are not symplectic if L 6= R
1
, but they are structure-
preserving, for diagonal L and R. The permutations
used for isolating the eigenvalues must, however, be
symplectic.
As described above, the elementary scaling matri-
ces in (4) are, for k = ℓ, · ·· ,n,
L
k
= bl
diag(I
k1
,l
k
,I
n1
,r
k
,I
nk
),
R
k
= bl diag(I
k1
,r
k
,I
n1
,l
k
,I
nk
), (11)
where bl
diag(X,Y,...) denotes a bloc-diagonal ma-
trix with diagonal blocks X, Y, etc. Note that
L
1
= P
1
···P
1
L
1
···L
1
n
. (12)
Formula (12) can be used to apply, from the left, the
back balancing transformations to a given 2n×m ma-
trix. Since L
1
has a formula similar to R , to obtain
R , one can apply to I
2n
the transformations in (12),
with scaling factors l
k
and r
k
replaced by 1/r
k
and
1/l
k
, respectively. There is no need for an algorithm
to apply the right back transformation.
A special case of application of the balancing
transformations appears when computing the solu-
tion of algebraic Riccati equations (AREs) using the
skew-Hamiltonian/Hamiltonian approach, see (Ben-
ner et al., 2013a), and the references therein. This
approach uses a basis for the stable right deflating
subspace of the matrix pencil αS βH , where S
and H are suitably defined in terms of the dynamic
system and performance index matrices. For conve-
nience, assume that S = I
2n
, which is the case for stan-
dard continuous-time dynamic systems, e.g., when
the control weighting matrix of the performance index
is well-conditioned. (The general case will be briefly
described below.) Let
e
U =
e
U
T
1
e
U
T
2
T
be a ba-
sis of the stable right deflating subspace of α
e
S β
e
H ,
corresponding to the balanced pencil. The stabilizing
solution of the balanced ARE is given by
e
X =
e
U
2
e
U
1
1
.
Balancing Skew-Hamiltonian/Hamiltonian Pencils - With Applications in Control Engineering
179
Note that
e
U is related to a basis, U, of the stable right
deflating subspace of the original pencil, αS βH , by
the transformation
e
U = bl
diag(R
1
,L
1
)U, hence
U
T
1
U
T
2
T
:= U = bl
diag(R,L)
e
U. Therefore,
the stabilizing solution of the original ARE can be
computed as follows
X = U
2
U
1
1
= L
e
U
2
e
U
1
1
R
1
. (13)
Formula (13) allows to represent and use the solution
X in a factored form, which may be useful for nu-
merical reasons. If balancing also involves symplec-
tic permutations of the form (5) only, then the right
hand side in (13) becomes PL
e
U
2
e
U
1
1
R
1
P
T
, where
P denotes here the product of the 1 permutations
performed. Such factorization, is, however, not possi-
ble if J -permutations are also needed, in which case
it is necessary to apply the balancing transformations
before solving the set of linear systems giving
e
X.
If S is a general matrix, the order of S and H is
2(n+ p), where n is the dimension of the state vector,
and p may be chosen as p = m/2 (i.e., p = m/2, if
m is even, and p = (m+ 1)/2, otherwise), with m the
number of system inputs. In this case, the computa-
tions are similar, but
e
U
i
and U
i
, i = 1,2, refer to the
lines 1 : n and n + p + 1 : 2n+ p of the matrix bases
e
U and U, respectively.
3 IMPLEMENTATION ISSUES
The developed balancing algorithm operates only on
the matrices A, D, E, C, V, and W, and preserves
the pencil structure. Moreover, the pairs of skew-
Hermitian matrices, D and E, and Hermitian matri-
ces, V andW, are stored compactly in two n× (n+1)
arrays,
DE
and
VW
. Specifically, the lower triangle of
E and the upper triangle of D are concatenated along
their main diagonals, and similarly for W and V:
DE
ij
= e
ij
, j = 1 : n, i j,
DE
i, j+1
= d
ij
, j = 1 : n, i j,
VW
ij
= w
ij
, j = 1 : n, i j,
VW
i, j+1
= v
ij
, j = 1 : n, i j. (14)
This way, the storage requirements are reduced by
4n
2
2n memory locations, in comparison with the
general algorithm in (Anderson et al., 1999), which
needs 8n
2
storage space for S and H . Note that in
the real case, the diagonal elements of D and E, d
ii
and e
ii
, i = 1 : n, which should be zero, by definition,
are not stored and not used. However, in the complex
case, d
ii
and e
ii
, i = 1 : n, should have zero real parts,
while the imaginary parts might be nonzero; more-
over, the diagonal elements of V and W should have
zero imaginary parts, but possibly nonzero real parts.
The implementation takes care of the compact storage
scheme when computing P
T
DP, ..., P
T
WP.
Consider a structured pencil, (S ,H ), and assume
that S = I
2n
. The order of magnitude of the dif-
ferences in size of the nonzero elements in the H
matrix can be huge. For instance, for the CM6
and CM6
IS examples from the COMPl
e
ib collec-
tion (Leibfritz and Lipinski, 2003), the maximum
and minimum absolute values of the nonzero en-
tries are about 2.53· 10
5
and 4.9407 · 10
324
, hence
their ratio is not representable in the usual double
precision representation, and it is evaluated as Inf
(). The standard balancing algorithm implemented
in the LAPACK Library subroutine
DGGBAL
(for ma-
trix pencils) returns scaling factors covering a very
large range of magnitudes, namely with maximum
and minimum values 10
159
and 10
47
, respectively,
for both left and right scaling factors. No eigenvalue
can be isolated. The scaling transformation matrices,
L and R , have condition numbers with values 10
206
.
The 1-norms of scaled matrices,
e
H and
e
S, are about
10
228
and 10
184
, respectively, while the 1-norms of
the original H and S are 6.9123· 10
5
and 1. Com-
puting the eigenvalues or deflating subspaces using
the balanced matrices would return results very far
from the true ones. CM6 and CM6
IS are not the only
examples in the COMPl
e
ib collection which produce
numerical troubles for eigenproblem-related compu-
tations. The CM5 and CM5
IS examples have also
a huge ratio, 10
279
, between the magnitudes of their
elements. Other very large ratios are 10
141
, for CM4
and CM4
IS, 10
72
, for CM3 and CM3 IS, or 10
37
, for
CM2 and CM2 IS.
The proposed algorithm uses an adaptation of
the basic LAPACK procedure for finding the scaling
factors, but optionally limits the range of their
variation, possibly via an outer loop. Specifically,
the user can set a threshold value, τ; if τ 0,
the entries whose absolute values are smaller than
τM
0
, where M
0
= max(kH (s,s)k
1
,kS(s,s)k
1
), with
s := : n n + : 2n, are considered negligible,
and do not count for computing the scaling factors.
(For the CM6 and CM6 IS examples, = 1, and
setting τ = 10
20
, all entries with magnitude less
than about 7 · 10
15
will be taken as zero by the
procedure.) If τ < 0 on entry, an outer loop over a set
of values τ
i
> 0 will enable to select a set of scaling
factors which, if possible, ensure the reduction of a
desired norm-related measure for the scaled matri-
ces. For τ = 1, this measure is the minimum of
max
i
(kH
i
(s,s)k
1
/kS
i
(s,s)k
1
,kS
i
(s,s)k
1
/kH
i
(s,s)k
1
)
where H
i
(s,s) and S
i
(s,s) are the scaled submatrices
corresponding to the threshold τ
i
. This strategy
ICINCO 2016 - 13th International Conference on Informatics in Control, Automation and Robotics
180
tries to balance H and S, but also have comparable
1-norms. For τ = 2, the same measure is used, but if
max(k
e
H (s,s)k
1
,k
e
S(s,s)k
1
) > cM
0
and t > T, where
c and T are given constants (c possibly larger than
1), and t is the maximum ratio of the scaling factors
found (the maximum of the condition numbers of
L and R ), then the scaling factors are set to 1 and
a warning indicator is set; here, the matrices with
tilde accent are the solution of the above norm ratio
reduction problem. This approach avoids to obtain
scaled matrices with too large norms, compared to
the given ones, and also limits the range of the scaling
factors. For τ = 3, the measure used is the smallest
product of norms, min
i
(kH
i
(s,s)k
1
kS
i
(s,s)k
1
), over
the sequence of τ
i
values tried, while for τ = 4,
the condition numbers of the scaling transformations
are additionally supervised, and the scaling factors
are set to 1 if the “optimal” scaling has a condition
number larger than T. This tends to reduce the
1-norms of both matrices. Finally, if τ = 10
k
, the
condition numbers of the acceptable scaling matrices
are bounded by 10
k
.
4 NUMERICAL RESULTS
An extensive testing has been performed to assess the
performance of the new balancing solver for skew-
Hamiltonian/Hamiltonian matrix pencils. Computa-
tion of eigenvalues, as well as solution of AREs, have
been considered as main applications. This section
summarizes part of the results.
The COMPl
e
ib collection (Leibfritz and Lipin-
ski, 2003) is taken here as a benchmark for eigen-
value computations. All 151 problems with n < 2000
have been tried. For testing purposes, the balancing
algorithms have been used in combination with the
MATLAB eigensolvers
eig
(H ), when S = I
2n
, and
eig
(H ,S), when S is general, or, for small-size prob-
lems, with the symbolic solvers
eig(vpa
(H )) (with
default number of digits, i.e., 32), or
eig(sym
(H )).
Usually,
eig(vpa
(H )) and
eig(sym
(H )) have been
used for examples with n 200 and n 30, respec-
tively. Examples with larger orders need significant
CPU times for symbolic solvers, e.g., over one hour
for CM4 or CM4
IS (n = 240) using
vpa
, or much
longer for
sym
. But a comparison between the eigen-
values computed by
eig(vpa
(H )) and
eig(sym
(H ))
has shown a very good agreement between their re-
sults. The maximum absolute difference between
the relative errors of
eig
(H ) in comparison with
eig(vpa
(H )) and
eig(sym
(H )) for 100 COMPl
e
ib
examples with n 30 is 3.0292· 10
28
. (Examples
JE2 and JE3, with n = 21 and n = 24, respectively,
have been excluded, since
eig(sym
(H )) needed un-
reasonably large CPU time.) However, the rela-
tive errors of
eig
(H ,S), when S = I
2n
, compared
to
eig(vpa
(H )) and
eig(sym
(H )), for example PAS
have been 1.2798· 10
6
, and 1.4432 · 10
6
, respec-
tively. Omitting PAS, the norm of the differences
between relative errors of
eig
(H ,S) compared to
eig(vpa
(H )) and
eig(sym
(H )) was 5.8491· 10
25
.
There are 18 examples with ratios between the
magnitude of the largest and smallest moduli of
their elements larger than 10
16
. Besides exam-
ples CM2* CM6*, this set also includes AC10,
BDT2, PAS, TL, CDP, NN16, ISS1, and ISS2 ex-
amples. Most of these are difficult examples for
any balancing algorithm, and for generalized eigen-
solvers using standard balancing. For instance, for
CDP example, even the structure-exploiting skew-
Hamiltonian/Hamiltonian solver with standard bal-
ancing, for S = I
2n
, returns eigenvalues with a relative
error of 0.19172 (when compared to
eig
(H )), and
of 0.24788 (when compared to
eig(vpa
(H )), while
the same solver without scaling delivers relative er-
rors 1.0476 · 10
14
and 1.3064 · 10
14
, respectively.
These values are close to the error obtained using
eig
(H ), 4.9248· 10
15
, compared to
eig(vpa
(H )).
The structure-exploiting solver with balancing option
τ < 0 returns a relative error of 5.4838 · 10
15
(when
compared to
eig
(H )), and of 3.9091· 10
15
(when
compared to
eig(vpa
(H )), even more accurate than
eig
(H ).
Figure 1 shows the relative errors of eigenval-
ues computed using
eig
(Hl,Sl), where Hl :=
e
H , Sl
:=
e
S, with LAPACK balancing (using DGGBAL)
and
eig
(H ,S ), in comparison with
eig
(H ) for 151
COMPl
e
ib examples. The error is infinite for CDP,
CM3, CM4, CM3
IS, and CM4 IS examples. There
are also several large errors using LAPACK balanc-
ing. However, for many examples, balancing im-
proves the accuracy of the computed eigenvalues.
Similar results have been obtained using the struc-
tured balancing solver with option τ = 0. Specifically,
except for ve examples, the differences between the
relative errors of eigenvalues computed using LA-
PACK balancing and the new balancing solver with
option τ = 0 have been less than about 3.64· 10
3
.
Figure 2 shows the relative errors of eigenvalues
computed using
eig
(Hl,Sl) with the new balancing
solver with option τ = 1 and
eig
(H ,S), in compari-
son with
eig
(H ) for 151 COMPl
e
ib examples. There
are no infinite errors. The errors are smaller for many
examples than when using LAPACK balancing. Note
that the ordinate axes have different scales. Note also
that
eig
(H ,S ) has often larger errors than
eig
(Hl,Sl)
with the new solver. The behaviorfor τ= 3 has been
Balancing Skew-Hamiltonian/Hamiltonian Pencils - With Applications in Control Engineering
181
0 20 40 60 80 100 120 140 160
Example #
10
-20
10
-15
10
-10
10
-5
10
0
10
5
Relative errors
eig(H,S) and eig(Hl,Sl), after LAPACK balancing compared to eig(H)
eig(H,S)
eig(Hl,Sl)
Figure 1: Relative errors of eigenvalues computed using
eig
(Hl,Sl) with LAPACK balancing and
eig
(H,S), com-
pared to
eig
(H), for COMPl
e
ib examples.
0 20 40 60 80 100 120 140 160
Example #
10
-20
10
-15
10
-10
10
-5
10
0
Relative errors
Structured solver (τ = -1) and eig(H,S) compared to eig(H)
eig(H,S)
new solver
Figure 2: Relative errors of eigenvalues computed us-
ing
eig
(Hl,Sl) with new balancing solver, τ = 1, and
eig
(H,S), compared to
eig
(H), for COMPl
e
ib examples.
almost the same (the maximum differencebetween er-
rors for the two options has been of order 10
9
).
Figure 3 illustrates the relative errors of eigen-
values of αS + βH , computed using
eig
and
the SLICOT structured eigensolver MB04BD (see
www.slicot.org), without balancing, in comparison
with
eig
(H ). The structured eigensolver is more ac-
curate than
eig
for almost all examples and compa-
rable with
eig
for the remaining examples. Figure 4
uses similarly the structured eigensolver applied to
e
S
and
e
H , obtained with balancing option τ = 1. It can
provide even more accurate results than
eig
(H ,S ),
see Fig. 2. For comparison, Fig. 5 shows results for
option τ = 0. Clearly, standard balancing does not
provide good results in all cases, even using a struc-
tured solver.
The remaining of this section considers the per-
formance of the new balancing solver, in combina-
tion with the structured SLICOT routine MB03LD,
to compute the deflating subspaces of skew-
Hamiltonian/Hamiltonian pencils, for solving AREs
0 20 40 60 80 100 120 140 160
Example #
10
-20
10
-15
10
-10
10
-5
10
0
Relative errors
Structured eigensolver on (H,S), no balancing, and eig(H,S) compared to eig(H)
eig(H,S)
structured eigensolver
Figure 3: Relative errors of eigenvalues computed by
SLICOT structured eigensolver without balancing and
eig
(H,S), compared to
eig
(H), for COMPl
e
ib examples.
0 20 40 60 80 100 120 140 160
Example #
10
-20
10
-15
10
-10
10
-5
10
0
Relative errors
Structured eigensolver (τ = -1) and eig(H,S) compared to eig(H)
eig(H,S)
structured eigensolver
Figure 4: Relative errors of eigenvalues computed by
SLICOTstructured eigensolver with balancing, τ = 1, and
eig
(H,S), compared to
eig
(H), for COMPl
e
ib examples.
0 20 40 60 80 100 120 140 160
Example #
10
-20
10
-15
10
-10
10
-5
10
0
10
5
Relative errors
Structured eigensolver (τ = 0) and eig(H,S) compared to eig(H)
eig(H,S)
structured eigensolver
Figure 5: Relative errors of eigenvalues computed by
SLICOT structured eigensolver with balancing, τ = 0, and
eig
(H,S), compared to
eig
(H), for COMPl
e
ib examples.
from the SLICOT CAREX collection (Abels and
Benner, 1999). This combination is referred to as
skHH
in the legends of the following figures. Both
pencils of order 2n (with S = I
2n
) and extended pen-
cils of order 2(n + p) have been tried. The com-
ICINCO 2016 - 13th International Conference on Informatics in Control, Automation and Robotics
182
0 2 4 6 8 10 12 14 16 18
Example #
10
-20
10
-15
10
-10
10
-5
10
0
Relative errors
Relative errors to care solution using balancing
skHH
Figure 6: Relative errors of the stabilizing ARE solutions
computed by SLICOT structured eigensolver (for extended
pencil), with balancing option τ = 0, compared to
care
, for
examples from the SLICOT CAREX collection.
0 2 4 6 8 10 12 14 16 18
Example #
10
-16
10
-14
10
-12
10
-10
10
-8
10
-6
10
-4
10
-2
Relative errors
Relative errors to care solution using balancing
skHH
Figure 7: Relative errors of the stabilizing ARE solutions
computed by SLICOT structured eigensolver (for 2n-order
pencil), with option τ = 1, compared to
care
, for exam-
ples from the SLICOT CAREX collection.
0 5 10 15
Example #
10
-30
10
-25
10
-20
10
-15
10
-10
10
-5
10
0
Relative errors
Relative errors using balancing
skHH
care
Figure 8: Relative errors of the stabilizing ARE solutions
computed by
care
and SLICOT structured eigensolver (for
extended pencil), with option τ = 0, compared to exact so-
lution, for examples from the SLICOT CAREX collection.
puted solutions are compared with the exact ones,
when known, or with the solutions returned by the
MATLAB function
care
. Note that
care
uses a
special balancing procedure, and Hamiltonian ma-
trices, instead of pencils, for all examples, except
for CAREX example 2.2 with default parameter, for
which the control weighting matrix tends to be sin-
gular. This creates an advantage for
care
over the
new solver. However, both solvers produce compa-
rable results, when the balancing option τ = 1 is
used. See (Sima and Benner, 2015b) for a compari-
son between a Hamiltonian-based solver and
care
for
CAREX examples.
Figure 6 shows the relative errors of the stabiliz-
ing solution of the AREs computed using SLICOT
structured eigensolver (for extended pencil), with bal-
ancing option τ = 0, in comparison with
care
, while
Fig. 7 plots similarly the relative errors when option
τ = 1 is used. Moreover,Fig. 8 and Fig. 9 present in
the same manner the relative errors compared to the
known, exact solutions. The balancing option τ = 1
ensures better results than the standard balancing (for
τ = 0). Finally, Fig. 10 plots the relative residuals for
care
and the structured solver.
0 5 10 15
Example #
10
-30
10
-25
10
-20
10
-15
10
-10
10
-5
10
0
Relative errors
Relative errors using balancing
skHH
care
Figure 9: Relative errors of the stabilizing ARE solutions
computed by
care
and SLICOT structured eigensolver (for
2n-order pencil), with option τ = 1, compared to exact so-
lution, for examples from the SLICOT CAREX collection.
0 5 10 15 20 25 30 35
Example #
10
-25
10
-20
10
-15
10
-10
10
-5
10
0
Relative errors
Relative residuals using balancing
skHH
care
Figure 10: Relative residuals of the stabilizing ARE so-
lutions computed by
care
and SLICOT structured eigen-
solver (for 2n-order pencil), with option τ = 1, for exam-
ples from the SLICOT CAREX collection.
Balancing Skew-Hamiltonian/Hamiltonian Pencils - With Applications in Control Engineering
183
5 CONCLUSIONS
A new structure-preserving balancing technique for
skew-Hamiltonian/Hamiltonian matrix pencils is pre-
sented. Symplectic (J -)permutations and equivalence
scaling transformations are used. Several enhance-
ments are described, which avoid a large increase of
the norms of the pencil matrices, and/or of the con-
dition numbers of the scaling transformations, which
can appear when using the standard balancing pro-
cedure. The numerical results show a good perfor-
mance of the new technique in comparison with state-
of-the-art solvers. Tens of examples from well-known
benchmark collections have been investigated.
ACKNOWLEDGEMENTS
The NICONET support is highly acknowledged.
REFERENCES
Abels, J. and Benner, P. (1999). CAREX—A collection
of benchmark examples for continuous-time algebraic
Riccati equations (Version 2.0). SLICOT Working
Note 1999-14.
Anderson, E., Bai, Z., Bischof, C., Blackford, S., Dem-
mel, J., Dongarra, J., Du Croz, J., Greenbaum, A.,
Hammarling, S., McKenney, A., and Sorensen, D.
(1999). LAPACK Users’ Guide: Third Edition. SIAM,
Philadelphia.
Benner, P. (2001). Symplectic balancing of Hamiltonian
matrices. SIAM J. Sci. Comput., 22(5):1885–1904.
Benner, P., Byers, R., Losse, P., Mehrmann, V., and
Xu, H. (2007). Numerical solution of real skew-
Hamiltonian/Hamiltonian eigenproblems. Technical
report, Technische Universit¨at Chemnitz, Chemnitz.
Benner, P., Byers, R., Mehrmann, V., and Xu, H. (2002).
Numerical computation of deflating subspaces of
skew Hamiltonian/Hamiltonian pencils. SIAM J. Ma-
trix Anal. Appl., 24(1):165–190.
Benner, P., Kressner, D., and Mehrmann, V. (2005). Skew-
Hamiltonian and Hamiltonian eigenvalue problems:
Theory, algorithms and applications. In Proceedings
of the Conference on Applied Mathematics and Scien-
tific Computing, 3–39. Springer-Verlag, Dordrecht.
Benner, P., Sima, V., and Voigt, M. (2012a). L
-norm com-
putation for continuous-time descriptor systems us-
ing structured matrix pencils. IEEE Trans. Automat.
Contr., AC-57(1):233–238.
Benner, P., Sima, V., and Voigt, M. (2012b). Robust and
efficient algorithms for L
-norm computations for de-
scriptor systems. In 7th IFAC Symposium on Robust
Control Design (ROCOND’12), 189–194. IFAC.
Benner, P., Sima, V., and Voigt, M. (2013a). FOR-
TRAN 77 subroutines for the solution of skew-
Hamiltonian/Hamiltonian eigenproblems. Part I: Al-
gorithms and applications. SLICOT Working Note
2013-1.
Benner, P., Sima, V., and Voigt, M. (2013b). FOR-
TRAN 77 subroutines for the solution of skew-
Hamiltonian/Hamiltonian eigenproblems. Part II: Im-
plementation and numerical results. SLICOT Working
Note 2013-2.
Bruinsma, N. A. and Steinbuch, M. (1990). A fast algo-
rithm to compute the H
-norm of a transfer function.
Systems Control Lett., 14(4):287–293.
Jiang, P. and Voigt, M. (2013). MB04BV A FOR-
TRAN 77 subroutine to compute the eigenvectors as-
sociated to the purely imaginary eigenvalues of skew-
Hamiltonian/Hamiltonian matrix pencils. SLICOT
Working Note 2013-3.
Kressner, D. (2005). Numerical Methods for General and
Structured Eigenvalue Problems. Springer-Verlag,
Berlin.
Laub, A. J. (1979). A Schur method for solving algebraic
Riccati equations. IEEE Trans. Automat. Contr., AC
24(6):913–921.
Leibfritz, F. and Lipinski, W. (2003). Description of the
benchmark examples in COMPl
e
ib. Technical report,
Department of Mathematics, University of Trier, D–
54286 Trier, Germany.
MATLAB (2016). MATLAB
R
Primer. R2016a. The Math-
Works, Inc., 3 Apple Hill Drive, Natick, MA.
Mehrmann, V. (1991). The Autonomous Linear Quadratic
Control Problem. Theory and Numerical Solution.
Springer-Verlag, Berlin.
Paige, C. and Van Loan, C. F. (1981). A Schur decomposi-
tion for Hamiltonian matrices. Lin. Alg. Appl., 41:11–
32.
Sima, V. (1996). Algorithms for Linear-Quadratic Opti-
mization. Marcel Dekker, Inc., New York.
Sima, V. and Benner, P. (2015a). Pitfalls when solving
eigenproblems with applications in control engineer-
ing. In Proceedings of the 12th International Con-
ference on Informatics in Control, Automation and
Robotics (ICINCO-2015), 21-23 July, 2015, Colmar,
France, volume 1, 171–178. SciTePress.
Sima, V. and Benner, P. (2015b). Solving SLICOT bench-
marks for continuous-time algebraic Riccati equa-
tions by Hamiltonian solvers. In Proceedings of the
2015 19th International Conference on System The-
ory, Control and Computing (ICSTCC 2015), October
14-16, 2015, Cheile Gradistei, Romania, 1–6. IEEE.
Van Loan, C. F. (1984). A symplectic method for approx-
imating all the eigenvalues of a Hamiltonian matrix.
Lin. Alg. Appl., 61:233–251.
Ward, R. C. (1981). Balancing the generalized eigenvalue
problem. SIAM J. Sci. Stat. Comput., 2:141–152.
ICINCO 2016 - 13th International Conference on Informatics in Control, Automation and Robotics
184