A Robust Hermitian and Skew-Hermitian Based Multiplicative Splitting Iterative Method for the Continuous Sylvester Equation

Abstract

For solving the continuous Sylvester equation, a class of Hermitian and skew-Hermitian based multiplicative splitting iteration methods is presented. We consider two symmetric positive definite splittings for each coefficient matrix of the continuous Sylvester equations, and it can be equivalently written as two multiplicative splitting matrix equations. When both coefficient matrices in the continuous Sylvester equation are (non-symmetric) positive semi-definite, and at least one of them is positive definite, we can choose Hermitian and skew-Hermitian (HS) splittings of matrices A and B in the first equation, and the splitting of the Jacobi iterations for matrices A and B in the second equation in the multiplicative splitting iteration method. Convergence conditions of this method are studied in depth, and numerical experiments show the efficiency of this method. Moreover, by numerical computation, we show that multiplicative splitting can be used as a splitting preconditioner and induce accurate, robust and effective preconditioned Krylov subspace iteration methods for solving the continuous Sylvester equation.

1. Introduction

Matrix equations arise in a number of problems of scientific computations and engineering applications, such as in control theory [1,2], model reduction [3], signal processing [4] and image processing [5], and many researchers focus on the matrix equations [6,7,8,9,10,11]. Nowadays, the continuous Sylvester equation is possibly the most famous and the most broadly employed linear matrix equation [5,6,11,12,13,14,15,16,17,18,19]. It is given as

$$AX+XB=C,$$

(1)

where $A\in {\mathbb{R}}^{n\times n}$, $B\in {\mathbb{R}}^{m\times m}$ and $C\in {\mathbb{R}}^{n\times m}$ are defined matrices and $X\in {\mathbb{R}}^{n\times m}$ is an unknown matrix. A Lyapunov equation is a special case with $m=n,B=A^T$ and $C=C^T$. Here and in the subsequent sections, $W^T$ is used to denote the transpose of the matrix $W\in {\mathbb{R}}^{n\times n}$. The continuous Sylvester Equation (1) has a unique solution if and only if A and $-B$ have no common eigenvalues, which will be assumed throughout this paper.

The continuous Sylvester equation has been studied extensively in recent decades as a key pipeline and has received wide applications in areas such as control theory, differential systems and signal processing. As a result, a number of algorithms focusing on the dimension reduction of the equation have been proposed [20]. The resulting algorithms can be extended to more general forms of the continuous algebraic Lyapunov equations and extended Sylvester equations, as well as to the Markov jump system. For the continuous Sylvester equation, various iterative methods are given in the literature. For instance, splitting the matrix equation to find an approximate solution of the continuous Sylvester equation has shown proper convergence, and is simpler than the corresponding multiplicative splitting and thusg has relatively lower computational complexity [20,21,22,23,24]. Another class of solutions focuses on overcoming the lack of stability via the combination of two split equations and corresponding new iterative methods to approximate the solution of the continuous Sylvester equation via the modified Richardson technique [23,25,26].

As a summary, if the equation can be expressed as a simple, symmetric, regular singular equation, it can exploit the mixed structure to solve the equation by finding a specific transformation to the standard form or degenerative decomposition for large-scale experiments. Toward a solution, different techniques like regular deflating linearization, the Hamiltonian boundary value method, dealing with fundamental subspaces which do not relate to eigenvalues of the M-matrix, utilizing standard polynomial eigenproblem algorithms without computation of any special splitting and implementing the Hessenberg matrix to convert the homogeneous equivalent matrix to a special block structure have been reported in the literature [25,27,28,29].

In general, the dimensions of A and B may be orders of magnitude different, and this fact is key in selecting the most appropriate numerical solution strategy [11]. For solving general Sylvester equations of small size, we use some methods classified such as direct methods. Among these direct methods are the Bartels–Stewart [30] and the Hessenberg–Schur [31] methods, which consist of transforming coefficient matrices A and B into triangular or Hessenberg form by an orthogonal similarity transformation and then solving the resulting system directly by a back-substitution process. When the coefficient matrices A and B are large and sparse, iterative methods are often the methods of choice for solving the Sylvester equation (1) efficiently and accurately. Many iterative methods have been developed for solving matrix equations, such as the alternating direction implicit (ADI) method [32], the Krylov subspace-based algorithms [13,33,34], the Hermitian and skew-Hermitian splitting (HSS) method, the inexact variant of HSS (IHSS) iteration method [6] and the nested splitting conjugate gradient (NSCG) method [17,35] and the nested splitting CGNR (NS-CGNR) method [18].

In order to study the numerical methods, we often rewrite the continuous Sylvester equation (1) as th following linear system of equations

$$\mathcal{A}x=c,$$

