稳定双共轭梯度法

算法步骤

无预处理稳定双共轭梯度法

1. ${\displaystyle {\boldsymbol {r}}_{0}={\boldsymbol {b}}-{\boldsymbol {Ax}}}$
2. 任意选择向量 ${\displaystyle {\boldsymbol {\hat {r}}}_{0}\;}$ 使得 ${\displaystyle ({\boldsymbol {\hat {r}}}_{0},{\boldsymbol {r}}_{0})\neq 0\;}$，例如，${\displaystyle {\boldsymbol {\hat {r}}}_{0}={\boldsymbol {r}}_{0}\;}$
3. ${\displaystyle \rho _{0}=\alpha =\omega _{0}=1\;}$
4. ${\displaystyle {\boldsymbol {v}}_{0}={\boldsymbol {p}}_{0}={\boldsymbol {0}}}$
5. ${\displaystyle i=1,2,3,\ldots }$
1. ${\displaystyle \rho _{i}=({\boldsymbol {\hat {r}}}_{0},{\boldsymbol {r}}_{i-1})\;}$
2. ${\displaystyle \beta =(\rho _{i}/\rho _{i-1})(\alpha /\omega _{i-1})\;}$
3. ${\displaystyle {\boldsymbol {p}}_{i}={\boldsymbol {r}}_{i-1}+\beta ({\boldsymbol {p}}_{i-1}-\omega _{i-1}{\boldsymbol {v}}_{i-1})}$
4. ${\displaystyle {\boldsymbol {v}}_{i}={\boldsymbol {Ap}}_{i}}$
5. ${\displaystyle \alpha =\rho _{i}/({\boldsymbol {\hat {r}}}_{0},{\boldsymbol {v}}_{i})\;}$
6. ${\displaystyle {\boldsymbol {s}}={\boldsymbol {r}}_{i-1}-\alpha {\boldsymbol {v}}_{i}}$
7. ${\displaystyle {\boldsymbol {t}}={\boldsymbol {As}}}$
8. ${\displaystyle \omega _{i}=({\boldsymbol {t}},{\boldsymbol {s}})/({\boldsymbol {t}},{\boldsymbol {t}})}$
9. ${\displaystyle {\boldsymbol {x}}_{i}={\boldsymbol {x}}_{i-1}+\alpha {\boldsymbol {p}}_{i}+\omega _{i}{\boldsymbol {s}}}$
10. ${\displaystyle {\boldsymbol {x}}_{i}}$ 足够精确则退出
11. ${\displaystyle {\boldsymbol {r}}_{i}={\boldsymbol {s}}-\omega _{i}{\boldsymbol {t}}}$

预处理稳定双共轭梯度法

1. ${\displaystyle {\boldsymbol {r}}_{0}={\boldsymbol {b}}-{\boldsymbol {Ax}}}$
2. 任意选择向量 ${\displaystyle {\boldsymbol {\hat {r}}}_{0}\;}$ 使得 ${\displaystyle ({\boldsymbol {\hat {r}}}_{0},{\boldsymbol {r}}_{0})\neq 0\;}$，例如，${\displaystyle {\boldsymbol {\hat {r}}}_{0}={\boldsymbol {r}}_{0}\;}$
3. ${\displaystyle \rho _{0}=\alpha =\omega _{0}=1\;}$
4. ${\displaystyle {\boldsymbol {v}}_{0}={\boldsymbol {p}}_{0}={\boldsymbol {0}}}$
5. ${\displaystyle i=1,2,3,\ldots }$
1. ${\displaystyle \rho _{i}=({\boldsymbol {\hat {r}}}_{0},{\boldsymbol {r}}_{i-1})\;}$
2. ${\displaystyle \beta =(\rho _{i}/\rho _{i-1})(\alpha /\omega _{i-1})\;}$
3. ${\displaystyle {\boldsymbol {p}}_{i}={\boldsymbol {r}}_{i-1}+\beta ({\boldsymbol {p}}_{i-1}-\omega _{i-1}{\boldsymbol {v}}_{i-1})}$
4. ${\displaystyle {\boldsymbol {y}}={\boldsymbol {K}}^{-1}{\boldsymbol {p}}_{i}}$
5. ${\displaystyle {\boldsymbol {v}}_{i}={\boldsymbol {Ay}}}$
6. ${\displaystyle \alpha =\rho _{i}/({\boldsymbol {\hat {r}}}_{0},{\boldsymbol {v}}_{i})\;}$
7. ${\displaystyle {\boldsymbol {s}}={\boldsymbol {r}}_{i}-\alpha {\boldsymbol {v}}_{i}}$
8. ${\displaystyle {\boldsymbol {z}}={\boldsymbol {As}}}$
9. ${\displaystyle {\boldsymbol {t}}={\boldsymbol {K}}^{-1}{\boldsymbol {z}}}$
10. ${\displaystyle \omega _{i}=({\boldsymbol {K}}_{1}^{-1}{\boldsymbol {t}},{\boldsymbol {K}}_{1}^{-1}{\boldsymbol {s}})/({\boldsymbol {K}}_{1}^{-1}{\boldsymbol {s}},{\boldsymbol {K}}_{1}^{-1}{\boldsymbol {s}})}$
11. ${\displaystyle {\boldsymbol {x}}_{i}={\boldsymbol {x}}_{i-1}+\alpha {\boldsymbol {y}}+\omega _{i}{\boldsymbol {z}}}$
12. ${\displaystyle {\boldsymbol {x}}_{i}}$ 足够精确则退出
13. ${\displaystyle {\boldsymbol {r}}_{i}={\boldsymbol {s}}-\omega _{i}{\boldsymbol {t}}}$

