An Enhancement of the Accuracy of the BiCGStab Method for Solving Linear Systems with Single or Multiple Right-Hand Sides

Bouyghf, F.

doi:https://doi.org/10.1155/2023/8078760

Computational and Mathematical Methods

On this page

Abstract Introduction Conclusion Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2023 | Article ID 8078760 | https://doi.org/10.1155/2023/8078760

An Enhancement of the Accuracy of the BiCGStab Method for Solving Linear Systems with Single or Multiple Right-Hand Sides

F. Bouyghf^1,2

Academic Editor: Pedro Alonso

Received04 May 2023

Revised17 Aug 2023

Accepted23 Aug 2023

Published13 Sept 2023

Abstract

In this paper, we present a technique to improve the convergence of the biconjugate gradient stabilized (BiCGStab) method. This method was developed by Van der Vorst for solving nonsymmetric linear systems with a single right-hand side. The global and block versions of the BiCGStab method have been proposed for solving nonsymmetric linear systems with multiple right-hand sides. Using orthogonal projectors to minimize the residual norm in each step, we get an enhancement of the convergence of each version of the BiCGStab method. The considered methods are BiCGStab, global BiCGStab, and block BiCGStab methods, noted, respectively, as Gl-BiCGStab and Bl-BiCGStab. To show the performance of our enhanced algorithms, we compare them with the standard, global, and block versions of the well-known generalized minimal residual method (GMRES).

1. Introduction

The aim of the BiCGStab method studied in this paper is to solve the following nonsymmetric linear system where is a nonsingular matrix in and the vectors and are in . Problems such as (1) occur in most applications of scientific computing, engineering applications, and Navier-Stokes equations in computational fluid dynamics and structural mechanic computations based on finite element analysis. If the order of the matrix is small, we can solve (1) using direct methods, but if is large, direct methods can be prohibitively expensive both in terms of memory and time. So, iterative methods become appealing. These methods include the Krylov subspace methods. It is true that current methods give us interesting computational performance, but it is important to realize that there is no single method that can solve every linear system with desirable accuracy. Because it all depends on the conditioning number of the problem matrix and some characteristics of the matrix, i.e., its symmetric and whether it is positive definite or not. Many experiments are given in the recent book [1] to show this remark. However, a number of factors can influence the choice of the method, including the conditioning of the matrix and the number of nonzero values. In literature, there are many problems of linear systems with a sparse matrix which are obtained from real applications (heat transfer, fluid flow, mass transport, etc.) by using a numerical strategy for solving partial differential equations (finite difference, finite volume, and finite element methods) [2, 3]. The difference between all the existing methods for solving linear systems is in the level of accuracy, turnaround time, and storage. For system (1) with a symmetric positive definite matrix, we use the conjugate gradient method (CG). If is a symmetric matrix, we usually use the MINRES. Finally, the nonsymmetric case can be solved using GMRES, developed as the most popular and the most optimal in terms of precision but suffers from storage problems. The nonsymmetric case can also be solved using some short-recurrence Krylov subspace methods. We find for example the BiCG, CGS, and BiCGStab methods. These methods are derived from the extension of the CG in the nonsymmetric case. In this work, we focus on the nonsymmetric case, especially on the BiCGStab method as the most stable method compared with the BiCG and CGS methods. For more details about all the Krylov subspace methods, see [4] and the current book [1].

The BiCGStab method has been developed in [5] for solving (1). The global version of the BiCGStab method has been given in [6], and the block version of this method has been given in [7] for solving linear systems with several right-hand sides. Using orthogonal projectors and a technique given in [8, 9] to improve the convergence of the IDR method [10, 11], and some Krylov methods (BiCG and CGS), we give an enhancement of the convergence of the BiCGStab method.

