Biconjugate gradient method

From Wikipedia, the free encyclopedia

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 $<var>A</var> *$ .

The algorithm

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^{*}-x_{0}^{*}\,A^{T}$
$p_{0}\leftarrow M^{{-1}}r_{0}\,$
$p_{0}^{*}\leftarrow r_{0}^{*}M^{{-1}}\,$
for do
1. $\alpha _{k}\leftarrow {r_{k}^{*}M^{{-1}}r_{k} \over p_{k}^{*}Ap_{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 Ap_{k}\,$
5. $r_{{k+1}}^{*}\leftarrow r_{k}^{*}-\overline {\alpha _{k}}\cdot p_{k}^{*}\,A$
6. $\beta _{k}\leftarrow {r_{{k+1}}^{*}M^{{-1}}r_{{k+1}} \over 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-Ax_{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 $\overline {\alpha }$ is the complex conjugate.

Unpreconditioned version of the algorithm

Choose initial guess $x_{0}\,$ ,
$r_{0}\leftarrow b-A\,x_{0}\,$
${\hat {r}}_{0}\leftarrow {\hat {b}}-{\hat {x}}_{0}A^{T}$
$p_{0}\leftarrow r_{0}\,$
${\hat {p}}_{0}\leftarrow {\hat {r}}_{0}\,$
for do
1. $\alpha _{k}\leftarrow {{\hat {r}}_{k}r_{k} \over {\hat {p}}_{k}Ap_{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 Ap_{k}\,$
5. ${\hat {r}}_{{k+1}}\leftarrow {\hat {r}}_{k}-\alpha _{k}\cdot {\hat {p}}_{k}A^{T}$
6. $\beta _{k}\leftarrow {{\hat {r}}_{{k+1}}r_{{k+1}} \over {\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^{[citation needed]} (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-Ax_{j}\right),$

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

where $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,$

with

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

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

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}}\left(Ax_{k}-b\right)$ , where

$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

$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 $<var>A</var> - 1$ 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}^{*}Ap_{j}=r_{i}^{*}M^{{-1}}r_{j}=0$ for $i\neq j$ .

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

References

Fletcher, R. (1976). "Conjugate gradient methods for indefinite systems". In Watson, G. Alistair. 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.
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

v t e Numerical linear algebra

Key concepts	Floating point Numerical stability BLAS

Problems	Matrix multiplication Matrix decompositions Linear equations Sparse problems

Hardware	CPU cache TLB Cache-oblivious algorithm SIMD Multiprocessing

Software	Specialized libraries General purpose software

This article is issued from Wikipedia. The text is available under the Creative Commons Attribution/Share Alike; additional terms may apply for the media files.