# An Alternating Direction Implicit Method for Solving Projected Generalized Continuous-Time Sylvester Equations

Authors: Yujian Zhou; Liang Bao; Yiqin Lin
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

We present the generalized low-rank alternating direction implicit method and the low-rank cyclic Smith method to solve projected generalized continuous-time Sylvester equations with low-rank right-hand sides. Such equations arise in control theory including the computation of inner products andℍ2 norms and the model reduction based on balanced truncation for descriptor systems. The requirements of these methods are moderate with respect to both computational cost and memory. Numerical experiments presented in this paper show the effectiveness of the proposed methods.


---

## Body

## 1. Introduction

In this paper we consider the projected generalized continuous-time Sylvester equation of the form(1)AXC+DXB+Pl,1EPr,2=0,X=Pr,1XPl,2,
where A,D∈ℝn×n, B,C∈ℝm×m, E∈ℝn×m, and X∈ℝn×m are the sought-after solution. Here, Pr,1 and Pr,2 are the spectral projectors onto the right deflating subspaces corresponding to the finite eigenvalues of the pencils λD-A and λC-B, respectively; Qr,1=I-Pr,1 and Qr,2=I-Pr,2 are the spectral projectors onto the right deflating subspaces corresponding to the eigenvalue at infinity. The spectral projectors onto the left deflating subspaces corresponding to the finite eigenvalues of λD-A and λC-B are denoted by Pl,1 and Pl,2, respectively, while Ql,1=I-Pl,1 and Ql,2=I-Pl,2 are the spectral projectors onto the left deflating subspaces corresponding to the eigenvalue at infinity. We assume that the matrices C and D are singular, but the pencils λD-A and λC-B are regular; that is, det(λD-A) and det(λC-B) are not identically zero. Under the assumption, the pencils λD-A and λC-B have the Weierstrass canonical forms: there exist nonsingular n×n matrices W1,T1 and m×m matrices W2,T2 such that
(2)D=W1[I00N(A)]T1,A=W1[J(A)00I]T1,C=W2[I00N(B)]T2,B=W2[J(B)00I]T2,
where J(A), J(B), N(A), and N(B) are block diagonal matrices with each diagonal block being a Jordan block. The eigenvalues of J(A) and J(B) are the finite eigenvalues of the pencils λD-A and λC-B, respectively. N(A) and N(B) correspond to the eigenvalue at infinity. The index γ of nilpotency of N(A) is called the index of the pencil λD-A. Using (2), Pl,1, Pl,2, Pr,1, and Pr,2 can be expressed as
(3)Pl,1=W1[I000]W1-1,Pl,2=W2[I000]W2-1,Pr,1=T1-1[I000]T1,Pr,2=T2-1[I000]T2.A number of numerical solution methods have been proposed for the standard/generalized Lyapunov and Sylvester equations. Two classical direct methods are the Bartels-Stewart method [1] and the Hammarling method [2]. These methods need to compute the real Schur forms/generalized real Schur forms of the underlying matrices/matrix pencils by means of the QR/QZ algorithm [3]. Besides direct methods, we mention, among several iterative methods, the Smith method [4], the alternating direction implicit iteration (ADI) method [5, 6], the Smith(l) method [7], the low-rank Smith methods [7–9], the Cholesky factor-alternating direction implicit method [10–13], and the (generalized) matrix sign function method [14–25].Asm=n, B=AT, and C=DT, (1) is referred to as the projected generalized continuous-time Lyapunov equation, which arises in stability analysis and control design problems for descriptor systems including the characterization of controllability and observability properties, balanced truncation model order reduction, determining the minimal and balanced realizations, and computing H2 and Hankel norms; see [26–31] and the references therein. If the pencil λD-A is c-stable, that is, all its finite eigenvalues have negative real part, then the projected generalized Lyapunov equation has a unique solution for each E, and if, additionally, E is symmetric and positive semidefinite, then the solution X is symmetric and positive semidefinite; see, for example, [32] for details. Recently, several numerical methods have been proposed in the literature for solving the projected generalized Lyapunov equation. In [33], two direct methods, the generalized Bartels-Stewart method and the generalized Hammarling method, were proposed for the projected generalized Lyapunov equation. The generalized Hammarling method is designed to obtain the Cholesky factor of the solution. These two methods are based on the generalized real Schur form of the pencil λD-A and require 𝒪(n3) flops and 𝒪(n2) memory. Iterative methods to solve the projected generalized Lyapunov equation have also been proposed. Stykel [31] extended the ADI method and the Smith method to the projected equation. Moreover, their low-rank versions were also presented, which could be used to compute low-rank approximations to the solution. These methods are especially suitable for large sparse equations with low-rank E. Another iterative method for the projected generalized Lyapunov equation is the modified generalized matrix sign function method [25]. Unlike the classical generalized matrix sign function method, the variant converges quadratically independent of the index of the underlying matrix pencil; see [25] for more details.The projected generalized Sylvester equation (1) also plays an important role in control theory for descriptor systems.LetH(s)=G(sD-A)-1F and H~(s)=G~(sC-B)-1F~ be two c-stable systems. As shown in [34], H(s) and H~(s) can be decomposed into
(4)H(s)=Hsp(s)+P(s),H~(s)=H~sp(s)+P~(s),
where
(5)Hsp(s)=G(sD-A)-1Pl,1F,P(s)=G(I-Pr,1)(sD-A)-1(I-Pl,1)F,H~sp(s)=G~(sC-B)-1Pl,2F,P~(s)=G~(I-Pr,2)(sC-B)-1(I-Pl,2)F~.
Here, Hsp(s),H~sp(s) are called the strictly proper parts of H(s),H~(s), while P(s),P~(s) are called the polynomial parts of H(s),H~(s), respectively. Define the ℍ𝕃2 inner product of H(s) and H~(s) by
(6)〈H(s),H~(s)〉ℍ𝕃2=〈Hsp(s),H~sp(s)〉ℍ2+〈P(s),P~(s)〉𝕃2=12π∫-∞∞‍trace(Hsp(iω)¯H~spT(iω))dω+12π∫02π‍trace(P(eiω)¯P~T(eiω))dω,
where the notation trace(·) denotes the trace of a matrix. Then, the ℍ𝕃2 norm of H(s) induced by the ℍ𝕃2 inner product is
(7)∥H(s)∥ℍ𝕃2=∥Hsp(s)∥ℍ22+∥P(s)∥𝕃22=〈Hsp(s),Hsp(s)〉ℍ2+〈P(s),P(s)〉𝕃2.
For F,GT∈ℝn,  F~,G~T∈ℝm, it can be shown that 〈Hsp(s),H~sp(s)〉ℍ2 could be expressed as
(8)〈Hsp(s),H~sp(s)〉ℍ2=trace(GXF~),
where X is the solution of the projected generalized continuous-time Sylvester equation (1) with E=FG~. Therefore, the norm ∥Hsp(s)∥ℍ2 can be computed by
(9)∥Hsp(s)∥ℍ2=trace(GXF),
where X is the solution of the following projected generalized continuous-time Sylvester equation:
(10)AXD+DXA+Pl,1FGPr,1=0,X=Pr,1XPl,1.We now present the application of the projected generalized continuous-time Sylvester equation in model reduction of c-stable descriptor systems. Let𝒢pc and 𝒢po denote the proper controllability and observability Gramians of the continuous-time descriptor system H(s)=G(sD-A)-1F, respectively. For the definitions of these Gramians and their applications in model reduction of H(s), the reader is referred to [30]. The proper Gramians 𝒢pc and 𝒢po are the unique symmetric, positive semidefinite solutions of of the projected generalized continuous-time controllability and observability Lyapunov equations associated with the system H(s), respectively.There exists another kind of Gramians, called cross-Gramians; for a standard square system, see [35]. We now extend the definition to square descriptor systems. The proper cross-Gramian for a c-stable square descriptor system H(s)=G(sD-A)-1F is the solution X of the projected generalized continuous-time Sylvester equation (10). It can been shown that the proper cross-Gramian X satisfies
(11)(XD)2=𝒢pcDT𝒢poD.
By making use of the relation above and a similar approach as proposed in [35], we can design a cross-Gramian-based model reduction method for the square descriptor system H(s).In this paper, we propose two iterative methods for solving the projected generalized Sylvester equation (1). This work presented here is an extension of [31]. Numerical experiments show the effectiveness of the proposed methods.We note that the perturbation analysis for (1) has been presented in [36]. The following theorem, for example, [37], gives sufficient conditions for the existence, uniqueness, and analytic formula of the solution of the projected generalized continuous-time Sylvester equation (1).Theorem 1.
LetλD-A and λC-B be regular pencils with finite eigenvalues {μ1,μ2,…,μn1} and {ν1,ν2,…,νm1} counted according to their multiplicities, respectively. Then, the projected generalized continuous-time Sylvester equation (1) has a unique solution for every E if μi+νj≠0 for any i=1,2,…,n1 and j=1,2,…,m1.
Moreover, ifλD-A and λC-B are c-stable, that is, all their finite eigenvalues have negative real part, then X can be expressed as
(12)X=12π∫-∞∞‍(-iωD-A)-1Pl,1EPr,2(iωC-B)-1dω.Throughout this paper, we adopt the following notations. The square identity and zero matrices are denoted byI and 0, respectively. The spaces of m×n real matrices are denoted by ℝm×n. The 2-norm and the Frobenius matrix norm are denoted by ∥·∥2 and ∥·∥F, respectively. The superscript “·T” denotes the transposition of a vector or a matrix. ρ(A) is the spectral radius of square matrix A, and κ(A)=∥A∥2∥A-1∥2 is the spectral condition number of A. The open left and right half-plane are denoted by ℂ- and ℂ+, respectively. We will also adopt MATLAB-like convention to access the entries of vectors and matrices. For a matrix X, X(i,j) is X’s (i,j)th entry; X’s submatrices X(k:l,i:j), X(k:l,:), and X(:,i:j) consist of intersections of row k to row l and column i to column j, row k to row l, and column i to column j, respectively.The remainder of the paper is organized as follows. In Section2, we propose the generalized low-rank alternating direction implicit method for the solution of (1). In Section 3, we discuss the choice of shift parameters. In Section 4, we show how to solve the projected Sylvester equations by the low-rank cyclic Smith method. Section 5 is devoted to some numerical tests. Some concluding remarks are given in the last section.


