Sparsity Preserving Optimal Control of Discretized
PDE Systems

arXiv:1801.05194v3 [math.OC] 20 Feb 2018

Aleksandar Habera,∗, Michel Verhaegenb
a Department

of Engineering Science and Physics, City University of New York, College of
Staten Island, New York 10314, USA
b Delft Center for Systems and Control, Delft University of Technology, 2628 CD Delft, The
Netherlands

Abstract
We focus on the problem of optimal control of large-scale systems whose models
are obtained by discretization of partial differential equations using the Finite
Element (FE) or Finite Difference (FD) methods. The motivation for studying this pressing problem originates from the fact that the classical numerical
tools used to solve low-dimensional optimal control problems are computationally infeasible for large-scale systems. Furthermore, although the matrices of
large-scale FE or FD models are usually sparse banded or highly structured,
the optimal control solution computed using the classical methods is dense and
unstructured. Consequently, it is not suitable for efficient centralized and distributed real-time implementations. We show that the a priori (sparsity) patterns of the exact solutions of the generalized Lyapunov equations for FE and
FD models are banded matrices. The a priori pattern predicts the dominant
non-zero entries of the exact solution. We furthermore show that for wellconditioned problems, the a priori patterns are not only banded but also sparse
matrices. On the basis of these results, we develop two computationally efficient
methods for computing sparse approximate solutions of generalized Lyapunov
equations. Using these two methods and the inexact Newton method, we show
that the solution of the generalized Riccati equation can be approximated by
a banded matrix. This enables us to develop a novel computationally efficient
optimal control approach that is able to preserve the sparsity of the control law.
We perform extensive numerical experiments that demonstrate the effectiveness
of our approach.
Keywords: Finite Element Methods; Optimal control; Large-scale systems;
Riccati equation; Lyapunov equation.

∗ Corresponding

author
Email address: aleksandar.haber@gmail.com; aleksandar.haber@csi.cuny.edu
(Aleksandar Haber)

Preprint submitted to Computer Methods in Applied Mechanics and EngineeringFebruary 21, 2018

1. Introduction
Large-scale systems with the dynamics described by Partial Differential
Equations (PDEs) are ubiquitous in engineering and science. This comes from
the fact that the majority of the fundamental physical laws and processes are
mathematically described by PDEs. On the other hand, the problem of optimal control of large-scale systems is a long-standing and pressing problem
whose solution is important for efficient, safe, and low-cost operation of various
mechanical, electrical and physical systems [1, 2, 3, 4, 5, 6, 7, 8].
The focus of this paper is on the systems whose dynamics can be represented
in the descriptor state-space form [9, 10, 11]: E ẋ = Ax + Bu, where x ∈ Rn
is the state vector, u ∈ Rm is the control input vector, E, A ∈ Rn×n , and B ∈
Rn×m are the system matrices. The models where E = I are referred to as nondescriptor state-space models. In this paper, we mainly focus on the descriptor
state-space models resulting from the discretization of linear PDEs using the
classical Finite Element (FE) or Finite Difference (FD) methods. Under the
term the classical FE method, we understand a Galerkin FE method with locally
supported basis functions and with structured or unstructured meshes composed
of the triangular or rectangular 2D elements (tetrahedral and hexahedral 3D
elements). FE models of various classes of linear PDEs can be written in the
descriptor state-space form. For example, the general form of FE models used
for the dynamical analysis is [12]: M1 s̈ + M2 ṡ + M3 s = B1 u, where s, ṡ, and s̈
are the displacement, velocity, and acceleration vectors, and M1 , M2 , and M3
are the mass, damping and stiffness matrices, respectively. The term B1 u is the
vector of applied loads. The corresponding descriptor state-space model is




 
 
I
0
0
I
0
s
ẋ =
x+
u, x =
.
(1)
0 M1
−M3 −M2
B1
ṡ
Due to the local support of the basis functions, and taking into consideration
the fact that the node degrees in the classical FE methods are relatively small,
the FE element matrices Mi are usually sparse [12, 13]. Furthermore, after
suitable permutations (relabeling of the nodes), they usually have a banded
structure [14]. This furthermore implies that the matrices A and E of the
descriptor state-space model are usually sparse banded matrices or its partitions
are sparse banded matrices. For example, in the case of (1), the non-zero blocks
of the matrices E and A are banded matrices after suitable permutations. For
clarity and brevity, the methods proposed in this paper are developed under the
assumption that the matrices A are E are sparse banded matrices (as in the
case of the discretized heat equation described below), and as we explain later
on, the developed methods can be easily generalized to the case where the block
partitions of A and E are sparse banded matrices (as in the case of (1)).
To illustrate the structure of the FE matrices, we generate a complex geometry 2D domain shown in Fig. 1(a). Such a domain can for example represent
a cross-section of a building or of an arbitrary mechanical element. We assume
that the heat actuators are distributed across the domain, and that we are
able to observe the temperature change at predefined points. We assume that
2

the temperature dynamics of such an object is described by the heat equation
with Dirichlet boundary conditions and we introduce the FE mesh shown in
Fig. 1(a). We use tent-shaped linear basis functions. The pattern of the matrix
A (the matrix E has a similar pattern) of the descriptor FE state-space model
is shown in Fig. 1(b). In Fig. 1(c) we show the pattern of the matrix A, after
applying the reverse Cuthill-McKee permutation algorithm [14]. Similar banded
matrix structures appear when the domain is three-dimensional, as well as in
the models resulting from the FD methods [15].
(a)

(b)

1

1

56

112

168

1

56

112

168

56
112

(c)

168
1
56
112
168

(d)

Figure 1: (a) A 2D domain and a FE mesh for the discretization of the heat equation. Patterns
of (b) the matrix A and (c) the matrix obtained by permuting A using the reverse CuthillMcKee algorithm. (d) A surface plot of the absolute values of the entries of the solution of
the Ricatti equation.

The classical numerical algorithms [16, 17] for solving optimal control problems have O(n3 ) computational and O(n2 ) memory complexities [18, 19, 20],
and consequently, they are computationally infeasible for large-scale systems.
Furthermore, these algorithms produce centralized control solutions which require all sensor data collected from the system to be transferred to a centralized
computing unit that calculates and sends control signals to system’s actuators. Such an approach usually fails because in the case of a large number of
actuators and sensors, an enormous amount of information needs to be transferred and processed in real-time. This implies that in the case of large-scale
systems, control has to be performed in a distributed/decentralized manner
3

requiring only local computation and communication between sensors and actuators [21, 22, 23, 24, 25]. This further implies that the matrices of the optimal
controllers need to be sparse and highly structured, mathematically reflecting
the need for distributed communication and computation.
The problem of designing sparse control matrices has been considered by
several authors, see for example [26, 27, 28] and follow-up approaches. Due to
their O(n3 ) computational and O(n2 ) memory complexities such methods are
computationally infeasible for large-scale systems. Several other methods with
lower computational and memory complexities have been developed [29, 30].
However, these methods are applicable only to a narrow class of systems, such
as systems composed of strings of subsystems [29], or decomposable systems
[30]. On the other hand, it has been observed that in the case where the number
of system’s inputs and outputs is relatively small compared to n, the solution
of the algebraic Riccati equation necessary to compute the Linear-Quadratic
(LQ) control law, can be factored into a product of low-rank matrices [31, 32,
33]. Such an LQ control law can be computed and implemented relatively
efficiently. The applicability of these methods is restricted mainly because largescale systems usually have a large number of inputs and outputs. Moreover, the
computed LQ feedback matrix is dense and thus it is difficult to implement the
resulting controller in distributed or decentralized manners.
The main computational bottleneck in solving the optimal control problems
originates from the high complexity of the numerical methods for solving the
Lyapunov and Riccati matrix equations [18, 34]. In this paper, the Lyapunov
equation defined for descriptor state-space models is referred to as the Generalized Lyapunov (GL) equation, whereas the one defined for the non-descriptor
state-space models (where E = I), is referred to as the Non-Generalized Lyapunov (NGL) equation. In the same spirit we refer to Riccati equations as
Generalized Riccati (GR) and Non-Generalized Riccati (NGR) equations. Recently, it has been shown that the solutions of the NGL equations with symmetric banded matrices belong to the class of localized matrices [35, 36]. Loosely
speaking, localized matrices are characterized by the fact that a small number
of matrix entries is significantly larger than other entries, and the dominant
entries are usually grouped in localized matrix regions [37, 38, 39, 40, 36, 41].
Furthermore, the entries are decaying away from these regions. Off-diagonally
decaying matrices are typical examples of localized matrices [42, 43, 44]. An
example of an off-diagonally decaying matrix is shown in Fig. 1(d).
The localization and spatially decaying phenomena of the solutions of distributed control problems for infinite-dimensional systems have been analyzed
in [45, 46, 47]. It has been shown that the solutions of the operator NGL and
NGR equations are spatially decaying. The operator NGL and NGR equations
are the infinite-dimensional counterparts of finite-dimensional NGL and NGR
equations [48, 49, 51]. The extension of these approaches to finite-dimensional
systems and to GL and GR equations, is not straightforward. Furthermore,
these approaches do not take into consideration the influence of the condition
numbers of the system matrices on the spatially decaying phenomenon. Finally, they do not deal with the pressing problem of efficiently computing the
4

