STRUCTURE-PRESERVING ALGORITHMS FOR DISCRETE-TIME
ALGEBRAIC MATRIX RICCATI EQUATIONS
Vasile Sima
National Institute for Research & Development in Informatics, 8-10 Bd. Mares¸al Averescu, Bucharest, Romania
Keywords:
Linear-quadratic Optimization, Optimization Algorithms, Structure-preserving Algorithms.
Abstract:
Structure-preserving algorithms for solving discrete-time algebraic matrix Riccati equations are presented.
The proposed techniques extract the stable deflating subspaces for extended, inverse-free symplectic matrix
pencils. The algorithms are based on skew-Hamiltonian/Hamiltonian pencils derived by an extended Cayley
transformation, which only involves matrix additions and subtractions. The structure-preserving approach
has the potential to avoid the numerical difficulties which are encountered for a traditional, non-structured
solution, returned by the currently available software tools.
1 INTRODUCTION
Consider the continuous-time algebraic Riccati equa-
tion (CARE),
0 = Q+ A
H
XE + E
H
XA− E
H
XWXE
and the discrete-time algebraic Riccati equation
(DARE),
E
H
XE = Q + A
H
XA− A
H
XB(R+ B
H
XB)
−1
B
H
XA,
where A,E,W,Q ∈ C
n×n
, B ∈ C
n×m
, R ∈ C
m×m
,
W = W
H
, Q = Q
H
, R = R
H
, E is nonsingular, and
W := BR
−1
B
H
. More general equations are obtained
by replacing E
H
XB and A
H
XB above by L + E
H
XB
and L+ A
H
XB, respectively, where L ∈ C
n×m
. The
real case is obtained by replacing C by R, and the
conjugate-transpose operator H by the transpose op-
erator T.
In applications, usually the stabilizing solution X
∗
is required, hence, e.g., for DARE, λE − (A − B(R +
B
H
X
∗
B)
−1
B
H
X
∗
A) is a (Schur) stable matrix pencil,
i.e., Λ(A−B(R+ B
H
X
∗
B)
−1
B
H
X
∗
A,E) ⊂ C
−
:= {z ∈
C : |z| < 1}, where Λ(M) denotes the spectrum of a
matrix or pencil M.
CAREs and DAREs arise in many applications,
such as, stabilization and linear-quadratic regulator
problems, Kalman filtering, linear-quadratic Gaus-
sian (H
2
-) optimal control problems, computation of
(sub)optimal H
∞
controllers, model reduction tech-
niques based on stochastic, positive or bounded real
LQG balancing, factorization procedures for transfer
functions (here, usually E 6= I
n
).
There are several basic approaches for solving al-
gebraic Riccati equations (AREs):
1. Treat an ARE as a nonlinear system of equations
using Newton’s method (with line search).
2. Use the connection to Hamiltonian eigenproblem.
The second approach for CARE, with E = I
n
, is based
on the identity

A −W
−Q −A
H

I
n
X

=

I
n
X

(A−WX),
hence, if X is stabilizing, then Λ(A−WX) = Λ(H) ∩
C
−
, where H is the first matrix in the above formula,
and C
−
:= {z ∈ C : ℜ(z) < 0}. Consequently, the
columns of [ I
n
X
T
]
T
span the stable invariant sub-
space of the Hamiltonian matrix H. Therefore, it is
possible to compute the stable H-invariant subspace
via eigendecomposition or block-Schur factorization,
U
−1
HU =
"
H
11
H
12
0 H
22
#
, U =
"
U
11
U
12
U
21
U
22
#
,
and the solution is given by X = U
21
U
−1
11
.
If R is ill-conditioned, it is advisable to use ex-
tended matrix pencils, for better accuracy (Bender
and Laub, 1987a; Bender and Laub, 1987b; Lancaster
and Rodman, 1995; Mehrmann, 1991; Van Dooren,
1981):
– extended pencil for CARE:
N − λM =


A 0 B
Q A
H
L
L
H
B
H
R


− λ


E 0 0
0 −E
H
0
0 0 0


;
187
Sima V. (2010).
STRUCTURE-PRESERVING ALGORITHMS FOR DISCRETE-TIME ALGEBRAIC MATRIX RICCATI EQUATIONS.
In Proceedings of the 7th International Conference on Informatics in Control, Automation and Robotics, pages 187-192
DOI: 10.5220/0003000101870192
Copyright
c
 SciTePress