The rest of this paper is organized as follows: In the next section, we recall the algorithm of the BiCGStab method [5]. Then, we propose an improvement in the convergence of this algorithm using orthogonal projectors. A partial and full improvement of the BiCGStab method is proposed and will be called PEnha-BiCGStab(k) and FEnha-BiCGStab, respectively. In Section 3, we focus on the solution of linear systems with multiple right-hand sides. We will recall the global version of the BiCGStab [6], which will be called the global BiCGStab (Gl-BiCGStab) method. We will also propose two improvements to this method, partial and full improvements, which will be called partial and full enhancements of the global BiCGStab and denoted by PEnha-Gl-BiCGStab(k) and FEnha-Gl-BiCGStab, respectively. In Section 4, we will recall the block version of the BiCGStab (Bl-BiCGStab) method. We will also propose two improvements to this method, partial and full improvements, which will be called partial and full enhancements of the block BiCGStab and denoted by PEnha-Bl-BiCGStab(k) and FEnha-Bl-BiCGStab, respectively.

In Section 5, we will present some numerical experiments to compare the proposed algorithms with the well-known GMRES [12], the global GMRES (Gl-GMRES) [13], and the block GMRES (Bl-GMRES) [14] methods.

Throughout this article, all vectors and matrices are assumed to be real, and the following notation is used. First, represents the transpose of any matrix . For two vectors and in the inner product is , with the Euclidean norm. In the block and global cases, we consider two matrices and in . The inner product is defined by , where denotes the trace of a square matrix . Moreover, the associated norm is the Frobenius norm noted . We denote by the identity matrix of order .

2. BiCGStab Method and Its Enhancement

The biconjugate gradient stabilized (BiCGStab) algorithm has been developed for solving nonsymmetric linear systems (1). This algorithm has been given from the conjugate gradient squared (CGS) algorithm of Sonneveld [15], which is obtained from the biconjugate gradient (BiCG) algorithm (see [4, 16]). This last algorithm was obtained by using the Lanczos biorthogonalization (see [4, 16]). All these methods are Krylov subspace methods for solving linear systems. So, in this section, we give some theoretical background and some preliminary results. For an initial guess , we associate a residual vector

Definition 1. We define the Krylov subspace of order associated to the matrix and vector by

Classical Krylov subspace methods compute the approximate solution and its corresponding residual such that where is a polynomial of degree . Let be the minimal polynomial of the matrix with respect to the vector , i.e.,

If the considered method converges after iterations, the polynomial matrix can be written as follows:

The way in which the coefficients of these two polynomials are calculated gives us the different variants of the Krylov subspace methods. For example, in the GMRES method with and is the pseudoinverse of the matrix . The BiCG method is obtained if we consider the following choice: with .

We remark that is an orthogonal projector. Then, the associated residual vector of the GMRES method is defined by an orthogonal projector, hence the optimal property of this method.

The CGS method is developed to avoid the calculation of the transpose of the matrix in the BiCG method, and then the residual associated to the CGS method is given as follows:

The CGS algorithm is based on squaring the residual polynomial, and, in cases of irregular convergence, this may lead to a substantial buildup of rounding errors or possibly even overflow. The BiCGStab algorithm is a variation of the CGS method which was developed to remedy this difficulty. Instead of seeking a method which delivers a residual vector of the form (9), BiCGStab produces iterates whose residual vectors are of the form

in which, as before, is the residual polynomial associated with the BiCG algorithm and is a new polynomial which is defined recursively at each step with the goal of “stabilizing” or “smoothing” the convergence behavior of the original algorithm. Specifically, is defined by the simple recurrence

with is determined by minimizing the norm of the residual. Based on this fact of minimization, our idea of improving the convergence of the BiCGStab method comes from the fact that we can give another form to these polynomials to improve convergence, using only the data that we have to keep the same storage and computation time.

The BiCGStab method is summarized by the following algorithm.

1. guess initial vector;
2. , , ;
3. for ;
4. ;
5. ;
6. ;
7. ;
8. ;
9. ;
10. ;
11. ;
12. ;
13. end for.

