Exploiting the Composite Step Strategy to the Biconjugate 𝐴-Orthogonal Residual Method for Non-Hermitian Linear Systems

Abstract

The Biconjugate 𝐴-Orthogonal Residual (BiCOR) method carried out in finite precision arithmetic by means of the biconjugate 𝐴-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 𝐴-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.

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 form 𝐴𝑥=𝑏. 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 𝐴-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 𝐴-orthonormalization methods [7].

As observed from numerous numerical experiments carried out with the Lanczos biconjugate 𝐴-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 𝐴 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 𝐴-Orthogonal Residual (BiCOR) method, which is the basic underlying variant of the Lanczos biconjugate 𝐴-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 𝐴-orthogonal Residual Squared (CORS) method and the Biconjugate 𝐴-Orthogonal Residual Stabilized (BiCORSTAB) method.

By comparison of the biconjugate 𝐴-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 𝐴-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 biconjugate 𝐴-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 “𝑇” the transpose of a vector, or matrix. For a non-Hermitian matrix 𝐴=(𝑎_𝑖𝑗)_𝑁×𝑁∈ℂ^𝑁×𝑁, the Hermitian conjugate of 𝐴 is denoted as 𝐴^𝐻≡𝐴^𝑇=𝑎_{𝑗𝑖𝑁×𝑁}. (1) The standard Hermitian inner product of two complex vectors 𝑢,𝑣∈ℂ^𝑁 is defined as ⟨𝑢,𝑣⟩=𝑢^𝐻𝑣=^𝑁_𝑖=1𝑢_𝑖𝑣_𝑖. (2) The nested Krylov subspace of dimension 𝑡 generated by 𝐴 from 𝑣 is of the form 𝒦_𝑡(𝐴,𝑣)=span𝑣,𝐴𝑣,𝐴²𝑣,…,𝐴^𝑡−1𝑣. (3) In addition, 𝑒_𝑖 denotes the 𝑖th 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 biconjugate 𝐴-orthonormalization procedure [7] is first briefly recalled as in Algorithm 1, which can ideally build up a pair of biconjugate 𝐴-orthonormal bases for the dual Krylov subspaces 𝒦_𝑚(𝐴,𝑣₁) and 𝒦_𝑚(𝐴^𝐻,𝑤₁), where 𝑣₁ and 𝑤₁ are chosen initially to satisfy certain conditions.

(1) Choose 𝑣1, 𝜔1, such that 𝜔1,𝐴𝑣1=1 (2) Set 𝛽1=𝛿1≡0, 𝜔0=𝑣0≡𝟎∈ℂ𝑁 (3) for  𝑗=1,2,…  do (4)  𝛼𝑗=𝜔𝑗,𝐴(𝐴𝑣𝑗) (5)  𝑣𝑗+1=𝐴𝑣𝑗−𝛼𝑗𝑣𝑗−𝛽𝑗𝑣𝑗−1 (6)  𝜔𝑗+1=𝐴𝐻𝜔𝑗−𝛼𝑗𝜔𝑗−𝛿𝑗𝜔𝑗−1 (7)  𝛿𝑗+1=𝜔𝑗+1,𝐴𝑣𝑗+11/2 (8)  𝛽𝑗+1=𝜔𝑗+1,𝐴𝑣𝑗+1/𝛿𝑗+1 (9)  𝑣𝑗+1=𝑣𝑗+1/𝛿𝑗+1 (10)   𝜔𝑗+1=𝜔𝑗+1/𝛽𝑗+1 (11) end for

Algorithm 1 

Biconjugate 𝐴-orthonormalization procedure.

Observe that the above algorithm is possible to have Lanczos-type breakdown whenever 𝛿_𝑗+1 vanishes while 𝑤_𝑗+1 and 𝐴𝑣_𝑗+1 are not equal to 𝟎∈ℂ^𝑁 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 Algorithm 1 are reviewed [7] in the following proposition for the preparation of the theoretical basis of the composite step method.

Proposition 1. If Algorithm 1 proceeds 𝑚 steps, then the right and left Lanczos-type vectors 𝑣_𝑗,𝑗=1,2,…,𝑚 and 𝑤_𝑖,𝑖=1,2,…,𝑚 form a biconjugate 𝐴-orthonormal system in exact arithmetic, that is, 𝜔_𝑖,𝐴𝑣_𝑗=𝛿_𝑖,𝑗,1≤𝑖,𝑗≤𝑚. (4)
Furthermore, denote by 𝑉_𝑚=[𝑣₁,𝑣₂,…,𝑣_𝑚] and 𝑊_𝑚=[𝑤₁,𝑤₂,…,𝑤_𝑚] the 𝑁×𝑚 matrices and by 𝑇_𝑚 the extended tridiagonal matrix of the form 𝑇_𝑚=𝑇_𝑚𝛿_𝑚+1𝑒^𝑇_𝑚, (5) where 𝑇_𝑚=𝛼₁𝛽₂𝛿₂𝛼₂𝛽₃𝛿⋱⋱⋱_𝑚−1𝛼_𝑚−1𝛽_𝑚𝛿_𝑚𝛼_𝑚, (6) whose entries are the coefficients generated during the algorithm implementation, and in which 𝛼₁,…,𝛼_𝑚,𝛽₂,…,𝛽_𝑚 are complex while 𝛿₂,…,𝛿_𝑚 are positive. Then with the biconjugate 𝐴-orthonormalization procedure, the following four relations hold: 𝐴𝑉_𝑚=𝑉_𝑚𝑇_𝑚+𝛿_𝑚+1𝑣_𝑚+1𝑒^𝑇_𝑚,𝐴^𝐻𝑊_𝑚=𝑊_𝑚𝑇^𝐻_𝑚+𝛽_𝑚+1𝜔_𝑚+1𝑒^𝑇_𝑚,𝑊^𝐻_𝑚𝐴𝑉_𝑚=𝐼_𝑚,𝑊^𝐻_𝑚𝐴²𝑉_𝑚=𝑇_𝑚. (7)(8)(9)(10)

