Consider the large Sylvester matrix equation
where , and , with . If we define the operator as
Then the Sylvester matrix equation (1) can be written as
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
where denotes the Kronecker product operator, and
denotes the vectorize operator defined as (in MATLAB notation)
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 , 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 . For example, Mohseni Moghadam et al.  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 , 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 . 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 
and the associated -norm of is defined as
For two matrices , the -inner product is defined as 
where denotes the trace of a matrix, and represents the transpose of the matrix . It can be verified that 
Also, the -norm associated with this inner product is
Next we introduce a useful product that will be used latter:  Let and , where . Then elements of the matrix is defined as
It was shown that . 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 .
The -step weighted global Arnoldi process
Input: , , , a positive diagonal matrix and an integer number .
Compute . If then stop.
The weighted global Arnoldi process constructs a -orthonormal basis , i.e.,
for the matrix Krylov subspace
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
where . Define
then (6) can be rewritten as
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
where . The corresponding residual is
where is the first canonical basis vector in . Note that the residual is -orthogonal to , i.e.,
where , and “” means orthgonal with respect to the “” inner product.
where we used . Thus,
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.
A restarted weighted global GMRES algorithm for large Sylvester matrix equations (W-GLGMRES)
Input: , , . Choose the initial guess, , an initial positive diagonal matrix D and an integer , and the tolerance .
Output: The approximation .
Compute and ,
Run Algorithm 1 with the initial block vector to obtain the matrices .
Solving for .
Compute and . If , then stop, else continue.
Set and update the weighting matrix according to Options 1–3, and go to Step 3.
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
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
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
If is nonsingular, (19) is equivalent to
Then we define the “weighted harmonic Ritz vector” as , and the corresponding harmonic residual is
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
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
so we have . Let and , we obtain from (21) that
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 Algorithm3. Then we orthonormalize against to get , and let .
As a result,
where we used and are collinear with each other; see Theorem 3. Therefore,
Define , where , so we have
If we denote
then we have from (23) that
and there is a matrix such that
The condition yields
where is generally a dense matrix. In conclusion, we obtain
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
That is, , and a global Arnoldi-like relation still holds after restarting.
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
We are ready to present the main algorithm of this paper.
A weighted global GMRES with deflation for large Sylvester matrix equations (W-GLGMRES-D)
Input: . Choose an initial guess a positive diagonal matrix , the positive integer number and a convergence tolerance .
Output: The approximation .
Compute , and
Run the weighted global Arnoldi algorithm to obtain and .
Solve for .
Compute , and .
If , then stop, else continue.
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.
Extend the vectors to length with zero entries, then orthonormalize the vector against the columns of to form Set .
Compute and , note that is generally a full matrix. Update according to Options 1–3.
Run the weighted global Arnoldi algorithm from index to to obtain and , where the th block is the last columns of