Low-Rank Methods for Solving Discrete-Time Projected Lyapunov Equations

Abstract

In this paper, we consider the numerical solution of large-scale discrete-time projected Lyapunov equations. We provide some reasonable extensions of the most frequently used low-rank iterative methods for linear matrix equations, such as the low-rank Smith method and the low-rank alternating-direction implicit (ADI) method. We also consider how to reduce complex arithmetic operations and storage when shift parameters are complex and propose a partially real version of the low-rank ADI method. Through two standard numerical examples from discrete-time descriptor systems, we will show that the proposed low-rank alternating-direction implicit method is efficient.

1. Introduction

Solving linear matrix equations is a very important topic in control theory. Such equations include Lyapunov equations and Sylvester equations. Let $E,A\in {\mathbb{R}}^{n\times n}$, where E is singular. Assume that the pencil $\lambda E-A$ is d-stable. That is, the moduli of all finite eigenvalues of the pencil $\lambda E-A$ are less than 1. According to [1], there exist nonsingular $n\times n$ matrices $W,T$ that transform $E,A$ into a Weierstrass canonical form, i.e.,

$$E=W\left[\begin{array}{cc}I& 0\\ 0& N\end{array}\right]T,A=W\left[\begin{array}{cc}J& 0\\ 0& I\end{array}\right]T$$

(1)

with $J\in {\mathbb{R}}^{n_f\times n_f}$ and $N\in {\mathbb{R}}^{n_{\infty }\times n_{\infty }}$. Define the left and right spectral projection matrices $P_l$, $P_r$ by

$$P_l=W\left[\begin{array}{cc}I& 0\\ 0& 0\end{array}\right]W^{-1},P_r=T^{-1}\left[\begin{array}{cc}I& 0\\ 0& 0\end{array}\right]T.$$

(2)

In this paper, we focus on the numerical solution of the discrete-time projected Lyapunov equation

$$EX{E}^T-AX{A}^T=P_{l}B{B}^T{P}_l^T,X=P_{r}X{P}_r^T,$$

(3)

where $X\in {\mathbb{R}}^{n\times n}$ is the solution, $B\in {\mathbb{R}}^{n\times m},m\ll n$. Here and in the following, the superscript T denotes the transpose of a vector or a matrix. Since $\lambda E-A$ is d-stable, (3) has a symmetric positive semi-definite solution; see, for example, ref. [2]. The discrete-time Lyapunov equation is also called the Stein equation in the literature. By using the Kronecker product [3], the first equation in (3) can be formulated as $(E\otimes E-A\otimes A)\mathrm{vec}\left(X\right)=\mathrm{vec}\left(P_{l}B{B}^T{P}_l^T\right)$, where $\mathrm{vec}\left(X\right)={[x_1^T,x_2^T,\cdots ,x_n^T]}^T$, and $x_i$ is the i-th column of X.

The Stein equation with nonsingular E plays an essential role in discrete-time dynamical systems, including stability analysis and control [4,5,6], model reduction [7,8,9,10,11,12], solutions of discrete-time algebraic Riccati equations (by Newton’s method) in optimal control [13], and the restoration of images [14]. In contrast, the projected Stein equation arises in the balanced truncation model reduction [2] of discrete-time descriptor systems. In the positive real and bounded real balanced truncation model reduction of discrete-time descriptor systems, we also need to solve a projected Stein equation at each iteration step of Newton’s method for projected Riccati equations.

In the past few decades, many researchers have focused on constructing numerically robust algorithms for the standard Stein equation, i.e., with E being the identity matrix. For example, a standard direct method was provided in [15], which is a direct extension of the well-known Bartels–Stewart algorithm [16] for continuous-time Lyapunov equations $AX+X{A}^T=Q$ to the standard Stein equation. Hammarling [17] proposed a variant of the Bartels–Stewart algorithm for both the continuous-time and discrete-time cases. This variant is named the Hammarling method in the literature and aims to compute the Cholesky factor of the solution, which is desired in the balanced truncation model order reduction of discrete-time systems. The Hammarling algorithm was further improved in [18,19] by using a rank-2 updating formula. These approaches are based on the real Schur decomposition, require a computational complexity of $O\left(n^3\right)$ flops, and thus are only suitable for small to moderately sized problems.

It is known that the solution matrix of a continuous-time Lyapunov equation has a low numerical rank in cases where it has a low-rank right-hand side; see, e.g., ref. [20]. Specifically, the low numerical rank means that the singular values of the solution X decay very rapidly. Penzl [21] showed theoretically that the singular values of the solution decay exponentially for the continuous-time Lyapunov equation with a symmetric coefficient matrix and a low-rank right-hand side. Baker, Embree, and Sabino [22] considered the nonsymmetric case, and they explained that a larger departure from normality probably means a faster decay of singular values. The fact that the solution has a rapid decay of singular values and can be well approximated by its low-rank factorization now enables the use of numerous iterative methods that seek accurate low-rank approximations to the solution. These iterative methods include the low-rank Smith method [23,24], the Cholesky factor alternating-direction implicit (ADI) method [25], the (generalized) matrix sign function method [26], and the extended Krylov subspace method [27], to name a few. For the continuous-time projected Lyapunov equation, Stykel [28] extended the low-rank ADI method and the low-rank Smith method to compute low-rank approximations to the solution. In [13], the ADI method was extended to discrete-time Lyapunov equations and was further improved to compute the real factors in [29] by utilizing the technique that was proposed in [30] for continuous-time Lyapunov equations.

In recent years, numerical methods for continuous-time Lyapunov equations have been further considered. In [31], a class of low-rank iterative methods is proposed by using Runge–Kutta integration methods. It is shown that a special instance of this class of methods is equivalent to the low-rank ADI method. In [32], Benner, Palitta, and Saak further improved the low-rank ADI. They used the extended Krylov subspace method to solve the shifted linear system at each iteration. It is shown that by using only a single subspace, all the shifted linear systems can be solved to achieve a prescribed accuracy. In [33], the authors considered the inexact rational Krylov subspace method and low-rank iteration, in which a shifted linear system of equations is solved inexactly. In [34,35], the choice of shift parameters is considered, and some selection techniques are proposed to achieve a fast convergence for the low-rank ADI method.

In this paper, we first transform the projected generalized Stein Lyapunov Equation (3) to an equivalent projected standard Stein equation and then extend the low-rank Smith method to the projected standard equation. After this, we extend the low-rank ADI method to (3) and propose how to compute the real low-rank factor by following the idea in [29,30]. We also consider the choice of ADI shift parameters. Finally, through two standard numerical examples from discrete-time descriptor systems, we show the efficiency of the proposed low-rank ADI method.

The main contributions of this paper include the following:

• The low-rank ADI method is extended to solve the discrete-time projected Lyapunov equation.
• A partially real low-rank ADI algorithm is proposed.
• Two numerical examples are presented to demonstrate the good convergence of the low-rank ADI method.

The low-rank ADI method is one of the most commonly used iterative methods for solving linear matrix equations. It has a good convergence curve, although the shift parameters are not optimal. Moreover, it always produces a low-rank positive semi-definite approximate solution for Lyapunov equations, which is desired for some applications, such as the balanced truncation model order reduction. In contrast, the Krylov subspace method cannot guarantee the generation of the positive semi-definite solution. The main drawback of the low-rank ADI method is that it requires selecting shifts and solving one linear system of equations for each iteration.

The rest of the paper proceeds as follows. In Section 2, we reformulate (3) and propose the low-rank Smith method. In Section 3, we extend the real version of the low-rank ADI method for the projected Stein equation. Section 4 is devoted to two numerical examples. Finally, conclusions are given in Section 5.

2. Low-Rank Smith Method

The Smith method [36] was originally proposed for solving the continuous-time Sylvester equation $AX+X\tilde{A}=C$. First, the continuous-time equation is equivalently transformed to a discrete-time equation via a Cayley transform. Then, the Smith iteration is derived from the series representation of the solution. For the projected Stein Equation (3), the Smith method can be applied directly without the Cayley transform.

Due to the singularity of E, its inverse does not exist. We use the $\left\{2\right\}$-inverse $E^{-}$ of E, which is defined by

$$E^{-}=P_r{(E{P}_r+A(I-P_r))}^{-1}={(P_{l}E+(I-P_l)A)}^{-1}{P}_l=T^{-1}\left[\begin{array}{cc}I& 0\\ 0& 0\end{array}\right]W^{-1},$$

