# On an integrated Krylov-ADI solver for large-scale Lyapunov equations

Peter Benner<EMAIL_ADDRESS>Davide Palitta
<EMAIL_ADDRESS>Jens Saak<EMAIL_ADDRESS>
###### Abstract

One of the most computationally expensive steps of the low-rank ADI method for
large-scale Lyapunov equations is the solution of a shifted linear system at
each iteration. We propose the use of the extended Krylov subspace method for
this task. In particular, we illustrate how a single approximation space can
be constructed to solve all the shifted linear systems needed to achieve a
prescribed accuracy in terms of Lyapunov residual norm. Moreover, we show how
to fully merge the two iterative procedures in order to obtain a novel,
efficient implementation of the low-rank ADI method, for an important class of
equations. Many state-of-the-art algorithms for the shift computation can be
easily incorporated into our new scheme, as well. Several numerical results
illustrate the potential of our novel procedure when compared to an
implementation of the low-rank ADI method based on sparse direct solvers for
the shifted linear systems.

###### keywords:

Lyapunov equations, low-rank ADI, extended Krylov method, shifted linear
systems.

15A24, 65F10, 65F30, 15A06

Extended Krylov subspace methods (EKSM) and the low-rank alternating
directions implicit (LR-ADI) iteration have been competing methods for the
solution of large-scale algebraic Lyapunov equations. In this paper, we make
an important step towards a new method merging them into a combined procedure
that inherits advantages from both worlds.

##  1 Introduction

The low-rank alternating direction implicit (LR-ADI) [42, 54] method is one of
the state-of-the-art methods for the numerical solution of large-scale
Lyapunov equations [19, 65]. This linear matrix equation can be encountered in
many applications: control and system theory [66, 34], especially some model
reduction techniques for dynamical systems [15, 3], but also discretization of
certain partial differential equations (PDEs) [71], and many more.

We consider Lyapunov equations of the form

$AXE^{\mkern-1.5mu\mathsf{T}}+EXA^{\mkern-1.5mu\mathsf{T}}+BB^{\mkern-1.5mu\mathsf{T}}=0,$
(1)

where $A,E\in\mathbb{R}^{n\times n}$, and $B\in\mathbb{R}^{n\times q}$, $q\ll
n$. Moreover, $E$ is supposed to be symmetric positive definite (SPD) and the
matrix pencil $(A,E)$ to be asymptotically stable, i.e., its spectrum is
contained in the open left half plane $\mathbb{C}_{-}$, which guarantees that
a unique solution $X$ exists, it is symmetric positive semidefinite [53].

A special case of equation (1) is attained whenever $E=I$, namely the equation
of interest is

$AX+XA^{\mkern-1.5mu\mathsf{T}}+BB^{\mkern-1.5mu\mathsf{T}}=0.$ (2)

Oftentimes the coefficient matrix $E$ possesses a structured sparsity pattern.
For instance, it is (block) diagonal when the matrices stem from a finite
element discretization that uses mass-lumping. In this case, we can easily
transform equation (1) and obtain an equation of the form (2). This can, for
example, be achieved by simply pre- and post-multiplying (1) by
$E^{-\frac{1}{2}}$ to potentially preserve symmetry of $A$. For the sake of
simplicity, we thus focus on equation (2) in the following.

In case of very large problem dimensions, the solution $X$ cannot be stored
since this matrix is, in general, dense. However, it is well-known that its
singular values quickly decay to zero under suitable assumptions, see, e.g.,
[55, 5, 33, 13], so that accurate low-rank approximations
$ZZ^{\mkern-1.5mu\mathsf{T}}\approx X$, $Z\in\mathbb{R}^{n\times t}$, $t\ll
n$, can be constructed. The efficient computation of the low-rank factor $Z$
is the task of LR-ADI and of all other _low-rank methods_. See, e.g., the
survey papers [65, 19] for further details on different low-rank methods for
linear matrix equations.

It is well-known that the convergence rate of the LR-ADI method is strictly
connected to the selection of some parameters
${\\{p_{i}\\}}_{i=1,\ldots,j}\subset\mathbb{C}_{-}$ called _shifts_ ††††††We
only consider _proper_ sets of shifts, namely ${\\{p_{i}\\}}_{i=1,\ldots,j}$
is closed with respect to complex conjugation.. The computation of effective
shifts is a highly non-trivial task and it has been a rather active research
topic in the last decades. Many strategies are available in the literature and
these can be divided into two categories: _Offline_ routines [54, 60, 73],
where the shifts are computed a-priori, before LR-ADI starts and then,
potentially, cyclically reused, and _online_ schemes [12, 37], where the
shifts are computed on the fly within the iterative procedure. The name shifts
for the values $p_{j}$ comes from the fact that in each LR-ADI iteration we
need to solve a shifted linear system with a coefficient matrix of the form
$A+p_{j}E$, or $A+p_{j}I$, in case of (1), or (2), respectively.  Notice that
since $(A,E)$ (or $A$ in case of (2)) is asymptotically stable and
$\\{p_{j}\\}\subset\mathbb{C}_{-}$, all the linear systems involved in the LR-
ADI scheme are well defined.

In Algorithm 1 we report an implementation of the LR-ADI scheme for the
solution of (1). Notice that Algorithm 1 is designed to drastically reduce the
amount of complex arithmetic calculations. Indeed, even though $A$ and $B$ in
(2) are real, the shifts $p_{j}$ are often complex if $A$ is nonsymmetric, so
that complex arithmetic may occur. See [11], [36, Chapter 4] and references
therein for details and derivations.

input : $A\in\mathbb{R}^{n\times n}$ stable, $B\in\mathbb{R}^{n\times q}$,
max. iteration count $j_{\max}$,
normalized residual bound $\varepsilon>0$.

output : $Z_{j}\in\mathbb{R}^{n\times jq}$,
$Z_{j}Z_{j}^{\mkern-1.5mu\mathsf{T}}=X_{j}\approx X$ approximate solution to
(2).

1 Set $W_{0}=B$, $Z_{0}=[]$, $j=1$, and select $p_{1}\in\mathbb{C}_{-}$

2 while  _$\|W_{j-1}^{*}W_{j-1}\|_{F}\geqslant\varepsilon\|B^{*}B\|_{F}$ and
$j\leqslant j_{\max}$_ do

3 Solve $(A+p_{j}I)S_{j}=W_{j-1}$

4 Set $W_{j}=W_{j-1}-2\text{Re}(p_{j})S_{j}$

5 if  _$\text{{Im}}(p_{j})\neq 0$_ then

6 Set $\beta=\text{Re}(p_{j})/\text{Im}(p_{j})$ and $p_{j+1}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{p}_{j}$

7 Set $S_{j+1}=\sqrt{2(\beta^{2}+1)}\text{Im}(S_{j})$

8 Set
$W_{j+1}=W_{j-1}-4\text{Re}(p_{j})(\text{Re}(S_{j})+\beta\text{Im}(S_{j}))$

9 Set $S_{j}=\sqrt{2}(\text{Re}(S_{j})+\beta\text{Im}(S_{j}))$

10 Set
$Z_{j+1}=[Z_{j-1},\sqrt{-2\text{Re}(p_{j})}S_{j},\sqrt{-2\text{Re}(p_{j+1})}S_{j+1}]$

11 Set $j=j+1$

12

13 else

14 Set $Z_{j}=[Z_{j-1},\sqrt{-2\text{Re}(p_{j})}S_{j}]$

15

16 Choose the next shift $p_{j+1}\in\mathbb{C}_{-}$

17 Set $j=j+1$

18

Algorithm 1 LR-ADI for Lyapunov equations.

One of the most computationally expensive steps of Algorithm 1 is the solution
of the shifted linear systems with $q$ right-hand sides in line 1. Such a job
has to be carried out at each LR-ADI iteration. In this contribution, we
propose to employ state-of-the-art block Krylov subspace methods for this
task. In particular, for equation (2), we illustrate how to efficiently reuse
the approximation space employed at the $j$-th LR-ADI iteration and utilize it
also in the next one. To this end, it is crucial that the right-hand side of
the linear system we need to solve at the $(j+1)$-st iteration can be
represented in terms of the basis of the subspace employed in the previous
iteration. This simple but critical observation lets us design a novel,
efficient procedure that can lead to noticeable savings in the running time
for the solution of (2). Indeed, all the LR-ADI steps can be completely merged
into the Krylov routine so that the LR-ADI iteration is only implicitly
performed. Moreover, also the LR-ADI shift computation can be incorporated
into the framework proposed in this paper.

The following is a synopsis of the paper. Section 2 is devoted to recalling
the general (block) Krylov subspace framework for shifted linear systems. In
particular, some details about the extended Krylov subspace method presented
in [64] are given in Section 2.1. In Section 3 we present the main
contribution of the paper and we show how to fully merge the LR-ADI iteration
into the projection method adopted for the linear system solution. The
selection of effective shifts is crucial for attaining a fast convergence in
terms of number of LR-ADI iterations, and numerous strategies have been
proposed in the literature to accomplish this task; see, e.g., [71, 54, 60,
16, 73, 12, 37, 58]. In Section 5 we illustrate how many of these routines can
be integrated into our novel framework with no additional cost. The potential
of our strategy is depicted in Section 6, where several numerical results are
reported. We close the paper with our conclusions in Section 7.

Throughout the paper, we adopt the following notation. The matrix inner
product is defined as $\langle
X,Y\rangle_{F}\colon=\operatorname*{trace}(Y^{\mkern-1.5mu\mathsf{T}}X)$ so
that the induced norm is $\|X\|_{F}^{2}=\langle X,X\rangle_{F}$. The Kronecker
product is denoted by $\otimes$ whereas $I_{n}$ and $O_{n\times m}$ denote the
identity matrix of order $n$ and the $n\times m$ zero matrix, respectively.
Only one subscript is used for a square zero matrix, i.e., $O_{n\times
n}=O_{n}$, and the subscript is omitted whenever the dimension of $I$ and $O$
is clear from the context. Moreover, $e_{i}$ is the $i$-th basis vector of the
canonical basis of $\mathbb{R}^{n}$. The brackets $[\cdot]$ are used to
concatenate matrices of conformal dimensions. In particular, a MATLAB-like
notation is adopted and $[M,N]$ denotes the matrix obtained by putting $M$ on
the left of $N$ whereas $[M;N]$ the one obtained by putting $M$ on top of $N$,
i.e.,
$[M;N]={[M^{\mkern-1.5mu\mathsf{T}},N^{\mkern-1.5mu\mathsf{T}}]}^{\mkern-1.5mu\mathsf{T}}$.
If $w\in\mathbb{R}^{n}$, $\text{diag}(w)$ denotes the $n\times n$ diagonal
matrix whose $i$-th diagonal entry corresponds to the $i$-th component of $w$.
Given $X\in\mathbb{C}^{n\times m}$, we write
$X=\text{Re}(X)+\imath\text{Im}(X)$, where $\text{Re}(X)$ and $\text{Im}(X)$
are its real and imaginary parts, respectively, and $\imath$ is the imaginary
unit. The complex conjugate of $X$ is denoted by $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{X}=\text{Re}(X)-\imath\text{Im}(X)$.

##  2 Block Krylov methods for shifted linear systems

The literature about the numerical solution of shifted linear systems by
Krylov subspace methods is rather vast. Indeed, sequences of shifted linear
systems arise in many applications belonging to different research areas like
control theory [23, 41], wave propagation problems [8], mechanical systems
[27], quantum chromodynamics [32], and many more.

This algebraic problem is trickier than it looks and many researchers have
contributed to its understanding providing important insights on its
properties and designing efficient, robust algorithms for its solution. Here
is an incomplete list of contributions on numerical schemes for sequences of
shifted linear systems and their analysis [62, 28, 29, 48, 70, 67, 68, 7].

In this section, we consider sequences of shifted linear systems of the form

$(A+p_{j}I)S_{j}=W,\quad W\in\mathbb{R}^{n\times q},$ (3)

where the right-hand side $W$ does not depend on the index $j$, even though,
in line 1 of Algorithm 1, $W_{j-1}$ does change at every LR-ADI iteration. In
Section 3 we show how to adapt the machinery, presented here, to the case of
linear systems of the form $(A+p_{j}I)S_{j}=W_{j-1}$, arising within the LR-
ADI scheme.

