1. bookVolume 6 (2021): Issue 2 (July 2021)
Journal Details
License
Format
Journal
First Published
01 Jan 2016
Publication timeframe
2 times per year
Languages
English
access type Open Access

The Incomplete Global GMERR Algorithm for Solving Sylvester Equation

Published Online: 03 Feb 2021
Page range: 1 - 6
Received: 27 Mar 2020
Accepted: 02 Dec 2020
Journal Details
License
Format
Journal
First Published
01 Jan 2016
Publication timeframe
2 times per year
Languages
English
Abstract

In this paper, the Incomplete Global GMERR algorithm and the Global GMERR algorithm are used to solve the Sylvester equation. The numerical experiment is given to compare the CPU run time and the number of iterations of the two methods.

Keywords

MSC 2010

Introduction

In control and communication theory, Sylvester matrix equation AX+XB=C AX + XB = C has important applications, where A ∈ ℝn×n, B ∈ ℝs×s, C,X ∈ ℝn×s. Many scholars have studied and explored its solution [1]. For (1), if B = AT, then the Lyapunov matrix equation is obtained. When n, s is very small, some direct algorithms such as the Hessenberg-Schur algorithm [2] can be used. When the matrix order is large, we can use Krylov subspace method to solve [3], [4], [5].

We have proposed the Incomplete Global GMERR algorithm [6], the Global GMERR algorithm [7]. In this paper, we use them to solve (1).

Let λ (A) and λ (−B) be the set of eigenvalues of A and −B, respectively. When and only when λ (A) ∩ λ (−B) = φ, (1) has a unique solution. Then solve (1) as follows.

The equation (1) is equivalent to A˜X˜=C˜ \tilde A\tilde X = \tilde C where à is stated by the following equation A˜=IsA+BTIs \tilde A = {I_s} \otimes A + {B^T} \otimes {I_s} where ⊗ represents the Kronecker product MN = [mi,jN], and X˜=(x1,1,...,xn,1,x1,2,...,xn,2,...,x1,s,...,xn,s)TnsC˜=(c1,1,...,cn,1,c1,2,...,cn,2,...,c1,s,...,cn,s)Tns \matrix{ \tilde{X}=(x_{1,1},...,x_{n,1},x_{1,2},...,x_{n,2},...,x_{1,s},...,x_{n,s})^T\in\, \mathbb{R}^{n s} \cr \tilde{C}=(c_{1,1},...,c_{n,1},c_{1,2},...,c_{n,2},...,c_{1,s},...,c_{n,s})^T\in\, \mathbb{R}^{n s}}

Thus, we can transform (2) with a direct algorithm. In this case, the amount of computation and the amount of storage increase as the number of iterations increases. Therefore, we hope that the problem (1) can be solved directly, rather than through the above transformation. On the basis of the block Krylov subspace method, some methods [4] are proposed to solve large Sylvester equation (sn). Each iteration must solve a low rank Sylvester equation. The Global FOM and GMRES algorithm [3] were proposed by Jbilou et al. Salkuyeh D K and Toutounian F [4] proposed GFSL method and GGSL method to solve some large Sylvester equations.

In this paper, we give new algorithms: the Global GMERRSL algorithm and Incomplete Global GMERRSL algorithm.

Applied to Sylvester Equation

We give a corresponding rule: XA(X)=AX+XB X\rightarrow\mathbb{A}(X)=AX+XB where 𝔸(ℝn×s → ℝn×s) is an operatr. then 𝔸T(ℝn×s → ℝn×s) has the following form: XAT(X)=ATX+XBT X\rightarrow\mathbb{A}^T(X)=A^TX+XB^T and then (1) can be showed as follows A(X)=C \mathbb{A}(X)=C

On the basis of the Global Arnoldi method [6], we give a new method to approximate the solution of (4). Give the initial estimate X0, then R0 = C − 𝔸(X). The algorithm solves (4), and the approximation Xm obtained in step m is satisfied XmX0=ZmATKm(AT,R0) X_m-X_0=Z_m\in\mathbb{A}^TK_m(\mathbb{A}^T,R_0) and Rm=CA(Xm)Km(AT,R0) R_m=C-\mathbb{A}(X_m)\perp K_m(\mathbb{A}^T,R_0) where 𝔸kR0 = 𝔸(𝔸(k−1)R0)

For the above definition, the following lemma [4] tells us that the Krylov subspace Km(𝔸,R0) is the same as the Krylov subspace Km(A,R0).

Lemma 1

Let 𝔸 be defined as described above, assuming that R0 is full rank, then Km(𝔸,R0) = Km(A,R0), where m ≥ 1.

