Biconjugate gradient stabilized method

In the following sections, (x,y) = x^T y denotes the dot product of vectors. 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. ρ₀ = (r̂₀, r₀)
4. p₀ = r₀
5. For i = 1, 2, 3, …
   1. v = Ap_i−1
   2. α = ρ_i−1/(r̂₀, v)
   3. h = x_i−1 + αp_i−1
   4. s = r_i−1 − αv
   5. If h is accurate enough, i.e., if s is small enough, then set x_i = h and quit
   6. t = As
   7. ω = (t, s)/(t, t)
   8. x_i = h + ωs
   9. r_i = s − ωt
   10. If x_i is accurate enough, i.e., if r_i is small enough, then quit
   11. ρ_i = (r̂₀, r_i)
   12. β = (ρ_i/ρ_i−1)(α/ω)
   13. p_i = r_i + β(p_i−1 − ωv)

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. ρ₀ = (r̂₀, r₀)
4. p₀ = r₀
5. For i = 1, 2, 3, …
   1. y = K −1
      2 K −1
      1 p_i−1
   2. v = Ay
   3. α = ρ_i−1/(r̂₀, v)
   4. h = x_i−1 + αy
   5. s = r_i−1 − αv
   6. If h is accurate enough then x_i = h and quit
   7. z = K −1
      2 K −1
      1 s
   8. t = Az
   9. ω = (K −1
      1 t, K −1
      1 s)/(K −1
      1 t, K −1
      1 t)
   10. x_i = h + ωz
   11. r_i = s − ωt
   12. If x_i is accurate enough then quit
   13. ρ_i = (r̂₀, r_i)
   14. β = (ρ_i/ρ_i−1)(α/ω)
   15. p_i = r_i + β(p_i−1 − ωv)

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.

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).

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)r₀ 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.

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).