# A Flexible Global GCRO-DR Method for Shifted Linear Systems and General Coupled Matrix Equations

Authors: Jing Meng; Xian-Ming Gu; Wei-Hua Luo; Liang Fang
Journal: Journal of Mathematics
 (2021)

---

## Abstract

In this paper, we mainly focus on the development and study of a new global GCRO-DR method that allows both the flexible preconditioning and the subspace recycling for sequences of shifted linear systems. The novel method presented here has two main advantages: firstly, it does not require the right-hand sides to be related, and, secondly, it can also be compatible with the general preconditioning. Meanwhile, we apply the new algorithm to solve the general coupled matrix equations. Moreover, by performing an error analysis, we deduce that a much looser tolerance can be applied to save computation by limiting the flexible preconditioned work without sacrificing the closeness of the computed and the true residuals. Finally, numerical experiments demonstrate that the proposed method illustrated can be more competitive than some other global GMRES-type methods.


---

## Body

## 1. Introduction

In the current study, we focus on the solution of sequences of shifted linear systems of the form(1)Ai−σiℓIxiℓ=biℓ,fori=1,2,…,s;ℓ=1,2,…,pi,where Aii=1s⊂ℂn×n are large nonsingular matrices of order n, the right-hand sides biℓℓ=1pi⊂ℂn are vectors of the length n, and σiℓℓ=1pi⊂ℂ are called shifts. Such shifted linear systems (1) arise in many realistic applications, for example, electromagnetic scattering [1], QCD [2–5], dynamics of structures [6], complex network [7–9], and model reduction [10].Note that problem (1) is reduced to the solution of a sequence of linear systems with the form(2)A˜ixi=bi,fori=1,2,…,s,when ℓ=1, where the coefficient matrix A˜i=Ai−σiI and the right-hand sides bi change from one system to the next. One approach to more efficient solution of such systems is to exploit Krylov subspace recycling; that is, recycling a judiciously selected subspace of the search space typically reduces the number of iterations substantially. There are several algorithms that take advantage of recycling for subsequent systems as well as for restarts of generalized minimal residual (GMRES) type methods, such as the generalized conjugate residual method with inner orthogonalization (GCRO) [11], GCRO with deflated restarting (GCRO-DR) [12], GMRES with deflated restarting (GMRES-DR) [13] (only for a single system), and GCRO with optimal truncation (GCROT) [14].It is also noted that, for eachi, we substantially solve a family of shifted linear systems with the form(3)A−σℓIxℓ=bℓ,forℓ=1,2,…,pi.Krylov subspace solvers [15] are particularly attractive for the solutions of (3) because of its shift-invariance property [16]. That is, for any shift σ, the m-th Krylov subspace satisfies(4)KmA,b=KmA+σI,b˜,as long as b˜=βb. This property has led to many approaches for solving such linear systems (3) simultaneously by generating only one subspace. Some solvers are built upon the (Bi-)Lanczos recurrences [17–21] and some others use the Arnoldi-like iteration [16, 22–28]. However, the shift-invariance property will be lost when using the general preconditioning, such as Jacobi preconditioning, SSOR preconditioning, ILU preconditioning, AMG preconditioning, and AINV preconditioning. Furthermore, when biℓℓ=1pi are unrelated, we cannot employ the shift-invariance property. It means that all residuals are not simultaneously minimized whiling maintaining collinearity, for example, the GMRES-type methods [16, 23, 29].Therefore, we propose an alternative without exploiting the shift-invariance. In fact, the solution of (3) is equivalent to solve a Sylvester equation [30],(5)AX+XT=B,where T=diagσ1,σ2,…,σpi∈ℂpi×pi and B=b1,b2,…,bpi∈ℂn×pi. Therefore, (1) can also be reformulated as a family of Sylvester equations with the slowly changing coefficient matrices(6)AiXi+XiTi=Bi,i=1,2,…s.For simplicity of the presentation, a linear operatorA is defined as follows:(7)A:X⟶AX+XT.Applying the linear operatorA, we can rewrite (1) as(8)AiXi=Bi.Moreover, we also assume that the linear operatorA is invertible; that is,(9)Ai−1Bi=Xi.Obviously, the linear systems (8) can be seen as generalized version of (2). Therefore, based on the special structure of (1), we propose a new generalized global version of the GCRO-DR method that not only combines with any type of the preconditioning but also allows one to retain important spectral information generated during the solution of the i-th linear system and then to exploit such information to accelerate convergence of the iterations when solving the subsequent i+1-th system. Moreover, the proposed method can be combined with any form of preconditioning and does not require the right-hand sides to maintain collinearity.For simplicity, we denoteAX by AX. In the following, we utilize the operator A to present our flexible global GCRO-DR algorithms for solving (1).The paper is organized as follows. A global version of GCRO-DR to solve (1) is derived in detail in Section 2. We also apply the flexible preconditioning operator at each iteration to improve the global GCRO-DR algorithm. By performing an error analysis, we deduce that a much looser tolerance can be applied to save computation by limiting the flexible preconditioned work without sacrificing the closeness of the computed and the true residuals. Meanwhile, we also propose a global version of the FGMRES-DR method and then analyze the relationship between two flexible global methods. In Section 2.5, the flexible global GCRO-DR method is applied to solve the general coupled matrix equations. Section 3 numerically demonstrates the effectiveness of the novel method. Finally, some conclusions and remarks are summarized in Section 4.


## 2. Flexible Global GCRO-DR Method

GCRO-DR is a well-established Krylov subspace solver that paired with a harmonic Ritz vector recycling strategy [12]. GCRO [11] is a nested Krylov method, which uses GCR as the outer solver and the restarted GMRES as an inner method. In [12], Parks et al. propose GCRO-DR that exploits deflated restarting within the framework of GCRO for a family of slowly changing linear equations. Because of its superior numerical behavior, we introduce a global version of GCRO-DR to solve (1). Here, we not only employ the subspace recycling at each cycle but also apply the flexible preconditioning operator at each iteration. The proposed method is referred to as flexible global GCRO-DR (FGL-GCRO-DR).


### 2.1. Flexible Global Arnoldi Process