## 2. Alternating Direction Implicit Method

The alternating direction implicit (ADI) method was first introduced by Peaceman and Rachford [38] to solve linear systems arising from the discretization of elliptic boundary value problems and then used in [5–7, 12, 13, 31] to solve Lyapunov or Sylvester matrix equations. In this section, we consider how to generalize the ADI method for iteratively solving (1).We always assume that the pencilsλD-A and λC-B are c-stable; that is, all their finite eigenvalues have negative real part. Hence, the matrices A and B are nonsingular. Multiplying the first equation in (1) on the left by A-1 and on the right by B-1, we get the following projected standard Sylvester equation:
(13)(A-1D)X+X(CB-1)+A-1Pl,1EPr,2B-1=0,X=Pr,1XPl,2.
The iterates Xk of the ADI iteration for (13) are usually generated by the alternating solution of two linear systems with multiple right-hand sides:
(14)(A-1D-βkI)Xk-(1/2)=Xk-1(CB-1+βkI)+A-1Pl,1EPr,2B-1,Xk(CB-1+αkI)=(A-1D-αkI)Xk-(1/2)-A-1Pl,1EPr,2B-1,
where X0=0 and the shift parameters {αi}i=1k and {βi}i=1k are elements of ℂ- and ℂ+, respectively. These two equations are equivalent to the following single iteration step:
(15)Xk=(A-1D-βkI)-1(A-1D-αkI)Xk-1×(CB-1+βkI)(CB-1+αkI)-1+(βk-αk)(A-1D-βkI)-1A-1Pl,1EPr,2B-1×(CB-1+αkI)-1.
We have from
(16)Pr,1(A-1D-βkI)-1=(A-1D-βkI)-1Pr,1,Pr,1(A-1D-αkI)=(A-1D-αkI)Pr,1,(CB-1+αkI)-1Pl,2=Pl,2(CB-1+αkI)-1,(CB-1+βkI)Pl,2=Pl,2(CB-1+βkI)
that Xk=Pr,1XkPl,2, that is, the second equation in (1) and (13) is satisfied exactly. We can rewrite iteration (15) as
(17)Xk=(D-βkA)-1(D-αkA)Xk-1×(C+βkB)(C+αkB)-1+(βk-αk)(D-βkA)-1Pl,1EPr,2(C+αkB)-1.
Let X denote the exact solution of (1). Then it is easily verified that the error matrix X-Xk obeys the recursion:
(18)X-Xk=(D-βkA)-1(D-αkA)(X-Xk-1)×(C+βkB)(C+αkB)-1=⋯=𝒜kXℬk,
where
(19)𝒜k=Pr,1(D-βkA)-1(D-αkA)⋯(D-β1A)-1(D-α1A),ℬk=(C+β1B)(C+α1B)-1⋯(C+βkB)(C+αkB)-1Pl,2.
By using the Weierstrass canonical form (2)–(3), we have
(20)𝒜k=T1-1[Jk(A)000]T1,ℬk=W2[Jk(B)000]W2-1
with
(21)Jk(A)=(I-βkJ(A))-1(I-αkJ(A))⋯(I-β1J(A))-1(I-α1J(A)),Jk(B)=(I+β1J(B))(I+α1J(B))-1⋯(I+βkJ(B))(I+αkJ(B))-1.
This implies that if {αj}j=1k contains all finite eigenvalues (multiple eigenvalues counted by their algebraic multiplicities) of A-1D or if {βj}j=1k contains all finite eigenvalues of -B-1C, then X-Xk≡0. This is due to the Cayley-Hamilton theorem [39], which states that p(A)≡0 for A’s characteristic polynomial p(λ)=det(λI-A) and q(B)≡0 for q(λ)=det(λI-B). From (18), we can see that the parameters {αi}i=1k and {βi}i=1k should be chosen to achieve minαi,βi∥𝒜k∥2·∥ℬk∥2; we will discuss the details of the choice of parameters in the next section.About thekth approximate solution Xk of the ADI method, by using (18)–(20), we have the following estimate.Theorem 2.
Assume that the pencilsλD-A and λC-B are in Weierstrass form (2) with J(A) and J(B) being diagonal. Then
(22)∥X-Xk∥2⩽κ(T1)κ(W2)ρ(𝒜k)ρ(ℬk)∥X∥2.Note that the solutionX is explicitly computed by the ADI iteration (17), so the storage requirement is 𝒪(mn). One should notice that in many cases the storage requirement is the limiting factor rather than the amount of computation. We note that low-rank schemes are the only existing methods that can effectively solve large-scale Lyapunov/Sylvester equations. Assume that the low-rank right-hand side E has the factored form E=FG with F∈ℝn×r and G∈ℝr×m. Instead of explicitly forming the solution X, the low-rank method computes and stores approximate solutions in low-rank factored form. If the numerical rank l of X is much smaller than min{m,n}, that is, l≪min{m,n}, then the storage is reduced from 𝒪(mn) to 𝒪(ml) or 𝒪(nl).The key idea in the low-rank version of the ADI iteration is to rewrite the iterationXk in (17) as an outer product:
(23)Xk=YkHkZk.
This is always possible since starting with the initial guess X0=0n×m. The low-rank alternating direction implicit (LR-ADI) method is based on (17). Replacing Xk-1 with Yk-1Hk-1Zk-1, (17) can be reformulated in terms of the low-rank factors as
(24)Xk=[(D-βkA)-1Pl,1F(D-βkA)-1(D-αkA)Yk-1]×[(βk-αk)IHk-1]×[GPr,2(C+αkB)-1Zk-1(C+βkB)(C+αkB)-1]≡YkHkZk,
where
(25)Yk=[(D-βkA)-1Pl,1F(D-βkA)-1(D-αkA)Yk-1],Hk=[(βk-αk)IHk-1],Zk=[GPr,2(C+αkB)-1Zk-1(C+βkB)(C+αkB)-1].
From the fact that Y0, H0, and Z0 are all zero matrices, it can be seen that Yk is n×kr, Hk is kr×kr, and Zk is kr×m. Thus the rank of Xk is no more than kr. Since the order of the ADI parameters {αi}i=1k and {βi}i=1k is not important, the ordering of {αi}i=1k and {βi}i=1k can be reversed. As shown in [12], we have the following iterative scheme:
(26)Yk=[Y(1)Y(2)⋯Y(k)],withY(1)=(D-β1A)-1Pl,1F,Y(i+1)=(D-βi+1A)-1(D-αiA)Y(i)(i=1,2,…,k-1),Zk=[Z(1)Z(2)⋯Z(k)],withZ(1)=GPr,2(C+α1B)-1,Z(i+1)=Z(i)(C+βiB)(C+αi+1B)-1(i=1,2,…,k-1),
and so we have
(27)Xk=YkHkZk
with
(28)Hk=diag((β1-α1)Ir,…,(βk-αk)Ir).
The algorithm is summarized as in Algorithm 1.Algorithm 1:The generalized LR-ADI method for the projected generalized Sylvester equation (1) with E=FG.
Input:  A,D∈ℝn×n, B,C∈ℝm×m, F∈ℝn×r and G∈ℝr×m with λD-A
andλC-B are c-stable. The number of ADI steps k and the ADI shifts
{α1,α2,…,αk} and {β1,β2,…,βk} computed by Algorithm 2.
Output: Y, H and Z such that X=YHZ is an approximate solution of the projected
generalized Sylvester equation (1) with E=FG.
(1) LetY(:,1:r)=(D-β1A)-1Pl,1F, Z(1:r,:)=GPr,2(C+α1B)-1 and H(1:r,1:r)=(β1-α1)Ir.
(2)    Fori=1,2,…
Y(:,ir+1:(i+1)r)=(D-βi+1A)-1(D-αiA)Y(:,(i-1)r+1:ir),
Z(ir+1:(i+1)r,:)=Z((i-1)r+1:ir,:)(C+βiB)(C+αi+1B)-1,
H(ir+1:(i+1)r,ir+1:(i+1)r)=(βi+1-αi+1)Ir.
End ForWe make a few comments on Algorithm1.(1)
Since the computation of the Weierstrass canonical form is sensitive under small perturbations, we should make use of the generalized real Schur factorization to compute the spectral projectors; see, for example, [33]. For large-scale problems, the computation of the spectral projectors by the generalized real Schur factorization may be very expensive. However, in some applications the spectral projectors can be expressed in explicit form by using the special block structure of the matrices A, D, B, and C; see numerical examples in this paper or [31].(2)
If the number of shift parametersk is smaller than the number of iterations required to obtain a prescribed tolerance, then we reuse these parameters in a cyclic manner.


