查看“︁稳定双共轭梯度法”︁的源代码

在[[数值线性代数]]中，'''稳定双共轭梯度法'''（{{lang-en|'''Biconjugate gradient stabilized method'''}}，通常简称为{{abbr|'''BiCGSTAB'''}}）是一种由荷兰数学家 H. A. van der Vorst 提出的用于数值求解非对称[[线性方程组]]的[[迭代方法]]。它是[[双共轭梯度法]]（BiCG）的一个变种，比双共轭梯度法本身以及诸如[[共轭梯度平方法]]（CGS）等其他变种有更快速和更平滑的收敛性。它是一种 [[Krylov 子空间]]方法。

==算法步骤==
===无预处理稳定双共轭梯度法===
要求解线性方程组 <math>\boldsymbol{Ax}=\boldsymbol{b}</math>，稳定双共轭梯度法从初始解 <math>\boldsymbol{x}_0</math> 开始按以下步骤迭代：

# <math>\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}</math>
# 任意选择向量 <math>\boldsymbol{\hat{r}}_0\;</math> 使得 <math>(\boldsymbol{\hat{r}}_0,\boldsymbol{r}_0) \neq0\;</math>，例如，<math>\boldsymbol{\hat{r}}_0=\boldsymbol{r}_0\;</math>
# <math>\rho_0=\alpha=\omega_0=1\;</math>
# <math>\boldsymbol{v}_0=\boldsymbol{p}_0=\boldsymbol{0}</math>
# 对 <math>i=1,2,3,\ldots</math>
## <math>\rho_i=(\boldsymbol{\hat{r}}_0,\boldsymbol{r}_{i-1})\;</math>
## <math>\beta=(\rho_i/\rho_{i-1})(\alpha/\omega_{i-1})\;</math>
## <math>\boldsymbol{p}_i=\boldsymbol{r}_{i-1}+\beta(\boldsymbol{p}_{i-1}-\omega_{i-1}\boldsymbol{v}_{i-1})</math>
## <math>\boldsymbol{v}_i=\boldsymbol{Ap}_i</math>
## <math>\alpha=\rho_i/(\boldsymbol{\hat{r}}_0,\boldsymbol{v}_i)\;</math>
## <math>\boldsymbol{s}=\boldsymbol{r}_{i-1}-\alpha\boldsymbol{v}_i</math>
## <math>\boldsymbol{t}=\boldsymbol{As}</math>
## <math>\omega_i=(\boldsymbol{t},\boldsymbol{s})/(\boldsymbol{t},\boldsymbol{t})</math>
## <math>\boldsymbol{x}_i=\boldsymbol{x}_{i-1}+\alpha\boldsymbol{p}_i+\omega_i\boldsymbol{s}</math>
## 若 <math>\boldsymbol{x}_i</math> 足够精确则退出
## <math>\boldsymbol{r}_i=\boldsymbol{s}-\omega_i\boldsymbol{t}</math>

===预处理稳定双共轭梯度法===
[[预处理]]通常被用来加速迭代方法的收敛。要使用预处理子 <math>\boldsymbol{K}=\boldsymbol{K}_1\boldsymbol{K}_2\approx\boldsymbol{A}</math> 来求解线性方程组 <math>\boldsymbol{Ax}=\boldsymbol{b}</math>，预处理稳定双共轭梯度法从初始解 <math>\boldsymbol{x}_0</math> 开始按以下步骤迭代：

# <math>\boldsymbol{r}_0=\boldsymbol{b}-\boldsymbol{Ax}</math>
# 任意选择向量 <math>\boldsymbol{\hat{r}}_0\;</math> 使得 <math>(\boldsymbol{\hat{r}}_0,\boldsymbol{r}_0) \neq0\;</math>，例如，<math>\boldsymbol{\hat{r}}_0=\boldsymbol{r}_0\;</math>
# <math>\rho_0=\alpha=\omega_0=1\;</math>
# <math>\boldsymbol{v}_0=\boldsymbol{p}_0=\boldsymbol{0}</math>
# 对 <math>i=1,2,3,\ldots</math>
##<math>\rho_i=(\boldsymbol{\hat{r}}_0,\boldsymbol{r}_{i-1})\;</math>
## <math>\beta=(\rho_i/\rho_{i-1})(\alpha/\omega_{i-1})\;</math>
## <math>\boldsymbol{p}_i=\boldsymbol{r}_{i-1}+\beta(\boldsymbol{p}_{i-1}-\omega_{i-1}\boldsymbol{v}_{i-1})</math>
## <math>\boldsymbol{y}=\boldsymbol{K}^{-1}\boldsymbol{p}_i</math>
## <math>\boldsymbol{v}_i=\boldsymbol{Ay}</math>
## <math>\alpha=\rho_i/(\boldsymbol{\hat{r}}_0,\boldsymbol{v}_i)\;</math>
## <math>\boldsymbol{s}=\boldsymbol{r}_i-\alpha\boldsymbol{v}_i</math>
## <math>\boldsymbol{z}=\boldsymbol{As}</math>
## <math>\boldsymbol{t}=\boldsymbol{K}^{-1}\boldsymbol{z}</math>
##<math>\omega_i=(\boldsymbol{K}_1^{-1}\boldsymbol{t},\boldsymbol{K}_1^{-1}\boldsymbol{s})/(\boldsymbol{K}_1^{-1}\boldsymbol{t},\boldsymbol{K}_1^{-1}\boldsymbol{t})</math>
## <math>\boldsymbol{x}_i=\boldsymbol{x}_{i-1}+\alpha\boldsymbol{y}+\omega_i\boldsymbol{z}</math>
## 若 <math>\boldsymbol{x}_i</math> 足够精确则退出
## <math>\boldsymbol{r}_i=\boldsymbol{s}-\omega_i\boldsymbol{t}</math>

