# Numerical Solution of Sylvester Equation Based on Iterative Predictor-Corrector Method

Authors: Ovgu Cidar Iyikal
Journal: Journal of Mathematics
 (2022)

---

## Abstract

The inspiration of the study concerns an iterative predictor-corrector method with order of convergencep=45 for computing the inverse of the coefficient matrix Λ=In⊗A+BT⊗Im, which is obtained by the Sylvester equation AX+XB=C. The numerical solutions of three examples by predictor-corrector algorithm are given. The final numerical results also support the applicability, fast convergency, and high accuracy of the method for finding matrix inverses.


---

## Body

## 1. Introduction

Consider the Sylvester matrix equations of the formAX+XB=C where A∈ℝm×m, B∈ℝn×n, C∈ℝm×n are the constant coefficient matrices and X∈ℝm×n is the unknown solution matrix. These equations play important roles in matrix decompositions of the eigenvalues, in control theory, in reduction models, in numerical solutions of the matrix differential Riccati equations, and in image processing (see [1–6]).Solution of Sylvester matrix equations follows the transformation to the linear formΛcolX=colC, where Λ=In⊗A+BT⊗Im is an mn×mn matrix [7]. So, higher-order dimensions of Λ can create problems. These problems can be avoided with iterative algorithms for finding matrix inverses. An iterative optimization-based method on partial swarm theory to solve Sylvester equation was given in [8]. Iterative gradient-based algorithm for solving the equations by minimizing specific criterion functions was presented in [7]. A preconditioned iterative gradient-based method is obtained by selection of two auxiliary matrices by using the Newton method, and it can be acceptable as a generalization of the iterative splitting method to solve system of linear equations in [9]. Later, a preconditioned positive-definite skew-Hermitian iteration method of splitting for continuous Sylvester equations with positive-definite/semi-definite matrices was presented in [10]. The implicit iteration alternating direction for solving the continuous Sylvester equation, where the coefficient matrices are taking as positive semi-definite matrices, and at least one of them to be positive definite is given in [11]. The implicit iteration alternating direction method for solving the continuous Sylvester equation, when the coefficient matrices are taking as positive semidefinite matrices (at least one of them to be positive definite) is given in [12]. A method based on the sign function iteration for solving large-scale Sylvester equation was presented in [13]. A transformation method such as Hessenberg–Schur algorithm was also used to solve Sylvester equation by reducing the coefficient matrices A and B to triangular form in [14]. A numerical Arnoldi-based method for solving Sylvester equation when A is sparse and large for partial pole-assignment problem for large matrices was given in [15]. A method consisting of orthogonal reduction of the coefficient matrix A of the Sylvester matrix equation to a block-upper-Hessenberg form was given in [16].Other than these methods, direct methods generally take time and require a lot of storage. Special interest is to use an iterative method for computing the Moore–Penrose inverseΛ^†∈ℂn2×n1 of the matrix Λ^∈ℂn1×n2. This method is derived by the Schulz method of 2nd order. A pth order iterative method for computing the Moore–Penrose inverse was studied in [17, 18] for the orders p≥2, and it is known as the hyperpower iterative method. Let Λ^∈ℂn1×n2, and we surveyed on iterative methods by which Vk presents the approximate Moore–Penrose inverse of Λ^ at the kth iteration. It is usually presented as(1)Vk+1=Vk∑j=0p−1Tkj,k=0,1,…,for n1≤n2. Tk=I−Λ^Vk performs p matrix-by-matrix multiplications per step. There are different representations of (1) which differ in the calculation of the power sum ∑j=0p−1Tkj. If the nested loops are used through the factorizations, then the computational effort in the pth order method decreases and the number of matrix-by-matrix multiplications and additions in the polynomial is reduced (see [19, 20]). Most recently, in [21], for a given integer t≥1, two different classes of iterative methods for finding matrix inverses of square nonsingular matrices were proposed and they were used as approximate inverse preconditioners for solving linear systems. Among these classes of methods, class 1 converges with order p=3∗2t+1 and class 2 converges with order p=5∗2t−1, requiring 2t+4 and 3t+4 matrix-by-matrix multiplications per step, respectively. The method of orders p=7,11,15,19 to approximate matrix inversion was applied for approximating the Schur complement matrices. The given algorithm was applied to solve Fredholm integral equations of first kind in [22]. Another recursive approach for constructing incomplete block-matrix factorization of the M-matrices by iterative two-step method for the approximation of the inverse of pivoting the diagonal block-matrices at each stage of the recursion was given in [23].In this study, an iterative predictor-corrector method of orderp=45 is used for computing the inverse of nonzero coefficient matrix Λ∈ℝmn×mn of the linear form ΛcolX=colC of the Sylvester equations. It was established in [24]. The method requires 5 matrix-by-matrix multiplications per step for the predictor step and 5 matrix-by-matrix multiplications per step for the corrector step. Then, to verify the predictor-corrector method’s fast convergence, experimental analysis is conducted for the three examples, whereas for large dimensional matrices, the computation of the Moore–Penrose inverse is costly when the direct method as singular value decomposition is used. Therefore, fast converging iterative methods for approximating the Moore–Penrose inverses that are also computationally efficient are essential. The need to compute Moore-Penrose inverse for the solution of the Sylvester equation leads to one of the application area. As stated before, these equations appear in control theory and other important real-life phenomena.


