穩定雙共軛梯度法

在數值線性代數中，穩定雙共軛梯度法（英語：Biconjugate gradient stabilized method，通常簡稱為BiCGSTAB）是一種由荷蘭數學家 H. A. van der Vorst 提出的用於數值求解非對稱線性方程組的迭代方法。它是雙共軛梯度法（BiCG）的一個變種，比雙共軛梯度法本身以及諸如共軛梯度平方法（CGS）等其他變種有更快速和更平滑的收斂性。它是一種 Krylov 子空間方法。

算法步驟

無預處理穩定雙共軛梯度法

要求解線性方程組 ${\boldsymbol {Ax}}={\boldsymbol {b}}$，穩定雙共軛梯度法從初始解 ${\boldsymbol {x}}_{0}$ 開始按以下步驟迭代：

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

預處理穩定雙共軛梯度法

預處理通常被用來加速迭代方法的收斂。要使用預處理子 ${\boldsymbol {K}}={\boldsymbol {K}}_{1}{\boldsymbol {K}}_{2}\approx {\boldsymbol {A}}$ 來求解線性方程組 ${\boldsymbol {Ax}}={\boldsymbol {b}}$，預處理穩定雙共軛梯度法從初始解 ${\boldsymbol {x}}_{0}$ 開始按以下步驟迭代：

1. ${\boldsymbol {r}}_{0}={\boldsymbol {b}}-{\boldsymbol {Ax}}$
2. 任意選擇向量 ${\boldsymbol {\hat {r}}}_{0}\;$ 使得 $({\boldsymbol {\hat {r}}}_{0},{\boldsymbol {r}}_{0})\neq 0\;$，例如，${\boldsymbol {\hat {r}}}_{0}={\boldsymbol {r}}_{0}\;$
3. $\rho _{0}=\alpha =\omega _{0}=1\;$
4. ${\boldsymbol {v}}_{0}={\boldsymbol {p}}_{0}={\boldsymbol {0}}$
5. 對 $i=1,2,3,\ldots$
   1. $\rho _{i}=({\boldsymbol {\hat {r}}}_{0},{\boldsymbol {r}}_{i-1})\;$
   2. $\beta =(\rho _{i}/\rho _{i-1})(\alpha /\omega _{i-1})\;$
   3. ${\boldsymbol {p}}_{i}={\boldsymbol {r}}_{i-1}+\beta ({\boldsymbol {p}}_{i-1}-\omega _{i-1}{\boldsymbol {v}}_{i-1})$
   4. ${\boldsymbol {y}}={\boldsymbol {K}}^{-1}{\boldsymbol {p}}_{i}$
   5. ${\boldsymbol {v}}_{i}={\boldsymbol {Ay}}$
   6. $\alpha =\rho _{i}/({\boldsymbol {\hat {r}}}_{0},{\boldsymbol {v}}_{i})\;$
   7. ${\boldsymbol {s}}={\boldsymbol {r}}_{i}-\alpha {\boldsymbol {v}}_{i}$
   8. ${\boldsymbol {z}}={\boldsymbol {As}}$
   9. ${\boldsymbol {t}}={\boldsymbol {K}}^{-1}{\boldsymbol {z}}$
   10. $\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. ${\boldsymbol {x}}_{i}={\boldsymbol {x}}_{i-1}+\alpha {\boldsymbol {y}}+\omega _{i}{\boldsymbol {z}}$
   12. 若 ${\boldsymbol {x}}_{i}$ 足夠精確則退出
   13. ${\boldsymbol {r}}_{i}={\boldsymbol {s}}-\omega _{i}{\boldsymbol {t}}$

這個形式等價於將無預處理的穩定雙共軛梯度法應用於顯式預處理後的方程組

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

其中 ${\boldsymbol {\tilde {A}}}={\boldsymbol {K}}_{1}^{-1}{\boldsymbol {AK}}_{2}^{-1}$，${\boldsymbol {\tilde {x}}}={\boldsymbol {K}}_{2}{\boldsymbol {x}}$，${\boldsymbol {\tilde {b}}}={\boldsymbol {K}}_{1}^{-1}{\boldsymbol {b}}$。換句話說，左預處理和右預處理都可以通過這個形式實施。

推導

雙共軛梯度法的多項式形式

在雙共軛梯度法中，搜索方向 ${\boldsymbol {p}}_{i}$ 和 ${\boldsymbol {\hat {p}}}_{i}$ 以及殘量 ${\boldsymbol {r}}_{i}$ 和 ${\boldsymbol {\hat {r}}}_{i}$ 通過以下遞推關係更新：

${\boldsymbol {p}}_{i}={\boldsymbol {r}}_{i}+\beta _{i}{\boldsymbol {p}}_{i-1}{\text{,}}$

${\boldsymbol {\hat {p}}}_{i}={\boldsymbol {\hat {r}}}_{i}+\beta _{i}{\boldsymbol {p}}_{i-1}{\text{,}}$

${\boldsymbol {r}}_{i}={\boldsymbol {r}}_{i-1}-\alpha _{i}{\boldsymbol {Ap}}_{i}{\text{,}}$

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

常數 $\alpha _{i}\;$ 和 $\beta _{i}\;$ 取值為

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

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

其中 $\rho _{i}=({\boldsymbol {\hat {r}}}_{i-1},{\boldsymbol {r}}_{i-1})$，使得殘量和搜索方向分別滿足雙正交性和雙共軛性，也就是對於 $i\neq j$，

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

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

不難證明，

${\boldsymbol {r}}_{i}=P_{i}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\text{,}}$

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

