GMRES

在數學上，**廣義最小殘量方法**(一般簡稱 **GMRES**) 是一個求解線性方程組 數值解的疊代方法。這個方法利用在 Krylov 子空間中有著最小殘量的向量來逼近解。Arnoldi 疊代方法被用來求解這個向量。

GMRES 方法由 Yousef Saad 和 Martin H. Schultz 在 1986 年提出。

## GMRES 方法

需要求解的線性方程組記為

　
: $Ax=b\,$ 。

假設矩陣 _A_ 是 _n_ $\times $ _n_ 階的可逆的。進一步，假設 _b_ 是標準化的，即 ||_b_|| = 1 (在這篇文章中，||・|| 是 Euclidean 範數)。

這個問題的 _m_ 階 Krylov 子空間是

　
: $K_{m}=\operatorname {span} \,\{b,Ab,A^{2}b,\ldots ,A^{m-1}b\}\,$ 。

GMRES 通過使得殘量 _Ax_m − _b_ 最小的向量 _x_m ∈ _K_m 來逼近 _Ax_ = _b_ 的精確解。

但是，向量 _b_, _Ab_, …, _A_m−1_b_ 幾乎是線性相依的。因此，用 Arnoldi 疊代方法求得的這組 _K_m 的標準正交基

　
: $q_{1},q_{2},\ldots ,q_{m}\,$

來取代上面的那組基。所以，向量 _x_m ∈ _K_m 寫成 _x_m = _Q_m_y_m，其中 _y_m ∈ **R**m 且 _Q_m 是由 _q_1, …, _q_m 組成的 _n_ $\times $ _m_ 矩陣。

Arnoldi 過程也產生一個 (_m_+1) $\times $ _m_ 階上 Hessenberg 矩陣 ${\tilde {H}}_{m}$ 滿足

　
: $AQ_{m}=Q_{m+1}{\tilde {H}}_{m}\,$ 。

因為 $Q_{m}$ 是正交的，我們有

　
: $\|Ax_{m}-b\|=\|{\tilde {H}}_{m}y_{m}-\beta e_{1}\|\,$ ，

其中

　
: $e_{1}=(1,0,0,\ldots ,0)\,$

是 **R**m+1 的標準基的第一個向量，並且

　
: $\beta =\|b-Ax_{0}\|\,$ ，

其中 $x_{0}$ 是初始向量 (通常是零向量)。因此，求使得殘量

　
: $r_{m}={\tilde {H}}_{m}y_{m}-\beta e_{1}$ 。

的範數最小的 $x_{m}$ 。這是一個 _m_ 階線性最小平方問題。

這就是 GMRES 方法。在疊代的每一步中：

1. 做一步 Arnoldi 疊代方法；
2. 尋找使得 ||_r_m|| 最小的 $y_{m}$ ；
3. 計算 $x_{m}=Q_{m}y_{m}$ ；
4. 如果殘量不夠小，重複以上過程。

在每一步疊代中，必須計算一次矩陣向量積 _Aq_m。對於一般的 _n_ 階稠密矩陣，這要計算複雜度大約 2_n_2 次浮點運算。但是對於稀疏矩陣，這個計算複雜度能減少到 O (_n_)。進一步，關於矩陣向量積，在 _m_ 次疊代中能進行 O (_m_ _n_) 次浮點運算。

## 收斂性

第 _m_ 次疊代獲得在 Krylov 子空間 _K_m 下的最小殘量。因為每個子空間包含於下一個子空間中，所以殘量單調遞減。在第 _n_ 次疊代後，其中 _n_ 是矩陣 _A_ 的階數，Krylov 空間 _K_n 是完整的 **R**n。因此，GMRES 方法達到精確解。然而，問題在於：在極少的幾次疊代後 (相對於 _n_)，向量 _x_m 幾乎已經是精確解的一個很好的逼近。

但是，在一般情況下這是不會發生的。事實上，Greenbaum，Pták 和 Strakoš的理論說明了對於每一個單調減少的序列 _a_1, …, _a_n−1, _a_n = 0 ，能夠找到一個矩陣 _A_ 對於所有 _m_ 滿足 ||_r_m|| = _a_m ，其中 _r_m 是上面所定義的殘量。特別的，有可能找到一個矩陣，使得前 _n_− 1 次疊代的殘量一直保持為常數，而只在最後一次疊代時達到零。