## 2. Sylvester Matrix Equation

Consider the Sylvester matrix equations of the form(2)AX+XB=C,where A∈ℝm×m, B∈ℝn×n, and C∈ℝm×n are constant matrices and X∈ℝm×n is the matrix of unknowns. Linear system (2) has a unique solution given as(3)colX=Λ−1colC,if and only if λiA+λjB≠0 for any i and j [7], where ⊗ presents the Kronecker product for(4)Λ=In⊗A+BT⊗Im∈ℝmn×mn,(5)colX=x1Tx2T…xmnTT,(6)colC=c1Tc2T…cmnTT.In practice, ifBT=A, then equation (2) becomes continuous time Lyapunov equation. For numerical solution of(7)ΛcolX=colC,the following iterative predictor-corrector method is used to find inverse of Λ in (4) for solving (7).


## 3. The Iterative Predictor-Corrector Method for Approximating Moore–Penrose Inverse of the Matrix

LetΛ^∈ℂn1×n2, I denote the unit matrix, Λ^∗ denote the conjugate transpose of Λ^, RΛ^ denote the range of Λ^, and PRΛ^ denote the orthogonal projection of RΛ^. The following matrix-valued functions are defined:(8)ΩTk=Tk+Tk2,(9)ΨTk=I+Tk2,(10)ΓTk=Tk4,where Tk=I−Λ^Vk. We propose the predictor-corrector method in [24] as follows.Predictor-corrector method:(11)In:V0=αΛ^∗,G:Tk=I−Λ^Vk,ΦTk=ΨTkΩTk,P:Vk+1/2=VkI+ΦTk,Tk+1/2=I−Λ^Vk+1/2,C:Vk+1=Vk+1/2I+ΦTk+1/2I+ΓTk+1/2,k=0,1,2,…,where In represents the initial step, G represents the generator, P represents the predictor at the fractional step k+1/2, and C represents the corrector for approximating the Moore–Penrose inverse Vk+1 of Λ^ at the k+1th step in (11). The predictor-corrector method given in (11) is effective in computational aspects and is useful for n1≤n2. If n1>n2, then dual version of the predictor-corrector iterative algorithm (11) can be used.Theorem 1 (see [24]).
LetΛ^∈ℂn1×n2and let the nonzero eigenvalues ofΛ^Λ^∗ satisfy λ1Λ^∗Λ^≥λ2Λ^∗Λ^≥…≥λrΛ^∗Λ^>0 where the real scalar α satisfies 0<α<2/λ1Λ^∗Λ^. Then, the sequence Vk+1 obtained by the predictor-corrector method (11) converges to the Moore–Penrose inverse Λ^† in the L2 norm when k⟶∞ with order of convergence p=45 and the asymptotic convergence factor is [21] ACF=10/ln45≈2.627. Moreover, the following error estimate is valid:(12)Λ^†−Vk+12≤αR0245k+1Λ^∗21−R02,whose degree of practicality increases in the given order and degree of precision decreases in the same order and R0=PRΛ^−αΛ^Λ^∗.


