# Exploiting the Composite Step Strategy to the BiconjugateA-Orthogonal Residual Method for Non-Hermitian Linear Systems

Authors: Yan-Fei Jing; Ting-Zhu Huang; Bruno Carpentieri; Yong Duan
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

The BiconjugateA-Orthogonal Residual (BiCOR) method carried out in
finite precision arithmetic by means of the biconjugate A-orthonormalization
procedure may possibly tend to suffer from two sources of numerical
instability, known as two kinds of breakdowns, similarly to those of the
Biconjugate Gradient (BCG) method. This paper naturally exploits the
composite step strategy employed in the development of the composite
step BCG (CSBCG) method into the BiCOR method to cure one of the
breakdowns called as pivot breakdown. Analogously to the CSBCG method,
the resulting interesting variant, with only a minor modification to the
usual implementation of the BiCOR method, is able to avoid near pivot
breakdowns and compute all the well-defined BiCOR iterates stably on the
assumption that the underlying biconjugate A-orthonormalization procedure
does not break down. Another benefit acquired is that it seems to be a
viable algorithm providing some further practically desired smoothing of
the convergence history of the norm of the residuals, which is justified
by numerical experiments. In addition, the exhibited method inherits
the promising advantages of the empirically observed stability and fast
convergence rate of the BiCOR method over the BCG method so that it
outperforms the CSBCG method to some extent.


---

## Body

## 1. Introduction

The computational cost of many simulations with integral equations or partial differential equations that model the process is dominated by the solution of systems of linear equations of the formAx=b. The field of iterative methods for solving linear systems has been observed as an explosion of activity spurred by demand due to extraordinary technological advances in engineering and sciences in the past half century [1]. Krylov subspace methods belong to one of the most widespread and extensively accepted techniques for iterative solution of today’s large-scale linear systems [2]. With respect to “the greatest influence on the development and practice of science and engineering in the 20th century” as written by Dongarra and Sullivan [3], Krylov subspace methods are considered as one of the “Top Ten Algorithms of the Century.”Many advances in the development of Krylov subspace methods have been inspired and made by the study of even more effective approaches to linear systems. Various methods differ in the way they extract information from Krylov spaces. A basis for the underlying Krylov subspace to form iterative recurrences involved is usually constructed based on the state-of-the-art Arnoldi orthogonalization procedure or Lanczos biorthogonalization procedure; see recent excellent and thorough review articles by Simoncini and Szyld [2] and by Philippe and Reichel [4] or monographs such as those of Greenbaum [5] and Saad [6]. A novel Lanczos-type biconjugate A-orthonormalization procedure has recently been established to give birth to a new family of efficient short-recurrence methods for large real nonsymmetric and complex non-Hermitian systems of linear equations, named as the Lanczos biconjugate A-orthonormalization methods [7].As observed from numerous numerical experiments carried out with the Lanczos biconjugateA-orthonormalization methods, it has been numerically demonstrated that this family of solvers shows competitive convergence properties, is cheap in memory as it is derived from short-term vector recurrences, is parameter-free, and does not require a symmetric preconditioner like other methods, if A is symmetric indefinite [8–10]. The efficiency of the BiCOR method is highlighted by the performance profiles on a set of 14 sparse matrix problems of size up to 1.5 M unknowns shown in [10].On the other hand, this family of solvers is often faced with apparently irregular convergence behaviors appearing as “spikes” in the convergence history of the norm of the residuals, possibly leading to substantial build-up of rounding errors and worse approximate solutions, or possibly even overflow [7, 9]. Therefore, it is quite necessary to tackle their irregular convergence properties to obtain more stabilized variants so as to improve the accuracy of the desired numerical physical solutions. It is emphasized that our main attention in this paper is focused on the straightforward natural enhancement of the Biconjugate A-Orthogonal Residual (BiCOR) method, which is the basic underlying variant of the Lanczos biconjugate A-Orthonormalization methods [7]. The improvement of the BiCOR method will simultaneously result in analogous improvements for the other two evolving variants known as the Conjugate A-orthogonal Residual Squared (CORS) method and the Biconjugate A-Orthogonal Residual Stabilized (BiCORSTAB) method.By comparison of the biconjugateA-orthonormalization procedure [7] and the Lanczos biorthogonalization procedure [6], respectively, for the BiCOR method [7] and the Biconjugate Gradient (BCG) method [11] as well as their corresponding implementation algorithms, it is obviously observed that when carried out in finite precision arithmetic, the BiCOR method does tend to suffer from two possible sources of numerical instability, known as two kinds of breakdowns, exactly similarly to those of the BCG method. We still call such two kinds of breakdowns as Lanczos breakdown and pivot breakdown, which will be analyzed and discussed in the following sections. As far as reported, a large number of strategies designed to handle such two kinds of breakdowns have been proposed in the literature [12–16], including composite step techniques [11, 17–21] for pivot breakdown and look-ahead techniques [22–30] for Lanczos breakdown or both types. For comprehensive reviews and discussions about the two kinds of breakdowns mentioned above, refer to [2, 6, 16] and the references therein. As shown and analyzed detailedly by Bank and Chan [18, 19], the 2×2 composite step strategy employed in the development of the composite step biconjugate gradient method (CSBCG) can eliminate (near) pivot breakdowns of the BCG method caused by (near) singularity of principal submatrices of the tridiagonal matrix generated by the underlying Lanczos biorthogonalization procedure, assuming that Lanczos breakdowns do not occur. Furthermore, besides simplifying implementations a great deal in comparison with those existing look-ahead techniques, the composite step strategy used there is able to provide some smoothing of the convergence history of the norm of the residuals without involving any arbitrary machine-dependent or user-supplied tolerances.This paper revisits the composite step strategy taken for the CSBCG method [18, 19] to naturally stabilize the convergence performance of the BiCOR method from the pivot-breakdown treatment point of view. Therefore, suppose that the underlying biconjugate A-orthonormalization procedure does not break down; that is to assume that there is no occurrence of the so-called Lanczos breakdown encountered during algorithm implementations throughout. The resulting interesting algorithm, given a name as the composite step BiCOR (CSBiCOR) method, could be also considered as a relatively simple modification of the regular BiCOR method.The main objectives are twofold. First, the CSBiCOR method is devised to be able to avoid near pivot breakdowns and compute all the well-defined BiCOR iterates stably with only minor modifications inspired by the advantages of the CSBCG method over the BCG method. Second, the CSBiCOR method can reduce the number of spikes in the convergence history of the norm of the residuals to the greatest extent, providing some further practically desired smoothing behavior towards stabilizing the behavior of the BiCOR method when it has erratic convergence behaviors. Additional purpose is that the CSBiCOR method inherits the promising advantages of the empirically observed stability and fast convergence rate of the BiCOR method over the BCG method so that it outperforms the CSBCG method to some extent.The remainder of this work is organized as follows. A brief review of the biconjugateA-orthonormalization procedure is made and the factorization of the resulted nonsingular tridiagonal matrices in the aforementioned procedure is considered in Section 2 for the preparation of the theoretical basis and relevant background of the CSBiCOR method. And then the breakdowns of the BiCOR method are presented in Section 3 while the CSBiCOR method is derived in Section 4. Section 5 analyzes the properties of the CSBiCOR method as well as providing a best approximation result for the convergence of the CSBiCOR method. In order to conveniently and simply present the CSBiCOR method, some implementation issues on both algorithm simplification and stepping strategy are discussed detailedly in Section 6. The improved performance of the CSBiCOR method will be illustrated in Section 7 in the perspective that the CSBiCOR method could hopefully smooth the convergence history through the reduction of the number of spikes when the BiCOR method has irregular convergence behavior. Finally, concluding remarks are made in the last section.Throughout the paper, we denote the overbar “–”  the conjugate complex of a scalar, vector or matrix and the superscript “T” the transpose of a vector, or matrix. For a non-Hermitian matrix A=(aij)N×N∈ℂN×N, the Hermitian conjugate of A is denoted as
(1)AH≡A¯T=(a¯ji)N×N.
The standard Hermitian inner product of two complex vectors u,v∈ℂN is defined as
(2)〈u,v〉=uHv=∑i=1N‍u¯ivi.
The nested Krylov subspace of dimension t generated by A from v is of the form
(3)𝒦t(A,v)=span{v,Av,A2v,…,At-1v}.
In addition, ei denotes the ith column of the appropriate identity matrix.When it will be helpful, we will use the word “ideally” (or “mathematically”) to refer to a result that could hold in exact arithmetic ignoring effects of rounding errors, and “numerically” (or “computationally”) to a result of a finite precision computation.