– extended pencil for DARE:
N − λM =


A 0 B
Q −E
H
L
L
H
0 R


− λ


E 0 0
0 −A
H
0
0 −B
H
0


.
If

U
T
1
U
T
2
U
T
3

T
spans the stable deflating sub-
space of N − λM, then X
∗
= U
2
(EU
1
)
−1
. The feed-
back gain matrix for the linear-quadratic optimal reg-
ulator can be computed directly via G = U
3
U
−1
1
.
If R is nonsingular, E = I
n
, and L = 0, the above
pencils can be reduced to 2n× 2n Hamiltonian and
symplectic pencils, respectively, by removing the sub-
pencils with infinite eigenvalues (Paige and Van Loan,
1981; Pappas et al., 1980; Mehrmann, 1991). A pen-
cil N −λM is Hamiltonian if NJ M
H
= −MJ N
H
, and
it is symplectic if NJ N
H
= MJ M
H
, where
J :=

0 I
n
−I
n
0

.
The general pencils inherit most of the spectral prop-
erties of the corresponding reduced Hamiltonian or
symplectic pencils.
The pencils above have much structure, which
should be exploited in order to improve the numer-
ical properties of the Riccati solvers. The approach
we follow is to transform the discrete-time problem
to an equivalent continuous-time problem, and use
the newly developed skew-Hamiltonian/Hamiltonian
eigensolvers for the latter problem, suitably extended.
2 EQUIVALENCE OF PENCILS
IN CONTINUOUS-TIME AND
DISCRETE-TIME PROBLEMS
A block column permutation (and sign change) gives,
equivalently:
– extended pencil for CARE:
λ


0 E 0
−E
H
0 0
0 0 0


−


0 A B
A
H
Q L
B
H
L
H
R


;
– extended pencil for DARE:
λ


0 E 0
−A
H
0 0
−B
H
0 0


−


0 A B
−E
H
Q L
0 L
H
R


.
These pencils are special cases of the following block
structured C-type and D-type pencils (Xu, 2006):
λE
C
− A
C
= λ

0
e
F
−
e
F
H
0

−

0
e
G
e
G
H
e
D

, (1)
and
λE
D
− A
D
= λ

0 F
−G
H
0

−

0 G
−F
H
D

, (2)
respectively, where F,G,
e
F,
e
G ∈ C
n,q
, q = n+ m, and
D,
e
D ∈ C
q,q
are Hermitian, i.e., D = D
H
,
e
D =
e
D
H
.
These pencils have important spectral properties:
C-type: symmetry about ℜ(z) = 0, i.e., pairs (λ,−
¯
λ);
D-type: symmetry about |z| = 1, i.e., pairs (λ,
¯
λ
−1
).
An equivalence transformation between the C-
type and D-type pencils can be established starting
from the Cayley transformation, c : C ∪ {∞} → C ∪
{∞}, defined by
µ = c(λ) = (λ − 1)(λ + 1)
−1
; c(−1) = ∞, c(∞) = 1.
Specifically, the generalized Cayley transformation
for matrix pairs is given by
(F ,B) = c(E,A) = (A + E,A − E). (3)
Let
(
e
E,
e
A) := c(E
D
,A
D
).
Then, an eigenvalue pair (λ,
¯
λ
−1
) of λE
D
− A
D
is
transformed to an eigenvalue pair (µ,−¯µ) of λ
e
E −
e
A,
with µ = c(λ), −¯µ = c(
¯
λ
−1
).
Unfortunately, λ
e
E −
e
A has not the same block
structure as λE
C
− A
C
, and it cannot be put into the
continuous-time setting. This inconvenience can be
removed using the Cayley transformation followed by
a drop/add transformation (Xu, 2006):
(E
C
,A
C
) = t(E
D
,A
D
),
where t(·) = d(c(·)), and d corresponds to drop-
ping/adding D in the E part.
The Cayley plus drop/add transformation is sug-
gestively represented by the following t transforma-
tion diagram:
λE
D
− A
D
= λ

0 F
−G
H
0

−

0 G
−F
H
D

c ↓↑ c
−1
λ
e
E −
e
A = λ

0
e
F
−
e
F
H
D

−

0
e
G
e
G
H
D

drop D from
e
E ↓↑ add D to
e
E
λE
C
− A
C
= λ

0
e
F
−
e
F
H
0

