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.

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 [3–5].

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

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 .