## 2. Theoretical Basis for the CSBiCOR Method

For the sake of discussing the theoretical basis and relevant background of the CSBiCOR method, the biconjugateA-orthonormalization procedure [7] is first briefly recalled as in Algorithm 1, which can ideally build up a pair of biconjugate A-orthonormal bases for the dual Krylov subspaces 𝒦m(A,v1) and 𝒦m(AH,w1), where v1 and w1 are chosen initially to satisfy certain conditions.Algorithm 1:Biconjugate A-orthonormalization procedure.
(1) Choosev1, ω1, such that 〈ω1,Av1〉=1
(2) Setβ1=δ1≡0, ω0=v0≡0∈ℂN
(3)for  j=1,2,…  do
(4)αj=〈ωj,A(Avj)〉
(5)v^j+1=Avj-αjvj-βjvj-1
(6)ω^j+1=AHωj-α-jωj-δjωj-1
(7)δj+1=|〈ω^j+1,Av^j+1〉|1/2
(8)βj+1=〈ω^j+1,Av^j+1〉/δj+1
(9)vj+1=v^j+1/δj+1
(10)ωj+1=ω^j+1/β-j+1
(11)end forObserve that the above algorithm is possible to have Lanczos-type breakdown wheneverδj+1 vanishes while w^j+1 and Av^j+1 are not equal to 0∈ℂN appearing in line 8. In the interest of counteraction against such breakdowns, refer oneself to remedies such as the so-called look-ahead strategies [22–30] which can enhance stability while increasing cost modestly. But that is outside the scope of this paper and we will not pursue that here. For more details, please refer to [2, 6] and the references therein. Throughout the rest of the present paper, suppose there is no such Lanczos-type breakdown encountered during algorithm implementations because most of our considerations concern the exploration of the composite step strategy [18, 19] to handle the pivot breakdown occurring in the BiCOR method for solving non-Hermitian linear systems.Next, some properties of the vectors produced by Algorithm1 are reviewed [7] in the following proposition for the preparation of the theoretical basis of the composite step method.Proposition 1.
If Algorithm1 proceeds m steps, then the right and left Lanczos-type vectors vj,j=1,2,…,m and wi,i=1,2,…,m form a biconjugate A-orthonormal system in exact arithmetic, that is,
(4)〈ωi,Avj〉=δi,j,1≤i,j≤m.
Furthermore, denote byVm=[v1,v2,…,vm] and Wm=[w1,w2,…,wm] the N×m matrices and by Tm_ the extended tridiagonal matrix of the form
(5)Tm_=[Tmδm+1emT],
where
(6)Tm=[α1β2δ2α2β3⋱⋱⋱δm-1αm-1βmδmαm],
whose entries are the coefficients generated during the algorithm implementation, and in which α1,…,αm,β2,…,βm are complex while δ2,…,δm are positive. Then with the biconjugate A-orthonormalization procedure, the following four relations hold:
(7)AVm=VmTm+δm+1vm+1emT,(8)AHWm=WmTmH+β-m+1ωm+1emT,(9)WmHAVm=Im,(10)WmHA2Vm=Tm.It is well known that the Lanczos biorthogonalization procedure can mathematically build up a pair of biorthogonal bases for the dual Krylov subspaces𝒦m(A,v1) and 𝒦m(AH,w1), where v1 and w1 are initially selected such that 〈ω1,v1〉=1 [6]. While the biconjugate A-orthonormalization procedure presented as in Algorithm 1 can ideally build up a pair of biconjugate A-orthonormal bases for the dual Krylov subspaces 𝒦m(A,v1) and 𝒦m(AH,w1), where v1 and w1 are chosen initially to satisfy 〈ω1,Av1〉=1. Keeping in mind the above differences in terms of different Krylov subspace bases constructed by the two different underlying procedures, we will in parallel borrow some of the results of the CSBCG method [18] for establishing the following counterparts for the CSBiCOR method; see [18] for relevant proofs.It should be emphasized that the present CSBiCOR method will be derived in exactly the same way the CSBCG method is obtained [19]. Hence, the analysis and theoretical background for the CSBiCOR method in this section are the counterparts of the CSBCG method in a different context, that is, the biconjugate A-orthonormalization procedure as presented in Algorithm 1. They are included as follows for the purpose of theoretical completeness and can be proved in the same way for the corresponding results of the CSBCG method analyzed in [18].The following proposition states the tridiagonal matrixTm formed in Proposition 1 cannot have two successive singular leading principal submatrices if the coefficient matrix A is nonsingular.Proposition 2.
LetTk(1≤k≤m) be the upper left principal submatrices of the nonsingular tridiagonal matrix Tm appeared as in Proposition 1. Then Tk-1 and Tk cannot be both singular.It is noted that there may exist possible breakdown in the factorization without pivoting of the tridiagonal matrixTm formed in Proposition 1. The following results illustrate how to correct such problem with the occasional use of 2×2 block pivots.Proposition 3.
LetTm be the nonsingular tridiagonal matrix formed in Proposition 1. Then Tm can be factored as
(11)Tm=LmDmUm,
where Lm is unit lower block bidiagonal, Um is unit upper block bidiagonal, and Dm is block diagonal, with 1×1 and 2×2 diagonal blocks.With the above two propositions, one can prove the following result insuring that if the upper left principal submatricesTk(1≤k≤m) of Tm is singular, it can still be factored.Corollary 4.
Suppose one upper left principal submatrixTk(1≤k≤m) of the nonsingular tridiagonal matrix Tm formed in Proposition 1 is singular, but Tk-1 must be nonsingular as stated in Proposition 2. Then
(12)Tk=LkDkUk,
where Lk is unit lower block bidiagonal, Uk is unit upper block bidiagonal, and Dk is block diagonal, with 1×1 and 2×2 diagonal blocks, and in particular, the last block of Dk is the 1×1 zero matrix.From Corollary4, the singularity of Tk(1≤k≤m) can be recognized by only examining the last diagonal element (and next potential pivot) dk. If dk=0, then the 2×2 block
(13)[dkβk+1δk+1αk+1]
will be nonsingular and can be used as a 2×2 pivot, where βk+1,δk+1, and αk+1 are the corresponding elements of Tm as shown in Proposition 1.It can be concluded that appropriate combinations of1×1 and 2×2 steps for the BiCOR method can skip over the breakdowns caused by the singular principal submatrices of the tridiagonal matrix Tm generated by the underlying biconjugate A-orthonormalization procedure presented as in Algorithm 1 since the BiCOR method pivots implicitly with those submatrices. The next section will have a detailed investigation of breakdowns possibly occurring in the BiCOR method.