distributed controllers. An approach for computing low-rank solutions to the
operator NGL equations has been developed in [48]. However, this approach is
not suitable for distributed implementation.
In [35], we have shown that the decay rate of the entries of the solution of the
NGL equation is faster for smaller condition numbers of A. By exploiting this
localization phenomena, we have developed two efficient methods for computing
approximate sparse banded solutions of the NGL equations. The importance
of the results of [35] lies in the fact that the Riccati equation can be solved by
solving a series of Lyapunov equations. It should be expected that the solution
of the Riccati equation is localized and that it can be approximated by a banded
matrix in a computationally efficient manner. Consequently, it might be possible
to efficiently approximate the solution of the LQ optimal control problem by a
banded matrix. However, the results and the methods presented in [35] are not
applicable either to state-space models with non-symmetric matrices, or to GL
and GR equations. On the other hand, numerical evidence indicates that the
solutions of GL and GR equations are localized. As an illustration, Fig. 1(d)
shows the surface plot of the solution of the GR equation for the discretized heat
equation and for the geometry defined in panel (a). This important observation
motivates us to develop the methods presented in this paper.
In this paper, we show that the a priori (sparsity) patterns of the exact
solutions of the GL equations for FE or FD descriptor state-space models are
banded matrices. The a priori pattern predicts the dominant non-zero entries of
the exact solution [50]. We furthermore show that for well-conditioned problems
(the notion of a well-conditioned problem is introduced in Section 3), the a
priori patterns are not only banded but also sparse matrices. On the basis of
these results, we develop two computationally efficient methods for computing
sparse approximate solutions of GL equations. In contrast to [35], the methods
presented in this paper are applicable to descriptor state-space models with
non-symmetric matrices. Using the developed methods and the inexact Newton
method [52, 53], we show that the solution of the generalized Riccati equation
can be approximated by a banded matrix. This enables us to approximate the
feedback matrices of the LQ controllers by (sparse) banded matrices and to
develop a novel computationally efficient approach for optimal control of largescale systems. We perform extensive numerical experiments to demonstrate the
effectiveness of the proposed approach. The MATLAB codes and the models
used in our simulations can be found in [54].
The main focus of this paper is on the development of the algorithms and
on the investigation of their numerical performance. In problem formulations
that differ from ours, the influence of the approximation errors (originating
from approximately solving NGL equations) on the stability and convergence
properties of the inexact Newton method, has been studied in [52, 53]. In our
future work, we will use these results to analyze the stability and convergence
properties of the inexact Newton method used in this paper.
The paper is organized as follows. In Section 2 we formulate the control
problem. In Section 3 we present a method for determining the a priori pattern
of the solution of the Lyapunov equation, which is used in Section 4 to develop
5

the computationally efficient approximation methods. In Section 5 we present
numerical results, and finally, in Section 6 we present the conclusion and discuss
the future work.
2. Problem Formulation
2.1. Notation and Preliminaries
We define the notation used in this paper. The notation X = [xi,j ] denotes
a matrix whose (i, j) entry is xi,j . The exponential of the matrix X is denoted by exp (X). The Frobenius and 2-norms are denoted by kXkF and kXk2 ,
respectively. The maximal and minimal eigenvalues are denoted by λmax (X)
and λmin (X), respectively. Similarly, the maximal and minimal singular values
are denoted by σmax (X) and σmin (X), respectively. The condition number is
denoted by κ (X). Let Q = [qi,j ] ∈ Rn×n be an arbitrary matrix. The projection operator, projecting the matrix Q onto the pattern of an arbitrary matrix
X = [xi,j ] ∈ Rn×n , is defined by

 
qi,j ,
xi,j 6= 0
VX Q =
(2)
0,
xi,j = 0
 
where VX Q ∈ Rn×n .
2.2. Control Problem Formulation
For clarity, the methods of this paper are developed assuming a concrete
large-scale optimal control problem in which the system matrices E and A have
sparse banded structures (after suitable permutations). The generalization of
the developed methods to the descriptor state-space models (1) in which the
blocks of the system matrices E and A are sparse banded matrices is explained
in Remark 4.1 in Section 4.
Despite the fact that the matrices of FE or FD models obtained for example
using COMSOL or Ansys modeling software are sparse, they are unstructured,
and as such, they might not be suitable for the design of distributed control algorithms. This motivates us to first explain a procedure for transforming unstructured FE models into banded state-space models, more suitable for the efficient
control design. Consider a 2D domain with the FE mesh shown in Fig. 1(a). We
assume that heat actuators and temperature sensors are randomly distributed
over the domain. In practice, localized heat actuators based on the resistive
heaters can be placed at predefined positions, and the thermocouples can be
used as temperature sensors. The control objective is to keep the temperature
distribution constant using the feedback information provided by the sensors.
Similar control problems appear in various real-life applications, such as for
example in control of thermally actuated deformable mirrors used in adaptive
optics applications [55, 56, 57, 58, 59]. The temperature dynamics is governed
by the 2D heat equation with constant parameters. For simplicity, we assume
that the temperature is constant and equal to zero on the boundaries (Dirichlet

6

boundary conditions), for more details see Section 5. The discretized model has
the following form: E1 ż(t) = A1 z(t), where z(t) ∈ Rn is the vector of temperatures at the discretization nodes, and E1 , A1 ∈ Rn×n are the system matrices.
The pattern of the matrix A1 is shown in Fig. 1(b) (the matrix E1 has a similar
pattern). We permute the matrices of the state equation into a banded form.
The permutation of A1 is defined by A = P A1 P T , A ∈ Rn×n , where P ∈ Rn×n is
a permutation matrix of the reverse Cuthill-McKee algorithm [14]. The banded
pattern of A is shown in Fig. 1(c). The matrix P defines a coordinate transformation x(t) = P z(t), where x(t) ∈ Rn×n is the reordered state vector. Taking
into account that permutation matrices are orthogonal (P T P = P P T = I), we
have z(t) = P T x(t), and consequently, P E1 P T ẋ(t) = P A1 P T x(t). The last
equation can be written compactly E ẋ(t) = Ax(t), where E = P E1 P T , and
E ∈ Rn×n . The pattern of E is similar to the pattern of A. We modify the
state equation by introducing the control input and add an output equation.
The resulting state-space model is
E ẋ(t) = Ax(t) + Bu(t), y(t) = Cx(t),
where u ∈ Rm is the control input vector, y ∈ Rr is the observed output
vector, B ∈ Rn×m and C ∈ Rr×n are the system matrices. Assuming that
approximately 50% of the mesh points are directly observed and controlled, the
patterns of the matrices B and C are shown in Fig. 3 in Section 5. The methods
developed in this paper are applicable to systems with the matrices B and C
having similar diagonal-like patterns.
We are interested in finding the control input vector producing a desired
temperature profile. That is, we are interested in the set point tracking problem
[16, 17]. For presentation clarity, we will use a simple set point tracking problem
formulation [17]. Let ud be the control input that brings the system to the
desired state xd and to the desired output yd in the steady state. That is, the
input ud should satisfy Axd + Bud = 0 and yd = Cxd , either exactly or in the
least-squares sense. Next, we introduce the error variables ũ(t) = u(t) − ud ,
x̃(t) = x(t) − xd , and ỹ(t) = y(t) − yd . It is easy to verify that the error
variables satisfy:
˙
E x̃(t)
= Ax̃(t) + B ũ(t),
ỹ(t) = C x̃(t).
The LQ control problem is defined by [34]:
Z ∞
min
ỹ(t)T Qỹ(t) + ũ(t)T Rũ(t)dt,
ũ(t)

(3)
(4)

(5)

0

subject to (3) and (4),
where Q ∈ Rr×r and R ∈ Rm×m are the weighting matrices. We assume
that E is nonsingular, R is positive-definite, (E, A, B) strongly stabilizable, and
(E, A, C T QC) is strongly detectable. Furthermore, for simplicity and without
the loss of generality, we assume that the matrices Q and R are diagonal. Under
7

these assumptions the solution of the LQ control problem exists, it is unique,
and it is given by the feedback law [34]:
ũ(t) = −F x̃(t), F = R−1 B T ZE,

(6)

where Z is the unique positive-definite solution of the GR equation D [Z] = 0,
where
D [Z] = C T QC + E T ZA + AT ZE − E T ZBR−1 B T ZE.

(7)

Furthermore, the (closed loop) dynamics of the system obtained with this control
˙
E x̃(t)
= Āx̃(t), Ā = A − BF,

(8)

is asymptotically stable. We solve the Riccati equation D [Z] = 0 using the inexact Newton method [52, 53]. It can be easily shown that the Fréchet derivative
of D [Z] is
0

DZ [Y ] = E T Y A + AT Y E − E T ZBR−1 B T Y E − E T Y BR−1 B T ZE.
In the k-th iteration, the inexact Newton method solves
h
i
h
i
0
DẐ
Ẑk − Ẑk−1 = −D Ẑk−1 + Sk ,
k−1

(9)

(10)

where Ẑk is an unknown matrix, and the matrix Ẑk−1 is determined in the
previous iteration k − 1. In (10), Sk is a residual matrix, quantifying the ”inexactness” of the Newton step. The convergence properties of the inexact Newton
methods have been studied in [53, 52, 60]. When Sk = 0, the equation (10) becomes the classical Newton iteration. After some transformations, the equation
(10) becomes the GL equation
E T Ẑk Āk−1 + ĀTk−1 Ẑk E = Pk + Sk ,

(11)

where
T
Pk = −C T QC − F̂k−1
RF̂k−1 , Āk−1 = A − B F̂k−1 , F̂k−1 = R−1 B T Ẑk−1 E.
(12)