see, e.g., ref. [37]. For the generalized inverse of a singular matrix, the interested reader is referred to [38].

Multiplying the first equation in (3) from the left and right by $E^{-}$ and ${\left(E^{-}\right)}^T$ and using $E^{-}E=P_r$, $E^{-}{P}_l=E^{-}$, and $X=P_{r}X{P}_r^T$, we obtain

$$X-\left(E^{-}A\right)X{\left(E^{-}A\right)}^T=E^{-}B{\left(E^{-}B\right)}^T,X=P_{r}X{P}_r^T.$$

(4)

The unique solution of (4) can be formulated as

$$X=\sum _{j=0}^{\infty }{\left(E^{-}A\right)}^j{E}^{-}B{\left({\left(E^{-}A\right)}^j{E}^{-}B\right)}^T.$$

(5)

Since the pencil $\lambda E-A$ is d-stable, the spectrum radius $\rho \left(E^{-}A\right)$ of $E^{-}A$, which is defined by $\rho \left(E^{-}A\right)={max}_{\lambda \in \Lambda \left(E^{-}A\right)}\left|\lambda \right|$, satisfies $\rho \left(E^{-}A\right)<1$. So, the series converges, and the solution X is symmetric positive semi-definite, i.e., $X\ge 0$. This series representation of the solution implies that the numerical rank of X is much smaller than its dimension n if the norm of the powers of $E^{-}A$ decreases rapidly. In [39], Benner, Khoury, and Sadkane considered the solution of the Stein equation with $E=I_n$ and obtained an inequality that explicitly describes the decay of the singular values of the solution. For the projected Stein Equation (4), by following [39], we can obtain

$$\frac{{\sigma }_{jm+1}}{{\sigma }_1}\le {\parallel {\left(E^{-}A\right)}^j\parallel }^2,$$

where ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_n$ denote the singular values of the solution X. This result explicitly shows that the solution has a low numerical rank if the norm of the powers of $E^{-}A$ decreases rapidly.

We can now apply the Smith method [24,28] to (4). It is a fixed-point iteration and is expressed as

$$X_{j+1}=\left(E^{-}A\right)X_j{\left(E^{-}A\right)}^T+E^{-}B{\left(E^{-}B\right)}^T,X_0=0.$$

This iteration converges since $\rho \left(E^{-}A\right)<0$, and the iterations can be written as the partial sum

$$X_j=\sum _{i=0}^{j-1}{\left(E^{-}A\right)}^i{E}^{-}B{\left({\left(E^{-}A\right)}^i{E}^{-}B\right)}^T,j=1,2,\cdots .$$

(6)

We see from (6) that the iterations can be reformulated by a low-rank representation of Cholesky factors, i.e.,

$$X_j=Z_j{Z}_j^T$$

with $Z_j=[E^{-}B,\left(E^{-}A\right)E^{-}B,\cdots ,{\left(E^{-}A\right)}^{j-1}{E}^{-}B]$. It follows from (5) and (6) that the error matrix $X-X_j$ can be expressed as

$$\begin{array}{ccc}\hfill X-X_j& =& \sum _{i=j}^{\infty }{\left(E^{-}A\right)}^i{E}^{-}B{\left({\left(E^{-}A\right)}^i{E}^{-}B\right)}^T\hfill \\ & =& {\left(E^{-}A\right)}^j\left(\sum _{i=0}^{\infty }{\left(E^{-}A\right)}^i{E}^{-}B{\left({\left(E^{-}A\right)}^i{E}^{-}B\right)}^T\right){\left({\left(E^{-}A\right)}^j\right)}^T\hfill \\ & =& {\left(E^{-}A\right)}^{j}X{\left({\left(E^{-}A\right)}^j\right)}^T.\hfill \end{array}$$

Consequently, we can obtain the relative error bound

$$\frac{\parallel X-X_j\parallel }{\parallel X\parallel }\le {\parallel {\left(E^{-}A\right)}^j\parallel }^2.$$

This shows that $X_j$ converges linearly to the solution if the spectrum radius $\rho \left(E^{-}A\right)<1$, and X can be accurately approximately by the low-rank iteration $X_j=Z_j{Z}_j^T$ if the norm of the powers of $E^{-}A$ decreases rapidly. Note that the norm of the error matrix is not computable since the solution X is unknown. For large-scale problems, it is also difficult to accurately estimate the relative error bound $\parallel {\left(E^{-}A\right)}^j{\parallel }^2$.

For the residual matrix $R_j$ defined by

$$R_j=A{X}_j{A}^T+P_{l}B{B}^T{P}_l^T-E{X}_j{E}^T,$$

(7)

we have

$$\begin{array}{ccc}\hfill R_j& =& A\left(\sum _{i=0}^{j-1}{\left(E^{-}A\right)}^i{E}^{-}B{\left({\left(E^{-}A\right)}^i{E}^{-}B\right)}^T\right)A^T+P_{l}B{B}^T{P}_l^T\hfill \\ & & -E\left(\sum _{i=0}^{j-1}{\left(E^{-}A\right)}^i{E}^{-}B{\left({\left(E^{-}A\right)}^i{E}^{-}B\right)}^T\right)E^T\hfill \\ & =& E\left(\sum _{i=1}^j{\left(E^{-}A\right)}^i{E}^{-}B{\left({\left(E^{-}A\right)}^i{E}^{-}B\right)}^T\right)E^T+E\left(E^{-}{P}_{l}B{\left(E^{-}{P}_{l}B\right)}^T\right)E^T\hfill \\ & & -E\left(\sum _{i=0}^{j-1}{\left(E^{-}A\right)}^i{E}^{-}B{\left({\left(E^{-}A\right)}^i{E}^{-}B\right)}^T\right)E^T\hfill \\ & =& E\left({\left(E^{-}A\right)}^j{E}^{-}B{\left({\left(E^{-}A\right)}^j{E}^{-}B\right)}^T\right)E^T.\hfill \end{array}$$

So, the Frobenius matrix norm of $R_j$ can be easily computed.

The dimension of the low-rank factor $Z_j$ will increase by m in each iteration step. Hence, if the Smith iteration converges slowly, the number of columns of $Z_j$ will easily reach unmanageable levels of memory requirements. To reduce the dimension of $Z_j$, we will approximate it by using the rank-revealing QR decomposition (RRQR) [40]. Assume that $Z_j$ has the low numerical rank $r_j$ with a prescribed tolerance $\tau$. Consider the RRQR decomposition of $Z_j$ with column pivoting:

$$Z_j^T{\mathrm{\Pi }}_j=Q_j{\mathrm{\Omega }}_j,Q_j=[Q_j^{\left(1\right)},Q_j^{\left(2\right)}],{\mathrm{\Omega }}_j=\left[\begin{array}{cc}{\mathrm{\Omega }}_j^{\left(1\right)}& {\mathrm{\Omega }}_j^{\left(2\right)}\\ 0& {\mathrm{\Omega }}_j^{\left(3\right)}\end{array}\right],$$

where ${\mathrm{\Omega }}_j$ is an upper triangular matrix with ${\mathrm{\Omega }}_j^{\left(1\right)}\in {\mathbb{R}}^{r_j\times r_j}$ and $\parallel {\mathrm{\Omega }}_j^{\left(3\right)}\parallel <\tau$, $Q_j$ is orthogonal, ${\mathrm{\Pi }}_j$ is a permutation matrix, and $Q_j^{\left(1\right)}$ has $r_j$ columns. Then, $Z_j{Z}_j^T$ can be approximated by

$$\begin{array}{ccc}\hfill Z_j{Z}_j^T& \approx & {\mathrm{\Pi }}_j{\left[\begin{array}{cc}{\mathrm{\Omega }}_j^{\left(1\right)}& {\mathrm{\Omega }}_j^{\left(2\right)}\end{array}\right]}^T{\left(Q_j^{\left(1\right)}\right)}^T{Q}_j^{\left(1\right)}\left[\begin{array}{cc}{\mathrm{\Omega }}_j^{\left(1\right)}& {\mathrm{\Omega }}_j^{\left(2\right)}\end{array}\right]{\mathrm{\Pi }}_j={\tilde{Z}}_j{\tilde{Z}}_j^T,\hfill \end{array}$$

