Table of Contents Author Guidelines Submit a Manuscript
Advances in Mathematical Physics
Volume 2018, Article ID 1369707, 12 pages
https://doi.org/10.1155/2018/1369707
Research Article

Residual-Based Simpler Block GMRES for Nonsymmetric Linear Systems with Multiple Right-Hand Sides

1School of Mathematics and Computational Science, Xiangtan University, Xiangtan 411105, China
2College of Science, Hunan University of Science and Engineering, Yongzhou 425100, China
3Department of Mathematics, East China University of Science and Technology, Shanghai 200237, China

Correspondence should be addressed to Qinghua Wu; moc.361@hqwkcaj

Received 16 November 2017; Accepted 12 March 2018; Published 16 April 2018

Academic Editor: Soheil Salahshour

Copyright © 2018 Qinghua Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

We propose in this paper a residual-based simpler block GMRES method for solving a system of linear algebraic equations with multiple right-hand sides. We show that this method is mathematically equivalent to the block GMRES method and thus equivalent to the simpler block GMRES method. Moreover, it is shown that the residual-based method is numerically more stable than the simpler block GMRES method. Based on the deflation strategy proposed by Calandra et al. (2013), we derive a deflation strategy to detect the possible linear dependence of the residuals and a near rank deficiency occurring in the block Arnoldi procedure. Numerical experiments are conducted to illustrate the performance of the new method.

1. Introduction

In this paper, we consider iterative methods for solving a system of linear algebraic equations:where is a nonsingular matrix of order and and are rectangular matrices of dimension with . For solving such systems, the block GMRES [1] and its variants are very popular. Block GMRES is based on the block Arnoldi process and is formally fully analogous to the ordinary GMRES algorithm by Saad and Schultz [2].

The following notation is used throughout the paper. Subscripts denote the iteration index and superscripts distinguish between individual columns in a block. We denote by the Euclidean vector norm and the induced matrix norm and by the Frobenius norm. Moreover, for of rank , is the spectral condition number, where are the extremal singular values of .

Given an initial approximation to the solution of (1), letand then in analogy to the unblocked case, we build a sequence of iterates such thatwhere . Equation (3) is equivalent to minimizing every column of , that is,and also to the orthogonality conditionwhere is the orthogonality relation induced by the Euclidean inner product. Assume that is of the form , where is a basis of . Then, we obtain the th residual matrix

Central to the usual implementations of block GMRES is the block Arnoldi process [1], which can be used to construct the orthonormal basis of . In practice, the possible linear dependence of the residuals of the systems requires an explicit reduction of the number of right-hand sides. In [3], this was called deflation. If the block residual is nearly rank deficient, block GMRES should be implemented with deflation and there are various sophisticated rank-revealing QR factorizations. For details, see [3] and the references therein. We can write the nondeflated block Arnoldi process as shown in Algorithm 1.

Algorithm 1: Nondeflated block Arnoldi process.

From Algorithm 1, we obtain formally the ordinary Arnoldi relationwhere the matrix isIn the block Arnoldi algorithm, holds due to the QR factorizations, and when , where is a unit matrix and 0 is a zero matrix of order . This indicates that the whole process is equivalent to the one in which the block vectors are generated column by column using an ordinary modified Gram-Schmidt process.

From (6) and (7), we obtain the fundamental block GMRES relationwhere is the first column of the unit matrix (the size changes with ), is an upper triangular matrix obtained in Arnoldi’s initialization step, and is the “block coordinates” of with respect to the block Arnoldi basis.

Using (9), the least squares problem (3) is solved by recursive QR factorization of , updated by applying Givens rotations. Once the norm of the residual is small enough, the triangular system with the computed -factor is solved, and the approximate solution is computed. The detailed algorithm of block GMRES can be found in [35].

The block GMRES method with deflation at each iteration was proposed in [6]. And a deflation strategy was investigated to detect when a linear combination of approximate solutions is already known; for details, see [7]. In this paper, we deal with a different approach and compare the situation with deflation and without deflation. Let be a block basis of . Instead of building a block orthogonal basis of , we look for a block orthogonal basis of . As a special case, simpler block GMRES (SBGMRES) was proposed by Liu and Zhong [8], where . In this paper, we will consider the basis . We call this case the residual-based simpler block GMRES (RB-SBGMRES).

The paper is organized as follows. In Section 2, the RB-SBGMRES and SBGMRES algorithms are described (without deflation). In Section 3, some comparison between RB-SBGMRES and SBGMRES is established. In Section 4, the RB-SBGMRES method with deflation and the corresponding algorithm RB-SBGMRES-D (Table 2) are derived. In Section 5, these algorithms are compared using test matrices taken from the Matrix Market [9]. Conclusions are included in the last section.