在實驗中，GMRES 方法經常表現得很好。在特殊的情況下這能夠被證明。如果 _A_ 是正定的，則

　
: $\|r_{m}\|\leq \left(1-{\frac {\lambda _{\mathrm {min} }(A^{T}+A)}{2\lambda _{\mathrm {max} }(A^{T}+A)}}\right)^{m/2}\|r_{0}\|$ ，

其中 $\lambda _{\mathrm {min} }(M)$ 和 $\lambda _{\mathrm {max} }(M)$ 分別為矩陣 $M$ 的最小和最大特徵值。

如果 _A_ 是對稱的並且是正定的，則

　
: $\|r_{m}\|\leq \left({\frac {\kappa _{2}^{2}(A)-1}{\kappa _{2}^{2}(A)}}\right)^{m/2}\|r_{0}\|$ 。

其中 $\kappa _{2}(A)$ 記為 _A_ 在 Euclidean 範數下的條件數。

一般情況下，其中 _A_ 是非正定的，則

　
: $\|r_{m}\|\leq \inf _{p\in P_{m}}\|p_{m}(A)\|\leq \kappa _{2}(V)\inf _{p\in P_{m}}\max _{\lambda \in \sigma (A)}|p(\lambda )|\,$ ，

其中 _P_m 記為次數不超過 _m_ 且 _p_(0) = 1 的多項式的集合，_V_ 是 _A_ 的譜分解中的矩陣，而 σ(_A_) 是 _A_ 的譜。粗略的說，當 _A_ 的特徵值聚集在遠離原點的區域且 _A_ 離正規不太遠時，收斂速度較快。

所有的不等式只界定殘量，而不是實際誤差 (精確解和當前疊代 _x_m 之間的距離)。

## GMRES 方法的拓展 (Restarted GMRES)

同其他疊代方法一樣，為了加快收斂，GMRES 經常結合預處理方法。

疊代的開銷以 O (_m_2) 增長，其中 _m_ 是疊代次數。然而有時候，GMRES 方法在 _k_ 次疊代後重新開始，即 _x_k 又變回初始值。這樣的方法叫做 GMRES (_k_)。

## 與其他解法的比較

對於對稱矩陣，Arnoldi 疊代方法變成 Lanczos 疊代方法。對應的 Krylov 子空間方法叫做 Paige 和 Saunders 的最小殘量方法 (MinRes)。不像非對稱的情況，MinRes 方法由三項循環關係 (three-term recurrence relation) 給出，並且同 GMRES 一樣，使殘量的範數最小。而對於一般矩陣，Krylov 子空間方法不能由短的循環關係 (short recurrence relation) 給出。

另一類方法由非對稱 Lanczos 疊代方法給出，特別的是 BiCG 方法。這個利用了 three-term recurrence relation，但他們沒有達到最小的殘量，因此對於這些方法殘量不會單調遞減。收斂性是不能保證的。

第三類方法由 CGS 和 BiCGSTAB 給出。這些也由 three-term recurrence relation 給出 (因此，非最優)。而且可能過早的終止疊代了而沒有達到收斂的目的。這些方法的想法是合適的選擇疊代序列所產生的多項式。

對於所有矩陣，這三類方法都不是最好的；總有例使得一類方法好於另一類。因而，各種解法應該進行實際的試驗，來決定對於給定的問題哪一種是最優的。

## 求解最小平方問題

GMRES 方法的其中一部分是求解向量 $y_{m}$ 使得

　
: $\|{\tilde {H}}_{m}y_{m}-e_{1}\|\,$

最小。這個可以通過計算 QR 分解來實現：找到一個 (_m_+1) $\times $ (_m_+1) 階正交矩陣 Ωm 和一個 (_m_+1) $\times $ _m_ 上三角矩陣 ${\tilde {R}}_{m}$ 滿足

　
: $\Omega _{m}{\tilde {H}}_{m}={\tilde {R}}_{m}$ 。

三角矩陣的行數比列數多 1，所以它的最後一行由零組成。因此，它能被分解為

　
: ${\tilde {R}}_{m}={\begin{bmatrix}R_{m}\\0\end{bmatrix}}$ ，

其中 $R_{m}$ 是一個 _m_ $\times $ _m_ 階三角 (方) 矩陣。