where

$${\tilde{Z}}_j={\mathrm{\Pi }}_j{\mathrm{\Omega }}_j^T\left[\begin{array}{c}{I}_{r_j}\\ 0\end{array}\right].$$

The low-rank Smith method for solving the projected Stein Equation (3) is presented in Algorithm 1.

Algorithm 1 Low-rank Smith method

Input: $E,A,B$, $\epsilon$, $\tau$. Output: Z such that $Z{Z}^T$ is the approximate solution of (3)  Set $j=1$, $V_1=E^{-}B$, $Z_1=V_1$;  Compute the rank-revealing QR decomposition $$[Q_1,{\mathrm{\Pi }}_1,{\mathrm{\Omega }}_1,r_1]=\mathrm{RRQR}(Z_1^T,\tau ).$$  Update $Z_1$ by $$Z_1={\mathrm{\Pi }}_1{\mathrm{\Omega }}_1^T\left[\begin{array}{c}{I}_{r_1}\\ 0\end{array}\right].$$  While $\parallel {\left(E{V}_j\right)}^T\left(E{V}_j\right){\parallel }_F>\epsilon$ do  •  $j=j+1$.  •  $V_j=E^{-}A{V}_{j-1}$.  •  $Z_j=[Z_{j-1},V_j]$.  •  Compute the rank-revealing QR decomposition $$[Q_j,{\mathrm{\Pi }}_j,{\mathrm{\Omega }}_j,r_j]=\mathrm{RRQR}(Z_j^T,\tau ).$$  •  Update $Z_j$ by $$Z_j={\mathrm{\Pi }}_j{\mathrm{\Omega }}_j^T\left[\begin{array}{c}{I}_{r_j}\\ 0\end{array}\right].$$  End While

The main advantage of the Smith iteration (6) is that it is very simple and can be easily implemented. However, we should note that the iterations converge very slowly if the spectrum radius $\rho \left(E^{-}A\right)\approx 1$. This is a significant motivation for further improvement of the Smith method.

3. Low-Rank ADI Method

The ADI method was first introduced in [41] and then applied to solve continuous-time Lyapunov matrix equations in [42]. Recently, this method was extended to the Stein Equation (3) by Benner and Faßbender [13] and further improved in [29]; see also [43].

For the projected Stein Equation (3), by generalizing the ADI method, we iteratively compute approximations $X_j$, $j\ge 1$, of the solution X by following the iteration scheme

$$\begin{array}{ccc}\hfill ({\overline{\mu }}_{j}A-E)X_{j-1/2}{A}^T& =& E{X}_{j-1}({\overline{\mu }}_j{E}^T-A^T)-{\overline{\mu }}_j{P}_{l}B{B}^T{P}_l^T,\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill E{X}_j(E^T-{\mu }_j{A}^T)& =& (A-{\mu }_{j}E)X_{j-1/2}{A}^T+P_{l}B{B}^T{P}_l^T,\hfill \end{array}$$

(9)

where $0<|{\mu }_j|<1$ denotes suitable shift parameters. Note that, although the iteration can work with any initial guess $X_0$, we use only $X_0=0$ in the sequel.

Since the pencil $\lambda E-A$ is d-stable and $0<|{\mu }_j|<1$, the matrices ${\overline{\mu }}_{j}A-E$ and $E^T-{\mu }_j{A}^T$ are nonsingular. From (8), we obtain

$$X_{j-1/2}{A}^T={({\overline{\mu }}_{j}A-E)}^{-1}E{X}_{j-1}({\overline{\mu }}_j{E}^T-A^T)-{\overline{\mu }}_j{({\overline{\mu }}_{j}A-E)}^{-1}{P}_{l}B{B}^T{P}_l^T,$$

(10)

Then, these half steps in the ADI iteration are rewritten into single steps by substituting $X_{j-1/2}{A}^T$ into (9) by the expression (10) for $X_{j-1/2}{A}^T$; i.e., we arrive at the single-step iteration

$$\begin{array}{ccc}\hfill E{X}_j& =& (A-{\mu }_{j}E)X_{j-1/2}{A}^T{(E-{\mu }_{j}A)}^{-T}+P_{l}B{B}^T{P}_l^T{(E-{\mu }_{j}A)}^{-T}\hfill \\ & =& (A-{\mu }_{j}E){({\overline{\mu }}_{j}A-E)}^{-1}E{X}_{j-1}({\overline{\mu }}_j{E}^T-A^T){(E-{\mu }_{j}A)}^{-T}\hfill \\ & & -{\overline{\mu }}_j(A-{\mu }_{j}E){({\overline{\mu }}_{j}A-E)}^{-1}{P}_{l}B{B}^T{P}_l^T{(E-{\mu }_{j}A)}^{-T}\hfill \\ & & +P_{l}B{B}^T{P}_l^T{(E-{\mu }_{j}A)}^{-T}.\hfill \end{array}$$

(11)

Observe that

$$\begin{array}{ccc}& & -{\overline{\mu }}_j(A-{\mu }_{j}E){({\overline{\mu }}_{j}A-E)}^{-1}{P}_{l}B{B}^T{P}_l^T{(E-{\mu }_{j}A)}^{-T}+P_{l}B{B}^T{P}_l^T{(E-{\mu }_{j}A)}^{-T}\hfill \\ & =& \left(-{\overline{\mu }}_j(A-{\mu }_{j}E){({\overline{\mu }}_{j}A-E)}^{-1}+I\right)P_{l}B{B}^T{P}_l^T{(E-{\mu }_{j}A)}^{-T}\hfill \\ & =& \left(-{\overline{\mu }}_j(A-{\mu }_{j}E)+({\overline{\mu }}_{j}A-E)\right){({\overline{\mu }}_{j}A-E)}^{-1}{P}_{l}B{B}^T{P}_l^T{(E-{\mu }_{j}A)}^{-T}\hfill \\ & =& \left(1-|{\mu }_j{|}^2\right)E{({\overline{\mu }}_{j}A-E)}^{-1}{P}_{l}B{B}^T{P}_l^T{({\mu }_{j}A-E)}^{-T}.\hfill \end{array}$$

(12)

Hence, we obtain

$$\begin{array}{ccc}\hfill E{X}_j& =& (A-{\mu }_{j}E){({\overline{\mu }}_{j}A-E)}^{-1}E{X}_{j-1}(A^T-{\overline{\mu }}_j{E}^T){({\mu }_{j}A-E)}^{-T}\hfill \\ & & +\left(1-|{\mu }_j{|}^2\right)E{({\overline{\mu }}_{j}A-E)}^{-1}{P}_{l}B{B}^T{P}_l^T{({\mu }_{j}A-E)}^{-T}.\hfill \end{array}$$

(13)

Multiplying (13) from the left by $E^{-}$ and using $X_j=P_r{X}_j{P}_r^T$,

$$\begin{array}{ccc}\hfill P_r{(\mu A-E)}^{-1}& =& {(\mu A-E)}^{-1}{P}_l,\hfill \\ \hfill P_l(A-\mu E)& =& (A-\mu E)P_r\hfill \end{array}$$

for any $\mu$, we obtain

$$\begin{array}{ccc}\hfill X_j& =& {({\overline{\mu }}_{j}A-E)}^{-1}(A-{\mu }_{j}E)X_{j-1}(A^T-{\overline{\mu }}_j{E}^T){({\mu }_{j}A-E)}^{-T}\hfill \\ & & +\left(1-|{\mu }_j{|}^2\right){({\overline{\mu }}_{j}A-E)}^{-1}{P}_{l}B{B}^T{P}_l^T{({\mu }_{j}A-E)}^{-T}.\hfill \end{array}$$

(14)

One can easily verify that the solution X of the projected Stein Equation (3) is a fixed point of the single-step iteration (14). That is to say,

$$\begin{array}{ccc}\hfill X& =& {({\overline{\mu }}_{j}A-E)}^{-1}(A-{\mu }_{j}E)X(A^T-{\overline{\mu }}_j{E}^T){({\mu }_{j}A-E)}^{-T}\hfill \\ & & +\left(1-|{\mu }_j{|}^2\right){({\overline{\mu }}_{j}A-E)}^{-1}{P}_{l}B{B}^T{P}_l^T{({\mu }_{j}A-E)}^{-T}.\hfill \end{array}$$