## 3. Breakdowns of the BiCOR Method

Given an initial guessx0 to the non-Hermitian linear system Ax=b associated with the initial residual r0=b-Ax0, define a Krylov subspace ℒm≡AHspan(Wm)=AH𝒦m(AH,w1), where Wm is defined in Proposition 1, v1=r0/||r0||2 and w1 is chosen arbitrarily such that 〈w1,Av1〉≠0. But w1 is often chosen to be equal to Av1/||Av1||22 subjecting to 〈w1,Av1〉=1. It is worth noting that this choice for w1 plays a significant role in establishing the empirically observed superiority of the BiCOR method to the BiCR [31] method as well as to the BCG method [7]. Thus running Algorithm 1  m steps, we can seek an mth approximate solution xm from the affine subspace x0+𝒦m(A,v1) of dimension m, by imposing the Petrov-Galerkin condition
(14)b-Axm⊥ℒm,
which can be mathematically written in matrix formulation as
(15)(AHWm)H(b-Axm)=0.Analogously, anmth dual approximation xm* of the corresponding dual system AHx*=b* is sought from the affine subspace x0*+𝒦m(AH,w1) of dimension m by satisfying
(16)b*-AHxm*⊥A𝒦m(A,v1),
which can be mathematically written in matrix formulation as
(17)(AVm)H(b*-AHxm)=0,
where x0* is an initial dual approximate solution and Vm is defined in Proposition 1 with v1=r0/||r0||2.Consequently, the BiCOR iteratesxj’s can be computed by the coming Algorithm 2, which is just the unpreconditioned BiCOR method with the preconditioner M there taken as the identity matrix [7] and has been rewritten with the algorithmic scheme of the unpreconditioned BCG method as presented in [6, 19].Algorithm 2:Algorithm BiCOR.
(1) Computer0=b-Ax0 for some initial guess x0.
(2) Chooser0*=P(A)r0 such that 〈r0*,Ar0〉≠0, where P(t) is a polynomial in t.
(e.g.,r0*=Ar0).
(3) Setp0=r0, p0*=r0*, q0=Ap0, q0*=AHp0*, r^0=Ar0, ρ0=〈r0*,r^0〉.
(4)for  n=0,1,…  do
(5)σn=〈qn*,qn〉
(6)αn=ρn/σn
(7)xn+1=xn+αnpn
(8)rn+1=rn-αnqn
(9)xn+1*=xn*+α-npn*
(10)rn+1*=rn*-α-nqn*
(11)r^n+1=Arn+1
(12)ρn+1=〈rn+1*,r^n+1〉
(13)if  ρn+1=0, method fails
(14)βn+1=ρn+1/ρn
(15)pn+1=rn+1+βn+1pn
(16)pn+1*=rn+1*+β-n+1pn*
(17)qn+1=r^n+1+βn+1qn
(18)qn+1*=AHpn*
(19)check convergence; continue if necessary
(20)end forSuppose Algorithm2 runs successfully to step n, that is, σi≠0, ρi≠0, i=0,1,…,n-1. The BiCOR iterates satisfy the following properties [7].Proposition 5.
LetRn+1=[r0,r1,…,rn], Rn+1*=[r0*,r1*,…,rn*] and Pn+1=[p0,p1,…,pn], Pn+1*=[p0*,p1*,…,pn*]. One has the following.(1)
Range(Rn+1)=Range(Pn+1)=𝒦n+1(A,r0), Range(Rn+1*)=Range(Pn+1*)=𝒦n+1(AH,r0*).(2)
Rn+1*HARn+1is diagonal.(3)
Pn+1*HA2Pn+1is diagonal.Similarly to the breakdowns of the BCG method [19], it is observed from Algorithm 2 that there also exist two possible kinds of breakdowns for the BiCOR method:(1)
ρn≡〈rn*,r^n〉≡〈rn*,Arn〉=0 but rn* and Arn are not equal to 0∈ℂN appearing in line 14;(2)
σn≡〈qn*,qn〉≡〈AHpn*,Apn〉=0 appearing in line 6.Although the computational formulae for the quantities where the breakdowns reside are different between the BiCOR method and the BCG method, we do not have a better name for them. And therefore, we still call the two cases of breakdowns described above as Lanczos breakdown and pivot breakdown, respectively.The Lanczos breakdown can be cured using look-ahead techniques [22–30] as mentioned in the first section, but such techniques require a careful and sophisticated way so as to make them become necessarily quite complicated to apply. This aspect of applying look-ahead techniques to the BiCOR method demands further research.In this paper, we attempt to resort to the composite step idea employed for the CSBCG method [18, 19] to handle the pivot breakdown of the BiCOR method with the assumption that the underlying biconjugate A-orthonormalization procedure depicted as in Algorithm 1 does not break down; that is the situation where σn=0 while ρn≠0.


## 4. The Composite Step BiCOR Method

Suppose Algorithm2 comes across a situation where σn=0 after successful algorithm implementation up to step n with the assumption that ρn≠0, which indicates that the updates of xn+1,rn+1,xn+1*,rn+1* are not well defined. Taking the composite step idea, we will avoid division by σn=0 via skipping this (n+1)th update and exploiting a composite step update to directly obtain the quantities in step (n+2) with scaled versions of rn+1 and rn+1* as well as with the previous primary search direction vector pn and shadow search direction vector pn*. The following process for deriving the CSBiCOR method is the same as that of the derivation of the CSBCG method [19] except for the different underlying procedures involved to correspondingly generate different Krylov subspace bases.Analogously, define auxiliary vectorszn+1∈𝒦n+2(A,r0) and zn+1*∈𝒦n+2(AH,r0*) as follows:
(18)zn+1=σnrn+1=σnrn-ρnApn,(19)zn+1*=σ-nrn+1*=σ-nrn*-ρ-nAHpn*,
which are then used to look for the iterates xn+2∈x0+𝒦n+2(A,r0) and xn+2*∈x0*+𝒦n+2(AH,r0*) in step (n+2) as follows:
(20)xn+2=xn+[pn,zn+1]fn,xn+2*=xn*+[pn*,zn+1*]fn*,
where fn,fn*∈ℂ2. Correspondingly, the (n+2)th primary residual rn+2∈𝒦n+3(A,r0) and shadow residual rn+2*∈𝒦n+3(AH,r0*) are, respectively, computed as
(21)rn+2=rn-A[pn,zn+1]fn,(22)rn+2*=rn*-AH[pn*,zn+1*]fn*.
The biconjugate A-orthogonality condition between the BiCOR primary residuals and shadow residuals shown as Property (2) in Proposition 5 requires
(23)〈[pn*,zn+1*],Arn+2〉=0,〈[pn,zn+1],AHrn+2*〉=0,
combining with (21) and (22) gives rise to the following two 2×2 systems of linear equations for, respectively, solving fn and fn*:(24)[〈AHpn*,Apn〉〈AHpn*,Azn+1〉〈AHzn+1*,Apn〉〈AHzn+1*,Azn+1〉][fn(1)fn(2)]=[〈pn*,Arn〉〈zn+1*,Arn〉],(25)[〈Apn,AHpn*〉〈Apn,AHzn+1*〉〈Azn+1,AHpn*〉〈Azn+1,AHzn+1*〉][fn*(1)fn*(2)]=[〈Apn,rn*〉〈Azn+1,rn*〉].Similarly, the(n+2)th primary search direction vector pn+2∈𝒦n+3(A,r0) and shadow search direction vector pn+2*∈𝒦n+3(AH,r0*) in a composite step are computed with the following form:
(26)pn+2=rn+2+[pn,zn+1]gn,(27)pn+2*=rn+2*+[pn*,zn+1*]gn*,
where gn,gn*∈ℂ2.The biconjugateA2-orthogonality condition between the BiCOR primary search direction vectors and shadow search direction vectors shown as Property (3) in Proposition 5 requires
(28)〈[pn*,zn+1*],A2pn+2〉=0,〈[pn,zn+1],(AH)2pn+2*〉=0,
combining with (26) and (27) results in the following two 2×2 systems of linear equations for respectively solving gn and gn*:(29)[〈AHpn*,Apn〉〈AHpn*,Azn+1〉〈AHzn+1*,Apn〉〈AHzn+1*,Azn+1〉][gn(1)gn(2)]=-[〈AHpn*,Arn+2〉〈AHzn+1*,Arn+2〉],(30)[〈Apn,AHpn*〉〈Apn,AHzn+1*〉〈Azn+1,AHpn*〉〈Azn+1,AHzn+1*〉][gn*(1)gn*(2)]=-[〈Apn,AHrn+2*〉〈Azn+1,AHrn+2*〉].Therefore, it could be able to advance from stepn to step (n+2) to provide xn+2, rn+2, xn+2*, rn+2*, pn+2, pn+2* by solving the above four 2×2 linear systems represented as in (24), (25), (29), and (30). With an appropriate combination of 1×1 and 2×2 steps, the CSBiCOR method can be simply obtained with only a minor modification to the usual implementation of the BiCOR method. Before explicitly presenting the algorithm, some properties of the CSBiCOR method will be given in the following section.