Any Krylov routine  for (3) computes a numerical solution of the form
$S_{m}^{(j)}=S_{0}+V_{m}Y_{m}^{(j)}\approx S_{j}$,
$V_{m}=[\mathcal{V}_{1},\ldots,\mathcal{V}_{m}]\in\mathbb{R}^{n\times m\ell
q}$, $\ell\geqslant 1$‡‡‡‡‡‡The value of $\ell$ depends on the adopted
approximation space. It holds, $\ell=1$ for the polynomial Krylov subspace in
(4), whereas $\ell=2$ for the extended Krylov subspace in (5).,
$\mathcal{V}_{i}\in\mathbb{R}^{n\times\ell q}$, $i=1,\ldots,m$,
$Y_{m}^{(j)}\in\mathbb{C}^{m\ell q\times q}$, where the orthonormal columns of
$V_{m}$ span a suitable subspace $\mathcal{K}_{m}$, namely,
$\operatorname*{Range\,}(V_{m})=\mathcal{K}_{m}$, $S_{0}$ is an initial guess,
and the matrix $Y_{m}^{(j)}$ can be computed by imposing different conditions.
In particular, $Y_{m}^{(j)}$ is often computed by either imposing a Galerkin
condition on the residual or minimizing the residual norm. For the sake of
simplicity, we consider $S_{0}=O$ in the following.

One of the most common choices for the approximation space $\mathcal{K}_{m}$
is the block Krylov subspace

$\mathbf{K}_{m}^{\square}(A,W)=\operatorname*{Range\,}([W,AW,\ldots,A^{m-1}W]).$
(4)

See, e.g., [50, 61, 30] and the references therein for further details on the
block polynomial Krylov subspace $\mathbf{K}_{m}^{\square}(A,W)$ and related
methods.

However, Simoncini showed in [64] that the extended Krylov subspace [24]

$\mathbf{EK}_{m}^{\square}(A,W)=\operatorname*{Range\,}([W,A^{-1}W,AW,A^{-2}W,\ldots,A^{m-1}W,A^{-m}W]),$
(5)

can be a powerful alternative for the solution of (3) in many cases. For
instance, when $A$ is large and real while the $p_{j}$’s are complex. See also
Section 2.1.

The basis $V_{m}$ of both the polynomial and extended Krylov subspace can be
constructed by means of the (extended) Arnoldi process and the following
Arnoldi relation is fulfilled

$AV_{m}=V_{m}T_{m}+\mathcal{V}_{m+1}E_{m+1}^{\mkern-1.5mu\mathsf{T}}\underline{T}_{m},$
(6)

where
$\underline{T}_{m}=V_{m+1}^{\mkern-1.5mu\mathsf{T}}AV_{m}\in\mathbb{R}^{(m+1)\ell
q\times m\ell q}$, $T_{m}$ is its principal square submatrix, and
$E_{m+1}=e_{m+1}\otimes I_{\ell q}$. See, e.g., [56, 63].

The Arnoldi relation (6) is one of the most crucial tools in the solution of
(3) by Krylov methods. Indeed, it can be used to show the fundamental _shift-
invariance_ property of the Krylov subspaces (4) and (5), and the following
relation holds true

$(A+p_{j}I_{n})V_{m}=V_{m}(T_{m}+p_{j}I_{m\ell
q})+\mathcal{V}_{m+1}E_{m+1}^{\mkern-1.5mu\mathsf{T}}\underline{T}_{m}.$ (7)

See, e.g., [62, Equation (2.1)], [64, Equation (3.1)].

Equation (7) says that we can compute only one approximation space for solving
(3). In particular, the space constructed using $A$, i.e.,
$\mathbf{K}_{m}^{\square}(A,W)$ or $\mathbf{EK}_{m}^{\square}(A,W)$, can be
employed, by possibly being expanded, to solve all the shifted linear systems
in the sequence (3).

Polynomial Krylov subspace methods often need many iterations to achieve the
prescribed accuracy, so that a large subspace is constructed. This leads to an
increment in both the storage demand and the computational efforts of the
selected solution procedure. Different strategies have been developed to avoid
the construction of a too large subspace.

With the goal of achieving a fast convergence in terms of number of
iterations, the linear system (3) can be preconditioned, namely is transformed
into an equivalent problem with better spectral properties. However, designing
effective preconditioning operators for a sequence of shifted linear systems
is a difficult task and often highly problem dependent. Very sophisticated
schemes have been proposed in the literature. See, e.g., [46, 9, 20, 21, 4].

Restarted routines are an alternative solution. In this framework, the
approximation space $\mathcal{K}_{m}$ is expanded until it reaches a
prescribed maximum dimension. If the desired level of accuracy is not
achieved, the last computed basis block $\mathcal{V}_{m+1}$ is employed as
initial block in the construction of a new subspace
$\mathcal{K}_{m}^{\prime}$. This procedure is iterated until a stopping
criterion is fulfilled. See, e.g., [62, 29] and [30, Section 3.2.1]. However,
in our framework the LR-ADI shifts $p_{j}$’s are often computed on the fly
and, thus, are not all available at the same time. Therefore, to fully take
advantage of the computational efforts needed to solve the linear system
$(A+p_{j-1}I)S_{j-1}=W$, we would have to store all the bases computed during
the employed restarted Krylov procedure and use them to solve the $j$-th
linear system, as well. Unfortunately, this would destroy all the benefits in
terms of storage complexity gained from the restart-paradigm.

In [64], Simoncini showed that the employment of the extended Krylov subspace
(5), in place of (4), often leads to a faster convergence, in terms of
iterations, to the point that the constructed subspace is usually smaller than
the polynomial counterpart needed to reach the same level of accuracy. We,
thus, decide to use such an approximation space for the solution of the
shifted linear systems within the LR-ADI method and in the next section we
recall some details of the extended Krylov subspace method.

Notice that the faster convergence of the extended Krylov subspace (5) comes
with a toll. Indeed, at each iteration, a linear system with $A$ has to be
solved during the basis construction. Nevertheless, the increase in the
overall workload of the solution process can be limited in general. Indeed, if
we want to use a direct solver to invert $A$, for instance, the LU factors of
$A$ can be computed once and for all before the LR-ADI scheme starts. On the
other hand, if an iterative procedure is employed, analogously a single
preconditioner for $A$ has to be designed once.

As already mentioned, in the formulation (3) the right-hand side $W$ is fixed,
namely it does not depend on the shift index $j$. However, in line 1 of
Algorithm 1, the linear systems we need to solve are of the form

$(A+p_{j}I)S_{j}=W_{j-1}.$

At a first glance, having a nonconstant right-hand side does not allow for the
employment of the shifted Krylov framework we briefly described above. A
larger class of solvers, the so-called _recycling_ Krylov methods, seems more
appropriate. See, e.g., [67, 52, 70, 31, 69] for general sequences of shifted
linear systems, and [2, 1, 26] for some recycling Krylov techniques applied in
a model reduction context. However, in Section 3 we show that, in the LR-ADI
context for $j>1$, the residual factor $W_{j-1}$ belongs to the subspace
$\mathcal{K}_{m}$ employed in the solution of the $(j-1)$-st linear system
$(A+p_{j-1}I)S_{j-1}=W_{j-2}$. Along with the shift-invariance property of the
Krylov subspace, this observation allows us to utilize only one subspace for
the solution of all the shifted linear systems within the LR-ADI method. In
turn, as shown in Section 6, we can notably reduce the computational effort of
the overall procedure.

###  2.1 The extended Krylov subspace method for shifted linear systems

In this section, we recall the extended Krylov subspace method for shifted
linear systems presented in [64].

Given the sequence of shifted linear systems (3), the extended Krylov subspace
method computes a solution of the form $S_{m}^{(j)}=V_{m}Y_{m}^{(j)}$, where
the $2mq$ orthonormal columns of $V_{m}$ span the extended Krylov subspace
(5), whereas the $2mq\times q$ matrix $Y_{m}^{(j)}$ can be computed in
different manners.

For instance, $Y_{m}^{(j)}$ can be computed by imposing a Galerkin condition
on the residual $R_{m}^{(j)}=(A+p_{j}I)V_{m}Y_{m}^{(j)}-W$, namely by imposing
$V_{m}^{\mkern-1.5mu\mathsf{T}}R_{m}^{(j)}=0$. Thanks to the shifted Arnoldi
relation (7), it is easy to show that such a Galerkin condition is equivalent
to solving the projected linear systems

$(T_{m}+p_{j}I)Y_{m}^{(j)}=E_{1}\gamma,$ (8)

where $E_{1}=e_{1}\otimes I_{2q}$, and $\gamma\in\mathbb{R}^{2q\times q}$ is
such that $W=V_{1}\gamma$.

With $Y_{m}^{(j)}$ at hand, the Frobenius norm of the residual
$\|R_{m}^{(j)}\|_{F}$ can be computed at low cost, as

$\|R_{m}^{(j)}\|_{F}=\|E_{m+1}^{\mkern-1.5mu\mathsf{T}}\underline{T}_{m}Y_{m}^{(j)}\|_{F},$
(9)

following [64, Equation (3.2)].

Alternatively, following the discussion in [68, Section 4.1], the matrix
$Y_{m}^{(j)}$ can be computed also by minimizing the residual norm, i.e.,

$Y_{m}^{(j)}=\operatorname*{argmin}_{Y\in\mathbb{R}^{2mq\times
q}}\|(A+p_{j}I)V_{m}Y-W\|_{F}.$ (10)

Once again, thanks to the shifted Arnoldi relation (7), the minimization
problem in (10) simplifies, and we can compute $Y_{m}^{(j)}$ as

$Y_{m}^{(j)}=\operatorname*{argmin}_{Y\in\mathbb{R}^{2mq\times
q}}\|(\underline{T}_{m}+p_{j}[I_{2mq};O_{2q\times 2mq}])Y-E_{1}\gamma\|_{F}.$
(11)

Note the abuse of notation in (11): $E_{1}\in\mathbb{R}^{2(m+1)q\times 2q}$
whereas $E_{1}\in\mathbb{R}^{2mq\times 2q}$ in (8).

If $QP=\underline{T}_{m}+p_{j}[I_{2mq};O_{2q\times 2mq}]$ denotes the QR
factorization of $\underline{T}_{m}+p_{j}[I_{2mq};O_{2q\times 2mq}]$, and we
consider the following partition

$Q=[Q_{1},Q_{2}],\,Q_{1}\in\mathbb{R}^{2(m+1)q\times
2mq},\,Q_{2}\in\mathbb{R}^{2(m+1)q\times 2q},\;P=\begin{bmatrix}P_{1}\\\
O_{2q\times 2mq}\\\ \end{bmatrix},\,P_{1}\in\mathbb{R}^{2mq\times 2mq},$

then the matrix $Y_{m}^{(j)}$ in (11) can be computed as

$Y_{m}^{(j)}=P_{1}^{-1}Q_{1}^{\mkern-1.5mu\mathsf{T}}E_{1}\gamma,$ (12)

and the residual norm is given by

$\|R_{m}^{(j)}\|_{F}=\|Q_{2}^{\mkern-1.5mu\mathsf{T}}E_{1}\gamma\|_{F}.$ (13)

The overall procedure is summarized in Algorithm 2,  where $\Sigma$ contains
the indices of all yet unsolved systems, whereas $\Sigma_{C}$ contains the
indices of all the systems that have already been solved. The basis block
$\mathcal{V}_{m+1}$ can be computed by following [63]. This operation involves
both matrix-vector products and linear system solves with $A$. Moreover, the
basis $V_{m}$ is real whenever $A$ and $W$ are so. Complex arithmetic may
occur in the computation of $Y_{m}^{(j)}$, if $\text{Im}(p_{j})\neq 0$.

Notice that as soon as the $j$-th linear system has converged, namely the
related relative residual norm is sufficiently small, we stop solving the
$j$-th projected problem******Either (8) or (11).. Once all the linear systems
have converged, we terminate the iterative process.

To conclude, we would like to point out that, to the best of our knowledge,
this is the first time the minimal residual condition (11) is proposed within
the extended Krylov subspace method for shifted linear systems.

input : $A\in\mathbb{R}^{n\times n}$, $\\{p_{j}\\}\subset\mathbb{C}$,
$W\in\mathbb{R}^{n\times q}$, and normalized residual bound $\varepsilon>0$

output : $S_{m}^{(j)}\in\mathbb{R}^{n\times q}$,
$S_{m}^{(j)}=V_{m}Y_{m}^{(j)}\approx S_{j}$,  where $(A+p_{j}I)S_{j}=W_{j}$

1 Set $\beta=\|W\|_{F}$, $\Sigma=\\{1,\ldots,\max j\\}$,
$\Sigma_{C}=\emptyset$

2 Perform economy-size QR,
$[W,A^{-1}W]=[\mathcal{V}_{1}^{(1)},\mathcal{V}_{1}^{(2)}][\gamma,\theta]$,
$\gamma,\theta\in\mathbb{R}^{2q\times q}$

3 for  _$m=1,2,\dots$ , till convergence,_ do

4 Compute next basis block $\mathcal{V}_{m+1}$ as in [63] and set
$V_{m+1}=[V_{m},\mathcal{V}_{m+1}]$

5 Update $T_{m}=V_{m}^{\mkern-1.5mu\mathsf{T}}AV_{m}$ as in [63]