(15)

Consequently, from (14) and (15), we obtain the following recursive formulation for the error matrix between the solution X and the approximation $X_j$:

$$X-X_j={({\overline{\mu }}_{j}A-E)}^{-1}(A-{\mu }_{j}E)(X-X_{j-1})(A^T-{\overline{\mu }}_j{E}^T){({\mu }_{j}A-E)}^{-T},j\ge 1.$$

(16)

We see from (16) and $X_0=0$ that the error matrix $X-X_j$ can be written as

$$X-X_j=\left(\prod _{i=1}^j{({\overline{\mu }}_{i}A-E)}^{-1}(A-{\mu }_{i}E)\right)X{\left(\prod _{i=1}^j{({\overline{\mu }}_{i}A-E)}^{-1}(A-{\mu }_{i}E)\right)}^H.$$

(17)

3.1. Low-Rank Version of ADI Method

For continuous-time Lyapunov equations with a low-rank right-hand side, Li and White [25] proposed the state-of-the-art Cholesky factor ADI algorithm, which generates a low-rank approximation to the solution. This method is a significant improvement of the ADI method [42] and is very appropriate for large-scale continuous-time Lyapunov equations. The Cholesky factor ADI method is developed by exploiting the low-rank structure of the iterations and reordering the shifts; see [25] for the details. The low-rank ADI method is generalized to the Stein equation in [29].

In this section, we follow these ideas to deduce the low-rank ADI method for the projected Stein Equation (3). Suppose now that $X_j$ and $X_{j-1}$ are written in their factored forms: $X_j=Z_j{Z}_j^H$ and $X_{j-1}=Z_{j-1}{Z}_{j-1}^H$. Then, from (14), we obtain the following recursive relation for $Z_j$ and $Z_{j-1}$:

$$Z_j=\left[\begin{array}{cc}\sqrt{1-|{\mu }_j{|}^2}{({\overline{\mu }}_{j}A-E)}^{-1}{P}_{l}B& {({\overline{\mu }}_{j}A-E)}^{-1}(A-{\mu }_{j}E)Z_{j-1}\end{array}\right]$$

(18)

with $Z_0=0$.

From (18), we easily see that the dimension of $Z_j$ would increase by m at each iteration step. Therefore, the factor $Z_j$ of $X_j$ has $mj$ columns. Since the number of columns increases by m at each iteration, the number of systems of linear equations with matrices $\overline{\mu }A-E$, which need to be solved at each iteration in the low-rank ADI method (18), increases by m. So, this iteration (18) for the factors is not suitable for practical implementation. By making use of the trick in [25], $Z_j$ can be reformulated as

$$\begin{array}{ccc}\hfill Z_j& =& \left[\sqrt{1-|{\mu }_1{|}^2}{T}_1{P}_{l}B,\sqrt{1-|{\mu }_2{|}^2}{T}_2{S}_1{T}_1{P}_{l}B,\right.\hfill \\ & & \left.\cdots ,\sqrt{1-|{\mu }_j{|}^2}{T}_j{S}_{j-1}{T}_{j-1}\cdots S_2{T}_2{S}_1{T}_1{P}_{l}B\right],\hfill \end{array}$$

where

$$T_j={({\overline{\mu }}_{j}A-E)}^{-1},S_j=(A-{\mu }_{j}E).$$

This directly leads to an efficient low-rank ADI iteration scheme: Let $V_1={({\overline{\mu }}_{1}A-E)}^{-1}{P}_{l}B$, and $Z_1=\sqrt{1-|{\mu }_1{|}^2}{V}_1$. Then, for $j\ge 1$,

$$\begin{array}{ccc}\hfill {\tilde{V}}_j& =& (A-{\mu }_{j}E)V_j,\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill V_{j+1}& =& {({\overline{\mu }}_{j+1}A-E)}^{-1}{\tilde{V}}_j,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill Z_{j+1}& =& [Z_j,\sqrt{1-|{\mu }_{j+1}{|}^2}{V}_{j+1}].\hfill \end{array}$$

(21)

We now investigate the residual matrix $R_j$ corresponding to the j-th approximate solution $X_j$. In the following theorem, we show that $R_j$ has a low-rank factorization of rank at most m.

Theorem 1.

Let $X_j=Z_j{Z}_j^H$, where $Z_j$ is the j-th iteration generated by the low-rank ADI for the projected Stein Equation (3), and let the $n\times m$ matrix ${\tilde{V}}_j$ be defined by (19). Then, the residual matrix $R_j$, defined by (7), can be formulated as

$$R_j={\tilde{V}}_j{\tilde{V}}_j^H.$$

Proof. 

From (3), $P_{l}B{B}^T{P}_l^T=EX{E}^T-AX{A}^T$. Thus,

$$R_j=E(X-X_j)E^T-A(X-X_j)A^T.$$

(22)

By inserting (17) into (22), we obtain

$$\begin{array}{ccc}\hfill R_j& =& E\left(\prod _{i=1}^j{({\overline{\mu }}_{i}A-E)}^{-1}(A-{\mu }_{i}E)\right)X{\left(\prod _{i=1}^j{({\overline{\mu }}_{i}A-E)}^{-1}(A-{\mu }_{i}E)\right)}^H{E}^T\hfill \\ & & -A\left(\prod _{i=1}^j{({\overline{\mu }}_{i}A-E)}^{-1}(A-{\mu }_{i}E)\right)X{\left(\prod _{i=1}^j{({\overline{\mu }}_{i}A-E)}^{-1}(A-{\mu }_{i}E)\right)}^H{A}^T\hfill \\ & =& \left(\prod _{i=1}^j{({\overline{\mu }}_{i}A-E)}^{-1}(A-{\mu }_{i}E)\right)P_{l}B{B}^T{P}_l^T{\left(\prod _{i=1}^j{({\overline{\mu }}_{i}A-E)}^{-1}(A-{\mu }_{i}E)\right)}^H.\hfill \end{array}$$

From (19) and (20), and $V_1={({\overline{\mu }}_{1}A-E)}^{-1}{P}_{l}B$, it follows that

$${\tilde{V}}_j=\left(\prod _{i=1}^j{({\overline{\mu }}_{i}A-E)}^{-1}(A-{\mu }_{i}E)\right)P_{l}B.$$

Thus, $R_j={\tilde{V}}_j{\tilde{V}}_j^H$.    □

The following theorem states that $V_{j+1}$ and ${\tilde{V}}_{j+1}$ can be obtained without solving systems of linear equations once $V_j$ has been computed.

Theorem 2.

Assume that a proper set of shift parameters is used in the low-rank ADI iteration. For the two subsequent blocks $V_{j+1}$ and ${\tilde{V}}_{j+1}$ related to the pair of complex shifts $\{{\mu }_j,{\mu }_{j+1}\}$ with ${\mu }_{j+1}={\overline{\mu }}_j$, it holds that

$$\begin{array}{ccc}\hfill V_{j+1}& =& {\mu }_j\mathrm{Re}\left(V_j\right)+\frac{1}{{\overline{\mu }}_j}\left((1-|{\mu }_j{|}^2)\frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}-j\right)\mathrm{Im}\left(V_j\right),\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill {\tilde{V}}_{j+1}& =& \frac{1}{|{\mu }_j{|}^2}\left({\tilde{V}}_{j-1}+(1-|{\mu }_j{|}^4)E\mathrm{Re}\left(V_j\right)+(1-|{\mu }_j{{|}^2)}^2\frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}E\mathrm{Im}\left(V_j\right)\right).\hfill \end{array}$$

(24)

Moreover, after applying a pair of complex shifts $\{{\mu }_j,{\mu }_{j+1}\}$, the generated ${\tilde{V}}_{j+1}$ is real.

Proof. 

Although the proof is similar to that of Theorem 1 in [30], we include the proof for completeness. In [30], the proof is split into three cases concerning different possible (sub)sequences of shift parameters.

Here, we only consider the first case in [30] for illustration. Assume that ${\mu }_1,\cdots ,{\mu }_{j-1}$ are real, ${\mu }_j$ has a nonzero imaginary part, and ${\mu }_{j+1}={\overline{\mu }}_j$. In this case, obviously, from (19) and (20), $V_1,\cdots ,V_{j-1},$${\tilde{V}}_1,\cdots ,{\tilde{V}}_{j-1}$ are real, and $V_j$ is the first complex iteration. From (20), it follows that