## 5. Some Properties of the CSBiCOR Method

In order to show the use of2×2 steps is sufficient for CSBiCOR method to compute exactly all those well-defined BiCOR iterates stably provided that the underlying biconjugate A-orthonormalization procedure depicted as in Algorithm 1 does not break down, it is necessary to establish some lemmas similar to those for the CSBCG method [19].Lemma 6.
Using either a1×1 or a 2×2 step to obtain pn, pn*, the following hold:
(31)〈AHrn*,Apn〉=〈AHpn*,Arn〉=〈AHpn*,Apn〉=〈qn*,qn〉≡σn.Proof.
Ifpn,pn* are obtained via a 1×1 step, that is, pn=rn+βnpn-1, pn*=rn*+β-npn-1* as shown in lines 15 and 16 of Algorithm 2, then with the biconjugate A2-orthogonality condition between the BiCOR primary search direction vectors and shadow search direction vectors shown as Property (3) in Proposition 5, it follows that
(32)〈AHrn*,Apn〉=〈AH(pn*-β-npn-1*),Apn〉=〈AHpn*,Apn〉-βn〈AHpn-1*,Apn〉=〈AHpn*,Apn〉=〈qn*,qn〉≡σn,〈AHpn*,Arn〉=〈AHpn*,A(pn-βnpn-1)〉=〈AHpn*,Apn〉-βn〈AHpn*,Apn-1〉=〈AHpn*,Apn〉=〈qn*,qn〉≡σn.
Ifpn,pn* are obtained via a 2×2 step, that is, pn=rn+gn-2(1)pn-2+gn-2(2)zn-1, pn*=rn*+gn-2*(1)pn-2*+gn-2*(2)zn-1* as presented in (26) and (27), then according to the biconjugate A2-orthogonality conditions mentioned above for deriving pn,pn* with a 2×2 step, it follows that
(33)〈AHrn*,Apn〉=〈AH(pn*-gn-2*(1)pn-2*-gn-2*(2)zn-1*),Apn〉=〈AHpn*,Apn〉-gn-2*(1)¯〈AHpn-2*,Apn〉-gn-2*(2)¯〈AHzn-1*,Apn〉=〈AHpn*,Apn〉=〈qn*,qn〉≡σn,〈AHpn*,Arn〉=〈AHpn*,A(pn-gn-2(1)pn-2-gn-2(2)zn-1)〉=〈AHpn*,Apn〉-gn-2(1)〈AHpn*,Apn-2〉-gn-2(2)〈AHpn*,Azn-1〉=〈AHpn*,Apn〉=〈qn*,qn〉≡σn.With the above lemma, we can show the relationships among the four2×2 coefficient matrices of the linear systems represented as in (24), (25), (29), and (30).Lemma 7.
The four2×2 coefficient matrices of the linear systems represented as in (24), (25), (29), and (30) have the following properties and relationships.(1)
The coefficient matrices in (24) and (29) are identical, so are those matrices in (25) and (30).(2)
All the coefficient matrices involved are symmetric.(3)
The coefficient matrices in (24) and (25) (correspondingly in (29) and (30)) are Hermitian to each other.Proof.
The first relation is obvious by observation of the corresponding coefficient matrices in (24) and (29) (correspondingly in (25) and (30)).
For proving the second property, it suffices to show〈AHzn+1*,Apn〉=〈AHpn*,Azn+1〉. From the presentations of zn+1 and zn+1* as respectively defined in (18) and (19), and following from Lemma 6, we have
(34)〈AHzn+1*,Apn〉=〈AH(σ-nrn*-ρ-nAHpn*),Apn〉=σn〈AHrn*,Apn〉-ρn〈(AH)2pn*,Apn〉=σn2-ρn〈(AH)2pn*,Apn〉,〈AHpn*,Azn+1〉=〈AHpn*,A(σnrn-ρnApn)〉=σn〈AHpn*,Arn〉-ρn〈AHpn*,A2pn〉=σn2-ρn〈(AH)2pn*,Apn〉.
The third fact directly comes from the Hermitian property of the standard Hermitian inner product of two complex vectors defined at the end of the first section; see, for instance, [6].It is noticed that similar Hermitian relationships also hold for the corresponding matriced involved in (1), (2), (3) and (4) in [19] when the CSBCG method is applied in complex cases.Now, an alternative result stating that a2×2 step is always sufficient to skip over the pivot breakdown of the BiCOR method when σn=0, just as stated at the end of Section 2, can be shown in the next lemma.Lemma 8.
Suppose that the biconjugateA-orthonormalization procedure depicted as in Algorithm 1 underlying the BiCOR method does not breakdown; that is the situation where ρi≠0, i=0,1,…,n and 〈zn+1*,Azn+1〉≠0. If σn≡〈qn*,qn〉≡〈AHpn*,Apn〉=0, then the 2×2 coefficient matrices in (24), (25), (29), and (30) are nonsingular.Proof.
It suffices to show that ifσn≡〈qn*,qn〉≡〈AHpn*,Apn〉=0, then 〈AHzn+1*,Apn〉≠0. From Lemma 7 and (18), we have
(35)〈AHzn+1*,Apn〉=〈AHpn*,Azn+1〉,〈AHzn+1*,Apn〉=〈AHzn+1*,-zn+1ρn〉=-1ρn〈zn+1*,Azn+1〉≠0.The next proposition shows some properties of the CSBiCOR iterates similar to those of the BiCOR method presented in Proposition5 with some minor changes.Proposition 9.
For the CSBiCOR algorithm, letRn+1=[r0,r1,…,rn], Rn+1*=[r0*,r1*,…,rn*], where ri, ri* are replaced by zi, zi* appropriately if the ith step is a composite 2×2 step; and let Pn+1=[p0,p1,…,pn], Pn+1*=[p0*,p1*,…,pn*]. Assuming the underlying biconjugate A-orthonormalization procedure does not break down, one has the following properties.(1)
Range(Rn+1)=Range(Pn+1)=𝒦n+1(A,r0),Range(Rn+1*)=Range(Pn+1*)=𝒦n+1(AH,r0*).(2)
Rn+1*HARn+1is diagonal.(3)
Pn+1*HA2Pn+1=diag(Dk), whereDkis either of order1×1or2×2and the sum of the dimensions ofDk’s is(n+1).Proof.
The properties can be proved with the same frame for those of Theorem  4.4 in [19].Therefore, we can show with the above properties that the use of2×2 steps is sufficient for the CSBiCOR method to compute exactly all those well-defined BiCOR iterates stably provided that the underlying biconjugate A-orthonormalization procedure depicted as in Algorithm 1 does not break down.Proposition 10.
The CSBiCOR method is able to compute exactly all those well-defined BiCOR iterates stably without pivot breakdown assuming that the underlying biconjugateA-orthonormalization procedure does not break down.Proof.
By construction,xn+2=xn+fn(1)pn+fn(2)zn+1; by induction, xn∈x0+𝒦n(A,r0) and pn∈𝒦n+1(A,r0); by definition, zn+1∈𝒦n+2(A,r0). Then it follows that xn+2∈x0+𝒦n+2(A,r0). From Proposition 9, we have rn+2⊥AH𝒦n+2(AH,r0*). Therefore, xn+2 exactly satisfies the Petrov-Galerkin condition defining the BiCOR iterates.Before ending this section, we present a best approximation result for the CSBiCOR method, which uses the same framework as that for the CSBCG method [19]. For discussion about convergence results for Krylov subspaces methods, refer to [2, 6, 32, 33] and the references therein. To the best of our knowledge, this is the first attempt to study the convergence result for the BiCOR method [7].First let𝒱m=span(v1,v2,…,vm) and 𝒲m=span(w1,w2,…,wm) denote the Krylov subspaces generated by the biconjugate A-orthonormalization procedure after running Algorithm 1  m steps. The norms associated with the spaces 𝒱N≡ℂN and 𝒲N≡ℂN are defined as
(36)|∥v∥|r2=vHMrv,|∥w∥|l2=wHMlw,
where, v,w∈ℂN while Mr and Ml are symmetric positive definite matrices. Then we have the following proposition with the initial guess x0 of the linear system being 0∈ℂN. The case with a nonzero initial guess can be treated adaptively.Proposition 11.
Suppose that for allv∈𝒱N and for all w∈𝒲N, one has
(37)|wHA2v|≤Γ|∥v∥|r|∥w∥|l,
where, Γ is a constant independent of v and w. Furthermore, assuming that for those steps in the composite step Biconjugate A-Orthogonal Residual (CSBiCOR) method in which we compute an approximation xm, we have
(38)infv∈𝒱m|∥v∥|r=1supw∈𝒲m|∥w∥|l≤1wHA2v≥γk≥γ>0.
Then,
(39)|∥x-xm∥|r≤(1+Γγ)infv∈𝒱m|∥x-v∥|r.Proof.
From the Petrov-Galerkin condition (15), we have for w∈𝒲m,
(40)(AHw)H(Ax-Axm)=wHA2(x-xm)=0,
yielding with arbitrary v∈𝒱m(41)wHA2(xm-v)=wHA2(x-v).
It is noted that xm-v∈𝒱m because of the assumption of the zero initial guess. Then taking the sup of both sides of the above equation for all |∥w∥|l≤1 with the inf-sup condition (38) to bound the left-hand side and the continuity assumption (37) to bound the right-hand side results in
(42)γk|∥xm-v∥|r≤Γ|∥x-v∥|r|∥w∥|l≤Γ|∥x-v∥|r,
which yields
(43)|∥x-xm∥|r≤|∥x-v∥|r+|∥xm-v∥|r≤(1+Γγ)|∥x-v∥|r
with the triangle inequality. The final assertion (39) follows immediately since v∈𝒱m is arbitrary.