### 3.1. Error Analysis

The minimum norm of the solution to the linear system:(13)ΛcolX=colC,Λ∈ℝmn×mn,colC∈ℝm×n,where colC∈RΛ with RΛ denotes the range of Λ, colX∈ℝm×n denotes the unknown solution, and rankΛ=r. The general least squares solution of (13) or minimum norm least squares solution is the solution of the problem(14)minx∈ScolX2,S=colX∈ℝm×n∣ΛcolX−colC2=min.System (14) has a unique solution obtained by the equation colX=Λ†colC, and it is called the pseudoinverse solution [25].IfΛ is nonsingular matrix, then Λ†=Λ−1. Moreover, the condition number of Λ is defined by condΛ=Λ2Λ†2.


### 3.2. Algorithm for Approximating Inverse of MatrixΛ to Solve Sylvester Equation

LetcolXk+1=Vk+1colC be the approximate solution of (14), where Vk+1 is the inverse of Λ obtained by the predictor-corrector method given in (11) at the k+1−th iteration by using initial approximation V0=αΛT, and α satisfies 0<α<2/λ1ΛTΛ where λ1ΛTΛ≥λ2ΛTΛ≥…≥λrΛTΛ>0. Therefore, the residual error at the corresponding step is τk+1=colC−ΛcolXk+1. Given a predescribed accuracy ε>0, we propose the following algorithm based on Algorithm A in [24] that uses the predictor-corrector method (11) and approximates the pseudoinverse for Λ (See Algorithm 1).Algorithm 1: Predictor-corrector algorithm
Choose an initial matrixV0=αΛT under the assumption 0<α<2/λ1ΛTΛ where λ1ΛTΛ≥λ2ΛTΛ≥…≥λrΛTΛ>0.Step 1. Let k=0. Evaluate colX0=V0colC and τ0=colC−ΛcolX0; then, calculate τ02/colC2.Do Step 2–Step 6 untilτk2/colC2≤ε.Step 2. Evaluate Tk=I−ΛVk.Step 3. Apply P, the predictor step of the predictor-corrector method (11), to find Vk+1/2.Step 4. Evaluate Tk+1/2=I−ΛVk+1/2.Step 5. Apply C, the corrector step of the predictor-corrector method (11), to find Vk+1.Step 6. Evaluate colXk+1=Vk+1colC and τk+1=colC−ΛcolXk+1, respectively. Then, calculate τk+12.Letk=k+1.Step 7. If k is the performed iteration number, then the approximate solution is yk=colXk that satisfies τk2/colC2≤ε.


## 3.1. Error Analysis

The minimum norm of the solution to the linear system:(13)ΛcolX=colC,Λ∈ℝmn×mn,colC∈ℝm×n,where colC∈RΛ with RΛ denotes the range of Λ, colX∈ℝm×n denotes the unknown solution, and rankΛ=r. The general least squares solution of (13) or minimum norm least squares solution is the solution of the problem(14)minx∈ScolX2,S=colX∈ℝm×n∣ΛcolX−colC2=min.System (14) has a unique solution obtained by the equation colX=Λ†colC, and it is called the pseudoinverse solution [25].IfΛ is nonsingular matrix, then Λ†=Λ−1. Moreover, the condition number of Λ is defined by condΛ=Λ2Λ†2.


## 3.2. Algorithm for Approximating Inverse of MatrixΛ to Solve Sylvester Equation