6 Compute $Y_{m}^{(j)}$ for $j\in\Sigma\diagdown\Sigma_{C}$ as in (8) or (12)

7 Compute $\|R_{m}^{(j)}\|_{F}$ accordingly as in (9) or (13)

8

9 if  _$\|R_{m}^{(j)}\|_{F}/\beta <\varepsilon$_ then

10 Set $\Sigma_{C}=\Sigma_{C}\cup\\{j\\}$

11

12 if  _$\Sigma\diagdown\Sigma_{C}=\emptyset$_ then

13 Break and go to 2

14

15Set $S_{m}^{(j)}=V_{m}Y_{m}^{(j)}$

Algorithm 2 Extended Krylov subspace method for shifted linear systems.

##  3 Merging the two iterative procedures

In this section we show how the LR-ADI iteration and the extended Krylov
subspace method for shifted linear systems can be merged together into a
novel, efficient iterative procedure for the solution of (2).

As already mentioned, in the sequence of shifted linear systems in line 1 of
Algorithm 1, also the right-hand side $W_{j-1}$ depends on the current LR-ADI
iteration $j$. Therefore, at a first glance, we seemingly have to build a new
subspace at each iteration $j$, by employing the current $W_{j-1}$ as initial
block. However, in the following theorem we show that $W_{j-1}$ belongs to the
subspace constructed to solve the $(j-1)$-st linear system so that such a
space can be used, by being possibly expanded, also in the solution of the
subsequent linear system.

######  Theorem 3.1.

Let $S_{j}=V_{m_{j}}Y_{m_{j}}$, $j\geqslant 1$,
$\operatorname*{Range\,}(V_{m_{j}})=\mathbf{EK}_{m_{j}}^{\square}(A,B)$ for
certain $m_{j}\geqslant 0$. Then

$\operatorname*{Range\,}(W_{j})\subseteq\mathbf{EK}_{m_{j}}^{\square}(A,B).$

###### Proof.

We are going to show the statement by induction on $j$.

The first linear system to be solved within the LR-ADI method is
$(A+p_{1}I)S_{1}=B$ and the extended Krylov subspace
$\mathbf{EK}_{m_{1}}^{\square}(A,B)$ can be employed to this end. The computed
solution is of the form $S_{1}=V_{m_{1}}Y_{m_{1}}$, $m_{1}>0$, where
$\operatorname*{Range\,}(V_{m_{1}})=\mathbf{EK}_{m_{1}}^{\square}(A,B)$ and
$Y_{m_{1}}\in\mathbb{C}^{2m_{1}q\times q}$. It is thus easy to show that
$W_{1}=B-2\text{Re}(p_{1})S_{1}=V_{m_{1}}(E_{1}\gamma-2\text{Re}(p_{1})Y_{m_{1}})$
is such that
$\operatorname*{Range\,}(W_{1})\subseteq\mathbf{EK}_{m_{1}}^{\square}(A,B)$.

We now assume the statement holds for a certain $j-1\geqslant 1$, and we show
it holds for $j$ as well. Since $S_{j}=V_{m_{j}}Y_{m_{j}}$ by assumption and
$\operatorname*{Range\,}(W_{j-1})\subseteq\mathbf{EK}_{m_{j-1}}^{\square}(A,B)$
by inductive hypothesis, namely we can write
$W_{j-1}=V_{m_{j-1}}\Upsilon_{j-1}$ for a certain
$\Upsilon_{j-1}\in\mathbb{R}^{2m_{j-1}q\times q}$, we have

$\displaystyle W_{j}=$ $\displaystyle
W_{j-1}-2\text{Re}(p_{j})S_{j}=V_{m_{j-1}}\Upsilon_{j-1}-2\text{Re}(p_{j})V_{m_{j}}Y_{m_{j}}$
$\displaystyle=$ $\displaystyle
V_{m_{j}}\left([\Upsilon_{j-1};O_{2(m_{j}-m_{j-1})q\times
q}]-2\text{Re}(p_{j})Y_{m_{j}}\right).$

Therefore,
$\operatorname*{Range\,}(W_{j})\subseteq\mathbf{EK}_{m_{j}}^{\square}(A,B)$. ∎

Theorem 3.1 shows that $W_{j}$ is exactly represented in
$\mathbf{EK}_{m_{j}}^{\square}(A,B)$. This means that the latter subspace can
be still employed for the computation of $S_{j+1}$ by being possibly expanded.
Indeed, no components of $W_{j}$ are annihilated when either the Galerkin or
the minimal residual condition is imposed. In the following corollary we show
how to easily write down the projected problems (8) and (11) along with the
corresponding residual norm computation.

######  Corollary 3.1.

Assume the prerequisites of Theorem 3.1 hold. If a Galerkin condition is
imposed for the computation of $S_{j}=V_{m_{j}}Y_{m_{j}}$, then the matrix
$Y_{m_{j}}$ amounts to the solution of the projected linear system

$(T_{m_{j}}+p_{j}I_{2m_{j}q})Y_{m_{j}}=[\Upsilon_{j-1};O_{2(m_{j}-m_{j-1})q\times
q}],$ (14)

where $\Upsilon_{j-1}\in\mathbb{R}^{2m_{j-1}q\times q}$ is such that
$W_{j-1}=V_{m_{j-1}}\Upsilon_{j-1}$, $m_{j-1}\leqslant m_{j}$. The related
residual norm can be computed by

$\|R_{m_{j}}\|_{F}=\|E_{m_{j}+1}^{\mkern-1.5mu\mathsf{T}}\underline{T}_{m_{j}}Y_{m_{j}}\|_{F}.$
(15)

Similarly, if a minimal residual norm condition is imposed, we have

$Y_{m_{j}}=\operatorname*{argmin}_{Y\in\mathbb{C}^{2m_{j}q\times
q}}\|(\underline{T}_{m_{j}}+p_{j}[I_{2m_{j}q};O_{2q\times
2m_{j}q}])Y-[\Upsilon_{j-1};O_{2(m_{j}-m_{j-1}+1)q\times q}]\|_{F},$ (16)

so that

$\|R_{m_{j}}\|_{F}=\|Q_{2}^{\mkern-1.5mu\mathsf{T}}[\Upsilon_{j-1};O_{2(m_{j}-m_{j-1}+1)q\times
q}]\|_{F},$ (17)

where the $2q$ orthonormal columns of $Q_{2}$ are a basis of the kernel of
$\underline{T}_{m_{j}}+p_{j}[I_{2m_{j}q};O_{2q\times 2m_{j}q}]$.

###### Proof.

Since $W_{j-1}=V_{m_{j-1}}\Upsilon_{j-1}$ and we look for a solution
$S_{j}=V_{m_{j}}Y_{m_{j}}$ to $(A+p_{j}I)S_{j}=W_{j-1}$, we can write

$\begin{array}[]{rcl}R_{m_{j}}&=&(A+p_{j}I)S_{j}-W_{j-1}=(A+p_{j}I)V_{m_{j}}Y_{m_{j}}-V_{m_{j-1}}\Upsilon_{j-1}\\\
&=&V_{m_{j}}\left((T_{m_{j}}+p_{j}I_{2m_{j}q})Y_{m_{j}}-[\Upsilon_{j-1};O_{2(m_{j}-m_{j-1})q\times
q}]\right)+\mathcal{V}_{m_{j}+1}E_{m_{j}+1}^{\mkern-1.5mu\mathsf{T}}\underline{T}_{m_{j}}\\\
&=&V_{m_{j}+1}\left((\underline{T}_{m_{j}}+p_{j}[I_{2m_{j}q};O_{2q\times
2m_{j}q}])Y_{m_{j}}-[\Upsilon_{j-1};O_{2(m_{j}-m_{j-1}+1)q\times
q}]\right).\end{array}$

If a Galerkin condition is imposed, namely
$V_{m_{j}}^{\mkern-1.5mu\mathsf{T}}R_{m_{j}}=0$, then $Y_{m_{j}}$ is the
solution of the linear system in (14) and the related residual norm
$\|R_{m_{j}}\|_{F}$ can be computed as in (15).

Similarly, if a minimal residual condition is imposed, $Y_{m_{j}}$ solves the
minimization problem (16) and the $\|R_{m_{j}}\|_{F}$ fulfills (17). ∎

Once $S_{j}=V_{m_{j}}Y_{m_{j}}$ is computed, namely the related residual norm
$\|R_{m_{j}}\|_{F}$ is sufficiently small, we proceed with the remaining LR-
ADI operations.

We would like to point out that the expression of $W_{j}$, i.e.,
$W_{j}=V_{m_{j}}\Upsilon_{j}$, can be exploited for the Lyapunov residual norm
as well. Indeed,

$\|W_{j}^{*}W_{j}\|_{F}=\|\Upsilon_{j}^{*}\Upsilon_{j}\|_{F}.$ (18)

This means that also the computation of the Lyapunov residual norm can be
carried out by manipulating small matrices of dimension $2m_{j}q\times q$.
Similarly, the solution $Z_{j}$ can be assembled at the very end of the LR-ADI
procedure once the residual norm in (18) is sufficiently small. Indeed,

$\begin{array}[]{rcl}Z_{j}&=&[Z_{j-1},\sqrt{-2\text{Re}(p_{j})}S_{j}]=[\sqrt{-2\text{Re}(p_{1})}S_{1},\sqrt{-2\text{Re}(p_{2})}S_{2},\ldots,\sqrt{-2\text{Re}(p_{j})}S_{j}]\\\
&&\\\
&=&[\sqrt{-2\text{Re}(p_{1})}V_{m_{1}}Y_{m_{1}},\sqrt{-2\text{Re}(p_{2})}V_{m_{2}}Y_{m_{2}},\ldots,\sqrt{-2\text{Re}(p_{j})}V_{m_{j}}Y_{m_{j}}]\\\
&&\\\ &=&V_{m_{j}}[[Y_{m_{1}};O_{2(m_{j}-m_{1})q\times
q}],[Y_{m_{2}};O_{2(m_{j}-m_{2})q\times q}],\ldots,Y_{m_{j}}]\\\
&&\cdot(\sqrt{-2\text{diag}(\text{Re}(p_{1}),\ldots,\text{Re}(p_{j}))}\otimes
I_{q}).\end{array}$ (19)

The overall procedure combining the LR-ADI iteration with the extended Krylov
subspace method for shifted linear systems is depicted in Algorithm
3*†*†*†Many subscripts have been removed to make the algorithm more readable..

input : $A\in\mathbb{R}^{n\times n}$, $A$ stable, $B\in\mathbb{R}^{n\times
q}$, max. inner and outer iteration count $m_{\max}$, $j_{\max}$, normalized
residual bound $\varepsilon_{\mathtt{out}}>0$

output : $Z_{j}\in\mathbb{R}^{n\times jq}$,
$Z_{j}Z_{j}^{\mkern-1.5mu\mathsf{T}}\approx X$ approximate solution to (2).

1 Set $\nu=\|B^{*}B\|_{F}$, $m_{0}=0$, $Y_{0}=O$, and select
$p_{1}\in\mathbb{C}_{-}$

2 Perform economy-size QR,
$[B,A^{-1}B]=[\mathcal{V}_{1}^{(1)},\mathcal{V}_{1}^{(2)}][\gamma,\theta]$,
$\gamma,\theta\in\mathbb{R}^{2q\times q}$

3 Set $\Upsilon_{0}=E_{1}\gamma$, $m=j=1$, and select
$\varepsilon_{\mathtt{inn}}^{(1)}$

4 while  _$m\leqslant m_{\max}$ and $j\leqslant j_{\max}$_ do

5 Compute next basis block $\mathcal{V}_{m+1}$ as in [63] and set
$V_{m+1}=[V_{m},\mathcal{V}_{m+1}]$

6 Update $T_{m}=V_{m}^{\mkern-1.5mu\mathsf{T}}AV_{m}$ as in [63]

7 Compute $Y_{m}$ as in (14) or (16)

8 Compute $\|R_{m}\|_{F}$ accordingly as in (15) or (17)

9

10 if  _$\|R_{m}\|_{F}\leqslant\varepsilon_{\mathtt{inn}}^{(j)}$_ then

11 $\mathtt{flag\\_noexpand}=1$

12 while  _$\mathtt{flag\\_noexpand}$_ do

13 Set $m_{j}=m$ and $Y_{m_{j}}=Y_{m}$

14 Set $\Upsilon_{j}=\Upsilon_{j-1}-2\text{Re}(p_{j})Y_{m_{j}}$

15 if  _$\text{{Im}}(p_{j})\neq 0$_ then

16 Set $m_{j+1}=m_{j}$, $\beta=\text{Re}(p_{j})/\text{Im}(p_{j})$, and
$p_{j+1}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{p}_{j}$

