Krylov 子空間法──線性方程的數值解法 (二):GMRES 與 FOM

本文的閱讀等級:高級

A 為一 n\times n 階實矩陣。對於非零向量 \mathbf{v}\in\mathbb{R}^n

\mathcal{K}_j(A,\mathbf{v})=\hbox{span}\{\mathbf{v},A\mathbf{v},\ldots,A^{j-1}\mathbf{v}\}

稱為 Krylov 子空間。設 m(t) 是次數最小的多項式使得 m(A)\mathbf{v}=\mathbf{0},稱為 \mathbf{v} 相對 A 的最小多項式。Krylov 子空間應用於線性方程的數值解法建立於下列基礎 (見“Krylov 子空間法──線性方程的數值解法 (一):Arnoldi 與 Lanczos 算法”):

  1. k=\deg m(t)。若 j\le k,則 Krylov 序列 \{\mathbf{v},A\mathbf{v}_,\ldots,A^{j-1}\mathbf{v}\} 為一線性獨立集,就有 \dim\mathcal{K}_j(A,\mathbf{v})=j。若 j>k,則 \mathcal{K}_{j}(A,\mathbf{v})=\mathcal{K}_k(A,\mathbf{v})
  2. j\le k,Arnoldi 算法可求得 Krylov 子空間 \mathcal{K}_j(A,\mathbf{v}) 的一組單範正交基底 (orthonormal basis) \{\mathbf{q}_1,\ldots,\mathbf{q}_j\}
  3. Arnoldi 算法給出 Arnoldi 關係式:

    AQ_j=Q_jH_j+h_{j+1,j}\mathbf{q}_{j+1}\mathbf{e}_j^T

    其中 Q_j=\begin{bmatrix}  \mathbf{q}_1&\cdots&\mathbf{q}_j  \end{bmatrix}n\times j 階矩陣滿足 Q_j^TQ_j=I_jH_jj\times j 階上 Hessenberg 矩陣 (見“特殊矩陣 (19):Hessenberg 矩陣”),\mathbf{e}_j^T=(0,\ldots,0,1)j 維單位向量。因為 \mathbf{q}_{j+1}\perp\mathbf{q}_i1\le i\le j,上式左乘 Q_j^T,可得 Q_j^TAQ_j=H_j

 
給定一可逆矩陣 A\mathbf{b}\neq\mathbf{0},何以 Krylov 子空間可用於計算線性方程 A\mathbf{x}=\mathbf{b} 的解?令 \mathbf{x}_\ast=A^{-1}\mathbf{b} 為精確解,\mathbf{x}_0 為初始猜測解,\mathbf{z}=\mathbf{x}_\ast-\mathbf{x}_0,且 \mathbf{r}_0=\mathbf{b}-A\mathbf{x}_0,稱為殘差 (residual)。因為 A\mathbf{z}=A\mathbf{x}_\ast-A\mathbf{x}_0=\mathbf{b}-A\mathbf{x}_0=\mathbf{r}_0,線性方程 A\mathbf{x}=\mathbf{b} 等價於 A\mathbf{z}=\mathbf{r}_0\mathbf{r}_0\neq 0 (否則 \mathbf{x}_0 即為精確解)。Cayley-Hamilton 定理保證 A^{-1} 可表示為 A^{n-1},\ldots,A,I 的線性組合 (見“利用 Cayley-Hamilton 定理計算矩陣函數”):

A^{-1}=c_{n-1}A^{n-1}+\cdots+c_1A+c_0I

推論線性方程 A\mathbf{z}=\mathbf{r}_0 的精確解為

\mathbf{z}_\ast=A^{-1}\mathbf{r}_0=c_{n-1}A^{n-1}\mathbf{r}_0+\cdots+c_1A\mathbf{r}_0+c_0\mathbf{r}_0