It is well known that the Lanczos biorthogonalization procedure can mathematically build up a pair of biorthogonal bases for the dual Krylov subspaces 𝒦_𝑚(𝐴,𝑣₁) and 𝒦_𝑚(𝐴^𝐻,𝑤₁), where 𝑣₁ and 𝑤₁ are initially selected such that ⟨𝜔₁,𝑣₁⟩=1 [6]. While the biconjugate 𝐴-orthonormalization procedure presented as in Algorithm 1 can ideally build up a pair of biconjugate 𝐴-orthonormal bases for the dual Krylov subspaces 𝒦_𝑚(𝐴,𝑣₁) and 𝒦_𝑚(𝐴^𝐻,𝑤₁), where 𝑣₁ and 𝑤₁ are chosen initially to satisfy ⟨𝜔₁,𝐴𝑣₁⟩=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 𝐴-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 matrix 𝑇_𝑚 formed in Proposition 1 cannot have two successive singular leading principal submatrices if the coefficient matrix 𝐴 is nonsingular.

Proposition 2. Let 𝑇_𝑘(1≤𝑘≤𝑚) be the upper left principal submatrices of the nonsingular tridiagonal matrix 𝑇_𝑚 appeared as in Proposition 1. Then 𝑇_𝑘−1 and 𝑇_𝑘 cannot be both singular.

It is noted that there may exist possible breakdown in the factorization without pivoting of the tridiagonal matrix 𝑇_𝑚 formed in Proposition 1. The following results illustrate how to correct such problem with the occasional use of 2×2 block pivots.

Proposition 3. Let 𝑇_𝑚 be the nonsingular tridiagonal matrix formed in Proposition 1. Then 𝑇_𝑚 can be factored as 𝑇_𝑚=𝐿_𝑚𝐷_𝑚𝑈_𝑚, (11) where 𝐿_𝑚 is unit lower block bidiagonal, 𝑈_𝑚 is unit upper block bidiagonal, and 𝐷_𝑚 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 submatrices 𝑇_𝑘(1≤𝑘≤𝑚) of 𝑇_𝑚 is singular, it can still be factored.

Corollary 4. Suppose one upper left principal submatrix 𝑇_𝑘(1≤𝑘≤𝑚) of the nonsingular tridiagonal matrix 𝑇_𝑚 formed in Proposition 1 is singular, but 𝑇_𝑘−1 must be nonsingular as stated in Proposition 2. Then 𝑇_𝑘=𝐿_𝑘𝐷_𝑘𝑈_𝑘, (12) where 𝐿_𝑘 is unit lower block bidiagonal, 𝑈_𝑘 is unit upper block bidiagonal, and 𝐷_𝑘 is block diagonal, with 1×1 and 2×2 diagonal blocks, and in particular, the last block of 𝐷_𝑘 is the 1×1 zero matrix.

From Corollary 4, the singularity of 𝑇_𝑘(1≤𝑘≤𝑚) can be recognized by only examining the last diagonal element (and next potential pivot) 𝑑_𝑘. If 𝑑_𝑘=0, then the 2×2 block 𝑑_𝑘𝛽_𝑘+1𝛿_𝑘+1𝛼_𝑘+1 (13) will be nonsingular and can be used as a 2×2 pivot, where 𝛽_𝑘+1,𝛿_𝑘+1, and 𝛼_𝑘+1 are the corresponding elements of 𝑇_𝑚 as shown in Proposition 1.

It can be concluded that appropriate combinations of 1×1 and 2×2 steps for the BiCOR method can skip over the breakdowns caused by the singular principal submatrices of the tridiagonal matrix 𝑇_𝑚 generated by the underlying biconjugate 𝐴-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 guess 𝑥₀ to the non-Hermitian linear system 𝐴𝑥=𝑏 associated with the initial residual 𝑟₀=𝑏−𝐴𝑥₀, define a Krylov subspace ℒ_𝑚≡𝐴^𝐻span(𝑊_𝑚)=𝐴^𝐻𝒦_𝑚(𝐴^𝐻,𝑤₁), where 𝑊_𝑚 is defined in Proposition 1, 𝑣₁=𝑟₀/||𝑟₀||₂ and 𝑤₁ is chosen arbitrarily such that ⟨𝑤₁,𝐴𝑣₁⟩≠0. But 𝑤₁ is often chosen to be equal to 𝐴𝑣₁/||𝐴𝑣₁||²₂ subjecting to ⟨𝑤₁,𝐴𝑣₁⟩=1. It is worth noting that this choice for 𝑤₁ 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  𝑚 steps, we can seek an 𝑚th approximate solution 𝑥_𝑚 from the affine subspace 𝑥₀+𝒦_𝑚(𝐴,𝑣₁) of dimension 𝑚, by imposing the Petrov-Galerkin condition 𝑏−𝐴𝑥_𝑚⟂ℒ_𝑚, (14) which can be mathematically written in matrix formulation as 𝐴^𝐻𝑊_𝑚^𝐻𝑏−𝐴𝑥_𝑚=0. (15)

Analogously, an 𝑚th dual approximation 𝑥^∗_𝑚 of the corresponding dual system 𝐴^𝐻𝑥^∗=𝑏^∗ is sought from the affine subspace 𝑥^∗₀+𝒦_𝑚(𝐴^𝐻,𝑤₁) of dimension 𝑚 by satisfying 𝑏^∗−𝐴^𝐻𝑥^∗_𝑚⟂𝐴𝒦_𝑚𝐴,𝑣₁, (16) which can be mathematically written in matrix formulation as 𝐴𝑉_𝑚^𝐻𝑏^∗−𝐴^𝐻𝑥_𝑚=0, (17) where 𝑥^∗₀ is an initial dual approximate solution and 𝑉_𝑚 is defined in Proposition 1 with 𝑣₁=𝑟₀/||𝑟₀||₂.

Consequently, the BiCOR iterates 𝑥_𝑗’s can be computed by the coming Algorithm 2, which is just the unpreconditioned BiCOR method with the preconditioner 𝑀 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].