## 3. Choice of Shift Parameters

The selection of good parameters is vitally important to the successful application of the generalized LR-ADI iteration. The rate of convergence is dominated by spectral radii of matrices𝒜k and ℬk given by (20). The minimization of these spectral radii with respect to the parameters {αj}j=1k and {βj}j=1k leads to the generalized ADI minimax problem
(29)min{α1,…,αk}∈ℂ-{β1,…,βk}∈ℂ+maxx∈𝔼fy∈𝔽f∏j=1k‍|(1-αjx)(1-βjy)(1-βjx)(1-αjy)|,
where 𝔼f and 𝔽f denote two sets of finite eigenvalues of the pencils λD-A and λC+B, respectively. In practice since the eigenvalues of the pencils λD-A and λC+B are unknown and expensive to compute, 𝔼f and 𝔽f are often replaced by domains that contain the eigenvalues of λD-A and λC+B, respectively. Noting that A and B are nonsingular, this problem is equivalent to find {αj}j=1k and {βj}j=1k such that
(30)min{α1,…,αk}∈ℂ-{β1,…,βk}∈ℂ+maxx∈𝔼∖{0}y∈𝔽∖{0}∏j=1k‍|(x-αj)(y-βj)(x-βj)(y-αj)|,
where 𝔼 and 𝔽 denote spectrums of the matrices A-1D and -B-1C, respectively.To compute the suboptimal ADI shift parameters for the standard Lyapunov equation withB=AT and C=D=I, a heuristic algorithm has been proposed in [7]. It chooses suboptimal ADI parameters from a set R, which is taken to be the union of the Ritz values of A and the reciprocals of the Ritz values of A-1, obtained by two Arnoldi processes, with A and A-1. As shown in [7], the Ritz values obtained by the Arnoldi process with A tend to be located near the “outer” eigenvalues, that is, the eigenvalues near the convex hull of the spectrum. In particular, the eigenvalues of large magnitude are usually approximated well. In contrast, they are generally poor approximations to the eigenvalues near the origin. Therefore, one computes the reciprocals of the Ritz values obtained by the Arnoldi process with A-1 to approximate the eigenvalues near the origin.In [12], this idea has been extended to the case for Sylvester equations. A heuristic procedure which is easy to implement has been proposed in [12]. It can also be naturally extended to the generalized problem (30). Here, we need to compute the largest and smallest (in modulus) nonzero approximate eigenvalues of A-1D and -B-1C, respectively. Noting that C and D are assumed to be singular, inverses of A-1D and -B-1C do not exist. In [31], Stykel proposed a strategy to overcome this difficulty. Define
(31)P1=Pr,1(DPr,1-AQr,1)-1=(Pl,1D-Ql,1A)-1Pl,1=T1-1[I000]W1-1,P2=Pr,2(CPr,2-BQr,2)-1=(Pl,2C-Ql,2B)-1Pl,2=T2-1[I000]W2-1,
where T1, T2, W1, and W2 are the transformation matrices as in (2). A simple calculation gives that
(32)P1A=T1-1[J(A)000]T1,-P2B=T2-1[-J(B)000]T2.
Then it is clear that the reciprocals of the smallest nonzero eigenvalues of A-1D and -B-1C are the largest eigenvalues of P1A and -P2B, respectively. Thus, we can run two Arnoldi processes with the matrices P1A and -P2B to compute the smallest nonzero eigenvalues of A-1D and -B-1C, respectively. In some applications, similar to the projectors Pr,1, Pl,1, Pr,2 and Pl,2, the matrices P1, and P2 can be also obtained in explicit form by making use of the special block structure of the matrices A, D, B, and C. The algorithm for choosing {αj}j=1k and {βj}j=1k is summarized as in Algorithm 2.Algorithm 2:ADI parameters by Ritz values.
Input:  A,D∈ℝn×n, B,C∈ℝm×m with λD-A and λC-B c-stable, F∈ℝn×r
andG∈ℝr×m.
Output: ADI parameters {α1,α2,…,αk} and {β1,β2,…,βk}.
(1)  Run Arnoldi process withA-1D on F to give the set 𝔼A+ of Ritz values,
(2)  Run Arnoldi process withP1A on F to give the set 𝔼A- of Ritz values,
(3)𝔼=𝔼A+∪(1/𝔼A-),
(4)  Run Arnoldi process with-B-1C on GT to give the set 𝔽B+ of Ritz values,
(5)  Run Arnoldi process with-P2B on GT to give the set 𝔽B- of Ritz values,
(6)𝔽=(𝔽B+)∪(1/𝔽B-),
(7) Set{α1,β1}=argminα∈𝔼β∈𝔽maxx∈𝔼y∈𝔽|(x-α)(y-β)(x-β)(y-α)|,
(8)    Fori=2,3,…,k
Set{αi,βi}=argminα∈𝔼′β∈𝔽′maxx∈𝔼y∈𝔽|(x-α)(y-β)(x-β)(y-α)|∏j-1i-1|(x-αj)(y-βj)(x-βj)(y-αj)|,
where𝔼′ is 𝔼 with α1,…,αi-1 deleted, and similarly for 𝔽′,
End ForFor more details about Algorithm2, the interesting reader is referred to [12].