因此,\mathbf{z}_\ast\in \mathcal{K}_{n}(A,\mathbf{r}_0)=\mathcal{K}_k(A,\mathbf{r}_0),其中 k\mathbf{r}_0 相對於 A 的最小多項式的次數。透過 Arnoldi 算法找到 \mathcal{K}_k(A,\mathbf{r}_0) 的一組單範正交基底可解出 \mathbf{z}_\ast,也就得到 A\mathbf{x}=\mathbf{b} 的精確解 \mathbf{x}_\ast=\mathbf{x}_0+\mathbf{z}_\ast。然而,數值計算存在誤差,最小多項式的次數 k 幾不可得。現實上,我們可以逐次計算近似解 \mathbf{x}_j=\mathbf{x}_0+\mathbf{z}_j,其中 \mathbf{z}_j\in\mathcal{K}_j(A,\mathbf{r}_0) 滿足某個優化條件。這種運用 Krylov 子空間的線性方程的數值解法,統稱 Krylov 子空間法。本文介紹採用 Arnoldi 算法的廣義最小殘量法 (generalized minimal residual method,簡稱 GMRES) 和滿正交法 (full orthogonalization method,簡稱 FOM)。

 
考慮 Krylov 子空間 \mathcal{K}_j=\hbox{span}\{\mathbf{r}_0,A\mathbf{r}_0,\ldots,A^{j-1}\mathbf{r}_0\}\hbox{rank}\mathcal{K}_j=j。在不造成混淆的情況下,以下用 \mathcal{K}_j 表示 \mathcal{K}_j(A,\mathbf{r}_0)。Arnoldi 算法得到 \mathcal{K}_j 的一組單範正交基底 \{\mathbf{q}_1,\ldots,\mathbf{q}_j\},其中 \mathbf{q}_1=\mathbf{r}_0/\Vert\mathbf{r}_0\Vert。Arnoldi 關係式可表示為

AQ_j=Q_{j+1}\tilde{H}_j

其中 \tilde{H}_j=\begin{bmatrix}  H_j\\  h_{j+1,j}\mathbf{e}_j^T  \end{bmatrix}。GMRES 是一個限定於 Krylov 子空間 \mathcal{K}_j 的最小平方法。在 Arnoldi 算法的每一次迭代,GMRES 計算 A\mathbf{z}=\mathbf{r}_0 的最佳近似解,具體地說,尋找 \mathbf{z}\in\mathcal{K}_j 使得殘量 \Vert\mathbf{r}_0-A\mathbf{z}\Vert 最小。因為 \mathcal{K}_j=\hbox{span}\{\mathbf{q}_1,\ldots,\mathbf{q}_j\}=C(Q_j),其中 C(Q_j) 表示 Q_j 的行空間 (column space),可設 \mathbf{z}=Q_j\mathbf{y}\mathbf{y}\in\mathbb{R}^j。使用 Arnoldi 關係式,

\displaystyle  \begin{aligned}  \mathbf{r}_0-A\mathbf{z}&=\mathbf{r}_0-AQ_j\mathbf{y}\\  &=\Vert\mathbf{r}_0\Vert\mathbf{q}_1-Q_{j+1}\tilde{H}_j\mathbf{y}\\  &=Q_{j+1}\left(\Vert\mathbf{r}_0\Vert\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y}\right),  \end{aligned}

其中 \tilde{\mathbf{e}}_1=(1,0,\ldots,0)^Tj+1 維單位向量。因為 Q_{j+1}^TQ_{j+1}=I_{j+1}

\displaystyle  \begin{aligned}  \Vert\mathbf{r}_0-A\mathbf{z}\Vert^2&=\left\|\mathbf{r}_0-AQ_j\mathbf{y}\right\|^2\\  &=\left\| Q_{j+1}\left(\Vert\mathbf{r}_0\Vert\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y}\right)\right\|^2\\  &= \left(\Vert\mathbf{r}_0\Vert\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y}\right)^T Q_{j+1}^T Q_{j+1}\left(\Vert\mathbf{r}_0\Vert\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y}\right)\\  &=\left\| \Vert\mathbf{r}_0\Vert\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y}\right\|^2.  \end{aligned}

