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

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 [ M. Embree, R. B. Morgan and H. V. Nguyen, 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.

## Authors

• 1 publication
• 1 publication
• 12 publications
11/25/2021

### The intrinsic Toeplitz structure and its applications in algebraic Riccati equations

In this paper we propose a new algorithm for solving large-scale algebra...
11/01/2020

### Highly accurate decoupled doubling algorithm for large-scale M-matrix algebraic Riccati equations

We consider the numerical solution of large-scale M-matrix algebraic Ric...
11/03/2020

### Decoupled Structure-Preserving Doubling Algorithm with Truncation for Large-Scale Algebraic Riccati Equations

In Guo et al, arXiv:2005.08288, we propose a decoupled form of the struc...
08/20/2019

### A generalized optimal fourth-order finite difference scheme for a 2D Helmholtz equation with the perfectly matched layer boundary condition

A crucial part of successful wave propagation related inverse problems i...
03/21/2020

### Parameter robust preconditioning for multi-compartmental Darcy equations

In this paper, we propose a new finite element solution approach to the ...
04/09/2021

### Householder orthogonalization with a non-standard inner product

Householder orthogonalization plays an important role in numerical linea...
08/20/2020

### Efficient and Accurate Algorithms for Solving the Bethe-Salpeter Eigenvalue Problem for Crystalline Systems

Optical properties of materials related to light absorption and scatteri...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Consider the large Sylvester matrix equation

 AX+XB=C, (1)

where , and , with . If we define the operator as

 A : Rn×s⟶Rn×s, X→ AX=AX+XB. (2)

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

 AX=C. (3)

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

 (Is⊗A+BT⊗In)vec(X)=vec(C), (4)

where denotes the Kronecker product operator, and

denotes the vectorize operator defined as (in MATLAB notation)

 vec(X)=[X(:,1);X(:,2);…;X(:,s)],

and is the -th column of the matrix . The linear equation systems (4) have unique solution if and only if the matrix 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 -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 -orthonormal basis for the Krylov subspace, where 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 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 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 be a diagonal matrix with , and let be given, then the -inner product with respect to two vectors is defined as [36]

 (u,v)D=vTDu =n∑i=1diuivj,

and the associated -norm of is defined as

 ∥u∥D =√(u,u)D,∀u∈Rn.

For two matrices , the -inner product is defined as [26]

 (Y,Z)D= trace(ZTDY),

where denotes the trace of a matrix, and represents the transpose of the matrix . It can be verified that [26]

 (Y,Z)D=(vec(Y),vec(Z))Is⊗D.

Also, the -norm associated with this inner product is

 ∥Y∥D =√(Y,Y)D,∀ Y∈Rn×n.

Next we introduce a useful product that will be used latter: [26] Let and , where . Then elements of the matrix is defined as

 (AT⋄DB)ij=(Ai,Bj)D,i=1,2,…,p, j=1,2,…,l. (5)

It was shown that [26]. Let be an initial block vector that is -orthogonal, that is, orthonormal with respect to the -inner product. The following algorithm presents an -step weighted global Arnoldi process [26].

###### Algorithm 1

The -step weighted global Arnoldi process

1. Input: , , , a positive diagonal matrix and an integer number .

2. Output:

3. Compute ,

4. for

5. Compute

6. for

7. end

8. Compute . If then stop.

9. Compute

10. end

The weighted global Arnoldi process constructs a -orthonormal basis , i.e.,

 (Vi,Vj)D =δij,∀i,j=1,2,…,m,

for the matrix Krylov subspace

 Km(A,V1) := span{V1,AV1,…,Am−1V1} = {m−1∑i=1αiAi−1V1|αi∈R,i=1,2,…,m−1}.

Let be a quasi-upper Hessenberg matrix whose nonzeros entries are defined by Algorithm 1, and the matrix is obtained from the matrix by deleting its last row. Note that the matrix is -orthogonal. With the help of Definition 2, we obtain the following relations

 [AV1,AV2,…,AVm] = Vm(Hm⊗Is)+hm+1,mVm+1(eTm⊗Is), (6) = Vm+1(¯Hm⊗Is), (7) ¯Hm = VTm+1⋄DAVm, (8)

where . Define

 AVm≡[AV1,AV2,…,AVm], (9)

then (6) can be rewritten as

 AVm=Vm+1(¯Hm⊗Is). (10)

We are in a position to consider the weighted global GMRES algorithm for solving (1). Let be the initial guess, and the initial residual be . In the weighted global GMRES algorithm, we construct an approximate solution of the form

 Xm=X0+Vm(ywm⊗Is), (11)

where . The corresponding residual is

 Rm = C−AXm−XmB (12) = (C−AX0−X0B)−(AVm(ywm⊗Is)+Vm(ywm⊗Is)B), = R0−AVm(ywm⊗Is),

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

 Rm = βV1−Vm+1(¯Hm⊗Is)(ywm⊗Is), (13) = Vm+1((βe1−¯Hmywm)⊗Is),

where is the first canonical basis vector in . Note that the residual is -orthogonal to , i.e.,

 Rm=C−AXm⊥DAKm(A,R0), (14)

where , and “” means orthgonal with respect to the “” inner product.

In order to compute , we have from (13) and (14) that

 miny∥Rm∥D = miny∥Vm+1((βe1−¯Hmy)⊗Is)∥D = miny∥(βe1−¯Hmy)T(Vm+1)TD(Vm+1)(βe1−¯Hmy)∥2 = miny∥βe1−¯Hmy∥2,

where we used . Thus,

 ywm=argminyw∈Rm∥βe1−¯Hmyw∥2, (15)

or equivalently,

 ¯HTm¯Hmywm=¯HTmβe1. (16)

As was pointed out in [9, 14, 36], the optimal choice of 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 , which could be updated during iterations:

Option 1:

Let , where is the -th column of residual matrix . Then we define , where stands for the absolute value of .

Option 2:

Similarly, let , then we define .

Option 3:

Use the mean of the block residual , i.e, .

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)