−

0
e
G
e
G
H
e
D

where
e
F := G + F,
e
G := G − F,
e
D = D. It is worth
mentioning that the t transformation involves matrix
additions and subtractions only.
Only regular pencils are considered in the sequel.
A pencil λE − A is regular if E and A are square and
det(γE − A) 6= 0 for some γ ∈ C. A necessary regu-
larity condition is: if the C-type and D-type pencils
of order n+ q are regular, then
q− rank
b
D ≤ n ≤ q,
ICINCO 2010 - 7th International Conference on Informatics in Control, Automation and Robotics
188

where
b
D =
e
D and
b
D = D, for C-type and D-type pen-
cils, respectively.
The relation between the eigen-structure of λE −
A and λF −B, (F ,B) = t(E,A), can be summarized
as follows (see, e.g., (Mehrmann, 1991; Xu, 2006)):
(i) λE − A is regular if and only if (iff) λF − B is
regular.
(ii) λ ∈ Λ(E,A) iff µ = c(λ) ∈ Λ(F ,B), and λ and µ
have the same geometric, partial, and algebraic multi-
plicities.
(iii) If λE −A is regular, then, R
λ
= R
µ
, L
λ
= L
µ
,µ =
c(λ), where R
x
and L
x
are the right and left deflating
subspaces corresponding to eigenvalue(s) x.
The C-type pencil (1) is skew-Hermitian/Hermi-
tian, i.e., E
H
C
= −E
C
,A
H
C
= A
C
, and it has the fol-
lowing main eigen-structure properties:
(i) λ ∈ Λ(E
C
,A
C
) iff −
¯
λ ∈ Λ(E
C
,A
C
), and λ and −
¯
λ
have the same geometric, partial, and algebraic multi-
plicities.
(ii) R
λ
= L
−
¯
λ
and L
λ
= R
−
¯
λ
.
(iii) U is a basis matrix of a right deflating subspace
of λE − A corresponding to λS− T iff U is a basis
matrix of a left deflating subspace corresponding to
λ(−S
H
) − T
H
.
The eigenvalue pairing (λ,−
¯
λ) does not hold for
λ with ℜ(λ) = 0, since then λ = −
¯
λ. But for such an
eigenvalue, R
λ
= L
λ
. This also holds for λ = ∞.
The regular D-type pencil (2) has the following
main eigen-structure properties (Mehrmann, 1991;
Xu, 2006):
(i) Nonzero finite eigenvalues come in pairs (λ,
¯
λ
−1
),
and λ,
¯
λ
−1
have the same geometric, partial, and al-
gebraic multiplicities.
(ii) span U = R
λ
, span V = L
λ
iff span
b
V = R
¯
λ
−1
,
span
b
U = L
¯
λ
−1
, where span X denotes the subspace
spanned by the columns of the matrix X, and
U =

U
1
U
2

, V =

V
1
V
2

∈ C
n+q,ℓ
,
E
D
UT = A
D
U, S
H
V
H
E
D
= V
H
A
D
, (4)
for T, S ∈ C
ℓ,ℓ
with Λ(T) = Λ(S
H
) = {λ}, λ 6= 0 (with
algebraic multiplicity ℓ), and
b
U =

U
1
T
U
2

,
b
V =

V
1
S
V
2

,
E
D
b
VS
−1
= A
D
b
V, T
−H
b
U
H
E
D
=
b
U
H
A
D
. (5)
Moreover, det V
H
E
D
U 6= 0 iff det
b
U
H
E
D
b
V 6= 0.
(iii) U =

U
T
1
U
T
2

T
is a basis matrix of a right
deflating subspace (left deflating subspace) of λE
D
−
A
D
corresponding to T ∈ C
p,p
nonsingular, iff
b
U =

(U
1
T)
T
U
T
2