If the residual term is zero, then Ẑk becomes Zk , that exactly solves
E T Zk Āk−1 + ĀTk−1 Zk E = Pk .

(13)

That is, the GR equation is actually solved by solving the GL equations for every
iteration k of the Newton method. The inexact Newton method is summarized
in Algorithm 1.
The initial guess in Algorithm 1 should be chosen as a scaled identity matrix
or as a sparse banded matrix. This algorithm ensures that the matrices Ẑk are
always banded (a method for selecting the structure and the bandwidth of the
8

approximate solution in the second step of Algorithm 1 is presented in Section 3).
Due to this fact, and taking into consideration that the matrices A and E are
banded, and that the matrices Q, R, and R−1 are diagonal, the approximation
of the feedback control matrix, denoted by F̂k , as well as the matrix Pk , are also
banded matrices. The matrices C and B, with the patterns shown in Fig. 3,
as well as the matrices C and B with similar diagonal-like patterns, do not
fundamentally change the banded structure of the matrices Ẑk , F̂k , and Pk .
Algorithm 1 Inexact Newton Algorithm
Input: The matrices E, A, B, C, Q, and R.
Output: Banded approximate solution Ẑ of the GR equation.
1. Choose an initial guess Ẑ0 .
2. For k = 1, 2, . . . , Nmax , compute a banded matrix Ẑk that approximately
solves the GL equation (up to a residual Sk in (11)):
E T Zk Āk−1 + ĀTk−1 Zk E = Pk ,
where Nmax is the maximal number of iterations. Terminate the loop if
the convergence tolerance is satisfied.

The main research question is: Under which conditions does the GL equation
admit banded approximate solutions? In addition, can the bandwidth and the
structure that capture the most dominant entries of the exact solution to the
GL equation be predicted a priori? The main challenge is to develop efficient
methods for computing banded approximate solutions to the GL equations. For
example, the numerical evidence shows that the exact solutions of the GL and
GR equations for the FE heat equation model belong to the class of localized
matrices (Fig. 1(d) shows the solution of the GR equation), implying that the
GL and GR equations admit banded approximate solutions. However, despite
the fact that similar numerical observations can be made for other types of
systems [55], in order to develop efficient numerical algorithms, we need analysis
methods that can provide us with deeper understanding of these phenomena.
In the next section we will present a method for computing the a priori
pattern of the exact solution Zk . The a priori pattern predicts the dominant
non-zero entries of the exact solution [50]. This method enables us to show that
the banded approximate solution exists due to the fact that the system matrices
are banded. For well-conditioned problems (the notion of a well-conditioned
problem is introduced in Section 3), the a priori pattern is not only banded
but also a sparse matrix, implying that the most dominant entries of the exact
solution belong to a sparse banded matrix pattern. These results are used in
Section 4 to develop computationally efficient methods for approximating the
solutions of the GL equations. The developed methods are used in the second
step of Algorithm 1.

9

3. A Priori Patterns
In this section we present a method for computing the a priori patterns of
the exact solutions to the GL equations. The a priori pattern will be used to
significantly reduce the dimensionality and the complexity of the approximation
problem and to develop the efficient approximation methods in Section 4. These
reductions will be accomplished by eliminating or neglecting the entries of the
exact solution that do not belong to the a priori pattern.
In the analysis presented here, we are mainly interested in the non-zero
structure of the a priori pattern, but not in the exact numerical values of its
entries. Furthermore, we ignore possible zero entries created by coincidence in
the numerical values of the matrix entries. This is a standard assumption that
is equivalent to the assumption that all the entries of the matrices used for the a
priori pattern analysis are positive, see for example [50] and references therein
(see also Remark 3.2).
The fact that the numerical evidence shows that the matrix Zk belongs to
the class of localized matrices, might lead to the conclusion that we can simply
assume a banded a priori pattern of Zk , with a completely dense bandwidth
region, and proceed with the development of the approximation algorithms.
However, the pattern bandwidth is difficult to guess. Since we are dealing with
large-scale problems, a priori patterns with completely dense bandwidth regions
and large bandwidths will not necessarily produce significant computational savings. In many cases, a sparser banded matrix would be a much better choice
that could lead to significant computational savings. In [35] we have developed
a method for predicting the a priori patterns of approximate solutions of NGL
equations. Here, we modify and generalize such a method to GL equations.
This method is summarized in Algorithm 2.
Applying the vec (·) operator to (13), we obtain
M z = p,
where M ∈ Rn

2

×n2

2

(14)

2

, z ∈ Rn , and p ∈ Rn are defined as follows
M = ĀTk−1 ⊗ E T + E T ⊗ ĀTk−1 ,

(15)

z = vec (Zk ) , p = vec (Pk ) .

(16)

Due to the fact that the dimension of the system (14) is extremely large, the
matrix M cannot be formed explicitly. Moreover, it might not be possible to
compute the solution z = M −1 p by using the direct solution methods [61] or
by using preconditioned iterative techniques [62]. Despite this obstacle, the fact
that M −1 can be formally expressed as the sum of powers of M can help us gain
more insight into the structure of the a priori pattern. Similarly to [35], we use
the Newton-Schulz iteration to analyze the a priori pattern. The Newton-Schulz
method approximates M −1 using the following iteration Xl+1 = Xl (2I − M Xl ),
where l = 0, 1, 2, . . . , is the iteration index, and Xl is an approximation of
M −1 . Let the approximation error of the Newton-Schulz iteration be quantified

10

by El = I − M Xl . It can be easily shown that if X0 =

2
σmin
(M ) M T , then
2

kEl k2 ≤

κ (M ) − 1
2

κ (M ) + 1


2
2/ σmax
(M ) +

!2l
.

(17)

Using the Netwton-Schulz iteration and the aforementioned initial guess, by
back substitution it can be shown that the pattern of Xl is determined by the
pattern of the following matrix

w  T
M ,
(18)
Xl = I + M T M + . . . + M T M
where w = 2l − 1. From (17) we see that if the matrix M is well-conditioned,
then for a relatively small l we obtain a good approximation to M −1 . This is
a direct consequence of the quadratic convergence rate of the Newton-Schulz
iteration described by (17). This, together with (18), implies that for a wellconditioned M , we can get a relatively good guess of the a priori pattern of
M −1 by summing only a few powers of M T M and by multiplying the resulting
sum with M T . Since the matrices A and E are sparse banded, the matrix M
is sparse (multi) banded. All this implies that for a well-conditioned M , the
a priori pattern that captures the most dominant entries of M −1 is a sparse
(multi) banded matrix. In the general case, the fill-in of the matrix Xl will
increase with l or consequently, with w. However, as long as w  n, the matrix
Xl will be a sparse (multi) banded matrix.
In practice, the matrix M cannot be formed explicitly due to its large-scale
nature. Consequently, we cannot compute the a priori patterns of the vector z
and the matrix Zk . This problem can be resolved by reversing the vec operation
used to form (14). Since the approximate pattern of M −1 is determined by Xl ,
from (14) and (18), we have
w
z = g + M T M g + . . . + M T M g,
(19)
where z is a vector whose non-zero entries determine an a priori pattern of
z, and g = M T p. From (19), we can reconstruct the a priori pattern of Zk .
Namely, by reversing the action of the vec (·) operator on (19), we obtain
Z = G1 + G2 + . . . + Gw+1 ,

(20)

where Gi , i = 1, 2, . . . , w + 1, is the matrix corresponding to (M T M )i−1 g term
in (19), and the matrix Z is the a priori pattern of Zk . The matrix terms Gi
can be computed recursively as follows. Let us start with the first entry in the
sum g = M T p. From the Kronecker sum structure of M given in (15), we can
conclude that the matrix form of this vector term, denoted by G1 , becomes
G1 = EPk ĀTk−1 + Āk−1 Pk E T .
11

(21)

By using the basic properties of the Kronecker product, it can be shown that
the vector M T M g has the following matrix form:


G2 = E E T G1 Āk−1 + ĀTk−1 G1 E ĀTk−1 + Āk−1 E T G1 Āk−1 + ĀTk−1 G1 E E T .
(22)
Using a similar reasoning, it can be shown that the matrix Gi+1 corresponding
i
to the vector M T M g can be computed recursively (see Remark 3.2)


Gi+1 = E E T Gi Āk−1 + ĀTk−1 Gi E ĀTk−1 + Āk−1 E T Gi Āk−1 + ĀTk−1 Gi E E T .
(23)
Our numerical experience shows that in order to get a good estimate of the a
priori pattern, an identity matrix has to be added to (20), resulting in
Z w = I + G1 + G2 + . . . + Gw+1 .

(24)

Using (24) we can predict the a priori pattern of the matrix Zk . The parameter w is a user choice. The procedure for computing the a priori pattern in
summarized in Algorithm 2.
Algorithm 2 Computation of the a priori pattern Z w
Input: the matrices Āk−1 , E and Pk and a relatively small positive integer w
(see Remark 3.1 and Remark 3.2).
Output: a priori pattern Z w .
1. Compute G1 = EPk ĀTk−1 + Āk−1 Pk E T .
2. For every i = 1, 2, . . . , w, compute recursively

Gi+1 = E E T Gi Āk−1 + ĀTk−1 Gi E ĀTk−1