$$({\overline{\mu }}_{j}A-E)V_j={\tilde{V}}_{j-1}.$$

Now, splitting ${\mu }_j$ and $V_j$ into their real and imaginary parts reveals

$$\begin{array}{ccc}\hfill {\tilde{V}}_{j-1}& =& (\mathrm{Re}\left({\mu }_j\right)A-E)\mathrm{Re}\left(V_j\right)+\mathrm{Im}\left({\mu }_j\right)A\mathrm{Im}\left(V_j\right)\hfill \\ & & +j[(\mathrm{Re}\left({\mu }_j\right)A-E)\mathrm{Im}\left(V_j\right)-\mathrm{Im}\left({\mu }_j\right)A\mathrm{Re}\left(V_j\right)].\hfill \end{array}$$

Since ${\tilde{V}}_{j-1}$ is real, then

$$\begin{array}{ccc}\hfill (\mathrm{Re}\left({\mu }_j\right)A-E)\mathrm{Im}\left(V_j\right)-\mathrm{Im}\left({\mu }_j\right)A\mathrm{Re}\left(V_j\right)& =& 0,\hfill \end{array}$$

which leads to

$$\begin{array}{ccc}\hfill A\mathrm{Re}\left(V_j\right)& =& \frac{1}{\mathrm{Im}\left({\mu }_j\right)}(\mathrm{Re}\left({\mu }_j\right)A-E)\mathrm{Im}\left(V_j\right),\hfill \\ \hfill E{V}_j& =& E\mathrm{Re}\left(V_j\right)+jE\mathrm{Im}\left(V_j\right)\hfill \end{array}$$

(25)

$$\begin{array}{ccc}& =& (E-j\mathrm{Im}\left({\mu }_j\right)A)\mathrm{Re}\left(V_j\right)+j\mathrm{Re}\left({\mu }_j\right)A\mathrm{Im}\left(V_j\right).\hfill \end{array}$$

(26)

From (19) and (20), it follows that

$$\begin{array}{ccc}\hfill {\tilde{V}}_j& =& (A-{\mu }_{j}E)V_j=(A-{\mu }_{j}E){({\overline{\mu }}_{j}A-E)}^{-1}{\tilde{V}}_{j-1}\hfill \\ & =& \frac{1}{{\overline{\mu }}_j}({\overline{\mu }}_{j}A-E+E-|{\mu }_j{|}^{2}E){({\overline{\mu }}_{j}A-E)}^{-1}{\tilde{V}}_{j-1}\hfill \\ & =& \frac{1}{{\overline{\mu }}_j}{\tilde{V}}_{j-1}+\frac{1-|{\mu }_j{|}^2}{{\overline{\mu }}_j}E{({\overline{\mu }}_{j}A-E)}^{-1}{\tilde{V}}_{j-1}\hfill \\ & =& \frac{1}{{\overline{\mu }}_j}{\tilde{V}}_{j-1}+\frac{1-|{\mu }_j{|}^2}{{\overline{\mu }}_j}E{V}_j.\hfill \end{array}$$

(27)

From (20), we obtain

$$\begin{array}{ccc}\hfill {\overline{\mu }}_j{V}_{j+1}& =& {\overline{\mu }}_j{({\overline{\mu }}_{j+1}A-E)}^{-1}{\tilde{V}}_j\hfill \\ & =& {({\mu }_{j}A-E)}^{-1}\left({\tilde{V}}_{j-1}+(1-|{\mu }_j{|}^2)E{V}_j\right)\hfill \\ & =& {\overline{V}}_j+(1-|{\mu }_j{|}^2){({\mu }_{j}A-E)}^{-1}E{V}_j\hfill \\ & =& |{\mu }_j{|}^2\mathrm{Re}\left(V_j\right)-j\mathrm{Im}\left(V_j\right)\hfill \\ & & +(1-|{\mu }_j{|}^2){({\mu }_{j}A-E)}^{-1}(({\mu }_{j}A-E)\mathrm{Re}\left(V_j\right)+E{V}_j).\hfill \end{array}$$

(28)

From (25) and (26), we obtain

$$\begin{array}{ccc}& & ({\mu }_{j}A-E)\mathrm{Re}\left(V_j\right)+E{V}_j\hfill \\ & =& ({\mu }_{j}A-E)\mathrm{Re}\left(V_j\right)+(E-j\mathrm{Im}\left({\mu }_j\right)A)\mathrm{Re}\left(V_j\right)+j\mathrm{Re}\left({\mu }_j\right)A\mathrm{Im}\left(V_j\right)\hfill \\ & =& \mathrm{Re}\left({\mu }_j\right)A\mathrm{Re}\left(V_j\right)+j\mathrm{Re}\left({\mu }_j\right)A\mathrm{Im}\left(V_j\right)\hfill \\ & =& \frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}(\mathrm{Re}\left({\mu }_j\right)A-E)\mathrm{Im}\left(V_j\right)+j\mathrm{Re}\left({\mu }_j\right)A\mathrm{Im}\left(V_j\right)\hfill \\ & =& \frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}({\mu }_{j}A-E)\mathrm{Im}\left(V_j\right).\hfill \end{array}$$

By inserting (29) into (28), we obtain

$$\begin{array}{ccc}\hfill {\overline{\mu }}_j{V}_{j+1}& =& |{\mu }_j{|}^2\mathrm{Re}\left(V_j\right)+\left((1-|{\mu }_j{|}^2)\frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}-j\right)\mathrm{Im}\left(V_j\right),\hfill \end{array}$$

which leads to (23).

From (27) and ${\mu }_{j+1}={\overline{\mu }}_j$, it follows that

$$\begin{array}{ccc}\hfill {\tilde{V}}_{j+1}& =& \frac{1}{{\mu }_j}{\tilde{V}}_j+\frac{1-|{\mu }_j{|}^2}{{\mu }_j}E{V}_{j+1}\hfill \\ & =& \frac{1}{{\mu }_j}\left(\frac{1}{{\overline{\mu }}_j}{\tilde{V}}_{j-1}+\frac{1-|{\mu }_j{|}^2}{{\overline{\mu }}_j}E{V}_j\right)\hfill \\ & & +\frac{1-|{\mu }_j{|}^2}{{\mu }_j}E\left({\mu }_j\mathrm{Re}\left(V_j\right)+\frac{1}{{\overline{\mu }}_j}\left((1-|{\mu }_j{|}^2)\frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}-j\right)\mathrm{Im}\left(V_j\right)\right).\hfill \\ & =& \frac{1}{|{\mu }_j{|}^2}{\tilde{V}}_{j-1}+\frac{1-|{\mu }_j{|}^2}{|{\mu }_j{|}^2}E{V}_j+(1-|{\mu }_j{|}^2)E\mathrm{Re}\left(V_j\right)\hfill \\ & & -j\frac{1-|{\mu }_j{|}^2}{|{\mu }_j{|}^2}E\mathrm{Im}\left(V_j\right)+\frac{(1-|{\mu }_j{{|}^2)}^2}{|{\mu }_j{|}^2}\frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}E\mathrm{Im}\left(V_j\right)\hfill \\ & =& \frac{1}{|{\mu }_j{|}^2}\left({\tilde{V}}_{j-1}+(1-|{\mu }_j{|}^4)E\mathrm{Re}\left(V_j\right)+(1-|{\mu }_j{{|}^2)}^2\frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}E\mathrm{Im}\left(V_j\right)\right).\hfill \end{array}$$

Obviously, ${\tilde{V}}_{j+1}$ is real.    □

3.2. Dealing with Complex Shifts

From $V_j={({\overline{\mu }}_{j}A-E)}^{-1}{\tilde{V}}_{j-1}$, it follows that if ${\mu }_j$ is a complex shift, then the complex $V_j$ will be added to the low-rank factor $Z_j$. In this case, complex arithmetic operations and storage are introduced into the process such that a complex low-rank factor $Z_j$ is generated in the end. From the numerical point of view, it is undesirable to use complex arithmetic operations in the iteration since operations and storage will increase.

