NUMERICAL LINEAR ALGEBRA WITH APPLICATIONS

Numer. Linear Algebra Appl. 2008; 15:747754

Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nla.621

EDITORIAL

Large-scale matrix equations of special type

GUEST EDITOR: Peter Benner, ,

Fakultat fur Mathematik, TU Chemnitz, 09107 Chemnitz, Germany

SUMMARY
This issue contains papers dealing with a variety of linear and nonlinear matrix equations arising in
systems and control theory, model reduction, and in various other areas of applied mathematics, economics,
engineering, and the sciences. Special emphasis is given to the numerical treatment of large-scale problems.
Copyright q 2008 John Wiley & Sons, Ltd.

KEY WORDS: matrix equation; control theory; systems theory; model reduction; Sylvester equation;
Lyapunov equation; algebraic Riccati equation; differential Riccati equation

This special issue of Numerical Linear Algebra with Applications is devoted to Large-scale matrix
equations of special type. The term matrix equations is generally used for situations where the
solution of a mathematical equation is given as a matrix (as opposed to a single vector). In this
sense, over some field F (usually, F {R, C}),

AX = C for given A Fnm , C Fn p (1)

with unknown X Fm p is a simple linear matrix equation. Such equations are not the topic of
this issue as (1) can be solved by usual methods for linear systems of equations or linear least-
squares problems, maybe with some exploitation of the fact that (1) represents p systems of linear
equations (or least-squares problems) sharing the same coefficient matrix A. Thus, (1) does not
pose any extra challenge compared with standard linear systems of equations.
The term special type in the title of this special issue reflects the fact that the contained papers
deal with matrix equations having certain structures arising in diverse applications. The simplest
equation considered here is the Sylvester equation

AX +XB = C for given A Fnn , B Fmm , C Fnm (2)

Correspondenceto: Peter Benner, Fakultat fur Mathematik, TU Chemnitz, 09107 Chemnitz, Germany.

E-mail:[email protected]

Professur Mathematik in Industrie und Technik.

Copyright q 2008 John Wiley & Sons, Ltd.

748 EDITORIAL

which arises in many decoupling problems as the solution matrix X Fnm can be used for block
diagonalization when embedded into a similarity transformation as follows:





In 0 B 0 In 0 B 0
= (3)
X Im C A X Im 0 A

This property can be exploited, e.g. for decoupling slow and fast states or stiff and non-stiff
components in linear(ized) dynamical systems, see e.g. [1, 2]. In numerical linear algebra, one of
the most elementary applications of Sylvester equations is the computation of an eigenvector basis
from a given Schur decomposition of a matrix [3, Section 7.6]. Among the many more applications
of (2), the ones in systems and control theory are probably the most commonly encountered (see
below). There the symmetric version of (2) with n = m, B = A , C = C , called Lyapunov equation,
also plays a major role.
Several nonlinear matrix equations have similar areas of applications as the Sylvester equation.
These are, for example, (non-symmetric) algebraic Riccati equations (AREs)

AX +XB+XDX = C for given A Fnn , B Fmm , C Fnm , D Fmn (4)

Their solutions X Fnm can for example be used for block triangularization using a similarity
transformation as follows:





In 0 B D In 0 B +DX D
= (5)
X Im C A X Im 0 (A +XD)
   
An important consequence of (5) is that the columns of IXn form an invariant subspace of CB AD
.
This observation forms the basis for solving a large variety of application problems by the use
of AREs. Similar to Sylvester and Lyapunov equations, the probably most prominent applications
of (4) can be found in systems and control theory, where a symmetric version (n = m, B = A ,
C = C , D = D ) of (4), commonly called continuous-time algebraic Riccati equation (CARE),
is omnipresent. In addition, in control theory for discrete-time linear systems, rational matrix
equations like the discrete-time ARE can be found.
Besides algebraic matrix equations, differential or difference matrix equations are encountered
in many application fields, including the already mentioned ones. For example, Riccati differential
equations have the form
d
X = AX +XB+XDX C for given A Fnn , B Fmm , C Fnm , D Fmn (6)
dt
where now, X X (t) is time dependent and also the coefficient matrices A, B, C, D may be matrix-
valued functions of t with suitable smoothness properties. These equations have already been the
topic of a monograph several decades ago [4], but remain an interesting and challenging research
topic. Sylvester and Lyapunov differential equations can be derived from (6) in an obvious way. In

