A weighted global GMRES algorithm with deflation for solving large Sylvester matrix equations

TL;DR

This paper introduces a novel weighted global GMRES algorithm with deflation for efficiently solving large Sylvester matrix equations, combining residual-based weighting and deflated restarting strategies.

Contribution

It proposes a new weighted global GMRES method with deflation for large Sylvester equations, improving efficiency without changing inner products, and provides theoretical and numerical validation.

Findings

The new algorithm effectively solves large Sylvester equations.

It reduces computational costs compared to existing methods.

Numerical examples demonstrate improved performance and accuracy.

Abstract

The solution of large scale Sylvester matrix equation plays an important role in control and large scientific computations. A popular approach is to use the global GMRES algorithm. In this work, we first consider the global GMRES algorithm with weighting strategy, and propose some new schemes based on residual to update the weighting matrix. Due to the growth of memory requirements and computational cost, it is necessary to restart the algorithm efficiently. The deflation strategy is popular for the solution of large linear systems and large eigenvalue problems, to the best of our knowledge, little work is done on applying deflation to the global GMRES algorithm for large Sylvester matrix equations. We then consider how to combine the weighting strategy with deflated restarting, and propose a weighted global GMRES algorithm with deflation for solving large Sylvester matrix equations. Theoretical analysis is given to show why the new algorithm works effectively. Further, unlike the weighted GMRES-DR presented in [{\sc M. Embree, R. B. Morgan and H. V. Nguyen}, {\em Weighted inner products for GMRES and GMRES-DR}, (2017), arXiv:1607.00255v2], we show that in our new algorithm, there is no need to change the inner product with respect to diagonal matrix to that with non-diagonal matrix, and our scheme is much cheaper. Numerical examples illustrate the numerical behavior of the proposed algorithms.

Tables

Table 1. Table 1: Example 1: Performance of W-GLGMRES with different choices of D 𝐷 D

			$n=22500$			$n=40000$
$s$	$m$	$D$	iter	res.norm	CPU	iter	res.norm	CPU
		$I$	271	9.6980e-07	1.7603e+03	502	9.8712e-07	3.6526e+03
		$D_1$	225	8.9007e-07	1.4532e+03	361	9.0881e-07	2.6656e+03
25	10	$D_2$	196	8.9191e-07	1.2564e+03	375	9.9170e-07	2.7972e+03
		$D_3$	149	8.9245e-07	964.4450	247	9.2564e-07	1.8723e+03
		$I$	287	9.6370e-07	1.0593e+03	511	9.9187e-07	2.6007e+03
		$D_1$	224	8.6859e-07	819.1613	350	8.0319e-07	1.7281e+03
16	10	$D_2$	221	9.5079e-07	828.4745	278	9.0354e-07	1.2168e+03
		$D_3$	147	8.9043e-07	540.1904	253	9.1782e-07	1.1863e+03
		$I$	135	9.9273e-07	1.0023e+03	229	9.4550e-07	2.5741e+03
		$D_1$	93	9.2973e-07	690.4734	164	9.6980e-07	1.8403e+03
16	15	$D_2$	85	8.8796e-07	612.2571	138	9.1576e-07	1.5791e+03
		$D_3$	77	8.9375e-07	558.9338	125	8.8805e-07	1.4288e+03

Table 2. Table 2: Summary of of the test matrices for Example 2–Example 4.

Matrix	n	nnz	Density	Density	Application area
saylr4	3564	22316	0.0017	real unsymmetric	Oil reservoir modeling
add32	4960	19848	0.0008	real unsymmetric	Electronic circuit design
sherman4	1104	3786	0.0031	real unsymmetric	Oil reservior modeling
sherman2	1080	23094	0.0198	real unsymmetric	Oil reservoir modeling
pde2961	2961	14585	0.0231	real unsymmetric	Partial differential equation

Table 3. Table 3: Example 2: Numerical results of W-GLGMRES with different weighting stratiges.

		$m=10$			$m=20$
Problem	Weighting Strategy	iter	res.norm	CPU	iter	res.norm	CPU
$A=sherman4$	[14]	102	1.8165e-06	997.0492	39	1.7936e-06	1.0142e+03
$B=fdm(cos\left(xy\right),e^{y^{2}x},100)$	[31]	143	7.9507e-07	1.4222e+03	44	9.6308e-07	1.8846e+03
$n=1104,s=400$	$D_3$	38	6.1240e-07	371.8751	16	8.9240e-07	374.3088
$A=add32$	[14]	25	8.0562e-06	1.6571e+03	12	3.4655e-06	1.3972e+03
$B=fdm(sin\left(xy\right),e^{xy},10)$	[31]	40	8.5137e-07	2.2993e+03	24	8.0102e-07	2.7319e+03
$n=4960,s=400$	$D_3$	26	8.1566e-07	957.1628	10	6.9810e-07	934.2591
$A=saylr4$	[14]	84	1.4897e-06	3.8107e+03	34	1.2236e-06	2.9038e+03
$B=fdm(sin\left(xy\right),e^{xy},10)$	[31]	115	9.2482e-07	4.9978e+03	40	6.45505e-07	4.7578e+03
$n=3564,s=400$	$D_3$	42	9.0481e-07	1.7269e+03	16	7.1080e-07	1.8709e+03
$A=sherman2$	[14]	31	7.7570e-07	996.7216	13	7.4740e-07	843.9966
$B=fdm(sin\left(xy\right),e^{xy},10)$	[31]	48	8.0974e-07	1.8602e+03	16	6.9319e-07	2.5471e+03
$n=1080,s=400$	$D_3$	20	1.9301e-07	544.5839	12	2.0916e-07	668.4175
$A=pde2961$	[14]	32	1.3552e-06	1.0066e+03	12	1.4062e-06	879.7433
$B=fdm(sin\left(xy\right),e^{xy},10)$	[31]	47	8.5788e-07	1.5451e+03	16	8.6097e-07	1.1724e+03
$n=2961,s=400$	$D_3$	28	6.6487e-07	878.394	11	9.8081e-07	795.3090

Table 4. Table 4: Example 3: Numerical results of the four algorithms.

	Problem	$(m,k)$	Algorithm	iter	res.norm	CPU
		(20,-)	GLGMRES	60	9.2014e-07	871.0590
		(20,-)	W-GLGMRES	18	7.8295e-07	387.4749
	$A=sherman4$	(20,10)	GLGMRES-D	7	5.8933e-07	146.6579
1		(20,10)	W-GLGMRES-D	4	6.1324e-13	79.3758
	$n=1104,s=400$	(20,15)	GLGMRES-D	7	5.5915e-08	103.9108
		(20,15)	W-GLGMRES-D	5	1.7006e-14	68.0248
		(30,20)	GLGMRES-D	8	3.0370e-08	243.1432
		(30,20)	W-GLGMRES-D	5	5.1292e-14	179.8224
		(20,-)	GLGMRES	14	6.7001e-07	1.1432e+03
		(20,-)	W-GLGMRES	11	2.9613e-07	1.2481e+03
	$A=add32$	(20,10)	GLGMRES-D	7	2.0322e-08	581.1866