T
is a basis matrix of a left deflat-
ing subspace (right deflating subspace) of λE
D
− A
D
corresponding to T
−H
.
(iv) If 0 ∈ Λ(E
D
,A
D
) with algebraic multiplicity
ℓ
0
, then ∞ ∈ Λ(E
D
,A
D
) with algebraic multiplicity
greater than or equal to ℓ
0
.
(v) The formulas for the relations between the basis
matrices of right/left deflating subspace for λ = 0 or
∞ are more complicated than for the case λ 6= 0,∞.
The sizes of the submatrices depend on the algebraic
multiplicities of 0 ∈ Λ(G
H
,F
H
) and 0 ∈ Λ(F,G).
The eigenvalue pairing (λ,
¯
λ
−1
) does not hold for
|λ| = 1, since then λ =
¯
λ
−1
. But for such an eigen-
value, U in (4) is a basis matrix of R
λ
iff
b
U in (5) is a
basis matrix of L
λ
.
Eigenvalues 0 and ∞ are paired in a weak sense,
since the algebraic multiplicity of ∞ may be greater
than or equal to the algebraic multiplicity of 0, and
R
0
and L
0
are only related to certain subspaces of L
∞
and R
∞
, respectively.
The equivalence relation between D-type and C-
type pencils is shown below.
Assume (E
C
,A
C
) = t(E
D
,A
D
) and that λE
D
−A
D
(or λE
C
− A
C
) is regular. Then,
(i) λ ∈ Λ(E
D
,A
D
), λ 6= −1,∞, iff µ = c(λ) ∈
Λ(E
C
,A
C
), µ 6= ∞,1, and λ and µ have the same geo-
metric, partial, and algebraic multiplicities.
(ii) span U = R
D
λ
, span V = L
D
λ
iff span
e
U = R
C
µ
,
span
e
V = L
C
µ
, where the superscript C or D refers
to (1) or (2), respectively, U and V satisfy (4) for
T,S ∈ C
ℓ,ℓ
with Λ(T) = Λ(S
H
) = {λ} (with algebraic
multiplicity ℓ), and
e
U =

U
1
(I + T)
2U
2

,
e
V =

V
1
(I + S)
2V
2

,
E
C
e
U
e
T = A
C
e
U,
e
S
H
e
V
H
E
C
=
e
V
H
A
C
,
where
e
T = c(T),
e
S = c(S), Λ(
e
T) = Λ(
e
S
H
) = {µ}, µ =
c(λ). Moreover, det V
H
E
D
U 6= 0 iff det
e
V
H
E
C
e
U 6= 0.
(iii) If −1 ∈ Λ(E
D
,A
D
), with algebraic multiplic-
ity ℓ
−1
, then ∞ ∈ Λ(E
C
,A
C
), with algebraic multi-
plicity greater than or equal to ℓ
−1
. Suppose also
−1 ∈ Λ(G
H
,F
H
), with algebraic multiplicity r
1
. Let
U =

U
11
U
12
0 U
22

∈ C
n+q,ℓ
−1
, U
11
∈ C
n,r
1
,
E
D
UT = A
D
U, rank E
D
U = ℓ
−1
,
T =

T
11
T
12
0 T
22

∈ C
ℓ
−1
,ℓ
−1
, T
11
∈ C
r
1
,r
1
,
with Λ(T) = {−1}. IfU is a basis matrix of R
D
−1
, then
the columns of
e
U =

2U
11
U
12
(T
22
+ I)
0 2U
22

STRUCTURE-PRESERVING ALGORITHMS FOR DISCRETE-TIME ALGEBRAIC MATRIX RICCATI EQUATIONS
189

span an ℓ
−1
-dimensional(right and left) deflating sub-
space of λE
C
− A
C
corresponding to eigenvalue ∞.
(iv) Let ℓ
−1
, ℓ
0
, and ℓ
∞
be the algebraic multiplic-
ities of the eigenvalues −1,0, ∞ ∈ Λ(E
D
,A
D
) and
e
ℓ
1
,
e
ℓ
∞
the algebraic multiplicities of the eigenvalues
1,∞ ∈ Λ(E
C
,A
C
). Then, ℓ
0
=
e
ℓ
1
and
e
ℓ
∞
= ℓ
∞
− ℓ
0
+ ℓ
−1
, ℓ
∞
=
e
ℓ
∞
− ℓ
−1
+
e
ℓ
1
.
Specifically, with t, ∞ ∈ Λ(E
C
,A
C
) comes from
−1 ∈ Λ(E
D
,A
D
) (with multiplicity ℓ
−1
) and ∞ ∈
Λ(E
D
,A
D
) (with multiplicity ℓ
∞
− ℓ
0
).
If −1 ∈ Λ(E
C
,A
C
) (i.e., 0 ∈ Λ(E
D
,A
D
)), then 1 ∈
Λ(E
C
,A
C
), and it comes from ∞ ∈ Λ(E
D
,A
D
). Only
part of ∞ ∈ Λ(E
D
,A
D
) is transformed into 1, to match
−1 ∈ Λ(E
C
,A
C
).
3 DEFLATING SUBSPACES FOR
SKEW-HAMILTONIAN/
HAMILTONIAN PENCILS
The structure-preserving algorithms and software are
more advanced for CAREs, based on deflating sub-
spaces for skew-Hamiltonian/Hamiltonian pencils.
Extensions of the HAPACK approach are currently
under development. In the sequel, the pencils λM −N
will be represented in the numerically better form
αM − βN, with λ = α/β (possibly ∞).
Since the structured algorithms for skew-Hamil-
tonian/Hamiltonian pencils work on problems with
even size, a basic idea is to embed the matrix pen-
cil, adding k ≥ 0 fictitious controls, so that m + k is
even. The solution of the optimal control problem
corresponding to CARE, hence to
αE
c
−βA
c
= α


