Biconjugate gradient method

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

$Ax=b.$

Unlike the conjugate gradient method, this algorithm does not require the matrix $A$ to be self-adjoint, but instead one needs to perform multiplications by the conjugate transpose A^*.

The Algorithm

1. Choose initial guess $x_0$, two other vectors $x_0^{*}$ and $b^{*}$ and a preconditioner $M$
2. $r_0\leftarrow b-A{x}_0$
3. $r_0^{*}\leftarrow b^{*}-x_0^{*}{A}^{*}$
4. $p_0\leftarrow M^{-1}{r}_0$
5. $p_0^{*}\leftarrow r_0^{*}{M}^{-1}$
6. for $k=0,1,\dots$ do
   1. ${\alpha }_k\leftarrow \frac{r_k^{*}{M}^{-1}{r}_k}{p_k^{*}A{p}_k}$
   2. $x_{k+1}\leftarrow x_k+{\alpha }_k\cdot p_k$
   3. $x_{k+1}^{*}\leftarrow x_k^{*}+\overline{{\alpha }_k}\cdot p_k^{*}$
   4. $r_{k+1}\leftarrow r_k-{\alpha }_k\cdot A{p}_k$
   5. $r_{k+1}^{*}\leftarrow r_k^{*}-\overline{{\alpha }_k}\cdot p_k^{*}{A}^{*}$
   6. ${\beta }_k\leftarrow \frac{r_{k+1}^{*}{M}^{-1}{r}_{k+1}}{r_k^{*}{M}^{-1}{r}_k}$
   7. $p_{k+1}\leftarrow M^{-1}{r}_{k+1}+{\beta }_k\cdot p_k$
   8. $p_{k+1}^{*}\leftarrow r_{k+1}^{*}{M}^{-1}+\overline{{\beta }_k}\cdot p_k^{*}$

In the above formulation, the computed $r_k$ and $r_k^{*}$ satisfy

$r_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 systems

$Ax=b,$

$x^{*}{A}^{*}=b^{*};$

$x^{*}$ is the adjoint, and $\bar{\alpha }$ is the complex conjugate.

Unpreconditioned version of the algorithm

1. Choose initial guess $x_0$,
2. $r_0\leftarrow b-A{x}_0$
3. ${\hat{r}}_0\leftarrow \hat{b}-{\hat{x}}_{0}A$
4. $p_0\leftarrow r_0$
5. ${\hat{p}}_0\leftarrow {\hat{r}}_0$
6. for $k=0,1,\dots$ do
   1. ${\alpha }_k\leftarrow \frac{{\hat{r}}_k{r}_k}{{\hat{p}}_{k}A{p}_k}$
   2. $x_{k+1}\leftarrow x_k+{\alpha }_k\cdot p_k$
   3. ${\hat{x}}_{k+1}\leftarrow {\hat{x}}_k+{\alpha }_k\cdot {\hat{p}}_k$
   4. $r_{k+1}\leftarrow r_k-{\alpha }_k\cdot A{p}_k$
   5. ${\hat{r}}_{k+1}\leftarrow {\hat{r}}_k-{\alpha }_k\cdot {\hat{p}}_{k}A$
   6. ${\beta }_k\leftarrow \frac{{\hat{r}}_{k+1}{r}_{k+1}}{{\hat{r}}_k{r}_k}$
   7. $p_{k+1}\leftarrow r_{k+1}+{\beta }_k\cdot p_k$
   8. ${\hat{p}}_{k+1}\leftarrow {\hat{r}}_{k+1}+{\beta }_k\cdot {\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}(b-A{x}_j),$

$x_k^{*}:=x_j^{*}+(b^{*}-x_j^{*}A)P_k{A}^{-1},$

where $j<k$ using the related projection

$P_k:={\mathbf{𝐮}}_k{\left({\mathbf{𝐯}}_k^{*}A{\mathbf{𝐮}}_k\right)}^{-1}{\mathbf{𝐯}}_k^{*}A,$

with

${\mathbf{𝐮}}_k=[u_0,u_1,\dots ,u_{k-1}],$

${\mathbf{𝐯}}_k=[v_0,v_1,\dots ,v_{k-1}].$

These related projections may be iterated themselves as

$P_{k+1}=P_k+(1-P_k)u_k\otimes \frac{v_k^{*}A(1-P_k)}{v_k^{*}A(1-P_k)u_k}.$

A relation to Quasi-Newton methods is given by $P_k=A_k^{-1}A$ and $x_{k+1}=x_k-A_{k+1}^{-1}(A{x}_k-b)$, where

$A_{k+1}^{-1}=A_k^{-1}+(1-A_k^{-1}A)u_k\otimes \frac{v_k^{*}(1-A{A}_k^{-1})}{v_k^{*}A(1-A_k^{-1}A)u_k}.$

The new directions

$p_k=(1-P_k)u_k,$

$p_k^{*}=v_k^{*}A(1-P_k)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(1-P_k)A^{-1}{r}_j,$

$r_k^{*}=r_j^{*}(1-P_k)$

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 A⁻¹ 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_i^{*}A{p}_j=r_i^{*}{M}^{-1}{r}_j=0$ for $i\ne j$.
• if $P_{j^{\prime }}$ is a polynomial with $\mathrm{deg}\left(P_{j^{\prime }}\right)+j<k$, then $r_k^{*}{P}_{j^{\prime }}\left(M^{-1}A\right)u_j=0$. The algorithm thus produces projections onto the Krylov subspace.
• if $P_{i^{\prime }}$ is a polynomial with $i+\mathrm{deg}\left(P_{i^{\prime }}\right)<k$, then $v_i^{*}{P}_{i^{\prime }}\left(A{M}^{-1}\right)r_k=0$.

See also

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

References

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