This work is licensed under the Creative Commons Attribution 4.0 International License.
Introduction
In control and communication theory, Sylvester matrix equation
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
\tilde A\tilde X = \tilde C
where à is stated by the following equation
\tilde A = {I_s} \otimes A + {B^T} \otimes {I_s}
where ⊗ represents the Kronecker product M ⊗ N = [mi,jN], and
\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 (s ≪ n). 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:
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:
X\rightarrow\mathbb{A}^T(X)=A^TX+XB^T
and then (1) can be showed as follows
\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
X_m-X_0=Z_m\in\mathbb{A}^TK_m(\mathbb{A}^T,R_0)
and
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 R0is full rank, then Km(𝔸,R0) = Km(A,R0), where m ≥ 1.
It is easy to get from Lemma 1\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
X_m-X_0=Z_m\in\mathbb{A}^TK_m(A^T,R_0)R_m=C-\mathbb{A}(X_m)\perp K_m(A^T,R_0)
By (7), we have
{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
{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
<V_i,R_0-\mathbb{A}(A^TU_m*y_m)>=0,\ i=1,2,...,m,
then
\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
\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
{y_m} = (y_m^{(1)},y_m^{(2)},...,y_m^{(m)}{)^T}
.
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.