因此,線性方程 A\mathbf{z}=\mathbf{r}_0 限定於 \mathbf{z}\in\mathcal{K}_j 的最小殘量近似解為 \mathbf{z}_j=Q_j\mathbf{y}_j,其中 \mathbf{y}_jAQ_j\mathbf{y}=\mathbf{r}_0\tilde{H}_j\mathbf{y}=\Vert\mathbf{r}_0\Vert\tilde{\mathbf{e}}_1 的最小平方解,滿足

\displaystyle\begin{aligned}  \mathbf{y}_j&=\arg \min_{\mathbf{y}}\left\|\mathbf{r}_0-AQ_j\mathbf{y}\right\|\\  &=\arg \min_{\mathbf{y}}\left\|\Vert\mathbf{r}_0\Vert\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y}\right\|.\end{aligned}

線性方程 A\mathbf{x}=\mathbf{b} 的最小殘量近似解則為

\mathbf{x}_j=\mathbf{x}_0+\mathbf{z}_j=\mathbf{x}_0+Q_j\mathbf{y}_j

殘差為

\displaystyle\begin{aligned}  \mathbf{b}-A\mathbf{x}_j&=\mathbf{b}-A(\mathbf{x}_0+Q_j\mathbf{y}_j)\\  &=\mathbf{b}-A\mathbf{x}_0-AQ_j\mathbf{y}_j\\  &=\mathbf{r}_0-AQ_j\mathbf{y}_j\\  &=\Vert\mathbf{r}_0\Vert\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y}_j.  \end{aligned}

 
我們將 GMRES 算法整理於下:對於一 n\times n 階可逆矩陣 A 和非零向量 \mathbf{b}\in\mathbb{R}^n,設定初始解 \mathbf{x}_0 (一般設為零)。

GMRES


  1. \mathbf{r}_0\leftarrow\mathbf{b}-A\mathbf{x}_0\beta\leftarrow\Vert\mathbf{r}_0\Vert\mathbf{q}_1\leftarrow \mathbf{r}_0/\beta
  2. 對於 j=1,2,\ldots,n
  3.     執行 Arnoldi 算法的第 j 次迭代,得到 \tilde{H}_jQ_j
  4.     \mathbf{y}_j\leftarrow\arg\min_{\mathbf{y}}\left\|\beta\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y}\right\|
  5.     \mathbf{x}_j\leftarrow \mathbf{x}_0+Q_j\mathbf{y}_j
  6.     若 \Vert\mathbf{b}-A\mathbf{x}_j\Vert/\Vert\mathbf{b}\Vert<\epsilon,停止計算。

當 Arnoldi 算法停止計算 (表示 j\mathbf{r}_0 相對於 A 的最小多項式次數),則 \mathbf{x}_j 為精確解 (若不考慮數值計算誤差)。因為 \mathcal{K}_j\subset\mathcal{K}_{j+1}\min_{\mathbf{z}\in\mathcal{K}_{j+1}}\Vert\mathbf{r}_0-A\mathbf{z}\Vert\le\min_{\mathbf{z}\in\mathcal{K}_j}\Vert\mathbf{r}_0-A\mathbf{z}\Vert,也就是說 GMRES 具有單調收斂性。為控制收斂,設相對殘量 \Vert\mathbf{b}-A\mathbf{x}_j\Vert/\Vert\mathbf{b}\Vert 的容許誤差為 \epsilon。在許多實際應用上,GMRES 的收斂表現良好 (針對對稱正定矩陣的收斂性分析請見附註[1])。

 
下面說明如何計算 \mathbf{y}_j 使最小化 \left\|\beta\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y}\right\|,其中 \tilde{H}_j(j+1)\times j 階矩陣,\beta=\Vert\mathbf{r}_0\Vert。計算 \tilde{H}_j 的 QR 分解,可得一 (j+1)\times (j+1) 階正交矩陣 (orthogonal matrix) \Omega_j\Omega_j^T=\Omega^{-1},以及 (j+1)\times j 階矩陣 \tilde{R}_j 使得