It is easy to get from Lemma 1 Km(AT,R0)=Km(AT,R0)ATKm(AT,R0)=ATKm(AT,R0) \matrix{ K_m(\mathbb{A}^T,R_0)=K_m(A^T,R_0) \cr \mathbb{A}^TK_m(\mathbb{A}^T,R_0)=A^TK_m(A^T,R_0)} then, (5) and (6) are equivalent to XmX0=ZmATKm(AT,R0) X_m-X_0=Z_m\in\mathbb{A}^TK_m(A^T,R_0) Rm=CA(Xm)Km(AT,R0). R_m=C-\mathbb{A}(X_m)\perp K_m(A^T,R_0)

By (7), we have Xm=X0+ATUm*ym, {X_m} = {X_0} + {A^T}{U_m}*{y_m}, where Um = [V1,V2,...,Vm] is a set of orthogonal bases of Km(𝔸T, R0) produced by the Global Arnoldi algorithm, and * is a vector product defined as Um*α=i=1mαiVi, {U_m}*\alpha = \sum\limits_{i = 1}^m {\alpha _i}{V_i}\;, where α = (α1,α2,...,αm) ∈ ℝm.

We define inner product < X,Y >= tr(XTY), where X,Y ∈ ℝn×s, and tr(XTY) is the trace of XTY. Then, based on (8), we have <Vi,R0A(ATUm*ym)>=0,i=1,2,...,m, <V_i,R_0-\mathbb{A}(A^TU_m*y_m)>=0,\ i=1,2,...,m, then <Vi,R0AATUm*ymATUm*ymB>=0,i=1,2,...,m,<Vi,AATUm*ym+ATUm*ymB>=<Vi,R0>,i=1,2,...,m. \matrix{ { < {V_i},{R_0} - A{A^T}{U_m}*{y_m} - {A^T}{U_m}*{y_m}B > = 0,\;i = 1,2,...,m,} \cr { < {V_i},A{A^T}{U_m}*{y_m} + {A^T}{U_m}*{y_m}B > = < {V_i},{R_0} > ,\;i = 1,2,...,m.} \cr }

Therefore, we have j=1m<Vi,AATVj+ATVjB>ym(j)=<Vi,R0>,i=1,2,...,m. \sum\limits_{j = 1}^m < {V_i},A{A^T}{V_j} + {A^T}{V_j}B > y_m^{(j)} = < {V_i},{R_0} > ,\;i = 1,2,...,m. where ym=(ym(1),ym(2),...,ym(m))T {y_m} = (y_m^{(1)},y_m^{(2)},...,y_m^{(m)}{)^T} .

Then, let V1 = R0/ ||R0 ||F, when i = 1, (tr(V1TAATV1+V1TATV1B),tr(V1TAATV2+V1TATV2B),...,tr(V1TAATVm+V1TATVmB))*ym=R0F, (tr(V_1^TA{A^T}{V_1} + V_1^T{A^T}{V_1}B),tr(V_1^TA{A^T}{V_2} + V_1^T{A^T}{V_2}B),...,tr(V_1^TA{A^T}{V_m} + V_1^T{A^T}{V_m}B))*{y_m} = \parallel {R_0}{\parallel _F}, when i = 2,3,...,m, (tr(ViTAATV1+ViTATV1B),tr(ViTAATV2+ViTATV2B),...,tr(ViTAATVm+ViTATVmB))*ym=0. (tr(V_i^TA{A^T}{V_1} + V_i^T{A^T}{V_1}B),tr(V_i^TA{A^T}{V_2} + V_i^T{A^T}{V_2}B),...,tr(V_i^TA{A^T}{V_m} + V_i^T{A^T}{V_m}B))*{y_m} = 0. Then, we get ym by the following matrix equation: (tr(V1TAATV1+V1TATV1B)tr(V1TAATV2+V1TATV2B)...tr(V1TAATVm+V1TATVmB)tr(V2TAATV1+V2TATV1B)tr(V2TAATV2+V2TATV2B)...tr(V2TAATVm+V2TATVmB)tr(VmTAATV1+VmTATV1B)tr(VmTAATV2+VmTATV2B)...tr(VmTAATVm+VmTATVmB))ym=(R0F00) \matrix{ {\left( {\matrix{ {tr(V_1^TA{A^T}{V_1} + V_1^T{A^T}{V_1}B)} & {tr(V_1^TA{A^T}{V_2} + V_1^T{A^T}{V_2}B)} & {...} & {tr(V_1^TA{A^T}{V_m} + V_1^T{A^T}{V_m}B)} \cr {tr(V_2^TA{A^T}{V_1} + V_2^T{A^T}{V_1}B)} & {tr(V_2^TA{A^T}{V_2} + V_2^T{A^T}{V_2}B)} & {...} & {tr(V_2^TA{A^T}{V_m} + V_2^T{A^T}{V_m}B)} \cr \vdots & \vdots & \ddots & \vdots \cr {tr(V_m^TA{A^T}{V_1} + V_m^T{A^T}{V_1}B)} & {tr(V_m^TA{A^T}{V_2} + V_m^T{A^T}{V_2}B)} & {...} & {tr(V_m^TA{A^T}{V_m} + V_m^T{A^T}{V_m}B)}\cr } } \right){y_m} = \left( {\matrix{ {\parallel {R_0}{\parallel _F}} \cr 0 \cr \vdots \cr 0 \cr } } \right)} \hfill \cr }