+ Āk−1 E T Gi Āk−1 + ĀTk−1 Gi E E T .
3. Compute Z w = I + G1 + G2 + . . . + Gw+1 .
For a relatively small value of w (w  n), the a priori pattern can be
computed very quickly. Namely, since the matrices A and E are sparse banded,
the a priori pattern can be computed with the O(n) computational and memory
complexities. Consequently, our method can be effectively used for large-scale
problems. Due to the fact that the matrices A, E, and Pk are banded, the matrix
Z w is also banded (for w  n), however, its fill-in increases by increasing w, see
also Remark 3.1. On the other hand, we have shown that for a well-conditioned
M , w is generally small, implying that the a priori pattern of Zk is a sparse
banded matrix. That is, for a well-conditioned M , the dominant entries of
the exact solution Zk belong to the sparse banded pattern of the matrix Z w .
Furthermore, as the condition number of M increases, according to the NewtonSchultz error analysis, the parameter w needs to be increased in order to ensure
12

that the a priori pattern captures the dominant entries of the exact solution. This
implies that the fill-in of the a priori pattern Z w increases with the condition
number of M .
The approximation problems for which the matrix M is well-conditioned are
referred to as the well-conditioned problems. Obviously, the condition number
of M is closely related to the spectrum or singular values of the matrices Āk−1
and E. However, to the best of our knowledge, establishing such a relation is a
difficult problem. In the case of non-descriptor state-space models (E = I) with
symmetric Āk = A, it is easy to show that the condition number of M is equal
to the condition number of A [35]. For small-size systems, we can explicitly
form the matrix M and compute its condition number.
Remark 3.1. Although the parameter w is fixed and set by the user, the fill-in
of the a priori pattern will increase with the increase of the iteration k of the
inexact Newton method, summarized in Algorithm 1. This is due to the fact
that the matrix Āk−1 depends on the previous solution of the GL equation Ẑk−1 ,
and in every step k we need to recompute the a priori pattern. After a certain
number of iterations k, the matrix Āk−1 of the inexact Newton method becomes
dense, implying a dense a priori pattern. To prevent this, the a priori pattern
should be computed only for small values of k, and then such an a priori pattern
should be fixed and kept constant in further iterations of Algorithm 1.
Remark 3.2. As stated at the beginning of this section, we are mainly interested in the non-zero structure of the a priori pattern, but not in the exact
numerical values of its entries. Although the probability of zeros created by coincidence in the numerical values of the matrix entries is very small, it is still
a good practice to provide a safeguard against this effect. Furthermore, since
in the derivation of the a priori pattern, we ignore the scalar coefficients of
the Newton-Schulz iteration, the entries of the matrices Gi can rapidly grow or
less likely, decrease with i, creating numerical instabilities or zeros produced by
numerical underflow. Good safeguards against these effects are to operate with
the matrices obtained by taking the absolute values (taken element-wise) of the
matrices Āk−1 , E, and Pk , and after every iteration i, to transform the matrices
Gi into binary matrices according to their non-zero patterns. That is, the (p, q)
entry of the binary matrix is equal to 1 if the corresponding (p, q) entry of Gi is
non-zero, otherwise, it is zero.
Next, we show the numerical performance of the pattern prediction method.
Consider the FE model of the heat equation described in Introduction and
illustrated in Fig. 1 (for more details see Section 5). We assume that w = 1 in
(24), and that the matrices Āk−1 and Pk are computed for Ẑk−1 = I. Figure
2(a) shows the predicted pattern of Ẑk whereas Fig. 2(b) shows the absolute
values of the entries of the exact solution Zk . The exact solution is computed
and visually represented using the lyap (·) and imagesc (·) MATLAB functions.
From Fig. 2 we see that the a priori pattern is able to accurately predict the
most dominant entries of the exact solution. It is of interest to compute the
relevant condition numbers, since according to our theory, they determine the
13

localization degree of the exact solution and the density of the a priori pattern.
The condition numbers of E T and ĀTk−1 are, respectively, 142 and 18, whereas
the condition number of M is 166. The matrix M is well-conditioned, implying
that the most dominant entries of the exact solution are grouped around the
main diagonal, which is numerically confirmed in Fig. 2.
(a)
1

1

42

84

126

148

(b)
1

42

42

84

84

126

126

168

168

1

42

84

126

148

5
-5
-15
-25
-35

Figure 2: The prediction of the pattern of the exact solution of the Lyapunov equation. (a)
The pattern predicted by Algorithm 2. (b) The absolute values of the entries of the exact
solution in the logarithmic scale. The model is obtained using the FE discretization of the
heat equation for the domain shown in Fig. 1.

4. Methods for Approximating the Solution of the Lyapunov Equation
On the basis of the derived a priori pattern, in this section we present two
methods for approximating the solution of the GL equation. The first method,
based on the solution of the reduced normal equations, is a simple generalization
of the method presented in [35]. For completeness, we briefly summarize this
approach in Method 1. However, although the computational complexity of
this method is quite low, and for well-conditioned problems it scales almost
with O(n), the memory complexity is still significant, regardless of the fact that
it still lower than the O(n2 ) memory complexity of the MATLAB lyap (·) solver.
This motivates us to develop a second method based on the Faber expansion of
the matrix exponential and on the gradient projection method. This method
is a generalization of the Chebyshev expansion method presented in [35], which
is only applicable to symmetric matrices and NGL equations. The memory
complexity of the second method is significantly lower than the first one (for
very sparse problems it is O(n)), however, the computational complexity is a
bit higher, for more details, see Section 5. The second approach is summarized
in Method 2 at the end of this section.
4.1. First Approximation Method
Consider the system of equations (14). Using the a priori pattern of Zk
computed using Algorithm 2, we can significantly reduce the size of this system
and efficiently solve it in the least-squares sense. Let us define zw = vec (Z w ).
The zero entries of the vector zw determine a zero-entry pattern of the vector
14

that approximates the exact solution z in (14). The dimension of the vector z
should be reduced by eliminating these zero entries. Let the reduced vector be
denoted by z1 ∈ Rn1 . Consequently, the corresponding columns of the matrix
M should be eliminated. More precisely, if the ith entry of zw is zero, then the
ith column of the matrix M should be eliminated. This column reduction can
also produce zero rows of the matrix M , and we also assume that such rows
are eliminated. Let the reduced matrix be denoted by M1 ∈ Rn2 ×n1 , where
n2 ≤ n2 (usually n2 < n2 ), and n1  n2 , and let the vector p after possible
row reductions be denoted by p1 ∈ Rn2 . Once the system of equations has
been reduced, the approximate solution can be found by solving the following
least-squares problem
2

min kp1 − M1 z1 k2 .
z1

(25)

The solution, denoted by ẑ1 , can be found by solving the system of normal
equations
M1T M1 ẑ1 = M1T p1 ,

(26)

using for example the Conjugate Gradient Least Squares (CGLS) method [63].
The computational and memory complexities of the first method are thoroughly
analyzed and numerically simulated in [35], and for brevity we omit this analysis. In brief, one iteration of the CGLS costs approximately nnz (M1 )+3n1 +2n2
flops, where nnz (·) denotes the number of non-zero entries determined by the
parameter w. If M is well-conditioned, w can be chosen as a small number,
and the solution can be computed with (almost) O(n) computational complexity. However, special attention needs to be given to the implementation of the
first method. Namely, a naive implementation based on a direct formation of
the matrix M using for example, the MATLAB kron (·) function, followed by a
reduction of this matrix, is not possible for a large n due to the high memory
demand. Instead, the non-zero entries of the reduced matrix M1 should be calculated without forming the matrix M . This itself requires some computational
time, however, our numerical experience shows that with proper implementation
in the MATLAB sparse matrix toolbox, this time is shorter than the time necessary to solve the least-squares problem. Also, although for well-conditioned
problems the memory complexity scales almost linearly with n, the reduced matrix M1 for a large n still occupies significant memory space, see Section 5 for
more details.
Method 1 Approximation of the GL equation using the reduced least-squares
problem
1. Compute the a priori pattern Z w using Algorithm 2.
2. On the basis of the non-zero structure of Z w , form the reduced matrix M1
and solve the least-squares problem (25) using the CGLS.

15

4.2. Second Approximation Method
The approximate solution of the GL equation (13) can be found by solving
the following optimization problem:
2
F

min Pk − E T Zk Āk−1 − ĀTk−1 Zk E
Zk

, subject to Zk ∈ ΩZ w ,

(27)

where ΩZ w is a set of all matrices with the pattern defined by the pattern of
the matrix Z w computed by Algorithm 2. The solution of (27) can be obtained
using the gradient projection method [64]:


Z̃ki+1 = VZ w Z̃ki − δ i Nki ,
(28)
where Z̃ki is the approximate solution at the i-th iteration, VZ w (·) is the projection operator defined in (2) (this operator projects the argument matrix on the
pattern of Z w ), Nki is the gradient matrix, and δ i is the step size. Expressing
the cost function (27) as a sum of matrix traces, and by differentiating such an
expression, it can be shown that
Nki = −2ERki ĀTk−1 − 2Āk−1 Rki E T , Rki = Pk − E T Z̃ki Āk−1 − ĀTk−1 Z̃ki E. (29)
To define the step size δ i , we first introduce
 
J Z̃ki = Pk − E T Z̃ki Āk−1 − ĀTk−1 Z̃ki E

2
F



 
, Z̃ki δ = VZ w Z̃ki − δNki .

(30)

Then, the step size is determined by the Armijo rule along the projection arc,
defined by δ i = ζ gi δ, where δ > 0, ζ ∈ (0, 1), and gi equals the first nonnegative
integer g satisfying
h  i

 
T
 