We first present a global Arnoldi process used in the flexible setting to create an F-orthonormal basis for the global Krylov subspace in Algorithm1. We define ℳj as the preconditioning operator used at the j-th iteration, and the action ℳj on rectangular matrix V is denoted by ℳjV.The flexible global Arnoldi process, recursively, generates anF-orthonormal basis Vm+1=V1,V2,…,Vm+1 of the global Krylov subspace spanned by the rectangular matrices V1,AZ1,…,AZm−1. At the end of the m-th iteration, we obtain a matrix relation (flexible global Arnoldi relation); that is,(10)AZm=Vm+1H¯m⊗Ip,where Zm=ℳ1V1,…,ℳmVm∈ℂn×mp, Vm+1∈ℂn×m+1p, and(11)H¯m=h1,1h1,2……h1,mh2,1h2,2⋮h3,2⋱⋮⋱⋱⋮hm−1,mhm,mhm+1,m=Hmhm+1,memT∈ℂm+1×m,is an upper Hessenberg matrix.


### 2.2. The FGL-GCRO-DR Algorithm

Suppose thatZk=Z1,…,Zk and Vk=V1,…,Vk∈ℂn×kp are given matrices satisfying(12)AZk=VkIk⊗Ip,Vk,VkF=Ik,where Zi,i=1,2,…,k is a search vector, recycled from either the previous linear system or the previous cycle, and 0≤k≤m. By minimizing the residual F-norm,(13)minX∈X0+ℛZkR0−AXF,over the search space X0+ℛZk, Xk and Rk∈ℂn×p satisfy(14)Xk=X0+ZkVk,R0F⊗Ip,Rk=R0−VkVk,R0F⊗Ipsuch thatRk⊥FℛVk.Based upon what is proven in [11, 14] and using the framework discussed in [12, 32], the m−k steps of the flexible global Arnoldi process are carried out, with the starting rectangular matrix Vk+1=Rk/Rk, while holding the orthogonality to Vk, to eventually obtain(15)AZk+1,Zk+2,…,Zm−VkVk,AZk+1,Zk+2,…,Zm⊗Ip=Vk+1,Vk+2,…,Vm+1H¯m−k⊗Ip.SinceRk⊥FℛVk, that is, Vk+1⊥FVk, the following flexible global Arnoldi-like relation holds:(16)AZk,Zk+1,…,Zm=V1,V2,…,Vm+1Ik⊗IpVk,AZk+1,…,Zm⊗Ip0H¯m−k⊗Ip,AZm=V1,V2,…,Vm+1⋅IkVk,AZk+1,…,Zm0H¯m−k⊗Ip,AZm=VmH¯m⊗Ip,Vm+1,Vm+1F=Im+1.At the end, we find the approximate solutionXm=X0+Zmy∗⊗Ip∈ℂn×p by minimizing the norm R0−AZmy⊗IpF, where y∗ satisfies(17)y∗=argminy∈ℂmR0−AZmy⊗IpF=argminy∈ℂmβ0−H¯my2,with β0=Vm+1,R0F.To improve the convergence of iterative solvers, the smallest magnitude eigenvalues will be deflated. Following the idea of [12, 32], we use harmonic Ritz vectors with those eigenvalues to construct approximate invariant subspace and then the subspace for next cycles is retained. Firstly, we compute k harmonic Ritz vectors y˜k=Wmgk⊗Ip, such that(18)AZmWmHFWmgi⊗Ip−θiWmgi⊥FAZm,where Wm∈ℂn×mp, and each harmonic Ritz pair θi,gi satisfies(19)H¯mHH¯mgi=θiH¯mHVm+1HWmgi,1≤i≤k.It is noted that the harmonic residual vectors satisfy(20)AZmWmHFWmgi⊗Ip−θiWmgi∈ℛVm+1.Observe that the residualRm=R0−AZmy∗⊗Ip=Vm+1y∗⊗Ip∈ℛVm+1⊥FℛAZm. So, the residual and the harmonic residuals are in the same p-dimensional space ℛAZm⊥F∩ℛVm+1. Then we can characterize the relationship by the following formula:(21)AZmGk⊗Ip=WmGk⊗Ip,Rm⋅diagθ1,…,θka1×k⊗Ip,where a:,i satisfies AZmgi⊗Ip−θiVmgi⊗Ip=a1,iRm, 1≤i≤k, and Gk=g1,…,gk∈ℂm×k.Due to the relations (10) and (21) can be also expressed as(22)Vm+1H¯m⊗IpGk⊗Ip=Vm+1H¯mGk⊗Ip=WmGk⊗Ip,Rm⋅diagθ1,…,θka1×k⊗Ip.Next, we compute the reduced QR-factorization ofH¯mGk=QR with Q∈ℂm+1×k and R∈ℂk×k and then define Zk, Vk∈ℂn×k as(23)Zk=ZmGkR−1⊗Ip,Vk=Vm+1Q⊗Ip,such that AZk=Vk. Since Rm is F-orthogonal to the subspace ℛAZm and Vm+1 is F-orthonormal, we have(24)Vk,RmF=Vm+1Q⊗IpR⊗Ip,RmF=Vm+1H¯mGk⊗Ip,RmF=AZmGk⊗Ip,RmF=0,Vk,VkF=Vm+1Q⊗Ip,Vm+1Q⊗IpF=QHQ=Ik.Therefore, at the end of the loop, an updateZk is generated, scaled such that Vk=AZk has F-orthonormal columns; the previous global Krylov subspace basis is dropped, and then we restart. Details of the FGL-GCRO-DR method are summarized in Algorithm 2.Algorithm 1: The FGL-GCRO-DR method.
Input:A∈ℂn×n, B∈ℂn×p. Choose m, k, tol and initial guesses X0.Output:Xm such that RmF/R0F≤tol.(1)
R0=B−AX0(2)
WhileRmF>R0F×toldo(3)
Ifm==1(4)
γ=R0F, V1=R0/γ(5)
Apply Algorithm2 to generate Zm, H¯m, Vm+1 such that AZm=Vm+1H¯m⊗Ip, and setWm=Vm(6)
Else(7)
Applym−k steps of the Algorithm 2 with the starting rectangular matrix Vk+1=Rk/RkF, such that AZk+1,…,Zm−VkVk,AZk+1,…,Zm⊗Ip=Vm−k+1H¯m−k⊗Ip with Zj=ℳj−1Vj(8)
Endif(9)
y∗=argminy∈ℂmγ−H¯my2 with γ=Vm+1,R0F(10)
SetXm=X0+Zmy∗⊗Ip, Rm=B−AXm(11)
Compute thek smallest eigenpairs θi,gi of H¯mHH¯mgi=θiH¯mHVm+1HWmgi,1≤i≤k,(12)
QR=H¯mGk(13)
Vk=Vm+1Q(14)
Zk=ZmGkR−1(15)
Wk=WmGkR−1(16)
Wm=Wk,Vk:,k+1p:mp(17)
EndwhileAlgorithm 2: Flexible global Arnoldi process.
γ=R0F and set V1=R0/γ;forj=1,…,m, doZj=ℳjVjV^=AZjfori=1,…,j, dohij=Vi,V^F, i=1,…,j;V^=V^−Vihij;end forhj+1,j=V^F and set Vj+1=V^/hj+1,jend forIn addition, the major advantage of FGL-GCRO-DR is that it can retain important spectral information generated during the solution of thei-th linear system and then exploit such information to accelerate convergence of the iterations when solving the subsequent i+1-th system. Then we detail how to recycle subspace when solving sequences of slowly changing linear systems (8) as follows (Algorithm 3).Algorithm 3: FGL-GCRO-DR for solving sequences of slowly changing linear systems.
(1)
Suppose thatVki, Zki and Wki are defined from solving i-th linear system linear and satisfy AiZki=Vki with Vki,Vki=Ik.(2)
R0i+1=Bi+1−Ai+1X0i+1(3)
QR=Ai+1Zki(4)
Vki+1=Q(5)
Zki+1=ZkiR−1⊗Ip(6)
Wki+1=WkiR−1⊗Ip(7)
X0i+1=X0i+1+Zki+1Vki+1,R0i+1⊗Ip(8)
R0i+1=R0i+1−Vki+1Vki+1,R0i+1⊗Ip(9)
Applym−k steps of Algorithm 2 with the starting rectangular matrix Vk+1i+1=Rk/RkF, such that Ai+1Zm−ki+1−Vki+1Vki+1,Ai+1Zm−ki+1⊗Ip=Vm−k+1i+1H¯m−k⊗Ip with Zji+1=ℳj−1Vji+1(10)
y∗=argminy∈ℂmβ−H¯my2 with β=Vm+1i+1,R0i+1F(11)
Xmi+1=Xm−1i+1+Zmi+1y∗⊗Ip(12)
Rmi+1=Bi+1−Ai+1Xmi+1(13)
Wmi+1=Wki+1,Wmi+11:n,k+1p:mp