Therefore, the restarting GMERRSL Algorithm is given as following:

Algorithm 1 : The Global GMERRSL

Step 1 Set restart step number m, then 2 ≤ qm, Give the tolerance tol, initial guess X0(n × s). Then R0 = CAX0X0B. Set V1 = R0/ ||R0 ||F;

Step 2 Find a set of orthogonal basis of Km(AT, V1) : V1, V2,...,Vm produced by the Global Arnoldi algorithm;

Step 3 Get ym by solving the matrix equation equation (9);

Step 4 Compute Xm = X0 + ATUm * ym;

Step 5 If ||R0 ||F=||BAX0 ||Ftol, stop;

if not, let X0 = Xm, go to step 1.

Similarly, by the article [7], we can solve the Sylvester equations with the Incomplete Global GMERR algorithm as follows:

Algorithm 2 The Incomplete Global GMERRSL

Step 1 Set the restart step number m, then 2 ≤ qm, Give the tolerance tol, initial guess X0(n × s). Then R0 = CAX0X0B. Set V1 = R0/ ||R0 ||F;

Step 2 Find a set of orthogonal basis of Km(AT, V1) : V1,V2,...,Vm produced by the Incomplete Global Arnoldi algorithm;

Step 3 Get ym by solving the following equation (tr(V1TAATV1+V1TATV1B)tr(V1TAATV2+V1TATV2B)...tr(V1TAATVm+V1TATVmB)tr(V2TAATV1+V2TATV1B)tr(V2TAATV2+V2TATV2B)...tr(V2TAATVm+V2TATVmB)tr(VmTAATV1+VmTATV1B)tr(VmTAATV2+VmTATV2B)...tr(VmTAATVm+VmTATVmB))ym=(R0F00tr(Vq+2TR0)tr(VmTR0)) \matrix{ {\left( {\matrix{ {tr(V_1^TA{A^T}{V_1} + V_1^T{A^T}{V_1}B)} & {tr(V_1^TA{A^T}{V_2} + V_1^T{A^T}{V_2}B)} & {...} & {tr(V_1^TA{A^T}{V_m} + V_1^T{A^T}{V_m}B)} \cr {tr(V_2^TA{A^T}{V_1} + V_2^T{A^T}{V_1}B)} & {tr(V_2^TA{A^T}{V_2} + V_2^T{A^T}{V_2}B)} & {...} & {tr(V_2^TA{A^T}{V_m} + V_2^T{A^T}{V_m}B)} \cr \vdots & \vdots & \ddots & \vdots \cr {tr(V_m^TA{A^T}{V_1} + V_m^T{A^T}{V_1}B)} & {tr(V_m^TA{A^T}{V_2} + V_m^T{A^T}{V_2}B)} & {...} & {tr(V_m^TA{A^T}{V_m} + V_m^T{A^T}{V_m}B)} \cr } } \right){y_m} = \left( {\matrix{ {\parallel {R_0}{\parallel _F}} \cr 0 \cr \vdots \cr 0 \cr {tr(V_{q + 2}^T{R_0})} \cr \vdots \cr {tr(V_m^T{R_0})} \cr } } \right)} \hfill \cr }

Step 4 Compute Xm = X0 + ATUm * ym;

Step 5 If ||R0 ||F=||BAX0 ||Ftol, stop;

if not, let X0 = Xm, then, go to step 1.

Numerical Examples

The example is computed with MATLAB, we take A as a n = 1000 order Toeplitz matrix, (310.5310.5310.5310.5313) \left( {\matrix{ 3 & 1 & {0.5} & {} & {} & {} & {} \cr {} & 3 & 1 & {0.5} & {} & {} & {} \cr {} & {} & 3 & 1 & {0.5} & {} & {} \cr {} & {} & {} & \ddots & \ddots & \ddots & {} \cr {} & {} & {} & {} & 3 & 1 & {0.5} \cr {} & {} & {} & {} & {} & 3 & 1 \cr {} & {} & {} & {} & {} & {} & 3 \cr } } \right)

B is the s = 10 order and s = 100 order Toeplitz matrix, respectively. And the right-hand matrix C = rand(n, s) is the random matrix. The Global GMERR (25) (GL-GMERR) and the Incomplete Global GMERR (25)(IGL-GMERR) are used to solve the two Sylvester equations. The restart step number is 25. Let Initial matrix X0 = 0.