E 0 0
0 −E
H
0
0 0 0


−β


A 0 B
Q A
H
L
L
H
B
H
R


,
is unchanged for k new controls, with
e
B = 0
n×k
,
e
R = I
k
, and D replaced by block-diag(D,
e
R), with
D :=

Q L
L
H
R

.
Partition, with ℓ = (m + k)/2, B
i
∈ C
n×ℓ
, L
i
∈
C
n×ℓ
, R
ij
∈ C
ℓ×ℓ
, i, j = 1,2,

B
e
B

=

B
1
B
2

,


Q L 0
L
H
R 0
0 0
e
R


=


Q L
1
L
2
L
H
1
R
11
R
12
L
H
2
R
21
R
22


.
Reordering the variables and equations, the following
skew-Hamiltonian/Hamiltonian pencil is obtained
αE
e
c
− βA
e
c
= α



E 0 0 0
0 0 0 0
0 0 E
H
0
0 0 0 0



− β




A B
1
0 B
2
L
H
2
R
H
12
B
H
2
R
22
−Q −L
1
−A
H
−L
2
−L
H
1
−R
11
−B
H
1
−R
12




. (6)
Let αS − βH be skew-Hamiltonian/Hamiltonian, i.e.,
(SJ )
H
= −SJ , (H J )
H
= H J . For some cases,
including in linear-quadratic optimization applica-
tions, S is given in a factored form, the so-called
skew-Hamiltonian Cholesky factorization, defined by
S = J Z
H
J
T
Z. For instance, in (6) with S = E
e
c
,
Z =



I
n
0 0 0
0 I
ℓ
0 0
0 0 E
H
0
0 0 0 0



.
A skew-Hamiltonian matrix having such a factoriza-
tion is said to be J -semidefinite.
Some properties of skew-Hamiltonian/Hamilto-
nian pencils are summarized below (Benner et al.,
2002):
(i) Real skew-Hamiltonian matrices are J -semide-
finite.
(ii) The structure of skew-Hamiltonian/Hamiltonian
matrix pencils is preserved under J -congruence:
αS − βH → J Y
H
J
T
(αS − βH )Y ,
for Y nonsingular.
(iii) A skew-Hamiltonian matrix S of order 2n is J -
semidefinite (J -definite) iff ıJ S has at most (exactly)
n positive and at most (exactly) n negative eigenval-
ues, where ı := (−1)
1/2
.
(iv) If S is skew-Hamiltonian (H is Hamiltonian) and
there is Y nonsingular, such that
J Y
H
J
T
SY =

S
11
S
12
0 S
H
11

(J Y
H
J
T
H Y =

H
11
H
12
0 −H
H
11

)
with S
11
,S
12
(H
11
,H
12
) ∈ C
n×n
, then S (ıH ) is J -
semidefinite.
(v) Let αS − βH be regular skew-Hamiltonian/
Hamiltonian with ν pairwise distinct, nonzero finite
eigenvalues ıα
i
, of algebraic multiplicity p
i
, and as-
sociated right deflating subspace Q
i
, i = 1 : ν; let p
∞
,
Q
∞
, be defined similarly for eigenvalue ∞. The fol-
lowing statements are equivalent:
ICINCO 2010 - 7th International Conference on Informatics in Control, Automation and Robotics
190

(a) There exists a nonsingular matrix Y , such that
J Y
H
J
T
(αS − βH )Y =
α

S
11
S
12
0 S
H
11

− β

H
11
H
12
0 −H
H
11