2		(20,10)	W-GLGMRES-D	4	2.7931e-09	348.0351
	$n=4960,s=400$	(20,15)	GLGMRES-D	7	3.5901e-08	549.6019
		(20,15)	W-GLGMRES-D	4	1.2383e-09	309.2734
		(30,20)	GLGMRES-D	6	4.6590e-08	991.0881
		(30,20)	W-GLGMRES-D	4	3.4711e-07	700.5107
		(20,-)	GLGMRES	41	7.9900e-07	1.9829e+03
		(20,-)	W-GLGMRES	18	6.1862e-07	1.5938e+03
	$A=saylr4$	(20,10)	GLGMRES-D	7	2.1766e-07	445.2339
3		(20,10)	W-GLGMRES	4	9.2134e-07	257.8798
	$n=3564,s=400$	(20,15)	GLGMRES-D	7	6.9319e-10	424.5378
		(20,15)	W-GLGMRES-D	6	2.3769e-09	346.7902
		(30,20)	GLGMRES-D	7	5.1598e-08	687.8864
		(30,20)	W-GLGMRES-D	5	9.8619e-07	584.5981
		(20,-)	GLGMRES	$\uparrow 2500$	$†$	-
		(20,-)	W-GLGMRES	$\uparrow 2500$	$†$	-
	$A=sherman2$	(20,10)	GLGMRES-D	5	1.8630e-07	101.9935
4		(20,10)	W-GLGMRES-D	4	5.2809e-08	96.6050
	$n=1080,s=400$	(20,15)	GLGMRES-D	6	8.8742e-09	92.2590
		(20,15)	W-GLGMRES-D	5	5.9891e-11	88.5462
		(30,20)	GLGMRES-D	6	3.1034e-10	211.2737
		(30,20)	W-GLGMRES-D	5	8.0189e-12	210.1942
		(20,-)	GLGMRES	16	7.8910e-07	865.9652
		(20,-)	W-GLGMRES	11	5.1991e-07	606.6930
	$A=pde2961$	(20,10)	GLGMRES-D	7	9.6839e-09	411.2545
5		(20,10)	W-GLGMRES-D	4	5.0377e-08	238.1900
	$n=2961,s=400$	(20,15)	GLGMRES-D	6	3.7837e-07	348.5320
		(20,15)	W-GLGMRES-D	5	3.7809e-08	289.2056
		(30,20)	GLGMRES-D	6	1.0925e-08	576.2694
		(30,20)	W-GLGMRES-D	4	8.1073e-10	391.9157

Table 5. Table 5: Example 4: A comparison of flexible GLGMRES and flexible W-GLGMRES.

		$m=20$
Problem	Algorithm	iter	res.norm	CPU
	GLGMRES	60	9.2014e-07	844.0590
$A=sherman4$	W-GLGMRES	18	7.8295e-07	398.4749
$B=fdm(\cos \left(xy\right),e^{y^{2}x},100)$
$n=1104,s=400$	F-GLGMRES	10	3.1852e-07	798.0550
	WF-GLGMRES	6	9.4672e-07	521.1000
	GLGMRES	14	6.7001e-07	1.1432e+03
$A=add32$	W-GLGMRES	11	2.9613e-07	1.0521e+03
$B=fdm(\sin \left(xy\right),e^{xy},10)$
$n=4960,s=400$	F-GLGMRES	3	7.9173e-07	812.2611
	WF-GLGMRES	2	6.8148e-11	520.8404
	GLGMRES	41	7.9900e-07	1.9829e+03
$A=saylr4$	W-GLGMRES	18	6.1862e-07	1.3438e+03
$B=fdm(\sin \left(xy\right),e^{xy},10)$
$n=3564,s=400$	F-GLGMRES	5	4.9185e-08	1.2897e+03
	WF-GLGMRES	4	6.0365e-09	1.0238e+03
	GLGMRES	$\uparrow 2500$	$†$	-
$A=sherman2$	W-GLGMRES	$\uparrow 2500$	$†$	-
$B=fdm(\sin \left(xy\right),e^{xy},10)$
$n=1080,s=400$	F-GLGMRES	$\uparrow 2500$	$†$	-
	WF-GLGMRES	$\uparrow 2500$	$†$	-
	GLGMRES	16	7.8910e-07	755.9652
$A=pde2961$	W-GLGMRES	11	5.1991e-07	661.4130
$B=fdm(\sin \left(xy\right),e^{xy},10)$
$n=2961,s=400$	F-GLGMRES	4	7.9803e-08	777.3842
	WF-GLGMRES	3	6.1608e-08	573.7325

Table 6. Table 6: Example 4: A comparison of flexible W-GLGMRES-D and flexible GLGMRES-D.

		$m=20,k=15$
Problem	Algorithm	iter	res.norm	CPU
$A=sherman4$	F-GLGMRES-D	7	8.39783e-08	513.5553
$B=fdm(\cos \left(xy\right),e^{y^{2}x},100)$
$n=1104,s=400$	WF-GLGMRES-D	5	8.0266e-07	370.8139
$A=add32$	F-GLGMRES-D	4	1.0656e-08	1.2597e+03
$B=fdm(\sin \left(xy\right),e^{xy},10)$
$n=4960,s=400$	WF-GLGMRES-D	3	8.7016e-12	928.9309
$A=saylr4$	F-GLGMRES-D	5	5.9416e-08	1.0732e+03
$B=fdm(\sin \left(xy\right),e^{xy},10)$
$n=3564,s=400$	WF-GLGMRES-D	3	2.1611e-07	640.1040
$A=sherman2$	F-GLGMRES-D	$\uparrow 2500$	$†$	-
$B=fdm(\sin \left(xy\right),e^{xy},10)$
$n=1080,s=400$	WF-GLGMRES-D	$\uparrow 2500$	$†$	-
$A=pde2961$	F-GLGMRES-D	5	2.5480e-08	887.0873
$B=fdm(\sin \left(xy\right),e^{xy},10)$
$n=2961,s=400$	WF-GLGMRES-D	3	4.3535e-07	530.0198

Equations

$$AX+XB=C,$$

$$\mathcal{A}$$

$$X\to$$

$$\mathcal{A}X=C.$$

$$\big{(}I_{s}\otimes A+B^{T}\otimes I_{n}\big{)}vec(X)=vec(C),$$

$$vec(X)=[X(:,1);X(:,2);\dots ;X(:,s)],$$

$$(\mathbf{u},\mathbf{v}{)}_D={\mathbf{v}}^{T}D\mathbf{u}$$

