# The MGHSS for Solving Continuous Sylvester EquationAX+XB=C

Authors: Yu-Ye Feng; Qing-Biao Wu; Xue-Na Jing
Journal: Complexity
 (2021)

---

## Abstract

This paper proposes the modified generalization of the HSS (MGHSS) to solve a large and sparse continuous Sylvester equation, improving the efficiency and robustness. The analysis shows that the MGHSS converges to the unique solution ofAX + XB = C unconditionally. We also propose an inexact variant of the MGHSS (IMGHSS) and prove its convergence under certain conditions. Numerical experiments verify the efficiency of the proposed methods.


---

## Body

## 1. Introduction

This paper focuses on solving the continuous Sylvester equation defined as(1)AX+XB=C.Firstly, we assumeA, B, and C are large and sparse matrices, and A∈ℂm×m,B∈ℂn×n,andC∈ℂm×n, respectively; then, equation (1) is a large sparse equation. Then, we assume that both A and B are positive semidefinite, and at least one of them is positive definite and at least one of A and B is non-Hermitian.Under the above assumptions, it is sufficient to prove that equation (1) has a unique solution [1]. When B=A∗ and C is Hermitian, the continuous Sylvester equation (1) is a special case of the continuous Lyapunov equation [2]. A∗ indicates the conjugate transpose of A. The continuous Sylvester equation (1) has numerous applications in many fields, such as in control and system theory [3], signal processing [4] and image restoring [5], and stability of linear systems [6]. Many authors considered such a linear matrix equation problem and concentrated on accelerating the HSS iteration on the continuous Sylvester equation (1) [7–10], which was first proposed in [11].By using Kronecker Product, equation (1) is rewritten as the following linear equation:(2)Ax=c,where A=In⊗A+BT⊗Im, and x and c are the columns of X and C, respectively. Im represents the identity matrix whose order is m, and ⊗ is the Kronecker Product. However, when the size of the linear equation (2) is large, it is ill-conditioned to solve it directly.Before the appearance of the HSS, direct algorithms were usually used to solve the continuous Sylvester equations, such as Hessenberg–Schur and Bartels–Stewart methods [1, 12]. However, they were only applicable to small-sized continuous Sylvester equations. For large and sparse continuous Sylvester equations, iteration methods were used, such as the gradient-based algorithm [13–18]. Such an iteration method has been studied in recent years, taking advantage of the low-rank and sparsity of right-hand C in equation (1).In 2011, Bai proposed the HSS to solve large sparse continuous Sylvester equation [11]. Since then, many HSS-based iteration methods [19–23] have been widely studied and achieved certain results in solving the continuous Sylvester equation. In the same direction of the research, this paper presents a modified GHSS method to solve the continuous Sylvester equations. Besides, there are numerical research studies which focus on solving complex Sylvester matrix equation with large size, based on the HSS method for solving (1) which is proposed in [11]. Modified HSS (MHSS) iteration method [24] and preconditioned MHSS (PMHSS) method [9] for solving complex Sylvester matrix equation were presented, and then, the generalized MHSS (GMHSS) method [10] is also based on the MHSS iteration method by parameterizing it. In recent years, some neural network methods for time-varying complex Sylvester equation were proposed [25, 26]. Many methods are updated to solve various types of Sylvester equation. In this paper, we focus on solving continuous Sylvester equation with non-Hermitian and positive definite/semidefinite matrices.Firstly, the Hermitian and skew-Hermitian splitting method in equation (1) is used [27]. Let the Hermitian part of V be HV=1/2V+V∗ and the skew-Hermitian part of V be SV=1/2V−V∗.HSS [2]: the following equations are computed with an initial matrix x0 until xk, where k=0,1,2… is converged:(3)αI+Hxk+1/2=αI−Sxk+b,αI+Sxk+1=αI−Hxk+1/2+b.where α>0 is a constant.HSS for solving continuous Sylvester equations [11]: Xk+1∈ℂm×n is computed with an initial matrix x0 through the following equations until Xk+1k=0∞⊆ℂm×n satisfies the stopping criterion:(4)αI+HAXk+1/2+Xk+1/2βI+HB=αI−SAXk+XkβI−SB+C,αI+SAXk+1+Xk+1βI+SB=αI−HAXk+1/2+Xk+1/2βI−HB+C,where α>0andβ>0 are constants. The iteration matrix for solving continuous Sylvester equation is Mγ=γI+S−1γI−HγI+H−1γI−S and γ=α+β, and it converges to the exact solution: H=I⊗HA+HBT⊗I and S=I⊗SA+SBT⊗I.GHSS [28]: similar to HSS, the following equations are computed until xk is converged:(5)αI+Gxk+1/2=αI−K−Sxk+b,αI+S+Kxk+1=αI−Gxk+1/2+b,where α>0 is a constant.GHSS for solving continuous Sylvester equations [8]: in GHSS, to solve continuous Sylvester equations, Xk+1∈ℂm×n, where k=0,1,2… is computed through the following scheme until Xk+1k=0∞⊆ℂm×n satisfies the stopping criterion:(6)αI+GAXk+1/2+Xk+1/2βI+GB=αI−SA−KAXk+XkβI−SB−KB+C,αI+SA+KAXk+1+Xk+1βI+SB+KB=αI−GAXk+1/2+Xk+1/2βI−GB+C,where α>0andβ>0 are constants. The iteration matrix is ℳγ=γI+S+K−1γI−GγI+G−1γI−S−K and G=I⊗GA+GBT⊗I, K=I⊗KA+KBT⊗I, and S=I⊗SA+SBT⊗I. The convergence factor is given by the spectral radius ρℳγ of the matrix ℳγ, bounded by ρℳγ≤σγ:=maxλi∈λGAmaxμi∈λGBγ−λi+μi/γ+λi+μi, and the optimal parameter γ∗=λminλmax, where λmin=λminGA+λminGB and λmax=λmaxGA+λmaxGB.As an extension of those iterative schemes, this paper proposes the modified GHSS to solve the continuous Sylvester equations and proves its convergence. Section2 presents the detailed MGHSS and analyzes its convergence for the continuous Sylvester equation. IMGHSS iteration is described in Section 3. In Section 4, we take two examples into experiments from previous experiments in other HSS-based iteration methods. The results show that the proposed MGHSS is more effective in both the iteration and runtime. Section 5 concludes this paper with several thoughts.In the remainder of this paper, especially in the proof of the convergence property of MGHSS, we rewrite the continuous Sylvester equation (1) as the linear-vector form firstly. When the vector sequence ykk=0∞⊆ℂn2 converges the vector y∈ℂn2, it is equivalent as the corresponding matrix sequence Ykk=0∞⊆ℂn×n converges to the corresponding matrix Y∈ℂn×n, where the corresponding vector y is the concatenated columns of the correspondent matrix Y.