(1) Compute 𝑟0=𝑏−𝐴𝑥0 for some initial guess 𝑥0. (2) Choose 𝑟∗0=𝑃(𝐴)𝑟0 such that 𝑟∗0,𝐴𝑟0≠0, where 𝑃(𝑡) is a polynomial in 𝑡.    (e.g., 𝑟∗0=𝐴𝑟0). (3) Set 𝑝0=𝑟0, 𝑝∗0=𝑟∗0, 𝑞0=𝐴𝑝0, 𝑞∗0=𝐴𝐻𝑝∗0, 𝑟0̂=𝐴𝑟0, 𝜌0=⟨𝑟∗0,𝑟0̂⟩. (4) for  𝑛=0,1,…  do (5)  𝜎𝑛=⟨𝑞∗𝑛,𝑞𝑛⟩ (6)  𝛼𝑛=𝜌𝑛/𝜎𝑛 (7)  𝑥𝑛+1=𝑥𝑛+𝛼𝑛𝑝𝑛 (8)  𝑟𝑛+1=𝑟𝑛−𝛼𝑛𝑞𝑛 (9)  𝑥∗𝑛+1=𝑥∗𝑛+𝛼𝑛𝑝∗𝑛 (10)   𝑟∗𝑛+1=𝑟∗𝑛−𝛼𝑛𝑞∗𝑛 (11)   𝑟̂𝑛+1=𝐴𝑟𝑛+1 (12)   𝜌𝑛+1=𝑟∗𝑛+1,𝑟̂𝑛+1 (13)   if  𝜌𝑛+1=0, method fails (14)   𝛽𝑛+1=𝜌𝑛+1/𝜌𝑛 (15)   𝑝𝑛+1=𝑟𝑛+1+𝛽𝑛+1𝑝𝑛 (16)   𝑝∗𝑛+1=𝑟∗𝑛+1+𝛽𝑛+1𝑝∗𝑛 (17)   𝑞𝑛+1=𝑟̂𝑛+1+𝛽𝑛+1𝑞𝑛 (18)   𝑞∗𝑛+1=𝐴𝐻𝑝∗𝑛 (19)  check convergence; continue if necessary (20) end for

Algorithm 2 

Algorithm BiCOR.

Suppose Algorithm 2 runs successfully to step 𝑛, that is, 𝜎_𝑖≠0, 𝜌_𝑖≠0, 𝑖=0,1,…,𝑛−1. The BiCOR iterates satisfy the following properties [7].

Proposition 5. Let 𝑅_𝑛+1=[𝑟₀,𝑟₁,…,𝑟_𝑛], 𝑅^∗_𝑛+1=[𝑟^∗₀,𝑟^∗₁,…,𝑟^∗_𝑛] and 𝑃_𝑛+1=[𝑝₀,𝑝₁,…,𝑝_𝑛], 𝑃^∗_𝑛+1=[𝑝^∗₀,𝑝^∗₁,…,𝑝^∗_𝑛]. One has the following.(1)Range(𝑅_𝑛+1)=Range(𝑃_𝑛+1)=𝒦_𝑛+1(𝐴,𝑟₀), Range(𝑅^∗_𝑛+1)=Range(𝑃^∗_𝑛+1)=𝒦_𝑛+1(𝐴^𝐻,𝑟^∗₀).(2)𝑅^∗𝐻_𝑛+1𝐴𝑅_𝑛+1is diagonal.(3)𝑃^∗𝐻_𝑛+1𝐴²𝑃_𝑛+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)𝜌_𝑛≡⟨𝑟^∗_𝑛,𝑟_𝑛̂⟩≡⟨𝑟^∗_𝑛,𝐴𝑟_𝑛⟩=0 but 𝑟^∗_𝑛 and 𝐴𝑟_𝑛 are not equal to 𝟎∈ℂ^𝑁 appearing in line 14;(2)𝜎_𝑛≡⟨𝑞^∗_𝑛,𝑞_𝑛⟩≡⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛⟩=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 𝐴-orthonormalization procedure depicted as in Algorithm 1 does not break down; that is the situation where 𝜎_𝑛=0 while 𝜌_𝑛≠0.

4. The Composite Step BiCOR Method

Suppose Algorithm 2 comes across a situation where 𝜎_𝑛=0 after successful algorithm implementation up to step 𝑛 with the assumption that 𝜌_𝑛≠0, which indicates that the updates of 𝑥_𝑛+1,𝑟_𝑛+1,𝑥^∗_𝑛+1,𝑟^∗_𝑛+1 are not well defined. Taking the composite step idea, we will avoid division by 𝜎_𝑛=0 via skipping this (𝑛+1)th update and exploiting a composite step update to directly obtain the quantities in step (𝑛+2) with scaled versions of 𝑟_𝑛+1 and 𝑟^∗_𝑛+1 as well as with the previous primary search direction vector 𝑝_𝑛 and shadow search direction vector 𝑝^∗_𝑛. 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 vectors 𝑧_𝑛+1∈𝒦_𝑛+2(𝐴,𝑟₀) and 𝑧^∗_𝑛+1∈𝒦_𝑛+2(𝐴^𝐻,𝑟^∗₀) as follows: 𝑧_𝑛+1=𝜎_𝑛𝑟_𝑛+1=𝜎_𝑛𝑟_𝑛−𝜌_𝑛𝐴𝑝_𝑛,𝑧^∗_𝑛+1=𝜎_𝑛𝑟^∗_𝑛+1=𝜎_𝑛𝑟^∗_𝑛−𝜌_𝑛𝐴^𝐻𝑝^∗_𝑛, (18)(19) which are then used to look for the iterates 𝑥_𝑛+2∈𝑥₀+𝒦_𝑛+2(𝐴,𝑟₀) and 𝑥^∗_𝑛+2∈𝑥^∗₀+𝒦_𝑛+2(𝐴^𝐻,𝑟^∗₀) in step (𝑛+2) as follows: 𝑥_𝑛+2=𝑥_𝑛+𝑝_𝑛,𝑧_𝑛+1𝑓_𝑛,𝑥^∗_𝑛+2=𝑥^∗_𝑛+𝑝^∗_𝑛,𝑧^∗_𝑛+1𝑓^∗_𝑛, (20) where 𝑓_𝑛,𝑓^∗_𝑛∈ℂ². Correspondingly, the (𝑛+2)th primary residual 𝑟_𝑛+2∈𝒦_𝑛+3(𝐴,𝑟₀) and shadow residual 𝑟^∗_𝑛+2∈𝒦_𝑛+3(𝐴^𝐻,𝑟^∗₀) are, respectively, computed as 𝑟_𝑛+2=𝑟_𝑛−𝐴𝑝_𝑛,𝑧_𝑛+1𝑓_𝑛,𝑟^∗_𝑛+2=𝑟^∗_𝑛−𝐴^𝐻𝑝^∗_𝑛,𝑧^∗_𝑛+1𝑓^∗_𝑛. (21)(22) The biconjugate 𝐴-orthogonality condition between the BiCOR primary residuals and shadow residuals shown as Property (2) in Proposition 5 requires 𝑝^∗_𝑛,𝑧^∗_𝑛+1,𝐴𝑟_𝑛+2𝑝=0,_𝑛,𝑧_𝑛+1,𝐴^𝐻𝑟^∗_𝑛+2=0, (23) combining with (21) and (22) gives rise to the following two 2×2 systems of linear equations for, respectively, solving 𝑓_𝑛 and 𝑓^∗_𝑛∶𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛𝐴^𝐻𝑝^∗_𝑛,𝐴𝑧_𝑛+1𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑝_𝑛𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑧_𝑛+1𝑓_𝑛⁽¹⁾𝑓_𝑛⁽²⁾=𝑝^∗_𝑛,𝐴𝑟_𝑛𝑧^∗_𝑛+1,𝐴𝑟_𝑛,𝐴𝑝_𝑛,𝐴^𝐻𝑝^∗_𝑛𝐴𝑝_𝑛,𝐴^𝐻𝑧^∗_𝑛+1𝐴𝑧_𝑛+1,𝐴^𝐻𝑝^∗_𝑛𝐴𝑧_𝑛+1,𝐴^𝐻𝑧^∗_𝑛+1𝑓_𝑛^∗(1)𝑓_𝑛^∗(2)=𝐴𝑝_𝑛,𝑟^∗_𝑛𝐴𝑧_𝑛+1,𝑟^∗_𝑛. (24)(25)