J Z̃ki − J Z̃ki ζ g δ ≥ σvec Nki vec Z̃ki − Z̃ki ζ g δ ,
(31)
and where σ > 0. The convergence of the gradient projection method can be
analyzed by recognizing that the cost function in (27) can be written in the
vectorized form, which corresponds to the least squares formulation of the solution to the linear system of equations (14). From there it follows that the
convergence rate is determined by the minimal and maximal singular values of
M , that is, it is determined by the condition number of M [64]. Provided that
the parameter w is much smaller than n, one step of the gradient projection
method can be implemented with O(n) computational and memory complexities. In the case of ill-conditioned problems, the convergence of the gradient
projection method can be slow. Consequently, for fast convergence, it is essential that it be initialized with a good initial guess of the solution. In the sequel
we present a method for generating such an initial guess. By first multiplying
−1
(13) from the left with E T
and then from the right with E −1 , we transform
the GL equation into the NGL equation:
Zk Āk−1 E −1 + E T

−1

ĀTk−1 Zk = E T
16

−1

Pk E −1 .

(32)

Assuming that Āk−1 E −1 is a stable matrix, the unique solution of (32) can be
expressed as follows (see Theorem 13.23 in [65]):
Z ∞

(33)
X =−
exp (tA) P exp tAT dt,
0

T
T
where A = E −1 ĀTk−1 , P = E −1 Pk E −1 , and where we rely upon the fact
T
−1
that E −1 = E T
. By inverting the matrix E, the sparsity of the problem
has been lost. To preserve the sparsity, we search for a sparse approximate
inverse of E. Such an inverse exists due to the following facts. Given that the
matrix E is banded, under some mild conditions its inverse will exhibit some
form of localization, or in the positive-definite case it will belong to the class
of off-diagonally decaying matrices, see for example [43, 38, 66, 37, 67]. In
the general case, when the matrix E is not strictly banded or positive-definite,
the sparse approximation of E −1 still exists if the matrix E is relatively wellconditioned. This follows from the convergence analysis of the Newton-Schulz
method, see Section 3. All these facts motivate us to search for a sparse banded
approximate inverse of E. This approximate inverse can be computed by solving
the following optimization problem [68, 69]
2

min kI − EEkF , subject to E ∈ ΩE .
E

(34)

where E is a sparse approximate inverse of E and ΩE is a set of matrices with
a given pattern. The problem (34) is highly parallelizable and can be split into
n independent low-dimensional least-squares problems that can be efficiently
solved using the SPAI solver [68, 69, 70, 71, 72]. In some situations, it might
be more appropriate to replace the term I − EE by I − EE in (34). In practice,
both forms can be easily computed and the most accurate solution should be
selected. The choice of the set ΩE depends on the pattern of E. Using arguments
analogous to those in Section 3, it can be shown that a good estimate of the
pattern of E is given by I + E + E 2 + . . . + E k1 , where k1 is a small positive
integer. Similarly to the analysis developed in Section 3, it can be shown that k1
is generally small for a well-conditioned E. By substituting in (33), the matrix
E −1 by the matrix E, that solves (34), we define the following approximation
Z ∞

X ≈ X1 , X1 = −
exp (tA1 ) P1 exp tAT1 dt,
(35)
0

where A1 = E T ĀTk−1 and P1 = E T Pk E. Because the matrices E, Āk−1 , and Pk
are banded, the matrices A1 and P1 are also banded.
Our goal is to approximate X1 by a banded matrix. Denote with λRS , λRL ,
and λIL the smallest real part, the largest real part, and the largest absolute
value of the imaginary part of eigenvalues of A1 . Because the matrix A1 is sparse
banded, its extreme eigenvalues can be estimated with O(n) complexity using
the ARPACK software (for more details see Section 4 of [73]). From Theorem

17

4.3 in [74], we have that the integral (35) can be approximated by:
X2 = −

q
X



ψωj exp t̃j A1 P1 exp t̃j AT1 ,

(36)

j=−q

where t̃j = ψtj , q is a sufficiently large integer, and
−1/2


3
,
, ωj = q + q exp −2jq −1/2
2|λRL |


 q


tj = log exp jq −1/2 + 1 + exp 2jq −1/2 .

ψ=

(37)

In [74] it has been shown that the approximation error
exponentially decreases
√ 
√
with q, that is, kX1 − X2 k2 ≤ KA1 kP1 k2 exp − q , where the constant KA1
depends on λRS , λRL and λIL .
To compute X2 , we need to approximate the matrix exponential by a banded
matrix. Consequently, in the sequel we develop an algorithm for approximating
the matrix exponential by a banded matrix. The developed approach is summarized in Algorithm 3. First of all, the justification of the existence of the banded
approximation of the matrix exponential follows from the following facts. If the
matrix t̃j A1 is diagonalizable for each t̃j ∈ [0, ∞], then its matrix exponential
is an off-diagonally decaying matrix [38]. In [41] these results are extended to
non-Hermitian, not necessarily diagonalizable matrices. Furthermore, our numerical results indicate that the matrix exponentials of a wide class of sparse
banded matrices belong to the class of localized matrices.
We use the Faber polynomials to approximate the matrix exponential [73].
Denote with λ̃RS,j , λ̃RL,j , and λ̃IL,j the smallest real part, the largest real part,
and the largest imaginary part of the eigenvalues of t̃j A1 (these scalars can be
computed directly from λRS , λRL , and λIL ). Next, define [73]:
r
q 2

 23
2
2
3
3
3
2
c
c
+
λ̃
+
c
+ λ̃2IL,j
λ̃
1,j
1,j
1,j
IL,j
IL,j
λ̃RL,j − λ̃RS,j
c1,j =
, c2,j =
,
2
2
 2
 4

2
4
λ̃RL,j + λ̃RS,j
3
3
3
3
c3,j = c1,j
+ λ̃IL,j
c1,j
− λ̃IL,j
, c4,j =
.
(38)
2
The matrix A2,j is defined as follows [73]:
A2,j = √


1
t̃j A1 − c4,j I .
c3,j

(39)


The truncated Faber series approximation of exp t̃j A1 is defined by [73]:
Kp,j =

p
X
l=0

18

al,j Fl,j ,

(40)

where p is a positive integer, Faber coefficients al,j are defined in Appendix, and
the Faber polynomials Fl,j are defined by:
√
F0,j = I, Fl,j = 2

c3,j
2c2,j

l
Tl,j (A2,j ) , l > 0,

(41)

and where Tl,j (A2,j ) are the Chebyshev polynomials of the first kind defined
by the recurrence relation:
T0,j = I, T1,j = A2,j , Tk+1,j = 2A2,j Tk,j − Tk−1,j , k = 1, 2, . . .

(42)

An upper bound on the approximation error of the truncated Faber series is
given in Section 3 of [73] and for brevity we omit it. Briefly speaking, the
approximation error is smaller if p is lager, however, larger values of p decrease
the number of non-zero elements of Kp,j . This is due to the fact that the
recurrence (42) implies that the matrix Kp,j in (40) can be expressed as the sum
of the first p powers of A2,j . If p is much smaller than n, then Kp is a sparse
banded matrix and it can be computed with O(n) complexity [38]. However,
if this condition is not satisfied, then the matrix Kp,j is no longer sparse and
consequently, it cannot be computed with O(n) complexity. To ensure that Kp,j
remains sparse for larger values of p, and to ensure that it can be computed with
O(n) complexity, we sparsify the recurrence relation (42) as follows:


T̃0,j = I, T̃1,j = A2,j , T̃k+1,j = VA2 2A2,j T̃k,j − T̃k−1,j , k = 1, 2, . . .
(43)
 
where VA2 · is a projection operator constructed for example by summing the
identity matrix with the first few powers of A2,j (or A1 ). In this way, we also
sparsify the Faber polynomials. Let the sparsified approximation of the matrix
exponential, obtained by using the sparsified Chebyshev polynomials (43), be
denoted by K̃p,j
K̃p,j =

p
X

al,j F̃l,j ,

(44)

l=0

where F̃l,j are the sparsified Faber polynomials computed by substituting (43)
in (41).

19

Algorithm 3 Computation
of the sparse banded approximation of the matrix

exponential exp t̃j A1 in (36)
Input: the matrices E, Āk−1 , parameter t̃j , and the approximation order p in
(40).
Output: the approximation of the matrix exponential K̃p,j .
1. Computation of the sparse approximate inverse E of the matrix E.
(a) Compute the matrix I +E +E 2 +. . .+E k1 that represents an a priori
pattern of E −1 . The small positive integer k1 is a user choice, and it
should be selected such that the a priory pattern is sparse (k1  n).
(b) For the computed a priori pattern, compute the approximate inverse
E by solving (34) using the SPAI solver.
2. Compute the matrix A1 = E T ĀTk−1 .
3. Compute λ̃RS,j , λ̃RL,j , and λ̃IL,j that represent the smallest real part, the
largest real part, and the largest imaginary part of the eigenvalues of t̃j A1 .
These scalars can be computed from the extreme parts of the eigenvalues
of A1 .
4. Compute the constants c1,j , c2,j , c3,j and c4,j defined in (38) and compute
the matrix A2,j defined in (39).
5. Compute the coefficients of the truncated Faber series expansion (the coefficients are defined in Appendix).
6. Computation of the sparsified Faber polynomials F̃l,j .
 