We will propose an improvement to the convergence of the BiCGStab method. Two enhancements of this method are studied: the first one will be called partial enhancement, denoted by PEnha-BiCGStab(k), and the second one will be called full enhancement, denoted by FEnha-BiCGStab. We propose to improve the convergence of the BiCGStab method by using the following well-known result.

Proposition 2. Consider the orthogonal projector where the rectangular matrix is a full rank matrix in and is its pseudoinverse (Moore-Penrose) (for more details of the pseudoinverse, see [17]). Applying the projector to any vector , we obtain a new vector, which we denote by

Then, we have

The proof of this proposition is given in [4] page 38.

To improve the convergence of an iterative method for solving linear systems, it is necessary to minimize and decrease the norm of its residual in as few iterations as possible. Then, by invoking Proposition 2 with the residual vector , we obtain an improvement in the accuracy and stability of the BiCGStab algorithm. Thus, we will apply an orthogonal projector to the residual of this method to obtain a new residual with a smaller norm. Furthermore, to avoid a storage problem, we use the pairs of vectors already calculated in the BiCGStab method to construct the orthogonal projector.

The partial enhancement of the convergence of the BiCGStab method (PEnha-BiCGStab(k)) is given by choosing equal to the last pairs of vectors , and by adding to line 10 in Algorithm 1 the following instructions (1)(2)(3)(4)(5)

The full enhancement of the convergence of the BiCGStab (FEnha-BiCGStab) method is defined by choosing equal to the all last pairs of vectors and by adding to line 10 in Algorithm 1 the following instructions (1)(2)(3)(4)(5)

3. Global BiCGStab Method and Its Enhancement

In this section, we consider the solution of large and sparse nonsymmetric systems with multiple right-hand sides of the form

where the coefficient matrix is a nonsingular real matrix of order , with . One class of solvers for solving problem (17) are the global methods, which are based on the use of a global projection process onto a matrix (global) Krylov subspace, including global the FOM and GMRES methods [13], global BCG and BiCGStab methods [10], and global Hessenberg and CMRH methods [18].

The other class is that of block solvers which are much more efficient when the matrix is relatively dense. The first block solvers are the block conjugate gradient (Bl-CG) and block biconjugate gradient (Bl-BiCG) methods proposed in [19];for nonsymmetric problems, the block generalized minimal residual (Bl-GMRES) algorithm [14], the block quasiminimum residual (Bl-QMR) algorithm [20], the block BiCGStab (Bl-BiCGStab) algorithm [10], and the block Lanczos method [7] have been developed.

In what follows, we recall the global biconjugate gradient stabilized (Gl-BiCGStab) algorithm.

1. guess initial matrix;
2. , , ;
3. for ;
4. ;
5. ;
6. ;
7. ;
8. ;
9. ;
10. ;
11. ;
12. ;
13. end for.

We will propose an improvement of the convergence of the Gl-BiCGStab method. Two enhancements of this method are studied: the first one will be called partial global enhancement, denoted by PEnha-Gl-BiCGStab(k), and the second one will be called full global enhancement, denoted by FEnha-Gl-BiCGStab. We propose to improve the convergence of the Gl-BiCGStab method by using the following well-known result.

Proposition 3. Consider the orthogonal projector where the rectangular matrix is a full rank matrix in and is its pseudoinverse (Moore-Penrose) (for more details of the pseudoinverse, see [17]). Applying the projector to any matrix , we obtain a new residual, which we denote by

Then, we have

The proof of this proposition is similar to that of the proposition in the standard case. By invoking Proposition 3 with the residual matrix , we obtain an improvement of the convergence Gl-BiCGStab algorithm. Thus, we will apply an orthogonal projector to the residual of this method and then to minimize its norm. To avoid a storage problem, we use the pairs of matrices already calculated in the Gl-BiCGStab method to construct the orthogonal projector.

The partial enhancement of the convergence of the Gl-BiCGStab method (PEnha-Gl-BiCGStab(k)) is given by choosing equal to the last pairs of matrices and by adding to line 10 in Algorithm 2 the following instructions (1)(2)(3)(4)(5)

