Biconjugate gradient method

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

${\displaystyle Ax=b.,円}${\displaystyle Ax=b.,円}

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

1. Choose initial guess ${\displaystyle x_{0},円}${\displaystyle x_{0},円}, two other vectors ${\displaystyle x_{0}^{*}}${\displaystyle x_{0}^{*}} and ${\displaystyle b^{*},円}${\displaystyle b^{*},円} and a preconditioner ${\displaystyle M,円}${\displaystyle M,円}
2. ${\displaystyle r_{0}\leftarrow b-A,円x_{0},円}${\displaystyle r_{0}\leftarrow b-A,円x_{0},円}
3. ${\displaystyle r_{0}^{*}\leftarrow b^{*}-x_{0}^{*},円A^{*}}${\displaystyle r_{0}^{*}\leftarrow b^{*}-x_{0}^{*},円A^{*}}
4. ${\displaystyle p_{0}\leftarrow M^{-1}r_{0},円}${\displaystyle p_{0}\leftarrow M^{-1}r_{0},円}
5. ${\displaystyle p_{0}^{*}\leftarrow r_{0}^{*}M^{-1},円}${\displaystyle p_{0}^{*}\leftarrow r_{0}^{*}M^{-1},円}
6. for ${\displaystyle k=0,1,\ldots }${\displaystyle k=0,1,\ldots } do
   1. ${\displaystyle \alpha _{k}\leftarrow {r_{k}^{*}M^{-1}r_{k} \over p_{k}^{*}Ap_{k}},円}${\displaystyle \alpha _{k}\leftarrow {r_{k}^{*}M^{-1}r_{k} \over p_{k}^{*}Ap_{k}},円}
   2. ${\displaystyle x_{k+1}\leftarrow x_{k}+\alpha _{k}\cdot p_{k},円}${\displaystyle x_{k+1}\leftarrow x_{k}+\alpha _{k}\cdot p_{k},円}
   3. ${\displaystyle x_{k+1}^{*}\leftarrow x_{k}^{*}+{\overline {\alpha _{k}}}\cdot p_{k}^{*},円}${\displaystyle x_{k+1}^{*}\leftarrow x_{k}^{*}+{\overline {\alpha _{k}}}\cdot p_{k}^{*},円}
   4. ${\displaystyle r_{k+1}\leftarrow r_{k}-\alpha _{k}\cdot Ap_{k},円}${\displaystyle r_{k+1}\leftarrow r_{k}-\alpha _{k}\cdot Ap_{k},円}
   5. ${\displaystyle r_{k+1}^{*}\leftarrow r_{k}^{*}-{\overline {\alpha _{k}}}\cdot p_{k}^{*},円A^{*}}${\displaystyle r_{k+1}^{*}\leftarrow r_{k}^{*}-{\overline {\alpha _{k}}}\cdot p_{k}^{*},円A^{*}}
   6. ${\displaystyle \beta _{k}\leftarrow {r_{k+1}^{*}M^{-1}r_{k+1} \over r_{k}^{*}M^{-1}r_{k}},円}${\displaystyle \beta _{k}\leftarrow {r_{k+1}^{*}M^{-1}r_{k+1} \over r_{k}^{*}M^{-1}r_{k}},円}
   7. ${\displaystyle p_{k+1}\leftarrow M^{-1}r_{k+1}+\beta _{k}\cdot p_{k},円}${\displaystyle p_{k+1}\leftarrow M^{-1}r_{k+1}+\beta _{k}\cdot p_{k},円}
   8. ${\displaystyle p_{k+1}^{*}\leftarrow r_{k+1}^{*}M^{-1}+{\overline {\beta _{k}}}\cdot p_{k}^{*},円}${\displaystyle p_{k+1}^{*}\leftarrow r_{k+1}^{*}M^{-1}+{\overline {\beta _{k}}}\cdot p_{k}^{*},円}

In the above formulation, the computed ${\displaystyle r_{k},円}${\displaystyle r_{k},円} and ${\displaystyle r_{k}^{*}}${\displaystyle r_{k}^{*}} satisfy

${\displaystyle r_{k}=b-Ax_{k},,円}${\displaystyle r_{k}=b-Ax_{k},,円}

${\displaystyle r_{k}^{*}=b^{*}-x_{k}^{*},円A^{*}}${\displaystyle r_{k}^{*}=b^{*}-x_{k}^{*},円A^{*}}

and thus are the respective residuals corresponding to ${\displaystyle x_{k},円}${\displaystyle x_{k},円} and ${\displaystyle x_{k}^{*}}${\displaystyle x_{k}^{*}}, as approximate solutions to the systems

${\displaystyle Ax=b,,円}${\displaystyle Ax=b,,円}

${\displaystyle x^{*},円A^{*}=b^{*},円;}${\displaystyle x^{*},円A^{*}=b^{*},円;}

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