2. Residual-Based Simpler Block GMRES

Suppose that is a basis of . The orthogonal basis of is thus obtained from the QR factorization of ; that is,where is an upper triangular matrix with order .

Due to the orthogonality property , the th residual matrix can be computed recursively aswhere .

Define . Since the columns of form a basis of , we can represent in the formwhere is the “block coordinates” of with respect to the block basis . Due to , , and , it follows thatHence, once the residual norm is small enough, we can solve this upper triangular system (13) and then compute the approximate solution .

We now present the RB-SBGMRES method without deflation as shown in Algorithms 2 and 3 and Table 1. For comparison, we also present the SBGMRES method proposed in [8].

Table 1: Computational cost of a cycle of RB-SBGMRES (or SBGMRES).
Table 2: Computational cost of a cycle of RB-SBGMRES-D (Algorithm 4).
Algorithm 2: RB-SBGMRES(m).
Algorithm 3: SBGMRES(m).

3. Comparison with SBGMRES

In [8], an equivalence between SBGMRES and classical block GMRES had been established. Algorithm 2 indicates that RB-SBGMRES is equivalent to block GMRES; that is, search a solution , such that .

On the other hand, for single right-hand side, it has been observed in [10] that the gap between the true residual and the updated residual can be strongly influenced by the conditioning of , which is the basis of , and the choice of the basis has an effect on the conditioning of the matrix [10]. Since we compute the coordinates of the correction in the basis by (13), the approximate solution becomes inaccurate as the conditioning of grows. Simpler GMRES [11] is, in general, less accurate than GMRES and is inherently unstable due to the choice of the basis . It is easy to formulate an analogous conclusion in the block case.

There is a theorem about condition number of in [8] stated as follows.

Theorem 1. For , one haswhere is the number of right-hand sides.

The condition number of may have an effect on the conditioning of the matrix . In the following, we formulate an analogous theorem on the condition number of .

Theorem 2. Suppose that steps of RB-SBGMRES have been taken and , , if and only if is a zero vector; then, one haswhere and is the number of right-hand sides.

Proof. Let be a QR factorization of . We have whereMoreover, we get .
From (11), we obtainThen, it follows thatSince the columns of are orthogonal, it follows thatOn the other hand,Thus, using the triangle inequality and the fact that the 1-norm is bounded from above by the 2-norm, we haveThe norm of can be expressed using (11) asWith , we haveThe proof then follows from .

Theorem 1 indicates that the conditioning of is inversely proportional to the actual relative norm of the residual. Once residuals become small, this will lead to the ill-conditioning of the matrices and the matrices , and SBGMRES will behave unstably after some initial residual reduction. However, from Theorem 2, the conditioning of is related to the intermediate decrease of the residual norms, not to the residual decrease with respect to the initial residual. For single right-hand sides, it has been observed that remains well conditioned while becomes ill-conditioned [10].

4. RB-SBGMRES with Deflation

When block Krylov subspace methods are used for the solution of linear systems of equations with multiple right-hand sides, the linear dependence of the residual of the systems may occur, and this is called deflation. Deflation may be possible at startup or in a later step. Sometimes, we need to incorporate a strategy for detecting when a linear combination of systems has approximately converged. Recently, Calandra et al. derived a deflation strategy to detect a near rank deficiency occurring in the block Arnoldi procedure in [7]. We provide a brief overview of the method.

Assume that the QR factorization of has been performed aswith having orthonormal columns and . To circumvent deflation at startup, the subspace decomposition at the beginning of the cycle is derived by finding a unitary matrix , such thatwith , , and .

The unitary matrix is determined by the singular value decomposition (SVD) of and set . Choose a relative deflation threshold and select singular values of such that .

Define and , for with orthonormal columns, and assume that the following block Arnoldi relation holds at the beginning of the th iteration:with , , , and . The th iteration of the deflated block Arnoldi procedure produces matrices and which satisfywith having the following block structure:

Defining , (28) can be reformulated asThe subspace decomposition is performed by finding a unitary matrix such thatHence, we obtain Defining as , this leads to which is the block Arnoldi relation required at the beginning of the next iteration.

From (31), the unitary matrix has the following matrix structure:Defining as and with , the unitary matrix is determined by the following steps (for details, see [7]): (1)SVD of , with .(2)QR decomposition of matrix .

In the case of simpler block GMRES method, the relationship , which is an important ingredient for the block GMRES method, cannot be established. We must find another strategy to perform the subspace decomposition.

Assume that the QR factorization of has been performed as

We compute analogously as in (26), which leads to the following formulas:with , and , . To build a block orthogonal basis of , we compute the QR factorization of :with , , , , , and , and is computed similarly to .

Assume that the block Arnoldi relationholds at the beginning of the th iteration of the deflated block Arnoldi procedure, with , , , , , and , where , , and are defined as follows: , , , and .