The full enhancement of the convergence of the Gl-BiCGStab (FEnha-Gl-BiCGStab) method is defined by choosing equal to the all last pairs of matrices

and by adding to line 10 in Algorithm 1 the following instructions (1)(2)(3)(4)(5)

4. Block BiCGStab Method and Its Enhancement

As for the Gl-BiCGStab method, we will propose an improvement of the convergence of the block BiCGStab method by applying Proposition 3. Two enhancements of this method are proposed: the first one will be called the block partial enhancement, denoted by PEnha-Bl-BiCGStab(k), and the second one will be called block full enhancement, denoted by FEnha-Bl-BiCGStab.

First, let us recall the block version of BiCGStab (Bl-BiCGStab).

1. guess initial matrix;
2. , , ;
3. for ;
4. ;
5. ;
6. ;
7. ;
8. ;
9. ;
10. ;
11. ;
12. ;
13. end for.

By invoking Proposition 3 with the residual vector , we obtain an improvement of the Bl-BiCGStab algorithm. Thus, we will apply an orthogonal projector to the residual of this method. To avoid a storage problem, we use the pairs of matrices already calculated in the Bl-BiCGStab method to construct the orthogonal projector.

The partial enhancement of the convergence of the Bl-BiCGStab method (PEnha-Bl-BiCGStab(k)) is given by choosing equal to the last pairs of matrices and by adding to line 10 in Algorithm 3 the following instructions (1)(2)(3)(4)(5)

Remark that is given by Algorithm3.

The full enhancement of the convergence of the Bl-BiCGStab (FEnha-Bl-BiCGStab) method is defined by choosing equal to the all last pairs of matrices and by adding to line 10 in Algorithm 3 the following instructions (1)(2)(3)(4)(5)

We notice that at each iteration, we compute two matrices and . The choice of the matrix is crucial because we have vectors in each matrix for the partial enhancement instead of like in the standard case. Then, the convergence is clear in these two cases. In other words, if the number of vectors that we use to construct the orthogonal projector is large, we obtain a clearer improvement in accuracy and stability, and then the convergence will be faster.

5. Numerical Examples

In this section, we consider the following convection-diffusion equation

where and . The discretization of this equation is done via centered finite differences with the standard 7-point stencil in three dimensions. (For more details about the 7-point stencil, see [21].) The obtained matrix is sparse. Then, we have used an example where the existing methods converge. We have shown numerically that the new methods give an improvement to this convergence. That is what we have shown theoretically. For all the examples, we choose , , and

The order of the system is .

In general, to compare two iterative methods in terms of convergence accuracy and stability, we need to compare the history of the norm of the residual and error vectors. Numerically, sometimes, even if the residual norm gives calculation results, the error norm does not. That is why it is necessary to check this also for the new methods too. We define the residual norm and the error norm in the standard case as follows:

For the global and block case, we use the following formulas:

We would also compare the execution time of each method, but in this case, the difference between the improved and unimproved methods is negligible (see Tables 1–3). To minimize the norm of the residual, we have used the matrices and vectors already computed at each iteration.

To illustrate the efficiency of our technique, we compare the BiCGStab and its enhancement methods for systems with single right-hand sides, given by Algorithm 1, with the GMRES methods. Then, we apply the classical BiCGStab and the new enhanced BiCGStab(k) (partial and full enhancements of BiCGStab), denoted by PEnha-BiCGStab(k) and FEnha-BiCGStab for and , we give the curves of residual norms and error norms. For these methods, the right-hand of the system is chosen as follows: where is the solution of the considered system, and the rand function creates a random -vector for , with coefficients uniformly distributed in [0, 1], and the initial guess was taken to be zero. In this case, the tests were stopped as soon as . Figures 1 and 2 illustrate the comparison of these algorithms for residual and error norms, respectively. Remark that the function randn can be also used, which creates a random matrix or vector with real random coefficients.