, (7)
where S
11
and H
11
are upper triangular, S
12
skew-
Hermitian, and H
12
Hermitian.
(b) There exists a unitary matrix Q , such that (7)
holds for Y replaced by Q .
(c) Q
H
k
J S Q
k
is congruent to a p
k
× p
k
copy of J ,
k = 1,2,. ..,ν; Q
H
∞
J H Q
∞
is congruent to a p
∞
× p
∞
copy of ıJ .
(vi) Factored version: Let αS − βH be a skew-
Hamiltonian/Hamiltonian pencil with nonsingular J -
semidefinite skew-Hamiltonian part S = J Z
H
J
T
Z. If
any of the equivalent statements above holds, then
there is a unitary matrix Q and a unitary symplectic
matrix U, such that
U
H
ZQ =

Z
11
Z
12
0 Z
22

,
J Q
H
J
T
H Q =

H
11
H
12
0 −H
H
11

,
where Z
11
, Z
T
22
, and H
11
are n× n upper triangular.
(vii) If ıH is also nonsingular J -semidefinite, i.e.,
ıH = J W
H
J
T
W , then
U
H
ZQ =

Z
11
Z
12
0 Z
22

,
U
H
W Q =

W
11
W
12
0 W
22

,
where Z
11
, Z
T
22
, W
11
, and W
T
22
are n× n upper triangu-
lar.
(viii) Factored version, real skew-Hamiltonian/skew-
Hamiltonian case: Let αS − βN be a real regu-
lar skew-Hamiltonian/skew-Hamiltonian pencil with
S = J Z
T
J
T
Z. Then, there is an orthogonal matrix Q
and an orthogonal symplectic matrix U, such that
U
T
ZQ =

Z
11
Z
12
0 Z
22

,
J Q
T
J
T
N Q =

N
11
N
12
0 N
T
11

,
where Z
11
, Z
T
22
are upper triangular, N
11
is upper
quasi-triangular, and N
12
= −N
T
12
. Moreover,
J Q
T
J
T
(αS − βN )Q =
α

Z
T
22
Z
11
Z
T
22
Z
12
− Z
T
12
Z
22
0 Z
T
11
Z
22

− β

N
11
N
12
0 N
T
11

is a J -congruent skew-Hamiltonian/skew-Hamilto-
nian matrix pencil.
Comments:
1. The result (viii) above is used to compute the struc-
tured Schur form of order 4n for a complex skew-
Hamiltonian/Hamiltonian pencil.
2. The periodic QZ algorithm is used.
3. Algorithms for eigenvalue reordering and deflating
subspace computation are available.
Below is a summary about the related software:
• Fortran and MATLAB software for eigenvalues and
deflating subspaces are under development.
• Both real and complex cases are considered.
• Factored or unfactored versions are covered.
• Auxiliary routines for problems (of even order) with
(quasi-)triangular structure are included.
• Optimized kernels for problems of order 2, 3, or 4,
called by the general solvers, are available.
To compute or reorder the eigenvalues, the com-
putations begin with an initial reduction, called gen-
eralized symplectic URV decomposition, whose fac-
tored version is defined as follows (Benner et al.,
2007):
Given a real skew-Hamiltonian/Hamiltonian 2n × 2n
pencil αT Z − βH , orthogonal matrices Q
1
, Q
2
and
orthogonal symplectic matrices U
1
, U
2
are deter-
mined, such that
Q
T
1
T U
1
=

T
11
T
12
0 T
22

,
U
T
2
ZQ
2
=

Z
11
Z
12
0 Z
22

,
Q
T
1
H Q
2
=

H
11
H
12
0 H
22

,
where T
11
, T
T
22
, Z
11
, Z
T
22
, and H
11
are upper triangular,
and H
T
22
is upper quasi-triangular. The matrices U
1
and U
2
are stored compactly (the first n rows only),
since, for i = 1,2,
U
i
=

U
i1
U
i2
−U
i2
U
i1

.
4 NUMERICAL RESULTS
This section presents some preliminary numerical re-
sults. These results have been obtained on a portable
Intel Dual Core computer at 2 GHz, with 2 GB
RAM, and relative machine precision ε ≈ 1.11 ×
10
−16
, using Windows XP (Service Pack 2) operat-
ing system, Intel Visual Fortran 11.1 compiler, and
MATLAB 7.8.0.347 (R2009a). The matrices
S =