Similarly, the (𝑛+2)th primary search direction vector 𝑝_𝑛+2∈𝒦_𝑛+3(𝐴,𝑟₀) and shadow search direction vector 𝑝^∗_𝑛+2∈𝒦_𝑛+3(𝐴^𝐻,𝑟^∗₀) in a composite step are computed with the following form: 𝑝_𝑛+2=𝑟_𝑛+2+𝑝_𝑛,𝑧_𝑛+1𝑔_𝑛,𝑝^∗_𝑛+2=𝑟^∗_𝑛+2+𝑝^∗_𝑛,𝑧^∗_𝑛+1𝑔^∗_𝑛, (26)(27) where 𝑔_𝑛,𝑔^∗_𝑛∈ℂ².

The biconjugate 𝐴²-orthogonality condition between the BiCOR primary search direction vectors and shadow search direction vectors shown as Property (3) in Proposition 5 requires 𝑝^∗_𝑛,𝑧^∗_𝑛+1,𝐴²𝑝_𝑛+2𝑝=0,_𝑛,𝑧_𝑛+1,𝐴^𝐻2𝑝^∗_𝑛+2=0, (28) combining with (26) and (27) results in the following two 2×2 systems of linear equations for respectively solving 𝑔_𝑛 and 𝑔^∗_𝑛∶𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛𝐴^𝐻𝑝^∗_𝑛,𝐴𝑧_𝑛+1𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑝_𝑛𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑧_𝑛+1𝑔_𝑛⁽¹⁾𝑔_𝑛⁽²⁾𝐴=−^𝐻𝑝^∗_𝑛,𝐴𝑟_𝑛+2𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑟_𝑛+2,𝐴𝑝_𝑛,𝐴^𝐻𝑝^∗_𝑛𝐴𝑝_𝑛,𝐴^𝐻𝑧^∗_𝑛+1𝐴𝑧_𝑛+1,𝐴^𝐻𝑝^∗_𝑛𝐴𝑧_𝑛+1,𝐴^𝐻𝑧^∗_𝑛+1𝑔_𝑛^∗(1)𝑔_𝑛^∗(2)=−𝐴𝑝_𝑛,𝐴^𝐻𝑟^∗_𝑛+2𝐴𝑧_𝑛+1,𝐴^𝐻𝑟^∗_𝑛+2. (29)(30)

Therefore, it could be able to advance from step 𝑛 to step (𝑛+2) to provide 𝑥_𝑛+2, 𝑟_𝑛+2, 𝑥^∗_𝑛+2, 𝑟^∗_𝑛+2, 𝑝_𝑛+2, 𝑝^∗_𝑛+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 of 2×2 steps is sufficient for CSBiCOR method to compute exactly all those well-defined BiCOR iterates stably provided that the underlying biconjugate 𝐴-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 a 1×1 or a 2×2 step to obtain 𝑝_𝑛, 𝑝^∗_𝑛, the following hold: 𝐴^𝐻𝑟^∗_𝑛,𝐴𝑝_𝑛=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑟_𝑛=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛=𝑞^∗_𝑛,𝑞_𝑛≡𝜎_𝑛. (31)