For global and block methods, the right hand of (17) is chosen as follows: the initial guess matrix equal to . The tests were stopped as soon as

For the global case, we compare the global BiCGStab (Gl-BiCGStab) and its enhancements, partial enhancement of global BiCGStab(k) (PEnha-Gl-BiCGStab(k)), and full enhancement of global BiCGStab (FEnha-Gl-BiCGStab) for , , and for , with the Gl-GMRES method. Figures 3 and 4 give this comparison of residual and error norms, respectively.

For the block case, we compare the block BiCGStab (Bl-BiCGStab) and its enhancements, partial enhancement of block BiCGStab(s) (PEnha-Bl-BiCGStab(k)), and full enhancement of block BiCGStab (FEnha-Bl-BiCGStab) for , , and for , with the Bl-GMRES method. Figures 5 and 6 show this comparison of residual and error norms, respectively.

5.1. BiCGStab Method

From the curves of Figures 1 and 2, we remark that the residual norm and the error norm of the enhanced algorithms decrease quickly compared with the existing methods. Furthermore, for , the convergence is clearly faster than the case when .

5.2. Global BiCGStab Method

In Figures 3 and 4, we observe that the enhanced solvers FEnha-Gl-BiCGStab and PEnha-Gl-BiCGStab(k) give the best result. In this example, we can remark also that for the PEnha-Gl-BiCGStab(k) is faster. A slight improvement of stability is also observed. For , the enhanced method is still better than the Gl-BiCGStab and Gl-GMRES methods.

5.3. Block BiCGStab Method

The curves clearly show that the accuracy of the block-enhanced methods is better than that of the Bl-BiCGStab and Bl-GMRES. And this is totally normal, as we have in the global and block cases, vectors to construct an orthogonal projector. The convergence is faster as the number of vectors increases. Hence, the importance of choosing matrices for constructing the orthogonal projectors in the global and block cases. The main aim of this paper is to improve the convergence of the BiCGStab method in terms of accuracy and stability while keeping the same storage properties and computation time. Tables show that when using this technique, the difference in computation time is negligible between the BiCGStab method and the improved BiCGStab method in all three cases. Here, we have also included the GMRES method as the most optimal, to show that even though we lose a little in terms of computation time, we still get a significant improvement in accuracy. So, in Tables 1–3, we focus on the comparison between the BiCGStab method and its improved version.

6. Conclusion

In this paper, we proposed a new technique to improve the convergence behavior of the BiCGStab method for the standard, global, and block cases. Using orthogonal projectors, we have proposed an enhancement of the convergence of the BiCGStab method. The orthogonal projectors are constructed using vectors and matrices already computed in each method to avoid storage problems and then keep the advantage of storage of BiCGStab in all cases. Numerically, we see that for all three cases, the enhanced algorithms of BiCGStab are more efficient, and they converge faster than the BiCGStab and GMRES methods with negligible deference to the turnaround time.

Data Availability

The data that support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The author declares that they have no conflicts of interest.

Acknowledgments

This work was supported by the Pure and Applied Mathematics Laboratory at the University of Littoral Côte d’Opale. Open Access funding is enabled and organized by COUPERIN CY23.

References