(2)

where the matrix $\mathcal{A}$ is of dimension $nm\times nm$ and is given by

$$\mathcal{A}=I_m\otimes A+B^T\otimes I_n,$$

(3)

where ⊗ denotes the Kronecker product $(A\otimes B=\left[a_{ij}B\right])$ and

$$\begin{array}{c}c=vec\left(C\right)={(c_{11},c_{21},\cdots ,c_{n1},c_{12},c_{22},\cdots ,c_{n2},\cdots ,c_{nm})}^T\hfill \\ x=vec\left(X\right)={(x_{11},x_{21},\cdots ,x_{n1},x_{12},x_{22},\cdots ,x_{n2},\cdots ,x_{nm})}^T.\hfill \end{array}$$

Of course, this is quite expensive and a numerically poor way to determine the solution X of the continuous Sylvester equation (1), as the linear system of Equation (2) is costly to solve and can be ill-conditioned.

Now, we recall some necessary notations and useful results, which will be used in the following section. In this paper, we use ${\lambda \left(M\right),\left|\right|M\left|\right|}_2$, ${\left|\right|M\left|\right|}_F$ and $I_n$ to denote the eigenvalue, the spectral norm, the Frobenius norm of a matrix $M\in {\mathbb{R}}^{n\times n}$ and the identity matrix with dimension n, respectively. Note that ${\left|\right|.\left|\right|}_2$ is also used to represent the 2-norm of a vector. For nonsingular matrix $\mathcal{B}$, we denote by $\kappa \left(\mathcal{B}\right)={\left|\right|\mathcal{B}\left|\right|}_2\left|\right|{\mathcal{B}}^{-1}{\left|\right|}_2$ its spectral condition number, and for a symmetric a positive definite matrix $\mathcal{B}$, we define the $\left|\right|·{\left|\right|}_{\mathcal{B}}$ norm of a vector $x\in {\mathbb{R}}^n$ as ${\left|\right|x\left|\right|}_{\mathcal{B}}=\sqrt{x^H\mathcal{B}x}$. Then, the induced $\left|\right|·{\left|\right|}_{\mathcal{B}}$ norm of a matrix $M\in {\mathbb{R}}^{n\times n}$ is defined as ${\left|\right|M\left|\right|}_{\mathcal{B}}=\left|\right|{\mathcal{B}}^{\frac{1}{2}}M{\mathcal{B}}^{-\frac{1}{2}}{\left|\right|}_2$. In addition, it holds that ${\left|\right|Mx\left|\right|}_{\mathcal{B}}\le {\left|\right|M\left|\right|}_{\mathcal{B}}{\left|\right|x\left|\right|}_{\mathcal{B}}$, ${\left|\right|M\left|\right|}_{\mathcal{B}}\le \sqrt{\kappa \left(\mathcal{B}\right)}{\left|\right|M\left|\right|}_2$ and ${\left|\right|I\left|\right|}_{\mathcal{B}}=1$, where I is the identity matrix. For any matrices $A=\left[a_{ij}\right]$ and $B=\left[b_{ij}\right]$, $A\otimes B$ denotes the Kronecker product, defined as $A\otimes B=\left[a_{ij}B\right]$. For the matrix $X=(x_1,x_2,\cdots ,x_m)\in {\mathbb{R}}^{n\times m}$, $vec\left(X\right)$ denotes the $vec$ operator, defined as $vec\left(X\right)={(x_1^T,x_2^T,\cdots ,x_m^T)}^T$. Moreover, for a matrix $M\in {\mathbb{R}}^{n\times n}$ and the vector $vec\left(M\right)\in {\mathbb{R}}^{nm}$, we have ${\left|\right|M\left|\right|}_F={\left|\right|vec\left(M\right)\left|\right|}_2$.

For matrix $\mathcal{A}\in {\mathbb{R}}^{n\times n}$, $\mathcal{A}=\mathcal{B}-\mathcal{C}$ is called a splitting of the matrix $\mathcal{A}$ if $\mathcal{B}$ is nonsingular. This splitting is a convergent splitting if $\rho \left({\mathcal{B}}^{-1}\mathcal{C}\right)<1$ and a contractive splitting if $\left|\right|{\mathcal{B}}^{-1}\mathcal{C}\left|\right|<1$ for some matrix norm.

The reminder of this paper is organized as follows. Section 2 presents our main contribution. In other words, the multiplicative splitting iteration (MSI) method for the continuous Sylvester equation and its convergence properties are studied deeply. Section 3 is devoted to an extensive numerical experiments with full comparison with other state-of-the-art methods in the literature. In Section 4, we address some challenges and suggestions for future work. Finally, we present our conclusions in Section 5.

2. Multiplicative Splitting Iterations

