- Email: [email protected]
A hybrid variant of the BiCOR method for a nonsymmetric linear system with a complex spectrum
Applied Mathematics Letters 26 (2013) 457–462
Contents lists available at SciVerse ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
A hybrid variant of the BiCOR method for a nonsymmetric linear system with a complex spectrum✩ Liang Zhao, Ting-Zhu Huang ∗ School of Math. Sci., University of Electronic Science and Technology of China, Chengdu, 611731, China
article
info
Article history: Received 18 September 2012 Received in revised form 22 November 2012 Accepted 22 November 2012 Keywords: Krylov subspace method BiCOR method A-orthonormalization procedure
abstract A biconjugate A-orthogonal residual method based on the biconjugate A-orthonormalization procedure has been proposed and denoted as BiCOR by Jing, Huang, et al. [Y.-F. Jing, T.-Z. Huang, Y. Zhang, L. Li, G.-H. Cheng, Z.-G. Ren, Y. Duan, T. Sogabe, and B. Carpentieri, Lanczos-type variants of the COCR method for complex nonsymmetric linear systems, J. Comput. Phys. 228 (2009) 6376–6394]. Here, a hybrid variant of BiCOR is presented to avoid generating the transpose of the system matrix and to optimize the convergence of a nonsymmetric linear system with a complex spectrum. © 2012 Elsevier Ltd. All rights reserved.
1. Introduction In nonsymmetric linear systems, i.e., Ax = b, with coefficient matrix A nonsymmetric, we expect to take advantage of the idea used in symmetric cases. BiCOR is an evolution of CR [1] based on a dual Krylov subspace and the A-orthogonal procedure [2] to generate a symmetric structure similar to the one derived from the Lanczos biorthogonalization procedure [1]. There are many similarities between BiCOR and BiCG [3,4]. From the numerical experiments in [2,5], we could easily conclude that BiCOR and its variants, like BiCORSTAB [2] and CORS [5], are better than traditional methods, including BiCG and its variants, for solving nonsymmetric linear systems. Furthermore, BiCORSTAB and CORS could be viewed as product-type variants of BiCOR, sharing the same strategies with BiCGSTAB [6] and CGS [7]. Also BiCOR and its variants will encounter similar problems suffered in BiCG, BiCGSTAB and CGS. For example, BiCORSTAB is established upon the orthogonal precondition that (AT ψi (AT )r0′ , φj (A)r0 ) = 0, if i ̸= j, where ψi and φj are polynomials with the subscripts i and j denoting corresponding degree, and r0 is the initial residual. However, the complex roots of ψi relating to the spectrum of the system matrix will impair both BiCGSTAB and BiCORSTAB severely (relevant information can be found in [6,8,9]). BiCORSTAB, indeed, involves similar troubles in numerical experiments, so we hope to utilize the idea that solves the same problem in BiCGSTAB to optimize BiCOR in this situation. Instead of setting polynomial ψi (t ) by a linear recurrence, we choose a quadratic recurrence formula to define the polynomial in partial steps of iteration as utilized in [9]. Different from CORS and BiCORSTAB, the parameters used to define the quadratic recurrence are computed and obtained with the information of the previous two steps of iteration. Corresponding numerical experiments testify that this method accelerates the convergence when the eigenvalues of the system matrix are complex. In [2], the author has supposed to utilize the ingenious idea of BiCGSTAB2 [9] to improve BiCOR and its variants; we will propose a method to testify Jing’s supposition in this paper.
✩ This research is supported by NSFC (61170311), Chinese Univ. Specialized Research Fund for the Doctoral Program (20110185110020).
∗
Corresponding author. E-mail addresses: [email protected] (L. Zhao), [email protected], [email protected] (T.-Z. Huang).
0893-9659/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.aml.2012.11.008
458
L. Zhao, T.-Z. Huang / Applied Mathematics Letters 26 (2013) 457–462
In Section 2, the original BiCOR and BiCORSTAB will be introduced briefly. The details of the hybrid variant of BiCOR and its algorithm are presented in Section 3. The numerical experiment results are given in Section 4, and we summarize this paper in Section 5. 2. BiCOR and BiCORSTAB The BiCOR method is a method for solving real nonsymmetric or complex non-Hermitian linear systems [2,5]. Given an initial guess x0 , for solving the following linear system Ax = b,
A ∈ Rn×n , b ∈ Rn ,
(1)
popular iterative Krylov methods, such as GMRES [10] and CG [11], all aim at extracting an approximate solution in the affine subspace, i.e., x0 + Km , where Km is formed as Km (A, r0 ) = span{r0 , Ar0 , A2 r0 , . . . , Am−1 r0 },
(2)
and r0 is the initial residual, which is obtained as r0 = b − Ax0 . Different from CR, a dual linear system AT x˜ = b˜
(3)
is implicitly solved by BiCOR with the Petrov–Galerkin approach [8]. Instead, we search for the approximation xm of (1) from x0 + Km (A, v1 ) satisfying b − Axm ⊥ AT Lm (A, w1 ), where Lm (A, w1 ) is defined as Lm (A, w1 ) = span{w1 , AT w1 , (AT )2 w1 , . . . , (AT )m−1 w1 },
(4)
and the approximation x˜ m of (3) is constructed over x˜ 0 + Lm (A, w1 ) satisfying b˜ − A x˜ m ⊥ AKm (A, v1 ). Generally, we select v1 and w1 as initial residuals in the original linear system (1) and its dual system (3), respectively. We keep updating a T
direction vector pi to modify the approximation xi at each step to make it approach the exact solution. And the direction vector depends on the residual vector of the current iteration and the last direction vector. The direction vectors and the residual vectors should obey the following propositions. Proposition 1. The direction vectors pi and p˜ j , along which to approach the exact solutions of (1) and (3), satisfy A2 -orthonormality:
(˜pj , A2 pi ) =
1, 0,
i = j, i ̸= j.
Proposition 2. The residuals in linear systems (1) and (3) are A-orthogonal [2], i.e., (˜ri , Ari ) = 0 (i ̸= j). The above two propositions have been proved in [5,2], so we do not state them again here. Similar to the CR method, the two-term recurrences xi+1 = xi + αi pi ,
(5)
ri+1 = ri − αi Api ,
(6)
pi+1 = ri+1 + βi pi
(7)
lead to the updates of solutions, residuals and directions. Furthermore, the dual linear system (3) is solved in the same way; that is, x˜ i+1 = x˜ i + αi p˜ i ,
(8)
r˜i+1 = r˜i − αi A p˜ i ,
(9)
T
p˜ i+1 = r˜i+1 + βi p˜ i .
(10) T
T
Considering the above two propositions, the orthogonalization conditions that ri+1 ⊥ A L and Api+1 ⊥ A L can be used to deduce the formulations of αi and βi . Corresponding details of the BiCOR algorithm can be found in [2,5]. Following the ingenious ideas of BiCGSTAB, Jing has introduced polynomial representations of residual vectors ri and direction vectors pi in [2] with ri and pi defined as follows: ri = ψi (A)φi (A)r0 ,
pi = φi (A)πi (A)r0 ,
r˜i = ψi (AT )φi (AT )˜r0 ,
p˜ i = φi (AT )πi (AT )˜r0 ,
where r0 and r˜0 are initial residual vectors in linear system (1) and (3), respectively. φi and ψi are both Lanczos-type polynomials with degrees according with a Krylov subspace (2). Algorithm 5 in [2] presents details of BiCORSTAB, and more details of BiCORSTAB can be found in the next section; the new method proposed in this paper could be regarded as a spread of BiCORSTAB combined with the idea of BiCGSTAB2 [9].
L. Zhao, T.-Z. Huang / Applied Mathematics Letters 26 (2013) 457–462
459
3. A hybrid variant of the BiCOR method In BiCORSTAB and BiCGSTAB, the residual vector is defined as a product of two polynomials φi (t ) and ψi (t ). Recalling the original BiCOR method in the previous section, we set the initial direction vector p0 = r0 , so that we can write the residual vector generated in the nth step as rn = φn (A)r0
and
r˜n = φn (AT )˜r0 ,
based on (6) and (9). Similarly, the direction vectors pi and p˜ i could also be expressed in polynomial form; that is, pn = πn (A)r0 and p˜ n = πn (AT )˜r0 . We can obtain CORS along with this idea. A new polynomial, then, is introduced to define the residual vectors and direction vectors in the form ri = ψi (A)φi (A)r0 ,
(11)
pi = ψi (A)πi (A)r0 ,
(12)
respectively, and analogously for the dual system (3). The new polynomial ψi (t ) is determined upon optimizing (fast and smooth) the convergence of the algorithm with a certain recurrence formula. In BiCORSTAB and BiCGSTAB, the authors both take the recurrence in a linear form:
ψi+1 (t ) = (1 − ωi t )ψi (t ),
(13)
where parameter ωi will be chosen according to different criteria. Different from the choice of the recurrence pattern in BiCORSTAB, we redefine the polynomial ψi (t ) in partial steps in a quadratic recurrence form as
ψi+1 (t ) = (1 − ζi )ψi−1 (t ) + (ζi + ηi t )ψi (t ),
(14)
where parameters ζi and ηi are extracted at each iteration step based on certain regulations. In this paper, we continue to take formula (13) as the pattern to define ψi in steps labeled by an odd number and use (14) in the other steps. Or in other words, we replace part of BiCORSTAB [2] with a new method and maintain the remaining segment unchanged, alternatively using (13) and (14) in order to weaken the adverse impact brought by complex roots of φi in form (13). Because of the introduction of two new parameters ζi and ηi , furthermore, they will be determined in the light of minimizing the norm of residual vector ri , which is described as an equation in two unknowns, i.e. ri = f (ζi , ηi ). We should notice that this method is primarily found on recurrences (6) and (7). Then, let us replace ri and pi with polynomials ri = φi (A)r0 and pi = πi (A)r0 , respectively, and multiply the new polynomials ψi−1 and ψi , considering the necessary information from the previous two steps, on both sides of equations (6) and (7). Here, we merely discuss the situation at the steps labeled with even numbers. We have
ψi−1 φi+1 = ψi−1 φi − αi Aψi−1 πi , ψi φi+1 = ψi φi − αi Aψi πi , ψi πi+1 = ψi φi+1 + βi ψi πi .
(15) (16) (17)
Next, we express (11) and (12) with the quadratic recurrence formula (14), which are presented as ri+1 = ψi+1 φi+1 r0 = [(1 − ζi )ψi−1 φi+1 + ζi ψi φi+1 + ηi Aψi φi+1 ]r0
(18)
pi+1 = ψi+1 πi+1 r0 = [ψi+1 φi+1 − βi+1 (1 − ζi )ψi−1 πi − βi+1 (ζi + ηi A)ψi πi ]r0 .
(19)
The next step is to determine the way to obtain parameters αi , βi , ζi and ηi . The four parameters can be divided into two classes, one depending on the original recurrences (5)–(7) including αi and βi , and the other depending on the quadratic recurrence (14), covering ζi and ηi . Based on the orthogonal properties of residuals and direction vectors, illustrated in Propositions 1 and 2, now we deduce the first class as follows:
(φi (AT )˜r0 , Aφi (A)r0 ) (˜ri , Ari ) = T T (A p˜ i , Api ) (A πi (AT )˜r0 , Aπi (A)r0 ) T (ψi (AT )˜r0 , Aφi (A)r0 ) (˜r0 , Aψi (A)φi (A)r0 ) (φi (A )˜r0 , Aφi (A)r0 ) = T = = . T T T (A φi (A )˜r0 , Aπi (A)r0 ) (A ψi (A )˜r0 , Aπi (A)r0 ) (˜r0 , A2 ψi (A)πi (A)r0 )
αi =
Recalling (11) and (12), and then defining ρi = (˜r0 , Ari ), we easily have the form of αi :
αi =
ρi . (˜r0 , A2 pi )
(20)
If we set the leading coefficients of polynomials ψi and φi to ai and bi , respectively, from (6) and (14), we have the following relations: bi+1 = −αi bi , and ai+1 = ηi ai . Hence, βi , differing from BiCORSTAB, is structured as
βi = −
ρi αi−1 × . ρi−1 ηi−1
(21)
460
L. Zhao, T.-Z. Huang / Applied Mathematics Letters 26 (2013) 457–462
Then we have obtained the formulas to compute the parameters of the first class, and the way to gain ζi and ηi will be presented in the following. If we denote (15)–(17) by fi = ψi−2 φi , si = ψi φi+1 , combining with (18), we have that ri+1 = Asi+1 ηi + (si+1 − fi+1 )ζi + fi+1 , so ζi and ηi are determined in order to minimize the norm of ri , i.e.,
∥ri+1 ∥2 = ∥Asi+1 ηi + (si+1 − fi+1 )ζi + fi+1 ∥2 ,
(22)
in which ζi and ηi are variables. The hybrid variant of BiCOR, named BiCORSTAB2 here, is presented below. Algorithm 1. BiCORSTAB2 1. initialization: select initial guess x0 and r0 = b − Ax0 ρ select r˜0 so that (˜r0 , Ar0 ) ̸= 0, set ρ0 = (˜r0 , Ar0 ), p0 = r0 , α0 = (˜r ,A02 r ) 0 0 2. for i = 0, 1, ... until convergence ρi 3. αi = (˜r ,ˆq ) 0 i 4. si+1 = ri − αi qi 5. sˆi+1 = rˆ − αi qˆ i 6. fi+1 = si − αi hˆ i 7. if at odd steps, we compute parameter ω in the same way as BiCORSTAB 8. ri+1 = si+1 − ωˆsi+1 9. xi+1 = xi + αi pi + ωsi+1 10. rˆ = Ari+1 11. ρi+1 = (˜r0 , ri+1 ) ρ 12. βi+1 = iρ+i 1 × αωi 13. pi+1 = ri+1 + βi+1 (pi − ωqi ) 14. qi+1 = rˆ + βi+1 (qi − ωˆqi ) 15. else we compute parameters ζ and η according to (22) 16. ri+1 = (1 − ζ )fi+1 + ζ si+1 + ηˆsi+1 17. xi+1 = (1 − ζ )(xi−1 + αi−1 pi−1 + αi hi ) + ζ (xi + αi pi ) − ηsi+1 18. rˆ = Ari+1 19. ρi+1 = (˜r0 , ri+1 ) ρ 20. βi+1 = − iρ+i 1 × αηi 21. pi+1 = ri+1 + βi+1 ((1 − ζ )hi + ζ pi + ηqi ) 22. qi+1 = rˆ − βi+1 ((1 − ζ )hˆ i + ζ qi + ηˆqi ) 23. end 24. qˆ i+1 = Aqi+1 25. hi+1 = si+1 − βi+1 pi 26. hˆ i+1 = sˆi+1 − βi+1 qi 27. end
We take advantage of suitable variables to replace the matrix–vector multiplications encountered throughout this method. Algorithm 1 describes such a method and the variables characterized with symbol ˆ are introduced to eliminate matrix–vector multiplications, following the same way proposed by Sogabe in [12] and used by Jing in [2] and Carpentieri in [5]. Numerical experiments will be presented in the following section.
4. Numerical experiments In this section, we take three numerical examples to test BiCORSTAB2. Here, two numerical examples used in [9] to testify BiCGSTAB2 will be utilized again because BiCORSTAB2 and BiCGSTAB2 share a similar executing method. The first two examples are banded Toeplitz matrices derived from the discretization of partial differential equations. BiCORSTAB2 is a method based on BiCOR and BiCORSTAB for making convergence stable and smooth, so we set BiCOR, CORS, and BiCORSTAB as reference solvers in our comparison. Moreover, we take cz628, borrowed in the MATLAB format from the University of Florida Sparse Matrix Collection provided by Davis [13], as another numerical example, which belongs to 2D/3D problems.
4 1
First, we give the first numerical example in the form A =
√
√
−2 4
..
.
−2 .. . 1
..
.
4 1
−2
, regarded as the system matrix with
4
spectrum located between 4 − 2i 2 and 4 + 2i 2. The right-hand side b in (1) is selected as a random vector. Fig. 1 displays ∥b−Ax ∥ the process of iteration, where the convergence history, setting log10 ∥b∥2i 2 here, with respect to iteration numbers. The first two examples are both of size 300. From Fig. 1, we can observe easily that BiCORSTAB2, faster and smoother, reaches a certain accuracy within fewer steps and without more fluctuations.
L. Zhao, T.-Z. Huang / Applied Mathematics Letters 26 (2013) 457–462
461
Fig. 1. Relative residual versus iterations under BiCOR, CORS, BiCORSTAB and BiCORSTAB2 on example 1.
Fig. 2. Relative residual versus iterations under BiCOR, CORS, BiCORSTAB and BiCORSTAB2 on example 2.
2 0 1
The second numerical example is presented as A =
1 2 0
1 2
1
..
..
.
1
.
0 1
..
.
2 0
, and we also set the right-hand side b as
1 2
a random vector. The performances of BiCORSTAB2 and its reference solvers in our comparison are illustrated in Fig. 2. We can see that its convergence is obviously faster than that of BiCORSTAB and the original BiCOR. Although the performance of CORS is similar to BiCORSTAB2, BiCORSTAB2 is a little faster than CORS in reaching a given accuracy (a two-step advantage: 39 steps in BiCORSTAB2 versus 41 steps in CORS), and the result of BiCORSTAB2, notably, is smoother than the one of CORS. Hence, BiCORSTAB2 keeps a more optimal result than the other three methods on the first two numerical examples. cz628, the third numerical example, is of size 628 × 628 with 6346 nonzero elements and a complex spectrum (group 2D/3D problems in the University of Florida Sparse Matrix Collection). Here we set the right-hand side b by b = Ae, where e is a vector of size n × 1 with elements all unity. Fig. 3 shows the consequences of BiCORSTAB2 and its reference solvers in our comparison. On the second numerical example, the performances of BiCORSTAB2 and CORS are similar while BiCORSTAB2 has a two-step advantage. But BiCORSTAB2 converges much faster than CORS on the third numerical example, and analogously for BiCOR. This is obviously illustrated in Fig. 3. The iteration process of BiCORSTAB is close to that of BiCORSTAB2 but it converges a little slower than BiCORSTAB2 for the same given accuracy.
462
L. Zhao, T.-Z. Huang / Applied Mathematics Letters 26 (2013) 457–462
Fig. 3. Relative residual versus iterations under BiCOR, CORS, BiCORSTAB and BiCORSTAB2 on example 3.
5. Summary In this paper, we recall basic concepts about BiCOR and one of its variants. To optimize the original BiCOR method, we take advantage of the idea of BiCGSTAB2 to overcome the problems encountered in nonsymmetric linear systems with a complex spectrum. The new method maintains part of BiCORSTAB and introduces a quadratic recurrence to define the polynomial ψi in other steps. The BiCORSTAB2 method, a hybrid variant of BiCOR, displays a better result than other three methods and is stable under different numerical examples. Furthermore, we expect to search out other variants of BiCOR to optimize it with ideas of the variants of BiCG. The original BiCOR method has shown us better numerical consequences than those of BiCG in [2,5], so we hope that the ingenious ideas of variants of BiCG will bring us more surprising inspirations to modify and optimize the BiCOR method. Acknowledgments The authors would like to express their great thankfulness to the referees and editor Prof. Alan Tucker for their very helpful suggestions for revising this manuscript. References [1] Y. Saad, Iterative Methods for Sparse Linear Systems, second ed., SIAM, Philadelphia, 2003. [2] Y.-F. Jing, T.-Z. Huang, Y. Zhang, L. Li, G.-H. Cheng, Z.-G. Ren, Y. Duan, T. Sogabe, B. Carpentieri, Lanczos-type variants of the COCR method for complex nonsymmetric linear systems, J. Comput. Phys. 228 (2009) 6376–6394. [3] R. Fletcher, Conjugate Gradient Methods for Indefinite Systems, in: Lecture Notes in Math., vol. 506, Springer-Verlag, Berlin, Heidelberg, New York, 1976, pp. 73–89. [4] C. Lanczos, Solution of systems of linear equations by minimized iterations, J. Res. Natl. Bur. Stand. 49 (1952) 33–53. [5] B. Carpentieri, Y.-F. Jing, T.-Z. Huang, The BICOR and CORS iterative algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Comput. 33 (2011) 3020–3036. [6] H.A. van der Vorst, Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 13 (1992) 631–644. [7] P. Sonneveld, CGS: a fast Lanczos-type solver for nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 10 (1989) 36–52. [8] Henk A. van der Vorst, Iterative Krylov Methods for Large Linear Systems, Cambridge University Press, Cambridge, UK, 2003. [9] M.H. Gutknecht, Variant of BiCGSTAB for matrices with complex spectrum, SIAM J. Sci. Comput. 14 (1993) 1020–1033. [10] Y. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear system, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869. [11] M.R. Hestenes, E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Natl. Bur. Stand. 49 (1952) 409–436. [12] T. Sogabe, Extensions of the conjugate residual method, Ph.D. Dissertation, University of Tokyo, 2006. [13] T. Davis, The University of Florida sparse matrix collection. NA Digest 97 (23) 1997. http://www.cise.ufl.edu/research/sparse/matrices.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
A hybrid variant of the BiCOR method for a nonsymmetric linear system with a complex spectrum✩ Liang Zhao, Ting-Zhu Huang ∗ School of Math. Sci., University of Electronic Science and Technology of China, Chengdu, 611731, China
article
info
Article history: Received 18 September 2012 Received in revised form 22 November 2012 Accepted 22 November 2012 Keywords: Krylov subspace method BiCOR method A-orthonormalization procedure
abstract A biconjugate A-orthogonal residual method based on the biconjugate A-orthonormalization procedure has been proposed and denoted as BiCOR by Jing, Huang, et al. [Y.-F. Jing, T.-Z. Huang, Y. Zhang, L. Li, G.-H. Cheng, Z.-G. Ren, Y. Duan, T. Sogabe, and B. Carpentieri, Lanczos-type variants of the COCR method for complex nonsymmetric linear systems, J. Comput. Phys. 228 (2009) 6376–6394]. Here, a hybrid variant of BiCOR is presented to avoid generating the transpose of the system matrix and to optimize the convergence of a nonsymmetric linear system with a complex spectrum. © 2012 Elsevier Ltd. All rights reserved.
1. Introduction In nonsymmetric linear systems, i.e., Ax = b, with coefficient matrix A nonsymmetric, we expect to take advantage of the idea used in symmetric cases. BiCOR is an evolution of CR [1] based on a dual Krylov subspace and the A-orthogonal procedure [2] to generate a symmetric structure similar to the one derived from the Lanczos biorthogonalization procedure [1]. There are many similarities between BiCOR and BiCG [3,4]. From the numerical experiments in [2,5], we could easily conclude that BiCOR and its variants, like BiCORSTAB [2] and CORS [5], are better than traditional methods, including BiCG and its variants, for solving nonsymmetric linear systems. Furthermore, BiCORSTAB and CORS could be viewed as product-type variants of BiCOR, sharing the same strategies with BiCGSTAB [6] and CGS [7]. Also BiCOR and its variants will encounter similar problems suffered in BiCG, BiCGSTAB and CGS. For example, BiCORSTAB is established upon the orthogonal precondition that (AT ψi (AT )r0′ , φj (A)r0 ) = 0, if i ̸= j, where ψi and φj are polynomials with the subscripts i and j denoting corresponding degree, and r0 is the initial residual. However, the complex roots of ψi relating to the spectrum of the system matrix will impair both BiCGSTAB and BiCORSTAB severely (relevant information can be found in [6,8,9]). BiCORSTAB, indeed, involves similar troubles in numerical experiments, so we hope to utilize the idea that solves the same problem in BiCGSTAB to optimize BiCOR in this situation. Instead of setting polynomial ψi (t ) by a linear recurrence, we choose a quadratic recurrence formula to define the polynomial in partial steps of iteration as utilized in [9]. Different from CORS and BiCORSTAB, the parameters used to define the quadratic recurrence are computed and obtained with the information of the previous two steps of iteration. Corresponding numerical experiments testify that this method accelerates the convergence when the eigenvalues of the system matrix are complex. In [2], the author has supposed to utilize the ingenious idea of BiCGSTAB2 [9] to improve BiCOR and its variants; we will propose a method to testify Jing’s supposition in this paper.
✩ This research is supported by NSFC (61170311), Chinese Univ. Specialized Research Fund for the Doctoral Program (20110185110020).
∗
Corresponding author. E-mail addresses: [email protected] (L. Zhao), [email protected], [email protected] (T.-Z. Huang).
0893-9659/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.aml.2012.11.008
458
L. Zhao, T.-Z. Huang / Applied Mathematics Letters 26 (2013) 457–462
In Section 2, the original BiCOR and BiCORSTAB will be introduced briefly. The details of the hybrid variant of BiCOR and its algorithm are presented in Section 3. The numerical experiment results are given in Section 4, and we summarize this paper in Section 5. 2. BiCOR and BiCORSTAB The BiCOR method is a method for solving real nonsymmetric or complex non-Hermitian linear systems [2,5]. Given an initial guess x0 , for solving the following linear system Ax = b,
A ∈ Rn×n , b ∈ Rn ,
(1)
popular iterative Krylov methods, such as GMRES [10] and CG [11], all aim at extracting an approximate solution in the affine subspace, i.e., x0 + Km , where Km is formed as Km (A, r0 ) = span{r0 , Ar0 , A2 r0 , . . . , Am−1 r0 },
(2)
and r0 is the initial residual, which is obtained as r0 = b − Ax0 . Different from CR, a dual linear system AT x˜ = b˜
(3)
is implicitly solved by BiCOR with the Petrov–Galerkin approach [8]. Instead, we search for the approximation xm of (1) from x0 + Km (A, v1 ) satisfying b − Axm ⊥ AT Lm (A, w1 ), where Lm (A, w1 ) is defined as Lm (A, w1 ) = span{w1 , AT w1 , (AT )2 w1 , . . . , (AT )m−1 w1 },
(4)
and the approximation x˜ m of (3) is constructed over x˜ 0 + Lm (A, w1 ) satisfying b˜ − A x˜ m ⊥ AKm (A, v1 ). Generally, we select v1 and w1 as initial residuals in the original linear system (1) and its dual system (3), respectively. We keep updating a T
direction vector pi to modify the approximation xi at each step to make it approach the exact solution. And the direction vector depends on the residual vector of the current iteration and the last direction vector. The direction vectors and the residual vectors should obey the following propositions. Proposition 1. The direction vectors pi and p˜ j , along which to approach the exact solutions of (1) and (3), satisfy A2 -orthonormality:
(˜pj , A2 pi ) =
1, 0,
i = j, i ̸= j.
Proposition 2. The residuals in linear systems (1) and (3) are A-orthogonal [2], i.e., (˜ri , Ari ) = 0 (i ̸= j). The above two propositions have been proved in [5,2], so we do not state them again here. Similar to the CR method, the two-term recurrences xi+1 = xi + αi pi ,
(5)
ri+1 = ri − αi Api ,
(6)
pi+1 = ri+1 + βi pi
(7)
lead to the updates of solutions, residuals and directions. Furthermore, the dual linear system (3) is solved in the same way; that is, x˜ i+1 = x˜ i + αi p˜ i ,
(8)
r˜i+1 = r˜i − αi A p˜ i ,
(9)
T
p˜ i+1 = r˜i+1 + βi p˜ i .
(10) T
T
Considering the above two propositions, the orthogonalization conditions that ri+1 ⊥ A L and Api+1 ⊥ A L can be used to deduce the formulations of αi and βi . Corresponding details of the BiCOR algorithm can be found in [2,5]. Following the ingenious ideas of BiCGSTAB, Jing has introduced polynomial representations of residual vectors ri and direction vectors pi in [2] with ri and pi defined as follows: ri = ψi (A)φi (A)r0 ,
pi = φi (A)πi (A)r0 ,
r˜i = ψi (AT )φi (AT )˜r0 ,
p˜ i = φi (AT )πi (AT )˜r0 ,
where r0 and r˜0 are initial residual vectors in linear system (1) and (3), respectively. φi and ψi are both Lanczos-type polynomials with degrees according with a Krylov subspace (2). Algorithm 5 in [2] presents details of BiCORSTAB, and more details of BiCORSTAB can be found in the next section; the new method proposed in this paper could be regarded as a spread of BiCORSTAB combined with the idea of BiCGSTAB2 [9].
L. Zhao, T.-Z. Huang / Applied Mathematics Letters 26 (2013) 457–462
459
3. A hybrid variant of the BiCOR method In BiCORSTAB and BiCGSTAB, the residual vector is defined as a product of two polynomials φi (t ) and ψi (t ). Recalling the original BiCOR method in the previous section, we set the initial direction vector p0 = r0 , so that we can write the residual vector generated in the nth step as rn = φn (A)r0
and
r˜n = φn (AT )˜r0 ,
based on (6) and (9). Similarly, the direction vectors pi and p˜ i could also be expressed in polynomial form; that is, pn = πn (A)r0 and p˜ n = πn (AT )˜r0 . We can obtain CORS along with this idea. A new polynomial, then, is introduced to define the residual vectors and direction vectors in the form ri = ψi (A)φi (A)r0 ,
(11)
pi = ψi (A)πi (A)r0 ,
(12)
respectively, and analogously for the dual system (3). The new polynomial ψi (t ) is determined upon optimizing (fast and smooth) the convergence of the algorithm with a certain recurrence formula. In BiCORSTAB and BiCGSTAB, the authors both take the recurrence in a linear form:
ψi+1 (t ) = (1 − ωi t )ψi (t ),
(13)
where parameter ωi will be chosen according to different criteria. Different from the choice of the recurrence pattern in BiCORSTAB, we redefine the polynomial ψi (t ) in partial steps in a quadratic recurrence form as
ψi+1 (t ) = (1 − ζi )ψi−1 (t ) + (ζi + ηi t )ψi (t ),
(14)
where parameters ζi and ηi are extracted at each iteration step based on certain regulations. In this paper, we continue to take formula (13) as the pattern to define ψi in steps labeled by an odd number and use (14) in the other steps. Or in other words, we replace part of BiCORSTAB [2] with a new method and maintain the remaining segment unchanged, alternatively using (13) and (14) in order to weaken the adverse impact brought by complex roots of φi in form (13). Because of the introduction of two new parameters ζi and ηi , furthermore, they will be determined in the light of minimizing the norm of residual vector ri , which is described as an equation in two unknowns, i.e. ri = f (ζi , ηi ). We should notice that this method is primarily found on recurrences (6) and (7). Then, let us replace ri and pi with polynomials ri = φi (A)r0 and pi = πi (A)r0 , respectively, and multiply the new polynomials ψi−1 and ψi , considering the necessary information from the previous two steps, on both sides of equations (6) and (7). Here, we merely discuss the situation at the steps labeled with even numbers. We have
ψi−1 φi+1 = ψi−1 φi − αi Aψi−1 πi , ψi φi+1 = ψi φi − αi Aψi πi , ψi πi+1 = ψi φi+1 + βi ψi πi .
(15) (16) (17)
Next, we express (11) and (12) with the quadratic recurrence formula (14), which are presented as ri+1 = ψi+1 φi+1 r0 = [(1 − ζi )ψi−1 φi+1 + ζi ψi φi+1 + ηi Aψi φi+1 ]r0
(18)
pi+1 = ψi+1 πi+1 r0 = [ψi+1 φi+1 − βi+1 (1 − ζi )ψi−1 πi − βi+1 (ζi + ηi A)ψi πi ]r0 .
(19)
The next step is to determine the way to obtain parameters αi , βi , ζi and ηi . The four parameters can be divided into two classes, one depending on the original recurrences (5)–(7) including αi and βi , and the other depending on the quadratic recurrence (14), covering ζi and ηi . Based on the orthogonal properties of residuals and direction vectors, illustrated in Propositions 1 and 2, now we deduce the first class as follows:
(φi (AT )˜r0 , Aφi (A)r0 ) (˜ri , Ari ) = T T (A p˜ i , Api ) (A πi (AT )˜r0 , Aπi (A)r0 ) T (ψi (AT )˜r0 , Aφi (A)r0 ) (˜r0 , Aψi (A)φi (A)r0 ) (φi (A )˜r0 , Aφi (A)r0 ) = T = = . T T T (A φi (A )˜r0 , Aπi (A)r0 ) (A ψi (A )˜r0 , Aπi (A)r0 ) (˜r0 , A2 ψi (A)πi (A)r0 )
αi =
Recalling (11) and (12), and then defining ρi = (˜r0 , Ari ), we easily have the form of αi :
αi =
ρi . (˜r0 , A2 pi )
(20)
If we set the leading coefficients of polynomials ψi and φi to ai and bi , respectively, from (6) and (14), we have the following relations: bi+1 = −αi bi , and ai+1 = ηi ai . Hence, βi , differing from BiCORSTAB, is structured as
βi = −
ρi αi−1 × . ρi−1 ηi−1
(21)
460
L. Zhao, T.-Z. Huang / Applied Mathematics Letters 26 (2013) 457–462
Then we have obtained the formulas to compute the parameters of the first class, and the way to gain ζi and ηi will be presented in the following. If we denote (15)–(17) by fi = ψi−2 φi , si = ψi φi+1 , combining with (18), we have that ri+1 = Asi+1 ηi + (si+1 − fi+1 )ζi + fi+1 , so ζi and ηi are determined in order to minimize the norm of ri , i.e.,
∥ri+1 ∥2 = ∥Asi+1 ηi + (si+1 − fi+1 )ζi + fi+1 ∥2 ,
(22)
in which ζi and ηi are variables. The hybrid variant of BiCOR, named BiCORSTAB2 here, is presented below. Algorithm 1. BiCORSTAB2 1. initialization: select initial guess x0 and r0 = b − Ax0 ρ select r˜0 so that (˜r0 , Ar0 ) ̸= 0, set ρ0 = (˜r0 , Ar0 ), p0 = r0 , α0 = (˜r ,A02 r ) 0 0 2. for i = 0, 1, ... until convergence ρi 3. αi = (˜r ,ˆq ) 0 i 4. si+1 = ri − αi qi 5. sˆi+1 = rˆ − αi qˆ i 6. fi+1 = si − αi hˆ i 7. if at odd steps, we compute parameter ω in the same way as BiCORSTAB 8. ri+1 = si+1 − ωˆsi+1 9. xi+1 = xi + αi pi + ωsi+1 10. rˆ = Ari+1 11. ρi+1 = (˜r0 , ri+1 ) ρ 12. βi+1 = iρ+i 1 × αωi 13. pi+1 = ri+1 + βi+1 (pi − ωqi ) 14. qi+1 = rˆ + βi+1 (qi − ωˆqi ) 15. else we compute parameters ζ and η according to (22) 16. ri+1 = (1 − ζ )fi+1 + ζ si+1 + ηˆsi+1 17. xi+1 = (1 − ζ )(xi−1 + αi−1 pi−1 + αi hi ) + ζ (xi + αi pi ) − ηsi+1 18. rˆ = Ari+1 19. ρi+1 = (˜r0 , ri+1 ) ρ 20. βi+1 = − iρ+i 1 × αηi 21. pi+1 = ri+1 + βi+1 ((1 − ζ )hi + ζ pi + ηqi ) 22. qi+1 = rˆ − βi+1 ((1 − ζ )hˆ i + ζ qi + ηˆqi ) 23. end 24. qˆ i+1 = Aqi+1 25. hi+1 = si+1 − βi+1 pi 26. hˆ i+1 = sˆi+1 − βi+1 qi 27. end
We take advantage of suitable variables to replace the matrix–vector multiplications encountered throughout this method. Algorithm 1 describes such a method and the variables characterized with symbol ˆ are introduced to eliminate matrix–vector multiplications, following the same way proposed by Sogabe in [12] and used by Jing in [2] and Carpentieri in [5]. Numerical experiments will be presented in the following section.
4. Numerical experiments In this section, we take three numerical examples to test BiCORSTAB2. Here, two numerical examples used in [9] to testify BiCGSTAB2 will be utilized again because BiCORSTAB2 and BiCGSTAB2 share a similar executing method. The first two examples are banded Toeplitz matrices derived from the discretization of partial differential equations. BiCORSTAB2 is a method based on BiCOR and BiCORSTAB for making convergence stable and smooth, so we set BiCOR, CORS, and BiCORSTAB as reference solvers in our comparison. Moreover, we take cz628, borrowed in the MATLAB format from the University of Florida Sparse Matrix Collection provided by Davis [13], as another numerical example, which belongs to 2D/3D problems.
4 1
First, we give the first numerical example in the form A =
√
√
−2 4
..
.
−2 .. . 1
..
.
4 1
−2
, regarded as the system matrix with
4
spectrum located between 4 − 2i 2 and 4 + 2i 2. The right-hand side b in (1) is selected as a random vector. Fig. 1 displays ∥b−Ax ∥ the process of iteration, where the convergence history, setting log10 ∥b∥2i 2 here, with respect to iteration numbers. The first two examples are both of size 300. From Fig. 1, we can observe easily that BiCORSTAB2, faster and smoother, reaches a certain accuracy within fewer steps and without more fluctuations.
L. Zhao, T.-Z. Huang / Applied Mathematics Letters 26 (2013) 457–462
461
Fig. 1. Relative residual versus iterations under BiCOR, CORS, BiCORSTAB and BiCORSTAB2 on example 1.
Fig. 2. Relative residual versus iterations under BiCOR, CORS, BiCORSTAB and BiCORSTAB2 on example 2.
2 0 1
The second numerical example is presented as A =
1 2 0
1 2
1
..
..
.
1
.
0 1
..
.
2 0
, and we also set the right-hand side b as
1 2
a random vector. The performances of BiCORSTAB2 and its reference solvers in our comparison are illustrated in Fig. 2. We can see that its convergence is obviously faster than that of BiCORSTAB and the original BiCOR. Although the performance of CORS is similar to BiCORSTAB2, BiCORSTAB2 is a little faster than CORS in reaching a given accuracy (a two-step advantage: 39 steps in BiCORSTAB2 versus 41 steps in CORS), and the result of BiCORSTAB2, notably, is smoother than the one of CORS. Hence, BiCORSTAB2 keeps a more optimal result than the other three methods on the first two numerical examples. cz628, the third numerical example, is of size 628 × 628 with 6346 nonzero elements and a complex spectrum (group 2D/3D problems in the University of Florida Sparse Matrix Collection). Here we set the right-hand side b by b = Ae, where e is a vector of size n × 1 with elements all unity. Fig. 3 shows the consequences of BiCORSTAB2 and its reference solvers in our comparison. On the second numerical example, the performances of BiCORSTAB2 and CORS are similar while BiCORSTAB2 has a two-step advantage. But BiCORSTAB2 converges much faster than CORS on the third numerical example, and analogously for BiCOR. This is obviously illustrated in Fig. 3. The iteration process of BiCORSTAB is close to that of BiCORSTAB2 but it converges a little slower than BiCORSTAB2 for the same given accuracy.
462
L. Zhao, T.-Z. Huang / Applied Mathematics Letters 26 (2013) 457–462
Fig. 3. Relative residual versus iterations under BiCOR, CORS, BiCORSTAB and BiCORSTAB2 on example 3.
5. Summary In this paper, we recall basic concepts about BiCOR and one of its variants. To optimize the original BiCOR method, we take advantage of the idea of BiCGSTAB2 to overcome the problems encountered in nonsymmetric linear systems with a complex spectrum. The new method maintains part of BiCORSTAB and introduces a quadratic recurrence to define the polynomial ψi in other steps. The BiCORSTAB2 method, a hybrid variant of BiCOR, displays a better result than other three methods and is stable under different numerical examples. Furthermore, we expect to search out other variants of BiCOR to optimize it with ideas of the variants of BiCG. The original BiCOR method has shown us better numerical consequences than those of BiCG in [2,5], so we hope that the ingenious ideas of variants of BiCG will bring us more surprising inspirations to modify and optimize the BiCOR method. Acknowledgments The authors would like to express their great thankfulness to the referees and editor Prof. Alan Tucker for their very helpful suggestions for revising this manuscript. References [1] Y. Saad, Iterative Methods for Sparse Linear Systems, second ed., SIAM, Philadelphia, 2003. [2] Y.-F. Jing, T.-Z. Huang, Y. Zhang, L. Li, G.-H. Cheng, Z.-G. Ren, Y. Duan, T. Sogabe, B. Carpentieri, Lanczos-type variants of the COCR method for complex nonsymmetric linear systems, J. Comput. Phys. 228 (2009) 6376–6394. [3] R. Fletcher, Conjugate Gradient Methods for Indefinite Systems, in: Lecture Notes in Math., vol. 506, Springer-Verlag, Berlin, Heidelberg, New York, 1976, pp. 73–89. [4] C. Lanczos, Solution of systems of linear equations by minimized iterations, J. Res. Natl. Bur. Stand. 49 (1952) 33–53. [5] B. Carpentieri, Y.-F. Jing, T.-Z. Huang, The BICOR and CORS iterative algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Comput. 33 (2011) 3020–3036. [6] H.A. van der Vorst, Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 13 (1992) 631–644. [7] P. Sonneveld, CGS: a fast Lanczos-type solver for nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 10 (1989) 36–52. [8] Henk A. van der Vorst, Iterative Krylov Methods for Large Linear Systems, Cambridge University Press, Cambridge, UK, 2003. [9] M.H. Gutknecht, Variant of BiCGSTAB for matrices with complex spectrum, SIAM J. Sci. Comput. 14 (1993) 1020–1033. [10] Y. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear system, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869. [11] M.R. Hestenes, E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Natl. Bur. Stand. 49 (1952) 409–436. [12] T. Sogabe, Extensions of the conjugate residual method, Ph.D. Dissertation, University of Tokyo, 2006. [13] T. Davis, The University of Florida sparse matrix collection. NA Digest 97 (23) 1997. http://www.cise.ufl.edu/research/sparse/matrices.