这个形式等价于将无预处理的稳定双共轭梯度法应用于显式预处理后的方程组
:<math>\boldsymbol{\tilde{A}\tilde{x}}=\boldsymbol{\tilde{b}}</math>，
其中 <math>\boldsymbol{\tilde{A}}=\boldsymbol{K}_1^{-1}\boldsymbol{AK}_2^{-1}</math>，<math>\boldsymbol{\tilde{x}}=\boldsymbol{K}_2\boldsymbol{x}</math>，<math>\boldsymbol{\tilde{b}}=\boldsymbol{K}_1^{-1}\boldsymbol{b}</math>。换句话说，左预处理和右预处理都可以通过这个形式实施。

==推导==
===双共轭梯度法的多项式形式===
在双共轭梯度法中，搜索方向 <math>\boldsymbol{p}_i</math>  和 <math>\boldsymbol{\hat{p}}_i</math> 以及残量 <math>\boldsymbol{r}_i</math> 和 <math>\boldsymbol{\hat{r}}_i</math> 通过以下递推关系更新：

:<math>\boldsymbol{p}_i=\boldsymbol{r}_i+\beta_i\boldsymbol{p}_{i-1}\text{,}</math>
:<math>\boldsymbol{\hat{p}}_i=\boldsymbol{\hat{r}}_i+\beta_i\boldsymbol{p}_{i-1}\text{,}</math>
:<math>\boldsymbol{r}_i=\boldsymbol{r}_{i-1}-\alpha_i\boldsymbol{Ap}_{i}\text{,}</math>
:<math>\boldsymbol{\hat{r}}_i=\boldsymbol{\hat{r}}_{i-1}-\alpha\boldsymbol{A}^{\mathrm{T}}\boldsymbol{\hat{p}}_i\text{.}</math>

常数 <math>\alpha_i\;</math> 和 <math>\beta_i\;</math> 取值为

:<math>\alpha_i=\rho_i/(\boldsymbol{\hat{p}}_i,\boldsymbol{Ap}_i)\text{,}</math>
:<math>\beta_i=\rho_i/\rho_{i-1}\text{,}\;</math>

其中 <math>\rho_i=(\boldsymbol{\hat{r}}_{i-1},\boldsymbol{r}_{i-1})</math>，使得残量和搜索方向分别满足双正交性和双共轭性，也就是对于 <math>i\neq  j</math>，

:<math>(\boldsymbol{\hat{r}}_i,\boldsymbol{r}_j)=0\text{,}</math>
:<math>(\boldsymbol{\hat{p}}_i,\boldsymbol{Ap}_j)=0\text{.}</math>

不难证明，

:<math>\boldsymbol{r}_i=P_i(\boldsymbol{A})\boldsymbol{r}_0\text{,}</math>
:<math>\boldsymbol{\hat{r}}_i=P_i(\boldsymbol{A}^{\mathrm{T}})\boldsymbol{\hat{r}}_0\text{,}</math>
:<math>\boldsymbol{p}_{i+1}=T_i(\boldsymbol{A})\boldsymbol{r}_0\text{,}</math>
:<math>\boldsymbol{\hat{p}}_{i+1}=T_i(\boldsymbol{A}^{\mathrm{T}})\boldsymbol{\hat{r}}_0\text{,}</math>

其中 <math>P_i(\boldsymbol{A})</math> 和 <math>T_i(\boldsymbol{A})</math> 是关于 <math>\boldsymbol{A}</math> 的 <math>i\;</math> 次多项式。这些多项式满足以下递推关系：