### 2.3. Comments on the FGL-GCRO-DR Method

To achieve better efficiency of the linear system solvers, we not only employ the subspace recycling at each cycle but also apply the flexible preconditioning operator at each iteration. As we know, the variable preconditioning can be construed as an inexact subspace method. If the matrixA is difficult to invert, we will resort to several steps of iterative method as flexible preconditioning, so we obtain in some way an approximate result,(25)Z^j≈A−1Vj,at j-th (inner) iteration for the equation AZj=Vj,j=1,…,m. In such a situation, what are the properties of the true residual obtained by using the flexible preconditioning? We will address this question in this section. Following [33, 34], we set(26)Pj=Vj−AZ^j,j=1,…,m,to be the residual that results after the iterative solver has been terminated; then(27)Z^j=A−1Vj−A−1Pj,j=1,…,m.Therefore, the perturbation in linear operation withA can be written explicitly as(28)AZ^j=AA−1Vj−AA−1Pj=I+PjVjHFVjF2Vj=I−EjVj.The perturbation implies that the error computed has the bound(29)Ej≤PjFVjF.In fact, the solution is obtained directly from the approximation subspaceℛZ^m; that is,(30)Xm=Z^my∗˜⊗Ip,indicating that the perturbation Ej does not affect the computation of the true residual. In other words, in exact arithmetic(31)Rm=R0−AZ^my˜∗⊗Ip=Vm+1βe1⊗Ip−Vm+1H¯my˜∗⊗Ip=Vm+1βe1−H¯my˜∗⊗Ip,and is orthogonal to ℛVm+1. As a consequence, the perturbation of the linear operations only influences the goodness of ℛZ^m and the accuracy of the approximate solution. Therefore, we deduce that, in the flexible setting, a dynamic stopping criterion or a much looser tolerance can be employed to save computation by limiting the flexible preconditioning work without sacrificing the closeness of the computed and the true residuals.


### 2.4. The Relationship between Two Flexible Global Methods with Subspace Recycling