$$\Vert \mathbf{u}{\Vert }_D$$

$$(Y,Z{)}_D=$$

$$(Y,Z)_{D}=\big{(}vec(Y),vec(Z)\big{)}_{I_{s}\otimes D}.$$

$$\Vert Y{\Vert }_D$$

$$(A^T{\diamond }_{D}B{)}_{ij}=(A_i,B_j{)}_D,i=1,2,\dots ,p,j=1,2,\dots ,l.$$

$$(V_i,V_j{)}_D$$

$${\mathcal{K}}_m(\mathcal{A},V_1)$$

$$[\mathcal{A}{V}_1,\mathcal{A}{V}_2,\dots ,\mathcal{A}{V}_m]$$

$$=$$

$${\bar{H}}_m$$

$$\mathcal{A}{\mathcal{V}}_m\equiv [\mathcal{A}{V}_1,\mathcal{A}{V}_2,\dots ,\mathcal{A}{V}_m],$$

$$\mathcal{A}{\mathcal{V}}_m={\mathcal{V}}_{m+1}({\bar{H}}_m\otimes I_s).$$

$$X_m=X_0+{\mathcal{V}}_m({\mathbf{y}}_m^w\otimes I_s),$$

$$R_m$$

$$R_m=C-\mathcal{A}{X}_m{\perp }_D\mathcal{A}{\mathcal{K}}_m(\mathcal{A},R_0),$$

$$\underset{\mathbf{y}}{\min }\Vert R_m{\Vert }_D$$

$${\mathbf{y}}_m^w={\mathrm{argmin}}_{{\mathbf{y}}^w\in {\mathbb{R}}^m}\Vert \beta {\mathbf{e}}_1-{\bar{H}}_m{\mathbf{y}}^w{\Vert }_2,$$

$${\bar{H}}_m^T{\bar{H}}_m{\mathbf{y}}_m^w={\bar{H}}_m^T\beta {\mathbf{e}}_1.$$

$${\hat{D}}_j=I_s\otimes D_j,j=1,2,3,$$

$$\mathcal{A}{\mathcal{V}}_m{Y}_i-{\theta }_i{\mathcal{V}}_m{Y}_i{\perp }_D(\mathcal{A}-\sigma I){\mathcal{K}}_m(\mathcal{A},V),i=1,\dots ,k,$$

$$\mathcal{A}{\mathcal{V}}_m=[\mathcal{A}{V}_1,\mathcal{A}{V}_2,\dots ,\mathcal{A}{V}_m]=[A{V}_1+V_{1}B,\dots ,A{V}_m+V_{m}B],$$

$$\mathcal{A}{\mathcal{V}}_m({\mathbf{g}}_i\otimes I_s)-{\theta }_i{\mathcal{V}}_m({\mathbf{g}}_i\otimes I_s)$$

$$\theta_{i}\big{(}(\mathcal{A}\mathcal{V}_{m})^{T}\diamond_{D}\mathcal{V}_{m}\big{)}{\bf g}_{i}=(\mathcal{A}\mathcal{V}_{m})^{T}\diamond_{D}(\mathcal{A}\mathcal{V}_{m}){\bf g}_{i}.$$

$${\theta }_i{H}_m^T{\mathbf{g}}_i={\bar{H}}_m^T{\bar{H}}_m{\mathbf{g}}_i.$$

Full text

1 Introduction

Consider the large Sylvester matrix equation

[TABLE]

where $A\in\mathbb{R}^{n\times n},B\in\mathbb{R}^{s\times s}$, $C\in\mathbb{R}^{n\times s}$ and $X\in\mathbb{R}^{n\times s}$, with $s\ll n$. If we define the operator $\mathcal{A}$ as

[TABLE]

Then the Sylvester matrix equation (1) can be written as

[TABLE]

The Sylvester matrix equation (1) plays an important role in control and communications theory, model reduction, image restoration, signal processing, filtering, decoupling techniques for ordinary partial differential equations, as well as the implementation of implicit numerical methods for ordinary differential equations; see [2, 3, 4, 5, 11, 13, 33], and the references there in. Also, (3) can be rewritten as the following large linear system

[TABLE]

where $\otimes$ denotes the Kronecker product operator, and $vec(X)$ denotes the vectorize operator defined as (in MATLAB notation)

[TABLE]

and $X(:,i)$ is the $i$-th column of the matrix $X\in\mathbb{R}^{n\times s}$. The linear equation systems (4) have unique solution if and only if the matrix $(I_{s}\otimes A)+(B^{T}\otimes I_{n})$ is nonsingular. Throughout this paper, we assume that the system (4) has a unique solution. However, the size of the linear equation systems (4) would be very huge. Therefore, we apply some iterative algorithms for solving (1) instead of (4).

There are some iterative algorithms based on the block or the matrix Krylov solvers for the solution of the Sylvester matrix equations, see, e.g. [1, 6, 11, 13, 15, 17, 19, 21, 27, 32, 33]. The main idea behind these algorithms is to exploit the global or block (extended) Arnoldi process to construct $F$-orthonormal or orthonormal bases for the matrix or block Krylov subspaces, respectively, and then apply some projection techniques to extract approximations.

In [9], Essai introduced a weighted Arnoldi process for solving large linear systems. The idea is to improve and accelerate the convergence rate of the standard algorithm by constructing a $D$-orthonormal basis for the Krylov subspace, where $D$ is called the weighting matrix and is generally a positive diagonal matrix. According to [8, 12], weighting strategy can improve the algorithm by alienating the eigenvalues that obstacle the convergence. The weighting strategy has been successfully developed for solving linear systems [14, 31], matrix equations [26, 28], and large eigenvalue problems [36]. For example, Mohseni Moghadam et al. [26] presented a weighted global FOM method for solving nonsymmetric linear systems. They used the Schur complement formula and a new matrix product, and gave some theoretical results to show rationality of the proposed algorithm. In [28], Panjeh Ali Beik et al. proposed weighted global FOM and weighted global GMRES algorithms for solving the general coupled linear matrix equations.

For the sake of the growth of memory requirements and computational cost, the global Krylov subspace algorithms will become impractical as the step of the global Arnoldi process proceeds. For Krylov subspace algorithms, one remedy is to use some restarting strategies [30]. A popular restarting strategy is the deflated restarting (also refer to as thick-restarting or deflation) strategy advocated in [18, 23, 24, 25, 34, 35, 36], in which the approximate eigenvectors are put firstly in the search subspace. Here “deflated restarting” (or deflation) means computing some approximate eigenvectors corresponding to some eigenvalues, and using them to “deflate” these eigenvalues from the spectrum of the matrix, to speed up the convergence of the iterative algorithm. The deflation strategy is popular for the solution of large linear systems [23, 24] and large eigenvalue problems [18, 25, 34, 35, 36], to the best of our knowledge, little work is done on applying the deflated restarting strategy on the global GMRES algorithm for large Sylvester matrix equations.