Proof. If 𝑝_𝑛,𝑝^∗_𝑛 are obtained via a 1×1 step, that is, 𝑝_𝑛=𝑟_𝑛+𝛽_𝑛𝑝_𝑛−1, 𝑝^∗_𝑛=𝑟^∗_𝑛+𝛽_𝑛𝑝^∗_𝑛−1 as shown in lines 15 and 16 of Algorithm 2, then with the biconjugate 𝐴²-orthogonality condition between the BiCOR primary search direction vectors and shadow search direction vectors shown as Property (3) in Proposition 5, it follows that 𝐴^𝐻𝑟^∗_𝑛,𝐴𝑝_𝑛=𝐴^𝐻𝑝^∗_𝑛−𝛽_𝑛𝑝^∗_𝑛−1,𝐴𝑝_𝑛=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛−𝛽_𝑛𝐴^𝐻𝑝^∗_𝑛−1,𝐴𝑝_𝑛=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛=𝑞^∗_𝑛,𝑞_𝑛≡𝜎_𝑛,𝐴^𝐻𝑝^∗_𝑛,𝐴𝑟_𝑛=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛−𝛽_𝑛𝑝_𝑛−1=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛−𝛽_𝑛𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛−1=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛=𝑞^∗_𝑛,𝑞_𝑛≡𝜎_𝑛. (32)
If 𝑝_𝑛,𝑝^∗_𝑛 are obtained via a 2×2 step, that is, 𝑝_𝑛=𝑟_𝑛+𝑔⁽¹⁾_𝑛−2𝑝_𝑛−2+𝑔⁽²⁾_𝑛−2𝑧_𝑛−1, 𝑝^∗_𝑛=𝑟^∗_𝑛+𝑔^∗(1)_𝑛−2𝑝^∗_𝑛−2+𝑔^∗(2)_𝑛−2𝑧^∗_𝑛−1 as presented in (26) and (27), then according to the biconjugate 𝐴²-orthogonality conditions mentioned above for deriving 𝑝_𝑛,𝑝^∗_𝑛 with a 2×2 step, it follows that 𝐴^𝐻𝑟^∗_𝑛,𝐴𝑝_𝑛=𝐴^𝐻𝑝^∗_𝑛−𝑔^∗(1)_𝑛−2𝑝^∗_𝑛−2−𝑔^∗(2)_𝑛−2𝑧^∗_𝑛−1,𝐴𝑝_𝑛=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛−𝑔^∗(1)_𝑛−2𝐴^𝐻𝑝^∗_𝑛−2,𝐴𝑝_𝑛−𝑔^∗(2)_𝑛−2𝐴^𝐻𝑧^∗_𝑛−1,𝐴𝑝_𝑛=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛=𝑞^∗_𝑛,𝑞_𝑛≡𝜎_𝑛,𝐴^𝐻𝑝^∗_𝑛,𝐴𝑟_𝑛=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛−𝑔⁽¹⁾_𝑛−2𝑝_𝑛−2−𝑔⁽²⁾_𝑛−2𝑧_𝑛−1=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛−𝑔⁽¹⁾_𝑛−2𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛−2−𝑔⁽²⁾_𝑛−2𝐴^𝐻𝑝^∗_𝑛,𝐴𝑧_𝑛−1=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛=𝑞^∗_𝑛,𝑞_𝑛≡𝜎_𝑛. (33)

With the above lemma, we can show the relationships among the four 2×2 coefficient matrices of the linear systems represented as in (24), (25), (29), and (30).

Lemma 7. The four 2×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 ⟨𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑝_𝑛⟩=⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑧_𝑛+1⟩. From the presentations of 𝑧_𝑛+1 and 𝑧^∗_𝑛+1 as respectively defined in (18) and (19), and following from Lemma 6, we have 𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑝_𝑛=𝐴^𝐻𝜎_𝑛𝑟^∗_𝑛−𝜌_𝑛𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛=𝜎_𝑛𝐴^𝐻𝑟^∗_𝑛,𝐴𝑝_𝑛−𝜌_𝑛𝐴^𝐻2𝑝^∗_𝑛,𝐴𝑝_𝑛=𝜎²_𝑛−𝜌_𝑛𝐴^𝐻2𝑝^∗_𝑛,𝐴𝑝_𝑛,𝐴^𝐻𝑝^∗_𝑛,𝐴𝑧_𝑛+1=𝐴^𝐻𝑝^∗_𝑛,𝐴𝜎_𝑛𝑟_𝑛−𝜌_𝑛𝐴𝑝_𝑛=𝜎_𝑛𝐴^𝐻𝑝^∗_𝑛,𝐴𝑟_𝑛−𝜌_𝑛𝐴^𝐻𝑝^∗_𝑛,𝐴²𝑝_𝑛=𝜎²_𝑛−𝜌_𝑛𝐴^𝐻2𝑝^∗_𝑛,𝐴𝑝_𝑛. (34)
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 a 2×2 step is always sufficient to skip over the pivot breakdown of the BiCOR method when 𝜎_𝑛=0, just as stated at the end of Section 2, can be shown in the next lemma.

Lemma 8. Suppose that the biconjugate 𝐴-orthonormalization procedure depicted as in Algorithm 1 underlying the BiCOR method does not breakdown; that is the situation where 𝜌_𝑖≠0, 𝑖=0,1,…,𝑛 and ⟨𝑧^∗_𝑛+1,𝐴𝑧_𝑛+1⟩≠0. If 𝜎_𝑛≡⟨𝑞^∗_𝑛,𝑞_𝑛⟩≡⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛⟩=0, then the 2×2 coefficient matrices in (24), (25), (29), and (30) are nonsingular.

Proof. It suffices to show that if 𝜎_𝑛≡⟨𝑞^∗_𝑛,𝑞_𝑛⟩≡⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑝_𝑛⟩=0, then ⟨𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑝_𝑛⟩≠0. From Lemma 7 and (18), we have 𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑝_𝑛=𝐴^𝐻𝑝^∗_𝑛,𝐴𝑧_𝑛+1,𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑝_𝑛=𝐴^𝐻𝑧^∗_𝑛+1𝑧,−^𝑛+1𝜌_𝑛1=−𝜌_𝑛𝑧^∗_𝑛+1,𝐴𝑧_𝑛+1≠0. (35)

The next proposition shows some properties of the CSBiCOR iterates similar to those of the BiCOR method presented in Proposition 5 with some minor changes.