As we know, FGCRO-DR [32] and FGMRES-DR [35] are two important methods in the general category of augmented (deflated) GMRES-type methods, which choose harmonic Ritz vectors as augmented (deflated) subspace. They emphasize on the solution of linear system Ax=b which has only one equation. The relationship between FGCRO-DR and FGMRES-DR has been well described by [32]. It is shown that they are algebraically equivalent as long as a collinearity condition is satisfied.In fact, FGL-GCRO-DR is a global generalization of FGCRO-DR. For the sake of descriptive integrality, we also present a global generalization of FGMRES-DR algorithm, referred to as FGL-GMRES-DR, in Algorithm 4. As for global versions, the FGL-GCRO-DR method is also closely related to the FGL-GMRES-DR approach in algebraical analysis, but their numerical behavior will be slightly different. Details will be shown in Section 3.Algorithm 4: FGL-GMRES-DR.
Input:A∈ℂn×n, B∈ℂn×p. Choose m, k, tol and initial guessesX0Output:Xm such that RmF/R0F≤tol(1)
R0=B−AX0, γ=R0F, V˜1=R0/γ(2)
Apply Algorithm2 to generate Z˜m, V˜m+1, H¯m such that AZ˜m=V˜m+1H¯m⊗Ip(3)
Solveminγ−H¯my2 for y∗. Set X0=X0+Z˜my∗⊗Im, R0=B−AX0(4)
whileRmF>R0F×toldo(5)
Compute thek smallest eigenpairs θi,gi of Hm+hm+1,mHhm+1,mHm−HememH.(6)
QR=Pk01×kγ−H¯my∗(7)
V˜k=V˜m+1Q(8)
Z˜k=Z˜mQ(9)
H¯k=QHH¯mQ1:m,1:k(10)
Applym−k steps of Algorithm 2 with the starting rectangular matrix V˜k+1=Rk/RkF, such that AZ˜m−k=V˜m−k+1H¯m−k⊗Ip with Z˜j=ℳj−1V˜j(11)
y∗=argminy∈ℂmγ−H¯my2 with γ=Vm+1,R0F(12)
Xm=Xm−1+Z˜my∗⊗Im(13)
Rm=B−AXm(14)
endwhileFurthermore, we compare the main calculation cost and storage capacity for each cycle of FGL-GCRO-DR and FGL-GMRES-DR in Table1. In practical applications, n≫m,k. Here, opM and opA denote floating-point calculation counts for the preconditioning application and the linear operations, respectively.Table 1
The main calculation cost and storage capacity of FGL-GCRO-DR and FGL-GMRES-DR.
Computation ofFGL-GMRES-DR (m, k)FGL-GCRO-DR (m, k)Vk+1:,1:k+1p2m+1k+1np2nm+1kp+npZk:,1:kp2nmkp2nmkpVm:,k+2p+1:m+1popAm−k+m−k2n+2k+5npopAm−k+m−k2n+2k+6npZm:,k+1p+1:mpopMm−kopMm−kTotal costC+m+k−1npCMemory requirement2m+3np3m+3npThe computational cost of harmonic Ritz vectors and “minγ−H¯my2” have not been induced in Table 1 since the cost is approximately Op3, much more less than n. C=opM+opAm−k+2m+1kp+p+2mkp+m−k2m+2k+6np is referred to as the total cost of FGL-GCRO-DR. One sees that the time complexities of two global algorithms are the same, that is, Om−kn2p, but FGL-GMRES-DR requires slightly more computation per cycle than FGL-GCRO-DR in terms of computational cost.


### 2.5. Application to the General Coupled Matrix Equations

FGL-GCRO-DR substantially belongs to the family of the global Krylov subspace method, which can be also applied to solve the general coupled matrix equations of the form(32)∑j=1lAi,1,jX1Bi,1,j+Ai,2,jX2Bi,2,j+⋯+Ai,l,jXlBi,l,j=Di,i=1,2,…,l,including the (coupled) Sylvester and the (coupled) Lyapunov matrix equations. The literature on solving (32) is particularly rich and still growing [36–53], due to the prominent role of such matrix equations in stability theory [54], perturbation analysis [55], control theory [56, 57], and so forth. In particular, the global subspace solvers perform better than other approaches for solving (32), especially for the coefficient matrices which are large and sparse; for more details, see [44–53, 58]. Therefore, we are interested in exploiting the FGL-GCRO-DR method for solving (32). To do so, the linear operator can be defined as follows:(33)ℱ:ℝm×n×⋯×ℝm×n⟶ℝm×pn,X=X1,X2,…,Xl⟶ℱX=A1X,A2X,…,AlX,where Ai=∑k=1l∑j=1lAi,k,jXkBi,k,j,i=1,2,…,l. Exploiting the operator ℱ, we can rewrite (32) as(34)ℱX=D,where D=D1,D2,…,Dl. The effectiveness of FGL-GCRO-DR for the solution of (32) is shown in Section 3.


## 2.1. Flexible Global Arnoldi Process

We first present a global Arnoldi process used in the flexible setting to create an F-orthonormal basis for the global Krylov subspace in Algorithm1. We define ℳj as the preconditioning operator used at the j-th iteration, and the action ℳj on rectangular matrix V is denoted by ℳjV.The flexible global Arnoldi process, recursively, generates anF-orthonormal basis Vm+1=V1,V2,…,Vm+1 of the global Krylov subspace spanned by the rectangular matrices V1,AZ1,…,AZm−1. At the end of the m-th iteration, we obtain a matrix relation (flexible global Arnoldi relation); that is,(10)AZm=Vm+1H¯m⊗Ip,where Zm=ℳ1V1,…,ℳmVm∈ℂn×mp, Vm+1∈ℂn×m+1p, and(11)H¯m=h1,1h1,2……h1,mh2,1h2,2⋮h3,2⋱⋮⋱⋱⋮hm−1,mhm,mhm+1,m=Hmhm+1,memT∈ℂm+1×m,is an upper Hessenberg matrix.


## 2.2. The FGL-GCRO-DR Algorithm