In this paper, we try to fill in this gap. As was pointed out in [8, 9, 12], the optimal choice of the weighting matrix $D$ in the weighted approaches is still an open problem and needs further investigation. We first apply the weighting strategy to the global GMRES algorithm, and present three new schemes to update the weighting matrix at each restart. To accelerate the convergence of the weighted global GMRES algorithm, we consider how to knit the deflation strategy together with it, and the key is that the Sylvester matrix equation can be rewritten as a linear system of the form (3) theoretically. The new algorithm can be understood as applying the deflation technique to remove some small eigenvalues of the matrix $(I_{s}\otimes A+B^{T}\otimes I_{n})$ at each restart. Theoretical results and numerical experiments show that the weighting strategy with deflation can produce iterations that give faster convergence than the conventional global GMRES algorithm, and a combination of these two strategies is more efficient and robust than its two original counterparts.

This paper is organized as follows. After presenting the weighted global GMRES algorithm for the solution of Sylvester matrix equations in section 2, the deflated version of this algorithm is established in section 3. Some numerical experiments confirm the superiority of our new algorithm over the conventional ones in section 4.

2 A weighted global GMRES algorithm for large Sylvester matrix equations

In this section, we recall some notations and definitions that will be used in this paper, and briefly introduce the weighted global Arnoldi process as well as the weighted global GMRES algorithm. Specifically, we propose three new schemes based on residual to update the weighting matrix during iterations.

The global generalized minimal residual (GLGMRES) algorithm is well-known for solving linear systems with multiple right-hand sides and for matrix equations [16, 27], which is an oblique projection method based on matrix Krylov subspace. Let us introduce the weighted global GMRES algorithm for Sylvester matrix equations. Let $D=diag(d_{1},d_{2},\ldots,d_{n})$ be a diagonal matrix with $d_{i}>0,~{}i=1,2,\ldots,n$, and let ${\bf u},{\bf v}\in\mathbb{R}^{n}$ be given, then the $D$-inner product with respect to two vectors is defined as [36]

[TABLE]

and the associated $D$-norm of $\bf u$ is defined as

[TABLE]

For two matrices $Y,Z\in\mathbb{R}^{n\times s}$, the $D$-inner product is defined as [26]

[TABLE]

where $trace(\cdot)$ denotes the trace of a matrix, and $Z^{T}$ represents the transpose of the matrix $Z$. It can be verified that [26]

[TABLE]

Also, the $D$-norm associated with this inner product is

[TABLE]

Next we introduce a useful product that will be used latter:

Definition 1**.**

[26]** Let $A=[A_{1},A_{2},\ldots,A_{p}]$ and $B=[B_{1},B_{2},\ldots,B_{l}]$, where $A_{i},B_{j}\in\mathbb{R}^{n\times s},~{}i=1,2,\ldots,p,~{}j=1,2,\ldots,l$. Then elements of the matrix $A^{T}\diamond_{D}B$ is defined as

[TABLE]

It was shown that $A^{T}\diamond_{D}B=A^{T}\diamond(DB)$ [26]. Let $V\in\mathbb{R}^{n\times s}$ be an initial block vector that is $D$-orthogonal, that is, orthonormal with respect to the $\diamond_{D}$-inner product. The following algorithm presents an $m$-step weighted global Arnoldi process [26].

Algorithm 1**.**

The $m$-step weighted global Arnoldi process

Input:* $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{s\times s}$, $V\in\mathbb{R}^{n\times s}$, a positive diagonal matrix $D$ and an integer number $m>0$.* 2. 2.

Output:* $\mathcal{V}_{m}=[V_{1},V_{2},\ldots,V_{m}],\bar{H}_{m}=[h_{i,j}].$* 3. 3.

Compute $\beta=\Arrowvert V\Arrowvert_{D}$, $V_{1}=V/\beta.$ 4. 4.

for $j=1,2,\ldots,m$ 5. 5.

Compute $W=\mathcal{A}V_{j}=AV_{j}+V_{j}B.$ 6. 6.

for $i=1,2,\ldots,j$ 7. 7.

$h_{ij}=(W,V_{i})_{D}.$** 8. 8.

$W=W-h_{ij}V_{i}.$** 9. 9.

end 10. 10.

Compute $h_{j+1,j}=\Arrowvert W\Arrowvert_{D}$. If $h_{j+1,j}=0$ then stop. 11. 11.

Compute $V_{j+1}=W/h_{j+1,j}.$ 12. 12.

end

The weighted global Arnoldi process constructs a $D$-orthonormal basis $\mathcal{V}_{m}=[V_{1},V_{2},\ldots,V_{m}]\in\mathbb{R}^{n\times ms}$, i.e.,

[TABLE]

for the matrix Krylov subspace

[TABLE]

Let \bar{H}_{m}=\Big{[}\begin{array}[]{c}H_{m}\\ h_{m+1,m}e_{m}^{T}\end{array}\Big{]}\in\mathbb{R}^{(m+1)\times m} be a quasi-upper Hessenberg matrix whose nonzeros entries $h_{ij}$ are defined by Algorithm 1, and the matrix $H_{m}$ is obtained from the matrix $\bar{H}_{m}$ by deleting its last row. Note that the matrix $\mathcal{V}_{m}$ is $D$-orthogonal. With the help of Definition 1, we obtain the following relations

[TABLE]

where $\mathcal{V}_{m+1}=[V_{1},V_{2},\ldots,V_{m+1}]\in\mathbb{R}^{n\times(m+1)s}$. Define

[TABLE]

then (6) can be rewritten as

[TABLE]

We are in a position to consider the weighted global GMRES algorithm for solving (1). Let $X_{0}\in\mathbb{R}^{n\times s}$ be the initial guess, and the initial residual be $R_{0}=C-AX_{0}-X_{0}B$. In the weighted global GMRES algorithm, we construct an approximate solution of the form

[TABLE]

where ${\bf y}^{w}_{m}\in\mathbb{R}^{m}$. The corresponding residual is

[TABLE]

here we used (1) and (9). Substituting (7) into (12), we arrive at

[TABLE]

where ${\bf e}_{1}$ is the first canonical basis vector in $\mathbb{R}^{m+1}$. Note that the residual is $D$-orthogonal to $\mathcal{A}\mathcal{K}_{m}(\mathcal{A},R_{0})$, i.e.,

[TABLE]

where $\mathcal{A}\mathcal{K}_{m}(\mathcal{A},R_{0})=span\{\mathcal{A}R_{0},\ldots,\mathcal{A}^{m}R_{0}\}$, and “$\perp_{D}$” means orthgonal with respect to the “$\diamond_{D}$” inner product.

In order to compute ${\bf y}^{w}_{m}$, we have from (13) and (14) that

[TABLE]

