Un article de Wikipédia, l'encyclopédie libre.
En mathématiques , plus spécifiquement en analyse numérique , la méthode du gradient biconjugué est un algorithme permettant de résoudre un système d'équations linéaires
A
x
=
b
.
{\displaystyle Ax=b.\,}
Contrairement à la méthode du gradient conjugué , cet algorithme ne nécessite pas que la matrice
A
{\displaystyle A}
soit auto-adjointe, en revanche, la méthode requiert des multiplications par la matrice adjointe
A
∗
{\displaystyle A^{*}}
.
Choisir
x
0
{\displaystyle x_{0}}
,
y
0
{\displaystyle y_{0}}
, un préconditionneur régulier
M
{\displaystyle M}
(on utilise fréquemment
M
−
1
=
1
{\displaystyle M^{-1}=1}
) et
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}}
et
s
k
=
c
−
A
∗
y
k
{\displaystyle s_{k}=c-A^{*}y_{k}}
sont le résidus);
β
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}\,}
.
La méthode est numériquement instable , mais on y remédie par la méthode stabilisée du gradient biconjugué (en) , et elle reste très importante du point de vue théorique : on définit l'itération par
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)}
et
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}
) en utilisant les projections suivantes :
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}
,
Avec
u
k
=
(
u
0
,
u
1
,
…
u
k
−
1
)
{\displaystyle \mathbf {u_{k}} =\left(u_{0},u_{1},\dots u_{k-1}\right)}
et
v
k
=
(
v
0
,
v
1
,
…
v
k
−
1
)
{\displaystyle \mathbf {v_{k}} =\left(v_{0},v_{1},\dots v_{k-1}\right)}
. On peut itérer les projections elles-mêmes, comme
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}}}
.
Les nouvelles directions de descente
d
k
:=
(
1
−
P
k
)
u
k
{\displaystyle d_{k}:=\left(1-P_{k}\right)u_{k}}
et
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}}
sont alors orthogonales aux résidus :
v
i
∗
r
k
=
f
i
∗
r
k
=
0
{\displaystyle v_{i}^{*}r_{k}=f_{i}^{*}r_{k}=0}
et
s
k
∗
u
j
=
s
k
∗
d
j
=
0
{\displaystyle s_{k}^{*}u_{j}=s_{k}^{*}d_{j}=0}
, qui satisfont aux mêmes
r
k
=
A
(
1
−
P
k
)
A
−
1
r
j
{\displaystyle r_{k}=A\left(1-P_{k}\right)A^{-1}r_{j}}
et
s
k
=
(
1
−
P
k
)
∗
s
j
{\displaystyle s_{k}=\left(1-P_{k}\right)^{*}s_{j}}
(
i
,
j
<
k
{\displaystyle i,j<k}
).
La méthode du gradient biconjugué propose alors le choix suivant :
u
k
:=
M
−
1
r
k
{\displaystyle u_{k}:=M^{-1}r_{k}}
et
v
k
:=
M
−
∗
s
k
{\displaystyle v_{k}:=M^{-*}s_{k}}
.
Ce choix particulier permet alors d'éviter une évaluation directe de
P
k
{\displaystyle P_{k}}
et
A
−
1
{\displaystyle A^{-1}}
, et donc augmenter la vitesse d'exécution de l'algorithme.
Si
A
=
A
∗
{\displaystyle A=A^{*}}
est auto-adjointe ,
y
0
=
x
0
{\displaystyle y_{0}=x_{0}}
et
c
=
b
{\displaystyle c=b}
, donc
r
k
=
s
k
{\displaystyle r_{k}=s_{k}}
,
d
k
=
f
k
{\displaystyle d_{k}=f_{k}}
, et la méthode du gradient conjugué produit la même suite
x
k
=
y
k
{\displaystyle x_{k}=y_{k}}
.
En dimensions finies
x
n
=
A
−
1
b
{\displaystyle x_{n}=A^{-1}b}
, au plus tard quand
P
n
=
1
{\displaystyle P_{n}=1}
: La méthode du gradient biconjugué rend la solution exacte après avoir parcouru tout l'espace et est donc une méthode directe.
La suite produite par l'algorithme est biorthogonale (en) :
f
i
∗
A
d
j
=
0
{\displaystyle f_{i}^{*}Ad_{j}=0}
et
s
i
∗
M
−
1
r
j
=
0
{\displaystyle s_{i}^{*}M^{-1}r_{j}=0}
pour
i
≠
j
{\displaystyle i\neq j}
.
SI
p
j
′
{\displaystyle p_{j'}}
est un polynôme avec
d
e
g
(
p
j
′
)
+
j
<
k
{\displaystyle deg\left(p_{j'}\right)+j<k}
, alors
s
k
∗
p
j
′
(
M
−
1
A
)
u
j
=
0
{\displaystyle s_{k}^{*}p_{j'}\left(M^{-1}A\right)u_{j}=0}
. L'algorithme est donc composé de projections sur des sous-espaces de Krylov ;
SI
p
i
′
{\displaystyle p_{i'}}
est un polynôme avec
i
+
d
e
g
(
p
i
′
)
<
k
{\displaystyle i+deg\left(p_{i'}\right)<k}
, alors
v
i
∗
p
i
′
(
A
M
−
1
)
r
k
=
0
{\displaystyle v_{i}^{*}p_{i'}\left(AM^{-1}\right)r_{k}=0}
.