${\displaystyle {\boldsymbol {{\tilde {A}}{\tilde {x}}}}={\boldsymbol {\tilde {b}}}}$

推导

双共轭梯度法的多项式形式

${\displaystyle {\boldsymbol {p}}_{i}={\boldsymbol {r}}_{i}+\beta _{i}{\boldsymbol {p}}_{i-1}{\text{,}}}$
${\displaystyle {\boldsymbol {\hat {p}}}_{i}={\boldsymbol {\hat {r}}}_{i}+\beta _{i}{\boldsymbol {p}}_{i-1}{\text{,}}}$
${\displaystyle {\boldsymbol {r}}_{i}={\boldsymbol {r}}_{i-1}-\alpha _{i}{\boldsymbol {Ap}}_{i}{\text{,}}}$
${\displaystyle {\boldsymbol {\hat {r}}}_{i}={\boldsymbol {\hat {r}}}_{i-1}-\alpha {\boldsymbol {A}}^{\mathrm {T} }{\boldsymbol {\hat {p}}}_{i}{\text{.}}}$

${\displaystyle \alpha _{i}=\rho _{i}/({\boldsymbol {\hat {p}}}_{i},{\boldsymbol {Ap}}_{i}){\text{,}}}$
${\displaystyle \beta _{i}=\rho _{i}/\rho _{i-1}{\text{,}}\;}$

${\displaystyle ({\boldsymbol {\hat {r}}}_{i},{\boldsymbol {r}}_{j})=0{\text{,}}}$
${\displaystyle ({\boldsymbol {\hat {p}}}_{i},{\boldsymbol {Ap}}_{j})=0{\text{.}}}$

${\displaystyle {\boldsymbol {r}}_{i}=P_{i}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\text{,}}}$
${\displaystyle {\boldsymbol {\hat {r}}}_{i}=P_{i}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0}{\text{,}}}$
${\displaystyle {\boldsymbol {p}}_{i+1}=T_{i}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\text{,}}}$
${\displaystyle {\boldsymbol {\hat {p}}}_{i+1}=T_{i}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0}{\text{,}}}$

${\displaystyle P_{i}({\boldsymbol {A}})=P_{i-1}({\boldsymbol {A}})-\alpha _{i}{\boldsymbol {A}}T_{i-1}({\boldsymbol {A}}){\text{,}}}$
${\displaystyle T_{i}({\boldsymbol {A}})=P_{i}({\boldsymbol {A}})-\beta _{i+1}T_{i-1}({\boldsymbol {A}}){\text{.}}}$

从双共轭梯度法导出稳定双共轭梯度 法

${\displaystyle {\boldsymbol {\tilde {r}}}_{i}=Q_{i}({\boldsymbol {A}})P_{i}({\boldsymbol {A}}){\boldsymbol {r}}_{0}}$

${\displaystyle Q_{i}({\boldsymbol {A}})P_{i}({\boldsymbol {A}}){\boldsymbol {r}}_{0}=({\boldsymbol {I}}-\omega _{i}{\boldsymbol {A}}){\bigl (}Q_{i-1}({\boldsymbol {A}})P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}-\alpha _{i}{\boldsymbol {A}}Q_{i-1}({\boldsymbol {A}})P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}{\text{.}}}$

${\displaystyle Q_{i}({\boldsymbol {A}})T_{i}({\boldsymbol {A}}){\boldsymbol {r}}_{0}=Q_{i}({\boldsymbol {A}})P_{i}({\boldsymbol {A}}){\boldsymbol {r}}_{0}-\beta _{i+1}({\boldsymbol {I}}-\omega _{i}{\boldsymbol {A}})Q_{i-1}({\boldsymbol {A}})T_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\text{.}}}$