where we used $\mathcal{V}_{m+1}^{T}\diamond_{D}\mathcal{V}_{m+1}=I_{m+1}$. Thus,

[TABLE]

or equivalently,

[TABLE]

As was pointed out in [9, 14, 36], the optimal choice of $D$ in the weighted approaches is still an open problem and needs further investigation. Some choices for the weighting matrix have been considered in, say, [14, 28, 31]. Also, to speed up the convergence rate, it was suggested to use a weighted inner product that changes at each restart [8, 12]. In this section, we propose three choices based on the residual $R_{m}$, which could be updated during iterations:

Option 1:

Let $\|R(:,t)\|_{2}=\max\{\|R_{m}(:,i)\|_{2},~{}1\leq i\leq s\}$, where $R_{m}(:,t)$ is the $t$-th column of residual matrix $R_{m}$. Then we define $D_{1}=diag\left(\frac{|R(:,t)|}{\|R(:,t)\|_{2}}\right)$, where $|R(:,t)|$ stands for the absolute value of $R(:,t)$.

Option 2:

Similarly, let $\|R(:,t)\|_{2}=\min\{\|R_{m}(:,i)\|_{2},~{}1\leq i\leq s\}$, then we define $D_{2}=diag\left(\frac{|R(:,t)|}{\|R(:,t)\|_{2}}\right)$.

Option 3:

Use the mean of the block residual $R_{m}$, i.e, $D_{3}=diag\left(|(\frac{\sum_{i=1}^{s}R_{m}(:,i)}{s})|\right)$.

Combining these weighting strategies with Algorithm 1, we have the following algorithm.

Algorithm 2**.**

A restarted weighted global GMRES algorithm for large Sylvester matrix equations (W-GLGMRES)

Input:* $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{s\times s}$, $C\in\mathbb{R}^{n\times s}$. Choose the initial guess, $X_{0}\in\mathbb{R}^{n\times s}$, an initial positive diagonal matrix D and an integer $m>0$, and the tolerance $tol>0$.* 2. 2.

Output:* The approximation $X_{m}$.* 3. 3.

Compute $R_{0}=C-AX_{0}-X_{0}B$ and $\beta=\Arrowvert R_{0}\Arrowvert_{D}$, $V=R_{0}/\beta.$ 4. 4.

Run Algorithm 1 with the initial block vector $V$ to obtain the matrices $\mathcal{V}_{m},\bar{H}_{m}$. 5. 5.

Solving $\min_{{\bf y}^{w}\in\mathbb{R}^{m}}\|\beta e_{1}-\bar{H}_{m}{\bf y}^{w}\|_{2}$ for ${\bf y}^{w}_{m}$. 6. 6.

Compute $X_{m}=X_{0}+\mathcal{V}_{m}({\bf y}^{w}_{m}\otimes I_{s})$ and $R_{m}=C-AX_{m}-X_{m}B$. If $\|R_{m}\|_{D}<tol$, then stop, else continue. 7. 7.

Set $X_{0}=X_{m}$ and update the weighting matrix $D$ according to Options 1–3, and go to Step 3.

Remark 2.1**.**

As was mentioned before, the Sylvester matrix equation (1) can be reformulated as the linear system (4). Thus, the three choices of $\{D_{j}\}_{j=1}^{3}$ for (1) can be understood as the weighted GMRES algorithm with the weighting matrices

[TABLE]

respectively, for solving the linear system (4). The theoretical results and discussions given in [8, 12] on weighted GMRES for large linear systems apply here trivially, and one refers to [8, 12] for why the weighted strategy can speed up the computation. This also interprets why the weighted strategy can improve the numerical performance of the standard global GMRES; see the numerical experiments made in Section 4.

3 A weighted global GMRES with deflation for large Sylvester matrix equations

In this section, we speed up the weighted global GMRES algorithm by using the deflated restarting strategy that is popular for large linear systems and large eigenvalue problems [18, 23, 24, 25, 35, 36]. In the first cycle of the weighted global GMRES algorithm with deflation, the standard weighted global GMRES algorithm is run. To apply the deflated restarting strategy, we need to compute $k~{}(1\leq k\leq m)$ weighted harmonic Ritz pairs. Let $\mathcal{V}_{m}$ be the $D$-orthonormal basis obtained from the “previous” cycle, we seek $k$ weighted harmonic Ritz pairs $(\theta_{i},Y_{i})$ that satisfy

[TABLE]

where

[TABLE]

and $Y_{i}={\bf g}_{i}\otimes I_{s}\in\mathbb{C}^{ms\times s}$ with ${\bf g}_{i}\in\mathbb{C}^{m\times 1}$. In this work, we want to deflate some smallest eigenvalues in magnitude, and a shift $\sigma=0$ is used throughout this paper.

From (10), we have that

[TABLE]

where \bar{I}_{m}=\left[\begin{array}[]{c}I_{m}\\ 0\end{array}\right]. By (17) and Definition 1, we can compute $(\theta_{i},{\bf g}_{i})$ via solving the following (small-sized) generalized eigenvalue problem

[TABLE]

From (10) and the fact that $\mathcal{V}_{m+1}^{T}\diamond_{D}\mathcal{V}_{m+1}=I_{m+1}$, we rewrite (18) as

[TABLE]

If $H_{m}$ is nonsingular, (19) is equivalent to

[TABLE]

Then we define the “weighted harmonic Ritz vector” as $Y_{i}={\bf g}_{i}\otimes I_{s}$, and the corresponding harmonic residual is

[TABLE]

where ${\bf\widetilde{r}}_{i}=(\bar{H}_{m}-\theta_{i}\bar{I}_{m}){\bf g}_{i},~{}i=1,\ldots,k$.

Remark 3.1**.**

In [7], a global harmonic Arnoldi method was proposed for computing harmonic Ritz pairs of large matrices. Here the difference is that our approach is based on the weighted projection techniques, and the method in [7] is a special case of ours as $D=I$.

In the following, we characterize the relationship between the weighted harmonic residuals and the residual from the weighted global GMRES algorithm. Let $X_{0}\in\mathbb{R}^{n\times s}$ be the initial guess and $R_{0}=C-AX_{0}-X_{0}B$ be the initial residual. After one cycle of the weighted global GMRES Algorithm, we have the approximate solution $X_{m}=X_{0}+\mathcal{V}_{m}({\bf y}^{w}_{m}\otimes I_{s})$, where ${\bf y}^{w}_{m}$ is defined in (15). The associated residual with respect to $X_{m}$ is

[TABLE]

where ${\bf{r}}_{m}\equiv\beta{\bf e}_{1}-\bar{H}_{m}{\bf y}^{w}_{m}$. The following result shows that ${\bf{r}}_{m}$ and $\{{\bf\widetilde{r}}_{i}\}_{i=1}^{k}$ are collinear with each other.

Theorem 2**.**