Proposition 9. For the CSBiCOR algorithm, let 𝑅_𝑛+1=[𝑟₀,𝑟₁,…,𝑟_𝑛], 𝑅^∗_𝑛+1=[𝑟^∗₀,𝑟^∗₁,…,𝑟^∗_𝑛], where 𝑟_𝑖, 𝑟^∗_𝑖 are replaced by 𝑧_𝑖, 𝑧^∗_𝑖 appropriately if the 𝑖th step is a composite 2×2 step; and let 𝑃_𝑛+1=[𝑝₀,𝑝₁,…,𝑝_𝑛], 𝑃^∗_𝑛+1=[𝑝^∗₀,𝑝^∗₁,…,𝑝^∗_𝑛]. Assuming the underlying biconjugate 𝐴-orthonormalization procedure does not break down, one has the following properties.(1)𝑅𝑎𝑛𝑔𝑒(𝑅_𝑛+1)=𝑅𝑎𝑛𝑔𝑒(𝑃_𝑛+1)=𝒦_𝑛+1(𝐴,𝑟₀), 𝑅𝑎𝑛𝑔𝑒(𝑅^∗_𝑛+1)=𝑅𝑎𝑛𝑔𝑒(𝑃^∗_𝑛+1)=𝒦_𝑛+1(𝐴^𝐻,𝑟^∗₀).(2)𝑅^∗𝐻_𝑛+1𝐴𝑅_𝑛+1 is diagonal.(3)𝑃^∗𝐻_𝑛+1𝐴²𝑃_𝑛+1=𝑑𝑖𝑎𝑔(𝐷_𝑘), where 𝐷_𝑘 is either of order 1×1 or 2×2 and the sum of the dimensions of 𝐷_𝑘’s is (𝑛+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 of 2×2 steps is sufficient for the CSBiCOR method to compute exactly all those well-defined BiCOR iterates stably provided that the underlying biconjugate 𝐴-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 biconjugate 𝐴-orthonormalization procedure does not break down.

Proof. By construction, 𝑥_𝑛+2=𝑥_𝑛+𝑓_𝑛⁽¹⁾𝑝_𝑛+𝑓_𝑛⁽²⁾𝑧_𝑛+1; by induction, 𝑥_𝑛∈𝑥₀+𝒦_𝑛(𝐴,𝑟₀) and 𝑝_𝑛∈𝒦_𝑛+1(𝐴,𝑟₀); by definition, 𝑧_𝑛+1∈𝒦_𝑛+2(𝐴,𝑟₀). Then it follows that 𝑥_𝑛+2∈𝑥₀+𝒦_𝑛+2(𝐴,𝑟₀). From Proposition 9, we have 𝑟_𝑛+2⟂𝐴^𝐻𝒦_𝑛+2(𝐴^𝐻,𝑟^∗₀). Therefore, 𝑥_𝑛+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 𝒱_𝑚=span(𝑣₁,𝑣₂,…,𝑣_𝑚) and 𝒲_𝑚=span(𝑤₁,𝑤₂,…,𝑤_𝑚) denote the Krylov subspaces generated by the biconjugate 𝐴-orthonormalization procedure after running Algorithm 1  𝑚 steps. The norms associated with the spaces 𝒱_𝑁≡ℂ^𝑁 and 𝒲_𝑁≡ℂ^𝑁 are defined as |‖𝑣‖|²_𝑟=𝑣^𝐻𝑀_𝑟𝑣,|‖𝑤‖|²_𝑙=𝑤^𝐻𝑀_𝑙𝑤, (36) where, 𝑣,𝑤∈ℂ^𝑁 while 𝑀_𝑟 and 𝑀_𝑙 are symmetric positive definite matrices. Then we have the following proposition with the initial guess 𝑥₀ of the linear system being 𝟎∈ℂ^𝑁. The case with a nonzero initial guess can be treated adaptively.

Proposition 11. Suppose that for all 𝑣∈𝒱_𝑁 and for all 𝑤∈𝒲_𝑁, one has 𝑤^𝐻𝐴²𝑣≤Γ|‖𝑣‖|_𝑟|‖𝑤‖|_𝑙, (37) where, Γ is a constant independent of 𝑣 and 𝑤. Furthermore, assuming that for those steps in the composite step Biconjugate 𝐴-Orthogonal Residual (CSBiCOR) method in which we compute an approximation 𝑥_𝑚, we have inf_{𝑣∈𝒱𝑚𝑣|‖‖|𝑟=1}sup_{𝑤∈𝒲𝑚𝑤|‖‖|𝑙≤1}𝑤^𝐻𝐴²𝑣≥𝛾_𝑘≥𝛾>0. (38) Then, 𝑥−𝑥_𝑚𝑟≤Γ1+𝛾inf_{𝑣∈𝒱𝑚}|‖𝑥−𝑣‖|_𝑟. (39)

Proof. From the Petrov-Galerkin condition (15), we have for 𝑤∈𝒲_𝑚, 𝐴^𝐻𝑤^𝐻𝐴𝑥−𝐴𝑥_𝑚=𝑤^𝐻𝐴²𝑥−𝑥_𝑚=0, (40) yielding with arbitrary 𝑣∈𝒱_𝑚𝑤^𝐻𝐴²𝑥_𝑚−𝑣=𝑤^𝐻𝐴²(𝑥−𝑣). (41) It is noted that 𝑥_𝑚−𝑣∈𝒱_𝑚 because of the assumption of the zero initial guess. Then taking the sup of both sides of the above equation for all |‖𝑤‖|_𝑙≤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) which yields 𝑥−𝑥_𝑚𝑟≤|‖𝑥−𝑣‖|_𝑟+𝑥_𝑚−𝑣_𝑟≤Γ1+𝛾|‖𝑥−𝑣‖|_𝑟 (43) with the triangle inequality. The final assertion (39) follows immediately since 𝑣∈𝒱_𝑚 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 four 2×2 linear systems represented as in (24), (25), (29), and (30).

Lemma 12. For the four 2×2 linear systems represented as in (24), (25), (29), and (30), the following relations hold:(1)⟨𝑝^∗_𝑛,𝐴𝑟_𝑛⟩=⟨𝐴𝑝_𝑛,𝑟^∗_𝑛⟩=⟨𝑟^∗_𝑛,𝐴𝑟_𝑛⟩=⟨𝑟^∗_𝑛,𝑟_𝑛̂⟩≡𝜌_𝑛 and ⟨𝑧^∗_𝑛+1,𝐴𝑟_𝑛⟩=⟨𝐴𝑧_𝑛+1,𝑟^∗_𝑛⟩=0,𝑓_𝑛=𝑓^∗_𝑛.(2)⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑟_𝑛+2⟩=⟨𝐴𝑝_𝑛,𝐴^𝐻𝑟^∗_𝑛+2⟩=0 and ⟨𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑟_𝑛+2⟩=⟨𝐴𝑧_𝑛+1,𝐴^𝐻𝑟^∗_𝑛+2⟩=−𝜌_𝑛+2/𝑓_𝑛⁽²⁾, where 𝜌_𝑛+2≡⟨𝑟^∗_𝑛+2,𝐴𝑟_𝑛+2⟩,𝑔_𝑛=𝑔^∗_𝑛.(3)⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑧_𝑛+1⟩=⟨𝐴𝑝_𝑛,𝐴^𝐻𝑧^∗_𝑛+1⟩=−𝜃_𝑛+1/𝜌_𝑛, where 𝜃_𝑛+1≡⟨𝑧^∗_𝑛+1,𝐴𝑧_𝑛+1⟩.