\Omega_j\tilde{H}_j=\tilde{R}_j

其中 \tilde{R}_j=\begin{bmatrix}  R_j\\  0  \end{bmatrix}R_jj\times j 階上三角矩陣 (見“Givens 旋轉於 QR 分解的應用”)。因為 \tilde{H}_j 具有特殊形態,QR 分解格外簡單。寫出

\displaystyle  \tilde{H}_{j+1}=\begin{bmatrix}  h_{11}&h_{12}& h_{13}&\cdots&h_{1j}&h_{1,j+1}\\  h_{21}&h_{22}&h_{23}&\cdots&h_{2j}&h_{2,j+1}\\  0&h_{32}&h_{33}&\cdots&h_{3j}&h_{3,j+1}\\  \vdots&\ddots&\ddots&\ddots&\vdots&\vdots\\  0&\cdots&0&h_{j,j-1}&h_{jj}&h_{j,j+1}\\  0&\cdots&\cdots&0&h_{j+1,j}&h_{j+1,j+1}\\  0&\cdots&\cdots&\cdots&0&h_{j+2,j+1}  \end{bmatrix}=\begin{bmatrix}  \tilde{H}_j&\mathbf{h}_{j+1}\\  0&h_{j+2,j+1}  \end{bmatrix}

其中 \mathbf{h}_{j+1}=(h_{1,j+1},\ldots,h_{j+1,j+1})^T。考慮分塊矩陣乘法

\begin{aligned}  \begin{bmatrix}  \Omega_j&0\\  0&1  \end{bmatrix}\tilde{H}_{j+1}&=\begin{bmatrix}  \Omega_j&0\\  0&1  \end{bmatrix}\begin{bmatrix}  \tilde{H}_j&\mathbf{h}_{j+1}\\  0&h_{j+2,j+1}  \end{bmatrix}=\begin{bmatrix}  \Omega_j\tilde{H}_j&\Omega_j\mathbf{h}_{j+1}\\  0&h_{j+2,j+1}  \end{bmatrix}\\  &=\begin{bmatrix}  \tilde{R}_j&\Omega_j\mathbf{h}_{j+1}\\  0&h_{j+2,j+1}  \end{bmatrix}=\begin{bmatrix}  R_j&\boldsymbol{\rho}_j\\  0&\rho_{j+1}\\  0&h_{j+2,j+1}  \end{bmatrix},\end{aligned}

其中 \boldsymbol{\rho}_j\Omega_j\mathbf{h}_{j+1} 的前 j 個元組成的向量,\rho_{j+1}\Omega_j\mathbf{h}_{j+1} 的最末一個元。為了產生上三角矩陣,我們希望消滅 h_{j+2,j+1}。設 Givens 旋轉矩陣

G_j=\begin{bmatrix}  I_j&0&0\\  0&c_j&-s_j\\  0&s_j&c_j  \end{bmatrix}

其中 c_j={\rho_{j+1}}\big/{\sqrt{\rho_{j+1}^2+h_{j+2,j+1}^2}}s_j=-{h_{j+2,j+1}}\big/{\sqrt{\rho_{j+1}^2+h_{j+2,j+1}^2}}。令 \Omega_{j+1}=G_j\begin{bmatrix}  \Omega_j&0\\  0&1  \end{bmatrix}。不難驗證 \Omega_{j+1}^T=\Omega_{j+1}^{-1},以及