Let $\widetilde{R}_{i}=\mathcal{V}_{m+1}({\bf\widetilde{r}}_{i}\otimes I_{s}),~{}i=1,\ldots,k$, be the weighted harmonic residuals and $R_{m}=\mathcal{V}_{m+1}({\bf{r}}_{m}\otimes I_{s})$ be the weighted global GMRES residual, respectively, where ${\bf\widetilde{r}}_{i}=(\bar{H}_{m}-\theta_{i}\bar{I}_{m}){\bf g}_{i}$ and $\ {\bf{r}}_{m}=\beta{\bf e}_{1}-\bar{H}_{m}{\bf y}^{w}_{m}.$ Then $\{{\bf\widetilde{r}}_{i}\}_{i=1}^{k}$ and ${\bf{r}}_{m}$ are collinear with each other.

Proof.

Note that both the weighted harmonic residuals and the residual of the weighted global GMRES algorithm are in $range(\mathcal{V}_{m+1})$, and they are both $D$-orthogonal to $\mathcal{A}\mathcal{V}_{m}$. Thus, there is an $s\times s$ matrix $T$, such that

[TABLE]

so we have $({\bf\widetilde{r}}_{i}\otimes I_{s})=({\bf{r}}_{m}\otimes I_{s})T$. Let ${\bf{r}}_{m}=(\eta_{1},\eta_{2},\ldots,\eta_{m+1})^{T}$ and ${\bf\widetilde{r}}_{i}=(\tau_{1},\tau_{2},\ldots,\tau_{m+1})^{T}$, we obtain from (21) that

[TABLE]

which implies that $\{{\bf\widetilde{r}}_{i}\}_{i=1}^{k}$ and ${\bf{r}}_{m}$ are collinear with each other. ∎

We are ready to consider how to apply the deflated restarting strategy to the weighted global GMRES algorithm, and show that a weighted global Arnoldi-like relation still holds after restarting. Let $G_{k}=[{\bf g}_{1},{\bf g}_{2},\ldots,{\bf g}_{k}]$, and let $G_{k}=Q_{k}\Gamma_{k}$ be the reduced QR factorization. We stress that both forming $G_{k}$ and computing the QR decomposition can be implemented in real arithmetics; see Step 9 of Algorithm 3. Then we orthonormalize ${\bf{r}}_{m}$ against \left[\begin{array}[]{c}Q_{k}\\ 0_{1\times k}\end{array}\right] to get ${\bf q}_{k+1}$, and let Q_{k+1}=\left[\begin{array}[]{c|c}\begin{array}[]{c}Q_{k}\\ 0_{1\times k}\end{array}&{\bf q}_{k+1}\end{array}\right].

Notice that both (6) and (7) hold in the first cycle, i.e.,

[TABLE]

and thus

[TABLE]

That is,

[TABLE]

We note that

[TABLE]

and

[TABLE]

As a result,

[TABLE]

where we used $\{{\bf\widetilde{r}}_{i}\}_{i=1}^{k}$ and ${\bf{r}}_{m}$ are collinear with each other; see Theorem 2. Therefore,

[TABLE]

and

[TABLE]

Define $\mathcal{V}_{k}^{new}\equiv\mathcal{V}_{m}(Q_{k}\otimes I_{s})=[V_{1}^{new},\ldots,V_{k}^{new}]$, where $V_{i}^{new}\in\mathbb{R}^{n\times s},~{}1\leq i\leq k$, so we have

[TABLE]

If we denote

[TABLE]

then we have from (26) that

[TABLE]

and there is a $(k+1)s\times ks$ matrix $P$ such that

[TABLE]

The condition $\mathcal{V}_{k+1}^{new^{T}}\diamond_{D}\mathcal{V}_{k+1}^{new}=I_{k+1}$ yields

[TABLE]

where $\bar{H}_{k}^{new}\equiv Q_{k+1}^{T}\bar{H}_{m}Q_{k}\in\mathbb{R}^{(k+1)\times k}$ is generally a dense matrix. In conclusion, we obtain

[TABLE]

Then, we run the weighted global Arnoldi Algorithm from index $(k+1)$ to $m$ with the last $s$ columns of $\mathcal{V}_{k+1}^{new}$ as the starting matrix to build a new $m$-step global Arnoldi relation. However, $D$ is updated in our new algorithm, see Step 11 of Algorithm 3. In other words, we will perform the remaining $m-k$ steps of global Arnoldi process with a new $D$ (denoted by $D^{new}$ that is from the “current” residual $R_{m}$) after deflated restarting.

Denote $\mathcal{V}_{m}=[\mathcal{V}_{k+1}^{new}~{}|~{}\mathcal{V}_{\perp}^{new}]$, where $\mathcal{V}_{\perp}^{new}=[V_{k+2}^{new},\ldots,V_{m}^{new}]$, then $\mathcal{V}_{\perp}^{new}$ is $D^{new}$-orthogonal and it is also $D^{new}$ orthogonal to $\mathcal{V}_{k+1}^{new}$. Denote by $D^{old}$ the weighting matrix obtained from the residual of the “previous” cycle, we define $\widetilde{D}$-orthogonality of $\mathcal{V}_{m}$ as follows

[TABLE]

That is, $\mathcal{V}_{m}^{T}\diamond_{\widetilde{D}}\mathcal{V}_{m}=I_{m}$, and a global Arnoldi-like relation still holds after restarting.

Remark 3.2**.**

In the weighted GMRES-DR presented in [8, pp.20], Embree et al. showed how to restart the weighted GMRES-DR algorithm with a change of inner product by using the Cholesky factorization. However, the new weighting matrix is non-diagonal any more in their strategy, and the computational cost will become much higher when computing the $D$-inner products with respect to non-diagonal matrices. Thanks to (31), we indicate that without changing inner products in the weighted and deflated restarting algorithms, a global Arnoldi-like relation is still hold. So there is no need to change the inner product with respect to diagonal matrix to that with non-diagonal matrix, and our scheme is cheaper.

In summary, we have the following theorem.

Theorem 3**.**

A global Arnoldi-like relation holds for the weighted global Arnoldi algorithm with deflation

[TABLE]

We are ready to present the main algorithm of this paper.

Algorithm 3**.**

A weighted global GMRES with deflation for large Sylvester matrix equations (W-GLGMRES-D)

Input:* $A\in\mathbb{R}^{n\times n},B\in\mathbb{R}^{s\times s},C\in\mathbb{R}^{n\times s}$. Choose an initial guess $X_{0}\in\mathbb{R}^{n\times s},$ a positive diagonal matrix $D$, the positive integer number $m$ and a convergence tolerance $tol>0$.* 2. 2.

Output:* The approximation $X_{m}$.* 3. 3.

Compute $R_{0}=C-AX_{0}-X_{0}B$, $\beta=\Arrowvert R_{0}\Arrowvert_{D}$ and $V_{1}=R_{0}/\beta.$ 4. 4.