Recently, some research has focused on how to deal with complex shift parameters in the Cholesky factor ADI method for continuous-time Lyapunov equations. In [24,25,44], a completely real formulation of the Cholesky factor ADI method is presented by concatenating steps associated with a pair of complex conjugate shift parameters into one step. Although this reformulation has the advantage that complex arithmetic operations and storage are avoided, systems of linear equations with matrices of the form $A^2+2\mathrm{Re}\left({\mu }_j\right)A+{\left|{\mu }_j\right|}^2{I}_n$ need to be solved in every two steps of the ADI method. This is a major drawback for the completely real formulation. Firstly, for large-scale problems, $A^2+2\mathrm{Re}\left({\mu }_j\right)A+{\left|{\mu }_j\right|}^2{I}_n$ may not preserve the original sparsity of A, and thus, sparse direct solvers cannot be applied to linear systems with such coefficient matrices. Secondly, from the perspective of numerical stability, it is undesirable to solve such linear equations since the condition number can be increased due to squaring. Iterative solvers such as Krylov subspace methods [45,46] can still be applied to linear systems with such coefficient matrices since they work with matrix–vector products only. However, it is known that the large condition number will deteriorate the efficiency of iterative solvers. In order to overcome these disadvantages in the completely real formulation and to avoid complex arithmetic and the storage of complex matrices as much as possible, Benner, Kürschner, and Saak [30] introduced a partially real reformulation of the Cholesky factor ADI method for continuous-time Lyapunov equations. They exploit the fact that the ADI shifts need to occur as a real number or as a pair of complex conjugate numbers. As a result, the resulting low-rank ADI method works with real low-rank factors $Z_j$; see also [47]. This idea is extended to the Stein equation in [29].

We now consider the generalization to obtain a partially real version of the low-rank ADI method for the projected Stein Equation (3) by investigating the blocks $V_j,V_{j+1}$, which are generated in the low-rank ADI with a pair of complex conjugate shifts ${\mu }_j,{\mu }_{j+1}={\overline{\mu }}_j$.

Define

$$\begin{array}{ccc}\hfill \widehat{Z}& =& \left[\begin{array}{cc}\sqrt{1-|{\mu }_j{|}^2}{V}_j& \sqrt{1-|{\mu }_{j+1}{|}^2}{V}_{j+1}\end{array}\right],\hfill \\ \hfill \tilde{Z}& =& \left[\begin{array}{cc}\mathrm{Re}\left(V_j\right)& \mathrm{Im}\left(V_j\right)\end{array}\right].\hfill \end{array}$$

From (19) and (23), with the help of the Kronecker product and ${\mu }_{j+1}={\overline{\mu }}_j$, we obtain

$$\widehat{Z}=\tilde{Z}(\mathrm{\Xi }\otimes I_m),$$

(29)

where

$$\begin{array}{ccc}\hfill \mathrm{\Xi }& =& \sqrt{1-|{\mu }_j{|}^2}\left[\begin{array}{cc}1& {\mu }_j\\ j& \frac{1}{{\overline{\mu }}_j}\left((1-|{\mu }_j{|}^2)\frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}-j\right)\end{array}\right].\hfill \end{array}$$

Direct calculation reveals

$$\begin{array}{ccc}\hfill \mathrm{\Xi }{\mathrm{\Xi }}^H& =& (1-|{\mu }_j{|}^2)\left[\begin{array}{cc}1+|{\mu }_j{|}^2& (1-|{\mu }_j{|}^2)\frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}\\ (1-|{\mu }_j{|}^2)\frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}& 1+\frac{1}{|{\mu }_j{|}^2}\left((1-|{\mu }_j{{|}^2)}^2\frac{{\left(\mathrm{Re}\left({\mu }_j\right)\right)}^2}{{\left(\mathrm{Im}\left({\mu }_j\right)\right)}^2}+1\right)\end{array}\right].\hfill \end{array}$$

The real symmetric positive definite matrix $\mathrm{\Xi }{\mathrm{\Xi }}^H$ has a unique Cholesky factorization given by

$$\mathrm{\Xi }{\mathrm{\Xi }}^H=L{L}^T,L=\left[\begin{array}{cc}{l}_1& 0\\ l_2& l_3\end{array}\right],$$

(30)

where

$$\begin{array}{ccc}\hfill l_1& =& \sqrt{1-|{\mu }_j{|}^4},\hfill \\ \hfill l_2& =& l_1^{-1}(1-|{\mu }_j{{|}^2)}^2\frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)},\hfill \\ \hfill l_3& =& \sqrt{(1-|{\mu }_j{|}^2)\left(1+\frac{1}{|{\mu }_j{|}^2}\left((1-|{\mu }_j{{|}^2)}^2\frac{{\left(\mathrm{Re}\left({\mu }_j\right)\right)}^2}{{\left(\mathrm{Im}\left({\mu }_j\right)\right)}^2}+1\right)\right)-l_2^2}.\hfill \end{array}$$

(31)

From (29) and the Cholesky factorization (30) of $\mathrm{\Pi }{\mathrm{\Pi }}^H$, we obtain a real low-rank expression for $\widehat{Z}{\widehat{Z}}^H$:

$$\widehat{Z}{\widehat{Z}}^H=\breve{Z}{\breve{Z}}^T$$

with

$$\breve{Z}=\left[\begin{array}{cc}{l}_1\mathrm{Re}\left(V_j\right)+l_2\mathrm{Im}\left(V_j\right)& l_3\mathrm{Im}\left(V_j\right)\end{array}\right].$$

In summary, we can propose a partially real low-rank ADI method for solving the projected Stein Equation (3), which is presented in Algorithm 2.

Algorithm 2 Real low-rank ADI method

Input: $E,A,B$, $\epsilon$, ${\mu }_j$ with $0<|{\mu }_j|<1$. Output: Z such that $Z{Z}^T$ is the approximate solution of (3)  Set $j=1$, ${\tilde{V}}_0=P_{l}B$, $Z=\left[\right]$;  While $\parallel {\tilde{V}}_{j-1}^T{\tilde{V}}_{j-1}{\parallel }_F>\epsilon$ do  •  Solve $V_j={({\overline{\mu }}_{j}A-E)}^{-1}{\tilde{V}}_{j-1}$ for $V_j$.  •  If   $\mathrm{Im}\left({\mu }_j\right)=0$ then  – $Z=[Z,\sqrt{1-|{\mu }_j{|}^2}{V}_j]$.  – ${\tilde{V}}_j=\frac{1}{{\mu }_j}\left({\tilde{V}}_{j-1}+(1-|{\mu }_j{|}^2)E{V}_j\right)$.  •  else  – Compute $l_1,l_2,l_3$ according to (31).  – Set $Z=[Z,l_1\mathrm{Re}\left(V_j\right)+l_2\mathrm{Im}\left(V_j\right),l_3\mathrm{Im}\left(V_j\right)]$.  – ${\tilde{V}}_{j+1}=\frac{1}{|{\mu }_j{|}^2}\left({\tilde{V}}_{j-1}+(1-|{\mu }_j{|}^4)E\mathrm{Re}\left(V_j\right)+(1-|{\mu }_j{{|}^2)}^2\frac{\mathrm{Re}\left({\mu }_j\right)}{\mathrm{Im}\left({\mu }_j\right)}E\mathrm{Im}\left(V_j\right)\right)$.  – $j=j+1$.  •  End If  •  $j=j+1$.  End While

3.3. Choosing the ADI Shift Parameters

We now consider how to compute appropriate shift parameters. These shifts are vitally important to the convergence rate of the ADI iteration.

For the ADI method for the projected Stein equation, the parameters ${\left\{{\mu }_j\right\}}_{j=1}^k$ should be chosen according to the minimax problem

$$\underset{\begin{array}{c}0<|{\mu }_j|<1,\\ j=1,2,\dots ,k\end{array}}{min}\underset{\begin{array}{c}t\in {\Lambda }_f\end{array}}{max}\prod _{j=1}^k\left|\frac{t-{\mu }_j}{{\overline{\mu }}_{j}t-1}\right|,$$

(32)

where ${\Lambda }_f$ denotes the set of finite eigenvalues of the pencil $\lambda E-A$. In practice, since the eigenvalues of the pencil $\lambda E-A$ are unknown and computationally expensive, ${\Lambda }_f$ is usually replaced by a domain containing a finite set of eigenvalues of $\lambda E-A$.

