Biconjugate gradient method explained

In mathematics, more specifically in numerical linear algebra, the biconjugate gradient method is an algorithm to solve systems of linear equations

A x = b. A to be self-adjoint, but instead one needs to perform multiplications by the conjugate transpose .

The Algorithm

1. Choose initial guess

x _{0} , two other vectors x ^{*} _{0} and b ^{*} and a preconditioner M r _{0} \leftarrow b - A x _{0} r ^{*} _{0} \leftarrow b ^{*} _{0} A ^{*} p _{0} \leftarrow M ^{-1} r _{0} p ^{*} _{0} \leftarrow r ^{*M} _{0} ^{-1}

1. for

k =0 , 1 , \ldots do\alpha _{k} \leftarrow {r ^{*} _{k} M ^{-1} r _{k} \over p ^{*} _{k} A p _{k}} x _{k+1} \leftarrow x _{k} + \alpha _{k} ⋅ p _{k} x ^{*} _{k +1} \leftarrow x ^{*} _{k} + \overline{\alpha _{k} ⋅} p ^{*} _{k} r _{k+1} \leftarrow r _{k} - \alpha _{k} ⋅ A p _{k} r ^{*} _{k +1} \leftarrow r ^{* -} _{k} \overline{\alpha _{k}} ⋅ p ^{*} _{k} A ^{*} \beta _{k} \leftarrow {r ^{*} _{k +1} M ^{-1} r _{k+1} \over r ^{*} _{k} M ^{-1} r _{k}} p _{k+1} \leftarrow M ^{-1} r _{k+1} + \beta _{k} ⋅ p _{k} p ^{*} _{k +1} \leftarrow r ^{*M} _{k +1} ^{-1} + \overline{\beta _{k} ⋅} p ^{*} _{k}

In the above formulation, the computed

r _{k} and r ^{*} _{k} satisfyr _{k} = b - A x _{k,} r ^{*} _{k} = b ^{*} - x ^{*} _{k} A ^{*}

and thus are the respective residuals corresponding to

x _{k} and x ^{*} _{k} , as approximate solutions to the systemsA x = b, x ^{*} A ^{*} = b ^{*;} x ^{*} is the adjoint, and \overline{\alpha} is the complex conjugate.

Unpreconditioned version of the algorithm

1. Choose initial guess

x _{0} ,r _{0} \leftarrow b - A x _{0} \hat{r} _{0} \leftarrow \hat{b} - \hat{x} ^{*} _{0 A} p _{0} \leftarrow r _{0} \hat{p} _{0} \leftarrow \hat{r} _{0}

1. for

k =0 , 1 , \ldots do\alpha _{k} \leftarrow {\hat{r} _{k} r _{k} \over \hat{p} _{k} A p _{k}} x _{k+1} \leftarrow x _{k} + \alpha _{k} ⋅ p _{k} \hat{x} _{k+1} \leftarrow \hat{x} _{k} + \alpha _{k} ⋅ \hat{p} _{k} r _{k+1} \leftarrow r _{k} - \alpha _{k} ⋅ A p _{k} \hat{r} _{k+1} \leftarrow \hat{r} _{k-} \alpha _{k} ⋅ \hat{p} _{k} A ^{*} \beta _{k} \leftarrow {\hat{r} _{k+1} r _{k+1} \over \hat{r} _{k} r _{k}} p _{k+1} \leftarrow r _{k+1} + \beta _{k} ⋅ p _{k} \hat{p} _{k+1} \leftarrow \hat{r} _{k+1} + \beta _{k} ⋅ \hat{p} _{k}

Discussion

The biconjugate gradient method is numerically unstable (compare to the biconjugate gradient stabilized method), but very important from a theoretical point of view. Define the iteration steps by

x _{k:=x} _{j+} P _{k} A ^{-1} \left(b - A x _{j} \right), x ^{*: =} _{k} x ^{* + \left(b} _{j} ^{*-} x ^{*} _{j} A \right) P _{k} A ^{-1} ,