Table 1 gives the results of the Sylvester equation (s = 10 and s = 100) by Gl-GMERR (25) and the IGl-GMERR (25). CPU represents the algorithm run time (seconds), l is the number of iterations, and the ratio is the run time ratio of IGL-GMERR(25) to the GL-GMERR (25) used at the same required tolerance. Obviously, when q < 25, the IGL-GMERR(25) algorithm is faster than the GL-GMERR(25), and the running time is shorter. When q = 25, the two are the same.

The Global GMERRSL (25) and the Incomplete Global GMERRSL (25) (s = 10 and s = 100)

s q cpu ratio l || R ||F
10 2 20.1400 0.7133 13 8.0989e-7
10 5 21.2500 0.7526 13 8.0989e-7
10 10 24.5940 0.8711 13 8.0989e-7
10 20 26.6250 0.9430 13 8.0978e-7
10 25 28.2340 1 14 4.1594e-7
100 2 374.0780 0.7039 16 7.6921e-7
100 5 405.8430 0.7637 16 7.6921e-7
100 10 446.2820 0.7039 16 7.6921e-7
100 20 492.4680 0.7039 16 7.6955e-7
100 25 531.4220 1 17 3.7744e-7
Conclusion

The Incomplete Global GMERR and the Global GMERR are used to solve the Sylvester equation respectively in this paper. The result: the Incomplete Global GMERR algorithm is more efficient. How to apply this method to other matrix equations and its convergence analysis need further research.

The Global GMERRSL (25) and the Incomplete Global GMERRSL (25) (s = 10 and s = 100)

s q cpu ratio l || R ||F
10 2 20.1400 0.7133 13 8.0989e-7
10 5 21.2500 0.7526 13 8.0989e-7
10 10 24.5940 0.8711 13 8.0989e-7
10 20 26.6250 0.9430 13 8.0978e-7
10 25 28.2340 1 14 4.1594e-7
100 2 374.0780 0.7039 16 7.6921e-7
100 5 405.8430 0.7637 16 7.6921e-7
100 10 446.2820 0.7039 16 7.6921e-7
100 20 492.4680 0.7039 16 7.6955e-7
100 25 531.4220 1 17 3.7744e-7

J. Z. Hearon, (1992), Nonsingular solutions of TABT = C. Linear Algebra Appl. 174:283–314. HearonJ. Z. 1992 Nonsingular solutions of TABT = C Linear Algebra Appl 174 283 314 Search in Google Scholar

R. H. Bartel and G W. Stewart, (1994), Algorithm 432:Solution of the matrix equation AX + XB = C. Circ Syst Signal Proc, 13:820–826. BartelR. H. StewartG W. 1994 Algorithm 432:Solution of the matrix equation AX + XB = C Circ Syst Signal Proc 13 820 826 Search in Google Scholar

K. Jbilou, A. Messaoudi and H. Sadok, (1999), Global Fom and GMRES algorithms for matrix equations. Appl Numer Math, 31:49–63. JbilouK. MessaoudiA. SadokH. 1999 Global Fom and GMRES algorithms for matrix equations Appl Numer Math 31 49 63 Search in Google Scholar

D. K. Salkuyeh and F. Toutounian, (2006), New approaches for solving large Sylvester equations. Appl Math Comput, 73:9–18. SalkuyehD. K. ToutounianF. 2006 New approaches for solving large Sylvester equations Appl Math Comput 73 9 18 Search in Google Scholar

Y. Q. Lin, (2004), Implicitly restarted global FOM and GMRES for nonsymmetric matrix equations and Sylvester equations. Appl Math Comput, 167:1004–1025. LinY. Q. 2004 Implicitly restarted global FOM and GMRES for nonsymmetric matrix equations and Sylvester equations Appl Math Comput 167 1004 1025 Search in Google Scholar

Y. H. Zheng, J. L. Li, D. X. Cheng and L. L. Lv, (2018), The Incomplete Global GMERR Algorithm to Solve AX = B. Journal of Computational Analysis and Applications, 24:760–772. ZhengY. H. LiJ. L. ChengD. X. LvL. L. 2018 The Incomplete Global GMERR Algorithm to Solve AX = B Journal of Computational Analysis and Applications 24 760 772 Search in Google Scholar

Y. H. Zheng, D. X. Cheng and X. H. Qian, (2012), Global GMERR algorithm for linear systems with multiple right-hand sides. Journal of Northeast Normal University (Natural Science Edition), 44:41–45. ZhengY. H. ChengD. X. QianX. H. 2012 Global GMERR algorithm for linear systems with multiple right-hand sides Journal of Northeast Normal University (Natural Science Edition) 44 41 45 Search in Google Scholar

Recommended articles from Trend MD

Plan your remote conference with Sciendo