Set , with , and define as . We generate by calculating the QR factorization of , which leads to the following formulas:with , , , , and , and is a unitary matrix, which is determined similarly to .

Compute and make it orthogonal to the columns of by an ordinary modified Gram-Schmidt process, such that with , and have the following block structure:where , , , and . Define as

We formulateby calculating the SVD of , instead of computing the SVD of . Note that ; we set and select singular values of such that . Once a near rank deficiency occurs in the block Arnoldi procedure, it will be reflected by the singular values of . We now present the RB-SBGMRES method with deflation as shown in Algorithms 4 and 5.

Algorithm 4: RB-SBGMRES-D(m).
Algorithm 5: Determination of (, , ) or (, , ).

Denote by and the number of Krylov directions keeping at the th iteration of the th cycle; Algorithm 4 indicates that and .

5. Numerical Experiments

In this section, RB-SBGMRES is tested and compared with SBGMRES. The test matrices were taken from the Matrix Market [9]. All computations were carried out using Matlab on a PC with the usual double precision, where the floating point relative accuracy is . In the following examples, we take as the zero matrix; thus, the initial residual matrix is , and we set , where is the order of the matrix . We plot the relative true norm of residualand condition number of , respectively, in Figures 1, 2, and 3. The condition number of is computed by Matlab internal function cond.

Figure 1: Example 1, the relative true residual norms versus the number of iterations.
Figure 2: Example 1, the condition numbers versus the number of iterations.
Figure 3: Example 1, the relative true residual norms versus the number of iterations.

Example 1. The matrix is fs1836 (, ) which is a matrix of order , and we have right-hand sides which are chosen as,, where is the th column of matrix , , and .

In Figures 1 and 2, we show the relative true norm of residual and condition number of for the SBGMRES and RB-SBGMRES, respectively. It is clearly seen from Figure 1 that the relative true norm of residual of SBGMRES may stagnate at a significantly higher level than that of RB-SBGMRES. The reason for this is that the condition number of the matrix of SBGMRES increases significantly faster than that of RB-SBGMRES as Figure 2 shows. In Figure 3, we show the relative true norm of residual for the RB-SBGMRES-D and the BFGMRES-S proposed in [7], with . It is clearly seen from the figure that RB-SBGMRES-D can compete with BFGMRES-S.

Example 2. The matrix is steam1 (, ) which is a matrix of order , and we have right-hand sides which are chosen as,,,,, where is the th column of matrix , , and is the order of matrix .

It is clear from Figure 4 that the condition number of the matrix of SBGMRES increases faster than that of RB-SBGMRES. Since the number of right-hand sides , the maximum iteration number is . From Figure 5, we can observe that the numerical performance of RB-SBGMRES is better than that of SBGMRES.

Figure 4: Example 2, the relative true residual norms versus the number of iterations.
Figure 5: Example 2, the condition numbers versus the number of iterations.

Figure 6 shows the relative true norm of residual for the RB-SBGMRES-D and BFGMRES-S, with . It is clearly seen that the performances of RB-SBGMRES-D and BFGMRES-S are almost same.

Figure 6: Example 2, the relative true residual norms versus the number of iterations.

Example 3. The matrix is add20 (, ) which is a matrix of order 2395, and we have right-hand sides which are chosen as ,,  for .

It is clear from Figure 7 that the numerical performance of RB-SBGMRES is better than that of SBGMRES, and from Figure 8, we see that the condition number of the matrix of SBGMRES increases fast while that of RB-SBGMRES remains at a significantly low level. Figure 9 shows the relative true norm of residual for the RB-SBGMRES-D and BFGMRES-S, with . The performances of RB-SBGMRES-D and BFGMRES-S are almost the same for Example 3.

Figure 7: Example 3, the relative true residual norms versus the number of iterations.
Figure 8: Example 3, the condition numbers versus the number of iterations.
Figure 9: Example 3, the relative true residual norms versus the number of iterations.

Example 4. The matrix is fs7601 (, ) which is a matrix of order 760, and we have right-hand sides which are chosen as,,,  for .

It is also obvious from Figure 10 that the RB-SBGMRES method is slightly more accurate than SBGMRES in this example, and Figure 11 also shows that the condition number of the matrix of SBGMRES increases faster than that of RB-SBGMRES. It is clear from Figure 12 that the performances of RB-SBGMRES-D and BFGMRES-S are almost the same for Example 4.

Figure 10: Example 4, the relative true residual norms versus the number of iterations.
Figure 11: Example 4, the condition numbers versus the number of iterations.
Figure 12: Example 4, the relative true residual norms versus the number of iterations.