Set ${\bf c}=[\beta;0_{m\times 1}].$ 5. 5.

Run the weighted global Arnoldi algorithm to obtain $\mathcal{V}_{m+1}$ and $\bar{H}_{m}$. 6. 6.

Solve $\min_{{\bf y}^{w}}\|{\bf c}-\bar{H}_{m}{\bf y}^{w}\|_{2}$ for ${\bf y}^{w}_{m}$. 7. 7.

Compute $X_{m}=X_{0}+\mathcal{V}_{m}({\bf y}^{w}_{m}\otimes I_{s})$, and $R_{m}=C-AX_{m}-X_{m}B$. 8. 8.

If $\|R_{m}\|_{D}<tol$, then stop, else continue. 9. 9.

*Compute the $k$ eigenpairs $(\theta_{i},{\bf g}_{i})$ with the smallest magnitude from (19) or (20). Set $G_{k}=[{\bf g}_{1},\ldots,{\bf g}_{k}]$: We first sperate the *$\{{\bf g}_{i}\}^{\prime}$s into real and imaginary parts if they are complex, to form the columns of $G_{k}\in\mathbb{R}^{m\times k}$. Both the real and the imaginary parts need to be included. Then we compute the reduced QR factorization of $G_{k}$ : $G_{k}=Q_{k}\Gamma_{k}$, where $Q_{k}=[{\bf q}_{1},\ldots,{\bf q}_{k}]$. Notice that both $G_{k}$ and $Q_{k}$ are real. 10. 10.

Extend the vectors ${\bf q}_{1},\ldots,{\bf q}_{k}$ to length $m+1$ with zero entries, then orthonormalize the vector ${\bf r}_{m}={\bf c}-\bar{H}_{m}{\bf y}^{w}_{m}$ against the columns of \left[\begin{array}[]{c}Q_{k}\\ 0\end{array}\right] to form ${\bf q}_{k+1}.$ Set Q_{k+1}=\left[\begin{array}[]{c|c}\begin{array}[]{c}Q_{k}\\ 0\end{array}&{\bf q}_{k+1}\end{array}\right]. 11. 11.

Compute $\mathcal{V}_{k+1}^{new}=\mathcal{V}_{m+1}(Q_{k+1}\otimes I_{s})$ and $\bar{H}_{k}^{new}=Q_{k+1}^{T}\bar{H}_{m}Q_{k}$, note that $\bar{H}_{k}^{new}$ is generally a full matrix. Update $D$ according to Options 1–3. 12. 12.

*Run the weighted global Arnoldi algorithm from index $(k+1)$ to $m$ to obtain $\mathcal{V}_{m+1}$ and $\bar{H}_{m}$, where the *$(k+1)$th block is the last $s$ columns of ${\mathcal{V}}_{k+1}^{new}$. 13. 13.

Set $X_{0}=X_{m},R_{0}=R_{m}$ and ${\bf c}=\mathcal{V}_{m+1}^{T}\diamond_{D}R_{0}$, and go to Step 6.

4 Numerical Experiments

In this section, we perform some numerical experiments to show the potential of our new algorithms for solving large Sylvester matrix equations. All the numerical examples were performed using MATLAB R2013b on PC-Pentium(R) with CPU 2.66 GHz and 4.00 of RAM. In all the algorithms, we set $X_{0}=0_{n\times s}$ to be the initial guess, and choose $C=sprand(n,s,s)$ as the right-hand side matrix, where $sprand(n,s,s)$ is the MATLAB command generating a random, $n$-by-$s$, sparse matrix with approximately $s\times n\times s$ uniformly distributed nonzero entries. Moreover, we use the stopping criterion

[TABLE]

and all the algorithms will be stopped as soon as the maximal iteration number $maxit=2500$ is reached. For all the algorithms, we consider comparisons in three aspects: the number of iterations (referred to iter), the runtime in terms of seconds (referred to CPU) and the “real” residual in terms of Frobenius norm (referred to res.norm) defined as

[TABLE]

where $X_{m}$ are the computed solutions from the algorithms. In the tables below, we denote by $m,k$ the number of global Arnoldi process and the number of harmonic Ritz block vectors added to the search subspace, respectively.

Example 1. In this example, we illustrate the numerical behavior of Algorithm 2 (W-GLGMRES) for different choices of $D$, and show efficiency of our three new weighting strategies proposed in Options 1–3. To this aim, we compare the performance of W-GLGMRES with the standard global GMRES algorithm (GLGMRES) proposed in [27].

The matrices $A$ and $B$ are obtained from the discretization of the operators [1]

[TABLE]

on the unit square $[0,1]\times[0,1]$ with homogeneous Dirichlet boundary conditions. In these operators, we have $f_{1,1}(x,y)=e^{x^{2}+y},f_{1,2}(x,y)=2xy,f_{2,1}(x,y)=\sin(x+2y),f_{2,2}(x,y)=e^{xy}$, $f_{3,1}(x,y)=\cos(xy)$ and $f_{3,2}(x,y)=xy.$ The dimensions of matrices $A$ and $B$ are $n=n_{0}^{2}$ and $s=s_{0}^{2}$, respectively. By using the command $fdm\_2d\_matrix$ from LYAPACK [29], we can extract the matrices $A=fdm(f_{1,1},f_{2,1},f_{3,1})$ and $B=fdm(f_{1,2},f_{2,2},f_{3,2})$.

We make use of three cases for $D$, i.e., $D_{1},D_{2}$ and $D_{3}$, which are proposed in Options 1–3. Note that they could be updated during the cycles. We also consider the case of $D=I_{n\times n}$ in which Algorithm 2 reduces to the standard GLGMRES algorithm for large Sylvester matrix equations [27]. Table 1 lists the numerical results for different choices of $m,s$ and $n$; and Figure 1 plots the convergence curves of the algorithms for $n=22500$ and 40000 as $s=16,m=15$.

From Table 1 and Figure 1, we observe that the three weighted GLGMRES algorithms need much fewer iterations and much less CPU time than the standard GLGMRES algorithm, and they reach about the same accuracy in terms of the “real” residual norm. More precisely, W-GLGMRES performs better than the standard GLGMRES algorithm, using $D_{1}$, $D_{2}$ or $D_{3}$ as the weighting matrix; and it works the best with $D_{3}$. All these demonstrate that the WGLGMRES algorithm has the potential to improve the convergence, and also it is more robust and efficient than the standard global GMRES algorithm.