Suppose thatZk=Z1,…,Zk and Vk=V1,…,Vk∈ℂn×kp are given matrices satisfying(12)AZk=VkIk⊗Ip,Vk,VkF=Ik,where Zi,i=1,2,…,k is a search vector, recycled from either the previous linear system or the previous cycle, and 0≤k≤m. By minimizing the residual F-norm,(13)minX∈X0+ℛZkR0−AXF,over the search space X0+ℛZk, Xk and Rk∈ℂn×p satisfy(14)Xk=X0+ZkVk,R0F⊗Ip,Rk=R0−VkVk,R0F⊗Ipsuch thatRk⊥FℛVk.Based upon what is proven in [11, 14] and using the framework discussed in [12, 32], the m−k steps of the flexible global Arnoldi process are carried out, with the starting rectangular matrix Vk+1=Rk/Rk, while holding the orthogonality to Vk, to eventually obtain(15)AZk+1,Zk+2,…,Zm−VkVk,AZk+1,Zk+2,…,Zm⊗Ip=Vk+1,Vk+2,…,Vm+1H¯m−k⊗Ip.SinceRk⊥FℛVk, that is, Vk+1⊥FVk, the following flexible global Arnoldi-like relation holds:(16)AZk,Zk+1,…,Zm=V1,V2,…,Vm+1Ik⊗IpVk,AZk+1,…,Zm⊗Ip0H¯m−k⊗Ip,AZm=V1,V2,…,Vm+1⋅IkVk,AZk+1,…,Zm0H¯m−k⊗Ip,AZm=VmH¯m⊗Ip,Vm+1,Vm+1F=Im+1.At the end, we find the approximate solutionXm=X0+Zmy∗⊗Ip∈ℂn×p by minimizing the norm R0−AZmy⊗IpF, where y∗ satisfies(17)y∗=argminy∈ℂmR0−AZmy⊗IpF=argminy∈ℂmβ0−H¯my2,with β0=Vm+1,R0F.To improve the convergence of iterative solvers, the smallest magnitude eigenvalues will be deflated. Following the idea of [12, 32], we use harmonic Ritz vectors with those eigenvalues to construct approximate invariant subspace and then the subspace for next cycles is retained. Firstly, we compute k harmonic Ritz vectors y˜k=Wmgk⊗Ip, such that(18)AZmWmHFWmgi⊗Ip−θiWmgi⊥FAZm,where Wm∈ℂn×mp, and each harmonic Ritz pair θi,gi satisfies(19)H¯mHH¯mgi=θiH¯mHVm+1HWmgi,1≤i≤k.It is noted that the harmonic residual vectors satisfy(20)AZmWmHFWmgi⊗Ip−θiWmgi∈ℛVm+1.Observe that the residualRm=R0−AZmy∗⊗Ip=Vm+1y∗⊗Ip∈ℛVm+1⊥FℛAZm. So, the residual and the harmonic residuals are in the same p-dimensional space ℛAZm⊥F∩ℛVm+1. Then we can characterize the relationship by the following formula:(21)AZmGk⊗Ip=WmGk⊗Ip,Rm⋅diagθ1,…,θka1×k⊗Ip,where a:,i satisfies AZmgi⊗Ip−θiVmgi⊗Ip=a1,iRm, 1≤i≤k, and Gk=g1,…,gk∈ℂm×k.Due to the relations (10) and (21) can be also expressed as(22)Vm+1H¯m⊗IpGk⊗Ip=Vm+1H¯mGk⊗Ip=WmGk⊗Ip,Rm⋅diagθ1,…,θka1×k⊗Ip.Next, we compute the reduced QR-factorization ofH¯mGk=QR with Q∈ℂm+1×k and R∈ℂk×k and then define Zk, Vk∈ℂn×k as(23)Zk=ZmGkR−1⊗Ip,Vk=Vm+1Q⊗Ip,such that AZk=Vk. Since Rm is F-orthogonal to the subspace ℛAZm and Vm+1 is F-orthonormal, we have(24)Vk,RmF=Vm+1Q⊗IpR⊗Ip,RmF=Vm+1H¯mGk⊗Ip,RmF=AZmGk⊗Ip,RmF=0,Vk,VkF=Vm+1Q⊗Ip,Vm+1Q⊗IpF=QHQ=Ik.Therefore, at the end of the loop, an updateZk is generated, scaled such that Vk=AZk has F-orthonormal columns; the previous global Krylov subspace basis is dropped, and then we restart. Details of the FGL-GCRO-DR method are summarized in Algorithm 2.Algorithm 1: The FGL-GCRO-DR method.
Input:A∈ℂn×n, B∈ℂn×p. Choose m, k, tol and initial guesses X0.Output:Xm such that RmF/R0F≤tol.(1)
R0=B−AX0(2)
WhileRmF>R0F×toldo(3)
Ifm==1(4)
γ=R0F, V1=R0/γ(5)
Apply Algorithm2 to generate Zm, H¯m, Vm+1 such that AZm=Vm+1H¯m⊗Ip, and setWm=Vm(6)
Else(7)
Applym−k steps of the Algorithm 2 with the starting rectangular matrix Vk+1=Rk/RkF, such that AZk+1,…,Zm−VkVk,AZk+1,…,Zm⊗Ip=Vm−k+1H¯m−k⊗Ip with Zj=ℳj−1Vj(8)
Endif(9)
y∗=argminy∈ℂmγ−H¯my2 with γ=Vm+1,R0F(10)
SetXm=X0+Zmy∗⊗Ip, Rm=B−AXm(11)
Compute thek smallest eigenpairs θi,gi of H¯mHH¯mgi=θiH¯mHVm+1HWmgi,1≤i≤k,(12)
QR=H¯mGk(13)
Vk=Vm+1Q(14)
Zk=ZmGkR−1(15)
Wk=WmGkR−1(16)
Wm=Wk,Vk:,k+1p:mp(17)
EndwhileAlgorithm 2: Flexible global Arnoldi process.
γ=R0F and set V1=R0/γ;forj=1,…,m, doZj=ℳjVjV^=AZjfori=1,…,j, dohij=Vi,V^F, i=1,…,j;V^=V^−Vihij;end forhj+1,j=V^F and set Vj+1=V^/hj+1,jend forIn addition, the major advantage of FGL-GCRO-DR is that it can retain important spectral information generated during the solution of thei-th linear system and then exploit such information to accelerate convergence of the iterations when solving the subsequent i+1-th system. Then we detail how to recycle subspace when solving sequences of slowly changing linear systems (8) as follows (Algorithm 3).Algorithm 3: FGL-GCRO-DR for solving sequences of slowly changing linear systems.
(1)
Suppose thatVki, Zki and Wki are defined from solving i-th linear system linear and satisfy AiZki=Vki with Vki,Vki=Ik.(2)
R0i+1=Bi+1−Ai+1X0i+1(3)
QR=Ai+1Zki(4)
Vki+1=Q(5)
Zki+1=ZkiR−1⊗Ip(6)
Wki+1=WkiR−1⊗Ip(7)
X0i+1=X0i+1+Zki+1Vki+1,R0i+1⊗Ip(8)
R0i+1=R0i+1−Vki+1Vki+1,R0i+1⊗Ip(9)
Applym−k steps of Algorithm 2 with the starting rectangular matrix Vk+1i+1=Rk/RkF, such that Ai+1Zm−ki+1−Vki+1Vki+1,Ai+1Zm−ki+1⊗Ip=Vm−k+1i+1H¯m−k⊗Ip with Zji+1=ℳj−1Vji+1(10)
y∗=argminy∈ℂmβ−H¯my2 with β=Vm+1i+1,R0i+1F(11)
Xmi+1=Xm−1i+1+Zmi+1y∗⊗Ip(12)
Rmi+1=Bi+1−Ai+1Xmi+1(13)
Wmi+1=Wki+1,Wmi+11:n,k+1p:mp


## 2.3. Comments on the FGL-GCRO-DR Method