1. Input: , , . Choose the initial guess, , an initial positive diagonal matrix D and an integer , and the tolerance .

2. Output: The approximation .

3. Compute and ,

4. Run Algorithm 1 with the initial block vector to obtain the matrices .

5. Solving for .

6. Compute and . If , then stop, else continue.

7. Set and update the weighting matrix 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 for (1) can be understood as the weighted GMRES algorithm with the weighting matrices

 ˆDj=Is⊗Dj, j=1,2,3,

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 weighted harmonic Ritz pairs. Let be the -orthonormal basis obtained from the “previous” cycle, we seek weighted harmonic Ritz pairs that satisfy

 AVmYi−θiVmYi⊥D(A−σI)Km(A,V),i=1,…,k, (17)

where

 AVm=[AV1,AV2,…,AVm]=[AV1+V1B,…,AVm+VmB],

and with . In this work, we want to deflate some smallest eigenvalues in magnitude, and a shift is used throughout this paper.

From (10), we have that

 AVm(gi⊗Is)−θiVm(gi⊗Is) = (AVm−θiVm)(gi⊗Is) = Vm+1((¯Hm−θi¯Im)⊗Is)(gi⊗Is),

where By (17) and Definition 2, we can compute via solving the following (small-sized) generalized eigenvalue problem

 θi((AVm)T⋄DVm)gi=(AVm)T⋄D(AVm)gi. (18)

From (10) and the fact that , we rewrite (18) as

 θiHTmgi=¯HTm¯Hmgi. (19)

If is nonsingular, (19) is equivalent to

 (Hm+h2m+1,mH−TmemeTm)gi=θigi. (20)

Then we define the “weighted harmonic Ritz vector” as , and the corresponding harmonic residual is

 ˜Ri=AVmYi−θiVmYi = Vm+1((¯Hm−θi¯Im)gi⊗Is) = Vm+1(˜ri⊗Is),i=1,…,k,

where .

###### 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 .

In the following, we characterize the relationship between the weighted harmonic residuals and the residual from the weighted global GMRES algorithm. Let be the initial guess and be the initial residual. After one cycle of the weighted global GMRES Algorithm, we have the approximate solution , where is defined in (15). The associated residual with respect to is

 Rm=Vm+1((βe1−¯Hmywm)⊗Is)=Vm+1(rm⊗Is),

where . The following result shows that and are collinear with each other. Let , be the weighted harmonic residuals and be the weighted global GMRES residual, respectively, where and Then and are collinear with each other. Note that both the weighted harmonic residuals and the residual of the weighted global GMRES algorithm are in , and they are both -orthogonal to . Thus, there is an matrix , such that

 Vm+1(˜ri⊗Is)=Vm+1(rm⊗Is)T, (21)