## 2. The Modified GHSS Iteration Method

This paper proposes a modified GHSS, and a direct solver can be used to solve each step of the inner iteration.Firstly, the GHSS splitsA and B into generalized Hermitian and skew-Hermitian parts [8]:(7)A=GA+SA+KA,B=GB+SB+KB,where SA and SB are the skew-Hermitian part of A and B, respectively, and HA=GA+KA and HB=GB+KB are the Hermitian part of A and B, respectively.With matrix splitting and GHSS [8], SA and SB are decomposed into two skew-Hermitian matrices: SA=RA+TA and SB=RB+TB. Then, A and B are rewritten:(8)A=αI+GA+RA−αI−TA−KA=αI+TA+KA−αI−GA−RA,B=βI+GB+RB−βI−TB−KB=βI+TB+KB−βI−GB−RB.Accordingly, equation (1) is rewritten in the matrix equation as follows:(9)αI+GA+RAX+XβI+GB+RB=αI−TA−KAX+XβI−TB−KB+C,αI+TA+KAX+XβI+TB+KB=αI−GA−RAX+XβI−GB−RB+C.It is known that no common eigenvalue exists betweenαI+GA+RA and −βI+GB+RA and between αI+TA+KA and −βI+TB+KB. Thus, there exist unique solutions in both two fixed-point matrix equations. Based on this, the MHSS iteration is conducted to solve equation (1).MGHSS:Xk+1∈ℂm×n is computed with an initial matrix x0 through equation (10). The process stops when Xk+1k=0∞⊆ℂm×n satisfies the stopping criterion:(10)αI+GA+RAXk+1/2+Xk+1/2βI+GB+RB=αI−TA−KAXk+XkβI−TB−KB+C,αI+TA+KAXk+1+Xk+1βI+TB+KB=αI−GA−RAXk+1/2+Xk+1/2βI−GB−RB+C.Lemma 1.
(see [11]). Let Mi and Ni denote two split of matrix A, where A=Mi−Nii=1,2. Then, a two-step iteration sequence Xk is defined as follows:(11)M1Xk+1/2B=N1XkB+C,M2Xk+1B=N2Xk+1/2B+C,where A,B,C∈ℂm×n and k=1,2…, x0 is an initial matrix. Then,(12)Xk+1=M2−1N2M1−1N1Xk+M2−1I+N2−1M1CB−1.This is rewritten as vector form as follows:(13)xk+1=I⊗M2−1N2M1−1N1xk+B−T⊗M2−1I+N2M1−1vecC.Furthermore, when the spectral radiusρI⊗M2−1N2M1−1N1<1, Xk converges to X∗∈ℂm×n for all X0∈ℂm×n.Lemma 2.
(see [29, 30]). Let A=H+S, where H=1/2A+A∗. When H is positive semidefinite and α⩾0,(14)αI−AαI+A−1≤1.WhenH is positive definite and α>0,(15)αI−AαI+A−1<1.Theorem 1.
LetA=G+K+R+T, where(16)G=In⊗GA+GBT⊗Im,K=In⊗KA+KBT⊗Im,R=In⊗RA+RBT⊗Im.T=In⊗TA+TBT⊗Im,γ=α+β,where R and T are skew-Hermitian matrices and G and K are symmetric positive semi-definite. The unique solution of (1) obtained by the MGHSS converges to the unique exact solution X∗∈ℂm×n when either G and K are symmetric positive definite.Proof.
By using the Kronecker product, we can rewrite (10) as follows:(17)γI+G+RvecXk+1/2=γI−T−KvecXk+vecC,γI+T+KvecXk+1=γI−G−RvecXk+1/2+vecC.
Then, it can reformulated as follows by eliminatingXk+1/2:(18)vecXk+1=ℳγvecXk+GγvecC,where(19)ℳγ=γI+T+K−1γI−G−RγI+G+R−1γI−T−K,Gγ=12αγI+T+K−1γI+G+R−1.
By a similar transformation of the components of the iteration matrixℳγ, we obtain(20)ℳγ^=γI−G−RγI+G+R−1γI−T−KγI+T+K−1.
DenoteA1=G+R and A2=T+K, then A1 and A2 are positive semidefinite, and clearly,(21)ℳγ^=γI−A1γI+A1−1γI−A2γI+A2−1.
From Lemma2, we have(22)γI−A1γI+A1−1≤1,γI−A2γI+A2−1≤1.
Respectively, ifA1 and A2 are positive definite, the above inequalities are strict. Thus, when either G and K is symmetric positive definite, we have(23)ρℳγ=ρℳ^γ≤γI−A1γI+A1−1γI−A2γI+A2−1<1,for any γ>0, completing the proof. Accordingly, the MGHSS unconditionally converges to the exact solution of equation (1).
From the results in Chapter 4 in [31], denote the inner product x,y=xTy. We know under the above definitions of A, G, K, R, T, A1, and A2, the following inequalities are satisfied:(24)Aiu,u≥mu,uandAiu,u≤Mu,u,i=1,2.
The upper bound on the spectral radius of the iteration matrixℳγ is minimized with the parameter γ, which is defined as follows:(25)γopt=1mM.It indicates that finding the optimal parameterγ is challenging but necessary because it relies on the spectral information of the selected iteration matrix.In the following section, the improved MGHSS, the inexact MGHSS (IMGHSS) is introduced.