2.1. Traditional MSI Method

Consider the linear system of Equation (2). Let $\mathcal{A}={\mathcal{M}}_i-{\mathcal{N}}_i$$(i=1,2)$ be two splittings of the coefficient matrix $\mathcal{A}$. The MSI method for solving the system of linear Equation (2) is defined as follows [36]:

MSI method for linear system of equations:

• Given an initial guess $x^{\left(0\right)}\in {\mathbb{R}}^n$,
• For $k=1,2,\cdots$ until convergence, do
  • $u^{(k+1)}={\mathcal{M}}_1^{-1}{\mathcal{N}}_1{x}^{\left(k\right)}+{\mathcal{M}}_1^{-1}c$
  • $x^{(k+1)}={\mathcal{M}}_2^{-1}{\mathcal{N}}_2{u}^{(k+1)}+{\mathcal{M}}_2^{-1}c$
• end

The MSI method can be equivalently written in the form

$$x^{(k+1)}={\mathcal{T}}_{msi}{x}^{\left(k\right)}+{\mathcal{G}}_{msi}c,k=0,1,2,\cdots$$

where ${\mathcal{T}}_{msi}={\mathcal{M}}_2^{-1}{\mathcal{N}}_2{\mathcal{M}}_1^{-1}{\mathcal{N}}_1$ and ${\mathcal{G}}_{msi}={\mathcal{M}}_2^{-1}{\mathcal{N}}_2{\mathcal{M}}_1^{-1}+{\mathcal{M}}_2^{-1}$. See [36] for more details.

2.2. MSI Method for the Sylvester Equation

Based on the MSI method proposed in [36], we obtain the MSI method for the continuous Sylvester equation. Let $A=M_i-N_i$ and $B=P_i-Q_i,(i=1,2)$ be two splittings of the matrices A and B, such that $M_i$ and $P_i,(i=1,2)$ are symmetric positive definite. The continuous Sylvester equation (1) can be equivalently written as the multiplicative splitting matrix equations