LetcolXk+1=Vk+1colC be the approximate solution of (14), where Vk+1 is the inverse of Λ obtained by the predictor-corrector method given in (11) at the k+1−th iteration by using initial approximation V0=αΛT, and α satisfies 0<α<2/λ1ΛTΛ where λ1ΛTΛ≥λ2ΛTΛ≥…≥λrΛTΛ>0. Therefore, the residual error at the corresponding step is τk+1=colC−ΛcolXk+1. Given a predescribed accuracy ε>0, we propose the following algorithm based on Algorithm A in [24] that uses the predictor-corrector method (11) and approximates the pseudoinverse for Λ (See Algorithm 1).Algorithm 1: Predictor-corrector algorithm
Choose an initial matrixV0=αΛT under the assumption 0<α<2/λ1ΛTΛ where λ1ΛTΛ≥λ2ΛTΛ≥…≥λrΛTΛ>0.Step 1. Let k=0. Evaluate colX0=V0colC and τ0=colC−ΛcolX0; then, calculate τ02/colC2.Do Step 2–Step 6 untilτk2/colC2≤ε.Step 2. Evaluate Tk=I−ΛVk.Step 3. Apply P, the predictor step of the predictor-corrector method (11), to find Vk+1/2.Step 4. Evaluate Tk+1/2=I−ΛVk+1/2.Step 5. Apply C, the corrector step of the predictor-corrector method (11), to find Vk+1.Step 6. Evaluate colXk+1=Vk+1colC and τk+1=colC−ΛcolXk+1, respectively. Then, calculate τk+12.Letk=k+1.Step 7. If k is the performed iteration number, then the approximate solution is yk=colXk that satisfies τk2/colC2≤ε.


## 4. Numerical Examples