Note:usually, block diagonalization is performed for block upper triangular matrices, which corresponds to a
transposed version of (3). We use the block lower triangular form in (3) in order to achieve coherence with (5).

Copyright q 2008 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. 2008; 15:747754
DOI: 10.1002/nla
 EDITORIAL 749

this context, the solutions of the algebraic matrix equations (2) and (4) represent the steady states
corresponding to (6).
As already indicated above, matrix equations have played a fundamental role in systems theory,
control, filtering, and estimation since the epochal discoveries of Kalman in the early sixties: the
CARE became known as the main tool to solve the linear-quadratic regulator problem and to find
the stationary Kalman filter, which together constitute the linear-quadratic Gaussian controller.
As new control methodologies like H2 - and H -control emerged in the early eighties, the major
work horse in computing the associated optimal or suboptimal controllers was still found to be the
CARE. In addition, linear matrix equations such as the Sylvester and Lyapunov equations occur in
many analysis and design methods for linear and nonlinear control systems, e.g. in observer design,
model reduction or stabilization, the H2 norm of a system is computed via solving a Lyapunov
equation, etc. There are numerous other applications of linear and nonlinear matrix equations. The
importance of the topic is emphasized in many textbooks and monographs that are completely
devoted to the theory and numerical solution of these matrix equations [58] or deal in large parts
with this topic [914], not counting the many texts on control theory and the numerous survey
papers.
As new design and modeling strategies evolve, such as the use of periodic systems, stochastic
systems or differential games, new aspects and variants of matrix equations are discovered. E.g. the
equations investigated in Damms paper in this issue are related to stochastic systems and play the
same role as Lyapunov equations play for classical (deterministic) linear time-invariant systems.
Periodic systems arise in various application areas describing, e.g. rotating objects like helicopter
blades, rolling tires or orbital dynamics as used in control of satellites. The formulations of the
usual control problems like linear-quadratic optimization over finite time horizons or H -control
for such systems lead to periodic Riccati differential equations of the form (6), where A, B, C, D, X
are periodic matrix functions, i.e. A(t) = A(t + T ), etc., for all t R and some fixed T >0. The
numerical solution of such equations is the topic of Vargas paper in this issue.
Besides special type, the term large-scale is part of the title of this issue. Solving matrix
equations tends to be a challenging large-scale problem even if n, m are relatively small due to the
fact that (2), (4), and (6) represent linear or nonlinear systems of equations with nm unknowns. A
direct approach for solving (2) employs the associated linear system of equations
((Im A)+(B T In ))vec(X ) = vec(C) (7)
where In , Im are identity matrices of size n, m, respectively, denotes the Kronecker (tensor)
product, and vec(M) is the vector obtained from stacking the columns of the matrix M on top of
each other (see e.g. [15]). The solution of (7) using Gaussian elimination would require 23 (nm)3
floating point operations. On a current desktop computer, this becomes already infeasible for n =
m = 100, which can be regarded as a fairly small problem size. Fortunately, there exist algorithms
that exploit the structure of (2) and have a complexity of O(max{n, m}3 ) only. Similar observations
hold for Lyapunov and Riccati equations. However, these algorithms are also of limited use if the
coefficient matrices of the matrix equation to be solved come from a discretized partial differential
equation (PDE). This issue arises, e.g. in linear-quadratic optimal control problems for PDEs, see
e.g. [1619]. In this application area, CAREs with n in the thousands up to the millions have
to be solved. There has been considerable progress in this direction, see e.g. [2025], but with
increasing challenges for PDE control problems, there is further need for more efficient and new
approaches for solving AREs/CAREs and also Sylvester and Lyapunov equations. The papers of
Baur, Damm, and Grasedyck in this issue address this need while the paper of Benner, Li, and

Copyright q 2008 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. 2008; 15:747754
DOI: 10.1002/nla
750 EDITORIAL