A heuristic algorithm [21,48] can calculate the suboptimal ADI shift parameters for standard Lyapunov or Sylvester equations. It selects suboptimal ADI parameters from a set $\mathrm{\Omega }$, which is taken as the union of Ritz values of A and the reciprocals of the Ritz values of $A^{-1}$, obtained by two Arnoldi processes, with A and $A^{-1}$.

The heuristic algorithm can also be naturally extended to the minimax problem (32). Since E is assumed to be singular, the inverse of E does not exist. However, it is clear that $E^{-}A$ has the same nonzero finite eigenvalues as the pencil $\lambda E-A$. Moreover, the reciprocals of the largest nonzero eigenvalues of $A^{-1}E$ are the smallest eigenvalues of $E^{-}A$. Thus, we can run one Arnoldi process with the matrix $A^{-}E$ to compute the smallest nonzero eigenvalues of $E^{-}A$.

The algorithm for computing ${\left\{{\mu }_j\right\}}_{j=1}^k$ is summarized in Algorithm 3. For more details about the implementation of this algorithm, the interested reader is referred to [24,28,48].

Algorithm 3 Choose ADI parameters

Input: $E,A\in {\mathbb{R}}^{n\times n}$ with $\lambda E-A$ being d-stable, $b\in {\mathbb{R}}^n$, $k_{+}$, $k_{-}$. Output: ADI parameters $\mathbb{P}$.  Run $k_{+}$ steps of the Arnoldi process with respect to $E^{-}A$ on b to obtain the set ${\mathrm{\Omega }}_{+}$ of Ritz values.  Run $k_{-}$ steps of the Arnoldi process with respect to $A^{-1}E$ on b to obtain the set ${\mathrm{\Omega }}_{-}$ of Ritz values.  Set $\mathrm{\Omega }={\mathrm{\Omega }}_{+}\cup (1/{\mathrm{\Omega }}_{-})$.  Set $${\mu }_1=arg\underset{\begin{array}{c}\mu \in \mathrm{\Omega }\end{array}}{min}\underset{\begin{array}{c}t\in \mathrm{\Omega }\end{array}}{max}\left|\frac{t-\mu }{\overline{\mu }t-1}\right|.$$  If $Im\left({\mu }_1\right)=0$, $\mathbb{P}=\mathbb{P}\cup \left\{{\mu }_1\right\}$, $j=1$; else $\mathbb{P}=\mathbb{P}\cup \{{\mu }_1,{\mu }_2={\overline{\mu }}_1\}$, $j=2$.  While $j<k$ do  •  Set $${\mu }_{j+1}=arg\underset{\begin{array}{c}\mu \in {\mathrm{\Omega }}^{\prime }\end{array}}{min}\underset{\begin{array}{c}t\in \mathrm{\Omega }\end{array}}{max}\left|\frac{t-\mu }{\overline{\mu }t-1}\right|\prod _{i=1}^j\left|\frac{t-{\mu }_i}{{\overline{\mu }}_{i}t-1}\right|,$$ where ${\mathrm{\Omega }}^{\prime }$ is $\mathrm{\Omega }$ with $\mathbb{P}$ deleted.  •  If $Im\left({\mu }_{j+1}\right)=0$, $\mathbb{P}=\mathbb{P}\cup \left\{{\mu }_{j+1}\right\}$, $j=j+1$; else $\mathbb{P}=\mathbb{P}\cup \{{\mu }_{j+1},{\mu }_{j+2}={\overline{\mu }}_{j+1}\}$, $j=j+2$.  End While

4. Numerical Examples

We provide two numerical examples to demonstrate the convergence performance of the LR-ADI method and the LR-Smith method for (3) in this section. Define the relative residual (RRes) as

$$\mathrm{RRes}\equiv \frac{\parallel A{X}_j{A}^T+P_{l}B{B}^T{P}_l^T-E{X}_j{E}^T{\parallel }_F}{\parallel P_{l}B{B}^T{P}_l^T{\parallel }_F},$$

where $X_j$ is generated by LR-Smith or LR-ADI.

For LR-ADI, we first compute $k=20$ shift parameters by making use of Algorithm 3 and then reuse these parameters in a circular manner if the number of shift parameters is less than the number of iterations required to achieve the specified tolerance. In LR-ADI, we solve the shift linear systems by the $LU$ factorization of the corresponding coefficient matrices.

All the numerical results are obtained by performing calculations on an Intel Core i7-8650U with CPU 1.90 GHz and RAM 16 GB.

4.1. Example 1

In this example, the differential algebraic equation (DAE) is

$$\left\{\begin{array}{cc}\hfill {\widehat{E}}_{11}\dot{x}\left(t\right)& ={\widehat{A}}_{11}x\left(t\right)+{\widehat{A}}_{12}p\left(t\right)+{\widehat{B}}_{1}u\left(t\right),\hfill \\ \hfill 0& ={\widehat{A}}_{21}x\left(t\right)+{\widehat{B}}_{2}u\left(t\right)\hfill \end{array}\right.$$

(33)

with ${\widehat{E}}_{11}=I$, ${\widehat{A}}_{11}={\widehat{A}}_{11}^T$, and ${\widehat{A}}_{21}={\widehat{A}}_{12}^T$. It comes from the spatial discretization of the 2D instationary Stokes equation [37].

In order to obtain discrete-time equations, we first use a semi-explicit Euler and a semi-implicit Euler method [49] with a timestep size $\Delta t$ to discretize the differential equation of (33). This leads to two difference equations:

$$\begin{array}{ccc}\hfill {\widehat{E}}_{11}{x}_{k+1}& =& ({\widehat{E}}_{11}+\Delta t{\widehat{A}}_{11})x_k+\Delta t{\widehat{A}}_{12}{p}_k+\Delta t{\widehat{B}}_1{u}_k,\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill ({\widehat{E}}_{11}-\Delta t{\widehat{A}}_{11})x_{k+1}& =& {\widehat{E}}_{11}{x}_k+\Delta t{\widehat{A}}_{12}{p}_k+\Delta t{\widehat{B}}_1{u}_k.\hfill \end{array}$$

(35)

Then, by averaging (34) and (35) and also discretizing the algebraic equation of (33), we obtain the final difference-algebraic equations

$$\underset{E}{\underbrace{\left[\begin{array}{cc}{E}_{11}& 0\\ 0& 0\end{array}\right]}}\left[\begin{array}{c}{x}_{k+1}\\ p_{k+1}\end{array}\right]=\underset{A}{\underbrace{\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& 0\end{array}\right]}}\left[\begin{array}{c}{x}_k\\ p_k\end{array}\right]+B{u}_k,$$

(36)

where $E_{11}={\widehat{E}}_{11}-\frac{\Delta t}{2}{\widehat{A}}_{11}$, $A_{11}={\widehat{E}}_{11}+\frac{\Delta t}{2}{\widehat{A}}_{11}$, $A_{12}=A_{21}^T=\Delta t{\widehat{A}}_{12}$, and $B=\Delta t{[{\widehat{B}}_1^T,{\widehat{B}}_2^T]}^T$.

Note that $E,A$ in (36) are sparse and have special block structures. Using these structures, the projectors $P_l$ and $P_r$ can be formulated as

$$\begin{array}{ccc}\hfill P_l& =& \left[\begin{array}{cc}{\mathrm{\Pi }}_l& -{\mathrm{\Pi }}_l{A}_{11}\mathrm{\Psi }\\ 0& 0\end{array}\right],\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill P_r& =& \left[\begin{array}{cc}{\mathrm{\Pi }}_r& 0\\ -\mathrm{\Phi }{A}_{11}{\mathrm{\Pi }}_r& 0\end{array}\right],\hfill \end{array}$$

(38)

where

$$\begin{array}{ccc}\hfill {\mathrm{\Pi }}_l& =& I-A_{12}{\left(A_{21}{E}_{11}^{-1}{A}_{12}\right)}^{-1}{A}_{21}{E}_{11}^{-1}\hfill \\ \hfill {\mathrm{\Pi }}_r& =& I-E_{11}^{-1}{A}_{12}{\left(A_{21}{E}_{11}^{-1}{A}_{12}\right)}^{-1}{A}_{21}=E_{11}^{-1}{\mathrm{\Pi }}_l{E}_{11},\hfill \\ \hfill \mathrm{\Phi }& =& {\left(A_{21}{E}_{11}^{-1}{A}_{12}\right)}^{-1}{A}_{21}{E}_{11}^{-1},\hfill \\ \hfill \mathrm{\Psi }& =& E_{11}^{-1}{A}_{12}{\left(A_{21}{E}_{11}^{-1}{A}_{12}\right)}^{-1}.\hfill \end{array}$$