In this section, all the calculations are carried out by Mathematica program in double precision for Examples1–3. The tables adopt the following notation:(15)TCP:total CPU time in seconds.Example 1.
Consider the following coefficient matrices for Sylvester equation (2) (see [7]):(16)A=112−4,B=11−11,C=310−12−8,with exact solution(17)X=1235.
Here,X0=ones2,2×10−6 (ones2,2 represent the 2×2 matrix of ones) and apply the predictor-corrector method to find the inverse of Λ in (4).
It is obtained as(18)Λ=21−102−30−11021012−3,and then, compute approximate solution of the Sylvester equation (2). In this example, the right side vector C is calculated by AX+XB. We apply the proposed predictor-corrector algorithm to obtain the inverse of Λ, by performing k iterations for an accuracy of ε≤1×10−6 for the corresponding Sylvester matrix equation (7). Table 1 presents L2 norm of the errors using the predictor-corrector algorithm and the condition numbers.
The best error of the approximate solution for Example1 obtained in [7] was 2.5974×10−4 after k=15 iterations. When we compare results in Table 1 with the literature [7], the predictor-corrector method is a more accurate, fast, and effective method.Table 1
L2 norm of the errors obtained by predictor-corrector algorithm and condition numbers for Example 1.
τ12/colC24.893274×10−5colX−V1colC2/colX25.814136×10−5colX−V1colC23.630926×10−4condA5.000000condΛ3.365853Example 2.
Sylvester equation (2) with the following coefficient matrices is considered:(19)A=triurandp,p,1+diagα+diagrandp,B=triurandt,t,1+diagα+diagrandp,where triu represents a tridiagonal matrix, rand represents random numbers, diag represents the diagonal matrix, and α is used for the weights of the diagonal entries [9]. For Example 2, the exact solution is(20)X=randp,t+speyep,t×2,with X0=onesp,t×10−6, where speye represents a sparse identity matrix and α=2. Stopping criterion is ε=1×10−9. Table 2 presents TCP, total iteration number k, and L2 norm errors for given values of the dimensions p and t.
It was already expected that the error would increase with the increase in the condition numbers as shown in Table2. Figure 1 illustrates the error err=τk2/colC2 for the test problems in Example 2 with respect to iteration number k for α=0.1,0.5,1,5 when p,t=10,10. The obtained numerical results are compared with the results given in [7] and those given in [9]. Hence, our results by predictor-corrector algorithm are more accurate and more faster than methods given in [7, 9].
According to Figure1, fast and accurate result can be obtained when α increases.Table 2
TCP, iteration numberk, and the relative errors obtained by predictor-corrector algorithm for Example 2.
p,tTCPkτk2/colC2condΛ5,30.0124.57222×10−172.3680310,60.1723.05786×10−163.8111015,102.2524.74074×10−164.0343020,1517.7922.54129×10−1012.992550,302167.1538.37652×10−1640.3054Figure 1
Error versus iteration numberk for Example 2 when p,t=10,10.Example 3.
Consider the following randomly generated coefficient matricesA, B, C∈ℝ6×6 [16] to solve Sylvester equation (2) via (7) by predictor-corrector algorithm with stopping criteria ε=1×10−15:(21)a11=0.846221417824,a12=0.202647357650,a13=0.891298966149,a21=0.525152496305,a22=0.672137468474,a23=0.762096833027,a31=0.921812970745,a32=0.176266144495,a33=0.935469699108,a41=0.738207245811,a42=0.405706213062,a43=0.916904439913,a51=0.138890881957,a52=0.198721742661,a53=0.272187924970,a61=0.20276521856,a62=0.603792479194,a63=0.198814267761,a14=0.456467665168,a15=0.821407164295,a16=0.615432348100,a24=0.018503643248,a25=0.444703364353,a26=0.791937037427,a34=0.410270206991,a35=0.057891304784,a36=0.813166497304,a44=0.893649530914,a45=0.352868132217,a46=0.009861300661,a54=0.015273927029,a55=0.445096432288,a56=0.465994341675,a64=0.746785676564,a65=0.931814578462,a66=0.418649467728.b11=0.546571151829,b12=0.522590349081,b13=0.875741899818,b21=0.444880204673,b22=0.880142207411,b23=0.737305988466,b31=0.694567240426,b32=0.172956141275,b33=0.136518742260,b41=0.621310130795,b42=0.979746896789,b43=0.011756687353,b51=0.794821080201,b52=0.271447258642,b53=0.893897966445,b61=0.956843448445,b62=0.252329346874,b63=0.199138067206,b14=0.298723012102,b15=0.582791681561,b16=0.579806873250,b24=0.661442576382,b25=0.423496256851,b26=0.760365009804,b34=0.284408589750,b35=0.515511752141,b36=0.529823116717,b44=0.469224285211,b45=0.333951479972,b46=0.640526498990,b54=0.064781122963,b55=0.432906596107,b56=0.209069404438,b64=0.988334938278,b65=0.225949868144,b66=0.379818370351.c11=0.838118445052,c12=0.709471392703,c13=−0.52465110045,c21=0.019639513865,c22=0.428892365341,c23=−0.157705064793,c31=0.681277161282,c32=0.304617366869,c33=−0.259097521287,c41=0.379481018028,c42=0.189653747547,c43=0.045079530894,c51=0.831796017610,c52=0.193431156405,c53=−0.195440134055,c61=0.502812883996,c62=0.682223223591,c63=0.128597436672,c14=0.883859755846,c15=0.920621353311,c16=0.337873490541,c24=0.906034170972,c25=0.486903569691,c26=0.307318127939,c34=0.356617663467,c35=0.936266830386,c36=0.135805928663,c44=0.742419230344,c45=0.836601295625,c46=0.670950124615,c54=0.592369870448,c55=0.589358440912,c56=0.57915246589,c64=0.975132462127,c65=0.918056131072,c66=0.352654884696.
Figure2 illustrates the error err=τk2/colC2 for the test problems in Example 3, and Table 3 presents L2 norm of the errors using the predictor-corrector algorithm and the condition numbers. The accuracy of the predictor-corrector algorithm is remarkable.Figure 2
Error versus iteration numberk for Example 3.Table 3
L2 norm of the errors obtained by predictor-corrector algorithm for Example 3.
τ52/colC2τ52condAcondΛ4.23589×10−151.50175×10−1470.801051.51


## 5. Concluding Remarks

An iterative predictor-corrector method of convergence orderp=45 is used for computing the inverse of a nonzero matrix Λ∈ℝmn×mn in (4) to solve Sylvester equation (2). It converges fast, and it is a highly accurate method. It can also be useful when Λ is rectangular matrix or an ill-conditioned matrix. The predictor-corrector iterative method may also be applied to precondition the coefficient matrices of the linear algebraic system of equations obtained by using finite difference method to solve boundary value problems, for example, for the solution of Laplace’s equation with singularities (see [26, 27]), for the solution of the heat equation on hexagonal grid (see [28]), and for the approximation of the derivatives of the solution of the heat equation (see [29, 30]). Moreover, the iterative method of predictor-corrector for finding matrix inverses may be applied to precondition the coefficient matrices of the system of equations obtained by using finite element method to solve parabolic partial differential equations (see [31]).



---

*Source: 6571126-2022-08-22.xml*