(a) Compute the pattern of the projection operator VA2 · of the sparsi2
fied Chebyshev polynomials in (43). The pattern is I +A2,j +. . . Ak2,j
,
where k2 is a small positive integer (k2  n) that should be chosen
such that the resulting matrix remains sparse.
(b) Compute the sparsified Chebyshev polynomials (43). Substitute such
polynomials in (41) and compute the sparsified Faber polynomials.
7. Compute the approximation of the matrix exponential (44).

Substituting K̃p,j in (36), we obtain:
X3 = −

q
X


T
ψωj K̃p,j P1 K̃p,j .

(45)

j=−q

The matrix (45) is used as an initial guess of the gradient projection method
(28). The initial guess (45) can be computed with O(n) computational and
memory complexities provided that q  n, p in (40) satisfies p  n, and that
the projection operator in (43) is chosen such that the projected matrix remains
sparse. In practice, we do not need a high accuracy of the initial guess, and
consequently, some inaccuracies can be tolerated for the benefit of computing
it with low computational complexity. The second approximation method is
summarized in Method 2.

20

Method 2 Approximation of the GL equation using the gradient projection
method and the sparsified Chebyshev polynomials
1. Compute the a priori pattern Z w using Algorithm 2.
2. Generate an initial guess for the gradient projection method.
(a) Choose a relatively small integer q in (45) (q  n). Compute P1 =
E T Pk E. Compute the parameters t̃j , ψ, and ωj using (37).
(b) For j = −q, −q+1, . . . , q, compute the matrix K̃p,j using Algorithm 3.
Compute the initial guess (45).
3. Set Z̃k0 = X3 , and propagate the gradient projection iteration (28) until
convergence or until the predefined maximal number of iterations has been
reached.
Remark 4.1. We have developed the approximation methods assuming sparse
banded structure of the system matrices E and A. Here we briefly explain the
generalization of the proposed methods to the case of model (1) where the matrices E and A are composed of sparse banded blocks. The simplest choice of
the initial guess in Algorithm 1 is a 2 × 2 block matrix (blocks corresponding to
the blocks of the matrices E and A in (1)), in which every block is an identity
matrix. Assuming such an initial guess, the a priori pattern computed using
Algorithm 2, is a 2 × 2 block matrix with (sparse) banded blocks. This is because
the pattern is computed by multiplying the matrices E, A, and Pk , whose blocks
are (sparse) banded matrices. The fill-in of these blocks will be increased after
every multiplication, however, as before, we assume that only a relatively small
number of multiplications is necessary, and consequently, the blocks will remain
sparse. Once this pattern has been computed, Method 1 and Method 2 can be
directly applied.
5. Numerical Experiments
The simulations are performed on a desktop computer with 16GB RAM and
Intel Xeon Processor E3-1245 v5. We use the FE model of the discretized heat
equation described in Introduction. The MATLAB codes and the models used
in our simulations can be found in [54]. The discretization domain is shown
in Fig. 1(a). We generate 4 models with an increasing number of nodes. As
the node number increases, the mesh becomes denser, and the state dimension
increases. We use the models having n = 168, 841, 3687 and 15379 states. The
matrix B ∈ Rn×bn/2c is constructed by randomly placing the localized actuators
on bn/2c nodes in the 2D domain in Fig. 1(a). Similarly, we construct the matrix
C ∈ Rbn/2c×n . The patterns of the matrices B and C are shown in Fig. 3. For
simplicity, the weighting matrices Q and R in (5) are chosen as Q = I and
R = I, where I is an identity matrix. The initial guess for the solution of the
Riccati equation in Algorithm 1 is Ẑ0 = 10I. This is a reasonable initial guess,
because it is expected that the exact solution is a localized matrix with the
localization region around the main diagonal.
21

(a)

1

1 21 42 63 84

1

(b)

21
56

42
63

112

84
1

56

112

168

168

Figure 3: The patterns of (a) the matrix B and (b) the matrix C.

The exact solutions of the Lyapunov and Riccati equations are computed
using the built-in MATLAB functions lyap (·) and dare (·), respectively. Due
to the O(n3 ) computational and O(n2 ) memory complexities of the algorithms
implemented in these functions, we are only able to compute the exact solutions
for the models for which n ≤ 3687. The accuracy of the proposed methods
for solving the Lyapunov equations is quantified by comparing it to the exact
solution
ek = Ẑk − Zk

F

/ kZk kF .

(46)

The convergence of the inexact Newton method is quantified using the Frobenius
norm of the residual
h i
vk = D Ẑk
,
(47)
F

where the Riccati residual D [·] is defined in (7). First, we show the accuracy
of the approximate solution of the GL equation as a function of the number
of nonzero entries of the a priori pattern Z w . The results of solving the GL
equation in the first iteration of the inexact Newton method (Algorithm 1) are
shown in Fig. 4(a), for w = 0, 1, 2, 3. The GL equation is solved using the first
method. The residual convergence tolerance of the CGLS method is 10−7 . The
a priori patterns are shown in Fig. 4(d). Figure 4(b) shows the convergence of
the inexact Newton method for solving the GR equation for w = 0, 1, 2. It can
be observed that the convergence is relatively stable and that the steady-state
residual decreases as the number of the nonzero entries of the a priori pattern
increases. This is a direct consequence of the fact that the accuracy of solving
the GL equation, shown in Fig. 4(a), increases with an increase of the number
of nonzero entries of the a priori patterns. Figure 4(c) shows the exact solution
of the GR equation. This figure confirms that the exact solution belongs to
the class of localized matrices and confirms that the a priori patterns shown in
Fig. 4(d) are able to accurately capture the most dominant entries of the exact
solution.

22

%

Figure 4: (a) The accuracy of solving the GL equation, measured by ek defined in (46), using
the first method for several a priori patterns. (b) The convergence of the inexact Newton
method, measured by vk defined in (47), for w = 0, 1, 2. (c) An image of the exact solution of
the GR equation. The colors of the pixels correspond to the logarithms of the absolute values
of matrix entries. (d)A priori patterns for several values of w.

(a) 1
1

210

420

630

841
-5

105

(b)

-15
w=0

210

-25

315

-35

420

-45

(c)
w=2

Figure 5: (a) The absolute values of the entries of the exact feedback matrix F in the logarithmic scale. The patterns of the computed feedback matrix for the a priori patterns determined
by (b) w = 0 and (c) w=2.

Figure 5(a) shows the feedback control matrix F computed on the basis of
the exact solution of the GR equation. Figures 5(b) and (c) show the patterns
of the approximate feedback control matrices F̂ computed for w = 0 and w = 2,
respectively. The percentages of non-zero entries of F̂ , for w = 0, 1, 2 and 3 are,
respectively, 5.7, 13.0, 21.4 and 31.3. It can be observed that the feedback control
matrices computed for w = 0 and w = 2 are able to accurately capture the most
dominant entries of the exact feedback matrix. We compare the performance
of the closed loop system controlled using the exact feedback matrix and the
approximate feedback matrix calculated using the first method. The results
23

are shown in Fig. 6. It can be observed that the loss of performance of the
approximate control law is almost negligible compared to the exact control law.
Notice that the presented results are generated for w = 0, implying that for
w > 0, performance will be even better.
100

||~xi ||2

w=0
true

10-10
10-20

0

20 40 60 80 100

i

Figure 6: The closed loop performance of the system for the exact feedback matrix (denoted
by ”true”) and for the approximate feedback matrix computed using the first method for
w = 0. The integer i is the discrete-time index used in the simulation of the continuous time
dynamics.

Next, we test the performance of the second method. We investigate the
convergence of the initial guess in Method 2 as a function of the parameter q in
(45). We solve the optimization problem (34) to approximate the matrix E −1 .
The a priori pattern of E −1 is chosen as I + E + E 2 + E 3 . Such an a priori
pattern has 3.1% of non-zero entries. The approximation error is kI − EEkF =
0.73. This is a relatively crude approximation, and consequently, we are not
expecting a high overall accuracy of the initial guess. Next, we compute the
matrix A1 . For the sparsified Chebyshev polynomials we use an a priori pattern
Q1 = I + A1 . The percentage of nonzero entries of Q1 is 4.7%. The order of the
Faber expansion is p = 30. The fill-in of the final initial guess X3 , defined by
(45), is 29%. We measure the accuracy of the initial guess by (46), where Ẑk is
substituted by X3 . The accuracy of the initial guess is shown in Fig. 7(a) as a
function of the parameter q. We initialize the gradient projection method with
such an initial guess. The convergence of the gradient projection method for
several a priori patterns is shown in Fig. 7(b). The final approximation errors
for w = 0, 1, 2 are, respectively, e = 0.32, 0.26, 0.25.

24

(a)

(b)

103

1

0.8

w=0
w=1
w=2

Ji

eq

102

0.6

101

0.4
0

5

10 15 20 25 30

0

q

5000

i

10000

Figure 7: (a) The accuracy of approximating the solution of the Lyapunov equation using
(45), as a function of q. (b) The convergence of the gradient projection method for several a
priori patterns. The residual J is defined in (30).

