A.12 GMRES 算法

A.12 GMRES 算法

考虑线性方程组

$$ \bar{A}\vec{x}=\vec{b}.\tag{A.90} $$

寻求以下形式的近似解：

$$ \vec{x}=\vec{x}_{0}+\vec{z},\tag{A.91} $$

其中，$ \vec{x}_{0} $ 为初始估计，$ \vec{z} $ 属于 Krylov 子空间

$$ \vec{z}\in\mathcal{K}_{m},\quad\mathcal{K}_{m}\equiv\operatorname{span}\left\{\vec{r}_{0},\;\bar{A}\vec{r}_{0},\;\bar{A}^{2}\vec{r}_{0},\;\ldots,\;\bar{A}^{m-1}\vec{r}_{0}\right\}\tag{A.92} $$

其中，$ \vec{r}_{0}=\vec{b}-\vec{A}\vec{x}_{0} $，m 为 K 的维数，也称搜索方向数。广义最小残差（GMRES）算法 [27] 确定 $ \vec{z} $，使残差的 2-范数，即

$$ \parallel\vec{b}-\bar{A}(\vec{x}_{0}+\vec{z})\parallel\tag{A.93} $$

达到最小。下面给出 GMRES 算法的具体步骤。

A.12.1 计算 $ K_{m} $ 的标准正交基

采用修正的 Gram–Schmidt 过程：

$$ \begin{aligned}\vec r_0&=\vec b-\bar A\vec x_0,\\\vec v_1&=\vec r_0/\|\vec r_0\|,\\&\text{对 }j=1,\ldots,m\text{ 循环：}\\\tilde{\vec v}_{j+1}&=\bar A\vec v_j,\\&\quad\text{对 }i=1,\ldots,j\text{ 循环：}\\h_{i,j}&=\tilde{\vec v}_{j+1}\cdot\vec v_i,\\\tilde{\vec v}_{j+1}&=\tilde{\vec v}_{j+1}-h_{i,j}\vec v_i,\\&\quad\text{结束内层循环；}\\h_{j+1,j}&=\|\tilde{\vec v}_{j+1}\|,\\\vec v_{j+1}&=\tilde{\vec v}_{j+1}/h_{j+1,j},\\&\text{结束外层循环。}\end{aligned} $$

其中，$ h_{i,j} $ 为上 Hessenberg 矩阵的系数，i 表示行，j 表示列。这里将矩阵扩展，加入元素 $ h_{j+1,j} $，因此维数为 $ (m+1) \times m $。

A.12.2 构造上 Hessenberg 矩阵

该矩阵具有以下近似三角形式：

$$ \bar{H}_{m}^{*}=\begin{bmatrix}h_{1,1}&h_{1,2}&\cdots&h_{1,m-1}&h_{1,m}\\ h_{2,1}&h_{2,2}&\cdots&h_{2,m-1}&h_{2,m}\\ 0&h_{3,2}&\ddots&\vdots&\vdots\\ \vdots&0&\ddots&h_{m-1,m-1}&h_{m-1,m}\\ \vdots&\vdots&\ddots&h_{m,m-1}&h_{m,m}\\ 0&0&\cdots&0&h_{m+1,m}\end{bmatrix}^{(m+1)\times m}.\tag{A.94} $$

下面将利用它建立并求解式（A.93）的残差最小化问题。

A.12.3 最小化残差

初始解 $ \vec{x}_{0} $ 的修正量定义为 [27]

$$ \vec{z}=\sum_{j=1}^{m} y_{j}\vec{ v}_{j},\tag{A.95} $$

其中，$ y_{j} $ 是以下矢量的分量：

$$ \vec{ y}=\left[ y_{1}, y_{2},\ldots, y_{m}\right]^{\mathrm{T}}.\tag{A.96} $$

此外，可以证明

$$ \bar{A}\bar{V}_{m}=\bar{V}_{m+1}\bar{H}_{m}^{*}\tag{A.97} $$

其中，

$$ \bar{V}_{m}=[\vec{ v}_{1},\vec{ v}_{2},\ldots,\vec{ v}_{m}]\tag{A.98} $$

是以 $ \vec{ v}_{j} $ 为各列的矩阵。引入记号

$$ \vec{e}=[\parallel\vec{r}_{0}\parallel,0,\ldots,0]^{\mathrm{T}},\tag{A.99} $$

其中，$ \vec{e} $ 具有 $ (m+1) $ 个元素。由式（A.99）可见，$ \vec{r}_{0}=\vec{b}-\vec{A}\vec{x}_{0}=\vec{V}_{m+1}\vec{e} $。因此，式（A.93）的残差为

$$ \begin{aligned}\parallel\vec{b}-\bar{A}(\vec{x}_{0}+\vec{z})\parallel&=\parallel\vec{r}_{0}-\bar{A}\left(\sum_{j=1}^{m} y_{j}\vec{v}_{j}\right)\parallel\\&=\parallel\vec{r}_{0}-\bar{A}\bar{V}_{m}\vec{ y}\parallel\\&=\parallel\bar{V}_{m+1}(\vec{e}-\bar{H}_{m}^{*}\vec{ y})\parallel\\&=\parallel\vec{e}-\bar{H}_{m}^{*}\vec{ y}\parallel.\end{aligned}\tag{A.100} $$