## 6. Implementation Issues

For the purpose of conveniently and simply presenting the CSBiCOR method, we first prove some relations simplifying the solutions of the four2×2 linear systems represented as in (24), (25), (29), and (30).Lemma 12.
For the four2×2 linear systems represented as in (24), (25), (29), and (30), the following relations hold:(1)
〈pn*,Arn〉=〈Apn,rn*〉¯=〈rn*,Arn〉=〈rn*,r^n〉≡ρn and 〈zn+1*,Arn〉=〈Azn+1,rn*〉=0,fn=f-n*.(2)
〈AHpn*,Arn+2〉=〈Apn,AHrn+2*〉=0 and 〈AHzn+1*,Arn+2〉=〈Azn+1,AHrn+2*〉¯=-ρn+2/fn(2), where ρn+2≡〈rn+2*,Arn+2〉,gn=g-n*.(3)
〈AHpn*,Azn+1〉=〈Apn,AHzn+1*〉¯=-θn+1/ρn, where θn+1≡〈zn+1*,Azn+1〉.Proof.
(1) If a1×1 step is taken, then pn*=rn*+β-npn-1*, giving that 〈pn*,Arn〉=〈rn*+β-npn-1*,Arn〉=〈rn*,Arn〉+βn〈pn-1*,Arn〉=〈rn*,Arn〉 because of the last-term vanishment due to Proposition 9. Analogously, the iteration pn=rn+βnpn-1 similarly gives 〈Apn,rn*〉=〈Arn+βnApn-1,rn*〉=〈Arn,rn*〉+β-n〈Apn-1,rn*〉=〈Arn,rn*〉=〈rn*,Arn〉¯ from Proposition 9 and the Hermitian property of the standard Hermitian inner product. If a 2×2 step is taken, then by construction as shown in (27) and Proposition 9, we have 〈pn*,Arn〉 = 〈rn*+gn-2*(1)pn-2*+gn-2*(2)zn-1*,Arn〉 = 〈rn*,Arn〉+gn-2*(1)¯〈pn-2*,Arn〉+gn-2*(2)¯〈zn-1*,Arn〉 = 〈rn*,Arn〉. Similarly, with Proposition 9 and the Hermitian property of the standard Hermitian inner product, it can be shown that 〈Apn,rn*〉= 〈A(rn+gn-2(1)pn-2+gn-2(2)zn-1),rn*〉 = 〈Arn,rn*〉+gn-2(1)¯〈Apn-2,rn*〉+gn-2(2)¯〈Azn-1,rn*〉 = 〈Arn,rn*〉 = 〈rn*,Arn〉¯. The equalities 〈zn+1*,Arn〉=〈Azn+1,rn*〉=0 directly follow from Proposition 9. Finally, because the coefficient matrices of the symmetric systems in (24) and (25) respectively governing fn and fn* are Hermitian to each other as shown in Lemma 7 while the corresponding right-hand sides are conjugate to each other as shown here, we can deduce that fn=f-n*.
(2) From Proposition9, it can directly verify that 〈AHpn*,Arn+2〉=〈Apn,AHrn+2*〉=0 with the fact that Apn∈𝒦n+2(A,r0) and AHpn*∈𝒦n+2(AH,r0*). Next, it can be noted that if σn=0, then fn(2)≠0. Otherwise, the first equation in the linear system (24) cannot hold for 〈AHpn*,Azn+1〉≠0 as proved in Lemma 8 and 〈pn*,Arn〉=ρn≠0 as shown in the above (1). Then by construction of (22) and Proposition 9, we have 〈AHzn+1*,Arn+2〉 = 〈(1/fn*(2))(rn*-rn+2*-fn*(1)AHpn*),Arn+2〉 = (1/fn*(2)¯)(〈rn*,Arn+2〉-〈rn+2*,Arn+2〉-fn*(1)¯〈AHpn*,Arn+2〉) = (-1/fn*(2)¯)〈rn+2*,Arn+2〉 = (-1/fn*(2)¯)ρn+2 = (-1/fn(2))ρn+2 with the fact that fn=f-n* proved in the above (1). By construction of (21) and Proposition 9, we can also show that 〈Azn+1,AHrn+2*〉=-ρn+2/fn(2)¯. Similarly to the relationships between fn and fn*, we can deduce that gn=g-n*, because the coefficient matrices of the symmetric systems in (29) and (30) respectively governing gn and gn* are Hermitian to each other as shown in Lemma 7 while the corresponding right-hand sides are conjugate to each other as shown here.
(3) By construction of (19) and Proposition 9, it follows that 〈AHpn*,Azn+1〉 = 〈(1/ρ-n)(σ-nrn*-zn+1*),Azn+1〉 = (1/ρn)(σn〈rn*,Azn+1〉-〈zn+1*,Azn+1〉) = (-1/ρn)〈zn+1*,Azn+1〉 = -θn+1/ρn, where θn+1≡〈zn+1*,Azn+1〉. Because it is already known that 〈AHzn+1*,Apn〉 = 〈AHpn*,Azn+1〉 as proved in Lemma 7, we have 〈Apn,AHzn+1*〉 = 〈AHzn+1*,Apn〉¯ = 〈AHpn*,Azn+1〉¯ = -θn+1/ρn¯. The proof is completed.As a result of Lemmas6 and 12, when a 2×2 step is taken, the solutions for fn and gn involved in the four linear systems represented as in (24), (25), (29), and (30) can be simplified only by solving the following two linear systems:
(44)[σn-θn+1ρn-θn+1ρnζn+1][fn(1)fn(2)]=[ρn0],[σn-θn+1ρn-θn+1ρnζn+1][gn(1)gn(2)]=[0ρn+2fn(2)],
where ζn+1=〈AHzn+1*,Azn+1〉.It should be emphasized that these above two systems have the same representations for the CSBCG method [19] but differ in detailed computational formulae for the corresponding quantities because of the different underlying procedures mentioned in Section 1. As a result, fn and gn can be explicitly represented as
(45)fn=(ζn+1ρn,θn+1)ρn2ηn,gn=(ρn+2ρn,σnρn+2θn+1),
where ηn≡σnζn+1ρn2-θn+12.Combining all the recurrences discussed above for either a1×1 step or a 2×2 step and taking the strategy of reducing the number of matrix-vector multiplications by introducing an auxiliary vector recurrence and changing variables adopted for the BiCOR method [7] as well as for the BiCR method [31], together lead to the CSBiCOR method. The pseudocode for the preconditioned CSBiCOR with a left preconditioner B can be represented by Algorithm 3, where sn+1 and zn+1 are used to respectively denote the unpreconditioned and the preconditioned residuals. It is obvious to note that the number of matrix-vector multiplications in this algorithm remains the same per step as in the BiCOR algorithm [7]. Therefore, the cost of the use of 2×2 composite steps is negligible. That is, 2×2 composite steps cost approximately twice as much as 1×1 steps just as stated for the CSBCG method [18]. Thus, from the point of view of choosing 1×1 or 2×2 steps, there is no significant difference with respect to algorithm cost. However, the presented algorithm is with the motivation of obtaining smoother and, hopefully, faster convergence behavior in comparison with the BiCOR method besides eliminating pivot breakdowns.Algorithm 3:Left preconditioned CSBiCOR method.
(1) Computer0=b-Ax0 for some initial guess x0.
(2) Chooser0*=P(A)r0 such that 〈r0*,Ar0〉≠0, where P(t) is a polynomial in t.
(e.g.,r0*=Ar0). Set p0=r0, p~0=r~0, q0=Ap0, q~0=AHp~0.
(3) Computeρ0=〈r~0,Ar0〉.
(4) Begin LOOP(n=0,1,2,…)
(5)σn=〈q~n,qn〉
(6)sn+1=σnrn-ρnqn
(7)s~n+1=σ-nr~n-ρ-nq~n
(8)yn+1=Asn+1
(9)y~n+1=AHs~n+1
(10)θn+1=〈s~n+1,yn+1〉
(11)ζn+1=〈y~n+1,yn+1〉
(12)if 1 × 1 step then
(13)αn=ρn/σn
(14)ρn+1=θn+1/σn2
(15)βn+1=ρn+1/ρn
(16)xn+1=xn+αnpn
(17)rn+1=rn-αnqn
(18)r~n+1=r~n-α-nq~n
(19)pn+1=sn+1/σn+βn+1pn
(20)qn+1=yn+1/σn+βn+1qn
(21)q~n+1=y~n+1/σ-n+β-n+1q~n
(22)n←n+1
(23)else
(24)δn=σnζn+1ρn2-θn+12
(25)αn=ζn+1ρn3/δn
(26)αn+1=θn+1ρn2/δn
(27)xn+2=xn+αnpn+αn+1sn+1
(28)rn+2=rn-αnqn-αn+1yn+1
(29)r~n+2=r~n-α-nq~n-α-n+1y~n+1
(30)solve  Bzn+2=rn+2
(31)solve  BHz~n+2=r~n+2
(32)z~n+2=Azn+2
(33)z~^n+2=A^Hz~n+2
(34)ρn+2=〈z~^n+2,rn+2〉
(35)βn+1=ρn+2/ρn
(36)βn+2=ρn+2σn/θn+1
(37)pn+2=zn+2+βn+1pn+βn+2sn+1
(38)qn+2=z^n+2+βn+1qn+βn+2yn+1
(39)q~n+2=z~^n+2+β-n+1q~n+β-n+2y~n+1
(40)n←n+2
(41)end if
(42) Check convergence; continue if necessary
(43) End LOOPFor the issue of deciding between1×1 and 2×2 updates, the heuristic based on the magnitudes of the residuals developed for the CSBCG method [18] is borrowed for the CSBiCOR method in order to choose the step size which maximizes numerical stability as well as to avoid overflow. The principle is to avoid “spikes” in the convergence history of the norm of the residuals. A 2×2 update will be chosen if the following circumstance satisfies
(46)∥rn+1∥>max{∥rn∥,∥rn+2∥}.For completeness of algorithm implementation, we recall the test procedure as follows. Defineπn+2≡ηnrn+2=ηnrn-ζn+1ρn3Apn-θn+1ρn2Azn+1. Then with the scaled versions of the corresponding quantities, the tests ∥rn+1∥≤∥rn∥ and ∥rn+2∥≤∥rn+1∥ can be, respectively, replaced by ∥zn+1∥≤|σn|∥rn∥ and |σn|∥πn+2∥≤|ηn|∥zn+1∥. Consequently, the test can be implemented with the code fragment shown in Algorithm 4.Algorithm 4
If||zn+1||≤|σn|||rn||, Then
1×1 Step
Else
||πn+2||=||ηnrn-ζn+1ρn3qn-θn+1ρn2yn+1||
If|σn|||πn+2||≤|ηn|||zn+1||, Then
2×2 Step
Else
1×1 Step
End If
End IfAnalogously to the effect obtained by the CSBCG method [19] in circumstances where (46) is satisfied, the use of 2×2 updates for the CSBiCOR method is also able to cut off spikes in the convergence history of the norm of the residuals possibly caused by taking two 1×1 steps in such circumstances. The resulting CSBiCOR algorithm inherits the promising properties of the CSBCG method [18, 19]. That is, the composite strategy employed can eliminate near pivot breakdowns while can provide some smoothing of the convergence history of the norm of the residuals without involving any arbitrary machine-dependent or user-supplied tolerances. However, it should be emphasized that the CSBiCOR method could only eliminate those spikes due to small pivots in the upper left principal submatrices Tk(1≤k≤m) of Tm formed in Proposition 1 but not all spikes. So that the residual norm of the CSBiCOR method does not decrease monotonically. By the way, learning from the mathematically equivalent but numerically different variants of the CSBCG method [18], we could also develop such kinds of variants of the CSBiCOR method, which will be taken into further account.