\begin{aligned}  \Omega_{j+1}\tilde{H}_{j+1}&=\begin{bmatrix}  I_j&0&0\\  0&c_j&-s_j\\  0&s_j&c_j  \end{bmatrix}\begin{bmatrix}  R_j&\boldsymbol{\rho}_{j}\\  0&\rho_{j+1}\\  0&h_{j+2,j+1}  \end{bmatrix}=\begin{bmatrix}  R_j&\boldsymbol{\rho}_{j}\\  0&r_{j+1,j+1}\\  0&0  \end{bmatrix}\\  &=\begin{bmatrix}  R_{j+1}\\  0  \end{bmatrix}=\tilde{R}_{j+1},  \end{aligned}

其中 r_{j+1,j+1}=\sqrt{\rho_{j+1}^2+h_{j+2,j+1}^2}。以上述的的迭代方式計算 \Omega_{j+1}R_{j+1} 僅耗費 \mathcal{O}(j) 個運算。

 
給定 QR 分解 \Omega_j\tilde{H}_j=\tilde{R}_j,殘量可表示為

\begin{aligned}  \left\|\beta\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y}\right\|  =\left\|\Omega_j(\beta\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y})\right\|  =\left\|\Omega_j(\beta\tilde{\mathbf{e}}_1)-\tilde{R}_j\mathbf{y}\right\|\\  =\left\|\begin{bmatrix}  \mathbf{g}_j\\  g_{j+1}  \end{bmatrix}-\begin{bmatrix}  R_j\\  0  \end{bmatrix}\mathbf{y}\right\|  =\left\|\begin{bmatrix}\mathbf{g}_j-R_j\mathbf{y}\\  g_{j+1}  \end{bmatrix}\right\|,  \end{aligned}

其中 \begin{bmatrix}  \mathbf{g}_j\\  g_{j+1}  \end{bmatrix}=\Omega_j(\beta\tilde{\mathbf{e}}_1)。所以,最小平方解滿足 R_j\mathbf{y}_j=\mathbf{g}_j。因為 R_j 是上三角矩陣,反向代入即可解出 \mathbf{y}_j。上式並給出最小殘量:

\left\|\mathbf{b}-A\mathbf{x}_j\right\|=\left\|\beta\tilde{\mathbf{e}}_1-\tilde{H}_j\mathbf{y}_j\right\|=\vert g_{j+1}\vert=\beta\left|(\Omega_j)_{j+1,1}\right|

 
GMRES 的運算量取決於 Arnoldi 算法以及最小平方算法。Arnoldi 算法的第 j 次迭代使用 \mathcal{O}(jn) 個運算,另加上矩陣─向量乘法,運算量為 \mathcal{O}(n^2)。利用 Givens 旋轉迭代計算 QR 分解與最小平方解僅耗用 \mathcal{O}(j) 個運算。為了控制 Arnoldi 算法的計算量並加速收斂,可設迭代次數為 m。給定一初始解 \mathbf{x}_0,待 Arnoldi 算法產生 \tilde{H}_mQ_m 後才計算最小平方解 \mathbf{x}_m,再設新的初始解為 \mathbf{x}_m,重啟計算程序。這個方法稱為 GMRES(m) 或重啟 (restarted) GMRES,如下:

GMRES(m)


  1. \mathbf{r}_0\leftarrow\mathbf{b}-A\mathbf{x}_0\beta\leftarrow \Vert\mathbf{r}_0\Vert\mathbf{q}_1\leftarrow \mathbf{r}_0/\beta
  2. 執行 Arnoldi 算法,得到 \tilde{H}_mQ_m
  3. \mathbf{y}_m\leftarrow\arg\min_{\mathbf{y}}\left\|\Vert\beta\Vert\tilde{\mathbf{e}}_1-\tilde{H}_m\mathbf{y}\right\|
  4. \mathbf{x}_m\leftarrow \mathbf{x}_0+Q_m\mathbf{y}_m
  5. \Vert\mathbf{b}-A\mathbf{x}_m\Vert/\Vert\mathbf{b}\Vert<\epsilon,停止計算。
  6. \mathbf{x}_0\leftarrow\mathbf{x}_m,返回步驟1。

 