Proof. (1) If a 1×1 step is taken, then 𝑝^∗_𝑛=𝑟^∗_𝑛+𝛽_𝑛𝑝^∗_𝑛−1, giving that ⟨𝑝^∗_𝑛,𝐴𝑟_𝑛⟩=⟨𝑟^∗_𝑛+𝛽_𝑛𝑝^∗_𝑛−1,𝐴𝑟_𝑛⟩=⟨𝑟^∗_𝑛,𝐴𝑟_𝑛⟩+𝛽_𝑛⟨𝑝^∗_𝑛−1,𝐴𝑟_𝑛⟩=⟨𝑟^∗_𝑛,𝐴𝑟_𝑛⟩ because of the last-term vanishment due to Proposition 9. Analogously, the iteration 𝑝_𝑛=𝑟_𝑛+𝛽_𝑛𝑝_𝑛−1 similarly gives ⟨𝐴𝑝_𝑛,𝑟^∗_𝑛⟩=⟨𝐴𝑟_𝑛+𝛽_𝑛𝐴𝑝_𝑛−1,𝑟^∗_𝑛⟩=⟨𝐴𝑟_𝑛,𝑟^∗_𝑛⟩+𝛽_𝑛⟨𝐴𝑝_𝑛−1,𝑟^∗_𝑛⟩=⟨𝐴𝑟_𝑛,𝑟^∗_𝑛⟩=⟨𝑟^∗_𝑛,𝐴𝑟_𝑛⟩ 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 ⟨𝑝^∗_𝑛,𝐴𝑟_𝑛⟩ = ⟨𝑟^∗_𝑛+𝑔^∗(1)_𝑛−2𝑝^∗_𝑛−2+𝑔^∗(2)_𝑛−2𝑧^∗_𝑛−1,𝐴𝑟_𝑛⟩ = ⟨𝑟^∗_𝑛,𝐴𝑟_𝑛⟩+𝑔^∗(1)_𝑛−2⟨𝑝^∗_𝑛−2,𝐴𝑟_𝑛⟩+𝑔^∗(2)_𝑛−2⟨𝑧^∗_𝑛−1,𝐴𝑟_𝑛⟩ = ⟨𝑟^∗_𝑛,𝐴𝑟_𝑛⟩. Similarly, with Proposition 9 and the Hermitian property of the standard Hermitian inner product, it can be shown that ⟨𝐴𝑝_𝑛,𝑟^∗_𝑛⟩= ⟨𝐴(𝑟_𝑛+𝑔⁽¹⁾_𝑛−2𝑝_𝑛−2+𝑔⁽²⁾_𝑛−2𝑧_𝑛−1),𝑟^∗_𝑛⟩ = ⟨𝐴𝑟_𝑛,𝑟^∗_𝑛⟩+𝑔⁽¹⁾_𝑛−2⟨𝐴𝑝_𝑛−2,𝑟^∗_𝑛⟩+𝑔⁽²⁾_𝑛−2⟨𝐴𝑧_𝑛−1,𝑟^∗_𝑛⟩ = ⟨𝐴𝑟_𝑛,𝑟^∗_𝑛⟩ = ⟨𝑟^∗_𝑛,𝐴𝑟_𝑛⟩. The equalities ⟨𝑧^∗_𝑛+1,𝐴𝑟_𝑛⟩=⟨𝐴𝑧_𝑛+1,𝑟^∗_𝑛⟩=0 directly follow from Proposition 9. Finally, because the coefficient matrices of the symmetric systems in (24) and (25) respectively governing 𝑓_𝑛 and 𝑓^∗_𝑛 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 𝑓_𝑛=𝑓^∗_𝑛.
(2) From Proposition 9, it can directly verify that ⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑟_𝑛+2⟩=⟨𝐴𝑝_𝑛,𝐴^𝐻𝑟^∗_𝑛+2⟩=0 with the fact that 𝐴𝑝_𝑛∈𝒦_𝑛+2(𝐴,𝑟₀) and 𝐴^𝐻𝑝^∗_𝑛∈𝒦_𝑛+2(𝐴^𝐻,𝑟^∗₀). Next, it can be noted that if 𝜎_𝑛=0, then 𝑓_𝑛⁽²⁾≠0. Otherwise, the first equation in the linear system (24) cannot hold for ⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑧_𝑛+1⟩≠0 as proved in Lemma 8 and ⟨𝑝^∗_𝑛,𝐴𝑟_𝑛⟩=𝜌_𝑛≠0 as shown in the above (1). Then by construction of (22) and Proposition 9, we have ⟨𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑟_𝑛+2⟩ = ⟨(1/𝑓_𝑛^∗(2))(𝑟^∗_𝑛−𝑟^∗_𝑛+2−𝑓_𝑛^∗(1)𝐴^𝐻𝑝^∗_𝑛),𝐴𝑟_𝑛+2⟩ = (1/𝑓_𝑛^∗(2))(⟨𝑟^∗_𝑛,𝐴𝑟_𝑛+2⟩−⟨𝑟^∗_𝑛+2,𝐴𝑟_𝑛+2⟩−𝑓_𝑛^∗(1)⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑟_𝑛+2⟩) = (−1/𝑓_𝑛^∗(2))⟨𝑟^∗_𝑛+2,𝐴𝑟_𝑛+2⟩ = (−1/𝑓_𝑛^∗(2))𝜌_𝑛+2 = (−1/𝑓_𝑛⁽²⁾)𝜌_𝑛+2 with the fact that 𝑓_𝑛=𝑓^∗_𝑛 proved in the above (1). By construction of (21) and Proposition 9, we can also show that ⟨𝐴𝑧_𝑛+1,𝐴^𝐻𝑟^∗_𝑛+2⟩=−𝜌_𝑛+2/𝑓_𝑛⁽²⁾. Similarly to the relationships between 𝑓_𝑛 and 𝑓^∗_𝑛, we can deduce that 𝑔_𝑛=𝑔^∗_𝑛, because the coefficient matrices of the symmetric systems in (29) and (30) respectively governing 𝑔_𝑛 and 𝑔^∗_𝑛 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 ⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑧_𝑛+1⟩ = ⟨(1/𝜌_𝑛)(𝜎_𝑛𝑟^∗_𝑛−𝑧^∗_𝑛+1),𝐴𝑧_𝑛+1⟩ = (1/𝜌_𝑛)(𝜎_𝑛⟨𝑟^∗_𝑛,𝐴𝑧_𝑛+1⟩−⟨𝑧^∗_𝑛+1,𝐴𝑧_𝑛+1⟩) = (−1/𝜌_𝑛)⟨𝑧^∗_𝑛+1,𝐴𝑧_𝑛+1⟩ = −𝜃_𝑛+1/𝜌_𝑛, where 𝜃_𝑛+1≡⟨𝑧^∗_𝑛+1,𝐴𝑧_𝑛+1⟩. Because it is already known that ⟨𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑝_𝑛⟩ = ⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑧_𝑛+1⟩ as proved in Lemma 7, we have ⟨𝐴𝑝_𝑛,𝐴^𝐻𝑧^∗_𝑛+1⟩ = ⟨𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑝_𝑛⟩ = ⟨𝐴^𝐻𝑝^∗_𝑛,𝐴𝑧_𝑛+1⟩ = −𝜃_𝑛+1/𝜌_𝑛. The proof is completed.