Penzl was part of the progress of the last decade as it is the reprint of an unpublished manuscript
from 1999, see below for details.
Another important driving force for many developments in the area of numerical methods for
Lyapunov and Sylvester equations in recent years is model reduction for linear, time-invariant
systems, see e.g. [2628]. The papers of Baur and Damm will find their applications particularly in
this area while the paper by Zhou and Sorensen directly addresses the model reduction problem.
In what follows, we give a short introduction of the contents of the papers in this issue.

Numerical solution of large-scale Lyapunov equations, Riccati equations, and linear-quadratic

optimal control problems by Peter Benner, Jing-Rebecca Li, and Thilo Penzl
In this paper, large-scale algebraic Lyapunov and Riccati equations are studied together with the
associated linear-quadratic optimal control problem. The main workhorse in this paper is an ADI-
based method to compute approximate low-rank Cholesky factors for the solutions of large-scale
Lyapunov equations. This is then used in the NewtonKleinman iteration [29] for solving CAREs.
Closely inspecting the resulting iteration, it is possible to compute the optimal feedback matrix
defining the optimal control function without constructing the solution of the CARE or a factor
of it. This is achieved by directly updating the feedback matrix during the ADI iteration. This
resembles a technique used in earlier work [20], with the difference that there, a different solver
for Lyapunov equations and a different version of Newtons method for CAREs is employed.
Notice that this is an unabridged, corrected reprint of an unpublished manuscript written by the
authors in 1999. The latest version of this work was finalized on 8 December 1999. Thilo Penzl
died in a tragic avalanche accident on 17 December 1999, and thus, work on the manuscript came
to an abrupt end. Nevertheless, this work was often cited and requested by people working on
large-scale matrix equations as it forms the basis for parts of the software package LyaPack [30].
This special issue offered the possibility to make this draft available in the open literature. No
developments since 2000 except for updated references are taken into account.
Owing to its nature as a reprint of a formerly unpublished draft, some parts of this work might
not appear to be new as they have been published in different papers since 2000. In particular, the
low-rank Cholesky factor ADI iteration was studied in much more detail in [31] and the solution
of CAREs and linear-quadratic optimal control problems was discussed in [23, 24], using parts of
this work.

Nonlinear multigrid for the solution of large scale Riccati equations in low-rank and H-matrix
format by Lars Grasedyck
The numerical solution of large-scale CAREs is of fundamental importance in linear-quadratic
optimal control problems for parabolic and hyperbolic PDEs as pointed out above. Often, the
discretization of the partial differential operator using finite elements or finite differences for
the spatial variables leads to a hierarchy of grids. This enables the usage of multigrid techniques
for the numerical solution of CAREs. This was already investigated in [21, 22] where a multigrid
solver for Lyapunov equations was embedded into the NewtonKleinman iteration for CAREs.
(Newtons method or its Kleinman version for CAREs [29] is based on solving a Lyapunov equation
in each Newton step.) In contrast to this, Grasedyck develops variants of nonlinear multigrid
methods that apply to the CARE directly. The two versions considered, Newton multigrid and
full nonlinear multigrid, are derived and tested for a finite-difference discretization of a parabolic
model control problem. The largest problem solved is of size n = 4 190 209, which corresponds

Copyright q 2008 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. 2008; 15:747754
DOI: 10.1002/nla
 EDITORIAL 751

to solving a system of nonlinear equations with more than 17 trillion unknowns or close to 9
trillion unknowns if symmetry is exploited! The method is also generalized to the case that C, D
are represented by hierarchical matrices rather than low-rank matrices.

On solving periodic Riccati equations by Andras Varga

Periodicity is encountered in many dynamical systems arising in the sciences and engineering. This
includes models of rotating or spinning objects as mentioned above. Control design for the resulting
systems leads in one way or another to the solution of matrix equations with periodic coefficient
matrices, either in a continuous setting to (6) with periodic matrix functions A(t), B(t), C(t), D(t)
or, after sampling/time-discretization of the underlying system, to discrete-time periodic Riccati
equations