A B
C A
H

, H =

D E
F −D
H

,
where A, B, C, D, E, F ∈ C
n×n
, have been initially
generated with MATLAB commands of the form
STRUCTURE-PRESERVING ALGORITHMS FOR DISCRETE-TIME ALGEBRAIC MATRIX RICCATI EQUATIONS
191

list
A = 10*rand(n)-5 + (10*rand(n)-5)*1i;
where rand is the uniform (0,1) random generator,
and 1i is the MATLAB notation for the purely imag-
inary unit, ı. Then, the B, C, E, and F matrices have
been transformed using the formulas
B := B− B
H
, B := B/2; C := C −C
H
,
E := E + E
H
, E := E/2; F := F + F
H
,
to become skew-Hermitian, and Hermitian, respec-
tively. Therefore, the pencil λS − H is skew-
Hamiltonian/Hamiltonian.
The order n took the values n = 100,200, ...,
800. For each order n ≤ 500, 10 problems have been
solved, and the means of the results are reported.
For larger n values, one problem has been solved for
each n. The generalized eigenvalues computed by a
structure-preserving algorithm have been compared
with those delivered by the standard QZ algorithm,
optimally implemented in the MATLAB function eig .
Fig. 1 presents the ratios of the mean CPU times,
in seconds, i.e., the speed-up factor of the structured
algorithm, in comparison with the standard algorithm.
0 200 400 600 800
1
2
3
4
5
6
7
n
Time eig/Time structured alg.
Comparison of CPU times for eig / structured algorithm
Figure 1: Ratios between the CPU times needed by
the MATLAB function eig and the structure-preserving
algorithm for randomly generated complex skew-
Hamiltonian/Hamiltonian pencils of order 2n.
The deviation from symmetry of the eigenval-
ues computed by eig has also been computed as
the difference between the vector of eigenvalues λ =
[λ
1
,λ
2
,... ,λ
2n
]
T
and a permutation of the elements
of the vector −
¯
λ, chosen so that the elements with
the same indices in the two vectors be as close as
possible. The largest norm has been 4 · 10
−10
, and
the smallest norm has been 1.90 · 10
−12
. The norms
should theoretically be 0.
5 CONCLUSIONS
Main issues related to the structure-preserving algo-
rithms for solving discrete-time algebraic matrix Ric-
cati equations are summarized. Stable deflating sub-
spaces for extended, inverse-free symplectic matrix
pencils, are computed. Algorithms based on skew-
Hamiltonian/Hamiltonian pencils derived by an ex-
tended Cayley transformation, which only involves
matrix additions and subtractions, are considered.
The preliminary results are encouraging.
ACKNOWLEDGEMENTS
The work was partially supported by the German Re-
search Foundation (DFG) and The MathWorks, Inc.
The collaboration with Peter Benner and Matthias
Voigt from TU Chemnitz is highly acknowledged.
REFERENCES
Bender, D. J. and Laub, A. J. (1987a). The linear-quadratic
optimal regulator for descriptor systems. IEEE Trans.
Automat. Contr., AC-32(8):672–688.
Bender, D. J. and Laub, A. J. (1987b). The linear-quadratic
optimal regulator for descriptor systems: Discrete-
time case. Automatica, 23(1):71–85.
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.
Lancaster, P. and Rodman, L. (1995). The Algebraic Riccati
Equation. Oxford University Press, Oxford.
Mehrmann, V. (1991). The Autonomous Linear Quadratic
Control Problem. Theory and Numerical Solution,
volume 163 of Lect. Notes in Control and Information
Sciences. Springer-Verlag, Berlin.
Paige, C. and Van Loan, C. F. (1981). A Schur decomposi-
tion for Hamiltonian matrices. Lin. Alg. Appl., 41:11–
32.
Pappas, T., Laub, A. J., and Sandell, N. R. (1980). On
the numerical solution of the discrete-time algebraic
Riccati equation. IEEE Trans. Automat. Contr., AC–
25(4):631–641.
Van Dooren, P. (1981). A generalized eigenvalue approach
for solving Riccati equations. SIAM J. Sci. Stat. Com-
put., 2(2):121–135.
Xu, H. (2006). On equivalence of pencils from discrete-time
and continuous-time control. Lin. Alg. Appl., 414:97–
124.
ICINCO 2010 - 7th International Conference on Informatics in Control, Automation and Robotics
192