Finally, we provide some insights into the computational and memory complexities of the proposed methods. We compare the complexity of the developed
methods to the complexity of the built-in MATLAB solver lyap (·). The computational and memory complexities of the proposed solvers for the GL equation
are shown in Fig. 8 for n = 168, 841, 3687, 15379. The exact solution obtained
using the lyap (·) function can only be computed for n ≤ 3687. For n = 15379,
the lyap (·) function cannot be used because the computer runs out of RAM,
and even if the computer had more RAM, it would take a significant amount
of time to compute such a solution. For the first method, the residual convergence tolerance of the CGLS method is 10−5 . We choose w = 1. For these
parameters, the accuracies of the converged solutions for n = 168, 841, 3687 are,
respectively, 4.4·10−4 , 0.02, 0.28, whereas the percentages of the non-zero entries
for n = 168, 841, 3687, 15379 are, respectively 24.5, 11.3, 3.9, 1.2. The accuracy
for n = 3687 degrades due to the fact that the a priori pattern generated for
w = 1 cannot capture the dominant entries as accurately as in the cases of
n = 168 and n = 841. This is most likely due to the increased condition number of the matrix M (which can only be computed for n = 168). To improve
the accuracy, a priori patterns w ≥ 1 need to be used. This will increase the
complexity shown in Fig. 8. However, this increase will not be significant for
w = 2 and w = 3, and the general growth trend of the complexity shown in the
figure will be preserved (for larger values of n).
10

4

10

2

10

0

10

10

9

10

-2

10

8

first
second
lyap

bytes

seconds

10

10

7

10

6

10

5

2

10

3

n

10

10

4

10

2

10

3

n

10

4

Figure 8: (a) The computational and (b) memory complexities of the developed methods and
the comparison with the built-in MATLAB solver lyap (·).

25

To generate the initial guess for the second method, we use the following
parameters p = 20 and q = 40. The gradient projection method is computed for w = 1 and computations are stopped when the iteration number
exceeds 4000. For such parameters the accuracies of the converged solutions for
n = 168, 841, 3687 are, respectively, e = 0.16, 0.15, 0.8. The accuracy of the final
solution can be improved by increasing both p and q together with the number
of iterations of the gradient projection method. In our future research, we will
investigate the optimal selection of these parameters that should produce high
accuracy, while at the same time preserving the computational efficiency. From
Fig. 8 we can observe that both developed approaches have lower computational
and memory complexities than the MATLAB lyap(·) function whose computational and memory complexities are O(n3 ) and O(n2 ), respectively. It can be
observed that the memory complexity of the second approximation method is in
the order of magnitude smaller than the memory complexity of the first method.
This is due to the fact that in the second method all operations are performed on
n×n sparse matrices. On the other hand, the first method has a lower computational complexity. Since it operates directly on sparse n×n matrices, the second
method can be used for problems whose dimensions are in the order of n = 105
(on a computer similar to ours that has at least 16GB RAM). It should be noted
that in the first method we are not using the preconditioning techniques. The
preconditioning techniques can significantly reduce the condition number of M1 ,
and consequently, they can significantly reduce the computational complexity
of the first method. Similarly, the convergence rate of the projected gradient method can be increased by using scaled gradient methods [64]. Finally,
the MATLAB implementation of our methods is relatively slow, and significant
computational savings can be achieved by implementing the methods in the C
or C++ programming languages.
6. Conclusions and Future Work
In this paper we have presented two computationally efficient methods for
computing (sparse) banded approximate solutions of the generalized Lyapunov
equations with banded matrices. These methods, together with the inexact
Newton method for solving the generalized Riccati equations, were used to derive a novel computationally efficient approach for optimal control of finite element or finite difference state-space models. The proposed approach is able
to compute a (sparse) banded feedback matrix of the linear-quadratic optimal
controller. The future research will be oriented toward the analysis of the influence of the approximation errors on the stability and performance of the inexact
Newton method for approximately solving the generalized Riccati equation.
The proposed control approach is based on the discretized PDE models. The
finite element discretization errors depend on the type of the basis functions,
as well as on the grid density. For coarse grids, the discretization errors might
degrade the performance of the optimal control algorithm. These discretization
errors can be seen as model uncertainties and their effect can be ameliorated
using the robust control techniques, see for example [1]. The future research
26

directions are to investigate the influence of the discretization errors on the control performance and to develop a method for selecting the weighting matrices
(the matrices Q and R in the cost function (5)), such that the robustness of the
method with respect to these errors is improved.
While preparing the final manuscript, several computationally efficient methods for approximating the solutions of the non-generalized Lyapunov equations
with sparse or low-rank matrices appeared [75, 76, 77]. In our future work,
we will investigate the possibility of generalizing these approaches to descriptor Lyapunov equations. Some of the ideas on which these methods have been
founded can be used to additionally increase the computational efficiency and
numerical stability of the methods developed in this paper. In order to develop
even faster methods for generating the initial guess in Method 2, we will further
investigate the possibilities and numerical performance of the matrix function
approximation algorithms [66, 78, 79].
7. Appendix
For presentation clarity, let us denote the coefficient al,j in (40) simply by
al , and the coefficients c2,j , c3,j , and c4,j simply by c2 , c3 , and c4 , respectively.
Then, the coefficients of the truncated Faber series (40) can be approximated
by [73]:
W −1
k
1 X
gk e−il2π W , gk = exp(sk ),
al =
W
k=0

sk =








c3
k
k
c3
+ c4 + i c2 −
,
cos 2π
sin
2π
c2 +
4c2
W
4c2
W
where W is a sufficiently large positive integer and i is the imaginary unit.
Acknowledgments
This work was supported by the European Research Council (ERC) (339681)
and PSC CUNY Award A (60156-00 48).
References
References
[1] K. Zhou, J. Doyle, K. Glover, Robust and Optimal Control, Prentice Hall,
1996.
[2] P. Benner, S. Dolgov, A. Onwunta, M. Stoll, Low-rank solvers for unsteady
Stokes–Brinkman optimal control problem with random data, Computer
Methods in Applied Mechanics and Engineering 304 (2016) 26–54.
27

[3] J. Sun, A. E. Motter, Controllability transition and nonlocality in network
control, Physical Review Letters 110 (2013) 208701.
[4] S. P. Cornelius, W. L. Kath, A. E. Motter, Realistic control of network
dynamics, Nature Communications 4 (2013) 1942.
[5] A. Haber, F. Molnar, A. E. Motter, State observation and sensor selection
for nonlinear networks, in press, IEEE Transactions on Control of Network
Systems, arXiv preprint arXiv:1706.05462 (2017).
[6] H. Lee, Optimal control for quasi-Newtonian flows with defective boundary conditions, Computer Methods in Applied Mechanics and Engineering
200 (33) (2011) 2498 – 2506.
[7] Y. Favennec, Y. Rouizi, D. Petit, On the use of reduced models obtained
through identification for feedback optimal control problems in a heat
convection-diffusion problem, Computer Methods in Applied Mechanics
and Engineering 199 (17) (2010) 1193–1201.
[8] A. Haber, R. Fraanje, M. Verhaegen, Linear computational complexity
robust ILC for lifted systems, Automatica 48 (6) (2012) 1102–1110.
[9] P. Kunkel, V. Mehrmann, The linear quadratic optimal control problem for
linear descriptor systems with variable coefficients, Mathematics of Control,
Signals and Systems 10 (3) (1997) 247–264.
[10] G. Duan, Analysis and Design of Descriptor Linear Systems, Springer, 2010.
[11] M. H. Verhaegen, P. Van Dooren, A reduced order observer for descriptor
systems, Systems & Control Letters 8 (1) (1986) 29–37.
[12] K.-J. Bathe, Finite Element Procedures, Prentice Hall, 1996.
[13] S. S. Rao, The Finite Element Method in Engineering, ButterworthHeinemann, 2017.
[14] W.-H. Liu, A. H. Sherman, Comparative analysis of the Cuthill-McKee and
the reverse Cuthill-McKee ordering algorithms for sparse matrices, SIAM
Journal on Numerical Analysis 13 (2) (1976) 198–213.
[15] O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, The Finite Element Method:
Its Basis and Fundamentals, Vol. 1, Butterworth-Heinemann, 2005.
[16] B. D. O. Anderson, J. B. Moore, Optimal Control: Linear Quadratic Methods, Courier Dover Publications, 2007.
[17] H. Kwakernaak, R. Sivan, Linear Optimal Control Systems, Wileyinterscience New York, 1972.
[18] V. Simoncini, Computational methods for linear matrix equations, SIAM
Review 58 (3) (2016) 377–441.
28

[19] Z. Gajic, M. T. J. Qureshi, Lyapunov Matrix Equation in System Stability
and Control, Courier Corporation, 2008.
[20] K. Kodra, M. Skataric, Z. Gajic, Finding Hankel singular values for singularly perturbed linear continuous-time systems, IET Control Theory &
Applications (2017).
[21] B. Bamieh, O. Paganini, M. A. Dahleh, Distributed control of spatially invariant systems, IEEE Transactions on Automatic Control 47 (2002) 1091–
1118.
[22] R. D’Andrea, G. Dullerud, Distributed control design for spatially interconnected systems, IEEE Transactions on Automatic Control 48 (2003)
1478–1495.
[23] M. R. Jovanović, B. Bamieh, On the ill-posedness of certain vehicular platoon control problems, IEEE Transactions on Automatic Control 50 (9)
(2005) 1307–1321.
[24] F. Lin, M. Fardad, M. R. Jovanović, Optimal control of vehicular formations with nearest neighbor interactions, IEEE Transactions on Automatic
Control 57 (9) (2012) 2203–2218.
[25] U. A. Khan, J. M. F. Moura, Distributing the Kalman filter for large-scale
systems, IEEE Transactions on Signal Processing 56 (10) (2008) 4919–4935.
[26] F. Lin, M. Fardad, M. R. Jovanović, Design of optimal sparse feedback gains
via the alternating direction method of multipliers, IEEE Transactions on
Automatic Control 58 (9) (2013) 2426–2431.
[27] F. Lin, M. Fardad, M. R. Jovanović, Augmented Lagrangian approach to
design of structured optimal state feedback gains, IEEE Transactions on
Automatic Control 56 (12) (2011) 2923–2929.
[28] S. Schuler, P. Li, J. Lam, F. Allgöwer, Design of structured dynamic outputfeedback controllers for interconnected systems, International Journal of
Control 84 (12) (2011) 2081–2091.
[29] J. Rice, M. Verhaegen, Distributed Control: A sequentially semi-separable
approach for heterogeneous linear systems, IEEE Transactions on Automatic Control 54(6) (2009) 1270–1283.
[30] P. Massioni, M. Verhaegen, Distributed control for identical dynamically
coupled systems: A decomposition approach, IEEE Transactions on Automatic Control 54(1) (2009) 124–135.
[31] P. Benner, J. R. Li, T. Penzl, Numerical solution of large-scale Lyapunov
equations, Riccati equations, and linear-quadratic optimal control problems, Numerical Linear Algebra Appl 15 (9) (2008) 755–777.