As a result of Lemmas 6 and 12, when a 2×2 step is taken, the solutions for 𝑓_𝑛 and 𝑔_𝑛 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: 𝜎_𝑛−𝜃_𝑛+1𝜌_𝑛−𝜃_𝑛+1𝜌_𝑛𝜁_𝑛+1𝑓_𝑛⁽¹⁾𝑓_𝑛⁽²⁾=𝜌_𝑛0,𝜎_𝑛−𝜃_𝑛+1𝜌_𝑛−𝜃_𝑛+1𝜌_𝑛𝜁_𝑛+1𝑔_𝑛⁽¹⁾𝑔_𝑛⁽²⁾=0𝜌_𝑛+2𝑓_𝑛⁽²⁾, (44) where 𝜁_𝑛+1=⟨𝐴^𝐻𝑧^∗_𝑛+1,𝐴𝑧_𝑛+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, 𝑓_𝑛 and 𝑔_𝑛 can be explicitly represented as 𝑓_𝑛=𝜁_𝑛+1𝜌_𝑛,𝜃_𝑛+1𝜌²_𝑛𝜂_𝑛,𝑔_𝑛=𝜌_𝑛+2𝜌_𝑛,𝜎_𝑛𝜌_𝑛+2𝜃_𝑛+1, (45) where 𝜂_𝑛≡𝜎_𝑛𝜁_𝑛+1𝜌²_𝑛−𝜃²_𝑛+1.

Combining all the recurrences discussed above for either a 1×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 𝐵 can be represented by Algorithm 3, where 𝑠_𝑛+1 and 𝑧_𝑛+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.

(1) Compute 𝑟0=𝑏−𝐴𝑥0 for some initial guess 𝑥0. (2) Choose 𝑟∗0=𝑃(𝐴)𝑟0 such that 𝑟∗0,𝐴𝑟0≠0, where 𝑃(𝑡) is a polynomial in 𝑡.    (e.g., 𝑟∗0=𝐴𝑟0). Set 𝑝0=𝑟0, 𝑝0=𝑟0̃, 𝑞0=𝐴𝑝0, 𝑞0=𝐴𝐻𝑝0. (3) Compute 𝜌0=𝑟0̃,𝐴𝑟0. (4) Begin LOOP (𝑛=0,1,2,…) (5) 𝜎𝑛=𝑞𝑛,𝑞𝑛 (6) 𝑠𝑛+1=𝜎𝑛𝑟𝑛−𝜌𝑛𝑞𝑛 (7) 𝑠̃𝑛+1=𝜎𝑛𝑟𝑛̃−𝜌𝑛𝑞𝑛 (8) 𝑦𝑛+1=𝐴𝑠𝑛+1 (9) 𝑦𝑛+1=𝐴𝐻𝑠̃𝑛+1 (10) 𝜃𝑛+1=𝑠̃𝑛+1,𝑦𝑛+1 (11) 𝜁𝑛+1=𝑦𝑛+1,𝑦𝑛+1 (12) if 1 × 1 step then (13)  𝛼𝑛=𝜌𝑛/𝜎𝑛 (14)  𝜌𝑛+1=𝜃𝑛+1/𝜎2𝑛 (15)  𝛽𝑛+1=𝜌𝑛+1/𝜌𝑛 (16)  𝑥𝑛+1=𝑥𝑛+𝛼𝑛𝑝𝑛 (17)  𝑟𝑛+1=𝑟𝑛−𝛼𝑛𝑞𝑛 (18)  𝑟̃𝑛+1=𝑟𝑛̃−𝛼𝑛𝑞𝑛 (19)  𝑝𝑛+1=𝑠𝑛+1/𝜎𝑛+𝛽𝑛+1𝑝𝑛 (20)  𝑞𝑛+1=𝑦𝑛+1/𝜎𝑛+𝛽𝑛+1𝑞𝑛 (21)  𝑞𝑛+1=𝑦𝑛+1/𝜎𝑛+𝛽𝑛+1𝑞𝑛 (22)  𝑛←𝑛+1 (23) else (24)  𝛿𝑛=𝜎𝑛𝜁𝑛+1𝜌2𝑛−𝜃2𝑛+1 (25)  𝛼𝑛=𝜁𝑛+1𝜌3𝑛/𝛿𝑛 (26)  𝛼𝑛+1=𝜃𝑛+1𝜌2𝑛/𝛿𝑛 (27)  𝑥𝑛+2=𝑥𝑛+𝛼𝑛𝑝𝑛+𝛼𝑛+1𝑠𝑛+1 (28)  𝑟𝑛+2=𝑟𝑛−𝛼𝑛𝑞𝑛−𝛼𝑛+1𝑦𝑛+1 (29)  𝑟̃𝑛+2=𝑟𝑛̃−𝛼𝑛𝑞𝑛−𝛼𝑛+1𝑦𝑛+1 (30)  solve  𝐵𝑧𝑛+2=𝑟𝑛+2 (31)  solve  𝐵𝐻𝑧𝑛+2=𝑟̃𝑛+2 (32)  𝑧𝑛+2=𝐴𝑧𝑛+2 (33)  𝑧𝑛+2=𝐴𝐻𝑧𝑛+2 (34)  𝜌𝑛+2=𝑧𝑛+2,𝑟𝑛+2 (35)  𝛽𝑛+1=𝜌𝑛+2/𝜌𝑛 (36)  𝛽𝑛+2=𝜌𝑛+2𝜎𝑛/𝜃𝑛+1 (37)