To achieve better efficiency of the linear system solvers, we not only employ the subspace recycling at each cycle but also apply the flexible preconditioning operator at each iteration. As we know, the variable preconditioning can be construed as an inexact subspace method. If the matrixA is difficult to invert, we will resort to several steps of iterative method as flexible preconditioning, so we obtain in some way an approximate result,(25)Z^j≈A−1Vj,at j-th (inner) iteration for the equation AZj=Vj,j=1,…,m. In such a situation, what are the properties of the true residual obtained by using the flexible preconditioning? We will address this question in this section. Following [33, 34], we set(26)Pj=Vj−AZ^j,j=1,…,m,to be the residual that results after the iterative solver has been terminated; then(27)Z^j=A−1Vj−A−1Pj,j=1,…,m.Therefore, the perturbation in linear operation withA can be written explicitly as(28)AZ^j=AA−1Vj−AA−1Pj=I+PjVjHFVjF2Vj=I−EjVj.The perturbation implies that the error computed has the bound(29)Ej≤PjFVjF.In fact, the solution is obtained directly from the approximation subspaceℛZ^m; that is,(30)Xm=Z^my∗˜⊗Ip,indicating that the perturbation Ej does not affect the computation of the true residual. In other words, in exact arithmetic(31)Rm=R0−AZ^my˜∗⊗Ip=Vm+1βe1⊗Ip−Vm+1H¯my˜∗⊗Ip=Vm+1βe1−H¯my˜∗⊗Ip,and is orthogonal to ℛVm+1. As a consequence, the perturbation of the linear operations only influences the goodness of ℛZ^m and the accuracy of the approximate solution. Therefore, we deduce that, in the flexible setting, a dynamic stopping criterion or a much looser tolerance can be employed to save computation by limiting the flexible preconditioning work without sacrificing the closeness of the computed and the true residuals.


## 2.4. The Relationship between Two Flexible Global Methods with Subspace Recycling

As we know, FGCRO-DR [32] and FGMRES-DR [35] are two important methods in the general category of augmented (deflated) GMRES-type methods, which choose harmonic Ritz vectors as augmented (deflated) subspace. They emphasize on the solution of linear system Ax=b which has only one equation. The relationship between FGCRO-DR and FGMRES-DR has been well described by [32]. It is shown that they are algebraically equivalent as long as a collinearity condition is satisfied.In fact, FGL-GCRO-DR is a global generalization of FGCRO-DR. For the sake of descriptive integrality, we also present a global generalization of FGMRES-DR algorithm, referred to as FGL-GMRES-DR, in Algorithm 4. As for global versions, the FGL-GCRO-DR method is also closely related to the FGL-GMRES-DR approach in algebraical analysis, but their numerical behavior will be slightly different. Details will be shown in Section 3.Algorithm 4: FGL-GMRES-DR.
Input:A∈ℂn×n, B∈ℂn×p. Choose m, k, tol and initial guessesX0Output:Xm such that RmF/R0F≤tol(1)
R0=B−AX0, γ=R0F, V˜1=R0/γ(2)
Apply Algorithm2 to generate Z˜m, V˜m+1, H¯m such that AZ˜m=V˜m+1H¯m⊗Ip(3)
Solveminγ−H¯my2 for y∗. Set X0=X0+Z˜my∗⊗Im, R0=B−AX0(4)
whileRmF>R0F×toldo(5)
Compute thek smallest eigenpairs θi,gi of Hm+hm+1,mHhm+1,mHm−HememH.(6)
QR=Pk01×kγ−H¯my∗(7)
V˜k=V˜m+1Q(8)
Z˜k=Z˜mQ(9)
H¯k=QHH¯mQ1:m,1:k(10)
Applym−k steps of Algorithm 2 with the starting rectangular matrix V˜k+1=Rk/RkF, such that AZ˜m−k=V˜m−k+1H¯m−k⊗Ip with Z˜j=ℳj−1V˜j(11)
y∗=argminy∈ℂmγ−H¯my2 with γ=Vm+1,R0F(12)
Xm=Xm−1+Z˜my∗⊗Im(13)
Rm=B−AXm(14)
endwhileFurthermore, we compare the main calculation cost and storage capacity for each cycle of FGL-GCRO-DR and FGL-GMRES-DR in Table1. In practical applications, n≫m,k. Here, opM and opA denote floating-point calculation counts for the preconditioning application and the linear operations, respectively.Table 1
The main calculation cost and storage capacity of FGL-GCRO-DR and FGL-GMRES-DR.
Computation ofFGL-GMRES-DR (m, k)FGL-GCRO-DR (m, k)Vk+1:,1:k+1p2m+1k+1np2nm+1kp+npZk:,1:kp2nmkp2nmkpVm:,k+2p+1:m+1popAm−k+m−k2n+2k+5npopAm−k+m−k2n+2k+6npZm:,k+1p+1:mpopMm−kopMm−kTotal costC+m+k−1npCMemory requirement2m+3np3m+3npThe computational cost of harmonic Ritz vectors and “minγ−H¯my2” have not been induced in Table 1 since the cost is approximately Op3, much more less than n. C=opM+opAm−k+2m+1kp+p+2mkp+m−k2m+2k+6np is referred to as the total cost of FGL-GCRO-DR. One sees that the time complexities of two global algorithms are the same, that is, Om−kn2p, but FGL-GMRES-DR requires slightly more computation per cycle than FGL-GCRO-DR in terms of computational cost.


## 2.5. Application to the General Coupled Matrix Equations

FGL-GCRO-DR substantially belongs to the family of the global Krylov subspace method, which can be also applied to solve the general coupled matrix equations of the form(32)∑j=1lAi,1,jX1Bi,1,j+Ai,2,jX2Bi,2,j+⋯+Ai,l,jXlBi,l,j=Di,i=1,2,…,l,including the (coupled) Sylvester and the (coupled) Lyapunov matrix equations. The literature on solving (32) is particularly rich and still growing [36–53], due to the prominent role of such matrix equations in stability theory [54], perturbation analysis [55], control theory [56, 57], and so forth. In particular, the global subspace solvers perform better than other approaches for solving (32), especially for the coefficient matrices which are large and sparse; for more details, see [44–53, 58]. Therefore, we are interested in exploiting the FGL-GCRO-DR method for solving (32). To do so, the linear operator can be defined as follows:(33)ℱ:ℝm×n×⋯×ℝm×n⟶ℝm×pn,X=X1,X2,…,Xl⟶ℱX=A1X,A2X,…,AlX,where Ai=∑k=1l∑j=1lAi,k,jXkBi,k,j,i=1,2,…,l. Exploiting the operator ℱ, we can rewrite (32) as(34)ℱX=D,where D=D1,D2,…,Dl. The effectiveness of FGL-GCRO-DR for the solution of (32) is shown in Section 3.