QR 分解能夠簡單的進行下去 (update)，從一步疊代到下一步疊代。因為每次的 Hessenberg 矩陣只在一行零元素和一列元素上有所不同：

　
: ${\tilde {H}}_{m+1}={\begin{bmatrix}{\tilde {H}}_{m}&h_{m}\\0&h_{m+1,m}\end{bmatrix}}$ ，

其中 _h_m = (_h_1_m_, … _h_mm)T。這意味著，Hessenberg 矩陣左乘上 Ωm 的擴大矩陣 (通過並上零元素和單位元素)，所得到的是類似於三角矩陣的矩陣：

　
: ${\begin{bmatrix}\Omega _{m}&0\\0&1\end{bmatrix}}{\tilde {H}}_{m+1}={\begin{bmatrix}R_{m}&r_{k}\\0&\rho \\0&\sigma \end{bmatrix}}$

這個矩陣可以三角化，如果 σ 為零。為了修正這個矩陣，需要進行 Givens 旋轉

　
: $G_{m}={\begin{bmatrix}I_{m-1}&0&0\\0&c_{m}&s_{m}\\0&-s_{m}&c_{m}\end{bmatrix}}$

其中

　
: $c_{m}={\frac {\rho }{\sqrt {\rho ^{2}+\sigma ^{2}}}}\quad {\mbox{and}}\quad s_{m}={\frac {\sigma }{\sqrt {\rho ^{2}+\sigma ^{2}}}}$ 。

通過這個 Givens 旋轉，我們構造

　
: $\Omega _{m+1}=G_{m}{\begin{bmatrix}\Omega _{m}&0\\0&1\end{bmatrix}}$ 。

事實上，

　
: $\Omega _{m+1}{\tilde {H}}_{m+1}={\begin{bmatrix}R_{m}&r_{m}\\0&r_{mm}\\0&0\end{bmatrix}}\quad {\text{其 中}}\quad r_{mm}={\sqrt {\rho ^{2}+\sigma ^{2}}}$

是一個三角矩陣。

給出了 QR 分解，最小值問題就容易解決了。注意到

　
: $\|{\tilde {H}}_{m}y_{m}-e_{1}\|=\|\Omega _{m}({\tilde {H}}_{m}y_{m}-e_{1})\|=\|{\tilde {R}}_{m}y_{m}-\Omega _{m}e_{1}\|$ 。

記 $\Omega _{m}e_{1}$ 為

　
: ${\tilde {g}}_{m}={\begin{bmatrix}g_{m}\\\gamma _{m}\end{bmatrix}}$

其中 _g_m ∈ **R**m 和 γm ∈ **R**，則

　
: $\|{\tilde {H}}_{m}y_{m}-e_{1}\|=\|{\tilde {R}}_{m}y_{m}-\Omega _{m}e_{1}\|=\left\|{\begin{bmatrix}R_{m}\\0\end{bmatrix}}y-{\begin{bmatrix}g_{m}\\\gamma _{m}\end{bmatrix}}\right\|$ 。

使得這個表達式最小的向量 _y_ 為

　
: $y_{m}=R_{m}^{-1}g_{m}$ 。

再一次，向量 $g_{m}$ 能夠簡單的進行下去 (update)。

## Example code

### Regular GMRES (python3)

## 註記

## 參考

* A. Meister, _Numerik linearer Gleichungssysteme_, 2nd edition, Vieweg 2005, ISBN 978-3-528-13135-7.
* Y. Saad, _Iterative Methods for Sparse Linear Systems_, 2nd edition, Society for Industrial and Applied Mathematics, 2003. ISBN 978-0-89871-534-7.
* Y. Saad and M.H. Schultz, "GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems", _SIAM J. Sci. Stat. Comput._, **7**:856-869, 1986. doi:10.1137/0907058.
* J. Stoer and R. Bulirsch, _Introduction to numerical analysis_, 3rd edition, Springer, New York, 2002. ISBN 978-0-387-95452-3.
* Lloyd N. Trefethen and David Bau, III, _Numerical Linear Algebra_, Society for Industrial and Applied Mathematics, 1997. ISBN 978-0-89871-361-9.
* Dongarra et al. , _Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods_ （頁面存檔備份，存於網際網路檔案館）, 2nd Edition, SIAM, Philadelphia, 1994
* https://github.com/J-N-ch/GMRES\_py\_restart （頁面存檔備份，存於網際網路檔案館）