## 3. The Inexact MGHSS Method

Unlike MGHSS (Section3) that solves the two fixed-point equations by direct algorithms, IMGHSS, presented in this section, iteratively solves the two subsystems. Similar to the IHSS [32, 33] for solving linear systems and IHSS for solving Sylvester equation [11], the process of the IMGHSS iteration scheme for solving a continuous Sylvester equation is as follows. Here, we denote ·F as Frobenius norm.IMGHSS:X0∈ℂm×n is an initial matrix. In the IMGHSS algorithm, the solution of equation (1) is derived as the following:(26)k=0,while (stopping condition = = false)(27)Rk=C−AXk−XkBis approximated to αI+GA+RAZk+ZkβI+GB+RA=Rk that is solved by the residual Pk=Rk−αI+GA+RAZk−βI+GB+RA of the iteration satisfies PkF≤ηkRkF:(28)Xk+1/2=Xk+Zk,(29)Rk+1/2=C−AXk+1/2−Xk+1/2B,is approximated to αI+TA+KAZk+1/2+Zk+1/2βI+TB+KB that is solved by the residual Qk+1/2=Rk+1/2−αI+TA+KAZk+1/2+Zk+1/2βI+TB+KB of the iteration which satisfies Qk+1/2F≤εkRk+1/2F:(30)Xk+1=Xk+1/2+Zk+1/2,(31)k=k+1,end.In the scheme of IMGHSS,εk and ηk are prescribed tolerances that control the accuracy of the inner iterations. In implements, the values of εkandηk do not necessarily decrease to zero when k increases so long as we choose suitable values of it, and we can also ensure the convergence of the IMGHSS.In [11], the convergence of the two-step iteration was explored. We analyze the convergence of IMGHSS as follows.Theorem 2.
LetXkk=0∞⊆ℂm×n denote an iteration sequence produced by IMGHSS and X∗∈ℂm×n denote the exact solution of equation (1). Then, under the assumption that the conditions of Theorem 1 are met, it holds(32)Xk+1−X∗S≤δγ+θρτkδ21+θηkδ2−1Xk−X∗S,k=1,2⋯,where the norm ·S is defined as follows:(33)YS=αI+TA+KAY+YβI+TB+KBF,for any matrix Y∈ℂm×n, and the constant ρ and θ are given by ρ=γI+T+KγI+G+R−12,θ=AγI+T+K−12, and δγ=δ1δ2, where δ1=γI−A1γI+A1−1 and δ2=γI−A2γI+A2−1.
In particular, whenδγ+θρεmaxδ21+θηmaxδ2−1<1, the iteration sequence Xkk=0∞⊆ℂm×n converges to X∗∈ℂm×n, where ηmax=maxkηk and εmax=maxkεk.Proof.
The IMGHSS can be rewritten in the matrix-vector form with the notations in Theorem1 and the Kronecker product as follows:(34)αI+G+Rzk=rk,xk+1/2=xk+zk,αI+T+Kzk+1/2=rk+1/2,xk+1=xk+1/2+zk+1/2,where rk=c−Axk and rk+1/2=c−Axk+1/2. zk is such that the residual pk=rk−αI+G+Rzk satisfies pk2≤ηkrk2, and zk+1/2 is such that the residual qk+1/2=rk+1/2−αI+T+Kzk+1/2 satisfies qk+1/22≤εk∥rk+1/22.
Equation (34) is the IMGHSS for solving (2), with A=G+K+R+T. Then, based on Theorem 2 in [11], we have(35)xk+1−x∗S≤δγ+θρεkδ21+θηkδ2−1xk−x∗S,k=1,2⋯,where the constants are given by(36)ρ=γI+T+KγI+G+R−12,θ=AγI+T+K−12,and δγ=γI−A1γI+A1−1γI−A2γI+A2−1andδ2=γI−A2γI+A2−1.
For a vectory∈ℂm×n that consists of the concatenated columns of Y, yS=αI+TA+KAY+YβI+TB+KBF.
Thus, the following is obtained:(37)Xk+1−X∗S≤δγ+θρεkδ21+θηkδ2−1Xk−X∗S,k=1,2⋯.
The proof is completed.According to Theorem2, it is important to choose a suitable value of the tolerance εkandηk to control the IMGHSS’s convergence. Still, analyzing the optimal tolerances is challenging.