## 3. Numerical Experiments

This section presents some numerical examples to illustrate the effectiveness of our algorithm for solving shifted linear system (1) and matrix equations (8) and compare the performance with the global GMRES m approach (referred to as GL-GMRES m) [31, 50, 59] and GL-GMRES-DR [60], with or without preconditioning. In each example, the flexible preconditioning is dominated by 5 steps of the global full GMRES [50].We consider three numerical illustrations here. The first two experiments are sequences of linear systems arising from lattice QCD [2–4, 61] and the finite difference scheme of the convection-diffusion-reaction equation [62], respectively. The last one is the general coupled Sylvester matrix equations. The comparisons are mainly made in two aspects: the execution time in seconds (T (s), for short) and number of matrix-vector products (mvps, for short).All the numerical experiments are carried out on MATLAB 2017b. We chooseX0=0n×p, m=20, and k=10. The stopping criteria are Bi−AiXmiF/BiF≤10−8,i=1,2,…,s for all the approaches, or when mvps exceed the set maximum of matrix-vector multiplications (MAXIT, for short). Moreover, the symbol † represents that the solver does not achieve the stopping criteria within MAXIT and RHSs=randn,p.Example 1.
We first consider two groups of complex Lattice QCD matrices from Matrix Market (please refer tohttp://math.nist.gov/Matrix-Market for details). One group has seven 3172×3172 matrices and the other has seven 49152×49152 matrices. For each matrix D, a critical parameter kc is provided such that I−k˜D with 0≤k˜≤kc is positive-real. For each set, we take Ai=kc+10−3I−Di as the sequence matrix, and Ti=0.01,0.02,0.03,1,2,3,4. The numerical computation is carried out with MAXIT=1000. The corresponding numerical results are reported in Tables 2 and 3. In addition, we also depict the residual histories for linear systems solved using iterative algorithms with or without flexible preconditioning.
As seen from Tables2 and 3 and Figures 1 and 2 , it is found that using deflation and flexible preconditioning is conducive to improving the convergence rate. FGL-GCRO-DR has a significant decrease in terms of mvps and T(s), which are similar to FGL-GMRES-DR, but performs slightly better. The GL-GMRES method is less effective and does not converge in some test problems. The above results show that FGL-GCRO-DR has the potential to improve the convergence, compared with other global Krylov subspace solvers.Table 2
Performance comparisons for3072×3072 QCD matrices.
ProblemsMethodGL-GCRO-DRFGL-GCRO-DRGL-GMRESmFGL-GMRESmGL-GMRES-DRFGL-GMRES-DRMatrix 1mvps780135940180†150T (s)6.4312.1664.1663.426†2.993Matrix 2mvps4318348010049090T (s)3.3281.6182.1291.7543.4641.741Matrix 3mvps50199800140710110T (s)4.4622.1083.7332.5234.0602.329Matrix 4mvps541100620140530110T (s)4.0902.1152.7222.5513.5252.127Matrix 5mvps265463406027060T (s)1.5731.0051.5601.0401.6181.223Matrix 6mvps803148†240820170T (s)5.0333.463†4.4145.4683.521Matrix 7mvps258503006026060T (s)2.0241.0041.3351.0641.7061.138Table 3
Performance comparisons for49152×49152 QCD matrices.
ProblemsMethodGL-GCRO-DRFGL-GCRO-DRGL-GMRESmFGL-GMRESmGL-GMRES-DRFGL-GMRES-DRMatrix 1mvps856162†220†180T (s)1.306e + 028.009e + 01†9.015e + 01†8.440e + 01Matrix 2mvps685139860180720150T (s)8.797e + 016.363e + 018.556e + 017.391e + 019.145e + 016.868e + 01Matrix 3mvps706131760180650140T (s)9.972e + 016.186e + 017.722e + 017.368e + 019.033e + 016.373e + 01Matrix 4mvps695150740180720160T (s)7.931e + 017.016e + 017.570e + 017.427e + 018.130e + 017.271e + 01Matrix 5mvps798110†140800120T (s)9.341e + 015.027e + 01†5.734e + 019.255e + 016.274e + 01Matrix 6mvps838115†180870140T (s)8.668e + 015.480e + 01†7.211e + 019.357e + 016.274e + 01Matrix 7mvps976115†160980130T (s)1.104e + 025.082e + 01†6.585e + 011.203e + 025.845e + 01Figure 1
The “.−” line and the “+−” line denote the convergence curves for the global method and the flexible global method for solving the s linear equations with n=3072 in sequence, respectively.Figure 2
The “.−” line and the “+−” line denote the convergence curves for the global method and the flexible global method for solving the s linear equations with n=49152 in sequence, respectively.Example 2.
For the second case, we assess the performance of global solvers on three-dimensional separable model convection-diffusion problem in the unit cube with the homogeneous Dirichlet boundary conditions [62]:(35)−ϵ△u+vT∇u−ru=f,inΩ,where ϵ=1 and r=0,100,400,600,800,1000. For obtaining the sequence matrices, we set v=0.1∗k,0.5∗i∗k,k with k=1,2,…,5 and i=−1 and then use the central finite difference scheme to discrete the convection-diffusion equation with the mesh size h=0.025 or h=0.02. Thus, two groups of complex matrices are obtained. One group has five 59319×59319 matrices and the other has five 117649×117649 matrices. As the definition of MAXIT, we choose MAXIT=500.
The numerical results reported in Tables4 and 5 and Figures 3 and 4 show that FGL-GMRES is slightly superior to FGL-GCRO-DR when the mesh size h=0.025. However, as the number of grid points increases, FGL-GCRO-DR enjoys the significant decrease in the number of mvps and T(s). Meanwhile, according to Tables 4 and 5, we can see that the convergence is speeded up remarkably by considering the flexible preconditioning. Moreover, even the size of grid points increases, and the average number of mvps for FGL-GCRO-DR does not increase.Table 4
Performance comparisons for 59319 convection-diffusion problem.
ProblemsMethodGL-GCRO-DRFGL-GCRO-DRGL-GMRESmFGL-GMRESmGL-GMRES-DRFGL-GMRES-DRMatrix 1mvps165293404017040T (s)2.781e + 011.735e + 013.693e + 011.693e + 012.811e + 011.826e + 01Matrix 2mvps172323404019040T (s)2.655e + 011.766e + 013.664e + 011.686e + 013.115e + 011.824e + 01Matrix 3mvps178313404018040T (s)2.907e + 011.713e + 013.718e + 011.688e + 012.946e + 011.828e + 01Matrix 4mvps182343404019040T (s)3.095e + 011.818e + 013.646e + 011.680e + 013.120e + 011.826e + 01Matrix 5mvps179323404020040T (s)2.747e + 011.747e + 013.658e + 011.675e + 013.271e + 011.822e + 01Table 5
Performance comparisons for 117649 convection-diffusion problem.
ProblemsMethodGL-GCRO-DRFGL-GCRO-DRGL-GMRESmFGL-GMRESmGL-GMRES-DRFGL-GMRES-DRMatrix 1mvps209354806025040T (s)7.486e + 013.630e + 011.018e + 024.998e + 018.223e + 013.604e + 01Matrix 2mvps21836†6024050T (s)7.660e + 014.006e + 01†4.988e + 017.876e + 014.536e + 01Matrix 3mvps21942†6024050T (s)7.665e + 014.490e + 01†5.013e + 018.117e + 014.541e + 01Matrix 4mvps242414806027050T (s)7.664e + 014.396e + 011.026e + 024.977e + 018.795e + 014.545e + 01Matrix 5mvps264444606023050T (s)8.102e + 014.351e + 019.859e + 014.985e + 017.514e + 014.534e + 01Figure 3
The “−” line and the “+−” line denote the convergence curves for the global method and the flexible global method for solving the s linear equations with n=59319 in sequence, respectively.Figure 4
The “−” line and the “+−” line denote the convergence curves for the global method and the flexible global method for solving the s linear equations with n=117649 in sequence, respectively.Example 3.
For the last test case, we are interested in solving the general coupled Sylvester matrix equations [45, 47](36)AX−YB=E,CX−YD=F,where(37)A=diagdiagrandn+2+triurandn,n,1,B=diagdiagrandn+1.75+triurandn,n,1,C=−diagdiagrandn+2+triurandn,n,1,D=diagdiagrandn+1.5+triurandn,n,n,E=randnn,n,F=randnn,n.
Here, we considern=50,60,70,80,90,100 different orders. As the definition of MAXIT, we choose MAXIT=500. The corresponding numerical results of mvps and T(s) are reported in Table 6 and the residual histories are depicted in Figures 5 and 6.
As seen from Table5 and Figures 5 and 6, it is noted that the flexible deflated GMRES-type methods (FGL-GCRO-DR and FGL-GMRES-DR) are found to be efficient in terms of performance test problems of different orders, while GL-GMRES m and FGL-GMRES m fail to converge in the given steps with or without variable preconditioning, while the parameter n is greater than 80. The convergence is improved remarkably by considering the eigenvalue deflation technique and the flexible preconditioning. One should note that FGL-GCRO-DR enjoys the significant decrease in the number of mvps and T(s) and performs slightly better than FGL-GMRES-DR method. In addition, when n increases, the average number of mvps and T(s) for FGL-GCRO-DR do not increase.Table 6
Performance comparisons for the general coupled Sylvester matrix equations.
ProblemsMethodGL-GCRO-DRFGL-GCRO-DRGL-GMRES (m)FGL-GMRES (m)GL-GMRES-DRFGL-GMRES-DRn  = 50mvps134372393713029T (s)4.060e − 011.336e − 012.988e − 011.396e − 019.266e − 013.421e − 01n  = 60mvps21741†5917346T (s)6.551e − 012.165e − 01†2.195e − 011.668e + 007.521e − 01n  = 70mvps24353†11523261T (s)1.042e + 003.929e − 01†6.168e − 013.112e + 001.354e + 00n  = 80mvps507100†507334138T (s)2.661e + 008.868e − 01†2.661e + 006.143e + 003.894e + 00n  = 90mvps507100††503175T (s)3.294e + 001.109e + 00††1.206e + 016.228e + 00n  = 100mvps†166†††228T (s)†2.326e + 00†††9.962e + 00Figure 5
The convergence curves for the global methods and the flexible global methods for solving general coupled Sylvester matrix equations.Figure 6
The convergence curves for the global methods and the flexible global methods for solving general coupled Sylvester matrix equations.


## 4. Conclusions

In this study, a new flexible global GCRO-DR method is proposed for addressing sequences of shifted linear systems and general coupled matrix equations. The proposed method does not require the right-hand sides to be related and can also be compatible with general preconditioning. By performing an error analysis, we deduce that a much looser tolerance can be applied to save computation by limiting the flexible preconditioned work without sacrificing the closeness of the computed and the true residuals. Meanwhile, we also propose a global version of the FGMRES-DR method and then analyze the relationship between two flexible global methods. Finally, numerical experiments show that FGL-GCRO-DR significantly reduces the mvps, saves the runtime, and enjoys faster convergence than some other global GMRES-type approaches on tough problems.As future work, we plan to investigate the numerical properties of FGL-GCRO-DR on large-scale realistic problems. Of interest are applications related to, for example, model adaptation in Gaussian process models [63] or stochastic finite element methods [64] in three dimensions where variable preconditioning using inexact solvers has to be usually considered. Another extension of this work is a distributed-memory parallel implementation and testing using MPI + OpenMP and MPI + CUDA. Investigation will focus on implementing, testing, and comparing the runtime of the introduced FGL-GCRO-DR on CPUs and GPUs with respect to other existent similar methods.



---

*Source: 5589582-2021-05-12.xml*