Biconjugate gradient stabilized method

In numerical linear algebra, the biconjugate gradient stabilized method, often abbreviated as BiCGSTAB, is an iterative method developed by H. A. van der Vorst for the numerical solution of nonsymmetric linear systems. It is a variant of the biconjugate gradient method (BiCG) and has faster and smoother convergence than the original BiCG as well as other variants such as the conjugate gradient squared method (CGS). It is a Krylov subspace method.

Algorithmic steps

Unpreconditioned BiCGSTAB

To solve a linear system Ax = b, BiCGSTAB starts with an initial guess x₀ and proceeds as follows:

1. r₀ = b − Ax₀
2. Choose an arbitrary vector r̂₀ such that (r̂₀, r₀) ≠ 0, e.g., r̂₀ = r₀
3. ρ₀ = α = ω₀ = 1
4. v₀ = p₀ = 0
5. For i = 1, 2, 3, …
   1. ρ_i = (r̂₀, r_i−1)
   2. β = (ρ_i/ρ_i−1)(α/ω_i−1)
   3. p_i = r_i−1 + β(p_i−1 − ω_i−1v_i−1)
   4. v_i = Ap_i
   5. α = ρ_i/(r̂₀, v_i)
   6. s = r_i−1 − αv_i
   7. t = As
   8. ω_i = (t, s)/(t, t)
   9. x_i = x_i−1 + αp_i + ω_is
   10. If x_i is accurate enough then quit
   11. r_i = s − ω_it

Preconditioned BiCGSTAB

Preconditioners are usually used to accelerate convergence of iterative methods. To solve a linear system Ax = b with a preconditioner K = K₁K₂ ≈ A, preconditioned BiCGSTAB starts with an initial guess x₀ and proceeds as follows:

1. r₀ = b − Ax₀
2. Choose an arbitrary vector r̂₀ such that (r̂₀, r₀) ≠ 0, e.g., r̂₀ = r₀
3. ρ₀ = α = ω₀ = 1
4. v₀ = p₀ = 0
5. For i = 1, 2, 3, …
   1. ρ_i = (r̂₀, r_i−1)
   2. β = (ρ_i/ρ_i−1)(α/ω_i−1)
   3. p_i = r_i−1 + β(p_i−1 − ω_i−1v_i−1)
   4. y = K⁻¹p_i
   5. v_i = Ay
   6. α = ρ_i/(r̂₀, v_i)
   7. s = r_i−1 − αv_i
   8. z = K⁻¹s
   9. t = Az
   10. ω_i = (K −1
       1 t, K −1
       1 s)/(K −1
       1 t, K −1
       1 t)
   11. x_i = x_i−1 + αy + ω_iz
   12. If x_i is accurate enough then quit
   13. r_i = s − ω_it

This formulation is equivalent to applying unpreconditioned BiCGSTAB to the explicitly preconditioned system

Ãx̃ = b̃

with à = K −1
1 AK −1
2 , x̃ = K₂x and b̃ = K −1
1 b. In other words, both left- and right-preconditioning are possible with this formulation.

Derivation

BiCG in polynomial form

In BiCG, the search directions p_i and p̂_i and the residuals r_i and r̂_i are updated using the following recurrence relations:

p_i = r_i−1 + β_ip_i−1,

p̂_i = r̂_i−1 + β_ip̂_i−1,

r_i = r_i−1 − α_iAp_i,

r̂_i = r̂_i−1 − α_iA^Tp̂_i.

The constants α_i and β_i are chosen to be

α_i = ρ_i/(p̂_i, Ap_i),

β_i = ρ_i/ρ_i−1

where ρ_i = (r̂_i−1, r_i−1) so that the residuals and the search directions satisfy biorthogonality and biconjugacy, respectively, i.e., for i ≠ j,

(r̂_i, r_j) = 0,

(p̂_i, Ap_j) = 0.

It is straightforward to show that

r_i = P_i(A)r₀,

r̂_i = P_i(A^T)r̂₀,

p_i+1 = T_i(A)r₀,

p̂_i+1 = T_i(A^T)r̂₀

where P_i(A) and T_i(A) are ith-degree polynomials in A. These polynomials satisfy the following recurrence relations:

P_i(A) = P_i−1(A) − α_iAT_i−1(A),

T_i(A) = P_i(A) − β_i+1T_i−1(A).

Derivation of BiCGSTAB from BiCG

It is unnecessary to explicitly keep track of the residuals and search directions of BiCG. In other words, the BiCG iterations can be performed implicitly. In BiCGSTAB, one wishes to have recurrence relations for

r̃_i = Q_i(A)P_i(A)r₀

where Q_i(A) = (I − ω₁A)(I − ω₂A)⋯(I − ω_iA) with suitable constants ω_j instead of r_i = P_i(A) in the hope that Q_i(A) will enable faster and smoother convergence in r̃_i than r_i.

It follows from the recurrence relations for P_i(A) and T_i(A) and the definition of Q_i(A) that

Q_i(A)P_i(A)r₀ = (I − ω_iA)(Q_i−1(A)P_i−1(A)r₀ − α_iAQ_i−1(A)T_i−1(A)r₀),

which entails the necessity of a recurrence relation for Q_i(A)T_i(A)r₀. This can also be derived from the BiCG relations:

Q_i(A)T_i(A)r₀ = Q_i(A)P_i(A)r₀ + β_i+1(I − ω_iA)Q_i−1(A)P_i−1(A)r₀.

Similarly to defining r̃_i, BiCGSTAB defines