In order to further verify that the condition number of the matrix of SBGMRES increases significantly faster than that of RB-SBGMRES, we compared a broad selection of different matrices from the Matrix Market and presented a comparison of the overall condition number trend. We select matrices, randomly, and set the same convergence threshold for two methods. We compare the iterations required to converge for two methods. It is easy to see from Table 3 that the condition number of the matrix of RB-SBGMRES is significantly smaller than that of SBGMRES and the number of iterations for the SBGMRES is slightly larger than that of RB-SBGMRES for most matrices.

Table 3

6. Conclusion

In this paper, we have proposed a minimum residual method mathematically equivalent to the block GMRES method for solving systems of linear equations with multiple right-hand sides. Numerical experiments show that, after some initial reduction, the relative true norms of residual SBGMRES may stagnate at a significantly higher level than that of RB-SBGMRES. This difference is clearly caused by the choice of the basis , which has an effect on the condition number of the matrix . Numerical experiments indicate that of RB-SBGMRES remains better-conditioned than of simpler block GMRES, which may become a very ill-conditioned triangular matrix. Since the coordinates of the correction in the basis are computed from (13), its error starts to diverge as grows and will become inaccurate. We see that the choice has a better numerical performance. In comparison with the case of deflation, we consider a deflation strategy to detect the possible linear dependence of the residuals of the systems and a near rank deficiency occurring in the block Arnoldi procedure for RB-SBGMRES method, which was later called RB-SBGMRES-D. Numerical experiments show that the performances of RB-SBGMRES-D and BFGMRES-S are almost the same.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Acknowledgments

This work is supported by the National Natural Science Foundation of China under Grants 11701170 and 11501193, the Natural Science Foundation of Hunan Province under Grants 2017JJ3092 and 2017JJ2102, the Key Program of the Scientific Research Foundation from Education Bureau of Hunan Province under Grant 09A033, the Scientific Research Foundation of Education Bureau of Hunan Province for Outstanding Young Scholars in University under Grant 10B038, and the Young Core Teacher Foundation of Hunan Province in University.

References

  1. Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, PA, USA, 2nd edition, 2003. View at Publisher · View at Google Scholar · View at MathSciNet
  2. Y. Saad and M. H. Schultz, “GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM Journal on Scientific and Statistical Computing, vol. 7, no. 3, pp. 856–869, 1986. View at Publisher · View at Google Scholar · View at MathSciNet
  3. M. H. Gutknecht, “Block krylov space methods for linear systems with multiple right-hand sides: an introduction,” in Modern Mathematical Models, Methods and Algorithms for Real World Systems, A. H. Siddiqi, I. S. Duff, and O. Christensen, Eds., pp. 420–447, Anamaya Publishers, New Delhi, India, 2007. View at Google Scholar
  4. V. Simoncini and E. Gallopoulos, “An iterative method for nonsymmetric systems with multiple right-hand sides,” SIAM Journal on Scientific Computing, vol. 16, no. 4, pp. 917–933, 1995. View at Publisher · View at Google Scholar · View at MathSciNet
  5. V. Simoncini and E. Gallopoulos, “Convergence properties of block GMRES and matrix polynomials,” Linear Algebra and its Applications, vol. 247, pp. 97–119, 1996. View at Publisher · View at Google Scholar · View at MathSciNet
  6. M. Robbé and M. Sadkane, “Exact and inexact breakdowns in the block GMRES method,” Linear Algebra and Its Applications, vol. 419, no. 1, pp. 265–285, 2006. View at Publisher · View at Google Scholar · View at MathSciNet
  7. H. Calandra, S. Gratton, R. Lago, X. Vasseur, and L. M. Carvalho, “A modified block flexible GMRES method with deflation at each iteration for the solution of non-Hermitian linear systems with multiple right-hand sides,” SIAM Journal on Scientific Computing, vol. 35, no. 5, pp. S345–S367, 2013. View at Publisher · View at Google Scholar · View at MathSciNet
  8. H. Liu and B. Zhong, “Simpler block GMRES for nonsymmetric systems with multiple right-hand sides,” Electronic Transactions on Numerical Analysis, vol. 30, pp. 1–9, 2008. View at Google Scholar · View at MathSciNet
  9. Matrix Market, http://math.nist.gov/matrixmarket.
  10. P. Jiránek, M. Rozložnk, and M. Gutknecht, “How to make simpler GMRES and GCR more stable,” SIAM Journal on Matrix Analysis and Applications, vol. 30, no. 4, pp. 1483–1499, 2008. View at Publisher · View at Google Scholar · View at MathSciNet
  11. H. F. Walker and L. Zhou, “A simpler GMRES,” Numerical Linear Algebra with Applications, vol. 1, no. 6, pp. 571–581, 1994. View at Publisher · View at Google Scholar · View at MathSciNet