Biconjugate gradient method

In mathematics, more specifically in numerical analysis, 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 $x_{0}$, $y_{0}$, a regular preconditioner $M$ (frequently $M^{-1}=1$ is used) and $c$;
2. $r_{0}\leftarrow b-Ax_{0},s_{0}\leftarrow c-A^{*}y_{0}\,$;
3. $d_{0}\leftarrow M^{-1}r_{0},f_{0}\leftarrow M^{-*}s_{0}\,$;
4. for $k=0,1,\dots \,$ do
5. $\alpha _{k}\leftarrow {s_{k}^{*}M^{-1}r_{k} \over f_{k}^{*}Ad_{k}}\,$;
6. $x_{k+1}\leftarrow x_{k}+\alpha _{k}d_{k}\,$ $\left(y_{k+1}\leftarrow y_{k}+{\overline {\alpha _{k}}}f_{k}\right)\,$;
7. $r_{k+1}\leftarrow r_{k}-\alpha _{k}Ad_{k}\,$, $s_{k+1}\leftarrow s_{k}-{\overline {\alpha _{k}}}A^{*}f_{k}\,$ ($r_{k}=b-Ax_{k}$ and $s_{k}=c-A^{*}y_{k}$ are the residuums);
8. $\beta _{k}\leftarrow {s_{k+1}^{*}M^{-1}r_{k+1} \over s_{k}^{*}M^{-1}r_{k}}\,$;
9. $d_{k+1}\leftarrow M^{-1}r_{k+1}+\beta _{k}d_{k}\,$, $f_{k+1}\leftarrow M^{-*}s_{k+1}+{\overline {\beta _{k}}}f_{k}\,$.

Discussion

The bicg method is numerically unstable, but very important from theoretical point of view: Define the iteration steps by $x_{k}:=x_{j}+P_{k}A^{-1}\left(b-Ax_{j}\right)$ and $y_{k}:=y_{j}+A^{-*}P_{k}^{*}\left(c-A^{*}y_{j}\right)$ ($j<k$) using the related projection

$P_{k}:=\mathbf {u_{k}} \left(\mathbf {v_{k}} ^{*}A\mathbf {u_{k}} \right)^{-1}\mathbf {v_{k}} ^{*}A$,

where $\mathbf {u_{k}} =\left(u_{0},u_{1},\dots u_{k-1}\right)$ and $\mathbf {v_{k}} =\left(v_{0},v_{1},\dots v_{k-1}\right)$. These related projections may be iterated themselves, as

$P_{k+1}=P_{k}+\left(1-P_{k}\right)u_{k}\otimes {v_{k}^{*}A\left(1-P_{k}\right) \over v_{k}^{*}A\left(1-P_{k}\right)u_{k}}$.

The new directions $d_{k}:=\left(1-P_{k}\right)u_{k}$ and $f_{k}:=\left(A\left(1-P_{k}\right)A^{-1}\right)^{*}v_{k}$ then are orthogonal to the residuums: $v_{i}^{*}r_{k}=f_{i}^{*}r_{k}=0$ and $s_{k}^{*}u_{j}=s_{k}^{*}d_{j}=0$, which themselves satisfy $r_{k}=A\left(1-P_{k}\right)A^{-1}r_{j}$ and $s_{k}=\left(1-P_{k}\right)^{*}s_{j}$($i,j<k$).

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

$u_{k}:=M^{-1}r_{k}$ and $v_{k}:=M^{-*}s_{k}$.

This special choice allows to avoid direct evaluations of $P_{k}$ and $A^{-1}$, and subsequently leads to the algorithm as stated above.

Properties

• If $A=A^{*}$ is self-adjoint, $y_{0}=x_{0}$ and $c=b$, then $r_{k}=s_{k}$, $d_{k}=f_{k}$, and the conjugate gradient method produces the same sequence $x_{k}=y_{k}$.

• In finite dimensions $x_{n}=A^{-1}b$, at the latest when $P_{n}=1$: The biconjugate gradient method returns the exact solution after iterating the full space and thus is a direct method.

• The sequences produced by the algorithm are biorthogonal: $f_{i}^{*}Ad_{j}=0$ and $s_{i}^{*}M^{-1}r_{j}=0$ for $i\neq j$.

• if $p_{j'}$ is a polynomial with $deg\left(p_{j'}\right)+j<k$, then $s_{k}^{*}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_{i}^{*}p_{i'}\left(AM^{-1}\right)r_{k}=0$.