G. Meurant and J. D. Tebbens, “Krylov methods for nonsymmetric linear systems-from theory to computation,” Springer Series in Computational Mathematics, vol. 57, 2020.
View at: Publisher Site | Google Scholar
I. S. Duff, R. G. Grimes, and J. G. Lewis, “Sparse matrix test problems,” ACM Transactions on Mathematical Software, vol. 15, no. 1, pp. 1–14, 1989.
View at: Publisher Site | Google Scholar
M. Market, Collection of Test Matrices, http://math.nist.gov/MatrixMarket/index.html.
Y. Saad, Iterative Methods for Sparse Linear System, SIAM, Philadelphia, Second edition, 2003.
View at: Publisher Site
H. A. Van Der Vorst, “Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems,” SIAM Journal on Scientific and Statistical Computing, vol. 13, no. 2, pp. 631–644, 1992.
View at: Publisher Site | Google Scholar
K. Jbilou, H. Sadok, and A. Tinzefte, “Oblique projection methods for linear systems with multiple right-hand sides,” Electronic Transactions on Numerical Analysis, vol. 20, pp. 119–138, 2005.
View at: Google Scholar
A. El Guennouni, K. Jbilou, and H. Sadok, “The block Lanczos method for linear systems with multiple right-hand sides,” Applied Numerical Mathematics, vol. 51, no. 2-3, pp. 243–256, 2004.
View at: Publisher Site | Google Scholar
F. Bouyghf, A. Messaoudi, and H. Sadok, “An enhancement of the convergence of IDR method,” ETNA–Electronic Transactions on Numerical Analysis, vol. 58, pp. 470–485, 2022.
View at: Publisher Site | Google Scholar
F. Bouyghf, A. Messaoudi, and H. Sadok, A unified approach to Krylov subspace methods for solving linear systems, Numerical Algorithms, 2023.
A. El Guennouni, K. Jbilou, and H. Sadok, “A block version of BiCGStab for linear systems with multiple right–hand sides,” Electronic Transactions on Numerical Analysis, vol. 16, pp. 129–142, 2003.
View at: Google Scholar
P. Sonneveld and M. B. Van Gijzen, “IDR(s): a family of simple and fast algorithms for solving large nonsymmetric systems of linear equations,” SIAM Journal on Scientific Computing, vol. 31, no. 2, pp. 1035–1062, 2009.
View at: Publisher Site | Google Scholar
Y. Saad and M. H. Schultz, “GMRES: a generalized minimal residual algorithm for solving non-symmetric linear systems,” SIAM Journal on Scientific and Statistical Computing, vol. 7, no. 3, pp. 856–869, 1986.
View at: Publisher Site | Google Scholar
K. Jbilou, A. Messaoudi, and H. Sadok, “Global FOM and GMRES algorithms for matrix equations,” Applied Numerical Mathematics, vol. 31, no. 1, pp. 49–63, 1999.
View at: Publisher Site | Google Scholar
B. Vital, Etude de quelques méthodes de résolution de problèmes linéaires de grande taille sur multiprocesseur, [Ph.D. thesis], Univérsité de Rennes, Rennes, France, 1990.
P. Sonneveld, “CGS: a fast Lanczos-type solver for nonsymmetric linear systems,” SIAM Journal on Scientific and Statistical Computing, vol. 10, no. 1, pp. 36–52, 1989.
View at: Publisher Site | Google Scholar
H. A. Van der Vorst, Iterative Krylov Methods for Large Linear Systems, Cambridge University Press, Cambridge, 2009.
A. Ben-Israel and T. N. E. Greville, Generalized Inverse, Theory and Application, Canadian Mathematical Society, Second edition, 2003.
M. El Heyouni, “The global Hessenberg and global CMRH methods for linear systems with multiple right-hand sides,” Numerical Algorithms, vol. 26, no. 4, pp. 317–332, 2001.
View at: Publisher Site | Google Scholar
D. O’Leary, “The block conjugate gradient algorithm and related methods,” Linear Algebra and its Applications, vol. 29, pp. 293–322, 1980.
View at: Publisher Site | Google Scholar
R. W. Freund and M. Malhotra, “A block QMR algorithm for non-Hermitian linear systems with multiple right-hand sides,” Linear Algebra and its Applications, vol. 254, no. 1-3, pp. 119–157, 1997.
View at: Publisher Site | Google Scholar
R. C. O’Reilly and J. M. Beck, “A family of large-stencil discrete Laplacian approximations in three-dimensions,” International Journal for Numerical Methods in Engineering, pp. 1–16, 2006.
View at: Google Scholar

Copyright

Copyright © 2023 F. Bouyghf. 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.

PDF Download Citation

Download other formats

Order printed copies

Views

248

Downloads

199

Citations