17 Set $Y_{m_{j+1}}=\sqrt{2(\beta^{2}+1)}\text{Im}(Y_{m_{j}})$

18 Set
$\Upsilon_{j+1}=\Upsilon_{j-1}-4\text{Re}(p_{j})(\text{Re}(Y_{m_{j}})+\beta\text{Im}(Y_{m_{j}}))$

19 Set $Y_{m_{j}}=\sqrt{2}(\text{Re}(Y_{m_{j}})+\beta\text{Im}(Y_{m_{j}}))$

20 Set $j=j+1$

21

22 if
_$\|\Upsilon_{j}^{*}\Upsilon_{j}\|_{F}\leqslant\nu\cdot\varepsilon_{\mathtt{out}}$_
then

23 Break and go to 3

24 Choose the next shift $p_{j+1}\in\mathbb{C}_{-}$

25 Set $j=j+1$

26 Compute $Y_{m}$ as in (14) or (16)

27 Compute $\|R_{m}\|_{F}$ accordingly as in (15) or (17)

28 Select $\varepsilon_{\mathtt{inn}}^{(j)}$

29 if  _$\|R_{m}\|_{F}\geqslant\varepsilon_{\mathtt{inn}}^{(j)}$_ then

30 $\mathtt{flag\\_noexpand}=0$

31

32

33

34 Set $m=m+1$

35

36

$Z_{j}=V_{m}[[Y_{m_{1}};O_{2(m_{j}-m_{1})q\times
q}],[Y_{m_{2}};O_{2(m_{j}-m_{2})q\times
q}],\ldots,Y_{m_{j}}](\sqrt{-2\text{diag}(\text{Re}(p_{1}),\ldots,\text{Re}(p_{j}))}\otimes
I_{q})$

Algorithm 3 LR-ADI-EKSM for Lyapunov equations.

As in Algorithm 1, if $\text{Im}(p_{j})\neq 0$, in lines 3 to 3 we set
$p_{j+1}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{p}_{j}$, and we follow the implementation suggested in [11, 36] to reduce
the amount of complex arithmetic. In particular, $Y_{m_{j+1}}$ can be obtained
from $Y_{m_{j}}$ without solving (14) or (16). Moreover, the adopted scheme
results in a real $Z_{j}$. See [11] and [36, Algorithm 4.3] for further
details.

######  Remark 4.

Theorem 3.1 shows that
$\text{Range}(W_{j})\subseteq\mathbf{EK}_{m_{j}}^{\square}(A,B)$ whenever
$W_{j}$ is updated as $W_{j}=W_{j-1}-2\text{Re}(p_{j})S_{j}$, namely whenever
all the employed shifts are real. In case of shifts with nonzero imaginary
part, the LR-ADI implementation we adopt sets

$W_{j+1}=W_{j-1}-4\text{Re}(p_{j})(\text{Re}(S_{j})+\beta\text{Im}(S_{j})).$

Therefore, we need to show that $W_{j+1}$ defined as above is still such that
$\text{Range}(W_{j+1})\subseteq\mathbf{EK}_{m_{j}}^{\square}(A,B)$. This can
be done by applying the same exact arguments used in the proof of Theorem 3.1.
In particular, the result follows by noticing that the basis $V_{m}$ is real,
as we assumed $A$ and $B$ to be real matrices, and that we can write

$\displaystyle W_{j+1}=$ $\displaystyle
W_{j-1}-4\text{Re}(p_{j})(\text{Re}(S_{j})+\beta\text{Im}(S_{j}))$
$\displaystyle=$ $\displaystyle
V_{m_{j}}\left([\Upsilon_{j-1};O_{2(m_{j}-m_{j-1})q\times
q}]-4\text{Re}(p_{j})(\text{Re}(Y_{m_{j}})+\beta\text{Im}(Y_{m_{j}})\right).$

Notice that two tolerances $\varepsilon_{\mathtt{inn}}^{(j)}$, and
$\varepsilon_{\mathtt{out}}$ are employed in Algorithm 3. In particular,
$\varepsilon_{\mathtt{out}}$ is used to assess the accuracy of the computed
solution in terms of the Lyapunov residual norm, whereas
$\varepsilon_{\mathtt{inn}}^{(j)}$ is employed to determine whether the
solution of the current linear system is sufficiently correct. In principle,
the user can provide a fixed value for the inner tolerance, i.e.
$\varepsilon_{\mathtt{inn}}^{(j)}\equiv\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\varepsilon}_{\mathtt{inn}}$ for all $j$. However, the theory developed in
[38] can be used to adaptively compute $\varepsilon_{\mathtt{inn}}^{(j)}$ as
the LR-ADI iterations proceed. The _relaxation_ strategy presented in [38,
Section 3] allows us to increase $\varepsilon_{\mathtt{inn}}^{(j)}$ as $j$
grows. Therefore, especially when $\varepsilon_{\mathtt{inn}}^{(j)}$ is rather
large, there is no need to expand the current extended Krylov subspace in
general. In all the results reported in Section 6, we employ such a strategy
and $\varepsilon_{\mathtt{inn}}^{(j)}$ is computed according to [38, Equation
(3.18b)]. See also [44] for similar results in case of Sylvester equations.

We would like to point out that the lines 3 to 3 in Algorithm 3 and the use of
the flag flag_noexpand are crucial to reduce the computational cost of the
overall procedure. Indeed, those lines are devoted to check whether the
current subspace already contains enough spectral information to solve the
current linear system. If this is the case, we do not expand the current space
avoiding unnecessary increments in the memory requirements and computational
efforts.

If $\mathbf{Y}:=[[Y_{m_{1}};O_{2(m_{j}-m_{1})q\times
q}],[Y_{m_{2}};O_{2(m_{2}-m_{1})q\times q}],\ldots,Y_{m_{j}}]$, (19) shows
that the numerical solution computed by the proposed LR-ADI implementation is
of the form

$Z_{j}Z_{j}^{\mkern-1.5mu\mathsf{T}}=-2V_{m_{j}}(\mathbf{Y}(\text{diag}(\text{Re}(p_{1}),\ldots,\text{Re}(p_{j}))\otimes
I_{q})\mathbf{Y}^{\mkern-1.5mu\mathsf{T}})V_{m_{j}}^{\mkern-1.5mu\mathsf{T}}.$
(20)

The right-hand side in (20) has the typical form of an approximate solution
computed by a projection method applied to (2). In particular, if the extended
Krylov subspace method (K-PIK) presented in [63] is applied to  solve (2), the
computed approximation is of the form
$X_{m}=V_{m}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>where the orthonormal columns of $V_{m}$ are a basis of
$\mathbf{EK}_{m}^{\square}(A,B)$ and $L_{m}$ is computed by imposing a
Galerkin condition on the residual matrix
$AV_{m}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}L_{m}}V_{m}^{\mkern-1.5mu\mathsf{T}}+V_{m}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>Therefore, the proposed LR-ADI implementation can be seen as a novel
projection method where the coefficients of the linear combination in terms of
the basis vectors that provides the approximate solution, namely the matrix
$\mathbf{Y}(\text{diag}(\text{Re}(p_{1}),\ldots,\text{Re}(p_{j}))\otimes
I_{q})\mathbf{Y}^{\mkern-1.5mu\mathsf{T}}$, is computed as outlined above and
not by imposing a Galerkin condition on the residual matrix. This perspective
may provide new insights on the relation between LR-ADI and K-PIK. However,
this is beyond the scope of this paper. Similar investigations, relating LR-
ADI and rational Krylov subspace methods have been reported in [25, 74, 75].

The expression (20) resembles the $LDL^{\mkern-1.5mu\mathsf{T}}$-form of the
LR-ADI solution. This formulation, while being more natural for projection-
based solvers, also turned out to be advantageous when LR-ADI is employed as
linear solver for differential matrix equations; see [39].

##  5 Shift computation

Many of the procedures, available in the literature, for the ADI shift
computation need the explicit construction of a basis of
$\operatorname*{Range\,}(Z_{j})$ or a subspace thereof. For instance, in [12]
the authors suggest to use, as shifts $p_{j}$, a subset of the Ritz values of
$A$ with respect to
$\mathcal{Z}_{j}=\operatorname*{Range\,}(\widetilde{Z}_{j})$, where
$\widetilde{Z}_{j}\in\mathbb{R}^{n\times h}$ consists of the last $h>0$
columns of $Z_{j}$ that have been orthogonalized with respect to each other.
However, (19) shows that Algorithm 3 provides us with a matrix $Z_{j}$ such
that
$\operatorname*{Range\,}(Z_{j})\subseteq\mathbf{EK}_{m_{j}}^{\square}(A,B)$ so
that the Ritz values of $A$ with respect to
$\mathbf{EK}_{m_{j}}^{\square}(A,B)$ can be employed as shifts. Moreover, in
standard LR-ADI implementations, one has to explicitly compute the projection
of $A$ onto $\mathcal{Z}_{j}$ increasing the computational efforts of the
overall procedure. In our approach, the projection of $A$ onto
$\mathbf{EK}_{m_{j}}^{\square}(A,B)$ is given for free as this amounts to
$T_{m_{j}}$ and no additional operations are required.

The observation above can be applied to many schemes for the shift
computation. In the following we give some details for the residual-
Hamiltonian-based shifts and the residual norm-minimizing shifts presented in
[37].

In [37, Section 2.1.3], at the $j$-th LR-ADI iteration, the Hamiltonian matrix
$\mathcal{H}_{j}=\begin{bmatrix}A^{\mkern-1.5mu\mathsf{T}}&O\\\
W_{j}W_{j}^{\mkern-1.5mu\mathsf{T}}&-A\\\ \end{bmatrix}$ is considered and its
projection onto $\mathcal{Z}_{j}$, namely
$\mathcal{\widetilde{H}}_{j}=\begin{bmatrix}{(\widetilde{Z}_{j}^{\mkern-1.5mu\mathsf{T}}A\widetilde{Z}_{j})}^{\mkern-1.5mu\mathsf{T}}&O\\\
\widetilde{Z}_{j}^{\mkern-1.5mu\mathsf{T}}W_{j}W_{j}^{\mkern-1.5mu\mathsf{T}}\widetilde{Z}_{j}&-\widetilde{Z}_{j}^{\mkern-1.5mu\mathsf{T}}A\widetilde{Z}_{j}\\\
\end{bmatrix}$, is constructed. In our case, we can easily construct the
projection of $\mathcal{H}_{j}$ onto $\mathbf{EK}_{m_{j}}^{\square}(A,B)$ and
this is given by

$\mathcal{\widetilde{H}}_{j}=\begin{bmatrix}T_{m_{j}}^{\mkern-1.5mu\mathsf{T}}&O\\\
\Upsilon_{j}^{\mkern-1.5mu\mathsf{T}}\Upsilon_{j}&-T_{m_{j}}\\\
\end{bmatrix}\in\mathbb{R}^{4m_{j}q\times 4m_{j}q}.$ (21)

With (21) at hand, we compute its stable eigenpairs
$\left(\lambda_{k},\begin{bmatrix}s_{k}\\\ t_{k}\\\ \end{bmatrix}\right)$,
$\text{Re}(\lambda_{k})<0$, $s_{k},t_{k}\in\mathbb{R}^{2m_{j}q}$, and  the
$(j+1)$-st residual-Hamiltonian-based shift $p_{j+1}$ is selected as the
eigenvalue $\lambda_{\widehat{k}}$ such that
$t_{\widehat{k}}=\operatorname*{argmax}\\{\|t_{k}\|\\}$.

For the computation of residual-norm-minimizing shifts, in [37, Section 3] a
rather involved optimization procedure is presented. In particular, the real
and imaginary parts of $p_{j+1}=\theta_{j+1}+\imath\xi_{j+1}$ are computed by
solving the following minimization problem

$[\theta_{j+1},\xi_{j+1}]=\operatorname*{argmin}_{\theta\in\mathbb{R}_{-},\xi\in\mathbb{R}}\|W_{j}-2\theta({(A+(\theta+\imath\xi)I)}^{-1})W_{j}\|^{2}.$
(22)

The objective function in (22) is expensive to evaluate, making the shift
computation often more expensive than a single LR-ADI iteration. To overcome
this issue, Kürschner proposes to employ smaller matrices $\widetilde{A}$ and
$\widetilde{W}_{j}$ in place of $A$ and $W_{j}$. Once again, $\widetilde{A}$
and $\widetilde{W}_{j}$ are the projection of $A$ and $W_{j}$ onto a suitable
subspace. This subspace is chosen to be
$\mathbf{EK}_{\ell}^{\square}(A,B)\cup\operatorname*{Range\,}(Z_{j})$ for a
certain, usually small, $\ell>0$. In our implementation,
$\mathbf{EK}_{\ell}^{\square}(A,B)\cup\operatorname*{Range\,}(Z_{j})\subseteq\mathbf{EK}_{m_{j}}^{\square}(A,B)$
if $\ell\leqslant m_{j}$. Therefore, we can set $\widetilde{A}=T_{m_{j}}$ and
$\widetilde{W}_{j}=\Upsilon_{j}$ for the approximation of
$[\theta_{j+1},\xi_{j+1}]$ (22).