:<math>P_i(\boldsymbol{A})=P_{i-1}(\boldsymbol{A})-\alpha_i\boldsymbol{A}T_{i-1}(\boldsymbol{A})\text{,}</math>
:<math>T_i(\boldsymbol{A})=P_i(\boldsymbol{A})-\beta_{i+1}T_{i-1}(\boldsymbol{A})\text{.}</math>

===从双共轭梯度法导出稳定双共轭梯度 法===
双共轭梯度法的残量和搜索方向不是必须显式跟踪的。换句话说，双共轭梯度法的迭代是可以隐式进行的。稳定双共轭梯度法中希望得到

:<math>\boldsymbol{\tilde{r}}_i=Q_i(\boldsymbol{A})P_i(\boldsymbol{A})\boldsymbol{r}_0</math>

的递推关系，其中 <math>Q_i(\boldsymbol{A})=(\boldsymbol{I}-\omega_1\boldsymbol{A})(\boldsymbol{I}-\omega_1\boldsymbol{A})\cdots(\boldsymbol{I}-\omega_i\boldsymbol{A})</math>，<math>\omega_j\;</math> 为适当选取的常数。以此代替 <math>\boldsymbol{r}_i=P_i(\boldsymbol{A})</math> 的目的是希望 <math>Q_i(\boldsymbol{A})</math>  可以使 <math>\boldsymbol{\tilde{r}}_i</math> 有比 <math>\boldsymbol{r}_i</math> 更快速和更平滑的收敛性。

根据 <math>P_i(\boldsymbol{A})</math> 和 <math>T_i(\boldsymbol{A})</math> 的递推关系以及 <math>Q_i(\boldsymbol{A})</math> 的定义，

:<math>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{.}</math>

于是还需要一条关于 <math>Q_i(\boldsymbol{A})T_i(\boldsymbol{A})\boldsymbol{r}_0</math>  的递推关系。这同样可以从双共轭梯度法的递推关系中导出：

:<math>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{.}</math>

类似于 <math>\boldsymbol{\tilde{r}}_i</math>，稳定双共轭梯度法定义

:<math>\boldsymbol{\tilde{p}}_{i+1}=Q_i(\boldsymbol{A})T_i(\boldsymbol{A})\boldsymbol{r}_0\text{.}</math>

写成向量形式，<math>\boldsymbol{\tilde{p}}_i</math> 和 <math>\boldsymbol{\tilde{r}}_i</math> 的递推关系就是

:<math>\boldsymbol{\tilde{p}}_i=\boldsymbol{\tilde{r}}_{i-1}+\beta_i(\boldsymbol{I}-\omega_{i-1}\boldsymbol{A})\boldsymbol{\tilde{p}}_{i-1}\text{,}</math>
:<math>\boldsymbol{\tilde{r}}_i=(\boldsymbol{I}-\omega_i\boldsymbol{A})(\boldsymbol{\tilde{r}}_{i-1}-\alpha_i\boldsymbol{A\tilde{p}}_i)\text{.}</math>

为了导出 <math>\boldsymbol{x}_i</math> 的递推关系，定义

:<math>\boldsymbol{s}_i=\boldsymbol{\tilde{r}}_{i-1}-\alpha_i\boldsymbol{A\tilde{p}}_i\text{.}</math>

于是 <math>\boldsymbol{\tilde{r}}_i</math> 的递推关系就可以写成

:<math>\boldsymbol{\tilde{r}}_i=\boldsymbol{\tilde{r}}_{i-1}-\alpha_i\boldsymbol{A\tilde{p}}_i-\omega_i\boldsymbol{As}_i\text{,}</math>

这对应于

:<math>\boldsymbol{x}_i=\boldsymbol{x}_{i-1}+\alpha_i\boldsymbol{\tilde{p}}_i+\omega_i\boldsymbol{s}_i\text{.}</math>

===确定稳定双共轭梯度法的常数===
现在只需确定双共轭梯度法的常数 <math>\alpha_i\;</math> 和 <math>\beta_i\;</math> 以及选择一个合适的 <math>\omega_i\;</math>。

在双共轭梯度法中，<math>\beta_i=\rho_i/\rho_{i-1}\;</math>， 其中

:<math>\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{.}</math>

由于稳定双共轭梯度法不显式跟踪 <math>\boldsymbol{\hat{r}}_i</math> 或 <math>\boldsymbol{r}_i</math>，<math>\rho_i\;</math> 不能立即用这条公式计算出来。但是，它可以和标量

:<math>\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})</math>