## 4. Smith Method

The Smith method [4] is derived from the projected discrete-time Sylvester equation
(33)𝒜Xℬ-X+(β-α)Pr,1ℱ𝒢Pl,2=0,X=Pr,1XPl,2,
with
(34)𝒜=(D-βA)-1(D-αA),ℬ=(C+βB)(C+αB)-1,ℱ=(D-βA)-1F,𝒢=G(C+αB)-1,
which is equivalent to (1) for any two parameters α∈ℂ- and β∈ℂ+. Under this assumption, the sequence {Xi}i=0∞ generated by the iteration
(35)X0=(β-α)Pr,1ℱ𝒢Pl,2,Xk=𝒜Xk-1ℬ+(β-α)Pr,1ℱ𝒢Pl,2
converges to the solution X.If all of the shifts in the generalized ADI iteration (17) are constant, that is, αk=α, βk=β (k=1,2,…), then the generalized ADI iteration reduces to the Smith method. In [7], Penzl illustrated that the ADI iteration with a single shift (Smith method) converges very slowly, while a moderate increase in the number of shifts s accelerates the convergence nicely. However, it is also observed that the speed of convergence is hardly improved by a further increase of s; see Table  2.1 in [7]. These observations lead to the idea of the cyclic Smith(s) iteration, a special case of ADI where s different shifts are used in a cyclic manner.The low-rank scheme based on the Smith(s) iteration was also introduced in [7]. The method is called the cyclic low-rank Smith method (LR-Smith(s)) and is a special case of LR-ADI where s number of shifts are reused in a cyclic manner. This idea can be generalized for (1). The generalized LR-Smith method consists of two steps. First the matrices Ys, Hs, and Zs are obtained by an s step generalized LR-ADI iteration. Then one solves the discrete-time Sylvester equation
(36)𝒜sXℬs-X+YsHsZs=0,
where 𝒜s and ℬs are as in (19) with k=s, respectively.The generalized LR-Smith method for solving the projected generalized Sylvester (1) is described as in Algorithm 3.Algorithm 3:The generalized LR-Smith(s) method for the projected generalized Sylvester equation (1) with E=FG.
Input: A,D∈ℝn×n, B,C∈ℝm×m, F∈ℝn×r and G∈ℝr×m. λD-A and λC-B
are c-stable.
Output: Y, H and Z such that X=YHZ is an approximate solution of the projected
generalized Sylvester equation (1) with E=FG.
(1)  ComputeYs, Hs and Zs by using Algorithm 1 with k=s.
(2)  SetY=Y-=Ys, H=H-=Hs and Z=Z-=Zs.
(3) Fori=1,2,…
Forj=1,2,…,s
Y-=(D-βjA)-1(D-αjA)Y-,
Z-=Z-(C+βjB)(C+αjB)-1,
End For
Y=[YPr,1Y-],Z=[ZZ-Pl,2],
H=[H00H-],
End For