## 7. Examples and Numerical Experiments

This section mainly focuses on the analysis of different numerical effects of the composite step strategy employed into the BiCOR method with, far from being exhaustive, a few typical circumstances of test problems as arising from electromagnetics, discretizations of 2D/3D physical domains, and acoustics, which are described in Table1. All of them are borrowed in the MATLAB format from the University of Florida Sparse Matrix Collection provided by Davis [34], in which the meanings of the column headers of Table 1 can be found. The effect analysis part may provide some suggestions that when to make use of the composite step strategy should make significant progress towards an exact solution according to the convergence history of the residual norms. It is with the hope that the stabilizing effect of the CSBiCOR method could make the residual norm plot become smoother and, hopefully, faster decreasing since far outlying iterates and residuals are avoided in the smoothed sequences.Table 1
Structures of the three sets of test problems.
Ex.
Group and name
id
# rows
#cols
Nonzeros
Sym
1
Dehghani/light_in_tissue
1873
29,282
29,282
406,084
0%
2
Kim/kim1
862
38,415
38,415
933,195
0%
3
HB/young1c
278
841
841
4,089
85%In a realistic setting, one would use a preconditioner. With appropriate preconditioning techniques, such as those existing well-established preconditioning methodologies [35–37] based on approximate inverses, all the following involved methods for numerical comparison are very attractive for solving relevant classes of non-Hermitian linear systems. But this is not the point we pursue here. All the experiments are performed without preconditioning techniques in default if without other clarification. That is, the preconditioner B in Algorithm 3 will be taken as the identity matrix except otherwise clarified. Refer to the survey by Benzi [38] and the book by Saad [6] on preconditioning techniques for improving the performance and reliability of Krylov subspace methods.For all these test problems below, the BiCOR method will be implemented using the same code as for the CSBiCOR method with the pivot test being modified to always choose 1-step update steps for the purpose of conveniently showing the stabilizing effect of the composite step strategy on the BiCOR method. So does for the BCG method as implementated in [18, 19].The experiments have been carried out with machine precision10-16 in double precision floating point arithmetic in MATLAB 7.0.4 with a PC-Pentium (R) D CPU 3.00 GHz, 1 GB of RAM. We make comparisons in four aspects: number of iterations (referred to as Iters), CPU consuming time in seconds (referred to as CPU), log10 of the updated and final true relative residual 2-norms defined, respectively, as log10∥rn∥2/∥r0∥2 and log10∥b-Axn∥2/∥r0∥2 (referred to as Relresand TRR). Iters takes the form of “*/*” recording the total number of iteration steps and the number of 2×2 iteration steps involved in the corresponding composite step methods. Numerical results in terms of Iters,CPU and TRR are reported by means of tables while convergence histories involved are shown in figures with Iters (on the horizontal axis) versus Relres (on the vertical axis). The stopping criterion used here is that the 2-norm of the residual be reduced by a factor (referred to as TOL) of the 2-norm of the initial residual, that is, ∥rn∥2/∥r0∥2<TOL, or when Iters exceeded the maximal iteration number (referred to as MAXIT). Here, we take MAXIT=500. All these tests are started with an initial guess equal to 0∈ℂn. Whenever the considered problem contains no right-hand side to the original linear system Ax=b, let b=Ae, where e is the n×1 vector whose elements are all equal to unity, such that x=(1,1,…,1)T is the exact solution. A symbol “*” is used to indicate that the method did not meet the required TOL before MAXIT or did not converge at all.


### 7.1. Example 1: Dehghani/Light_in_Tissue

The first example is one in which the BCG method and the BiCOR method have almost quite smooth convergence behaviors and as few as near breakdowns. The CSBiCOR method has exactly the same convergence behavior as the BiCOR method by investigating Table2 and Figure 1(b). Under such circumstances when the convergence history of the norm of the residuals is smooth enough so that the composite step strategy does not gain a lot for the CSBiCOR method, the BiCOR method is much preferred and there is no need to employ the composite step strategy because the MATLAB implementation may be probably a little penalizing for the CSBiCOR code with respect to CPU. By the way, from Figure 1(c), the CSBCG method works as predicted in [18, 19], smoothing a little the convergence history of the BCG method. In particular, the CSBCG method computes a subset of the BCG iterates by clipping three “spikes” in the BCG history by means of the composite 2×2 steps from the particular investigations given in the bottom first two plots of Figures 1(d), 1(e) and Table 2.Table 2
Comparison results of Example 1 withTOL=10-6.
Method
BCG
BiCOR
CSBCG
CSBiCOR
Iters
330
318
327/3
318/0
CPU
40.7807
39.0256
40.1300
39.2116
TRR
−6.0318
−6.0150
−6.0318
−6.0150Convergence histories of Example 1 withTOL=10-6.
(a)
(b)
(c)
(d)
(e)
(f)Another interesting observation is that the good properties of the BiCOR method in terms of fast convergence speed and smooth convergence behavior in comparison with the BCG method are instinctually inherited by the CSBiCOR method so that the CSBiCOR method outperforms the CSBCG method in aspects of Iters and CPU as well as smooth effect, as shown in Table2 and Figure 1(f).


### 7.2. Example 2: Kim/Kim1

This is a good test example to vividly demonstrate the efficacy of the composite step strategy as depicted in Figure2. The BCG method, as observed in Figure 2(a), takes on many local “spikes” in the convergence curve. In such a case, the CSBCG method does a wonderful attempt to smooth the convergence history of the norm of the residuals as shown in Figure 2(c), although not leading to a monotonic decreasing convergence behavior. So does the CSBiCOR method to the BiCOR method by observation of Figure 2(d). Moreover, the CSBiCOR method seems to take less CPU than the BiCOR method while keeping approximately the same TRR when the BiCOR method is implemented using the same code as for the CSBiCOR method and the pivot test is modified to always choose 1-step update steps. Finally, the CSBiCOR method is again shown to be competitive to the CSBCG method as seen in Table 3 and Figure 2(b).Table 3
Comparison results of Example 2 withTOL=10-6.
Method
BCG
BiCOR
CSBCG
CSBiCOR
Iters
248
182
158/90
129/56
CPU
53.8496
39.2517
43.3087
33.8528
TRR
−6.3680
−6.0909
−6.3680
−6.0483Convergence histories of Example 2 withTOL=10-6.
(a)
(b)
(c)
(d)


### 7.3. Example 3: HB/Young1c

In the third example, the BCG and BiCOR methods have small irregularities and do not converge superlinearly, as reflected in Figure3(a). The two methods above seem to have almost the same “asymptotical” speed of convergence while the BiCOR method seems a little smoother than the BCG method. From Figures 3(c) and 3(d), the composite step strategy seems to play a surprisingly good stabilizing effect to the BCG and BiCOR method. Furthermore, the CSBCG and CSBiCOR methods take less CPU, respectively, than the BCG and BiCOR methods as reported in Table 4. It is stressed that although the CSBiCOR method consumes a little more Iters and CPU than the CSBCG method, it presents a little more smooth convergence behavior than the latter method as displayed in Figure 3(b). This is still thanks to the promising advantages of the empirically observed stability and fast convergence rate of the BiCOR method over the BCG method.Table 4
Comparison results of Example 3 withTOL=10-6.
Method
BCG
BiCOR
CSBCG
CSBiCOR
Iters
210
208
158/52
171/41
CPU
0.6848
0.5974
0.5372
0.5573
TRR
−6.1155
−6.0984
−6.0389
−6.0017Convergence histories of Example 3 withTOL=10-6.
(a)
(b)
(c)
(d)


## 7.1. Example 1: Dehghani/Light_in_Tissue

The first example is one in which the BCG method and the BiCOR method have almost quite smooth convergence behaviors and as few as near breakdowns. The CSBiCOR method has exactly the same convergence behavior as the BiCOR method by investigating Table2 and Figure 1(b). Under such circumstances when the convergence history of the norm of the residuals is smooth enough so that the composite step strategy does not gain a lot for the CSBiCOR method, the BiCOR method is much preferred and there is no need to employ the composite step strategy because the MATLAB implementation may be probably a little penalizing for the CSBiCOR code with respect to CPU. By the way, from Figure 1(c), the CSBCG method works as predicted in [18, 19], smoothing a little the convergence history of the BCG method. In particular, the CSBCG method computes a subset of the BCG iterates by clipping three “spikes” in the BCG history by means of the composite 2×2 steps from the particular investigations given in the bottom first two plots of Figures 1(d), 1(e) and Table 2.Table 2
Comparison results of Example 1 withTOL=10-6.
Method
BCG
BiCOR
CSBCG
CSBiCOR
Iters
330
318
327/3
318/0
CPU
40.7807
39.0256
40.1300
39.2116
TRR
−6.0318
−6.0150
−6.0318
−6.0150Convergence histories of Example 1 withTOL=10-6.
(a)
(b)
(c)
(d)
(e)
(f)Another interesting observation is that the good properties of the BiCOR method in terms of fast convergence speed and smooth convergence behavior in comparison with the BCG method are instinctually inherited by the CSBiCOR method so that the CSBiCOR method outperforms the CSBCG method in aspects of Iters and CPU as well as smooth effect, as shown in Table2 and Figure 1(f).


## 7.2. Example 2: Kim/Kim1

This is a good test example to vividly demonstrate the efficacy of the composite step strategy as depicted in Figure2. The BCG method, as observed in Figure 2(a), takes on many local “spikes” in the convergence curve. In such a case, the CSBCG method does a wonderful attempt to smooth the convergence history of the norm of the residuals as shown in Figure 2(c), although not leading to a monotonic decreasing convergence behavior. So does the CSBiCOR method to the BiCOR method by observation of Figure 2(d). Moreover, the CSBiCOR method seems to take less CPU than the BiCOR method while keeping approximately the same TRR when the BiCOR method is implemented using the same code as for the CSBiCOR method and the pivot test is modified to always choose 1-step update steps. Finally, the CSBiCOR method is again shown to be competitive to the CSBCG method as seen in Table 3 and Figure 2(b).Table 3
Comparison results of Example 2 withTOL=10-6.
Method
BCG
BiCOR
CSBCG
CSBiCOR
Iters
248
182
158/90
129/56
CPU
53.8496
39.2517
43.3087
33.8528
TRR
−6.3680
−6.0909
−6.3680
−6.0483Convergence histories of Example 2 withTOL=10-6.
(a)
(b)
(c)
(d)


## 7.3. Example 3: HB/Young1c

In the third example, the BCG and BiCOR methods have small irregularities and do not converge superlinearly, as reflected in Figure3(a). The two methods above seem to have almost the same “asymptotical” speed of convergence while the BiCOR method seems a little smoother than the BCG method. From Figures 3(c) and 3(d), the composite step strategy seems to play a surprisingly good stabilizing effect to the BCG and BiCOR method. Furthermore, the CSBCG and CSBiCOR methods take less CPU, respectively, than the BCG and BiCOR methods as reported in Table 4. It is stressed that although the CSBiCOR method consumes a little more Iters and CPU than the CSBCG method, it presents a little more smooth convergence behavior than the latter method as displayed in Figure 3(b). This is still thanks to the promising advantages of the empirically observed stability and fast convergence rate of the BiCOR method over the BCG method.Table 4
Comparison results of Example 3 withTOL=10-6.
Method
BCG
BiCOR
CSBCG
CSBiCOR
Iters
210
208
158/52
171/41
CPU
0.6848
0.5974
0.5372
0.5573
TRR
−6.1155
−6.0984
−6.0389
−6.0017Convergence histories of Example 3 withTOL=10-6.
(a)
(b)
(c)
(d)


## 8. Concluding Remarks

We have presented a new interesting variant of the BiCOR method for solving non-Hermitian systems of linear equations. Our approach is naturally based on and inspired by the composite step strategy taken for the CSBCG method [18, 19]. It is noted that exact pivot breakdowns are rare in practice; however, near breakdowns could cause severe numerical instability. The resulting CSBiCOR method is both theoretically and numerically demonstrated to avoid near pivot breakdowns and compute all the well-defined BiCOR iterates stably with only minor modifications with the assumption that the underlying biconjugate A-orthonormalization procedure does not break down. Besides reducing the number of spikes in the convergence history of the norm of the residuals to the greatest extent, the CSBiCOR method could provide some further practically desired smoothing behavior towards stabilizing the behavior of the BiCOR method when it has erratic convergence behaviors. Additionally, the CSBiCOR method seems to be superior to the CSBCG method to some extent because of the inherited promising advantages of the empirically observed stability and fast convergence rate of the BiCOR method over the BCG method.Since the BiCOR method is the most basic variant of the family of Lanczos biconjugateA-orthonormalization methods, its improvement will analogously lead to similar improvements for the CORS and BiCORSTAB methods, which is under investigation.It is important to note the nonnegligible added complexity when deciding the composite step strategy. That is we have to make some extra computation before deciding to take a2×2 step. Therefore, when the 2×2 step will not be taken after the decision test, it will lead to extra computation as well as CPU consuming time. It seems critical for us to optimize the code to implement the CSBiCOR method.



---

*Source: 408167-2013-02-20.xml*