${\boldsymbol {p}}_{i+1}=T_{i}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\text{,}}$

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

其中 $P_{i}({\boldsymbol {A}})$ 和 $T_{i}({\boldsymbol {A}})$ 是關於 ${\boldsymbol {A}}$ 的 $i\;$ 次多項式。這些多項式滿足以下遞推關係：

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

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

從雙共軛梯度法導出穩定雙共軛梯度 法

雙共軛梯度法的殘量和搜索方向不是必須顯式跟蹤的。換句話說，雙共軛梯度法的迭代是可以隱式進行的。穩定雙共軛梯度法中希望得到

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

的遞推關係，其中 $Q_{i}({\boldsymbol {A}})=({\boldsymbol {I}}-\omega _{1}{\boldsymbol {A}})({\boldsymbol {I}}-\omega _{1}{\boldsymbol {A}})\cdots ({\boldsymbol {I}}-\omega _{i}{\boldsymbol {A}})$，$\omega _{j}\;$ 為適當選取的常數。以此代替 ${\boldsymbol {r}}_{i}=P_{i}({\boldsymbol {A}})$ 的目的是希望 $Q_{i}({\boldsymbol {A}})$ 可以使 ${\boldsymbol {\tilde {r}}}_{i}$ 有比 ${\boldsymbol {r}}_{i}$ 更快速和更平滑的收斂性。

根據 $P_{i}({\boldsymbol {A}})$ 和 $T_{i}({\boldsymbol {A}})$ 的遞推關係以及 $Q_{i}({\boldsymbol {A}})$ 的定義，

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

於是還需要一條關於 $Q_{i}({\boldsymbol {A}})T_{i}({\boldsymbol {A}}){\boldsymbol {r}}_{0}$ 的遞推關係。這同樣可以從雙共軛梯度法的遞推關係中導出：

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

類似於 ${\boldsymbol {\tilde {r}}}_{i}$，穩定雙共軛梯度法定義

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

寫成向量形式，${\boldsymbol {\tilde {p}}}_{i}$ 和 ${\boldsymbol {\tilde {r}}}_{i}$ 的遞推關係就是

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

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

為了導出 ${\boldsymbol {x}}_{i}$ 的遞推關係，定義

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

於是 ${\boldsymbol {\tilde {r}}}_{i}$ 的遞推關係就可以寫成

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

這對應於

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

確定穩定雙共軛梯度法的常數

現在只需確定雙共軛梯度法的常數 $\alpha _{i}\;$ 和 $\beta _{i}\;$ 以及選擇一個合適的 $\omega _{i}\;$。

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

$\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{.}}$

由於穩定雙共軛梯度法不顯式跟蹤 ${\boldsymbol {\hat {r}}}_{i}$ 或 ${\boldsymbol {r}}_{i}$，$\rho _{i}\;$ 不能立即用這條公式計算出來。但是，它可以和標量

${\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})$

關聯起來。由於雙正交性，${\boldsymbol {r}}_{i-1}=P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}$ 正交於 $U_{i-2}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0}$，其中 $U_{i-2}({\boldsymbol {A}}^{\mathrm {T} })$ 是關於 ${\boldsymbol {A}}^{\mathrm {T} }$ 的任意 $i-2\;$ 次多項式。因此在點積 ${\bigl (}P_{i-1}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0},P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}$ 和 ${\bigl (}Q_{i-1}({\boldsymbol {A}}^{\mathrm {T} }){\boldsymbol {\hat {r}}}_{0},P_{i-1}({\boldsymbol {A}}){\boldsymbol {r}}_{0}{\bigr )}$ 中只需考慮 $P_{i-1}({\boldsymbol {A}}^{\mathrm {T} })$ 和 $Q_{i-1}({\boldsymbol {A}}^{\mathrm {T} })$ 的最高次項。$P_{i-1}({\boldsymbol {A}}^{\mathrm {T} })$ 和 $Q_{i-1}({\boldsymbol {A}}^{\mathrm {T} })$ 的最高次項係數分別是 $(-1)^{i-1}\alpha _{1}\alpha _{2}\cdots \alpha _{i-1}$ 和 $(-1)^{i-1}\omega _{1}\omega _{2}\cdots \omega _{i-1}$。因此

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

於是

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

關於 $\alpha _{i}\;$ 的簡單公式可以類似地導出。在雙共軛梯度法中，

$\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{.}}$

類似於上面的情況，由於雙正交性和雙共軛性，在點積中只需考慮 $P_{i-1}({\boldsymbol {A}}^{\mathrm {T} })$ 和 $T_{i-1}({\boldsymbol {A}}^{\mathrm {T} })$ 的最高次項。$P_{i-1}({\boldsymbol {A}}^{\mathrm {T} })$ 和 $T_{i-1}({\boldsymbol {A}}^{\mathrm {T} })$ 的最高次項係數恰巧是相同的。因此，它們可以在公式中被同時替換為 $Q_{i-1}({\boldsymbol {A}}^{\mathrm {T} })$，於是

$\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{.}}$

最後，穩定雙共軛梯度法選擇 $\omega _{i}\;$ 使得 ${\boldsymbol {\tilde {r}}}_{i}=({\boldsymbol {I}}-\omega _{i}{\boldsymbol {A}}){\boldsymbol {s}}_{i}$ 的 2-範數作為 $\omega _{i}\;$ 的函數被最小化。這在

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

時達到，因此 $\omega _{i}\;$ 的最優值是

$\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. ISBN 978-0-898715-34-7. doi:10.2277/0898715342.