最后一步采用了 $ \bar{V}_{m+1} $ 的标准正交性，即 $ \vec{ v}_{i}^{T} \cdot \vec{ v}_{j} = 0 $ 在 $ i \neq j $ 时成立，$ \vec{ v}_{i}^{T} \cdot \vec{ v}_{j} = 1 $ 在 $ i = j $ 时成立。因此，残差最小化问题可简化为

$$ \min_{\vec{z}\in\mathcal{K}_{m}}\parallel\vec{b}-\bar{A}(\vec{x}_{0}+\vec{z})\parallel=\min_{\vec{ y}\in\mathcal{R}^{m}}\parallel\vec{e}-\bar{H}_{m}^{*}\vec{ y}\parallel.\tag{A.101} $$

该最小化问题可采用下面介绍的 Q–R 算法求解。

A.12.4 Q–R 算法

定义 $ \bar{R}_{m} = \bar{Q}_{m}\bar{H}_{m}^{*} $，其中

$$ \bar{Q}_{m}\stackrel{\mathrm{def}}{=}\bar{F}_{m}\bar{F}_{m-1}\ldots\bar{F}_{1}\tag{A.102} $$

为以下 Givens 旋转矩阵的乘积：

$$ \bar F_j=\begin{bmatrix}\bar I_{j-1}&0&0\\0&\begin{matrix}c_j&s_j\\-s_j&c_j\end{matrix}&0\\0&0&\bar I_{m-j}\end{bmatrix}^{(m+1)\times(m+1)}.\tag{A.103} $$

式（A.103）中，$ \bar{I}_{j} $ 为 j 维单位矩阵，$ c_{j} $ 与 $ s_{j} $（$ c_{j}^{2} + s_{j}^{2} = 1 $）分别表示旋转角的正弦与余弦。选择旋转，使 $ \bar{H}_{m}^{*} $ 变换为上三角矩阵 $ \bar{R}_{m} $，其维数为 $ (m + 1) \times m $，最后一行全为零。由于 $ \bar{Q}_{m}^{T} \bar{Q}_{m} = \bar{I} $，式（A.101）可写为

$$ \begin{aligned}\|\vec{e}-\bar{H}_{m}^{*}\vec{ y}\|=&\|\bar{Q}_{m}^{\mathrm{T}}(\bar{Q}_{m}\vec{e}-\bar{Q}_{m}\bar{H}_{m}^{*}\vec{ y})\|\\=&\|\vec{g}-\bar{R}_{m}\vec{ y}\|,\end{aligned}\tag{A.104} $$

其中，$ \vec{g}=\bar{Q}_{m}\vec{e} $ 为矢量 $ \vec{e} $（式（A.99））变换后的结果。$ \bar{R}_{m} $ 的最后一行全为零，因此只有项 $ g_{m+1} $ 在矢量的第 $ (m+1) $ 行非零，该矢量为 $ (\vec{g}-\bar{R}_{m}\vec{y}) $。将 $ (\vec{g}-\bar{R}_{m}\vec{y}) $ 的前 m 个分量记为 $ p_{j} $（$ j=1,\ldots,m $），则式（A.104）的范数为

$$ \|\vec{g}-\bar{R}_{m}\vec{ y}\|=\sqrt{g_{m+1}^{2}+\sum_{j=1}^{m}p_{j}^{2}}.\tag{A.105} $$

选择各分量 $ y_{j} $，组成 $ \vec{ y} $，使 $ p_{j}^{2}=0 $ 对所有 $ j=1,\ldots,m $ 成立，则式（A.101）的最小化问题得到

$$ \min_{\vec{ y}\in\mathcal{R}^{m}}\parallel\vec{g}-\bar{R}_{m}\vec{ y}\parallel=\left|g_{m+1}\right|.\tag{A.106} $$

各分量 $ y_{j} $ 可通过回代求解以下线性方程组获得：

$$ \begin{bmatrix}R_{1,1}&\cdots&R_{1,m-1}&R_{1,m}\\ 0&\ddots&\vdots&\vdots\\ \vdots&\ddots&R_{m-1,m-1}&R_{m-1,m}\\ 0&\cdots&0&R_{m,m}\end{bmatrix}\begin{bmatrix} y_{1}\\ \vdots\\ y_{m-1}\\ y_{m}\end{bmatrix}=\begin{bmatrix}g_{1}\\ \vdots\\ g_{m-1}\\ g_{m}\end{bmatrix}\tag{A.107} $$

已知 $ y_{j} $ 后，由式（A.95）得到方程组（A.90）的解。

需要注意，

$$ \begin{aligned}\parallel\vec{r}_{m}\parallel=&\parallel\vec{b}-\bar{A}(\vec{x}_{0}+\vec{z})\parallel\\=&\parallel\vec{e}-\bar{H}_{m}^{*}\vec{ y}\parallel\\=&\parallel\vec{g}-\bar{R}_{m}\vec{ y}\parallel\\=&\left|g_{m+1}\right|.\end{aligned}\tag{A.108} $$

这意味着，实际残差可容易地由 $ \left|g_{m+1}\right| $ 确定。