X k = Ck + ATk X k+1 Ak (ATk X k+1 Bk + Sk )(Dk + BkT X k+1 Bk )1 (ATk X k+1 Bk + Sk )T (8)

with given matrix sequences Ak , Bk , Ck , Dk , Sk of appropriate dimensions and periodicity Ak =

Ak+K , Bk = Bk+K , . . . for all k N0 and some fixed K N.
Shooting methods are an often-used approach to solve periodic differential equations; thus, Varga
develops a new multiple shooting-type algorithm for periodic differential Riccati equations by
employing suitable discretizations of the continuous-time problems. For the discrete-time periodic
Riccati equation (8) he suggests several methods based either on extensions of the periodic QZ
algorithm [32, 33] to non-square periodic pairs or on an extension of a quotientproduct swapping
and collapsing method [34, 35]. Numerical experiments illustrate the efficiency and accuracy of
the proposed algorithms.

Low rank solution of data-sparse Sylvester equations by Ulrike Baur

The Sylvester equation (2) with a factorized right-hand side C = FG, F, G being of low rank
r  n, m, arises in various applications of the control theory, including model reduction based on
the so-called cross-Gramian approach (see e.g. [36]) and observer design [9]. The approach taken
here combines the idea of using a particular iteration based on the sign function for this kind of
Sylvester equations with the formatted arithmetic for hierarchical matrices (see e.g. [37]). In this
way, very large-scale Sylvester equations with matrices A, B, F, G resulting from finite-element
or boundary-element discretizations of control problems for partial differential equations can be
solved efficiently.

Direct methods and ADI-preconditioned Krylov subspace methods for generalized Lyapunov equa-
tions by Tobias Damm
In this paper, linear matrix equations related to stochastic or bilinear control systems are considered.
These equations generalize algebraic Lyapunov equations in the sense that terms of the form
A j X ATj , j = 1, . . . ,
, are added to the Lyapunov part AX + X AT . These equations cannot be
handled with standard algorithms for Sylvester or Lyapunov equations. Vectorization similar to (7)
is one possibility to treat these equations, but leads to the above discussed difficulties regarding
complexity. In case the considered equation is a low-rank perturbation of a Lyapunov equation, it

A trillion in the short-scale naming convention used in American and modern British English refers to 1012 .

Copyright q 2008 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. 2008; 15:747754
DOI: 10.1002/nla
752 EDITORIAL

is proposed to exploit this low-rank structure so that the solution can be obtained by solving a finite
sequence of standard Lyapunov equations. Alternatively, an iterative approach based on convergent
regular operator splittings, employing Krylov subspace methods like GMRES(k) or BiCGStab with
an ADI-iteration as preconditioner, is suggested and tested for several numerical examples.

Approximate implicit subspace iteration with alternating directions for LTI system model reduction
by Yunkai Zhou and Dan C. Sorensen
This paper is actually not about solving matrix equations, but about avoiding them. In some
sense the explicit computation of the solution of Lyapunov equations (which define the symmetric
positive semidefinite system Gramians P, Q Rnn of stable, linear time-invariant (LTI) systems)
as required in algorithms for balanced truncation and related model reduction methods is partially
avoided by algorithms computing approximate low-rank factors of the solution as those suggested
in [31, 3840] and in many other related publications. This paper goes one step further and avoids
even computing these solution factors by approximating the dominant eigenspaces of the product
PQ of the system Gramians directly. The investigated Approximate Implicit Subspace Iteration
with Alternating Directions (AISIAD) framework for LTI system model reduction draws from
earlier ideas of balancing-free balanced truncation [41], which also uses the dominant eigenspaces
of P Q and the approximate dominant subspace algorithms proposed in [42]. Different versions of
the AISIAD approach result from using different updating techniques for the computed subspaces.
These are discussed and compared using several numerical examples.
We hope that the papers in this special issue on matrix equations will have an impact on the
awareness that in recent years, a variety of methods to solve large-scale matrix equations in control
theory, model reduction, and other areas of applied mathematics and computational engineering
have been developed and are now ready for use. Moreover, this special issue should further the
development of new and improved numerical algorithms for solving matrix equations. Certainly,
the ideas presented here can also be applied to related matrix equations and areas, and also these
opportunities should and will be explored in the future.

