In mathematics, more specifically in numerical linear algebra, the biconjugate gradient method is an algorithm to solve systems of linear equations
![{\displaystyle Ax=b.\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/83f1eec82741fe1371d5b10a1d8eb2360367c396)
Unlike the conjugate gradient method, this algorithm does not require the matrix
to be self-adjoint, but instead one needs to perform multiplications by the conjugate transpose A*.
The algorithm
- Choose initial guess
, two other vectors
and
and a preconditioner ![{\displaystyle M\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/aa66ea8156203e93f2dca12aca24538c2bdce761)
![{\displaystyle r_{0}\leftarrow b-A\,x_{0}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/645a7895cbe4795127dee1886a3f9709ea499b4e)
![{\displaystyle r_{0}^{*}\leftarrow b^{*}-x_{0}^{*}\,A}](https://wikimedia.org/api/rest_v1/media/math/render/svg/82f4c8d8f088924f16a9d01f0ffc5085a58c55d2)
![{\displaystyle p_{0}\leftarrow M^{-1}r_{0}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e51a2568076be0bed2cf9b2d3568f75a42266e20)
![{\displaystyle p_{0}^{*}\leftarrow r_{0}^{*}M^{-1}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5b425310928ae31fe214f183cb1153161eb4d3b1)
- for
do
![{\displaystyle \alpha _{k}\leftarrow {r_{k}^{*}M^{-1}r_{k} \over p_{k}^{*}Ap_{k}}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/022ee0d3fb5d2c6b8a2d9edf2a87577d4bc8ce3d)
![{\displaystyle x_{k+1}\leftarrow x_{k}+\alpha _{k}\cdot p_{k}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9ff6f21ecf9067251bc5e1f8b7089b064b87cfe7)
![{\displaystyle x_{k+1}^{*}\leftarrow x_{k}^{*}+{\overline {\alpha _{k}}}\cdot p_{k}^{*}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2d243b204d2d668bb2f76ca6dc48224f4eb3e352)
![{\displaystyle r_{k+1}\leftarrow r_{k}-\alpha _{k}\cdot Ap_{k}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d8f13fe1a781b1bb892f5a08f714b72fa0afef7e)
![{\displaystyle r_{k+1}^{*}\leftarrow r_{k}^{*}-{\overline {\alpha _{k}}}\cdot p_{k}^{*}\,A}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1041a8c48c2977741fe041a4bf66979f31d7f90a)
![{\displaystyle \beta _{k}\leftarrow {r_{k+1}^{*}M^{-1}r_{k+1} \over r_{k}^{*}M^{-1}r_{k}}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8725a7130ba943dbdfb5827a61fbf185c61f259c)
![{\displaystyle p_{k+1}\leftarrow M^{-1}r_{k+1}+\beta _{k}\cdot p_{k}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ba8c847bf583663824419fee8423f91f8983579c)
![{\displaystyle p_{k+1}^{*}\leftarrow r_{k+1}^{*}M^{-1}+{\overline {\beta _{k}}}\cdot p_{k}^{*}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d7f35e5ff4788071d771756edc554d9cf3fce1d2)
In the above formulation, the computed
and
satisfy
![{\displaystyle r_{k}=b-Ax_{k},\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/daac71fa9ccf99317571262f961e67d65929c99f)
![{\displaystyle r_{k}^{*}=b^{*}-x_{k}^{*}\,A}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5f613de2813f9abaa07503d6fcd87d74b1930bb5)
and thus are the respective residuals corresponding to
and
, as approximate solutions to the systems
![{\displaystyle Ax=b,\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ecff317a7ef55f43478aed6543005e9c6f66e49e)
![{\displaystyle x^{*}\,A=b^{*}\,;}](https://wikimedia.org/api/rest_v1/media/math/render/svg/733e8d5f6d7951d4329dd25e5514003f01566742)
is the adjoint, and
is the complex conjugate.
Unpreconditioned version of the algorithm
- Choose initial guess
,
![{\displaystyle r_{0}\leftarrow b-A\,x_{0}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/645a7895cbe4795127dee1886a3f9709ea499b4e)
![{\displaystyle {\hat {r}}_{0}\leftarrow {\hat {b}}-{\hat {x}}_{0}A}](https://wikimedia.org/api/rest_v1/media/math/render/svg/07074e54734f1ec6311229661431b557c29d4535)
![{\displaystyle p_{0}\leftarrow r_{0}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/778453e707906960f01eddc0e9761c1c69412c7b)
![{\displaystyle {\hat {p}}_{0}\leftarrow {\hat {r}}_{0}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e62055f08cc95ec1c86be57e28618768e3005358)
- for
do
![{\displaystyle \alpha _{k}\leftarrow {{\hat {r}}_{k}r_{k} \over {\hat {p}}_{k}Ap_{k}}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e9003309cd49d353c27ef22014807af56bdf3184)
![{\displaystyle x_{k+1}\leftarrow x_{k}+\alpha _{k}\cdot p_{k}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9ff6f21ecf9067251bc5e1f8b7089b064b87cfe7)
![{\displaystyle {\hat {x}}_{k+1}\leftarrow {\hat {x}}_{k}+\alpha _{k}\cdot {\hat {p}}_{k}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1726f4f79db7d4106384fb1ec6b2366bd42dc017)
![{\displaystyle r_{k+1}\leftarrow r_{k}-\alpha _{k}\cdot Ap_{k}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d8f13fe1a781b1bb892f5a08f714b72fa0afef7e)
![{\displaystyle {\hat {r}}_{k+1}\leftarrow {\hat {r}}_{k}-\alpha _{k}\cdot {\hat {p}}_{k}A}](https://wikimedia.org/api/rest_v1/media/math/render/svg/622c2e95caccabfcfd131e8e192d40ad484934f0)
![{\displaystyle \beta _{k}\leftarrow {{\hat {r}}_{k+1}r_{k+1} \over {\hat {r}}_{k}r_{k}}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/416bcfc683e0663d35b341ee7458a65e280ae377)
![{\displaystyle p_{k+1}\leftarrow r_{k+1}+\beta _{k}\cdot p_{k}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/beabfbd566b77dc65915c00efcc015a1424ced49)
![{\displaystyle {\hat {p}}_{k+1}\leftarrow {\hat {r}}_{k+1}+\beta _{k}\cdot {\hat {p}}_{k}\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2147dab1f496e9f9fcc299991491b2f975000a6b)
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
![{\displaystyle x_{k}:=x_{j}+P_{k}A^{-1}\left(b-Ax_{j}\right),}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1614ba2525d0ae8cd250e28654e6dc715b360983)
![{\displaystyle x_{k}^{*}:=x_{j}^{*}+\left(b^{*}-x_{j}^{*}A\right)P_{k}A^{-1},}](https://wikimedia.org/api/rest_v1/media/math/render/svg/493ffbdc5d3265aa409b210a32f87e82408c521b)
where
using the related projection
![{\displaystyle P_{k}:=\mathbf {u} _{k}\left(\mathbf {v} _{k}^{*}A\mathbf {u} _{k}\right)^{-1}\mathbf {v} _{k}^{*}A,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b2b1105cf87f2ab9aa9af9addee50086f124ab79)
with
![{\displaystyle \mathbf {u} _{k}=\left[u_{0},u_{1},\dots ,u_{k-1}\right],}](https://wikimedia.org/api/rest_v1/media/math/render/svg/356e9dd32012d4b25f4a3c78179554570098e78e)
![{\displaystyle \mathbf {v} _{k}=\left[v_{0},v_{1},\dots ,v_{k-1}\right].}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0f9338d3797abffd4007c6c2ab27aaed7e04e893)
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}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4990ef33691eb8ff46a537f82933b7882fa03074)
A relation to Quasi-Newton methods is given by
and
, 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}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/43da44535b6971984e59cde8c6460f45f508b0f6)
The new directions
![{\displaystyle p_{k}=\left(1-P_{k}\right)u_{k},}](https://wikimedia.org/api/rest_v1/media/math/render/svg/38bde397f420024d2c998f9a9fed2d2cb02d32cc)
![{\displaystyle p_{k}^{*}=v_{k}^{*}A\left(1-P_{k}\right)A^{-1}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8e67e0cd8c8333ecb103d37717f1980bd1bcb8ef)
are then orthogonal to the residuals:
![{\displaystyle v_{i}^{*}r_{k}=p_{i}^{*}r_{k}=0,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/13502cd02064035a0b797d0e175c324946ab2c6f)
![{\displaystyle r_{k}^{*}u_{j}=r_{k}^{*}p_{j}=0,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fc214533995a96ceb0121ea71fc608a1b1bb29e1)
which themselves satisfy
![{\displaystyle r_{k}=A\left(1-P_{k}\right)A^{-1}r_{j},}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8c99206844f87acb97b01fdf6b4667037b86cd9e)
![{\displaystyle r_{k}^{*}=r_{j}^{*}\left(1-P_{k}\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/008bc80cb97c2a3c377e8718b28be14f7832ae23)
where
.
The biconjugate gradient method now makes a special choice and uses the setting
![{\displaystyle u_{k}=M^{-1}r_{k},\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9cb5a9b9bd2916b3fb86e376536bee0feecf512f)
![{\displaystyle v_{k}^{*}=r_{k}^{*}\,M^{-1}.\,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a405678be1ba07be23f3f1457bec979236688524)
With this particular choice, explicit evaluations of
and A−1 are avoided, and the algorithm takes the form stated above.
Properties
- If
is self-adjoint,
and
, then
,
, and the conjugate gradient method produces the same sequence
at half the computational cost.
- The sequences produced by the algorithm are biorthogonal, i.e.,
for
.
- if
is a polynomial with
, then
. The algorithm thus produces projections onto the Krylov subspace.
- if
is a polynomial with
, then
.
See also
References
|
---|
Key concepts | |
---|
Problems | |
---|
Hardware | |
---|
Software | |
---|