1. Choose initial guess ${\displaystyle x_{0},円}${\displaystyle x_{0},円},
2. ${\displaystyle r_{0}\leftarrow b-A,円x_{0},円}${\displaystyle r_{0}\leftarrow b-A,円x_{0},円}
3. ${\displaystyle {\hat {r}}_{0}\leftarrow {\hat {b}}-{\hat {x}}_{0}A^{*}}${\displaystyle {\hat {r}}_{0}\leftarrow {\hat {b}}-{\hat {x}}_{0}A^{*}}
4. ${\displaystyle p_{0}\leftarrow r_{0},円}${\displaystyle p_{0}\leftarrow r_{0},円}
5. ${\displaystyle {\hat {p}}_{0}\leftarrow {\hat {r}}_{0},円}${\displaystyle {\hat {p}}_{0}\leftarrow {\hat {r}}_{0},円}
6. for ${\displaystyle k=0,1,\ldots }${\displaystyle k=0,1,\ldots } do
   1. ${\displaystyle \alpha _{k}\leftarrow {{\hat {r}}_{k}r_{k} \over {\hat {p}}_{k}Ap_{k}},円}${\displaystyle \alpha _{k}\leftarrow {{\hat {r}}_{k}r_{k} \over {\hat {p}}_{k}Ap_{k}},円}
   2. ${\displaystyle x_{k+1}\leftarrow x_{k}+\alpha _{k}\cdot p_{k},円}${\displaystyle x_{k+1}\leftarrow x_{k}+\alpha _{k}\cdot p_{k},円}
   3. ${\displaystyle {\hat {x}}_{k+1}\leftarrow {\hat {x}}_{k}+\alpha _{k}\cdot {\hat {p}}_{k},円}${\displaystyle {\hat {x}}_{k+1}\leftarrow {\hat {x}}_{k}+\alpha _{k}\cdot {\hat {p}}_{k},円}
   4. ${\displaystyle r_{k+1}\leftarrow r_{k}-\alpha _{k}\cdot Ap_{k},円}${\displaystyle r_{k+1}\leftarrow r_{k}-\alpha _{k}\cdot Ap_{k},円}
   5. ${\displaystyle {\hat {r}}_{k+1}\leftarrow {\hat {r}}_{k}-\alpha _{k}\cdot {\hat {p}}_{k}A^{*}}${\displaystyle {\hat {r}}_{k+1}\leftarrow {\hat {r}}_{k}-\alpha _{k}\cdot {\hat {p}}_{k}A^{*}}
   6. ${\displaystyle \beta _{k}\leftarrow {{\hat {r}}_{k+1}r_{k+1} \over {\hat {r}}_{k}r_{k}},円}${\displaystyle \beta _{k}\leftarrow {{\hat {r}}_{k+1}r_{k+1} \over {\hat {r}}_{k}r_{k}},円}
   7. ${\displaystyle p_{k+1}\leftarrow r_{k+1}+\beta _{k}\cdot p_{k},円}${\displaystyle p_{k+1}\leftarrow r_{k+1}+\beta _{k}\cdot p_{k},円}
   8. ${\displaystyle {\hat {p}}_{k+1}\leftarrow {\hat {r}}_{k+1}+\beta _{k}\cdot {\hat {p}}_{k},円}${\displaystyle {\hat {p}}_{k+1}\leftarrow {\hat {r}}_{k+1}+\beta _{k}\cdot {\hat {p}}_{k},円}

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

${\displaystyle x_{k}:=x_{j}+P_{k}A^{-1}\left(b-Ax_{j}\right),}${\displaystyle x_{k}:=x_{j}+P_{k}A^{-1}\left(b-Ax_{j}\right),}

${\displaystyle x_{k}^{*}:=x_{j}^{*}+\left(b^{*}-x_{j}^{*}A\right)P_{k}A^{-1},}${\displaystyle x_{k}^{*}:=x_{j}^{*}+\left(b^{*}-x_{j}^{*}A\right)P_{k}A^{-1},}

where ${\displaystyle j<k}${\displaystyle j<k} using the related projection

${\displaystyle P_{k}:=\mathbf {u} _{k}\left(\mathbf {v} _{k}^{*}A\mathbf {u} _{k}\right)^{-1}\mathbf {v} _{k}^{*}A,}${\displaystyle P_{k}:=\mathbf {u} _{k}\left(\mathbf {v} _{k}^{*}A\mathbf {u} _{k}\right)^{-1}\mathbf {v} _{k}^{*}A,}

with

${\displaystyle \mathbf {u} _{k}=\left[u_{0},u_{1},\dots ,u_{k-1}\right],}${\displaystyle \mathbf {u} _{k}=\left[u_{0},u_{1},\dots ,u_{k-1}\right],}

${\displaystyle \mathbf {v} _{k}=\left[v_{0},v_{1},\dots ,v_{k-1}\right].}${\displaystyle \mathbf {v} _{k}=\left[v_{0},v_{1},\dots ,v_{k-1}\right].}

These related projections may be iterated themselves as

${\displaystyle 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}}.}${\displaystyle 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}}.}