## 5. Numerical Experiments

In this section, we present some numerical examples to illustrate the performance of the generalized LR-ADI method for the projected generalized Sylvester equation (1). Let LR-ADI and LR-Smith(s) denote Algorithms 1 and 3, respectively. In the following examples, we compare the numerical behavior of LR-ADI with LR-Smith(s) with respect to the number of iterations (Iter for short), CPU time (CPU for short), and the relative residuals ∥Rk∥F/∥R0∥F (Res for short). All numerical experiments are performed on a PC with the usual double precision, where the floating point relative accuracy is 2.22×10-16.The stopping criterion for both methods is(37)∥Rk∥F∥R0∥F=∥AXkC+DXkB+Pl,1EPr,2∥F∥Pl,1EPr,2∥F⩽10-12.Let mvp denote the number of matrix vector products and let lss denote the number of linear system solvers at every iteration for the two methods. Table1 has a rough count of the expenses of LR-ADI and LR-Smith(s) at every iteration. Only the major expenses are considered. From Table 1, we can see that LR-Smith(s) needs more mvp and lss than LR-ADI at every iteration, so LR-ADI is more efficient than LR-Smith(s) for the computation cost at every iteration.Table 1
mvp and lss at every iteration step.
LR-ADI
LR-Smith(s)
mvp forD-αA
r
r
s
2
mvp for(C+βB)T
r
r
s
2
lss forD-βA
r
r
s
2
lss for(C+αB)T
r
r
s
2Example 1.
For the first experiment, we consider the 2D instationary Stokes equation that describes the flow of an incompressible fluid in a domain. The spatial discretization of this equation by the finite difference method on a uniform staggered grid leads to the descriptor systems(38)Dx˙(t)=Ax(t)+Fu(t),y(t)=G~x(t),(39)Cx˙(t)=Bx(t)+F~u(t),y(t)=Gx(t).
The matrix coefficients in (38) are given by
(40)D=[I000]∈ℝn×n,A=[A11A12A210]∈ℝn×n.
These matrices are sparse and have a special block structure. Using this structure, the projectors Pl,1 and Pr,1 onto the left and right deflating subspaces of the pencil λD-A can be computed as
(41)Pl,1=[Π-ΠA11A12(A21A12)-100],Pr,1=[Π0-(A21A12)-1A21A11Π0],
where Π=I-A12(A21A12)-1A21 is the orthogonal projector onto the kernel of A21 along the image of A12; see [40]. The matrices A12 and A21 have full rank, and the pencil λD-A is of index 2. Analogously, we can get the projectors Pl,2 and Pr,2 of the pencil λB-C. In our experiment the state space dimensions of the problems are n=2295 and m=1975, respectively. The matrix E in (1) is E=FG with F∈ℝ2295×1 and G∈ℝ1×1975.
See Figure1 for a graph of the convergence. The results in Table 2 clearly indicate that the LR-ADI method is more efficient than LR-Smith(3) for this example.Table 2
LR-ADI,k=12 versus LR-Smith(3).
LR-ADI
LR-Smith(3)
Iter
30
12
mvp forD-αA
30
108
mvp for(C+βB)T
30
108
lss forD-βA
30
108
lss for(C+αB)T
30
108
CPU
10.68
12.21
Res
8.93
×
1
0
-
13
3.72
×
1
0
-
13Figure 1
Comparison for Example1.Example 2.
For the second experiment, we consider a holonomically constrained damped mass-spring system withg masses as in [25]. The ith mass is connected to the (i+1)th mass by a spring and a damper and also to the ground by another spring and damper. Moreover, the first mass is connected to the last one by a rigid bar, and it can be influenced by a control. The vibration of this system is described by the descriptor system (38) with the matrices
(42)D=[I000M10000],A=[0I0K1H1-N1TN100],
where M1∈ℝg×g is the symmetric positive definite mass matrix, K1∈ℝg×g is the stiffness matrix, H1∈ℝg×g is the damping matrix, and N1 is the matrix of constraints. If N1∈ℝ1×g has full row rank, the pencil λD-A is of index 3, and the spectral projectors Pl,1 and Pr,1 can be computed as
(43)Pl,1=[Π0-ΠM1-1H1Q-ΠTH1(I-Π)ΠT-ΠT(K1+H1ΠM1-1H1)Q000],Pr,1=[Π00-ΠM1-1H1(I-Π)Π0-QT(K1Π-H1ΠM1-1H1(I-Π))QTH1Π0].
Here Q=M1-1N1T(N1M1-1N1T)-1 and Π=I-M1-1N1T(N1M1-1N1T)-1N1=I-QN1 are a projector onto the kernel of N1 along the image of M1-1N1T; see [41] for details. The matrix F has the form F=[e1,e2,eg-1]∈ℝn×3, where ej denotes the jth column of the identity matrix I2g+1. Analogously, we can get the projectors Pl,2 and Pr,2 of the pencil λB-C. In this experiment the state space dimensions of the problems are n=1261 and m=1161, respectively. The matrix E in (1) is E=FG with F∈ℝ1261×3 and G∈ℝ3×1161.
The computational results were reported in Figure2 and Table 3. We note that the iteration steps of the LR-Smith(6) method is 21 while the steps of LR-ADI is 42, but the CPU time of LR-Smith(6) method is much more than that of LR-ADI method. This is because the computation cost at every iteration of two methods is very different.Table 3
LR-ADI,k=12 versus LR-Smith(6).
LR-ADI
LR-Smith(6)
Iter
42
21
mvp forD-αA
126
2268
mvp for(C+βB)T
126
2268
lss forD-βA
126
2268
lss for(C+αB)T
126
2268
CPU
30.80
44.13
Res
7.16
×
1
0
-
13
7.69
×
1
0
-
13Figure 2
Comparison for Example2.


## 6. Conclusions

In this paper, we have proposed the generalized low-rank alternating direction implicit method and the generalized low-rank cyclic Smith method to solve the projected generalized Sylvester equations. Numerical experiments presented in this paper show the effectiveness of the proposed methods.



---

*Source: 387854-2013-09-26.xml*