接下來我們討論 GMRES 的涵義與推廣。令 A(\mathcal{K}_j) 代表 \mathcal{K}_jA 映射所得的子空間,A(\mathcal{K}_j)=\{A\mathbf{z}\vert\mathbf{z}\in\mathcal{K}_j\}。因此,A(\mathcal{K}_j)=C(AQ_j)\hbox{rank} A(\mathcal{K}_j)=\hbox{rank}\mathcal{K}_j=j。線性方程 AQ_j\mathbf{y}=\mathbf{r}_0 的最小平方解 \mathbf{y}_j 滿足正交性質,\mathbf{r}_0-AQ_j\mathbf{y}_j\perp C(AQ_j) (見“從線性變換解釋最小平方近似”)。GMRES 所設定的近似解滿足 \mathbf{x}_j-\mathbf{x}_0\in\mathcal{K}_j (即 C(Q_j)) 和 \mathbf{b}-A\mathbf{x}_j\perp A(\mathcal{K}_j) (即 C(AQ_j))。推廣這個結果,Krylov 子空間法可用下列兩個條件來定義:

  1. \mathbf{x}_j-\mathbf{x}_0\in \mathcal{K}_j
  2. \mathbf{b}-A\mathbf{x}_j\perp \mathcal{L}_j

其中子空間 \mathcal{K}_j\mathcal{L}_j 的維數皆等於 j,條件2稱為 Petrov-Galerkin 條件。因為這個緣故,Krylov 子空間法被歸類為投影法。令 P_jn\times j 階矩陣使得 C(P_j)=\mathcal{L}_j。矩陣 Q_jP_j 的線性獨立行向量分別組成 \mathcal{K}_j\mathcal{L}_j 的基底。根據條件1,設 \mathbf{x}_j=\mathbf{x}_0+Q_j\mathbf{y}_j。左零空間 N(P_j^T) 是行空間 C(P_j) 的正交補餘 (見“正交補餘與投影定理”),條件2說明 \mathbf{b}-A\mathbf{x}_j\in N(P^T_j)。合併上面兩式,

\displaystyle\begin{aligned}  \mathbf{0}&=P_j^T(\mathbf{b}-A\mathbf{x}_j)\\  &=P_j^T(\mathbf{b}-A\mathbf{x}_0-AQ_j\mathbf{y}_j)\\  &=P_j^T\mathbf{r}_0-P_j^T AQ_j\mathbf{y}_j.  \end{aligned}

Krylov 子空間法的近似解 \mathbf{y}_j 滿足正規方程 (normal equation)

\displaystyle  P_j^T AQ_j\mathbf{y}_j=P_j^T\mathbf{r}_0

A 是一 n\times n 階大尺寸矩陣時,Krylov 子空間法運用上式計算近似解不啻是一個令人讚嘆的作法,因為 P_j^TAQ_j 是一個「接近」Aj\times j 階小矩陣。

 
下面介紹另一種 Krylov 子空間法,稱作 FOM。類似 GMRES,FOM 也採用 Arnoldi 算法,設定的近似解 \mathbf{x}_j 滿足 \mathbf{x}_j-\mathbf{x}_0\in\mathcal{K}_j,但 Petrov-Galerkin 條件為 \mathbf{b}-A\mathbf{x}_j\perp\mathcal{K}_j。FOM 與 GMRES 的差異在於 FOM 定義 \mathcal{L}_j=\mathcal{K}_j,GMRES 定義 \mathcal{L}_j=A(\mathcal{K}_j)。對於一般矩陣 A,FOM 所優化的目標函數並不明確。若 A 是一對稱正定矩陣,我們可證明 FOM 尋找 \mathbf{z}\in\mathcal{K}_j 使最小化

f(\mathbf{z})=(\mathbf{x}_\ast-(\mathbf{x}_0+\mathbf{z}))^T A(\mathbf{x}_\ast-(\mathbf{x}_0+\mathbf{z}))