A relation to Quasi-Newton methods is given by ${\displaystyle P_{k}=A_{k}^{-1}A}${\displaystyle P_{k}=A_{k}^{-1}A} and ${\displaystyle x_{k+1}=x_{k}-A_{k+1}^{-1}\left(Ax_{k}-b\right)}${\displaystyle x_{k+1}=x_{k}-A_{k+1}^{-1}\left(Ax_{k}-b\right)}, where

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

The new directions

${\displaystyle p_{k}=\left(1-P_{k}\right)u_{k},}${\displaystyle p_{k}=\left(1-P_{k}\right)u_{k},}

${\displaystyle p_{k}^{*}=v_{k}^{*}A\left(1-P_{k}\right)A^{-1}}${\displaystyle p_{k}^{*}=v_{k}^{*}A\left(1-P_{k}\right)A^{-1}}

are then orthogonal to the residuals:

${\displaystyle v_{i}^{*}r_{k}=p_{i}^{*}r_{k}=0,}${\displaystyle v_{i}^{*}r_{k}=p_{i}^{*}r_{k}=0,}

${\displaystyle r_{k}^{*}u_{j}=r_{k}^{*}p_{j}=0,}${\displaystyle r_{k}^{*}u_{j}=r_{k}^{*}p_{j}=0,}

which themselves satisfy

${\displaystyle r_{k}=A\left(1-P_{k}\right)A^{-1}r_{j},}${\displaystyle r_{k}=A\left(1-P_{k}\right)A^{-1}r_{j},}

${\displaystyle r_{k}^{*}=r_{j}^{*}\left(1-P_{k}\right)}${\displaystyle r_{k}^{*}=r_{j}^{*}\left(1-P_{k}\right)}

where ${\displaystyle i,j<k}${\displaystyle i,j<k}.

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

${\displaystyle u_{k}=M^{-1}r_{k},,円}${\displaystyle u_{k}=M^{-1}r_{k},,円}

${\displaystyle v_{k}^{*}=r_{k}^{*},円M^{-1}.,円}${\displaystyle v_{k}^{*}=r_{k}^{*},円M^{-1}.,円}

With this particular choice, explicit evaluations of ${\displaystyle P_{k}}${\displaystyle P_{k}} and A⁻¹ are avoided, and the algorithm takes the form stated above.

• If ${\displaystyle A=A^{*},円}${\displaystyle A=A^{*},円} is self-adjoint, ${\displaystyle x_{0}^{*}=x_{0}}${\displaystyle x_{0}^{*}=x_{0}} and ${\displaystyle b^{*}=b}${\displaystyle b^{*}=b}, then ${\displaystyle r_{k}=r_{k}^{*}}${\displaystyle r_{k}=r_{k}^{*}}, ${\displaystyle p_{k}=p_{k}^{*}}${\displaystyle p_{k}=p_{k}^{*}}, and the conjugate gradient method produces the same sequence ${\displaystyle x_{k}=x_{k}^{*}}${\displaystyle x_{k}=x_{k}^{*}} at half the computational cost.
• The sequences produced by the algorithm are biorthogonal, i.e., ${\displaystyle p_{i}^{*}Ap_{j}=r_{i}^{*}M^{-1}r_{j}=0}${\displaystyle p_{i}^{*}Ap_{j}=r_{i}^{*}M^{-1}r_{j}=0} for ${\displaystyle i\neq j}${\displaystyle i\neq j}.
• if ${\displaystyle P_{j'},円}${\displaystyle P_{j'},円} is a polynomial with ${\displaystyle \deg \left(P_{j'}\right)+j<k}${\displaystyle \deg \left(P_{j'}\right)+j<k}, then ${\displaystyle r_{k}^{*}P_{j'}\left(M^{-1}A\right)u_{j}=0}${\displaystyle r_{k}^{*}P_{j'}\left(M^{-1}A\right)u_{j}=0}. The algorithm thus produces projections onto the Krylov subspace.
• if ${\displaystyle P_{i'},円}${\displaystyle P_{i'},円} is a polynomial with ${\displaystyle i+\deg \left(P_{i'}\right)<k}${\displaystyle i+\deg \left(P_{i'}\right)<k}, then ${\displaystyle v_{i}^{*}P_{i'}\left(AM^{-1}\right)r_{k}=0}${\displaystyle v_{i}^{*}P_{i'}\left(AM^{-1}\right)r_{k}=0}.

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

• Fletcher, R. (1976). "Conjugate gradient methods for indefinite systems". In Watson, G. Alistair (ed.). Numerical analysis : proceedings of the Dundee Conference on Numerical Analysis. Lecture Notes in Mathematics. Vol. 506. Springer. pp. 73–89. doi:10.1007/BFb0080116. ISBN 978-3-540-07610-0.
• 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.

• Numerical linear algebra
• Gradient methods

• Articles with short description
• Short description is different from Wikidata
• Articles lacking in-text citations from September 2013
• All articles lacking in-text citations
• All articles with unsourced statements
• Articles with unsourced statements from September 2009