7.2 GMRES方法 7.2.1 算法描述 GMRES 方法是目前求解非对称线性方程组的最常用算法之一. 在该算法中, “最佳近似解”的判别方法为“使得 ∥ r m ∥ 2 = ∥ b − A x ( m ) ∥ 2 \| r_m\| _2 = \| b - Ax^{(m)}\| _2 ∥ r m ∥ 2 = ∥ b − A x ( m ) ∥ 2 最小”, 即
x ( m ) = arg min x ∈ x ( 0 ) + K m ∥ b − A x ∥ 2 . (7.6) x ^ {(m)} = \arg \min _ {x \in x ^ {(0)} + \mathcal {K} _ {m}} \| b - A x \| _ {2}. \tag {7.6} x ( m ) = arg x ∈ x ( 0 ) + K m min ∥ b − A x ∥ 2 . ( 7.6 ) 下面我们就根据这个最优性条件来推导出GMRES方法
设迭代初始向量为 x ( 0 ) x^{(0)} x ( 0 ) ,则对任意向量 x ∈ x ( 0 ) + K m x\in x^{(0)} + \mathcal{K}_m x ∈ x ( 0 ) + K m ,可设 x = x ( 0 ) + V m y x = x^{(0)} + V_{m}y x = x ( 0 ) + V m y ,其中 y ∈ R m y\in \mathbb{R}^m y ∈ R m .于是有
r = b − A x = b − A ( x ( 0 ) − V m y ) = r 0 − A V m y = β v 1 − V m + 1 H m + 1 , m y = V m + 1 ( β e 1 − H m + 1 , m y ) , \begin{array}{l} r = b - A x \\ = b - A \left(x ^ {(0)} - V _ {m} y\right) \\ = r _ {0} - A V _ {m} y \\ = \beta v _ {1} - V _ {m + 1} H _ {m + 1, m} y \\ = V _ {m + 1} \left(\beta e _ {1} - H _ {m + 1, m} y\right), \\ \end{array} r = b − A x = b − A ( x ( 0 ) − V m y ) = r 0 − A V m y = β v 1 − V m + 1 H m + 1 , m y = V m + 1 ( β e 1 − H m + 1 , m y ) , 这里 β = ∥ r 0 ∥ 2 \beta = \| r_0\| _2 β = ∥ r 0 ∥ 2 由于 V m + 1 V_{m + 1} V m + 1 列正交,所以
∥ r ∥ 2 = ∥ V m + 1 ( β e 1 − H m + 1 , m y ) ∥ 2 = ∥ β e 1 − H m + 1 , m y ∥ 2 . \| r \| _ {2} = \| V _ {m + 1} (\beta e _ {1} - H _ {m + 1, m} y) \| _ {2} = \| \beta e _ {1} - H _ {m + 1, m} y \| _ {2}. ∥ r ∥ 2 = ∥ V m + 1 ( β e 1 − H m + 1 , m y ) ∥ 2 = ∥ β e 1 − H m + 1 , m y ∥ 2 . 于是最优性条件 (7.6) 就等价于
y ( m ) = arg min y ∈ R m ∥ β e 1 − H m + 1 , m y ∥ 2 . (7.7) y ^ {(m)} = \arg \min _ {y \in \mathbb {R} ^ {m}} \| \beta e _ {1} - H _ {m + 1, m} y \| _ {2}. \tag {7.7} y ( m ) = arg y ∈ R m min ∥ β e 1 − H m + 1 , m y ∥ 2 . ( 7.7 ) 这是一个最小二乘问题. 由于 H m + 1 , m H_{m+1,m} H m + 1 , m 是一个上 Hessenberg 矩阵, 且通常 m m m 不是很大, 所以我们可以用基于 Givens 变换的 QR 分解来求解. 下面就是 GMRES 方法的基本框架.
算法7.4.GMRES方法基本框架 1: 选取初值 x ( 0 ) x^{(0)} x ( 0 ) , 停机标准 ε > 0 \varepsilon > 0 ε > 0 , 以及最大迭代步数 IterMax
2: r 0 = b − A x ( 0 ) , β = ∥ r 0 ∥ 2 r_0 = b - Ax^{(0)},\beta = \| r_0\| _2 r 0 = b − A x ( 0 ) , β = ∥ r 0 ∥ 2 3: v 1 = r 0 / β v_{1} = r_{0} / \beta v 1 = r 0 / β 4: for j = 1 j = 1 j = 1 to IterMax do 5: w = A v j w = Av_{j} w = A v j 6: for i = 1 i = 1 i = 1 to j j j do % Arnoldi 过程 7: h i , j = ( v i , w ) h_{i,j} = (v_i, w) h i , j = ( v i , w ) 8: w = w − h i , j v i w = w - h_{i,j}v_i w = w − h i , j v i 9: end for 10: h j + 1 , j = ∥ w ∥ 2 h_{j + 1,j} = \| w\| _2 h j + 1 , j = ∥ w ∥ 2 11: if h j + 1 , j = 0 h_{j + 1,j} = 0 h j + 1 , j = 0 then 12: m = j , m = j, m = j , break 13: end if
14: v j + 1 = w / h j + 1 , j v_{j + 1} = w / h_{j + 1,j} v j + 1 = w / h j + 1 , j 15: relres = ||rj||2/β %相对残量 16: if relres < ε < \varepsilon < ε then % 检测是否收敛 17: m = j , m = j, m = j , break 18: end if 19: end for 20: 解最小二乘问题 (7.7), 得到 y ( m ) y^{(m)} y ( m ) 21: x ( m ) = x ( 0 ) + V m y ( m ) x^{(m)} = x^{(0)} + V_m y^{(m)} x ( m ) = x ( 0 ) + V m y ( m )
7.2.2 具体实施细节 需要解决的问题有:
(1) 如何计算残量 r j ≜ b − A x ( j ) r_j \triangleq b - Ax^{(j)} r j ≜ b − A x ( j ) 的范数? (2) 如何求解最小二乘问题 (7.7)?
这两个问题可以同时处理. 首先采用 QR 分解来求解最小二乘问题. 设 H m + 1 , m H_{m+1,m} H m + 1 , m 的 QR 分解为
H m + 1 , m = Q m + 1 ⊤ R m + 1 , m , H _ {m + 1, m} = Q _ {m + 1} ^ {\top} R _ {m + 1, m}, H m + 1 , m = Q m + 1 ⊤ R m + 1 , m , 其中 Q m + 1 ∈ R ( m + 1 ) × ( m + 1 ) Q_{m + 1}\in \mathbb{R}^{(m + 1)\times (m + 1)} Q m + 1 ∈ R ( m + 1 ) × ( m + 1 ) 是正交矩阵, R m + 1 , m ∈ R ( m + 1 ) × m R_{m + 1,m}\in \mathbb{R}^{(m + 1)\times m} R m + 1 , m ∈ R ( m + 1 ) × m 是上三角矩阵. 则
∥ β e 1 − H m + 1 , m y ∥ 2 = ∥ β Q m + 1 e 1 − R m + 1 , m y ∥ 2 = ∥ β q 1 − [ R m 0 ] y ∥ 2 , (7.8) \| \beta e _ {1} - H _ {m + 1, m} y \| _ {2} = \| \beta Q _ {m + 1} e _ {1} - R _ {m + 1, m} y \| _ {2} = \left\| \beta q _ {1} - \left[ \begin{array}{c} R _ {m} \\ 0 \end{array} \right] y \right\| _ {2}, \tag {7.8} ∥ β e 1 − H m + 1 , m y ∥ 2 = ∥ β Q m + 1 e 1 − R m + 1 , m y ∥ 2 = β q 1 − [ R m 0 ] y 2 , ( 7.8 ) 其中 R m ∈ R m × m R_{m}\in \mathbb{R}^{m\times m} R m ∈ R m × m 是非奇异上三角矩阵(这里假定 H m + 1 , m H_{m + 1,m} H m + 1 , m 不可约).所以问题(7.7)的解为
y ( m ) = β R m − 1 q 1 ( 1 : m ) , y ^ {(m)} = \beta R _ {m} ^ {- 1} q _ {1} (1: m), y ( m ) = β R m − 1 q 1 ( 1 : m ) , 且
∥ r m ∥ 2 = ∥ b − A x ( m ) ∥ 2 = ∥ β e 1 − H m + 1 , m y ( m ) ∥ 2 = β ⋅ ∣ q 1 ( m + 1 ) ∣ , \| r _ {m} \| _ {2} = \| b - A x ^ {(m)} \| _ {2} = \| \beta e _ {1} - H _ {m + 1, m} y ^ {(m)} \| _ {2} = \beta \cdot | q _ {1} (m + 1) |, ∥ r m ∥ 2 = ∥ b − A x ( m ) ∥ 2 = ∥ β e 1 − H m + 1 , m y ( m ) ∥ 2 = β ⋅ ∣ q 1 ( m + 1 ) ∣ , 其中 q 1 ( m + 1 ) q_{1}(m + 1) q 1 ( m + 1 ) 表示 q 1 q_{1} q 1 的第 m + 1 m + 1 m + 1 个分量
H m + 1 , m H_{m + 1,m} H m + 1 , m 的QR分解的递推计算方法由于 H m + 1 , m H_{m + 1,m} H m + 1 , m 是上Hessenberg矩阵,因此我们采用Givens变换
(1) 当 m = 1 m = 1 m = 1 时, H 21 = [ h 11 h 21 ] H_{21} = \begin{bmatrix} h_{11} \\ h_{21} \end{bmatrix} H 21 = [ h 11 h 21 ] , 构造 Givens 变换 G 1 G_{1} G 1 使得
G 1 H 21 = [ ∗ 0 ] = R 21 , 即 H 21 = G 1 T R 21 . G _ {1} H _ {2 1} = {\left[ \begin{array}{l} {*} \\ 0 \end{array} \right]} = R _ {2 1}, \quad \text {即} \quad H _ {2 1} = G _ {1} ^ {\mathsf {T}} R _ {2 1}. G 1 H 21 = [ ∗ 0 ] = R 21 , 即 H 21 = G 1 T R 21 . (2) 假定存在 G 1 , G 2 , … , G m − 1 G_{1}, G_{2}, \ldots, G_{m-1} G 1 , G 2 , … , G m − 1 , 使得
( G m − 1 … G 2 G 1 ) H m , m − 1 = R m , m − 1 , \left(G _ {m - 1} \dots G _ {2} G _ {1}\right) H _ {m, m - 1} = R _ {m, m - 1}, ( G m − 1 … G 2 G 1 ) H m , m − 1 = R m , m − 1 , 即
H m , m − 1 = ( G m − 1 … G 2 G 1 ) T R m , m − 1 ≜ Q m T R m , m − 1 . H _ {m, m - 1} = \left(G _ {m - 1} \dots G _ {2} G _ {1}\right) ^ {\mathsf {T}} R _ {m, m - 1} \triangleq Q _ {m} ^ {\mathsf {T}} R _ {m, m - 1}. H m , m − 1 = ( G m − 1 … G 2 G 1 ) T R m , m − 1 ≜ Q m T R m , m − 1 . 为了书写方便, 这里假定 G i G_{i} G i 的维数自动扩张, 以满足矩阵乘积的需要.
(3) 考虑 H m + 1 , m H_{m+1,m} H m + 1 , m 的 QR 分解. 易知
H m + 1 , m = [ H m , m − 1 h m 0 h m + 1 , m ] , 其 中 h m = [ h 1 m , h 2 m , … , h m m ] T . H _ {m + 1, m} = \left[ \begin{array}{c c} {{H _ {m, m - 1}}} & {{h _ {m}}} \\ {{0}} & {{h _ {m + 1, m}}} \end{array} \right], \quad \text {其 中} \quad h _ {m} = [ h _ {1 m}, h _ {2 m}, \ldots , h _ {m m} ] ^ {\mathsf {T}}. H m + 1 , m = [ H m , m − 1 0 h m h m + 1 , m ] , 其 中 h m = [ h 1 m , h 2 m , … , h mm ] T . 所以有
[ Q m 0 0 1 ] H m + 1 , m = [ R m , m − 1 Q m h m 0 h m + 1 , m ] = [ R m − 1 h ~ m − 1 0 h ^ m m 0 h m + 1 , m ] , \left[ \begin{array}{c c} Q _ {m} & 0 \\ 0 & 1 \end{array} \right] H _ {m + 1, m} = \left[ \begin{array}{c c} R _ {m, m - 1} & Q _ {m} h _ {m} \\ 0 & h _ {m + 1, m} \end{array} \right] = \left[ \begin{array}{c c} R _ {m - 1} & \tilde {h} _ {m - 1} \\ 0 & \hat {h} _ {m m} \\ 0 & h _ {m + 1, m} \end{array} \right], [ Q m 0 0 1 ] H m + 1 , m = [ R m , m − 1 0 Q m h m h m + 1 , m ] = R m − 1 0 0 h ~ m − 1 h ^ mm h m + 1 , m , 其中 h ~ m − 1 \tilde{h}_{m-1} h ~ m − 1 是 Q m h m Q_{m}h_{m} Q m h m 的前 m − 1 m-1 m − 1 个元素组成的向量, h ^ m m \hat{h}_{mm} h ^ mm 是 Q m h m Q_{m}h_{m} Q m h m 的最后一个元素. 构造Givens变换 G m G_{m} G m :
G m = [ I m − 1 0 0 0 c m s m 0 − s m c m ] ∈ R ( m + 1 ) × ( m + 1 ) , G _ {m} = \left[ \begin{array}{c c c} I _ {m - 1} & 0 & 0 \\ 0 & c _ {m} & s _ {m} \\ 0 & - s _ {m} & c _ {m} \end{array} \right] \in \mathbb {R} ^ {(m + 1) \times (m + 1)}, G m = I m − 1 0 0 0 c m − s m 0 s m c m ∈ R ( m + 1 ) × ( m + 1 ) , 其中
c m = h ^ m , m h ~ m , m , s m = h m + 1 , m h ~ m , m , h ~ m , m = h ^ m , m 2 + h m + 1 , m 2 . c _ {m} = \frac {\hat {h} _ {m , m}}{\tilde {h} _ {m , m}}, s _ {m} = \frac {h _ {m + 1 , m}}{\tilde {h} _ {m , m}}, \tilde {h} _ {m, m} = \sqrt {\hat {h} _ {m , m} ^ {2} + h _ {m + 1 , m} ^ {2}}. c m = h ~ m , m h ^ m , m , s m = h ~ m , m h m + 1 , m , h ~ m , m = h ^ m , m 2 + h m + 1 , m 2 . 令
Q m + 1 = G m [ Q m 0 0 1 ] , Q _ {m + 1} = G _ {m} \left[ \begin{array}{c c} Q _ {m} & 0 \\ 0 & 1 \end{array} \right], Q m + 1 = G m [ Q m 0 0 1 ] , 则
Q m + 1 H m + 1 , m = G m [ R m − 1 h ~ m − 1 0 h ^ m , m 0 h m + 1 , m ] = [ R m − 1 h ~ m − 1 0 h ~ m , m 0 0 ] ≜ R m + 1 , m . Q _ {m + 1} H _ {m + 1, m} = G _ {m} \left[ \begin{array}{c c} R _ {m - 1} & \tilde {h} _ {m - 1} \\ 0 & \hat {h} _ {m, m} \\ 0 & h _ {m + 1, m} \end{array} \right] = \left[ \begin{array}{c c} R _ {m - 1} & \tilde {h} _ {m - 1} \\ 0 & \tilde {h} _ {m, m} \\ 0 & 0 \end{array} \right] \triangleq R _ {m + 1, m}. Q m + 1 H m + 1 , m = G m R m − 1 0 0 h ~ m − 1 h ^ m , m h m + 1 , m = R m − 1 0 0 h ~ m − 1 h ~ m , m 0 ≜ R m + 1 , m . 所以可得 H m + 1 , m H_{m + 1,m} H m + 1 , m 的QR分解 H m + 1 , m = Q m + 1 T R m + 1 , m H_{m + 1,m} = Q_{m + 1}^{\mathsf{T}}R_{m + 1,m} H m + 1 , m = Q m + 1 T R m + 1 , m
由 H m , m − 1 H_{m,m - 1} H m , m − 1 的QR分解到 H m + 1 , m H_{m + 1,m} H m + 1 , m 的QR分解,我们需要
(1) 计算 Q m h m Q_{m}h_{m} Q m h m ,即将之前的 m − 1 m - 1 m − 1 个Givens变换作用到 H m + 1 , m H_{m + 1,m} H m + 1 , m 的最后一列的前 m m m 个元素上,所以我们需要保留所有的Givens变换; (2) 残量计算: ∥ r m ∥ 2 = ∣ β q 1 ( m + 1 ) ∣ = ∣ β Q m + 1 ( m + 1 , 1 ) ∣ \| r_{m}\|_{2} = |\beta q_{1}(m + 1)| = |\beta Q_{m + 1}(m + 1,1)| ∥ r m ∥ 2 = ∣ β q 1 ( m + 1 ) ∣ = ∣ β Q m + 1 ( m + 1 , 1 ) ∣ , 即
G m G m − 1 … G 2 G 1 ( β e 1 ) G _ {m} G _ {m - 1} \dots G _ {2} G _ {1} (\beta e _ {1}) G m G m − 1 … G 2 G 1 ( β e 1 ) 的最后一个分量的绝对值. 由于在计算 r m − 1 r_{m-1} r m − 1 时就已经计算出 G m − 1 ⋯ G 2 G 1 ( β e 1 ) G_{m-1} \cdots G_2 G_1(\beta e_1) G m − 1 ⋯ G 2 G 1 ( β e 1 ) , 因此这里只需做一次 Givens 变换即可;
(3) y ( m ) y^{(m)} y ( m ) 的计算:当相对残量满足精度要求时,需要计算 y ( m ) = R m − 1 q 1 ( 1 : m ) y^{(m)} = R_m^{-1}q_1(1:m) y ( m ) = R m − 1 q 1 ( 1 : m ) ,而 q 1 q_{1} q 1 即为 G m G m − 1 … G 2 G 1 ( β e 1 ) G_{m}G_{m - 1}\dots G_{2}G_{1}(\beta e_{1}) G m G m − 1 … G 2 G 1 ( β e 1 )
算法7.5.实用GMRES方法 20 : [ h i j h i + 1 , j ] = [ c i s i − s i c i ] [ h i j h i + 1 , j ] 2 0: \quad \left[ \begin{array}{c} h _ {i j} \\ h _ {i + 1, j} \end{array} \right] = \left[ \begin{array}{c c} c _ {i} & s _ {i} \\ - s _ {i} & c _ {i} \end{array} \right] \left[ \begin{array}{c} h _ {i j} \\ h _ {i + 1, j} \end{array} \right] 20 : [ h ij h i + 1 , j ] = [ c i − s i s i c i ] [ h ij h i + 1 , j ] 25 : τ = h j j / h j + 1 , j , s j = 1 / 1 + τ 2 , c j = s j τ 2 5: \quad \tau = h _ {j j} / h _ {j + 1, j}, s _ {j} = 1 / \sqrt {1 + \tau^ {2}}, c _ {j} = s _ {j} \tau 25 : τ = h jj / h j + 1 , j , s j = 1/ 1 + τ 2 , c j = s j τ 1: 选取初值 x ( 0 ) x^{(0)} x ( 0 ) , 停机标准 ε > 0 \varepsilon > 0 ε > 0 , 以及最大迭代步数 IterMax 2: r 0 = b − A x ( 0 ) , β = ∥ r 0 ∥ 2 r_0 = b - Ax^{(0)},\beta = \| r_0\| _2 r 0 = b − A x ( 0 ) , β = ∥ r 0 ∥ 2 3: if β / ∥ b ∥ 2 < ε \beta / \| b \|_2 < \varepsilon β /∥ b ∥ 2 < ε then 4: 停止计算, 输出近似解 x ( 0 ) x^{(0)} x ( 0 ) 5: end if 6: v 1 = r 0 / β v_{1} = r_{0} / \beta v 1 = r 0 / β 7: ξ = β e 1 \xi = \beta e_{1} ξ = β e 1 8: for j = 1 j = 1 j = 1 to IterMax do 9: w = A v j w = Av_{j} w = A v j 10: for i = 1 i = 1 i = 1 to j j j do % \% % Arnoldi过程 11: h i , j = ( v i , w ) h_{i,j} = (v_i,w) h i , j = ( v i , w ) 12: w = w − h i , j v i w = w - h_{i,j}v_{i} w = w − h i , j v i 13: end for 14: h j + 1 , j = ∥ w ∥ 2 h_{j + 1,j} = \| w\| _2 h j + 1 , j = ∥ w ∥ 2 15: if h j + 1 , j = 0 h_{j+1,j} = 0 h j + 1 , j = 0 then % 迭代中断 16: m = j m = j m = j , break 17: end if 18: v j + 1 = w / h j + 1 , j v_{j + 1} = w / h_{j + 1,j} v j + 1 = w / h j + 1 , j 19: f o r i = 1 \mathbf{for}i = 1 for i = 1 to j − 1 j - 1 j − 1 do % \% % 计算 G j − 1 … G 2 G 1 H j + 1 , j ( 1 : j , j ) G_{j - 1}\dots G_2G_1H_{j + 1,j}(1:j,j) G j − 1 … G 2 G 1 H j + 1 , j ( 1 : j , j ) 21: end for 22: if ∣ h j j ∣ > ∣ h j + 1 , j ∣ |h_{jj}| > |h_{j+1,j}| ∣ h jj ∣ > ∣ h j + 1 , j ∣ then % 构造 Givens 变换 G j G_j G j 23: τ = h j + 1 , j / h j j , c j = 1 / 1 + τ 2 , s j = c j τ \tau = h_{j + 1,j} / h_{jj},c_{j} = 1 / \sqrt{1 + \tau^{2}},s_{j} = c_{j}\tau τ = h j + 1 , j / h jj , c j = 1/ 1 + τ 2 , s j = c j τ 24: else 26: end if 27: h j j = c j h j j + s j h j + 1 , j h_{jj} = c_j h_{jj} + s_j h_{j+1,j} h jj = c j h jj + s j h j + 1 , j % \% % 计算 G j H j + 1 , j ( 1 : j , j ) G_j H_{j+1,j}(1:j,j) G j H j + 1 , j ( 1 : j , j ) 28: h j + 1 , j = 0 h_{j + 1,j} = 0 h j + 1 , j = 0 [ ξ j ξ j + 1 ] = [ c j s j − s j c j ] [ ξ j 0 ] \left[ \begin{array}{l}\xi_{j}\\ \xi_{j + 1} \end{array} \right] = \left[ \begin{array}{ll}c_{j} & s_{j}\\ -s_{j} & c_{j} \end{array} \right]\left[ \begin{array}{l}\xi_{j}\\ 0 \end{array} \right] [ ξ j ξ j + 1 ] = [ c j − s j s j c j ] [ ξ j 0 ] % 计算 G j ( β G j − 1 … G 2 G 1 e 1 ) G_{j}(\beta G_{j - 1}\dots G_{2}G_{1}e_{1}) G j ( β G j − 1 … G 2 G 1 e 1 ) 30: relres = |ξj+1|/β %相对残量 31: if relres < ε < \varepsilon < ε then
32: m = j , m = j, m = j , break 33: end if 34: end for 35: m = j m = j m = j 36: y ( m ) = H ( 1 : m , 1 : m ) \ ξ ( 1 : m ) y^{(m)} = H(1:m, 1:m) \backslash \xi(1:m) y ( m ) = H ( 1 : m , 1 : m ) \ ξ ( 1 : m ) % \% % 求最小二乘问题的解, 回代求解 37: x ( m ) = x ( 0 ) + V m y ( m ) x^{(m)} = x^{(0)} + V_m y^{(m)} x ( m ) = x ( 0 ) + V m y ( m ) 38: if relres < ε < \varepsilon < ε then 39: 输出近似解 x x x 及相关信息 40: else 41: 输出算法失败信息 42: end if
7.2.3 GMRES 方法的中断 在上面的GMRES方法中,当执行到某一步时有 h k + 1 , k = 0 h_{k + 1,k} = 0 h k + 1 , k = 0 ,则算法会中断(breakdown).如果出现这种中断,则我们就找到了精确解.
定理7.4 设 A ∈ R n × n A \in \mathbb{R}^{n \times n} A ∈ R n × n 非奇异且 r 0 ≠ 0 r_0 \neq 0 r 0 = 0 . 若 h i + 1 , i ≠ 0 , i = 1 , 2 , … , k − 1 h_{i+1,i} \neq 0, i = 1,2,\ldots,k-1 h i + 1 , i = 0 , i = 1 , 2 , … , k − 1 , 则 h k + 1 , k = 0 h_{k+1,k} = 0 h k + 1 , k = 0 当且仅当 x ( k ) x^{(k)} x ( k ) 是方程组的精确解. (不考虑舍入误差) (板书)
证明. 设 h k + 1 , k = 0 h_{k+1,k} = 0 h k + 1 , k = 0 , 则有
A V k = V k H k , y ( k ) = H k − 1 ( β e 1 ) . A V _ {k} = V _ {k} H _ {k}, \quad y ^ {(k)} = H _ {k} ^ {- 1} (\beta e _ {1}). A V k = V k H k , y ( k ) = H k − 1 ( β e 1 ) . 所以
∥ r k ∥ 2 = ∥ b − A x ( k ) ∥ 2 = ∥ b − A ( x ( 0 ) + V k y ( k ) ) ∥ 2 = ∥ r 0 − V k H k y ( k ) ∥ 2 = ∥ β v 1 − V k ( β e 1 ) ∥ 2 = 0. \begin{array}{l} \| r _ {k} \| _ {2} = \| b - A x ^ {(k)} \| _ {2} = \| b - A \left(x ^ {(0)} + V _ {k} y ^ {(k)}\right) \| _ {2} \\ = \| r _ {0} - V _ {k} H _ {k} y ^ {(k)} \| _ {2} = \| \beta v _ {1} - V _ {k} (\beta e _ {1}) \| _ {2} = 0. \\ \end{array} ∥ r k ∥ 2 = ∥ b − A x ( k ) ∥ 2 = ∥ b − A ( x ( 0 ) + V k y ( k ) ) ∥ 2 = ∥ r 0 − V k H k y ( k ) ∥ 2 = ∥ β v 1 − V k ( β e 1 ) ∥ 2 = 0. 反之,设 x ( k ) x^{(k)} x ( k ) 是精确解,则
0 = b − A x ( k ) = r 0 − V k + 1 H k + 1 , k y ( k ) = V k + 1 ( β e 1 − H k + 1 , k y ( k ) ) . 0 = b - A x ^ {(k)} = r _ {0} - V _ {k + 1} H _ {k + 1, k} y ^ {(k)} = V _ {k + 1} \left(\beta e _ {1} - H _ {k + 1, k} y ^ {(k)}\right). 0 = b − A x ( k ) = r 0 − V k + 1 H k + 1 , k y ( k ) = V k + 1 ( β e 1 − H k + 1 , k y ( k ) ) . 反证法, 假设 h k + 1 , k ≠ 0 h_{k+1,k} \neq 0 h k + 1 , k = 0 , 则 v k + 1 ≠ 0 v_{k+1} \neq 0 v k + 1 = 0 . 因此 V k + 1 V_{k+1} V k + 1 单位列正交, 故列满秩, 所以由上式可知
β e 1 − H k + 1 , k y ( k ) = 0. \beta e _ {1} - H _ {k + 1, k} y ^ {(k)} = 0. β e 1 − H k + 1 , k y ( k ) = 0. 由于 H k + 1 , k H_{k + 1,k} H k + 1 , k 是上Hessenberg矩阵,且 h i + 1 , i ≠ 0 , i = 1 , 2 , … , k . h_{i + 1,i}\neq 0,i = 1,2,\ldots ,k. h i + 1 , i = 0 , i = 1 , 2 , … , k . 通过向后回代求解可得 y ( k ) = 0 y^{(k)} = 0 y ( k ) = 0 于是 β = 0. \beta = 0. β = 0. 这与 r 0 ≠ 0 r_0\neq 0 r 0 = 0 矛盾.所以 h k + 1 , k = 0 h_{k + 1,k} = 0 h k + 1 , k = 0
7.2.4 带重启的GMRES方法 由于随着迭代步数的增加, GMRES 方法的每一步所需的运算量和存储量都会越来越大. 因此当迭代步数很大时, GMRES 方法就不太实用. 通常的解决方法就是重启, 即事先设定一个重启迭代步数 k k k , 如 k = 20 k = 20 k = 20 或 50 等等, 当 GMRES 达到这个迭代步数时仍不收敛, 则计算出方程组在 x ( 0 ) + K k x^{(0)} + \mathcal{K}_k x ( 0 ) + K k 中的最佳近似解 x ( k ) x^{(k)} x ( k ) , 然后令 x ( 0 ) = x ( k ) x^{(0)} = x^{(k)} x ( 0 ) = x ( k ) , 并重新开始新的 GMRES 迭代. 不断重复该过程, 直到收敛为止.
算法7.6.GMRES ( k ) (k) ( k ) :带重启的GMRES方法 1: 设定重启步数 k ( ≪ n ) k \ (\ll n) k ( ≪ n ) 2: 选取初值 x ( 0 ) x^{(0)} x ( 0 ) , 停机标准 ε > 0 \varepsilon > 0 ε > 0 , 以及最大迭代 3: r 0 = b − A x ( 0 ) , β = ∥ r 0 ∥ 2 r_0 = b - Ax^{(0)}, \beta = \| r_0\|_2 r 0 = b − A x ( 0 ) , β = ∥ r 0 ∥ 2 4: if β / ∥ b ∥ 2 < ε \beta / \| b \|_2 < \varepsilon β /∥ b ∥ 2 < ε then 5: 停止计算, 输出近似解 x = x ( 0 ) x = x^{(0)} x = x ( 0 ) 6: end if 7: for iter=1 to ceil(IterMax/k) do % 外循环 8: v 1 = r 0 / β v_1 = r_0 / \beta v 1 = r 0 / β 9: ξ = β e 1 \xi = \beta e_1 ξ = β e 1 10: for j = 1 j = 1 j = 1 to k k k do 11: 调用 GMRES 循环 12: end for 13: m = j m = j m = j 14: y ( m ) = H ( 1 : m , 1 : m ) \ ξ ( 1 : m ) y^{(m)} = H(1:m, 1:m) \backslash \xi(1:m) y ( m ) = H ( 1 : m , 1 : m ) \ ξ ( 1 : m ) 15: x ( m ) = x ( 0 ) + V m y ( m ) x^{(m)} = x^{(0)} + V_m y^{(m)} x ( m ) = x ( 0 ) + V m y ( m ) 16: if relres < ε \varepsilon ε then 17: break 18: end if 19: x ( 0 ) = x ( m ) x^{(0)} = x^{(m)} x ( 0 ) = x ( m ) % 重启 GMRES 20: r 0 = b − A x ( 0 ) , β = ∥ r 0 ∥ 2 r_0 = b - Ax^{(0)}, \beta = \| r_0\|_2 r 0 = b − A x ( 0 ) , β = ∥ r 0 ∥ 2 21: end for 22: if relres < ε \varepsilon ε then 23: 输出近似解 x ( m ) x^{(m)} x ( m ) 及相关信息 24: else 25: 输出算法失败信息 26: end if
带重启的GMRES方法需要注意的问题:
(1) 如何选取合适的重启步数 k k k ? 一般只能依靠经验来选取; (2) 不带重启的 GMRES 方法能保证算法的收敛性, 但带重启的 GMRES 方法却无法保证, 有时可能出现停滞现象 (stagnation).