p̃_i+1 = Q_i(A)T_i(A)r₀.

Written in vector form, the recurrence relations for p̃_i and r̃_i are

p̃_i = r̃_i−1 + β_i(I − ω_i−1A)p̃_i−1,

r̃_i = (I − ω_iA)(r̃_i−1 − α_iAp̃_i).

To derive a recurrence relation for x_i, define

s_i = r̃_i−1 − α_iAp̃_i.

The recurrence relation for r̃_i can then be written as

r̃_i = r̃_i−1 − α_iAp̃_i − ω_iAs_i,

which corresponds to

x_i = x_i−1 + α_ip̃_i + ω_is_i.

Determination of BiCGSTAB constants

Now it remains to determine the BiCG constants α_i and β_i and choose a suitable ω_i.

In BiCG, β_i = ρ_i/ρ_i−1 with

ρ_i = (r̂_i−1, r_i−1) = (P_i−1(A^T)r̂₀, P_i−1(A)r₀).

Since BiCGSTAB does not explicitly keep track of r̂_i or r_i, ρ_i is not immediately computable from this formula. However, it can be related to the scalar

ρ̃_i = (Q_i−1(A^T)r̂₀, P_i−1(A)r₀) = (r̂₀, Q_i−1(A)P_i−1(A)r₀) = (r̂₀, r_i−1).

Due to biorthogonality, r_i−1 = P_i−1(A)r₀ is orthogonal to U_i−2(A^T)r̂₀ where U_i−2(A^T) is any polynomial of degree i − 2 in A^T. Hence, only the highest-order terms of P_i−1(A^T) and Q_i−1(A^T) matter in the dot products (P_i−1(A^T)r̂₀, P_i−1(A)r₀) and (Q_i−1(A^T)r̂₀, P_i−1(A)r₀). The leading coefficients of P_i−1(A^T) and Q_i−1(A^T) are (−1)ⁱ⁻¹α₁α₂⋯α_i−1 and (−1)ⁱ⁻¹ω₁ω₂⋯ω_i−1, respectively. It follows that

ρ_i = (α₁/ω₁)(α₂/ω₂)⋯(α_i−1/ω_i−1)ρ̃_i,

and thus

β_i = ρ_i/ρ_i−1 = (ρ̃_i/ρ̃_i−1)(α_i−1/ω_i−1).

A simple formula for α_i can be similarly derived. In BiCG,

α_i = ρ_i/(p̂_i, Ap_i) = (P_i−1(A^T)r̂₀, P_i−1(A)r₀)/(T_i−1(A^T)r̂₀, AT_i−1(A)r₀).

Similarly to the case above, only the highest-order terms of P_i−1(A^T) and T_i−1(A^T) matter in the dot products thanks to biorthogonality and biconjugacy. It happens that P_i−1(A^T) and T_i−1(A^T) have the same leading coefficient. Thus, they can be replaced simultaneously with Q_i−1(A^T) in the formula, which leads to

α_i = (Q_i−1(A^T)r̂₀, P_i−1(A)r₀)/(Q_i−1(A^T)r̂₀, AT_i−1(A)r₀) = ρ̃_i/(r̂₀, AQ_i−1(A)T_i−1(A)r₀) = ρ̃_i/(r̂₀, Ap̃_i).

Finally, BiCGSTAB selects ω_i to minimize r̃_i = (I − ω_iA)s_i in 2-norm as a function of ω_i. This is achieved when

((I − ω_iA)s_i, As_i) = 0,

giving the optimal value

ω_i = (As_i, s_i)/(As_i, As_i).

Generalization

BiCGSTAB can be viewed as a combination of BiCG and GMRES where each BiCG step is followed by a GMRES(1) (i.e., GMRES restarted at each step) step to repair the irregular convergence behavior of CGS, as an improvement of which BiCGSTAB was developed. However, due to the use of degree-one minimum residual polynomials, such repair may not be effective if the matrix A has large complex eigenpairs. In such cases, BiCGSTAB is likely to stagnate as confirmed by numerical experiments.

One may expect that higher-degree minimum residual polynomials may better handle this situation. This gives rise to algorithms including BiCGSTAB2 and the more general BiCGSTAB(l). In BiCGSTAB(l), a GMRES(l) step follows every l BiCG steps. BiCGSTAB2 is equivalent to BiCGSTAB(l) with l = 2.

See also

• Biconjugate gradient method
• Conjugate gradient squared method
• Conjugate gradient method

References

• Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems". SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035. 
• Saad, Y. (2003). "§7.4.2 BICGSTAB". Iterative Methods for Sparse Linear Systems (2nd ed.). SIAM. pp. 231–234. doi:10.2277/0898715342. ISBN 978-0-89871-534-7. 
• ^ Gutknecht, M. H. (1993). "Variants of BICGSTAB for Matrices with Complex Spectrum". SIAM J. Sci. Comput. 14 (5): 1020–1033. doi:10.1137/0914062. 
• ^ Sleijpen, G. L. G.; Fokkema, D. R. (November 1993). "BiCGstab(l) for linear equations involving unsymmetric matrices with complex spectrum" (PDF). Electronic Transactions on Numerical Analysis (Kent, OH: Kent State University) 1: 11–32. ISSN 1068-9613. 

Numerical linear algebra Key concepts Floating point Numerical stability Problems Matrix multiplication (algorithms) Matrix decompositions Linear equations Sparse problems Hardware CPU cache TLB Cache-oblivious algorithm SIMD Multiprocessing Software BLAS Specialized libraries General purpose software