稳定双共轭梯度法

要求解线性方程组 ${\displaystyle {\boldsymbol {Ax}}={\boldsymbol {b}}}$，稳定双共轭梯度法从初始解 ${\displaystyle {\boldsymbol {x}}_{0}}$ 开始按以下步骤迭代：

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}}}$

预处理通常被用来加速迭代方法的收敛。要使用预处理子 ${\displaystyle {\boldsymbol {K}}={\boldsymbol {K}}_{1}{\boldsymbol {K}}_{2}\approx {\boldsymbol {A}}}$ 来求解线性方程组 ${\displaystyle {\boldsymbol {Ax}}={\boldsymbol {b}}}$，预处理稳定双共轭梯度法从初始解 ${\displaystyle {\boldsymbol {x}}_{0}}$ 开始按以下步骤迭代：

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 {t}},{\boldsymbol {K}}_{1}^{-1}{\boldsymbol {t}})}$
   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 {\tilde {A}}}={\boldsymbol {K}}_{1}^{-1}{\boldsymbol {AK}}_{2}^{-1}}$，${\displaystyle {\boldsymbol {\tilde {x}}}={\boldsymbol {K}}_{2}{\boldsymbol {x}}}$，${\displaystyle {\boldsymbol {\tilde {b}}}={\boldsymbol {K}}_{1}^{-1}{\boldsymbol {b}}}$。换句话说，左预处理和右预处理都可以通过这个形式实施。

在双共轭梯度法中，搜索方向 ${\displaystyle {\boldsymbol {p}}_{i}}$ 和 ${\displaystyle {\boldsymbol {\hat {p}}}_{i}}$ 以及残量 ${\displaystyle {\boldsymbol {r}}_{i}}$ 和 ${\displaystyle {\boldsymbol {\hat {r}}}_{i}}$ 通过以下递推关系更新：

${\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}\;}$ 和 ${\displaystyle \beta _{i}\;}$ 取值为

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

${\displaystyle \beta _{i}=\rho _{i}/\rho _{i-1}{\text{,}}\;}$

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

${\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}})}$ 和 ${\displaystyle T_{i}({\boldsymbol {A}})}$ 是关于 ${\displaystyle {\boldsymbol {A}}}$ 的 ${\displaystyle i\;}$ 次多项式。这些多项式满足以下递推关系：

${\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}})=({\boldsymbol {I}}-\omega _{1}{\boldsymbol {A}})({\boldsymbol {I}}-\omega _{1}{\boldsymbol {A}})\cdots ({\boldsymbol {I}}-\omega _{i}{\boldsymbol {A}})}$，${\displaystyle \omega _{j}\;}$ 为适当选取的常数。以此代替 ${\displaystyle {\boldsymbol {r}}_{i}=P_{i}({\boldsymbol {A}})}$ 的目的是希望 ${\displaystyle Q_{i}({\boldsymbol {A}})}$ 可以使 ${\displaystyle {\boldsymbol {\tilde {r}}}_{i}}$ 有比 ${\displaystyle {\boldsymbol {r}}_{i}}$ 更快速和更平滑的收敛性。

根据 ${\displaystyle P_{i}({\boldsymbol {A}})}$ 和 ${\displaystyle T_{i}({\boldsymbol {A}})}$ 的递推关系以及 ${\displaystyle Q_{i}({\boldsymbol {A}})}$ 的定义，

${\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}}$ 的递推关系。这同样可以从双共轭梯度法的递推关系中导出：

${\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 {r}}}_{i}}$，稳定双共轭梯度法定义

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

写成向量形式，${\displaystyle {\boldsymbol {\tilde {p}}}_{i}}$ 和 ${\displaystyle {\boldsymbol {\tilde {r}}}_{i}}$ 的递推关系就是

${\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 {x}}_{i}}$ 的递推关系，定义

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

于是 ${\displaystyle {\boldsymbol {\tilde {r}}}_{i}}$ 的递推关系就可以写成

${\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 \alpha _{i}\;}$ 和 ${\displaystyle \beta _{i}\;}$ 以及选择一个合适的 ${\displaystyle \omega _{i}\;}$。

在双共轭梯度法中，${\displaystyle \beta _{i}=\rho _{i}/\rho _{i-1}\;}$， 其中

${\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 {\boldsymbol {\hat {r}}}_{i}}$ 或 ${\displaystyle {\boldsymbol {r}}_{i}}$，${\displaystyle \rho _{i}\;}$ 不能立即用这条公式计算出来。但是，它可以和标量

${\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 {\boldsymbol {r}}_{i-1}=P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}}$ 正交于 ${\displaystyle U_{i-2}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0}}$，其中 ${\displaystyle U_{i-2}({\boldsymbol {A}}^{\mathrm {T} })}$ 是关于 ${\displaystyle {\boldsymbol {A}}^{\mathrm {T} }}$ 的任意 ${\displaystyle i-2\;}$ 次多项式。因此在点积 ${\displaystyle {\bigl (}P_{i-1}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0},P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}}$ 和 ${\displaystyle {\bigl (}Q_{i-1}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0},P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}}$ 中只需考虑 ${\displaystyle P_{i-1}({\boldsymbol {A}}^{\mathrm {T} })}$ 和 ${\displaystyle Q_{i-1}({\boldsymbol {A}}^{\mathrm {T} })}$ 的最高次项。${\displaystyle P_{i-1}({\boldsymbol {A}}^{\mathrm {T} })}$ 和 ${\displaystyle Q_{i-1}({\boldsymbol {A}}^{\mathrm {T} })}$ 的最高次项系数分别是 ${\displaystyle (-1)^{i-1}\alpha _{1}\alpha _{2}\cdots \alpha _{i-1}}$ 和 ${\displaystyle (-1)^{i-1}\omega _{1}\omega _{2}\cdots \omega _{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}\;}$ 的简单公式可以类似地导出。在双共轭梯度法中，

${\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 P_{i-1}({\boldsymbol {A}}^{\mathrm {T} })}$ 和 ${\displaystyle T_{i-1}({\boldsymbol {A}}^{\mathrm {T} })}$ 的最高次项。${\displaystyle P_{i-1}({\boldsymbol {A}}^{\mathrm {T} })}$ 和 ${\displaystyle T_{i-1}({\boldsymbol {A}}^{\mathrm {T} })}$ 的最高次项系数恰巧是相同的。因此，它们可以在公式中被同时替换为 ${\displaystyle Q_{i-1}({\boldsymbol {A}}^{\mathrm {T} })}$，于是

${\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 \omega _{i}\;}$ 使得 ${\displaystyle {\boldsymbol {\tilde {r}}}_{i}=({\boldsymbol {I}}-\omega _{i}{\boldsymbol {A}}){\boldsymbol {s}}_{i}}$ 的 2-范数作为 ${\displaystyle \omega _{i}\;}$ 的函数被最小化。这在

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

时达到，因此 ${\displaystyle \omega _{i}\;}$ 的最优值是

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