\mathbf{z}=Q_j\mathbf{y}。上面的約束優化問題等價於求 \mathbf{y}\in\mathbb{R}^j 使最小化二次型:

f(\mathbf{y})=(\mathbf{x}_\ast-(\mathbf{x}_0+Q_j\mathbf{y}))^T A(\mathbf{x}_\ast-(\mathbf{x}_0+Q_j\mathbf{y}))

因為 A^T=A,求導可得 (見“矩陣導數”,SV-10)

\displaystyle  \frac{\partial f}{\partial\mathbf{y}}=-2Q_j^T(A\mathbf{x}_\ast-A(\mathbf{x}_0+Q_j\mathbf{y}))=-2Q_j^T(\mathbf{b}-A(\mathbf{x}_0+Q_j\mathbf{y}))

近似解 \mathbf{x}_j=\mathbf{x}_0+Q_j\mathbf{y}_j 使上式為零,Q_j^T(\mathbf{b}-A\mathbf{x}_j)=\mathbf{0},即 \mathbf{b}-A\mathbf{x}_j\in N(Q_j^T),也就是說 \mathbf{b}-A\mathbf{x}_j\perp\mathcal{K}_j

 
寫出 FOM 的正規方程

Q_j^TAQ_j\mathbf{y}_j=Q_j^T\mathbf{r}_0

代入 Q_j^T AQ_j=H_j,可得 H_j\mathbf{y}_j=Q_j^T\mathbf{r}_0=\beta Q_j^T\mathbf{q}_1=\beta\mathbf{e}_1。所以,FOM 算出的近似解為 \mathbf{x}_j=\mathbf{x}_0+Q_j\mathbf{y}_j,其中 \mathbf{y}_j=H_j^{-1}(\beta\mathbf{e}_1)。因為 H_j 是上 Hessenberg 矩陣,高斯消去法將 H_j 化為上三角矩陣即可解出 \mathbf{y}_j。FOM 的殘量不需要再另行計算,

\begin{aligned}  \mathbf{b}-A\mathbf{x}_j&=\mathbf{r}_0-AQ_j\mathbf{y}_j\\  &=\mathbf{r}_0-(Q_jH_jH_j^{-1}(\beta{\mathbf{e}}_1)+h_{j+1,j}\mathbf{q}_{j+1}\mathbf{e}_j^T\mathbf{y}_j)\\  &=\mathbf{r}_0-\beta\mathbf{q}_1-h_{j+1,j}\mathbf{q}_{j+1}\mathbf{e}_j^T\mathbf{y}_j\\  &=-h_{j+1,j}(\mathbf{e}_j^T\mathbf{y}_j)\mathbf{q}_{j+1},  \end{aligned}

\Vert\mathbf{b}-A\mathbf{x}_j\Vert=h_{j+1,j}\vert \mathbf{e}^T_j\mathbf{y}_j\vert。FOM 的最小殘差 \mathbf{b}-A\mathbf{x}_j 平行於 \mathbf{q}_{j+1},推論 \mathbf{b}-A\mathbf{x}_j\perp \mathbf{b}-A\mathbf{x}_ii<j。這個性質來自於 FOM 的定義條件,\mathbf{x}_j-\mathbf{x}_0\in\mathcal{K}_j\mathbf{b}-A\mathbf{x}_j\perp\mathcal{K}_j

 
類似 GMRES 的計算步驟,FOM 算法如下:

FOM


  1. \mathbf{r}_0\leftarrow\mathbf{b}-A\mathbf{x}_0\beta\leftarrow \Vert\mathbf{r}_0\Vert\mathbf{q}_1\leftarrow \mathbf{r}_0/\beta
  2. 對於 j=1,2,\ldots,n
  3.     執行 Arnoldi 算法的第 j 次迭代,得到 \tilde{H}_jQ_j
  4.     \mathbf{y}_j\leftarrow H_j^{-1}(\beta{\mathbf{e}}_1)
  5.     \mathbf{x}_j\leftarrow \mathbf{x}_0+Q_j\mathbf{y}_j
  6.     若 \Vert\mathbf{b}-A\mathbf{x}_j\Vert/\Vert\mathbf{b}\Vert<\epsilon,停止計算。