${\displaystyle {\boldsymbol {\tilde {p}}}_{i+1}=Q_{i}({\boldsymbol {A}})T_{i}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\text{.}}}$

${\displaystyle {\boldsymbol {\tilde {p}}}_{i}={\boldsymbol {\tilde {r}}}_{i-1}+\beta _{i}({\boldsymbol {I}}-\omega _{i-1}{\boldsymbol {A}}){\boldsymbol {\tilde {p}}}_{i-1}{\text{,}}}$
${\displaystyle {\boldsymbol {\tilde {r}}}_{i}=({\boldsymbol {I}}-\omega _{i}{\boldsymbol {A}})({\boldsymbol {\tilde {r}}}_{i-1}-\alpha _{i}{\boldsymbol {A{\tilde {p}}}}_{i}){\text{.}}}$

${\displaystyle {\boldsymbol {s}}_{i}={\boldsymbol {\tilde {r}}}_{i-1}-\alpha _{i}{\boldsymbol {A{\tilde {p}}}}_{i}{\text{.}}}$

${\displaystyle {\boldsymbol {\tilde {r}}}_{i}={\boldsymbol {\tilde {r}}}_{i-1}-\alpha _{i}{\boldsymbol {A{\tilde {p}}}}_{i}-\omega _{i}{\boldsymbol {As}}_{i}{\text{,}}}$

${\displaystyle {\boldsymbol {x}}_{i}={\boldsymbol {x}}_{i-1}+\alpha _{i}{\boldsymbol {\tilde {p}}}_{i}+\omega _{i}{\boldsymbol {s}}_{i}{\text{.}}}$

确定稳定双共轭梯度法的常数

${\displaystyle \rho _{i}=({\boldsymbol {\hat {r}}}_{i-1},{\boldsymbol {r}}_{i-1})={\bigl (}P_{i-1}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0},P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}{\text{.}}}$

${\displaystyle {\tilde {\rho }}_{i}={\bigl (}Q_{i-1}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0},P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}={\bigl (}{\boldsymbol {\hat {r}}}_{0},Q_{i-1}({\boldsymbol {A}})P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}=({\boldsymbol {\hat {r}}}_{0},{\boldsymbol {r}}_{i-1})}$

${\displaystyle \rho _{i}=(\alpha _{1}/\omega _{1})(\alpha _{2}/\omega _{2})\cdots (\alpha _{i-1}/\omega _{i-1}){\tilde {\rho }}_{i}{\text{,}}}$

${\displaystyle \beta _{i}=\rho _{i}/\rho _{i-1}=({\tilde {\rho }}_{i}/{\tilde {\rho }}_{i-1})(\alpha _{i-1}/\omega _{i-1}){\text{.}}}$

${\displaystyle \alpha _{i}=\rho _{i}/({\boldsymbol {\hat {p}}},{\boldsymbol {Ap}}_{i})={\bigl (}P_{i-1}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0},P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}{\big /}{\bigl (}T_{i-1}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0},{\boldsymbol {A}}T_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}{\text{.}}}$

${\displaystyle \alpha _{i}={\bigl (}Q_{i-1}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0},P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}{\big /}{\bigl (}Q_{i-1}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0},{\boldsymbol {A}}T_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}={\tilde {\rho }}_{i}{\big /}{\bigl (}{\boldsymbol {\hat {r}}}_{0},{\boldsymbol {A}}Q_{i-1}({\boldsymbol {A}})T_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}={\tilde {\rho }}_{i}/({\boldsymbol {\hat {r}}}_{0},{\boldsymbol {A{\tilde {p}}}}_{i}){\text{.}}}$

${\displaystyle {\bigl (}({\boldsymbol {I}}-\omega _{i}{\boldsymbol {A}}){\boldsymbol {s}}_{i},{\boldsymbol {As}}_{i}{\bigr )}=0}$

${\displaystyle \omega _{i}=({\boldsymbol {As}}_{i},{\boldsymbol {s}}_{i})/({\boldsymbol {As}}_{i},{\boldsymbol {As}}_{i}){\text{.}}}$

参考文献

• Van der Vorst, H. A. Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing. 1992, 13: 631–644. doi:10.1137/0913035.
• Saad, Y. §7.4.2 BICGSTAB. Iterative Methods for Sparse Linear Systems 2nd. SIAM. 2003: 231–234. doi:10.2277/0898715342. ISBN 978-0-898715-34-7.