REFERENCES
1. Enright W. Improving the efficiency of matrix operations in the numerical solution of stiff ordinary differential
equations. ACM Transactions on Mathematical Software 1978; 4:127136.
2. Epton M. Methods for the solution of AXDBXC = E and its application in the numerical solution of implicit
ordinary differential equations. BIT 1980; 20:341345.
3. Golub G, Van Loan C. Matrix Computations (3rd edn). Johns Hopkins University Press: Baltimore, 1996.
4. Reid W. Riccati Differential Equations. Academic Press: New York, 1972.
5. Abou-Kandil H, Freiling G, Ionescu V, Jank G. Matrix Riccati Equations in Control and Systems Theory.
Birkhauser: Basel, Switzerland, 2003.
6. Bittanti S. Proceedings of the Workshop on the Riccati Equation in Control, Systems, and Signals, Como, Italy.
Pitagora Editrice: Bologna, Italy, 1989.
7. Gajic Z, Qureshi M. Lyapunov matrix equation in system stability and control. Mathematics in Science and
Engineering. Academic Press: San Diego, CA, 1995.
8. Lancaster P, Rodman L. The Algebraic Riccati Equation. Oxford University Press: Oxford, 1995.
9. Datta B. Numerical Methods for Linear Control Systems. Elsevier Academic Press: San Diego, CA, U.S.A.,
London, U.K., 2004.
10. Halanay A, Ionescu V. Time-varying Discrete Linear Systems: InputOutput Operators, Riccati Equations and
Disturbance Attenuation, Operator Theory: Advances and Applications, vol. 68. Birkhauser: Basel, Switzerland,
1994.

Copyright q 2008 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. 2008; 15:747754
DOI: 10.1002/nla
 EDITORIAL 753