Example 2.  In this example, we compare our weighting strategies with the ones proposed in [14, 31], and show effectiveness of our new strategies. In [14], Heyouni et al. considered the linear equation $AX=C$, and proposed a weighted matrix $D$ with elements $(D)_{i,j}=\sqrt{ns}|C|/\|C\|_{F}$, where $|C|$ is the matrix with absolute values of $C$. They then introduced a weighted strategy as $(X,Y)_{D}=tr(X^{T}(D\circ Y))$, where $D\circ Y$ denotes the Hadamard product of $D$ and $Y$. In [31], Saberi Najafi et al. suggested choosing the weighting matrix $D$ as a diagonal random matrix whose diagonal elements are uniformly and randomly chosen from $(0,2)$. In all the numerical examples from now on, we use $D_{3}$ as the weighting strategy in our new algorithms.

The test matrices are available from the Matrix Market Collection [22]. The names of these matrices, the size, the density of nonzeros elements and the type of the matrices are shown in the Table 2. Table 3 lists the iteration numbers, CPU time and residual norms of the approximations, obtained from running W-GLGMRES with three different weighting strategies. The results demonstrate that by using our weighting strategy $D_{3}$, the weighted global GMRES converges faster, and it needs fewer number of iterations and less CPU time than the other two strategies given in [14, 31]. In this example, the “Hadamard product” strategy [14] is better than the “randomized” strategy [31] according to iteration numbers and CPU time, while our new strategy based on the residual works the best. However, we find that the “real residual” norm res.norm computed from the “Hadamard product” strategy [14] may be larger than the desired tolerance $tol=10^{-6}$ in some cases, and it is obvious to see that our new strategy can cure this drawback very well. Indeed, the stopping criterion used is (34) in practical calculations, rather than (35). Figures 2 and 3 plot the convergence curves of the three algorithms. Again, they illustrate the effectiveness and efficiency of our new choice of weighting matrix.

Example 3. When $D=I$, Algorithm 3 reduces to the global GMRES algorithm with deflation, which is mathematically equivalent to the algorithm proposed in [20]. In this example, we try to show that the weighted global GMRES with deflation is more efficient than the global GMRES algorithm with deflation. To show the efficiency of Algorithm 3 (W-GLGMRES-D), we compare it with the global GMRES algorithm (GLGMRES), Algorithm 2 (W-GLGMRES), and the global GMRES algorithm with deflation (GLGMRES-D). In the first test problem, we use $B=fdm(\cos(xy),e^{y^{2}x},100)$, and for test problems 2–5, we use $B=fdm(\sin(xy),e^{xy},10)$.

Table 4 reports the results of the five test problems, where we use different $m$ and $k$ in W-GLGMRES-D and GLGMRES-D. It is seen that both W-GLGMRES-D and W-GLGMRES outperforms GLGMRES-D and GLGMRES in most cases, which illustrates the effectiveness of our weighting strategy. Furthermore, W-GLGMRES-D is superior to the other three algorithms in terms of iteration numbers, CPU time, and accuracy. Specifically, for the 4-th test problem, we see GLGMRES and W-GLGMRES fail to converge with in 2500 iterations, while the algorithms with deflation work quite well. This illustrates that the deflation strategy can improve convergence of the standard global GMRES algorithms for large Sylvester matrix equations. In addition, in Figures 4 and 5, we compare GLGMRES, WGLGMRES, GLGMRES-D and W-GLGMRES-D, where $m=20,k=10$. From Table 4 and Figures 4 –5, we conclude that applying deflation strategy with weighting technique leads to much better solutions.

**Example 4. **  In this example, we combine the weighted and deflated strategy with the flexible preconditoning strategy [30], and show the numerical behavior of the resulting algorithm. In the flexible preconditioned algorithms, the preconditioner may vary from one step to the next, for more details, refer to [30]. In this example, the flexible preconditioner consists of five steps of full GLGMRES for solving the linear systems with multiple right-hand sides in the inner iterations, and we make use of $D_{3}$ as the weighting strategy in our new algorithm.

We first compare GLGMRES, W-GLGMRES with their flexible preconditioning versions: flexible global GMRES (F-GLGMRES), and weighted flexible global GMRES (WF-GLGMRES). Table 5 reports the number of iterations, CPU time and residual norm of the four algorithms. It is seen that by combining the weighted and flexible strategies together, we pay fewer iterations and less CPU time, compared with the standard global GMRES algorithm and the weighted global GMRES algorithm, except for Sherman4. Indeed, for this problem, both F-GLGMRES and WF-GLGMRES use fewer iterations than W-GLGMRES, while the CPU time used for the two former algorithms is (much) more than W-GLGMRES. The reason is that in the flexible algorithms, one has to approximately solve $m=20$ linear systems with $s=400$ right-hand sides per cycle. Note that for Sherman2, both the standard and the weighted flexible global GMRES do not work within 2500 iterations. Figures 6–7 plot the convergence behavior of WF-GLGMRES and F-GLGMRES for the four problems.

Next, we compare the weighted flexible global GMRES with deflation with flexible global GMRES with deflation [10] for solving the five problems, where $m=20,k=15$ is used. Table 6 lists the numerical results. Figures 8–9 plot the convergence curves of the two algorithms during iterations. Again, it is obvious to see that the weighted algorithm is better than the standard one in terms of iteration numbers and CPU time. Compared with the numerical results given in Table 4, we find that the flexible and deflated algorithms often need fewer iterations than the deflated versions, however, the CPU time of the former can be much more than the latter. As we have pointed out before, this is due to the fact that the inner iterations bring us a large amount of computational overhead. How to reduce the high cost from inner iterations is beyond the scope of this paper, but deserves further investigation.

Moreover, the two flexible and deflated algorithm still do not work for Sherman2, just like the bare flexible algorithms. One reason is that only five steps of full GLGMRES for solving the linear systems in the inner iterations is not enough for this problem. Thus, we suggest to use deflated global GMRES when $s$, the number of columns of $C$, is large, say, more than one hundred. On the other hand, when $s$ is of medium-sized, we recommend to use the flexible and deflated global GMRES algorithm.

5 Conclusion

The global GMRES algorithm is popular for large Sylvester matrix equations. The weighting strategy can improve the algorithm by alienating the eigenvalues that obstacle the convergence. However, the optimal choice of the weighting matrix is still an open problem and needs further investigation. Moreover, due to the growth of memory requirements and computational cost, it is necessary to restart the algorithm efficiently.

The contribution of this work is three-fold. First, we present three new schemes based on residual to update the weighting matrix during iterations, and propose a weighted global GMRES algorithm. Second, we apply the deflated restarting strategy to the weighted algorithm, and propose a weighted global GMRES algorithm with deflation for solving large Sylvester matrix equations. Third, we show that in the weighted and deflated global GMRES algorithm, there is no need to change the inner product with respect to diagonal matrix to that with non-diagonal matrix, and our scheme is much cheaper than the one proposed in weighted GMRES-DR algorithm [8]. Further, we consider other acceleration technology such as the flexible preconditioning strategy. For the weighted flexible global GMRES algorithm with deflation, it is interesting to reduce the high cost from inner iterations, and it is definitely a part of our future work.