##  6 Numerical examples

In this section we illustrate the potential of the scheme we propose in this
paper. The two variants of the LR-ADI-EKSM method, we have illustrated in
Section 3, will be denoted by LR-ADI-EKSM(G) and LR-ADI-EKSM(MR). In
particular, in LR-ADI-EKSM(G) we solve the linear systems by imposing a
Galerkin condition, i.e., the matrix $Y$ is computed by solving the reduced
problem (14). In LR-ADI-EKSM(MR), $Y$ solves the least squares problem (16).

We test Algorithm 3 on different instances of (2) coming from the
discretization of certain PDEs, and we study how the computational cost of the
main steps of Algorithm 3 depends on the problem dimension $n$ and rank of the
right hand side $q$.

The results achieved by Algorithm 3 are also compared to the ones obtained by
running a standard implementation of the LR-ADI method. In particular, we
employed the MATLAB® function mess_lradi available in the M-M.E.S.S. package
[59]. Notice that mess_lradi is intended to be a _black-box_ routine so that
many checks and inspections are performed before the actual solution process
starts. This may increase the overall running time of mess_lradi. Therefore,
to have fair comparisons, we also report the results obtained by running a
standard implementation of LR-ADI where the overhead cost mentioned above is
not present. Such a routine is simply denoted by lradi in the tables that
follow.

For a better understanding, in Table 1 we summarize the adopted linear system
solver included in the tested routines for each of the numerical experiments
that follow. Similarly, in Table 1 we indicate whether a given scheme is
equipped with the relaxation strategy coming from [38] for the selection of
$\varepsilon_{\mathtt{inn}}$.

|  | Solver | Relaxation
---|---|---|---
Experiment 1 | LR-ADI-EKSM(G) | backslash | ✓
lradi | backslash | ✗
mess_lradi | backslash | ✗
Experiment 2 | LR-ADI-EKSM(MR) | PGMRES | ✓
lradi | PGMRES | ✓
mess_lradi | PGMRES | ✗
Experiment 3 | LR-ADI-EKSM(G) | backslash | ✓
mess_lradi | backslash | ✗
K-PIK | backslash | —
Experiment 4 | LR-ADI-EKSM(G) | backslash | ✓
mess_lradi | backslash | ✗
Table 1: Solver: solver employed for solving the linear systems with $A$ in
LR-ADI-EKSM and K-PIK and with $A+p_{j}I$ in lradi and mess_lradi. In the
column Relaxation we indicate whether a certain scheme is equipped with the
relaxation strategy proposed in [38].

For all experiments, the tolerance $\varepsilon_{\mathtt{out}}$ for the
relative residual norm is set to $10^{-8}$. Moreover, except for Experiment
6.3, we always employ the residual-Hamiltonian-based shifts presented in [37]
and computed as illustrated in Section 5.

All results were obtained by running MATLAB R2020b [47] on a standard
node*‡*‡*‡CPU: 2x Intel Xeon Skylake Silver 4110 @ 2.1 GHz, 8 cores per CPU.
RAM: 192 GB DDR4 ECC. of the Linux cluster mechthild hosted at the Max Planck
Institute for Dynamics of Complex Technical Systems in Magdeburg, Germany.

######  Example 6.1.

In the first experiment we consider a Lyapunov equation where

$A=I_{h}\otimes D_{h}+D_{h}\otimes I_{h},\quad
D_{h}=\text{tridiag}(1,-2,1)\in\mathbb{R}^{h\times h}.$

Therefore, $A\in\mathbb{R}^{n\times n}$, $n=h^{2}$, is symmetric and stable.
We first consider a matrix $B\in\mathbb{R}^{n\times q}$ with random entries
and unit norm, and in Table 2 we depict how the overall solution time
distributes among the main steps of our algorithm for different values of $n$
and $q$.

In both LR-ADI-EKSM(G) and LR-ADI-EKSM(MR), the linear systems with $A$
required for the basis construction are solved by means of the MATLAB sparse
direct solver “backslash”. In particular, $A$ is factorized once and for all
before the iterative procedures start so that only triangular systems are
actually solved during the basis construction. The computational time for the
factorization of $A$ is always included in the results that follow.

In this experiment, LR-ADI-EKSM(G) and LR-ADI-EKSM(MR) perform very similarly.
We thus report only the results achieved by the former.

|  |  | Basis | Projected Pr. | Shift | Etc | Total
---|---|---|---|---|---|---|---
$n$ | $q$ | It. | Time (s) | Time (s) | Time (s) | Time (s) | Time (s)
360 000 | 1 | 28 | 8.00 | 0.07 | 1.84 | 2.33 | 12.24
| 3 | 27 | 21.70 | 0.26 | 10.35 | 3.66 | 35.97
| 5 | 27 | 32.79 | 0.50 | 22.38 | 4.60 | 60.27
640 000 | 1 | 29 | 16.87 | 0.09 | 2.59 | 4.74 | 24.29
| 3 | 31 | 42.74 | 0.31 | 14.91 | 7.45 | 65.41
| 5 | 31 | 101.15 | 0.81 | 34.70 | 9.26 | 145.92
1 000 000 | 1 | 26 | 28.99 | 0.09 | 2.64 | 7.52 | 39.24
| 3 | 31 | 112.60 | 0.45 | 19.04 | 11.67 | 143.76
| 5 | 29 | 183.80 | 1.12 | 40.56 | 14.52 | 240.00
Table 2: Experiment 6.1. Computational timings devoted to the different main
steps of LR-ADI-EKSM(G) for different values of the problem size $n$ and rank
of the right-hand side $q$. Basis: Basis construction (Algorithm 3 — lines 3,
3, and 3). Projected Pr.: Computation of $Y$ (Algorithm 3 — lines 3, and 3).
Shift: Shifts computation (Algorithm 3 — line 3). Etc: Remainder of the
algorithm (e.g., factorization of $A$). It. indicates the number of ADI
iterations that have been implicitly performed.

As expected, the time devoted to the basis construction represents the
majority of the overall computational efforts. This is the usual case in
Krylov projection algorithms. This cost increases as $q$ grows. Indeed, a
larger subspace is computed making the basis construction, and in particular
the orthogonalization step, rather demanding. Having a large dimensional
approximation space leads to a more expensive shift computation, as well.

In Figure 1 (left $y$-axis) we illustrate how the dimension of the computed
extended Krylov subspace grows in terms of $j$ for $n=360\,000$ and different
values of $q$.

$0$$5$$10$$15$$20$$25$$30$$0$$100$$200$$300$$400$$500$$600$ADI index $j$Space
dimension
$2qm_{j}$$q=1$$q=3$$q=5$$10^{-9}$$10^{-8}$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$Normalized
residual norm $\frac{r_{j}}{r_{0}}$ Figure 1: Experiment 6.1. Dimension of the
constructed extended Krylov subspace and  computed normalized residual norms
as $j$ grows, i.e. the ADI progresses, for problem size $n=360\,000$.

In this experiment, we can notice that the subspace constructed to solve the
second shifted linear system, namely $(A+p_{2}I)S_{2}=W_{1}$, is a very rich
approximation space in terms of spectral information. Indeed, we need to only
slightly expand it to solve the subsequent linear systems  without
compromising the decrease in the Lyapunov residual norm; see Figure 1 (right
$y$-axis). This means that the majority of the computational efforts are
dedicated to solve the second linear system, and we can capitalize on them for
$j>2$ reducing the overall workload of the solution process. We would like to
mention that such a phenomenon is partially due to the adaptive selection of
the inner tolerance $\varepsilon_{\mathtt{inn}}^{(j)}$ coming from [38].

We now compare LR-ADI-EKSM(G) with the function mess_lradi of the M-M.E.S.S.
package [59], an abstract function handle based implementation of the LR-ADI,
and lradi, a plain matrix-based implementation of the same algorithm.

To this end, we make $B\in\mathbb{R}^{n}$ the normalized vector of all ones.
For having fair comparisons, we employ the shifts computed by the LR-ADI-
EKSM(G) in all the different implementations. This leads to a very similar
trend in the relative residual norm achieved by the routines even though the
shifted linear systems in mess_lradi and lradi are solved at very high
accuracy*§*§*§The MATLAB sparse direct solver “backslash” is employed for
solving $(A+p_{j}I)S_{j}=W_{j-1}$ for all $j$., whereas the relaxation
strategy of [38] is implemented in LR-ADI-EKSM(G). In Figure 2, we report the
relative difference between the relative residual norms computed by LR-ADI-
EKSM(G) and mess_lradi throughout all the necessary iterations $j$ for
different problem dimension $n$ along with the values of
$\epsilon^{(j)}_{\mathtt{inn}}$ we employed. In agreement with the results
presented in [38], we can notice that  the distance between the computed
relative residual norms is always rather moderate and smaller than
$\epsilon^{(j)}_{\mathtt{inn}}$*¶*¶*¶This is true for all the experiments we
ran except for $n=640\,000$, at the very last iteration where
$r_{27}^{\text{LR-ADI-EKSM}}\approx 1.3\times 10^{-8}$ whereas
$\epsilon_{\mathtt{inn}}^{(27)}\approx 1.6\times 10^{-8}$ . Very similar
results are obtained by comparing the residual norms attained by lradi in
place of mess_lradi.

figure2.tikzexternalize/figure2.pdf

Figure 2: Experiment 6.1. Relative gap $\left(\frac{|r_{j}^{\text{LR-ADI-
EKSM}}-r_{j}^{\texttt{mess\\_lradi}}|}{r_{j}^{\texttt{mess\\_lradi}}}\right)$
between the residual norms $r_{j}^{\text{LR-ADI-EKSM}}$ and
$r_{j}^{\texttt{mess\\_lradi}}$ computed by LR-ADI-EKSM(G) and mess_lradi,
respectively, as $j$ grows, i.e. ADI converges, and different problem sizes
$n$, together with the corresponding inner inexact solver tolerance
$\varepsilon^{(j)}_{\mathtt{inn}}$, denoted $\varepsilon_{n}$ to relate the
problem sizes.

We also compare the routines in terms of computation time. The results are
collected in Table 3. Since we employ the shifts computed within LR-ADI-
EKSM(G) also for mess_lradi and lradi, we do not consider the time devoted to
the shift computation when reporting the performances of LR-ADI-EKSM(G) in
Table 3.

|  | LR-ADI-EKSM(G) | lradi | mess_lradi
---|---|---|---|---
$n$ | It. | Time (s) | Time (s) | Time (s)
360 000 | 24 | 10.08 | 31.84 | 30.94
640 000 | 27 | 19.99 | 65.84 | 66.22
1 000 000 | 36 | 39.18 | 101.18 | 101.61
Table 3: Experiment 6.1. Computational timings achieved by LR-ADI-EKSM(G),
lradi, and mess_lradi for different problem sizes $n$. It. indicates the
number of ADI iterations that have been (implicitly) performed.

The results in Table 3 show that, for this experiment, our proposed scheme
combined with the relaxation strategy presented in [38] leads to a remarkable
speed-up of the solution process — up to 50% — when compared to a standard
implementation of the LR-ADI method.

######  Example 6.2.

In the second experiment, we consider a problem similar to [51, Example 6]. In
particular, the matrix $A$ comes from the centered finite difference
discretization of the 3-dimensional convection-diffusion operator

$\mathcal{L}(u)=-\zeta\Delta u+\mathbf{w}\cdot\nabla u,$

on the unit cube with zero Dirichlet boundary conditions. The convection
vector $\mathbf{w}$ is given by
$\mathbf{w}=(\phi_{1}(x)\psi_{1}(y)\pi_{1}(z),\,0,\,\pi_{3}(z))=((1-x^{2})yz,\,0,\,e^{z})$
whereas $\zeta>0$. By employing $h$ nodes in each direction, the
discretization phase leads to a matrix $A$ that can be written as

$A=(D_{h}+\Pi_{3}N^{\mkern-1.5mu\mathsf{T}})\otimes I_{h}\otimes
I_{h}+I_{h}\otimes D_{h}\otimes I_{h}+I_{h}\otimes I_{h}\otimes
D_{h}+\Pi_{1}\otimes\Psi_{1}\otimes\Phi_{1}N,$