where

j<k using the related projectionP _{k:=} u _{k} \left( v ^{*} _{k} A u _{k} \right) ^{-1} v ^{*} _{k} A,

with

u _{k=\left[u} _{0,} u _{1,} ... , u _{k-1} \right], v _{k=\left[v} _{0,} v _{1,} ... , v _{k-1} \right].

These related projections may be iterated themselves as

P _{k+1} = P _{k+} \left( 1- P _{k\right)} u _{k} ⊗ {v ^{*} _{k} A\left( 1- P _{k} \right) \over v ^{*} _{k} A\left( 1- P _{k} \right) u _{k}.}

A relation to Quasi-Newton methods is given by

P _{k=} A ^{-1} _{k} A and x _{k+1} = x _{k-} A ^{-1} _{k +1} \left(A x _{k} - b \right) , whereA ^{-1} _{k +1} = A ^{-1} _{k} + \left( 1- A ^{-1} _{k} A\right) u _{k} ⊗ {v ^{*} _{k} \left( 1- A A ^{-1} _{k} \right) \over v ^{*} _{k} A\left( 1- A ^{-1} _{k} A \right) u _{k}.}

The new directions

p _{k} = \left( 1- P _{k} \right) u _{k,} p ^{*} _{k} = v ^{*} _{k} A \left( 1- P _{k} \right) A ^{-1}

are then orthogonal to the residuals:

v ^{*} _{i} r _{k=} p ^{*} _{i} r _{k=0,} r ^{*} _{k} u _{j} = r ^{*} _{k} p _{j=} 0 ,

which themselves satisfy

r _{k=} A \left( 1- P _{k} \right) A ^{-1} r _{j,} r ^{* =} _{k} r ^{*} _{j} \left( 1- P _{k} \right)

where

i,j<k .

The biconjugate gradient method now makes a special choice and uses the setting

u _{k} = M ^{-1} r _{k,} v ^{*} _{k} = r ^{*} _{k} M ^{-1} .

With this particular choice, explicit evaluations of

P _{k} and are avoided, and the algorithm takes the form stated above.

Properties

• If

A = A ^{*} is self-adjoint, x ^{* =} _{0} x _{0} and b ^{*=b} , then r _{k=} r ^{*} _{k} , p _{k=} p ^{*} _{k} , and the conjugate gradient method produces the same sequence x _{k=} x ^{*} _{k} at half the computational cost.

• The sequences produced by the algorithm are biorthogonal, i.e.,

p ^{*Ap} _{j = r} ^{*M} ^{i} ^{-1} r _{j=0} for i ≠ j .

• if

P _{j'} is a polynomial with \deg\left(P _{j'} \right) + j<k , then r ^{*P} _{j'} \left(M ^{-1} A\right)u _{j=0} . The algorithm thus produces projections onto the Krylov subspace.

• if

P _{i'} is a polynomial with i + \deg\left(P _{i'} \right)<k , then v ^{*P} _{i'} \left(AM ^{-1} \right)r _{k=0} .

See also

• Biconjugate gradient stabilized method (BiCG-Stab)
• Conjugate gradient method (CG)
• Conjugate gradient squared method (CGS)

References

• R.. Fletcher. 1976. Conjugate gradient methods for indefinite systems. Numerical Analysis. 506. Lecture Notes in Mathematics. Springer Berlin / Heidelberg. 1617-9692. 978-3-540-07610-0. 73 - 89. 10.1007/BFb0080109. Watson. G. Alistair. free.
• Book: Press. WH. Teukolsky. SA. Vetterling. WT. Flannery. BP. 2007. Numerical Recipes: The Art of Scientific Computing. 3rd. Cambridge University Press. New York. 978-0-521-88068-8. Section 2.7.6. http://apps.nrbook.com/empanel/index.html?pg=87 .