so we have . Let and , we obtain from (21) that

 (rm⊗Is)T=⎛⎜ ⎜ ⎜ ⎜⎝η1Isη2Is⋮ηm+1Is⎞⎟ ⎟ ⎟ ⎟⎠T=⎛⎜ ⎜ ⎜ ⎜⎝η1Tη2T⋮ηm+1T⎞⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜⎝τ1Isτ2Is⋮τm+1Is⎞⎟ ⎟ ⎟ ⎟⎠,

which implies that and 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 , and let be the reduced QR factorization. We stress that both forming

and computing the QR decomposition can be implemented in real arithmetics; see Step 9 of Algorithm

3. Then we orthonormalize against to get , and let .

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

 AVm=Vm(Hm⊗Is)+hm+1,mVm+1(eTm⊗Is),

and thus

 AVm(gi⊗Is)=Vm(Hm⊗Is)(gi⊗Is)+hm+1,mVm+1(eTmgi⊗Is), i=1,…,k. (22)

That is,

 AVm(gi⊗Is)=[Vm Vm+1][Hmgihm+1,meTmgi]⊗Is.

We note that

 ˜ri=(¯Hm−θi¯Im)gi=[(Hm−θiI)gihm+1,meTmgi],

and

 [Hmgihm+1,meTmgi] = [Hmgi−θigi+θigihm+1,meTmgi] = θi[gi0]+˜ri, i=1,…,k.

As a result,

 [Hmgihm+1,meTmgi]∈span{[Gk0],rm}=span{Qk+1}, i=1,…,k,

where we used and are collinear with each other; see Theorem 3. Therefore,

 AVm(gi⊗Is)∈span{Vm+1(Qk+1⊗Is)},  i=1,…,k,

and

 AVm(Qk⊗Is)⊆span{Vm+1(Qk+1⊗Is)}. (23)

Define , where , so we have

 AVm(Qk⊗Is)=AVnewk=A[Vnew1,…,Vnewk].

If we denote

 Vnewk+1 = Vm+1(Qk+1⊗Is),

then we have from (23) that

 AVm(Qk⊗Is)⊆span{Vnewk+1},

and there is a matrix such that

 AVm(Qk⊗Is)=Vm+1(¯Hm⊗Is)(Qk⊗Is)=Vm+1(Qk+1⊗Is)P=Vnewk+1P.

The condition yields

 P = (Qk+1⊗Is)TVTm+1⋄DAVm(Qk⊗Is) = (Qk+1⊗Is)TVTm+1⋄DVm+1(¯Hm⊗Is)(Qk⊗Is), = QTk+1¯HmQk⊗Is=¯Hnewk⊗Is,

where is generally a dense matrix. In conclusion, we obtain

 AVnewk=Vnewk+1(¯Hnewk⊗Is). (24)

Then, we run the weighted global Arnoldi Algorithm from index to with the last columns of as the starting matrix to build a new -step global Arnoldi relation. However, is updated in our new algorithm, see Step 11 of Algorithm 3. In other words, we will perform the remaining steps of global Arnoldi process with a new (denoted by that is from the “current” residual ) after deflated restarting.

Denote , where , then is -orthogonal and it is also orthogonal to . Denote by the weighting matrix obtained from the residual of the “previous” cycle, we define -orthogonality of as follows

 VTm⋄˜DVm≡⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩VnewTk+1⋄DoldVnewk+1=Ik+1,VnewT⊥⋄DnewVnew⊥=Im−k−1,VnewT⊥⋄DnewVnewk+1=0. (25)

That is, , 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 -inner products with respect to non-diagonal matrices. Thanks to (25), 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.

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

 AVm = Vm(Hm⊗Is)+hm+1,mVm+1(eTm⊗Is), (26) = Vm+1(¯Hm⊗Is). (27)

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)

1. Input: . Choose an initial guess a positive diagonal matrix , the positive integer number and a convergence tolerance .

2. Output: The approximation .

3. Compute , and

4. Set

5. Run the weighted global Arnoldi algorithm to obtain and .

6. Solve for .

7. Compute , and .

8. If , then stop, else continue.

9. Compute the eigenpairs with the smallest magnitude from (19) or (20). Set : We first sperate the s into real and imaginary parts if they are complex, to form the columns of . Both the real and the imaginary parts need to be included. Then we compute the reduced QR factorization of : , where . Notice that both and are real.

10. Extend the vectors to length with zero entries, then orthonormalize the vector against the columns of to form Set .

11. Compute and , note that is generally a full matrix. Update according to Options 1–3.

12. Run the weighted global Arnoldi algorithm from index to to obtain and , where the th block is the last columns of