关联起来。由于双正交性，<math>\boldsymbol{r}_{i-1}=P_{i-1}(\boldsymbol{A})\boldsymbol{r}_0</math>  正交于 <math>U_{i-2}(\boldsymbol{A}^{\mathrm{T}})\boldsymbol{\hat{r}}_0</math>，其中 <math>U_{i-2}(\boldsymbol{A}^{\mathrm{T}})</math> 是关于 <math>\boldsymbol{A}^{\mathrm{T}}</math> 的任意 <math>i-2\;</math> 次多项式。因此在点积 <math>\bigl(P_{i-1}(\boldsymbol{A}^{\mathrm{T}})\boldsymbol{\hat{r}}_0,P_{i-1}(\boldsymbol{A})\boldsymbol{r}_0\bigr)</math> 和 <math>\bigl(Q_{i-1}(\boldsymbol{A}^{\mathrm{T}})\boldsymbol{\hat{r}}_0,P_{i-1}(\boldsymbol{A})\boldsymbol{r}_0\bigr)</math> 中只需考虑 <math>P_{i-1}(\boldsymbol{A}^{\mathrm{T}})</math> 和 <math>Q_{i-1}(\boldsymbol{A}^{\mathrm{T}})</math> 的最高次项。<math>P_{i-1}(\boldsymbol{A}^{\mathrm{T}})</math> 和 <math>Q_{i-1}(\boldsymbol{A}^{\mathrm{T}})</math> 的最高次项系数分别是 <math>(-1)^{i-1}\alpha_1\alpha_2\cdots\alpha_{i-1}</math>  和 <math>(-1)^{i-1}\omega_1\omega_2\cdots\omega_{i-1}</math>。因此

:<math>\rho_i=(\alpha_1/\omega_1)(\alpha_2/\omega_2)\cdots(\alpha_{i-1}/\omega_{i-1})\tilde{\rho}_i\text{,}</math>

于是

:<math>\beta_i=\rho_i/\rho_{i-1}=(\tilde{\rho}_i/\tilde{\rho}_{i-1})(\alpha_{i-1}/\omega_{i-1})\text{.}</math>

关于 <math>\alpha_i\;</math> 的简单公式可以类似地导出。在双共轭梯度法中，

:<math>\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{.}</math>

类似于上面的情况，由于双正交性和双共轭性，在点积中只需考虑 <math>P_{i-1}(\boldsymbol{A}^{\mathrm{T}})</math> 和 <math>T_{i-1}(\boldsymbol{A}^{\mathrm{T}})</math> 的最高次项。<math>P_{i-1}(\boldsymbol{A}^{\mathrm{T}})</math> 和 <math>T_{i-1}(\boldsymbol{A}^{\mathrm{T}})</math> 的最高次项系数恰巧是相同的。因此，它们可以在公式中被同时替换为 <math>Q_{i-1}(\boldsymbol{A}^{\mathrm{T}})</math>，于是

:<math>\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{.}</math>

最后，稳定双共轭梯度法选择 <math>\omega_i\;</math> 使得 <math>\boldsymbol{\tilde{r}}_i=(\boldsymbol{I}-\omega_i\boldsymbol{A})\boldsymbol{s}_i</math> 的 2-范数作为 <math>\omega_i\;</math> 的函数被最小化。这在

:<math>\bigl((\boldsymbol{I}-\omega_i\boldsymbol{A})\boldsymbol{s}_i,\boldsymbol{As}_i\bigr)=0</math>

时达到，因此 <math>\omega_i\;</math>  的最优值是

:<math>\omega_i=(\boldsymbol{As}_i,\boldsymbol{s}_i)/(\boldsymbol{As}_i,\boldsymbol{As}_i)\text{.}</math>

==相关主题==
* [[双共轭梯度法]]
* [[共轭梯度法]]

==参考文献==
* {{Cite journal en
| doi = 10.1137/0913035
| title = Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems
| year = 1992
| last = Van der Vorst | first = H. A.
| journal = SIAM Journal on Scientific and Statistical Computing
| volume = 13
| pages = 631–644
| author1 =
}}
* {{cite book en
| last = Saad
| first = Y.
| title =  Iterative Methods for Sparse Linear Systems
| url = https://archive.org/details/iterativemethods00saad_376
| edition = 2nd
| chapter =  §7.4.2 BICGSTAB
| pages = [https://archive.org/details/iterativemethods00saad_376/page/n250 231]–234
| year = 2003
| publisher =  SIAM
| isbn  = 978-0-898715-34-7
| doi =  10.2277/0898715342
}}

[[Category:数值线性代数]]
[[Category:带有伪代码示例的条目]]