where
$D_{h}=\zeta{(h-1)}^{2}\cdot\text{tridiag}(-1,2,-1)\in\mathbb{R}^{h\times h}$,
$N=-\frac{(h-1)}{2}\cdot\text{tridiag}(-1,0,1)\in\mathbb{R}^{h\times h}$, and
$\Phi_{i}$, $\Psi_{i}$, and $\Pi_{i}$ are diagonal matrices whose diagonal
entries correspond to the nodal values of the corresponding functions
$\phi_{i}$, $\psi_{i}$, and $\pi_{i}$. See [51] for further details.
$B\in\mathbb{R}^{n}$, $n=h^{3}$, is a vector with random entries.

Due to the 3D nature of the problem, the nonsymmetric linear systems with $A$
involved in the basis construction in LR-ADI-EKSM are solved by GMRES [57]. In
particular, we employ the GMRES implementation written by Lund et al [35],
namely the function bgmres in [45]. GMRES is stopped whenever the computed
relative residual norm gets smaller than $10^{-10}$.

It is well-known that (polynomial) Krylov methods for linear systems need to
be preconditioned to achieve a fast convergence in terms of number of
iterations. To this end, as suggested in [51], we employ the following
preconditioning operator when solving the linear systems with $A$,

$\mathcal{P}=(D_{h}+\Pi_{3}N^{\mkern-1.5mu\mathsf{T}})\otimes I_{h}\otimes
I_{h}+I_{h}\otimes D_{h}\otimes I_{h}+I_{h}\otimes I_{h}\otimes
D_{h}+\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\pi}_{1}I_{h}\otimes\Psi_{1}\otimes\Phi_{1}N,$

where $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\pi}_{1}$ is the mean value of the function $\pi_{1}$ in $[0,1]$. At each
GMRES iteration, we thus have to invert $\mathcal{P}$, namely we have to
compute $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{v}=\mathcal{P}^{-1}v$ for $v\in\mathbb{R}^{n}$. This operation is
performed by solving the Sylvester equation

$(D_{h}\otimes I_{h}+I_{h}\otimes D_{h}+\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\pi}_{1}\Psi_{1}\otimes\Phi_{1}N)\mathbf{\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{V}}+\mathbf{\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{V}}{(D_{h}+\Pi_{3}N^{\mkern-1.5mu\mathsf{T}})}^{\mkern-1.5mu\mathsf{T}}=\mathbf{V},$

where $\mathbf{\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{V}},\mathbf{V}\in\mathbb{R}^{h^{2}\times h}$ are such that
$\text{vec}(\mathbf{\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{V}})=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{v}$
and $\text{vec}(\mathbf{V})=v$. Since the coefficient matrices in the equation
above have moderate dimensions, the Bartels-Stewart method [6] is employed for
its solution and the Schur decompositions of the coefficient matrices are
computed once and for all before the iterative procedure starts. We always
employ a right preconditioning scheme in order to easily have access to the
actual residual norm.

Also, for the shifted linear systems with $A+p_{j}I$, within mess_lradi and
lradi, we employ preconditioned GMRES equipped with the preconditioning
operator $\mathcal{P}+p_{j}I$. Once again, this preconditioner is applied by
solving the Sylvester equation

$(D_{h}\otimes I_{h}+I_{h}\otimes D_{h}+\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\pi}_{1}\Psi_{1}\otimes\Phi_{1}N)\mathbf{\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{V}}+\mathbf{\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{V}}{(D_{h}+\Pi_{3}N^{\mkern-1.5mu\mathsf{T}}+p_{j}I_{h})}^{\mkern-1.5mu\mathsf{T}}=\mathbf{V}.$

Even though this is in general a better preconditioner for $A+p_{j}I$ compared
to $\mathcal{P}$, its application involves complex arithmetic whenever
$\text{Im}(p_{j})\neq 0$ with a consequent increment in the computational
efforts devoted to the preconditioning step.

For this experiment, lradi is equipped with the relaxation strategy presented
in [38].

Also for this experiment, LR-ADI-EKSM(G) and LR-ADI-EKSM(MR) perform very
similarly, with LR-ADI-EKSM(MR) achieving slightly better results in terms of
computational time. We thus report only the performance of LR-ADI-EKSM(MR).

The results are collected in Table 4 for different values of $n$ and $\zeta$.
In Table 4 we also report the number of shifts with nonzero imaginary part.

We would like to mention that we ran some experiments with mess_lradi where
the shifted linear systems were solved by means of the MATLAB sparse direct
solver “backslash” in place of preconditioned GMRES. However, for this example
the potentially higher accuracy of the direct solves did not benefit the
computation and the execution times we achieved with “backslash” could not
keep up with the ones reported for GMRES in Table 4. We, thus, decided to omit
them here.

|  |  |  | LR-ADI-EKSM(MR) | lradi | mess_lradi
---|---|---|---|---|---|---
$\zeta$ | $n$ | It. | #$\\{p_{j}\notin\mathbb{R}\\}$ | Time (s) | Time (s) | Time (s)
0.05 | 125 000 | 20 | 12 | 80.61 | 153.28 | 187.52
| 512 000 | 20 | 12 | 812.74 | 1 342.27 | 2 161.99
| 1 000 000 | 22 | 10 | 3 183.54 | 5 764.21 | 6 133.54
0.005 | 125 000 | 45 | 44 | 65.23 | 361.58 | 382.11
| 512 000 | 61 | 60 | 419.45 | 2 156.94 | 3 497.51
| 1 000 000 | 67 | 62 | 1 194.37 | 5 802.41 | 10 517.13
Table 4: Experiment 6.2. Computational timings achieved by LR-ADI-EKSM(MR),
lradi, and mess_lradi for different problem sizes $n$ and diffusivities
$\zeta$. It. indicates the number of ADI iterations that have been
(implicitly) performed.

From the results in Table 4 we can see that LR-ADI-EKSM(MR) is very
competitive and always achieves computational timings that are significantly
smaller than the ones required by mess_lradi. Thanks to the relaxation
procedure coming from [38], lradi performs better than mess_lradi.

The performance of all the tested routines is strictly related to the number
of complex shifts needed to converge. When this is sizable with respect to the
total number of iterations, many of the $n\times n$ linear systems $A+p_{j}I$
within mess_lradi and lradi involve complex arithmetic, whereas this is needed
only in the solution of the small dimensional least squares problem for the
computation of $Y$ in LR-ADI-EKSM(MR).

We notice that, for a fixed $n$, the computational time of LR-ADI-EKSM(MR)
decreases, in general, by reducing $\zeta$, even tough the number of LR-ADI
iterations that are implicitly performed increases. This is due to the
computational efforts required by the solution of the linear systems with $A$
during the basis construction. Indeed, for $\zeta=0.05$, many more GMRES
iterations are required than what is necessary for $\zeta=0.005$. In Figure 3,
we report the number of GMRES iterations needed to solve the linear system
with $A$ at each $m$, namely every time a new basis vector of the adopted
extended Krylov subspace needs to be computed.

figure3.tikzexternalize/figure3.pdf

Figure 3: Experiment 6.2. Number of GMRES iterations needed to solve the
linear systems with $A$ during the basis construction in LR-ADI-EKSM(MR) for
different values of the diffusivity $\zeta$ and $n=125\,000$.

A rather large number of GMRES iterations is required for solving the linear
systems with $A$ in case of $\zeta=0.05$ making the construction of the basis
of $\mathbf{EK}^{\square}_{m}(A,B)$ more demanding. On the other hand, few
GMRES iterations are sufficient to meet the prescribed accuracy for
$\zeta=0.005$ and the overall solution procedure turns out to be very
successful.

######  Example 6.3.

In this experiment we compare LR-ADI-EKSM also with K-PIK [63], since the two
routines construct the same subspace*∥*∥*∥The implementation of K-PIK we
employed will be available in the next M-M.E.S.S. release, along with other
projection methods for matrix equations. Such implementation is equivalent to
the one that can be found on Simoncini’s webpage,
http://www.dm.unibo.it/~simoncin/software.html.. We consider the thermal part
of the thermo-elastic modeling of a building-block of an experimental machine
tool given by the following heat equation

