From Wikipedia, the free encyclopedia
In mathematics , more specifically in numerical analysis , the biconjugate gradient method is an algorithm to solve systems of linear equations
A
x
=
b
.
{\displaystyle Ax=b.\,}
Unlike the conjugate gradient method , this algorithm does not require the matrix
A
{\displaystyle A}
to be self-adjoint, but instead one needs to perform multiplications by the conjugate transpose
A
∗
{\displaystyle A^{*}}
.
The algorithm
Choose
x
0
{\displaystyle x_{0}}
,
y
0
{\displaystyle y_{0}}
, a regular preconditioner
M
{\displaystyle M}
(frequently
M
−
1
=
1
{\displaystyle M^{-1}=1}
is used) and
c
{\displaystyle c}
;
r
0
←
b
−
A
x
0
,
s
0
←
c
−
A
∗
y
0
{\displaystyle r_{0}\leftarrow b-Ax_{0},s_{0}\leftarrow c-A^{*}y_{0}\,}
;
d
0
←
M
−
1
r
0
,
f
0
←
M
−
∗
s
0
{\displaystyle d_{0}\leftarrow M^{-1}r_{0},f_{0}\leftarrow M^{-*}s_{0}\,}
;
for
k
=
0
,
1
,
…
{\displaystyle k=0,1,\dots \,}
do
α
k
←
s
k
∗
M
−
1
r
k
f
k
∗
A
d
k
{\displaystyle \alpha _{k}\leftarrow {s_{k}^{*}M^{-1}r_{k} \over f_{k}^{*}Ad_{k}}\,}
;
x
k
+
1
←
x
k
+
α
k
d
k
{\displaystyle x_{k+1}\leftarrow x_{k}+\alpha _{k}d_{k}\,}
(
y
k
+
1
←
y
k
+
α
k
¯
f
k
)
{\displaystyle \left(y_{k+1}\leftarrow y_{k}+{\overline {\alpha _{k}}}f_{k}\right)\,}
;
r
k
+
1
←
r
k
−
α
k
A
d
k
{\displaystyle r_{k+1}\leftarrow r_{k}-\alpha _{k}Ad_{k}\,}
,
s
k
+
1
←
s
k
−
α
k
¯
A
∗
f
k
{\displaystyle s_{k+1}\leftarrow s_{k}-{\overline {\alpha _{k}}}A^{*}f_{k}\,}
(
r
k
=
b
−
A
x
k
{\displaystyle r_{k}=b-Ax_{k}}
and
s
k
=
c
−
A
∗
y
k
{\displaystyle s_{k}=c-A^{*}y_{k}}
are the residuums);
β
k
←
s
k
+
1
∗
M
−
1
r
k
+
1
s
k
∗
M
−
1
r
k
{\displaystyle \beta _{k}\leftarrow {s_{k+1}^{*}M^{-1}r_{k+1} \over s_{k}^{*}M^{-1}r_{k}}\,}
;
d
k
+
1
←
M
−
1
r
k
+
1
+
β
k
d
k
{\displaystyle d_{k+1}\leftarrow M^{-1}r_{k+1}+\beta _{k}d_{k}\,}
,
f
k
+
1
←
M
−
∗
s
k
+
1
+
β
k
¯
f
k
{\displaystyle 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
(
b
−
A
x
j
)
{\displaystyle x_{k}:=x_{j}+P_{k}A^{-1}\left(b-Ax_{j}\right)}
and
y
k
:=
y
j
+
A
−
∗
P
k
∗
(
c
−
A
∗
y
j
)
{\displaystyle y_{k}:=y_{j}+A^{-*}P_{k}^{*}\left(c-A^{*}y_{j}\right)}
(
j
<
k
{\displaystyle j<k}
) using the related projection
P
k
:=
u
k
(
v
k
∗
A
u
k
)
−
1
v
k
∗
A
{\displaystyle P_{k}:=\mathbf {u_{k}} \left(\mathbf {v_{k}} ^{*}A\mathbf {u_{k}} \right)^{-1}\mathbf {v_{k}} ^{*}A}
,
where
u
k
=
(
u
0
,
u
1
,
…
u
k
−
1
)
{\displaystyle \mathbf {u_{k}} =\left(u_{0},u_{1},\dots u_{k-1}\right)}
and
v
k
=
(
v
0
,
v
1
,
…
v
k
−
1
)
{\displaystyle \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
+
(
1
−
P
k
)
u
k
⊗
v
k
∗
A
(
1
−
P
k
)
v
k
∗
A
(
1
−
P
k
)
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}}}
.
The new directions
d
k
:=
(
1
−
P
k
)
u
k
{\displaystyle d_{k}:=\left(1-P_{k}\right)u_{k}}
and
f
k
:=
(
A
(
1
−
P
k
)
A
−
1
)
∗
v
k
{\displaystyle 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
{\displaystyle v_{i}^{*}r_{k}=f_{i}^{*}r_{k}=0}
and
s
k
∗
u
j
=
s
k
∗
d
j
=
0
{\displaystyle s_{k}^{*}u_{j}=s_{k}^{*}d_{j}=0}
, which themselves satisfy
r
k
=
A
(
1
−
P
k
)
A
−
1
r
j
{\displaystyle r_{k}=A\left(1-P_{k}\right)A^{-1}r_{j}}
and
s
k
=
(
1
−
P
k
)
∗
s
j
{\displaystyle s_{k}=\left(1-P_{k}\right)^{*}s_{j}}
(
i
,
j
<
k
{\displaystyle i,j<k}
).
The biconjugate gradient method now makes a special choice and uses the setting
u
k
:=
M
−
1
r
k
{\displaystyle u_{k}:=M^{-1}r_{k}}
and
v
k
:=
M
−
∗
s
k
{\displaystyle v_{k}:=M^{-*}s_{k}}
.
This special choice allows to avoid direct evaluations of
P
k
{\displaystyle P_{k}}
and
A
−
1
{\displaystyle A^{-1}}
, and subsequently leads to the algorithm as stated above.
Properties
If
A
=
A
∗
{\displaystyle A=A^{*}}
is self-adjoint ,
y
0
=
x
0
{\displaystyle y_{0}=x_{0}}
and
c
=
b
{\displaystyle c=b}
, then
r
k
=
s
k
{\displaystyle r_{k}=s_{k}}
,
d
k
=
f
k
{\displaystyle d_{k}=f_{k}}
, and the conjugate gradient method produces the same sequence
x
k
=
y
k
{\displaystyle x_{k}=y_{k}}
.
In finite dimensions
x
n
=
A
−
1
b
{\displaystyle x_{n}=A^{-1}b}
, at the latest when
P
n
=
1
{\displaystyle 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
∗
A
d
j
=
0
{\displaystyle f_{i}^{*}Ad_{j}=0}
and
s
i
∗
M
−
1
r
j
=
0
{\displaystyle s_{i}^{*}M^{-1}r_{j}=0}
for
i
≠
j
{\displaystyle i\neq j}
.
if
p
j
′
{\displaystyle p_{j'}}
is a polynomial with
d
e
g
(
p
j
′
)
+
j
<
k
{\displaystyle deg\left(p_{j'}\right)+j<k}
, then
s
k
∗
p
j
′
(
M
−
1
A
)
u
j
=
0
{\displaystyle s_{k}^{*}p_{j'}\left(M^{-1}A\right)u_{j}=0}
. The algorithm thus produces projections onto the Krylov subspace ;
if
p
i
′
{\displaystyle p_{i'}}
is a polynomial with
i
+
d
e
g
(
p
i
′
)
<
k
{\displaystyle i+deg\left(p_{i'}\right)<k}
, then
v
i
∗
p
i
′
(
A
M
−
1
)
r
k
=
0
{\displaystyle v_{i}^{*}p_{i'}\left(AM^{-1}\right)r_{k}=0}
.