11. Ionescu V, Oara C, Weiss M. Generalized Riccati Theory and Robust Control. A Popov Function Approach.
Wiley: Chichester, 1999.
12. Mehrmann V. The Autonomous Linear Quadratic Control Problem, Theory and Numerical Solution. Lecture
Notes in Control and Information Sciences, vol. 163. Springer: Heidelberg, 1991.
13. Petkov P, Christov N, Konstantinov M. Computational Methods for Linear Control Systems. Prentice-Hall:
Hertfordshire, U.K., 1991.
14. Sima V. Algorithms for Linear-Quadratic Optimization. Pure and Applied Mathematics, vol. 200. Marcel Dekker:
New York, NY, 1996.
15. Lancaster P, Tismenetsky M. The Theory of Matrices (2nd edn). Academic Press: Orlando, 1985.
16. Curtain R, Zwart H. An Introduction to Infinite-Dimensional Linear Systems Theory, Texts in Applied Mathematics,
vol. 21. Springer: New York, 1995.
17. Lasiecka I, Triggiani R. Control Theory for Partial Differential Equations: Continuous and Approximation Theories
I. Abstract Parabolic Systems. Cambridge University Press: Cambridge, U.K., 2000.
18. Lasiecka I, Triggiani R. Control theory for partial differential equations: continuous and approximation theories
II. Abstract hyperbolic-like systems over a finite time horizon. Encyclopedia of Mathematics and its Applications,
vol. 75. Cambridge University Press: Cambridge, 2000; 6451067.
19. Bensoussan A, Da Prato G, Delfour MC, Mitter SK. Representation and control of infinite dimensional systems.
Systems and Control: Foundations and Applications (2nd edn). Birkhauser Boston Inc.: Boston, MA, 2007.
20. Banks H, Ito K. A numerical algorithm for optimal feedback gains in high dimensional linear quadratic regulator
problems. SIAM Journal on Control and Optimization 1991; 29(3):499515.
21. Rosen I, Wang C. A multi-level technique for the approximate solution of operator Lyapunov and algebraic
Riccati equations. SIAM Journal on Numerical Analysis 1995; 32(2):514541.
22. Penzl T. A multi-grid method for generalized Lyapunov equations. Technical Report SFB393/97-24, Fakultat
fur Mathematik, TU Chemnitz, 09107 Chemnitz, FRG, 1997. Available from: http://www.tu-chemnitz.de/sfb393/
sfb97pr.html.
23. Benner P. Efficient algorithms for large-scale quadratic matrix equations. Proceedings in Applied Mathematics
and Mechanics 2002; 1(1):492495.
24. Benner P. Solving large-scale control problems. IEEE Control Systems Magazine 2004; 14(1):4459.
25. Morris K, Navasca C. Solution of algebraic Riccati equations arising in control of partial differential equations.
Control and Boundary Analysis. Lecture Notes in Pure and Applied Mathematics, vol. 240. Chapman & Hall/CRC:
Boca Raton, FL, 2005; 257280.
26. Antoulas A. Approximation of Large-Scale Dynamical Systems. SIAM Publications: Philadelphia, PA, 2005.
27. Benner P, Mehrmann V, Sorensen D. Dimension Reduction of Large-Scale Systems. Lecture Notes in Computational
Science and Engineering, vol. 45. Springer: Berlin/Heidelberg, Germany, 2005.
28. Benner P, Freund R, Sorensen D, Varga A. Special issue on order reduction of large-scale systems. Linear
Algebra and its Applications 2006; 415(23).
29. Kleinman D. On an iterative technique for Riccati equation computations. IEEE Transactions on Automatic
Control 1968; AC-13:114115.
30. Penzl T. Lyapack users guide. Technical Report SFB393/00-33, Sonderforschungsbereich 393 Numerische
Simulation auf massiv parallelen Rechnern, TU Chemnitz, 09107 Chemnitz, FRG, 2000. Available from:
http://www.tu-chemnitz.de/sfb393/sfb00pr.html.
31. Li JR, White J. Low rank solution of Lyapunov equations. SIAM Journal on Matrix Analysis and Applications
2002; 24(1):260280.
32. Bojanczyk A, Golub G, Van Dooren P. The periodic Schur decomposition. Algorithms and applications. In
Advanced Signal Processing Algorithms, Architectures, and Implementations III, Proceedings SPIE, San Diego,
CA, U.S.A., Luk F (ed.), vol. 1770, 1992; 3142.
33. Hench J, Laub A. Numerical solution of the discrete-time periodic Riccati equation. IEEE Transactions on
Automatic Control 1994; 39:11971210.
34. Benner P, Byers R. Evaluating products of matrix pencils and collapsing matrix products. Numerical Linear
Algebra with Applications 2001; 8(67):357380.
35. Benner P, Byers R. An arithmetic for matrix pencils: theory and new algorithms. Numerische Mathematik 2006;
103(4):539573.
36. Fernando K, Nicholson H. On a fundamental property of the cross-Gramian matrix. IEEE Transactions on
Circuits and Systems 1984; CAS-31(5):504505.
37. Grasedyck L, Hackbusch W. Construction and arithmetics of H-matrices. Computing 2003; 70(4):295334.

Copyright q 2008 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. 2008; 15:747754
DOI: 10.1002/nla
754 EDITORIAL

38. Jaimoukha I, Kasenally E. Krylov subspace methods for solving large Lyapunov equations. SIAM Journal on
Numerical Analysis 1994; 31:227251.
39. Benner P, Quintana-Ort E. Solving stable generalized Lyapunov equations with the matrix sign function. Numerical
Algorithms 1999; 20(1):75100.
40. Penzl T. A cyclic low rank Smith method for large sparse Lyapunov equations. SIAM Journal on Scientific
Computing 2000; 21(4):14011418.
41. Safonov M, Chiang R. A Schur method for balanced-truncation model reduction. IEEE Transactions on Automatic
Control 1989; AC34:729733.
42. Penzl T. Algorithms for model reduction of large dynamical systems. Linear Algebra and its Applications 2006;
415(23):322343 (Reprint of Technical Report SFB393/99-40, TU Chemnitz, 1999).

Copyright q 2008 John Wiley & Sons, Ltd. Numer. Linear Algebra Appl. 2008; 15:747754
DOI: 10.1002/nla