仿造 GMRES(m) 的步驟也可以得到 FOM(m)。

 
GMRES 與 FOM 係針對非對稱可逆矩陣的線性方程解法。如果 A 是一對稱矩陣,則 GMRES 所使用的 Arnoldi 算法簡化為 Lanczos 算法,此法稱為 MINRES (minimal residual method),對應 FOM 的算法稱為 SYMMLQ (symmetric LQ method)[2]。不難想見,MINRES/SYMMLQ 的運算量遠少於 GMRES/FOM。下文將介紹一個適用於對稱正定矩陣的 Krylov 子空間法,稱為共軛梯度法 (conjugate gradient method)。

 
註解
[1] 維基百科:Convergence of GMRES
[2] Iterative methods for symmetric Ax=b

繼續閱讀:
廣告
本篇發表於 線性代數專欄, 數值線性代數 並標籤為 , , , , 。將永久鏈結加入書籤。

8 Responses to Krylov 子空間法──線性方程的數值解法 (二):GMRES 與 FOM

  1. ss 說道:

    老師您好 這邊想請老師幫我看 我對這方法的理解是否有錯

    有些符號用文字代替 內容稍多 可能要麻煩老師

    gmres方法

    找出 Ax=b x近似解 其中A為nxn可逆矩陣 b為(nx1)非0矩陣

    計算部分

    迭代次數 j=1,2,..,n ,自訂起始值X0和誤差

    第j次演算

    先算算出 r0, Ar0,… ,A^j 共j+1個 (r0=b-AX0)

    再做Arnoldi-MGS正交化

    (qr分解,減少比一般GS得到的正交向量誤差,要算到q(j+1))

    得到Qj=(q1,q2,..,qj),R=Hj 其中Hj為Hessenberg矩陣

    並且有 AQj=Q(j+1)Hj’ 其中

    找出某個X(j)使||b-AX(j)|| 最小

    ||b-AX(j)|| = ||Be’1-Hj’y|| (這邊B=beta)

    對Hj’做Givens的QR分解 得到OjHj’=Rj’ 其中Oj為(j+1)乘(j+1)矩陣(這邊O=omega)

    做Givens的目的是要使Rj’變成上三角矩陣且最後一個row為0

    ||Be1′-Hj’y||=||Oj(Be’1-Hj’y)||
    (這邊可以直接乘Oj會相等是因為Oj為么正矩陣所以不影響norm計算?)

    =|g(j+1)|(Oj(Be1′)的最後一個元素)=B|(Oj)j+1,1| (這邊老師寫的看不懂)

    所以可以先找|g(j+1)|是否小於自訂誤差

    如果是 再計算Rjyj=gj 反代找出X(j)近似解

    想法
    Krylov序列得到的空間 Kj(A,r0)

    再做MGS正交化得到的標準正交基底{q1,..,qj}形成的空間與Kj(A,r0)相同

    要找||r0-Az||最小值

    由Cayley-Hamilton 定理知道A^-1是A^(n-1),…,A,I的線性組合

    從 qj 多一次迭代得到q(j+1) 得到的誤差會越小(單調收斂性)

    由最小多項式知道理論上一定找的到一個j<=n 使A^-1是A^(j-1),…,A,I的線性組合

    另外想問老師怎麼用投影概念解是這方法

    (內容好像有提到但不太懂)

    gmres(m)重啟

    直接計算gmres 第m次 得到的xm

    讓起始值為xm再跑一次gmres j=1,2…

  2. Meiyue Shao 說道:

    类似于GMRES/MINRES,FOM的对称版本是SYMMLQ。

  3. K 說道:

    周老師您好 我的論文要寫GMRES相關資料 可不可以引用您的部份文章

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s