Moreover,

$$\begin{array}{ccc}\hfill E^{-}& =& {(P_{l}E+(I-P_l)A)}^{-1}{P}_l=\left[\begin{array}{cc}{\mathrm{\Pi }}_r{E}_{11}^{-1}& -{\mathrm{\Pi }}_r{E}_{11}^{-1}{A}_{11}\mathrm{\Psi }\\ -\mathrm{\Phi }{A}_{11}{\mathrm{\Pi }}_r{E}_{11}^{-1}& \mathrm{\Phi }{A}_{11}{\mathrm{\Pi }}_r{E}_{11}^{-1}{A}_{11}\mathrm{\Psi }\end{array}\right],\hfill \\ \hfill E^{-}A& =& \left[\begin{array}{cc}{\mathrm{\Pi }}_r{E}_{11}^{-1}{A}_{11}{\mathrm{\Pi }}_r& 0\\ -\mathrm{\Phi }{A}_{11}{\mathrm{\Pi }}_r{E}_{11}^{-1}{A}_{11}{\mathrm{\Pi }}_r& 0\end{array}\right],\hfill \end{array}$$

see, for example, refs. [28,37] for details.

In this example, the timestep size $\Delta t$ is taken to be $0.05$, and ${[{\widehat{B}}_1^T,{\widehat{B}}_2^T]}^T\in {\mathbb{R}}^{n\times 2}$ is a matrix with each element being 1, except the $(1,2)$-th element, which is 0.

We first test a medium-size problem of order $n=1280$. Figure 1 illustrates the sparsity structure of the matrix $A_{21}{E}_{11}^{-1}{A}_{12}$. This may show that for larger problems, it is expensive to compute the $LU$ factorization of $A_{21}{E}_{11}^{-1}{A}_{12}$. For this experiment, as well as for the larger problems later, the final relative residual accuracy for both the LR-ADI method and the LR-Smith method was set to ${10}^{-8}$. The convergence curves of LR-ADI and LR-Smith are depicted in Figure 2. The ADI method reaches a relative residual of 6.2472 × ${10}^{-9}$ after 13 steps of iteration, while the Smith method has a relative residual of 9.1324 × ${10}^{-9}$ after 75 steps. From Figure 2, it is clear that LR-Smith is much slower than LR-ADI with respect to the number of iterations. From Section 2, we know that the convergence factor of LR-Smith is the spectrum radius $\rho \left(E^{-}A\right)$. If the spectrum radius $\rho \left(E^{-}A\right)\approx 1$, LR-Smith will converge very slowly. Note that the spectrum radius $\rho \left(E^{-}A\right)$ is the absolute value of the largest finite eigenvalues (in modulus) of the pencil $\lambda E-A$. For this medium-size problem, the largest finite eigenvalue is 0.9554, which verifies the slow convergence of LR-Smith. However, from

Table 1. Numerical results for Example 1.

n	Method	Iter	Cpu Time	RRes
1280	LR-ADI	13	0.27	6.2472 × ${10}^{-9}$
	LR-Smith	75	0.34	9.1324 × ${10}^{-9}$
3604	LR-ADI	13	2.1	5.8859 × ${10}^{-9}$
	LR-Smith	179	5.0	9.6261 × ${10}^{-9}$
7700	LR-ADI	16	14.8	6.8133 × ${10}^{-9}$
	LR-Smith	335	40.8	9.9569 × ${10}^{-9}$
14,559	LR-ADI	22	86.8	2.6304 × ${10}^{-9}$
	LR-Smith	564	227.9	9.9038 × ${10}^{-9}$

, we can see that LR-Smith is only slightly more expensive, with respect to execution time, than LR-ADI. The reason is that linear systems with coefficient matrices A and $\overline{\mu }A-E$ must be solved in the latter method.

We also tested problems with larger dimensions, namely, $n=3604$, $n=7700$, and n = 14,559. All the numerical results are reported in Table 1. As expected, the execution time becomes increasingly larger for both LR-Smith and LR-ADI as the problem dimension expands. Moreover, it is clear that both the number of iterations and the cpu time of LR-ADI are better than those of LR-Smith. However, it is interesting that the number of iterations for LR-Smith increases much more dramatically than for LR-ADI. This is the reason that LR-ADI is more favorable with respect to the cpu time for problems with larger dimensions. For illustration, we also present the convergence curves of the two methods for dimensions $n=7700$ and n = 14,559, respectively, in Figure 3 and Figure 4.

4.2. Example 2

In the second example, the DAE (33) is used to describe a holonomically constrained damped mass–spring system with g masses [28,37]. The system matrices have the following structures:

$$\begin{array}{ccc}\hfill {\widehat{E}}_{11}=\left[\begin{array}{cc}I& 0\\ 0& M\end{array}\right],& & {\widehat{A}}_{11}=\left[\begin{array}{cc}0& I\\ K& D\end{array}\right],\hfill \\ \hfill {\widehat{A}}_{12}=\left[\begin{array}{c}0\\ -N^T\end{array}\right],& & {\widehat{A}}_{21}=\left[\begin{array}{cc}N& 0\end{array}\right].\hfill \end{array}$$

Similarly, we discretize the DAE by the same method as in the first example to obtain difference-algebraic equations, in which the matrices $E,A$ have the same block structures as in the first example.

By setting $g=2000,6000$, 10,000, we obtain three DAEs of order $n=2g+1=4001$, 12,001, 200,001. The timestep size $\Delta t$ is taken to be $0.1$. All the elements of ${[{\widehat{B}}_1^T,{\widehat{B}}_2^T]}^T\in {\mathbb{R}}^{n\times 1}$ are set to 0, except the $(g+1)$-th element, which is 1.

In Figure 5, we first illustrate the convergence curves of LR-Smith and LR-ADI for the problem with the dimension $n=4001$. Oddly enough, we can obviously see that the relative residual of LR-Smith does not monotonously decrease. This is due to the rounding error, which can make the relative residual strangely increase if the convergence factor is almost equal to 1. LR-Smith did not converge even if it had run 1000 steps, which cost 45.6 s. However, LR-ADI, with 21 steps of iteration, had a relative residual of 1.2886 × ${10}^{-9}$ and a cpu time of 1.8 s.

Since LR-Smith does not converge for this example, we only report our numerical results of LR-ADI for all problems of different dimensions ($n=4001$, 12,001, 20,001) in

Table 2. Numerical results for Example 2.

n	Method	Iter	Cpu Time	RRes
4001	LR-ADI	21	1.8	1.2886 × ${10}^{-9}$
12,001	LR-ADI	22	5.3	3.0284 × ${10}^{-9}$
20,001	LR-ADI	26	9.1	2.3404 × ${10}^{-9}$

. We observe that LR-ADI is very fast for this example and, in particular, that the number of iterations does not increase considerably with the expansion of the problem dimension. Moreover, we also notice that the execution time for this example is much less than that for the first example. The reason is that, in the second example, $A_{21}{E}_{11}^{-1}{A}_{12}$ is a number, and $E,A$ are sparser than those in the first one, which makes it cheaper to compute the $LU$ factorization of $A_{21}{E}_{11}^{-1}{A}_{12},A$, and $\overline{\mu }A-E$ in this example. Finally, we show the convergence curve of LR-ADI for the problem with the dimension n = 20,001 in Figure 6.

5. Conclusions

We have proposed a low-rank Smith method and a low-rank ADI method for the solutions of large-scale projected Stein equations. Although the Smith iteration and ADI iteration are very common in the field of efficient numerical solutions of linear matrix equations such as Lyapunov equations and Sylvester equations, this is the first time that they are adapted to numerically solve large-scale projected Stein equations, which arise in the balanced truncation model reduction of discrete-time descriptor systems. We also present a partially real version of the low-rank ADI method, as some of the shift parameters are imaginary. Our numerical experiments seem to show that the low-rank ADI method is more competitive than the low-rank Smith method with respect to the execution time and the number of iterations for large-scale problems if the convergence factor of the low-rank Smith method is large.