29

[32] D. A. Bini, B. Iannazzo, B. Meini, Numerical Solution of Algebraic Riccati
Equations, SIAM, 2011.
[33] P. Benner, Solving large-scale control problems, IEEE Control Systems
24 (1) (2004) 44–59.
[34] V. L. Mehrmann, The Autonomous Linear Quadratic Control Problem:
Theory and Numerical Solution, Springer, 1991.
[35] A. Haber, M. Verhaegen, Sparse solution of the Lyapunov equation for
large-scale interconnected systems, Automatica 73 (2016) 256 – 268.
[36] V. Simoncini, The Lyapunov matrix equation. Matrix analysis from a computational perspective, arXiv preprint arXiv:1501.07564 (2015).
[37] M. Benzi, V. Simoncini, Decay bounds for functions of Hermitian matrices
with banded or Kronecker structure, SIAM Journal on Matrix Analysis and
Applications 36 (3) (2015) 1263–1282.
[38] M. Benzi, N. Razouk, Decay bounds and O(n) algorithms for approximating
functions of sparse matrices, Electronic Transactions on Numerical Analysis
(2007) 16–39.
[39] A. Haber, M. Verhaegen, Framework to trade optimality for local processing in large-scale wavefront reconstruction problems, Optics Letters 41 (22)
(2016) 5162–5165.
[40] A. Haber, M. Verhaegen, Moving horizon estimation for large-scale interconnected systems, IEEE Transactions on Automatic Control 58 (11)
(2013) 2834–2847.
[41] S. Pozza, V. Simoncini, Decay bounds for functions of banded nonHermitian matrices, arXiv preprint arXiv:1605.01595 (2016).
[42] M. Benzi, G. H. Golub, Bounds for the entries of matrix functions with
applications to preconditioning, BIT Numerical Mathematics 39 (3) (1999)
417–438.
[43] S. Demko, W. F. Moss, P. W. Smith, Decay rates for inverses of band
matrices, Mathematics of Computation 43 (168) (1984) 491–499.
[44] M. Benzi, Localization in matrix computations: Theory and applications,
Technical Report. Department of Mathematics and Computer Science,
Emory University, Atlanta.
[45] N. Motee, Q. Sun, Measuring sparsity in spatially interconnected systems,
in: Decision and Control (CDC), 2013 IEEE 52nd Annual Conference on,
IEEE, 2013, pp. 1520–1525.
[46] N. Motee, A. Jadbabaie, Distributed multi-parametric quadratic programming, IEEE Transactions on Automatic Control 54 (10) (2009) 2279–2289.
30

[47] N. Motee, Q. Sun, Sparsity and spatial localization measures for spatially
distributed systems, SIAM Journal on Control and Optimization 55 (1)
(2017) 200–235.
[48] M. R. Opmeer, T. Reis, W. Wollner, Finite-rank ADI iteration for operator
Lyapunov equations, SIAM Journal on Control and Optimization 51 (5)
(2013) 4084–4117.
[49] K. Ito, K. A. Morris, An approximation theory of solutions to operator
Riccati equations for H-infinity control, SIAM Journal on Control and Optimization 36 (1) (1998) 82–99.
[50] T. Huckle, Approximate sparsity patterns for the inverse of a matrix and
preconditioning, Applied Numerical Mathematics (30) (1999) 291–303.
[51] L. Grubišić, D. Kressner, On the eigenvalue decay of solutions to operator
Lyapunov equations, Systems & Control Letters 73 (2014) 42–47.
[52] P. Benner, M. Heinkenschloss, J. Saak, H. K. Weichelt, An inexact low-rank
Newton–ADI method for large-scale algebraic Riccati equations, Applied
Numerical Mathematics 108 (2016) 125–142.
[53] F. Feitzinger, T. Hylla, E. W. Sachs, Inexact Kleinman-Newton method
for Riccati equations, SIAM Journal on Matrix Analysis and Applications
31 (2) (2009) 272–288.
[54] A. Haber, MATLAB Codes for Computing Sparse Solutions of Lyapunov
Equations, http://aca.cabineer.com (2018).
[55] A. Haber, Estimation and Control of Large-Scale Systems with an Application to Adaptive Optics for EUV Lithography, Ph.D. thesis (2014).
[56] A. Haber, A. Polo, S. Ravensbergen, H. P. Urbach, M. Verhaegen, Identification of a dynamical model of a thermally actuated deformable mirror,
Optics Letters 38 (16) (2013) 3061–3064.
[57] R. Saathof, G. J. M. Schutten, J. W. Spronck, R. H. M. Schmidt, Actuation
profiles to form Zernike shapes with a thermal active mirror, Optics Letters
40 (2) (2015) 205–208.
[58] A. Haber, A. Polo, I. Maj, S. Pereira, H. Urbach, M. Verhaegen, Predictive
control of thermally induced wavefront aberrations, Optics Express 21 (18)
(2013) 21530–21541.
[59] A. Polo, A. Haber, S. F. Pereira, M. Verhaegen, H. P. Urbach, Linear phase
retrieval for real-time adaptive optics, Journal of the European Optical
Society-Rapid publications 8 (2013).
[60] J.-P. Chehab, M. Raydan, Inexact Newton‘s method with inner implicit
preconditioning for algebraic Riccati equations, Computational and Applied Mathematics 36 (2) (2017) 955–969.
31

[61] T. A. Davis, S. Rajamanickam, W. M. Sid-Lakhdar, A survey of direct
methods for sparse linear systems, Acta Numerica 25 (2016) 383–566.
[62] M. Benzi, Preconditioning techniques for large linear systems: A survey,
Journal of Computational Physics 182 (2) (2002) 418 – 477.
[63] Å. Björck, Numerical Methods for Least Squares Problems, SIAM, 1996.
[64] D. P. Bertsekas, Nonlinear Programming, Athena Scientific, 1999.
[65] A. J. Laub, Matrix Analysis for Scientists and Engineers, SIAM, 2005.
[66] M. Benzi, V. Simoncini, Approximation of functions of large matrices with
Kronecker structure, Numerische Mathematik 135 (1) (2017) 1–26.
[67] M. Benzi, P. Boito, Decay properties for functions of matrices over C∗ algebras, Linear Algebra and its Applications 456 (2014) 174–198.
[68] M. J. Grote, T. Huckle, Parallel preconditioning with sparse approximate
inverses, SIAM Journal on Scientific Computing 18 (3) (1997) 838–853.
[69] E. Chow, Y. Saad, Approximate inverse preconditioners via sparse-sparse
iterations, SIAM Journal on Scientific Computing 19 (3) (1998) 995–1023.
[70] W.-P. Tang, Toward an effective sparse approximate inverse preconditioner,
SIAM Journal on Matrix Analysis and Applications 20 (4) (1999) 970–986.
[71] E. Chow, A priori sparsity patterns for parallel sparse approximate inverse
preconditioners, SIAM Journal on Scientific Computing 21 (5) (2000) 1804–
1822.
[72] M. Benzi, M. Tuma, A sparse approximate inverse preconditioner for nonsymmetric linear systems, SIAM Journal on Scientific Computing 19 (3)
(1998) 968–994.
[73] L. Bergamaschi, M. Caliari, M. Vianello, Efficient approximation of the exponential operator for discrete 2D advection–diffusion problems, Numerical
Linear Algebra with Applications 10 (3) (2003) 271–289.
[74] L. Grasedyck, W. Hackbusch, B. N. Khoromskij, Solution of large scale
algebraic matrix Riccati equations by use of hierarchical matrices, Computing 70 (2) (2003) 121–165.
[75] S. Massei, D. Palitta, L. Robol, Solving rank structured Sylvester and
Lyapunov equations, arXiv preprint arXiv:1711.05493 (2017).
[76] D. Palitta, V. Simoncini, Numerical methods for large-scale Lyapunov
equations with symmetric banded data, arXiv preprint arXiv:1711.04187
(2017).

32

[77] D. Kressner, S. Massei, L. Robol, Low-rank updates and a divideand-conquer method for linear matrix equations, arXiv preprint
arXiv:1712.04349 (2017).
[78] N. Hale, N. J. Higham, L. N. Trefethen, Computing Aα , log(A), and related
matrix functions by contour integrals, SIAM Journal on Numerical Analysis
46 (5) (2008) 2505–2523.
[79] N. J. Higham, Functions of Matrices: Theory and Computation, SIAM,
2008.

33