$$\left\{\begin{array}{c}{M}_{1}U+U{P}_1=N_{1}X+X{Q}_1+C\hfill \\ M_{2}X+X{P}_2=N_{2}U+U{Q}_2+C\hfill \end{array}\right.$$

Under the assumption that $M_i$ and $P_i,(i=1,2)$ are symmetric positive definite, we easily know that there is no common eigenvalues between the matrices $M_i$ and $-P_i$, $(i=1,2)$, so that this two multiplicative splitting matrix equations have unique solutions for all given right-hand-side matrices.

Now, based on the above observations, we can establish the following multiplicative splitting iterations for solving the continuous Sylvester Equation (1):

MSI method for Sylvester equation:

• Given an initial guess $X^{\left(0\right)}\in {\mathbb{R}}^{m\times n}$,
• For $k=1,2,\cdots$ until convergence, do
  • Solve $M_1{U}^{(k+1)}+U^{(k+1)}{P}_1=N_1{X}^{\left(k\right)}+X^{\left(k\right)}{Q}_1+C$
  • Solve $M_2{X}^{(k+1)}+X^{(k+1)}{P}_2=N_2{U}^{(k+1)}+U^{(k+1)}{Q}_2+C$
• end

In special cases, when both coefficient matrices A and B in the Sylvester equation (1) are (non-symmetric) positive semi-definite, and at least one of them is positive definite, in the first equation in the MSI method we can choose Hermitian and skew-Hermitian (HS) splittings of matrices A and B, i.e., $A=H_A-S_A$ and $B=H_B-S_B$, where $H_A,S_A,H_B,S_B$ are the Hermitian and skew-Hermitian parts of A and B, respectively. Also, in the second equation in the MSI method, we consider the splitting of the Jacobi iterations [37] for matrices A and B, i.e., $A=D_A-N_A$ and $B=D_B-N_B$, where $D_A,N_A,D_B,N_B$ are the diagonal and non-diagonal parts of A and B, respectively. Therefore, we can rewrite this method as follows:

• Given an initial guess $X^{\left(0\right)}\in {\mathbb{R}}^{m\times n}$,
• For $k=1,2.\cdots$ until convergence, do
  • Solve system $H_A{U}^{(k+1)}+U^{(k+1)}{H}_B=S_A{X}^{\left(k\right)}+X^{\left(k\right)}{S}_B+C$
  • Solve system $D_A{X}^{(k+1)}+X^{(k+1)}{D}_B=N_A{U}^{(k+1)}+U^{(k+1)}{N}_B+C$
• end

Achieving two Sylvester equations that we can easily solve is our motivation for choosing these splittings. This is because the first system can be solved by the Sylvester conjugate gradient method [38], and the following routine can be used for direct solution of the second system:

Directly solution of matrix equation ${\mathit{D}}_{\mathit{A}}\mathit{X}+{\mathit{XD}}_{\mathit{B}}=\mathit{F}$:

• For $i=1:n$
  • For $j=1:m$
    • $x_{ij}=\frac{f_{ij}}{a_{ii}+b_{jj}}$
  • end
• end

where $a_{ii}$ and $b_{jj}$ are the diagonal elements of matrices A and B, respectively. Moreover, $F=N_{A}U+U{N}_B+C$, and U is the solution of the first equation in the MSI method for Sylvester equation.

2.3. Using Multiplicative Splitting as a Preconditioner

Given the fact that any matrix splitting can naturally induce a splitting preconditioner for the Krylov subspace methods (see [39]) in Section 3, by numerical computation, we show that multiplicative splitting can be used as a splitting preconditioner to induce accurate, robust and effective preconditioned Krylov subspace iteration methods for solving the continuous Sylvester equation.

2.4. Convergence Analysis

In the subsequent solution, we need the following lemmas.

Lemma 1

([36]). Let $B,C\in {\mathbb{R}}^{n\times n}$ be two Hermitian matrices. Then, $BC=CB$ if and only if B and C have a common set of orthonormal eigenvectors.

Lemma 2

([40]). Let $\mathcal{A}\in {\mathbb{R}}^{n\times n}$ be a symmetric positive definite matrix. Then, for all $x\in {\mathbb{R}}^n$, we have $\left|\right|{\mathcal{A}}^{\frac{1}{2}}{x\left|\right|}_2={\left|\right|x\left|\right|}_{\mathcal{A}}$ and

$$\sqrt{{\lambda }_{min}\left(\mathcal{A}\right)}{\left|\right|x\left|\right|}_{\mathcal{A}}\le {\left|\right|\mathcal{A}x\left|\right|}_2\le \sqrt{{\lambda }_{max}\left(\mathcal{A}\right)}{\left|\right|x\left|\right|}_{\mathcal{A}}.$$

Lemma 3

([41]). Suppose that $A,B\in {\mathbb{R}}^{n\times n}$ are two Hermitian matrices, and denote the minimum and the maximum eigenvalues of a matrix M with ${\lambda }_{min}\left(M\right)$ and ${\lambda }_{max}\left(M\right)$, respectively. Then,

$$\begin{array}{c}{\lambda }_{max}(A+B)\le {\lambda }_{max}\left(A\right)+{\lambda }_{max}\left(B\right),\hfill \\ {\lambda }_{min}(A+B)\ge {\lambda }_{min}\left(A\right)+{\lambda }_{min}\left(B\right).\hfill \end{array}$$

Lemma 4

([41]). Let $A,B\in {\mathbb{C}}^{n\times n}$, and λ and μ be the eigenvalues of A and B, and x and y be the corresponding eigenvectors, respectively. Then $\lambda \mu$ is an eigenvalue of $A\otimes B$ corresponding to the eigenvector $x\otimes y$.

Lemma 5.

Suppose that $\mathcal{A}=\mathcal{M}-\mathcal{N}$ is a splitting such that $\mathcal{M}$ is symmetric positive definite, with $\mathcal{M}=I_m\otimes M+P^T\otimes I_n$ and $\mathcal{N}=I_m\otimes N+Q^T\otimes I_n$. If

$${\theta }^3\frac{max\left|\lambda \right(N\left)\right|+max\left|\lambda \right(Q\left)\right|}{{\lambda }_{min}\left(M\right)+{\lambda }_{min}\left(P\right)}<1,$$

where $\theta =\sqrt{\frac{{\lambda }_{max}\left(M\right)+{\lambda }_{max}\left(P\right)}{{\lambda }_{min}\left(M\right)+{\lambda }_{min}\left(P\right)}}$, then $\left|\right|{\mathcal{M}}^{-1}\mathcal{N}{\left|\right|}_{\mathcal{M}}<1$.

Proof. 

By Lemmas 3 and 4, we have

$${\left|\right|\mathcal{M}\left|\right|}_2={\lambda }_{max}\left(\mathcal{M}\right)\ge {\lambda }_{min}\left(\mathcal{M}\right)\ge {\lambda }_{min}\left(M\right)+{\lambda }_{min}\left(P\right),$$

and

$${\left|\right|\mathcal{N}\left|\right|}_2=\underset{\lambda \in \Lambda \left(\mathcal{N}\right)}{max}\left|\lambda \left(\mathcal{N}\right)\right|\le max\left|\lambda \left(N\right)\right|+max\left|\lambda \left(Q\right)\right|.$$

Therefore, it follows that

$$\begin{array}{ccc}\hfill \left|\right|{\mathcal{M}}^{-1}\mathcal{N}{\left|\right|}_{\mathcal{M}}& \le & \sqrt{\kappa \left(\mathcal{M}\right)}\left|\right|{\mathcal{M}}^{-1}\mathcal{N}{\left|\right|}_2\hfill \\ & \le & \sqrt{\kappa \left(\mathcal{M}\right)}\left|\right|{\mathcal{M}}^{-1}{\left|\right|}_2{\left|\right|\mathcal{N}\left|\right|}_2\hfill \\ & \le & {\left(\kappa \left(\mathcal{M}\right)\right)}^{\frac{3}{2}}\frac{{\left|\right|\mathcal{N}\left|\right|}_2}{{\left|\right|\mathcal{M}\left|\right|}_2}\hfill \\ & \le & {\left(\kappa \left(\mathcal{M}\right)\right)}^{\frac{3}{2}}\frac{max\left|\lambda \right(N\left)\right|+max\left|\lambda \right(Q\left)\right|}{{\lambda }_{min}\left(M\right)+{\lambda }_{min}\left(P\right)}.\hfill \end{array}$$

Again, the use of Lemmas 3 and 4 implies that

$$\sqrt{\kappa \left(\mathcal{M}\right)}=\sqrt{\frac{{\lambda }_{max}\left(\mathcal{M}\right)}{{\lambda }_{min}\left(\mathcal{M}\right)}}\le \sqrt{\frac{{\lambda }_{max}\left(M\right)+{\lambda }_{max}\left(P\right)}{{\lambda }_{min}\left(M\right)+{\lambda }_{min}\left(P\right)}}=\theta .$$

(4)

So, we can write

$$\left|\right|{\mathcal{M}}^{-1}\mathcal{N}{\left|\right|}_{\mathcal{M}}\le {\theta }^3\frac{max\left|\lambda \right(N\left)\right|+max\left|\lambda \right(Q\left)\right|}{{\lambda }_{min}\left(M\right)+{\lambda }_{min}\left(P\right)}.$$

(5)

This clearly proves the lemma. □

Theorem 1.

Let $A\in {\mathbb{R}}^{n\times n}$ and $B\in {\mathbb{R}}^{m\times m}$ and consider two splittings $A=M_i-N_i$ and $B=P_i-Q_i(i=1,2)$ such that $M_i$ and $P_i,(i=1,2)$ are symmetric positive definite. Denote by $\mathcal{A}={\mathcal{M}}_i-{\mathcal{N}}_i(i=1,2)$ with ${\mathcal{M}}_i=I_m\otimes M_i+P_i^T\otimes I_n$ and ${\mathcal{N}}_i=I_m\otimes N_i+Q_i^T\otimes I_n(i=1,2)$, and assume that ${\mathcal{M}}_1{\mathcal{A}}^{-1}$ and ${\mathcal{M}}_2{\mathcal{A}}^{-1}$ are Hermitian matrices and ${\mathcal{M}}_1{\mathcal{A}}^{-1}{\mathcal{M}}_2={\mathcal{M}}_2{\mathcal{A}}^{-1}{\mathcal{M}}_1$. Then, the MSI method is convergent if ${\varrho }_1{\varrho }_2<1$, where

$${\varrho }_i={\theta }_i^3\frac{max|\lambda \left(N_i\right)|+max|\lambda \left(Q_i\right)|}{{\lambda }_{min}\left(M_i\right)+{\lambda }_{min}\left(P_i\right)},and{\theta }_i=\sqrt{\frac{{\lambda }_{max}\left(M_i\right)+{\lambda }_{max}\left(P_i\right)}{{\lambda }_{min}\left(M_i\right)+{\lambda }_{min}\left(P_i\right)}},(i=1,2).$$

Proof. 

By making use of the Kronecker product, we can rewrite the above-described MSI method in the following matrix–vector form:

$$\left\{\begin{array}{c}(I_m\otimes M_1+P_1^T\otimes I_n)u^{(k+1)}=(I_m\otimes N_1+Q_1^T\otimes I_n)x^{\left(k\right)}+c\hfill \\ (I_m\otimes M_2+P_2^T\otimes I_n)x^{(k+1)}=(I_m\otimes N_2+Q_2^T\otimes I_n)u^{(k+1)}+c\hfill \end{array}\right.$$

which can be arranged equivalently as

$$\left\{\begin{array}{c}{\mathcal{M}}_1{u}^{(k+1)}={\mathcal{N}}_1{x}^{\left(k\right)}+c\hfill \\ {\mathcal{M}}_2{x}^{(k+1)}={\mathcal{N}}_2{u}^{(k+1)}+c\hfill \end{array}\right.$$

which can be obtained by the following iteration method:

$$\left\{\begin{array}{c}{u}^{(k+1)}={\mathcal{M}}_1^{-1}{\mathcal{N}}_1{x}^{\left(k\right)}+{\mathcal{M}}_1^{-1}c\hfill \\ x^{(k+1)}={\mathcal{M}}_2^{-1}{\mathcal{N}}_2{u}^{(k+1)}+{\mathcal{M}}_2^{-1}c\hfill \end{array}\right.$$

(6)

Evidently, the above iteration scheme is the MSI method [36] for solving the system of linear Equation (2) with $\mathcal{A}={\mathcal{M}}_i-{\mathcal{N}}_i(i=1,2)$. The MSI iteration (6) can be neatly expressed as a stationary fixed-point iteration as follows:

$$x^{(k+1)}=\mathcal{T}{x}^{\left(k\right)}+\mathcal{G}c$$

with $\mathcal{T}={\mathcal{M}}_2^{-1}{\mathcal{N}}_2{\mathcal{M}}_1^{-1}{\mathcal{N}}_1$ and $\mathcal{G}={\mathcal{M}}_2^{-1}{\mathcal{N}}_2{\mathcal{M}}_1^{-1}+{\mathcal{M}}_2^{-1}$.

Because ${\mathcal{M}}_1{\mathcal{A}}^{-1}{\mathcal{M}}_2={\mathcal{M}}_2{\mathcal{A}}^{-1}{\mathcal{M}}_1$ is equivalent to the two matrices ${\mathcal{M}}_1{\mathcal{A}}^{-1}$ and ${\mathcal{M}}_2{\mathcal{A}}^{-1}$ being commutative, according to Lemma 1, we know that ${\mathcal{M}}_1{\mathcal{A}}^{-1}$ and ${\mathcal{M}}_2{\mathcal{A}}^{-1}$ have a common set of orthonormal eigenvectors. That is say, there exists a unitary matrix $\mathcal{Q}\in {\mathbb{R}}^{nm\times nm}$ and two diagonal matrices ${\Lambda }_i=\mathit{diag}({\lambda }_1^{\left(i\right)},{\lambda }_2^{\left(i\right)},\cdots ,{\lambda }_{nm}^{\left(i\right)}),i=1,2$, such that $\mathcal{Q}{\mathcal{M}}_i^{-1}\mathcal{A}{\mathcal{Q}}^{*}={\Lambda }_i,i=1,2$. Noticing that

$$\begin{array}{cc}\hfill \mathcal{T}& ={\mathcal{M}}_2^{-1}{\mathcal{N}}_2{\mathcal{M}}_1^{-1}{\mathcal{N}}_1\hfill \\ & ={\mathcal{M}}_2^{-1}({\mathcal{M}}_2-\mathcal{A}){\mathcal{M}}_1^{-1}({\mathcal{M}}_1-\mathcal{A})\hfill \\ & =(I-{\mathcal{M}}_2^{-1}\mathcal{A})(I-{\mathcal{M}}_1^{-1}\mathcal{A})\hfill \\ & =({\mathcal{Q}}^{*}\mathcal{Q}-{\mathcal{Q}}^{*}\mathcal{Q}{\mathcal{M}}_2^{-1}\mathcal{A}{\mathcal{Q}}^{*}\mathcal{Q})({\mathcal{Q}}^{*}\mathcal{Q}-{\mathcal{Q}}^{*}\mathcal{Q}{\mathcal{M}}_1\mathcal{A}{\mathcal{Q}}^{*}Q)\hfill \\ & =({\mathcal{Q}}^{*}\mathcal{Q}-{\mathcal{Q}}^{*}{\Lambda }_2\mathcal{Q})({\mathcal{Q}}^{*}\mathcal{Q}-{\mathcal{Q}}^{*}{\Lambda }_1\mathcal{Q})\hfill \\ & ={\mathcal{Q}}^{*}(I-{\Lambda }_2){\mathcal{Q}}^{*}\mathcal{Q}(I-{\Lambda }_1)\mathcal{Q}\hfill \\ & ={\mathcal{Q}}^{*}(I-{\Lambda }_2)(I-{\Lambda }_1)\mathcal{Q}\hfill \end{array}$$

by definition we have

$$\begin{array}{cc}\hfill \rho \left(\mathcal{T}\right)& \le \underset{1\le i,j\le nm}{max}\left|(1-{\lambda }_i^{\left(2\right)})(1-{\lambda }_j^{\left(1\right)})\right|\hfill \\ & \le \underset{1\le i\le nm}{max}|(1-{\lambda }_i^{\left(2\right)})|\underset{1\le j\le nm}{max}\left|(1-{\lambda }_j^{\left(1\right)})\right|\hfill \\ & =\rho (I-{\mathcal{M}}_2^{-1}\mathcal{A})\rho (I-{\mathcal{M}}_1^{-1}\mathcal{A})\hfill \\ & =\rho \left({\mathcal{M}}_2^{-1}{\mathcal{N}}_2\right)\rho \left({\mathcal{M}}_1^{-1}{\mathcal{N}}_1\right)\hfill \\ & \le \left|\right|{\mathcal{M}}_2^{-1}{\mathcal{N}}_2{\left|\right|}_{{\mathcal{M}}_2}\left|\right|{\mathcal{M}}_1^{-1}{\mathcal{N}}_1{\left|\right|}_{{\mathcal{M}}_1}\hfill \end{array}$$

Therefore, by Lemma 5 we have

$$\rho \left(\mathcal{T}\right)\le {\theta }_1^3\frac{max|\lambda \left(N_1\right)|+max|\lambda \left(Q_1\right)|}{{\lambda }_{min}\left(M_1\right)+{\lambda }_{min}\left(P_1\right)}{\theta }_2^3\frac{max|\lambda \left(N_2\right)|+max|\lambda \left(Q_2\right)|}{{\lambda }_{min}\left(M_2\right)+{\lambda }_{min}\left(P_2\right)}={\varrho }_1{\varrho }_2$$

and this completes the proof. □

3. Numerical Results

All numerical experiments presented in this section were computed in double precision with a number of MATLAB R2018a codes. All iterations are started from the initial zero matrix $X^{\left(0\right)}$ and terminated when the current iteration satisfies

$$\frac{\parallel R^{\left(k\right)}{\parallel }_F}{\parallel R^{\left(0\right)}{\parallel }_F}\le {10}^{-8},$$

where $R^{\left(k\right)}=C-A{X}^{\left(k\right)}-X^{\left(k\right)}B$ is the residual of the kth iterate. Also, we use the tolerance $\epsilon =0.01$ for inner iterations in corresponding methods. For each experiment, we report the number of iterations or the number of total outer iteration steps (shown as out-itr) and CPU time. In the tables, the norm of the residual is shown as res-norm.

The MSI method was compared with two Hermitian- and skew-Hermitian-based splitting methods such the NSCG [17] and the HSS [6], and two familiar iterative methods such the GMRES [34] and the BiCGSTAB [14] methods. Note that although the NS-CGNR method is a Hermitian- and skew-Hermitian-based splitting method, it works well for problems with a dominant skew-Hermitian part and not efficiently for problems with a dominant Hermitian part, see [18]. Therefore, it is not fair to compare this version of the MSI method with it.

Example 1.

For this example, we use the matrices

$$A=B=M+2rN+\frac{100}{{(n+1)}^2}I,$$

where $M,N\in {\mathbb{R}}^{n\times n}$ are the tridiagonal matrices given by

$$M=\mathrm{tridiag}(-1,2,-1)\mathrm{and}N=\mathrm{tridiag}(0.5,0,-0.5).$$

We consider $r=0.01$ and $n=m=256$ [6].

This class of problems may arise in the preconditioned Krylov subspace iteration methods used for solving the systems of linear equations resulting from the finite difference or Sinc–Galerkin discretization of various differential equations and boundary value problems [6].

We apply the iteration methods to this problem. The results are given in

Table 1. Results of the Example 1.

Method	Out-Itr	CPU Time	Res-Norm
MSI	7	10.54	2.0887 × 10−6
NSCG	7	8.19	2.1050 × 10−6
HSS	298	75.32	3.2107 × 10−6
GMRES(10)	151	40.56	3.0400 × 10−6
BiCGSTAB	255	17.91	2.4616 × 10−6

.

From the results presented in Table 1, it can be seen that for this problem, the MSI and the NSCG methods are more efficient than the other methods.

Example 2.

For the second experiment, consider $A=\mathrm{tridiag}(-2,4,-1)$ and $B=\mathrm{tridiag}(-1,4,-2)$ with dimensions $2048\times 2048$ and $128\times 128$, respectively.

This is a problem with a strong Hermitian part [42,43]. The numerical results for this problem are listed in

Table 2. Results of Example 2.

Method	Out-Itr	CPU Time	Res-Norm
MSI	7	5.66	4.3252 × 10−5
NSCG	9	6.18	7.7746 × 10−5
HSS	21	39.82	8.7083 × 10−5
GMRES(10)	3	2.84	2.3509 × 10−6
BiCGSTAB	14	2.67	9.7977 × 10−5

.

Regarding Table 2, it is obvious that although the MSI method is more effective versus the NSCG and the HSS methods, the GMRES and the BiCGSTAB methods are more effective than it.

Example 3.

We consider the continuous Sylvester Equation (1) with $n=m$ and the coefficient matrices

$$\left\{\begin{array}{c}A=\mathrm{diag}(1,2,\cdots ,n)+r{L}^T,\hfill \\ B=2^{-t}{I}_n+\mathrm{diag}(1,2,\cdots ,n)+r{L}^T+2^{-t}L,\hfill \end{array}\right.$$

where L is the strictly lower triangular matrix having ones in the lower triangle part [6]. Here, t is a problem parameter to be specified in actual computations.

The iteration methods were used for this problem, and the results are given in

Table 3. Results of Example 3.

Method	Out-Itr	CPU Time	Res-Norm
MSI	5	16.43	0.0029
NSCG	8	21.65	0.0070
HSS	99	326.71	0.0288
GMRES(10)	20	49.87	0.0027
BiCGSTAB	75	16.22	0.0028

. Moreover, we compare the convergence history of the iterative methods by residual norm decreasing in Figure 1.

In Table 3, we report the number of outer iterations (out-itr), the CPU time and the residual norm (res-norm) after convergence. For this example, we observe that the MSI method is superior to the other iterative methods in terms of the number of iterations and it is similar to the BiCGSTAB method in terms of CPU time. Comparing the convergence history of the iterative methods by residual norm decreasing shows that the MSI method converges more rapidly and smoothly than the BiCGSTAB method (see Figure 1).

Example 4.

For this example, we used the nonsymmetric sparse matrix SHERMAN3 of dimension $5005\times 5005$ with 20033 nonzero entries from the Harwell–Boeing collection [44] instead of the coefficient matrix A. For the coefficient matrix B, we used $B=\mathrm{tridiag}(-1,4,-2)$ of dimension $8\times 8$ [17].

We apply the iteration methods to this problem, and the results are given in

Table 4. Results of Example 4.

Method	Out-Itr	CPU Time	Res-Norm
MSI	34	78.437	1.57 × 10−4
NSCG	64	121.265	2.61 × 10−4
HSS	>5000	>1000	2.32
GMRES(10)	>5000	>1000	247.77
BiCGSTAB	†	†	NaN

. Moreover, we compare the convergence history of the iterative methods by residual norm decreasing in Figure 2.

In Table 4, we report the number of outer iterations (out-itr), the CPU time and the residual norm (res-norm) after convergence or in 5000 outer iterations. For this example, we observe that the MSI method is superior to the other iterative methods in terms of the number of iterations and CPU times, the NSCG method has an acceptable performance. Furthermore, the HSS and the GMRES methods have a very slow convergence rate ( see Figure 2). From the Table 4 and Figure 2, we observe that the BiCGSTAB method was diverged for this problem. Therefore, we use splitting of the MSI, the NSCG and the HSS method as the splitting preconditioner denoted as MSI-BiCGSTAB, NSCG-BiCGSTAB and HSS-BiCGSTAB, respectively. The results of the preconditioned method for these preconditioners are given in

Table 5. Results of the preconditioned BiCGSTAB for Example 4.

Method	Out-Itr	CPU Time	Res-Norm
BiCGSTAB	†	†	NaN
MSI-BiCGSTAB	2	113.51	9.0749 × 10−6
NSCG-BiCGSTAB	4	263.92	5.6612 × 10−6
HSS-BiCGSTAB	†	†	NaN

. In Table 4 and Table 5 dagger (†) and notation “>1000” show that no solution has been obtained after 5000 iteration or CPU time is more than 1000 s respectively.

The results in Table 5 show that the use of the MSI method as a preconditioner improves the results obtained by the corresponding methods.

4. Future Work

As an important challenge in the fractional mathematical field, the extension of this work for solving the multi-term fractional Sylvester equation with frequency is of interest because most of the existing methods are only useful for untangling the multi-term fractional differential equation. Therefore, investigation of capturing appropriate methods that are applicable to solve a specific type of multi-term fractional equation is attractive and meaningful [45,46,47,48].

5. Conclusions

In this paper, we have proposed an efficient iterative method for solving the continuous Sylvester equation $AX+XB=C$. This method employs two symmetric positive definite splittings of the coefficient matrices A and B and present a multiplicative splitting iteration method. The convergence conditions have been derived based on the iteration matrix.

We have compared the MSI method with well-known iterative methods such as the NSCG method, the HSS method, the BiCGSTAB method and the GMRES method for some problems. We have observed that, for these problems, the MSI method is more efficient than the other methods.

In summary, by focusing on the results presented in Table 1, Table 3, Table 4 and Table 5, one can observe that the MSI method is often superior to the other iterative methods. Moreover, the use of the multiplicative splitting as a preconditioner can induce an accurate and effective preconditioned BiCGSTAB method.