$\left\\{\begin{array}[]{rll}c_{p}\rho\frac{\partial T}{\partial
t}&=&\lambda\Delta T,\quad\text{in }\Omega,\\\ \lambda\frac{\partial
T}{\partial\textbf{n}}&=&f,\quad\text{on }\Gamma_{c}\subset\partial\Omega,\\\
\lambda\frac{\partial T}{\partial\textbf{n}}&=&\alpha(T_{ext}-T),\quad\text{on
}\Gamma_{ext}\subset\partial\Omega,\\\ T(0)&=&0.\\\ \end{array}\right.$ (23)

Figure 4: Experiment 6.3. Finite element grid of the machine frame indicated
on the CAD model of the full machine. (Source: DFG CRC/TR-96
(https://transregio96.de))

The discretization in space using the finite element method (here applying the
proprietary tool ANSYS*********https://www.ansys.com/) on the three-
dimensional domain, given by the machine frame indicated in Figure 4, leads to
the LTI system

$E\dot{T}=\left(A-\sum_{i=1}^{t}\alpha_{i}F_{i}\right)T+Bu(t).$ (24)

Here, $A$ represents the discretized Laplacian together with the Robin
boundary contributions from $\Gamma_{ext}$ and represented by $F_{i}$, while
$B$ results from the external control inputs (heats fluxes, e.g. induced by
the drive motors) on $\Gamma_{c}$. Note that the elastic part of the thermo-
elastic model can be encoded entirely in the output equation of the
corresponding dynamical system and is, thus, not relevant here [40]. The
algebraic problem resulting from this system amounts to a Lyapunov equation of
the form (1). However, due to mass lumping in ANSYS, the mass matrix $E$ is
diagonal and SPD. We can, thus, easily invert its square root and consider the
Lyapunov equation

$E^{-\frac{1}{2}}\left(A-\sum_{i=1}^{t}\alpha_{i}F_{i}\right)E^{-\frac{1}{2}}\widetilde{X}+\widetilde{X}E^{-\frac{1}{2}}{\left(A-\sum_{i=1}^{t}\alpha_{i}F_{i}\right)}^{\mkern-1.5mu\mathsf{T}}E^{-\frac{1}{2}}+E^{-\frac{1}{2}}BB^{\mkern-1.5mu\mathsf{T}}E^{-\frac{1}{2}}=0,\quad\widetilde{X}=(E^{\frac{1}{2}}XE^{\frac{1}{2}}).$

So, again, we can efficiently retract to a problem of the form (2). Once a
low-rank approximation $\widetilde{Z}\widetilde{Z}^{\mkern-1.5mu\mathsf{T}}$
to $\widetilde{X}$ is computed, the low-rank factor $Z$ such that
$ZZ^{\mkern-1.5mu\mathsf{T}}\approx X$ can be retrieved by performing
$Z=E^{-\frac{1}{2}}\widetilde{Z}$.

The actual machine frame in Figure 4 consists of several parts itself, which
are discretized separately. This leads to differently sized models of the
structure in (24). These are reflected by the rows of Table 5. Accordingly, we
solve the Lyapunov equation considering different configurations of the PDE
(23), respectively the LTI system in (24). In particular, this allows us to
vary the number of degrees of freedom employed in the discretization phase,
leading to different problem dimensions $n$, modify the Neumann boundary
conditions obtaining diverse matrices $F_{i}$, and consider different values
for the rank $q$ of $B$. Moreover, we set $\alpha_{i}=10$ for all
$i=1,\ldots,t$.

The results are collected in Table 5. It turns out that the Wachspress ADI
shifts [73, 42] are particularly effective for this experiment, since $A$ as
well as all the $F_{i}$ and thus
$E^{-\frac{1}{2}}\left(A-\sum_{i=1}^{t}\alpha_{i}F_{i}\right)E^{-\frac{1}{2}}$
are symmetric, i.e. the spectrum is real. These are the ideal circumstances
for Wachspress shifts. We, thus, employ those shifts in LR-ADI-EKSM(G) and
mess_lradi.

|  |  | LR-ADI-EKSM(G) | mess_lradi | K-PIK
---|---|---|---|---|---
$n$ | $t$ | $q$ | It. | $\text{dim}\left(\mathbf{EK}^{\square}\right)$ | Time (s) | It. | Time (s) | It. | $\text{dim}\left(\mathbf{EK}^{\square}\right)$ | Time (s)
4 813 | 1 | 23 | 54 | 644 | 2.73 | 54 | 7.93 | 15 | 736 | 8.13
13 551 | 2 | 5 | 53 | 430 | 8.89 | 53 | 18.09 | 43 | 440 | 24.93
25 872 | 1 | 10 | 63 | 1060 | 34.32 | 63 | 63.37 | 53 | 1080 | 97.74
Table 5: Experiment 6.3. Computational timings achieved by LR-ADI-EKSM(G),
K-PIK, and mess_lradi for different values of problem size $n$, number of
Robin boundary conditions $t$, and rank of the right-hand side $q$. It.
indicates the number of ADI/K-PIK iterations that have been (implicitly)
performed.

For this experiment, the LR-ADI method, either based on our new formulation or
on a standard scheme as the one in mess_lradi, turns out to be more efficient
in terms of computational time than K-PIK. Indeed, in spite of the smaller
number of iterations needed to converge, the large dimension of the extended
Krylov subspace constructed by K-PIK leads to a rather costly solution of the
projected equations. Also LR-ADI-EKSM(G) requires the construction of an
extended Krylov subspace whose dimension is similar to the one computed by
K-PIK. However, if
$\text{dim}\left(\mathbf{EK}^{\square}_{m}(E^{-\frac{1}{2}}\left(A-\sum_{i=1}^{t}\alpha_{i}F_{i}\right)E^{-\frac{1}{2}},E^{-\frac{1}{2}}B)\right)=2mq$,
the computational cost of solving the inner problems within LR-ADI-EKSM(G) is
$\mathcal{O}(4m^{2}q^{2})$ floating-point operations (FLOPs) whereas it
amounts to $\mathcal{O}(8m^{3}q^{3})$ FLOPs for K-PIK.

We conclude by mentioning that in this experiment we relied on the ease of
computing $E^{-\frac{1}{2}}$. However, it may happen that the mass matrix $E$
cannot be easily manipulated, e.g., it can be possibly singular, so that the
routine presented in this paper cannot be readily applied as we have done in
this experiment. We plan to extend the LR-ADI-EKSM framework to this more
challenging class of equations in the near future.

######  Example 6.4.

In the last experiment, we show that the proposed framework still needs some
further improvements to efficiently deal with generalized Lyapunov equations
of the form (1) where the mass matrix $E$ is not diagonal. To this end, we
consider the Steel Profile data set [49, 18] from the MORwiki repository [72].

We compute the observability Gramian of the system, namely the solution $X$ to
the equation

$A^{\mkern-1.5mu\mathsf{T}}XE+E^{\mkern-1.5mu\mathsf{T}}XA+C^{\mkern-1.5mu\mathsf{T}}C=0,$
(25)

where $A\in\mathbb{R}^{n\times n}$ is symmetric negative definite,
$C\in\mathbb{R}^{q\times n}$, $q=6$, and $E\in\mathbb{R}^{n\times n}$ is SPD
but not diagonal. See [17] fur further details on the model.

If $E=LL^{\mkern-1.5mu\mathsf{T}}$ denotes the Cholesky factorization of $E$,
we consider the transformed equation

$(L^{-1}A^{\mkern-1.5mu\mathsf{T}}L^{-\mkern-1.5mu\mathsf{T}})\widetilde{X}+\widetilde{X}(L^{-1}AL^{-\mkern-1.5mu\mathsf{T}})+L^{-1}C^{\mkern-1.5mu\mathsf{T}}CL^{-\mkern-1.5mu\mathsf{T}}=0,\quad\widetilde{X}=L^{\mkern-1.5mu\mathsf{T}}XL,$
(26)

and, due to symmetry of $A$, employ the extended Krylov subspace
$\mathbf{EK}_{m}^{\square}(L^{-1}AL^{-\mkern-1.5mu\mathsf{T}},L^{-1}C^{\mkern-1.5mu\mathsf{T}})$
as approximation space. Notice that the matrix
$L^{-1}AL^{-\mkern-1.5mu\mathsf{T}}$ does not need to be explicitly
constructed. See, e.g., [63, Example 5.4]. As before, once
$\widetilde{Z}\widetilde{Z}^{\mkern-1.5mu\mathsf{T}}\approx\widetilde{X}$ is
computed, we obtain a low-rank approximation to the original $X$ by performing
$Z=L^{-\mkern-1.5mu\mathsf{T}}\widetilde{Z}$.

In Table 6 we report the results achieved by LR-ADI-EKSM(G) and mess_lradi for
different values of $n$.

| LR-ADI-EKSM(G) | mess_lradi
---|---|---
$n$ | It. | $\text{dim}\left(\mathbf{EK}^{\square}\right)$ | $\text{rank}(X)$ | Time (s) | It. | $\text{rank}(X)$ | Time (s)
20 209 | 30 | 564 | 180 | 7.07 | 30 | 180 | 0.54
79 841 | 31 | 816 | 186 | 34.09 | 31 | 186 | 2.99
Table 6: Experiment 6.4. Computational timings achieved by LR-ADI-EKSM(G) and
mess_lradi for different values of the problem size $n$. The running time
devoted to the shift computation is not included. It. indicates the number of
ADI iterations that have been (implicitly) performed.

From the results in Table 6 we can readily see that the standard scheme of the
LR-ADI method implemented in mess_lradi is much faster than LR-ADI-EKSM(G).
This is due to the fact that the latter algorithm needs to construct a quite
large subspace to achieve the prescribed accuracy with a consequent increment
in the computational efforts of the overall procedure.

We also mention that the rank of the approximate solution computed by LR-ADI-
EKSM(G) is much lower than the dimension of the constructed subspace. We
believe that the transformation we performed in (26), and thus the employment
of
$\mathbf{EK}_{m}^{\square}(L^{-1}AL^{-\mkern-1.5mu\mathsf{T}},L^{-1}C^{T})$,
may lead to some spectral redundancy in the adopted approximation subspace and
a slower convergence of the method. On the other hand, mess_lradi is able to
deal with the original formulation (25) of the problem.

To address generalized equations of the form (1), we plan to study the
employment of different techniques within the Krylov LR-ADI framework we
presented in this paper. In particular, the use of nonstandard inner products
and (extended) generalized Krylov subspaces [43] will be explored.

##  7 Conclusions

A new formulation of the LR-ADI algorithm for large-scale standard Lyapunov
equations has been proposed. The computational core of the LR-ADI scheme
consists in the solution of a shifted linear system at each iteration. We
showed that the extended Krylov subspace method can be a valid candidate for
this task. In particular, we described how only one extended Krylov subspace
needs to be constructed to solve all the necessary linear systems required by
the LR-ADI method. The LR-ADI iteration has been completely merged into the
extended Krylov subspace method for shifted linear systems resulting in a
novel, efficient solution procedure. We also showed that many state-of-the-art
algorithms for the shift computation can be easily integrated into our new
scheme. Numerical results demonstrate the potential of our novel algorithm,
especially when this is equipped with the relaxation strategy proposed in
[38], and many complex shifts are needed to converge.

In future work we will consider more involved Lyapunov equations of the form
(1) that cannot be easily transformed into (2). While standard implementations
of the LR-ADI method naturally address such a scenario by solving linear
systems of the form $A+p_{j}E$, further care has to be taken to employ the
scheme we presented in this paper. Indeed, the shifted Arnoldi relation (7)
can no longer be exploited. The use of non-standard inner products and
generalized Krylov subspace methods [43] will be investigated.

The framework presented in this paper can be generalized to enhance other LR-
ADI-like algorithms for matrix equations. For instance, the LR-ADI method for
Sylvester equations [14], or LR-RADI schemes for Riccati equations [10, 22]
can be equipped with a procedure similar to the one we proposed here.

## Acknowledgments

The second author is member of the Italian INdAM Research group GNCS.

The work presented in this paper has been carried out when the second author
was affiliated with the Research Group Computational Methods in Systems and
Control Theory (CSC), Max Planck Institute for Dynamics of Complex Technical
Systems, Sandtorstraße 1, 39106 Magdeburg, Germany.

The _Steel Profile_ dataset is available in the MORwiki repository [72]. All
the other datasets generated during and/or analysed during the current study
are available from the corresponding author on reasonable request.

## References

  * [1] K. Ahuja, P. Benner, E. de Sturler, and L. Feng, Recycling BiCGSTAB with an application to parametric model order reduction, SIAM J. Sci. Comput., 37 (2015), pp. S429–S446, https://doi.org/10.1137/140972433.
  * [2] K. Ahuja, E. de Sturler, S. Gugercin, and E. R. Chang, Recycling BiCG with an application to model reduction, SIAM J. Sci. Comput., 34 (2012), pp. A1925–A1949, https://doi.org/10.1137/100801500.
  * [3] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems, vol. 6 of Adv. Des. Control, SIAM Publications, Philadelphia, 2005, https://doi.org/10.1137/1.9780898718713.
  * [4] H. Anzt, E. Chow, J. Saak, and J. Dongarra, Updating incomplete factorization preconditioners for model order reduction, Numer. Algorithms, 73 (2016), pp. 611–630, https://doi.org/10.1007/s11075-016-0110-2.
  * [5] J. Baker, M. Embree, and J. Sabino, Fast singular value decay for Lyapunov solutions with nonnormal coefficients, SIAM J. Matrix Anal. Appl., 36 (2015), pp. 656–668, https://doi.org/10.1137/140993867.
  * [6] R. H. Bartels and G. W. Stewart, Solution of the matrix equation ${AX}+{XB}={C}$: Algorithm 432, Comm. ACM, 15 (1972), pp. 820–826, https://doi.org/10.1145/361573.361582.
  * [7] M. Baumann and M. B. van Gijzen, Nested Krylov methods for shifted linear systems, SIAM J. Sci. Comput., 37 (2015), pp. S90–S112, https://doi.org/10.1137/140979927.
  * [8] M. Baumann and M. B. van Gijzen, Efficient iterative methods for multi-frequency wave propagation problems: A comparison study, Procedia Computer Science, 108 (2017), pp. 645 – 654, https://doi.org/10.1016/j.procs.2017.05.088.  International Conference on Computational Science, ICCS 2017, 12-14 June 2017, Zurich, Switzerland.
  * [9] S. Bellavia, V. De Simone, D. di Serafino, and B. Morini, Efficient preconditioner updates for shifted linear systems, SIAM J. Sci. Comput., 33 (2011), pp. 1785–1809, https://doi.org/10.1137/100803419.
  * [10] P. Benner, Z. Bujanović, P. Kürschner, and J. Saak, RADI: A low-rank ADI-type algorithm for large scale algebraic Riccati equations, Numer. Math., 138 (2018), pp. 301–330, https://doi.org/10.1007/s00211-017-0907-5.
  * [11] P. Benner, P. Kürschner, and J. Saak, Efficient handling of complex shift parameters in the low-rank Cholesky factor ADI method, Numer. Algorithms, 62 (2013), pp. 225–251, https://doi.org/10.1007/s11075-012-9569-7.
  * [12] P. Benner, P. Kürschner, and J. Saak, Self-generating and efficient shift parameters in ADI methods for large Lyapunov and Sylvester equations, Electron. Trans. Numer. Anal., 43 (2014), pp. 142–162, http://etna.mcs.kent.edu/volumes/2011-2020/vol43/abstract.php?vol=43&pages=142-162.
  * [13] P. Benner, P. Kürschner, and J. Saak, Frequency-limited balanced truncation with low-rank approximations, SIAM J. Sci. Comput., 38 (2016), pp. A471–A499, https://doi.org/10.1137/15M1030911.
  * [14] P. Benner, R.-C. Li, and N. Truhar, On the ADI method for Sylvester equations, J. Comput. Appl. Math., 233 (2009), pp. 1035–1045, https://doi.org/10.1016/j.cam.2009.08.108.
  * [15] P. Benner, V. Mehrmann, and D. C. Sorensen, Dimension Reduction of Large-Scale Systems, vol. 45 of Lect. Notes Comput. Sci. Eng., Springer-Verlag, Berlin/Heidelberg, Germany, 2005, https://doi.org/10.1007/3-540-27909-1.
  * [16] P. Benner, H. Mena, and J. Saak, On the parameter selection problem in the Newton-ADI iteration for large-scale Riccati equations, Electron. Trans. Numer. Anal., 29 (2008), pp. 136–149, https://etna.math.kent.edu/volumes/2001-2010/vol29/abstract.php?vol=29&pages=136-149.
  * [17] P. Benner and J. Saak, Linear-quadratic regulator design for optimal cooling of steel profiles, Tech. Report SFB393/05-05, Sonderforschungsbereich 393 Parallele Numerische Simulation für Physik und Kontinuumsmechanik, TU Chemnitz, D-09107 Chemnitz (Germany), 2005, http://nbn-resolving.de/urn:nbn:de:swb:ch1-200601597.
  * [18] P. Benner and J. Saak, A semi-discretized heat transfer model for optimal cooling of steel profiles, in Dimension Reduction of Large-Scale Systems, P. Benner, V. Mehrmann, and D. Sorensen, eds., vol. 45 of Lect. Notes Comput. Sci. Eng., Springer-Verlag, Berlin/Heidelberg, Germany, 2005, pp. 353–356, https://doi.org/10.1007/3-540-27909-1_19.
  * [19] P. Benner and J. Saak, Numerical solution of large and sparse continuous time algebraic matrix Riccati and Lyapunov equations: a state of the art survey, GAMM Mitteilungen, 36 (2013), pp. 32–52, https://doi.org/10.1002/gamm.201310003.
  * [20] M. Benzi and D. Bertaccini, Approximate inverse preconditioning for shifted linear systems, BIT, 43 (2003), pp. 231–244, https://doi.org/10.1023/A:1026089811044.
  * [21] D. Bertaccini, Efficient preconditioning for sequences of parametric complex symmetric linear systems, Electron. Trans. Numer. Anal., 18 (2004), pp. 49–64, https://etna.math.kent.edu/volumes/2001-2010/vol18/abstract.php?vol=18&pages=49-64.
  * [22] C. Bertram and H. Faßbender, Riccati ADI: Existence, uniqueness and new iterative methods, e-print 2004.11212, arXiv, 2020, https://arxiv.org/abs/2004.11212.  math.NA.
  * [23] B. N. Datta and Y. Saad, Arnoldi methods for large Sylvester-like observer matrix equations, and an associated algorithm for partial spectrum assignment, Linear Algebra Appl., 154–156 (1991), pp. 225 – 244, https://doi.org/10.1016/0024-3795(91)90378-A.
  * [24] V. Druskin and L. Knizhnerman, Extended Krylov subspaces: Approximation of the matrix square root and related functions, SIAM J. Matrix Anal. Appl., 19 (1998), pp. 755–771, https://doi.org/10.1137/S0895479895292400.
  * [25] V. Druskin, L. Knizhnerman, and V. Simoncini, Analysis of the rational Krylov subspace and ADI methods for solving the Lyapunov equation, SIAM J. Numer. Anal., 49 (2011), pp. 1875–1898, https://doi.org/10.1137/100813257.
  * [26] L. Feng, P. Benner, and J. G. Korvink, Subspace recycling accelerates the parametric macromodeling of MEMS, Internat. J. Numer. Methods Engrg., 94 (2013), pp. 84–110, https://doi.org/10.1002/nme.4449.
  * [27] A. Feriani, F. Perotti, and V. Simoncini, Iterative system solvers for the frequency analysis of linear mechanical systems, Comp. Meth. Appl. Mech. Eng., 190 (2000), pp. 1719 – 1739, https://doi.org/10.1016/S0045-7825(00)00187-0.
  * [28] R. W. Freund, Solution of shifted linear systems by quasi-minimal residual iterations, in Numerical Linear Algebra, De Gruyter, 1993, pp. 101–122, https://doi.org/10.1515/9783110857658.101.
  * [29] A. Frommer and U. Glässner, Restarted GMRES for shifted linear systems, SIAM J. Sci. Comput., 19 (1998), pp. 15–26, https://doi.org/10.1137/S1064827596304563.
  * [30] A. Frommer, K. Lund, and D. B. Szyld, Block Krylov subspace methods for functions of matrices, Electron. Trans. Numer. Anal., 47 (2017), pp. 100–126, https://etna.math.kent.edu/volumes/2011-2020/vol47/abstract.php?vol=47&pages=100-126.
  * [31] A. Gaul, Recycling Krylov subspace methods for sequences of linear systems – analysis and applications, PhD thesis, Technische Universität Berlin, 2014, https://doi.org/10.13140/2.1.4015.3284.
  * [32] U. Glässner, S. Güsken, T. Lippert, G. Ritzenhöfer, K. Schilling, and A. Frommer, How to compute Green’s functions for entire mass trajectories within Krylov solvers, Int. J. Mod. Phys. C, 7 (1996), pp. 635–644, https://doi.org/10.1142/S0129183196000533.
  * [33] L. Grasedyck, Existence of a low rank or $H$-matrix approximant to the solution of a Sylvester equation, Numer. Lin. Alg. Appl., 11 (2004), pp. 371–389, https://doi.org/10.1002/nla.366.
  * [34] R. E. Kalman, P. L. Falb, and M. A. Arbib, Topics in Mathematical System Theory, McGraw-Hill, New York, 1969.
  * [35] D. Kressner, K. Lund, S. Massei, and D. Palitta, Compress-and-restart block Krylov subspace methods for Sylvester matrix equations, Numer. Lin. Alg. Appl., 28 (2021), https://doi.org/doi.org/10.1002/nla.2339.
  * [36] P. Kürschner, Efficient Low-Rank Solution of Large-Scale Matrix Equations, PhD thesis, Otto-von-Guericke-Universität, Magdeburg, Germany, Apr. 2016, http://hdl.handle.net/11858/00-001M-0000-0029-CE18-2.  Shaker Verlag, ISBN 978-3-8440-4385-3.
  * [37] P. Kürschner, Approximate residual-minimizing shift parameters for the low-rank ADI iteration, Electron. Trans. Numer. Anal., 51 (2019), pp. 240–261, https://doi.org/10.1553/etna_vol51s240.
  * [38] P. Kürschner and M. Freitag, Inexact methods for the low rank solution to large scale Lyapunov equations, BIT, 60 (2020), pp. 1221–1259, https://doi.org/10.1007/s10543-020-00813-4.
  * [39] N. Lang, H. Mena, and J. Saak, On the benefits of the $LDL^{T}$ factorization for large-scale differential matrix equation solvers, Linear Algebra Appl., 480 (2015), pp. 44–71, https://doi.org/10.1016/j.laa.2015.04.006.
  * [40] N. Lang, J. Saak, and P. Benner, Model order reduction for systems with moving loads, at-Automatisierungstechnik, 62 (2014), pp. 512–522, https://doi.org/10.1515/auto-2014-1095.
  * [41] A. J. Laub, Numerical linear algebra aspects of control design computations, IEEE Trans. Autom. Control, 30 (1985), pp. 97–108, https://doi.org/10.1109/TAC.1985.1103900.
  * [42] J.-R. Li and J. White, Low rank solution of Lyapunov equations, SIAM J. Matrix Anal. Appl., 24 (2002), pp. 260–280, https://doi.org/10.1137/S0895479801384937.
  * [43] R.-C. Li and Q. Ye, A Krylov subspace method for quadratic matrix polynomials with application to constrained least squares problems, SIAM J. Matrix Anal. Appl., 25 (2003), pp. 405–428, https://doi.org/10.1137/S0895479802409390.
  * [44] Z. Liu., Y. Zhou., and Y. Zhang., On inexact alternating direction implicit iteration for continuous Sylvester equations, Numer. Lin. Alg. Appl., 27:e2320 (2020), https://doi.org/10.1002/nla.2320.
  * [45] K. Lund, S. Massei, and D. Palitta, CRKSM_MatEQ, 2020, https://gitlab.com/katlund/compress-and-restart-KSM.
  * [46] W.-H. Luo, T.-Z. Huang, L. Li, Y. Zhang, and X.-M. Gu, Efficient preconditioner updates for unsymmetric shifted linear systems, Comput. Math. Appl., 67 (2014), pp. 1643 – 1655, https://doi.org/10.1016/j.camwa.2014.03.005.
  * [47] The MathWorks, Inc., MATLAB, https://www.matlab.com.
  * [48] K. Meerbergen, The solution of parametrized symmetric linear systems, SIAM J. Matrix Anal. Appl., 24 (2003), pp. 1038–1059, https://doi.org/10.1137/S0895479800380386.
  * [49] Oberwolfach Benchmark Collection, Steel profile.  hosted at MORwiki – Model Order Reduction Wiki, 2005, http://modelreduction.org/index.php/Steel_Profile.
  * [50] D. P. O’Leary, The block conjugate gradient algorithm and related methods, Linear Algebra Appl., 29 (1980), pp. 293 – 322, https://doi.org/10.1016/0024-3795(80)90247-5.
  * [51] D. Palitta and V. Simoncini, Matrix-equation-based strategies for convection-diffusion equations, BIT, 56 (2016), pp. 751–776, https://doi.org/10.1007/s10543-015-0575-8.
  * [52] M. L. Parks, E. de Sturler, G. Mackey, D. D. Johnson, and S. Maiti, Recycling Krylov subspaces for sequences of linear systems, SIAM J. Sci. Comput., 28 (2006), pp. 1651–1674, https://doi.org/10.1137/040607277.
  * [53] T. Penzl, Numerical solution of generalized Lyapunov equations, Adv. Comp. Math., 8 (1997), pp. 33–48, https://doi.org/10.1023/A:1018979826766.
  * [54] T. Penzl, A cyclic low rank Smith method for large sparse Lyapunov equations, SIAM J. Sci. Comput., 21 (2000), pp. 1401–1418, https://doi.org/10.1137/S1064827598347666.
  * [55] T. Penzl, Eigenvalue decay bounds for solutions of Lyapunov equations: the symmetric case, Systems Control Lett., 40 (2000), pp. 139–144, https://doi.org/10.1016/S0167-6911(00)00010-4.
  * [56] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, USA, second ed., 2003, https://doi.org/10.1137/1.9780898718003.
  * [57] Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Statist. Comput., 7 (1986), pp. 856–869, https://doi.org/10.1137/0907058.
  * [58] J. Saak, Efficient Numerical Solution of Large Scale Algebraic Matrix Equations in PDE Control and Model Order Reduction, PhD thesis, Technische Universität Chemnitz, Chemnitz, Germany, July 2009, http://nbn-resolving.de/urn:nbn:de:bsz:ch1-200901642.
  * [59] J. Saak, M. Köhler, and P. Benner, M-M.E.S.S.-2.2 – the Matrix Equations Sparse Solvers library, Feb. 2022, https://doi.org/10.5281/zenodo.5938237.  See also: https://www.mpi-magdeburg.mpg.de/projects/mess.
  * [60] J. Sabino, Solution of Large-Scale Lyapunov Equations via the Block Modified Smith Method, PhD thesis, Rice University, Houston, Texas, June 2007, http://www.caam.rice.edu/tech_reports/2006/TR06-08.pdf.
  * [61] T. Schmelzer, Block Krlyov Methods for Hermitian Linear Systems, PhD thesis, University of Kaiserslautern, 2004.
  * [62] V. Simoncini, Restarted full orthogonalization method for shifted linear systems, BIT, 43 (2003), pp. 459–466, https://doi.org/10.1023/A:1026000105893.
  * [63] V. Simoncini, A new iterative method for solving large-scale Lyapunov matrix equations, SIAM J. Sci. Comput., 29 (2007), pp. 1268–1288, https://doi.org/10.1137/06066120X.
  * [64] V. Simoncini, Extended Krylov subspace for parameter dependent systems, Appl. Numer. Math., 60 (2010), pp. 550–560, https://doi.org/10.1016/j.apnum.2010.03.001.
  * [65] V. Simoncini, Computational methods for linear matrix equations, SIAM Rev., 38 (2016), pp. 377–441, https://doi.org/10.1137/130912839.
  * [66] E. D. Sontag, Mathematical Control Theory, Texts in Applied Mathematics, Springer-Verlag, New York, 2nd ed., 1998, https://doi.org/10.1007/978-1-4612-0577-7.
  * [67] K. M. Soodhalter, Block Krylov subspace recycling for shifted systems with unrelated right-hand sides, SIAM J. Sci. Comput., 38 (2016), pp. A302–A324, https://doi.org/10.1137/140998214.
  * [68] K. M. Soodhalter, Two recursive GMRES-type methods for shifted linear systems with general preconditioning, Electron. Trans. Numer. Anal., 45 (2016), pp. 499–523, https://etna.math.kent.edu/volumes/2011-2020/vol45/abstract.php?vol=45&pages=499-523.
  * [69] K. M. Soodhalter, E. de Sturler, and M. Kilmer, A survey of subspace recycling iterative methods, GAMM Mitteilungen, 43 (2020), https://doi.org/10.1002/gamm.202000016.
  * [70] K. M. Soodhalter, D. B. Szyld, and F. Xue, Krylov subspace recycling for sequences of shifted linear systems, Appl. Numer. Math., 81 (2014), pp. 105 – 118, https://doi.org/10.1016/j.apnum.2014.02.006.
  * [71] G. Starke, Optimal alternating directions implicit parameters for nonsymmetric systems of linear equations, SIAM J. Numer. Anal., 28 (1991), pp. 1431–1445, https://doi.org/10.1137/0728074.
  * [72] The MORwiki Community, MORwiki - Model Order Reduction Wiki.  http://modelreduction.org.
  * [73] E. L. Wachspress, The ADI Model Problem, Springer New York, 2013, https://doi.org/10.1007/978-1-4614-5122-8.
  * [74] T. Wolf, $\mathcal{H}_{2}$ Pseudo-Optimal Model Order Reduction, PhD thesis, Technische Universität München, Munich, Germany, 2015, https://d-nb.info/1064075568/34.
  * [75] T. Wolf and H. Panzer, The ADI iteration for Lyapunov equations implicitly performs $\mathcal{H}_{2}$ pseudo-optimal model order reduction, Internat. J. Control, 89 (2016), pp. 481–493, https://doi.org/10.1080/00207179.2015.1081985.