## 4. Numerical Analysis

The feasibility and efficiency of the MGHSS are verified in several examples in this section. The proposed method was compared with other methods in terms of the number of iteration steps (nis) and the computational time (t [sec]). The numerical analysis was conducted in Matlab on Intel dual-core CPU (2.5 GHz) and 8 GB RAM. Zero matrix was used as an initial guess, and the termination condition is defined as(38)C−AXk−XkBFCF≤10−6.Example 1.
The continuous Sylvester equation (1) with m=n is considered, and the matrices are as follows:(39)M=2.3−1−12.3−1⋱⋱⋱−12.3−1−12.3,N=00.5−0.500.5⋱⋱⋱−0.500.5−0.50,A=B=M+2rN+100/n+12I, where I represents the identity matrix.
For the practical iteration parameters of all those iteration methods, we takeα=β. Also, all the subproblems are exactly solved by the Bartels–Stewart method [3] in each step of the HSS, GHSS, and MGHSS. Tables 1 and 2 compare HSS and MGHSS in solving the continuous Sylvester equation (1) in terms of nis and tsec. The optimal values of αexpandβexpαexp=βexp were analyzed in Tables 3 and 4, respectively, for HSS/MGHSS and GHSS/MGHSS.
According to Table3 and Table 4, as the rank n of equation (2) is incremented, αexp and βexp of the HSS, GHSS, and MGHSS are all decreased. In Figure 1, the logarithm versus iteration of the HSS, GHSS, and MGHSS methods (n=128) are shown in (a) and (b) when r=0.1 and r=1.0, respectively. It shows the efficiency of the MGHSS method.Table 1
Comparison of HSS and MGHSS in terms ofnis and t.
rHSSMGHSS0.010.11.00.010.11.0nistnistnistnistnistnistn = 8120.002130.007110.007110.001110.00190.004n = 16200.012230.013150.017190.003200.004160.009n = 32440.076430.080320.055430.014410.074230.028n = 64861.556730.543440.332720.622680.494290.077n = 12817010.11411913.956691.547784.469794.870512.435n = 25633566.75817534.7348314.3787827.8897929.0827411.475Table 2
Comparison of GHSS and MGHSS in terms ofnis and t.
rGHSSMGHSS0.010.11.00.010.11.0nistnistnistnistnistnistn = 8110.003120.002100.004110.001110.00190.004n = 16190.012150.05170.048190.003200.004160.009n = 32160.121170.175290.065430.014410.074230.028n = 64270.864310.842410.834720.622680.494290.077n = 128442.754456.3475519.437784.469794.870512.435n = 2566830.5846731.462167125.6437827.8897929.0827411.475Table 3
The analysis of the optimal valuesαexpandβexp for HSS and MGHSS.
rHSSMGHSS0.010.11.00.010.11.0n = 82.752.552.452.652.482.61n = 161.261.261.481.251.32.53n = 320.640.651.150.620.631.78n = 640.320.320.930.270.31.71n = 1280.160.20.750.060.050.81n = 2560.830.140.650.010.040.71Table 4
The analysis of the optimal valuesαexpandβexp for GHSS and MGHSS.
rGHSSMGHSS0.010.11.00.010.11.0n = 82.052.12.352.652.482.61n = 160.260.361.281.251.32.53n = 320.340.351.050.620.631.78n = 640.320.360.830.270.31.71n = 1280.260.30.750.060.050.81n = 2560.30.340.550.010.040.71Figure 1
The logarithmRXk versus iteration of the HSS, GHSS, and MGHSS methods.
(a)(b)Example 2.
The continuous Sylvester equation (1) with m = n and the matrices
A=diag1,2,…,n+rLT and B=2−pI+diag1,2,…,n+rLT+2−pL, where L is the strictly lower triangular matrix with ones in the lower triangle part and p is a problem parameter specified in practical computations.
Table5 shows that the MGHSS outperforms the GHSS and HSS in solving the continuous Sylvester equation. In Table 6, the continuous Sylvester equation in Example 2 are solved by the IMGHSS and MGHSS iteration methods and the results show that IMGHSS is much better than the MGHSS. Here, we set εk=ηk=0,01,k=0,1,2,…, and use the ADI method as the inner iteration scheme.Table 5
Comparisons of HSS, GHSS, and MGHSS in terms ofnis and t.
HSSGHSSMGHSSnistnistnistn = 8170.014170.022170.013n = 16270.035210.087250.027n = 32350.608290.453340.048n = 64503.244343.231470.303n = 1287119.5394310.543643.87n = 25610030.765020.432789.095Table 6
Comparison of IMGHSS and MGHSS in terms ofnis and t.
IMGHSSMGHSSnistnistn = 880.011170.013n = 16190.021250.027n = 32240.029340.048n = 64270.256470.303n = 128331.123643.87


## 5. Conclusions

HSS-based methods have been widely used to solve the continuous Sylvester equations. In this paper, a modified generalization of the HSS method (MGHSS) is proposed. A preconditioner can also be taken for all of the generalizations of the HSS, although many researchers concentrated on the studies of the relations between parameters and the convergence property of each. Furthermore, we establish the IMGHSS as an efficient solver. The convergence of the MGHSS and IMGHSS were analyzed. Also, the efficiency and robustness of the proposed method were verified in several examples.



---

*Source: 6615728-2021-02-11.xml*