2/115 大规模稀疏线性方程组 Ax = b, A R n n, b R n. Krylov 子空间方法是子空间方法的成功代表 首选方法 Krylov 子空间方法. 基本思想 在一个维数较小的子空间 K R n 中寻找近似解. 这类方法也被看作是一种 投影方法, 即寻找真解在某个子空间中 的投影

Transcription

1 第四讲 Krylov 子空间方法 投影方法 Krylov 子空间与 Arnoldi 过程 一般线性方程组的 Krylov 子空间方法 对称线性方程的 Krylov 子空间方法 收敛性分析 基于双正交化过程的迭代方法 免转置迭代方法 正规方程的迭代方法

2 2/115 大规模稀疏线性方程组 Ax = b, A R n n, b R n. Krylov 子空间方法是子空间方法的成功代表 首选方法 Krylov 子空间方法. 基本思想 在一个维数较小的子空间 K R n 中寻找近似解. 这类方法也被看作是一种 投影方法, 即寻找真解在某个子空间中 的投影 ( 可以是正交投影, 也可以是斜投影 ) 如果没有特别注明, 本讲内容都是在实数域中讨论.

3 3/115 1 投影方法 设 K 是 R n 的一个子空间, 记 m dim(k) n 目标在 K 中寻找真解的一个 最佳 近似. 定解条件设置 m 个约束 : 要求残量满足 m 个正交性条件, 即 r = b A x L (4.1) 其中 x 是近似解, L 是另一个 m 维子空间. 不同的 L 对应不同的投影方法 m 维空间的向量有 m 个自由度, 因此需要加 m 的限制条件 Petrov-Galerkin 条件 Galerkin 条件 : L = K L 称为 K 称为 约束空间 搜索空间 当 L = K 时, 我们称为正交投影法, 否则称为斜投影法

4 4/115 投影方法的数学描述 find x K such that b A x L (4.2) 若给定初值 x (0) R n, 则改用仿射空间 x (0) + K, 即 find x x (0) + K such that b A x L. (4.3) 好的初值一般都包含有价值的信息 事实上, 如果将 x 写成 : x = x (0) + ˆx, 其中 ˆx K, 则 (4.3) 就等价于 find ˆx K such that r 0 Aˆx L, (4.4) 求解 x 等价于求解 ˆx 其中 r 0 = b Ax (0) 是初始残量. 这与 (4.2) 的形式是一样的.

5 5/115 如何计算近似解 设 {v 1, v 2,..., v m } 和 {w 1, w 2,..., w m } 分别是 K 和 L 一组基. 令 V = [v 1, v 2,..., v m ], W = [w 1, w 2,..., w m ]. W R n m, V R n m 由于 x x (0) + K, 因此存在向量 y R m 使得 x = x (0) + V y ˆx x x (0) = V y 由正交性条件 (4.4) 可知 r 0 AV y w i, i = 1, 2,..., m, 即 W AV y = W r 0. 若 W AV 非奇异, 则解得 y = (W AV ) 1 W r 0. 因此 在实际计算中, W AV 通常可以直接形成, 无需另外计算. x = x (0) + V (W AV ) 1 W r 0.

6 6/115 近似解 x 存在唯一的条件是矩阵 W AV 非奇异. 定理 1 如果 A, K 和 L 满足下面条件之一 (1) A 正定且 L = K; (2) A 非奇异且 L = AK, 证明留作练习 则矩阵 W AV 非奇异.

7 7/115 如何实施投影方法在实施投影方法时, 我们需要考虑下面两个问题 : 怎样选择搜索空间 K 和约束空间 L? 如果找到的近似解 x 达不到精度要求, 则需要更换搜索空间 K 和约束空间 L, 如何更新? 目前能够很好地解决上面两个问题的方法就是 Krylov 子空间方法

8 8/115 2 Krylov 子空间与 Arnoldi 过程 Krylov 子空间 设 A R n n, r R n, 我们称 K m (A, r) span{r, Ar,..., A m 1 r} R n Krylov 子空间方法就是在 Krylov 子空间中寻找近似解 是由 A 和 r 生成的 Krylov 子空间. 为书写方便, 通常简记为 K m. 几个基本性质 : Krylov 子空间是嵌套的, 即 : K 1 K 2 K m K m 的维数不超过 m { } K m (A, r) = x = p(a)r : p 为次数不超过 m 1 的多项式

9 9/115 约束子空间 L 的选取 目前被广泛采用的选取方法有 : L = K : 如 FOM, CG, SYMMLQ L = AK : 如 MINRES, GMRES 通过对这些方法进行变形与改进, 可以构造更多的子空间方法. L = K(A, r) : 如 BiCG

10 10/115 Arnoldi 过程 : 计算 K m 的一组正交基 算法 2.1 基于 Gram-Schmidt 正交化的 Arnoldi 过程 1: 给定非零向量 r, 计算 v 1 = r/ r 2 2: for j = 1, 2,..., m 1 do 3: w j = Av j 4: for i = 1, 2,..., j do 5: h ij = (w j, v i ) 6: end for 7: w j = w j j h ij v i i=1 8: h j+1,j = w j 2 9: if h j+1,j = 0 then 10: break 11: end if 12: v j+1 = w j /h j+1,j 13: end for 做计算时一定要注意任何可能存在的 中断.

11 11/115 如果到第 k (k < m) 步时有 h k+1,k = 0, 则算法将提前终止. 此时 Av k 必定可由 v 1, v 2,..., v k 线性表出 ( 不考虑浮点运算误差 ) 向量 v i 称为 Arnoldi 向量 需要注意的是, 在算法中, 我们是用 A 乘以 v j, 然后与之前的 Arnoldi 向量正交, 而不是计算 A j r. 事实上, 它们是等价的. 引理 1 设算法 2.1 不提前终止 ( 即 h k+1,k 0, k = 1, 2,..., m 1), 则 v j K j (A, r), j = 1, 2,..., m. K j (A, r) = span{v 1,..., v j } ( 板书 )

12 12/115 根据上面的引理, 我们可以立即得到下面的结论. 定理 2 如果算法 2.1 不提前终止, 则向量 v 1, v 2,..., v m 构成 K m 的一组标准正交基, 其中 K m (A, r) = span{r, Ar,..., A m 1 r}.

13 13/115 Arnoldi 过程的重要性质 由 Arnoldi 过程可知 h j+1,j v j+1 = Av j j h ij v i, 因此 Av j = h j+1,j v j+1 + j h ij v i i=1 = [v 1, v 2,..., v j+1 ] = V m+1 H m+1,m (:, j), h 1j h 2j. h j+1,j i=1 = [v 1,..., v j+1, v j+2,..., v m+1 ] h 1j. h j+1,j 0. 0 w j = w j j i=1 h ij v i v j+1 = w j /h j+1,j

14 14/115 其中 V m+1 = [v 1, v 2,..., v m+1 ], h 11 h 12 h 13 h 1,m 1 h 1,m h 21 h 22 h 23 h 2,m 1 h 2,m 0 h 32 h 33 h 3,m 1 h 2,m H m+1,m = 0 0 h 43 h 4,m 1 h 4,m h m,m 1 h m,m h m+1,m R (m+1) m w j = Av j for i = 1, 2,..., j do h ij = (w j, v i ) end for w j = w j j i=1 h j+1,j = w j 2 h ij v i 这里 h ij 是由 Arnoldi 过程所定义的.

15 15/115 定理 3 设 V m = [v 1, v 2,..., v m ], 则 AV m = V m+1 H m+1,m = V m H m + h m+1,m v m+1 e m (4.5) VmAV m = H m (4.6) 这两个公式在后面的算法推导过程中非常重要 其中 e m = [0,..., 0, 1] R m, H m = H m+1,m (1 : m, 1 : m) R m m, 即 H m 是由 H m+1,m 的前 m 行组成的 Hessenberg 矩阵. H m = = A V m V m+1 H m+1,m V m +

16 16/115 如果 Arnoldi 方法提前终止, 则我们可以得到一个不变子空间. 定理 4 如果算法 2.1 在第 k (k m) 步终止, 即 h k+1,k = 0, 则有 即 K k 是 A 的一个不变子空间. AV k = V k H k, 一旦出现不变子空间, Arnoldi 过程就会被中断, 这是好事还是坏事?

17 17/115 考虑到数值稳定性, 通常使用 修正的 Gram-Schmidt 过程 (MGS) 算法 2.2 基于 MGS 的 Arnoldi 过程 1: 给定非零向量 r 2: v 1 = r/ r 2 3: for j = 1, 2,..., m 1 do 4: w j = Av j 5: for i = 1, 2,..., j do 6: h ij = (w j, v i ) 7: w j = w j h ij v i 8: end for w j = Av j for i = 1, 2,..., j do h ij = (w j, v i ) end for w j = w j j i=1 h ij v i 9: h j+1,j = w j 2 10: if h j+1,j = 0 then 11: break 12: end if 13: v j+1 = w j /h j+1,j 14: end for

18 18/115 注记 基于 MGS 的 Arnoldi 过程与基于 CGS 的 Arnoldi 过程是等价的, 但 MGS 更稳定. 已有学者证明修正 Gram-Schmidt 过程是向后稳定的. 注记 但在某些极端情况下, 算法 2.2 仍不够可靠, 此时可以再做一次 MGS, 或者采用 Householder 变换.

19 19/115 3 一般线性方程组的 Krylov 子空间方法 给定一个迭代初始值 x (0) R n, 令 K m = span{r 0, Ar 0, A 2 r 0,..., A m 1 r 0 }, 其中 r 0 = b Ax (0) 为初始残量. 下面介绍分别基于 L = K m 和 L = AK m 的投影方法 : 完全正交方法 : L = K m 广义极小残量法 : L = AK m

20 20/ 完全正交方法 (FOM): L = K m 寻找 x x (0) + K m 满足 b A x K m 由 Arnoldi 过程可知 VmAV m = H m 且 Vmr 0 = Vm(βv 1 ) = βe 1, 其中 β = r 0 2, e 1 = [1, 0,..., 0] R m. 由 x x (0) + K m 可知, 存在向量 ỹ R m 使得 x = x (0) + V m ỹ. 由正交性条件 b A x K m 可得 Vm(b A x) = 0, 矩阵 V m 的列向量构成 K m 的一组基 与 K m 正交等价于与 K m 的一组基正交

21 21/115 即 0 = V m(b Ax (0) AV m ỹ) = V mr 0 V mav m ỹ = βe 1 H m ỹ. 如果矩阵 H m 非奇异, 则 ỹ = βh 1 m e 1, 所以 x = x (0) + βv m H 1 m e 1

22 22/115 算法 3.1 完全正交方法 (FOM) 1: 给定初值 x (0), 计算 r 0 = b Ax (0), β = r 0 2 和 v 1 = r 0 /β 2: for j = 1, 2,..., m 1 do 3: w j = Av j 4: for i = 1, 2,..., j do 5: h ij = (w j, v i ) 6: w j = w j h ij v i 7: end for 8: h j+1,j = w j 2 9: if h j+1,j = 0 then 10: m = j, break 11: end if 12: v j+1 = w j /h j+1,j 13: end for 14: 求解线性方程组 H m ỹ = βe 1 并计算近似解 x = x (0) + V m ỹ

23 23/115 注记 由于子空间 K m 的维数远远小于 n, 因此方程组 H m ỹ = βe 1 可 以通过直接法求解, 如列主元 Gauss 消去法. 事实上, 由于 H m 是一个上 Hessenberg 矩阵, 我们也可以采用基于 Givens 变换的 QR 分解来求解. 这样更稳定. 注记 FOM 方法的一个难点是如何选取 m. 如果 m 太小, 则近似解 x 可能达不到精度要求, 如果太大, 可能会多做一些无用功. 因此在实际计算时, 我们采用动态的方法, 即不断增长 m 的值. 当残量满足精度要求时就停止迭代. 因此, 在每一次迭代后, 我们都需要估计残量的范数.

24 24/115 残量估计 定理 5 设 x x (0) + K m 由算法 3.1 得到的近似解, 则残量表达式为 r b A x = h m+1,m (e mỹ) v m+1 其中 ỹ = βh 1 m e 1, e m = [0,..., 0, 1] R m. 因此, r 2 = h m+1,m e mỹ = h m+1,m ỹ(m) ( 板书 ) AV m = V m H m + h m+1,m v m+1 e m

25 25/115 注记 由定理 5 可知, 我们可以直接计算出残量 r 的大小, 而无需在每次迭代中计算近似解 x. 当残量满足精度要求时, 我们再计算近似解. 注记 在不考虑舍入误差的情况下, 当 m = n 时, 必有 r 2 = 0. 注记 如果 Arnoldi 过程提前终止, 即存在整数 k (k < n) 使得 h k+1,k = 0, 则由定理 5 可知, 此时 r = 0, 因此近似解 x 就是方程组的精确解.

26 26/115 算法 3.2 实用 FOM 方法 1: 给定初值 x (0) 和 ( 相对 ) 精度要求 ε > 0 2: 计算 r 0 = b Ax (0), β = r 0 2 和 v 1 = r 0 /β 3: for j = 1, 2,... do 4: w j = Av j 5: for i = 1, 2,..., j do 6: h ij = (w j, v i ) 7: w j = w j h ij v i 8: end for 9: h j+1,j = w j 2 10: 求解方程 H j ỹ = βe 1 11: if h j+1,j ỹ(j) /β < ε then % relative residual 12: break 13: end if 14: v j+1 = w j /h j+1,j 15: end for 16: 计算近似解 x = x (0) + V j ỹ

27 27/ 广义极小残量法 (GMRES): L = AK 广义极小残量法 (Generalized Minimum Residual) 可描述为 寻找 x x (0) + K 满足 b A x AK. (4.7) GMRES 是当前求解大规模非对称线性方程组的首选方法 定理 6 设 A R n n, K = K(A, r 0 ), 则 x 是问题 (4.7) 的解当且仅当 x 是下面的极小化问题的解 min b Ax 2, (4.8) x x (0) +K 即在仿射空间 x (0) + K 中, x 所对应的残量范数最小. ( 板书 ) 与 FOM 的推导过程不同, 我们需要用另外一种方式来推导 GMRES 这个最优性质是 GMRES 方法名称的由来

28 28/115 GMRES 的推导过程任意向量 x x (0) + K m 都可表示为 x = x (0) + V m y, 其中 y R m, 故 b Ax = b A(x (0) + V m y) = r 0 AV m y = βv 1 V m+1 H m+1,m y = V m+1 (βe 1 H m+1,m y), (4.9) 其中 e 1 = [1, 0,..., 0] R m+1. 又 V m+1 的列向量标准正交的, 所以 b Ax 2 = V m+1 (βe 1 H m+1,m y) 2 = βe 1 H m+1,m y 2. (4.10) 因此, 极小化问题 (4.8) 的解可表示为 x = x (0) + V m ỹ, (4.11)

29 29/115 其中 ỹ 是下面最小二乘问题的解 min y R m βe 1 H m+1,m y 2. (4.12) 注记 与 FOM 方法中求解一个线性方程组不同的是, 在 GMRES 方法中, 我们是求解一个最小二乘问题.

30 30/115 最小二乘求解 一般 m 不会很大, 因此可以使用 QR 方法来求解 (4.12). 设 H m+1,m = Q m+1r m+1,m 是 H m+1,m 的 QR 分解, 其中 Q m+1 R (m+1) (m+1) 是正交矩阵, R m+1,m R (m+1) m 是上三角矩阵. 于是 βe 1 H m+1,m y 2 = βe 1 Q m+1r m+1,m y 2 [ ] = βq m+1 e 1 R m+1,m y 2 = βq R m 1 y, 0 2 其中 q 1 是 Q m+1 的第一列, R m = R m+1,m (1 : m, :) R m m 表示 R m+1,m 的前 m 行. 所以 y 可以通过求解下面的上三角方程组获得 R m ỹ = βq 1 (1 : m), (4.13) 其中 q 1 (1 : m) 表示 q 1 的前 m 个元素组成的向量.

31 31/115 残量计算 与 FOM 方法类似, 我们也可以直接得到残量范数的表达式. 定理 7 设 x x (0) + K m 是由 GMRES 方法所得到的近似解, 则 r 2 = b A x 2 = β q 1 (m + 1) 其中 q 1 (m + 1) 表示 q 1 R m+1 的最后一个分量. r 2 = b A x 2 = βe 1 H m+1,m ỹ 2 [ ] = βq R m 1 ỹ 0 = β q 1 (m + 1). 2

32 32/115 算法 3.3 广义极小残量法 (GMRES) 1: 给定初值 x (0) 和 ( 相对 ) 精度要求 ε > 0 2: 计算 r 0 = b Ax (0) 和 β = r 0 2 3: v 1 = r 0 /β 4: for j = 1, 2,... do 5: w j = Av j 6: for i = 1, 2,..., j do % Arnoldi process 7: h ij = (w j, v i ) 8: w j = w j h ij v i 9: end for 10: h j+1,j = w j 2 11: if h j+1,j = 0 then 12: m = j, break 13: end if

33 33/115 14: v j+1 = w j /h j+1,j 15: if r 2 /β < ε then % check convergence 16: m = j, break 17: end if 18: end for 19: 求解上三角线性方程组 (4.13) 20: 计算近似解 x = x (0) + V m ỹ

34 34/115 算法的中断 在 GMRES 算法中, 如果第 k (k < n) 步时 h k+1,k = 0, 则提前终止, 且 AV k = V k H k, ỹ = H 1 k (βe 1). min βe 1 H k+1,k y 2 H k+1,k 最后一行全为 0 因此 r 2 = b A x 2 = b Ax (0) AV k ỹ 2 = r 0 V k H k ỹ 2 = βv 1 V k (βe 1 ) 2 = 0. 这意味着 x 就是原方程组的精确解. 因此, 这种提前终止是好事.

35 35/115 反之, 如果在第 k 步时有 r = b A x = 0, 则必有 h k+1,k = 0. 定理 8 设 A R n n 非奇异, 则 GMRES 方法在第 k 步提前终止的充要条件是近似解即为精确解. 证明. 留作练习.

36 36/ GMRES 方法的实施细节 我们注意到, 在 GMRES 方法中, 计算残量的范数时需要计算 H j+1,j 的 QR 分解, 工作量太大! 在实际计算中, 我们采用递推方式和 Givens 变换, 在 H j,j 1 的 QR 分解基础上, 通过一次 Givens 变换得到 H j+1,j 的 QR 分解, 从而大大降低运算量. 下面我们详细描述具体的递推计算过程. 第一步 [ 当 j = 1 时, 只需对 H 2,1 = h 11 h 21 ] 做一次 Givens 变换即可得 QR 分解

37 第二步 假定 H j,j 1 R j (j 1) 的 QR 分解为 H j,j 1 = (G j 1 G j 2 G 1 ) R j,j 1 = Q j [ R j 1 0 ] j (j 1) 其中 R j 1 R (j 1) (j 1) 是上三角矩阵, Q j = G 1 G 2 G j 1. 这里 G i 都是 Givens 变换. 为了确保矩阵乘积的相容性, 我们假定 G i 的维数会根据需要自动增大, 即在右下角添加单位阵. 第三步 考虑 H j+1,j 的 QR 分解. 易知 [ ] H j,j 1 h j H j+1,j =, 0 h j+1,j 其中 h j = [h 1,j, h 2,j,..., h j,j ] 是由 H j+1,j 最后一列的前 j 个元素构 37/115

38 38/115 成的向量. 于是 [ ] [ ] [ Q j 0 Q j 0 H j+1,j = = H j,j 1 h j 0 h j+1,j ] R j 1 hj 1 Q j h j = 0 ĥ j,j 0 h j+1,j 0 h j+1,j [ R j,j 1 其中 hj 1 是 Q j h j 的前 j 1 个元素构成的向量, ĥj,j 是 Q j h j 的最后一个元素. 下面我们利用 Givens 变换, 将 h j+1,j 化为 0. 计算 ], 构造 Givens 变换 h j,j = ĥ2 j,j + h2 j+1,j, c j = ĥj,j h j,j, s j = h j+1,j h j,j.

39 39/115 I 0 0 G j = 0 c j s j R (j+1) (j+1). 0 s j c j 令 则 Q j+1 H j+1,j = G j [ ] Q j 0 Q j+1 = G j R (j+1) (j+1), 0 1 R j 1 hj 1 0 ĥ j,j 0 h j+1,j = R j 1 hj 1 0 hj,j 0 0 R j+1,j. 所以 H j+1,j = Q j+1 R j+1,j 就是 H j+1,j 的 QR 分解.

40 40/115 注记 在上述过程中, 我们只需将 Givens 变换 G 1, G 2,..., G j 依次作用到 H j+1,j 最后一列上, 就可以得到 R j 的最后一列. 当残量充分小时, 我们再通过解方程 (4.13) 得到 ỹ. 由 Q m+1 = G m G m 1 G 1 可知, Q m+1 (βe 1 ) 可通过将 G 1, G 2,..., G m 依次作用到向量 βe 1 上得到. 在 GMRES 的具体实现过程. 为节省存储, 我们将 R m 存在 H m+1,m 中.

41 41/115 算法 3.4 实用的 GMRES 方法 1: 给定初值 x (0) 和 ( 相对 ) 精度要求 ε > 0 2: 计算 r 0 = b Ax (0) 和 β = r 0 2 3: v 1 = r 0 /β, ξ = βe 1 4: for j = 1, 2,... do 5: w j = Av j 6: for i = 1, 2,..., j do % Arnoldi process 7: h ij = (w j, v i ) 8: w j = w j h ij v i 9: end for 10: h j+1,j = w j 2 11: for i = 1, 2,..., j 1 do % Apply G j 1,..., G 1 to the last column of H j+1,j

42 42/115 12: [ 13: end for h i,j h i+1,j ] = [ 14: if h j+1,j = 0 then c i s i 15: set m = j and break 16: end if 17: v j+1 = w j /h j+1,j s i c i ] [ h i,j h i+1,j 18: if h j,j > h j+1,j then % Form the Givens rotation G j 19: c j = 1 1+τ 2, s j = c j τ where τ = hj+1,j h j,j 20: else 21: s j = 1 1+τ 2, c j = s j τ where τ = 22: end if ] hj,j h j+1,j 23: h j,j = c j h j,j + s j h j+1,j, h j+1,j = 0 % Apply G j to last column of H j+1,j

43 43/115 24: [ ξ j ξ j+1 ] = [ c j s j s j c j ] [ ] ξ j 0 % Apply G j to the right-hand side 25: if ξ j+1 /β < ε then % Check convergence: r 2 = ξ j+1 26: set m = j and break 27: end if 28: end for 29: 计算 ỹ = R 1 m ξ (m), 其中 R m = H(1 : m, 1 : m), ξ (m) = ξ(1 : m) 30: 计算近似解 x = x (0) + V m ỹ 注记 在 FOM 方法和 GMRES 方法中, 都是采用 MGS 来构造 Krylov 子空间的标准正交基, 在某些特殊情况下, 可能需要使用 Householder 变换来增加方法的稳定性.

44 44/115 注记 在不考虑舍入误差的情况下, 当 GMRES 方法迭代到 n 步时, 残量肯定为 0, 因此方法肯定收敛. 注记 对于大规模的线性方程组, 我们不可能迭代 n 步, 因为 GMRES 每一步迭代的运算量和存储量都会随着迭代步数的增加而增加, 所以当迭代步数增大时, 运算时间和存储量会大大增加, 这使得方法失去竞争性.

45 45/115 带重启的 GMRES 我们希望将迭代步数控制在一个能接受的范围内. 目前实现这一目标的通用做法是采用重启策略, 即给定一个最大迭代步数 m ( 通常远远小于 n, 如 20, 50 等 ) 如果 GMRES 迭代到 m 步后仍不收敛, 则计算出近似解 x (m) x (0) + K m. 然后将其作为新的迭代初值, 即令 x (0) = x (m), 重新启动 GMRES 方法. 该过程可以重复执行, 直到找到满意的近似解为止.

46 46/115 算法 3.5 GMRES(m): 带重启的 GMRES 方法 1: 给定正整数 m n, 初值 x (0) 和 ( 相对 ) 精度要求 ε > 0 2: 计算 r 0 = b Ax (0) 和 β = r 0 2 3: v 1 = r 0 /β 4: for j = 1, 2,..., m do 5: 执行 GMRES 方法 3.4 的第 5 步至第 27 步 6: end for 7: 计算 ỹ = R 1 m ξ (m), 其中 R m = H(1 : m, 1 : m), ξ (m) = ξ(1 : m) 8: 计算近似解 x = x (0) + V m ỹ 9: if ξ m+1 /β < ε then 10: stop 11: end if 12: 令 x (0) = x (m), 返回第 2 步

47 47/115 注记 加入重启技术可以节省运算量和存储量, 但随之也会带来一些负面效应. 如收敛速度会减慢, 这是因为在重启方法时, 一些有用信息会被丢弃. 如果 A 不是正定矩阵, 则可能还会出现方法停滞不前的情形, 即残量范数很难减小. 此时即使总迭代步数已经达到 n 步, 或者甚至超过 n 步, 但方法仍然不收敛. 另外, 怎样选取合适的重启步数 m, 目前仍然没有很好的办法, 经验很重要!

48 48/ 举例 由最优性可知, GMRES 残量范数 r k 2 递减, 但不一定严格递减. 在最坏的情况下, 可能会出现 r k 2 保持不变, 直到最后一步才降为 0. 例 1 考虑线性方程组 Ax = b, 其中 A = , b = a 1 a 2 a n 1 a n 具体求解过程参见 MATLAB 程序 GMRES_Example01.m. 下图绘出了 n = 64 时的 ( 相对 ) 残量范数下降曲线.

49 49/115

50 50/115 第二个例子是关于偏微分方程的. 例 2 考虑下面的带齐次 Dirichlet 边界条件的偏微分方程 : u(x, y) + u x (x, y) + u y (x, y) + u(x, y) = f(x, y), (4.14) (x, y) Ω [0, 1] [0, 1] u(x, 0) = u(x, 1) = u(0, y) = u(1, y) = 0, 0 < x, y < 1. 二阶导数用五点差分离散, 一阶导数用中心差分离散, 步长 h x = h y = 1/(N + 1). 可得具体离散格式为 4u i,j u i+1,j u i 1,j u i,j+1 u i,j 1 h 2 + u i+1,j u i 1,j 2h + u i,j+1 u i,j 1 2h + u i,j = f i,j

51 51/115 对应的离散线性方程组为 [I T + T I + h(i D) + h(d I) + h 2 I]u = h 2 f. 其中 T = tridiag( 1, 2, 1), D = tridiag( 1/2, 0, 1/2). 在实际测试中, 我们假设解析解为 u(x, y) = xy(1 x)(1 y). 由此可以求出右端项 f(x, y) = (3 2x)(1 y)y + (3 2y)(1 x)x + x(1 x)y(1 y). 用 GMRES 方法求解. 最大迭代步数设为 IterMax = N 2, 迭代终止条

52 52/115 件 r k 2 r 0 2 tol 我们画出 n = 32 时的相对残量下降曲线. MATLAB 程序见 GMRES_Example02_PDE.m

53 53/115

54 54/115 4 对称线性方程的 Krylov 子空间方法 4.1 Lanczos 过程 如果 A 是对称的, 则由 Arnoldi 过程的性质可知 H m = V mav m, 因此 H m 也对称. 又 H m 是上 Hessenberg 矩阵, 所以 H m 是一个对称三对角矩阵. 为了书写方便, 我们记 α j = h j,j, β j = h j,j+1, 并将 H m 改写为 T k, 即 α 1 β 1. β T k = βk 1 β k 1 α k

55 55/115 与 Arnoldi 过程类似, 我们有下面的性质 AV k = V k T k + β k v k+1 e k, V k AV k = T k. (4.15) (4.16) 考察关系式 (4.15) 两边的第 j 列可知 β j v j+1 = Av j α j v j β j 1 v j 1, j = 1, 2,..., 这个三项递推公式是 CG 方法的核心, 它使得每一步迭代的运算量与迭代步数无关. 其中 v 0 = 0 和 β 0 = 0. 基于这个三项递推公式, Arnoldi 过程可简化为 Lanczos 过程.

56 56/115 算法 4.1 Lanczos 过程 1: 给定非零向量 r, 令 v 0 = 0, β 0 = 0 2: v 1 = r/ r 2 3: for j = 1, 2,... do 4: w j = Av j β j 1 v j 1 5: α j = (w j, v j ) 6: w j = w j α j v j 7: β j = w j 2 8: if β j = 0 then 9: break 10: end if 11: v j+1 = w j /β j 12: end for

57 57/115 如果不考虑舍入误差, 则 Lanczos 过程生成的向量是相互正交的. 定理 9 设 v 1, v 2,..., v k,... 由 Lanczos 方法 4.1 生成, 则 1, i = j, (v i, v j ) = δ ij 0, i j, for i, j = 1, 2,.... (4.17) ( 板书 ) 注记 该结论说明, 只需正交化相邻的三个向量, 因此只需存三个向量, 而且每次迭代的运算量也固定不变. 这就是 Lanczos 过程的优点. 实际计算中, 由于舍入误差影响, Lanczos 过程生成的向量可能会逐渐失去正交性. 这时需采取一些补救措施, 如全部重新正交化或部分重新正交化.

58 58/ 共轭梯度法 (CG) 首先考虑 A R n n 是对称正定情形. 共轭梯度法 (Conjugate Gradient) 本质上与 FOM 是等价的, 即都是取 L = K. 但由于 A 对称, 因此可借助三项递推公式来简化运算. CG 方法是当前求解对称正定线性方程组的首选方法. 方法描述 设 x (0) 是初始值, 在仿射空间 x (0) + K k 中寻找近似解 x (k), 满足 x (k) x (0) + K k 且 b Ax (k) K k. (4.18) CG 方法就是根据这个性质来推导的.

59 方法推导 首先, 与 FOM 方法类似, 我们可得 x (k) = x (0) + V k y k, y k = T 1 k (βe(k) 1 ), (4.19) 其中 V k = [v 1, v 2,..., v k ], β = r 0 2, e (k) 1 = [1, 0,..., 0] R k. 因此需要求解一个以 T k 为系数矩阵的线性方程组. 由于 A 对称正定, 故 T k 也对称正定. 所以可以用 Cholesky 分解或 LDL 分解来求解. 设 T k 的 LDL 分解为 T k = L k D k L k. 于是可得 x (k) = x (0) +V k y k = x (0) +V k T 1 k (βe(k) 1 ) = x (0) +(V k L k 如果 x (k) 满足精度要求, 则计算结束. 否则我们需要计算 )(βd 1 k L 1 k e(k) 1 ). x (k+1) = x (0) +V k+1 T 1 k+1 (βe(k+1) 1 ) = x (0) +(V k+1 L k+1 )(βd 1 k+1 L 1 k+1 e(k+1) 1 ), 59/115

60 60/115 这里假定 T k+1 的 LDL 分解为 T k+1 = L k+1 D k+1 L k+1. 下面就考虑如何利用递推方式来计算 x (k), k = 1, 2,.... 记 P k V k L k = [ p 1, p 2,..., p k ] R n k, ỹ k βd 1 k L 1 k e(k) 1 = [η 1,..., η k ] R k. 我们首先证明下面的结论. 引理 2 下面的递推公式成立 : P k+1 V k+1 L k+1 = [ P k, p k+1 ], y k+1 βd 1 k+1 L 1 k+1 e(k+1) 1 = [ỹ k, η k+1], k = 1, 2,.... ( 板书 )

61 61/115 有了上面的性质, 我们就可以得到从 x (k) 到 x (k+1) 的递推公式 x (k+1) = P k+1 ỹ k+1 = [ P k, p k+1 ] [ ỹ k η k+1 ] = x (k) + η k+1 p k+1. 为了判断近似解的精度, 我们还需要计算残量. 通过直接计算可知 (4.20) r k+1 = b Ax (k+1) = b A(x (k) + η k+1 p k+1 ) = r k η k+1 A p k+1. (4.21)

62 62/115 另一方面, 我们有 r k = b Ax (k) = b A(x (0) + V k y k ) = r 0 AV k y k = βv 1 V k T k y k β k v k+1 e k y k = β k (e k y k)v k+1, 即 r k 与 v k+1 平行. 记 r k = τ k v k+1, k = 0, 1, 2,..., (4.22) 其中 τ 0 = β = r 0 2, τ k = β k (e k y k), k = 1, 2,.... 定义 p k = τ k 1 p k, k = 1, 2,....

63 63/115 所以 {p k } 满足下面的递推公式 : p k+1 = τ k p k+1 = τ k (v k+1 l k p k ) = r k + µ k p k, (4.23) p k+1 = l k p k + v k+1 r k = τ k v k+1 其中 µ k = l kτ k τ k 1, k = 1, 2,.... 于是我们就得到下面的递推公式 r k+1 = r k η k+1 A p k+1 = r k ξ k+1 Ap k+1 x (k+1) = x (k) + η k+1 p k+1 = x (k) + ξ k+1 p k+1, (4.24) (4.25) 其中 ξ k+1 = η k+1 τ k, k = 0, 1, 2,....

64 64/115 系数 ξ k+1 和 µ k 的计算方法 引理 3 下面的结论成立 : (1) r 1, r 2,..., r k,... 相互正交 ; (2) p 1, p 2,..., p k,... 相互 A 共轭 ( 或 A 正交 ), 即当 i j 时有 p i Ap j = 0. ( 板书 ) p k+1 = τ k p k+1 = τ k (v k+1 l k p k ) = r k + µ k p k, (4.26)

65 65/115 r k+1 = r k η k+1 A p k+1 = r k ξ k+1 Ap k+1 x (k+1) = x (k) + η k+1 p k+1 = x (k) + ξ k+1 p k+1, (4.27) (4.28) 在等式 (4.26) 两边同时左乘 p k+1 A 可得 p k+1 Ap k+1 = p k+1 Ar k + µ k p k+1 Ap k = r k Ap k+1. 再用 r k 左乘方程 (4.27) 可得 0 = r k r k+1 = r k r k ξ k+1 r k Ap k+1, 故 ξ k+1 = r k r k r k Ap k+1 = r k r k p k+1 Ap. (4.29) k+1

66 66/115 等式 (4.26) 两边同时左乘 p k A 可得 故 0 = p k Ap k+1 = p k Ar k + µ k p k Ap k, µ k = p k Ar k p k Ap k = r k Ap k p k Ap. (4.30) k 为了进一步减少运算量, 我们将上式进行改写. 用 r k+1 左乘方程 (4.27) 可得 r k+1 r k+1 = r k+1 r k ξ k+1 r k+1 Ap k+1 = ξ k+1 r k+1 Ap k+1, 故 ξ k+1 = r k+1 r k+1 r k+1 Ap. k+1

67 所以 ξ k = r k r k r k Ap, 即 r k Ap k = r k r k/ξ k, 代入 (4.30) 可得 k µ k = r k Ap k p k Ap k = r k r k p k Ap 1 = k ξ k r k r k p k Ap k p k Ap k r k 1 r k 1 = r k r k r k 1 r. k 1 (4.31) 综合公式 (4.26), (4.27), (4.28) 和 (4.29), (4.31) 即可得 CG 的迭代格式 : p k+1 = r k + µ k p k, 其中 µ k = r k r k r k 1 r, k 1 x (k+1) = x (k) + ξ k+1 p k+1, r k+1 = r k ξ k+1 Ap k+1, 其中 ξ k+1 = r k r k p k+1 Ap. k+1 注意, 以上递推公式是从 k = 1 开始. k = 0 时的公式需要另外推导. 67/115

68 k = 0 时的计算公式 首先, 由 p 1 的定义可知 p 1 = P 1 = V 1 L 1 = v 1. 因此 p 1 = τ 0 p 1 = βv 1 = r 0. 其次, 由 Lanczos 过程可知 T 1 = α 1 = v 1 Av 1. 注意到 β 2 = r 0 r 0, 于是 x (1) = x (0) + V 1 T 1 1 (βe 1 ) = x (0) + β v 1 Av 1 v 1 = x (0) + r 0 r 0 p 1Ap 1 p 1. 令 ξ 1 = r 0 r 0 p 1Ap 1, 则当 k = 0 时关于 x (k+1) 的递推公式仍然成立. ( 注 : 之前的 ξ k+1 计算公式 (4.29) 只对 k 1 有定义 ) 68/115

69 69/115 最后考虑残量. 易知 r 1 = b Ax (1) = b Ax (0) r 0 r 0 p 1Ap 1 Ap 1 = r 0 ξ 1 Ap 1, 即当 k = 0 时关于 r k+1 的递推公式也成立. 综合公式 (4.26), (4.27), (4.28) 和 (4.29), (4.31) 即可得 CG 方法 : 算法 4.2 共轭梯度法 (CG) 1: 给定初值 x (0), ( 相对 ) 精度要求 ε > 0 和最大迭代步数 IterMax 2: 计算 r 0 = b Ax (0) 和 β = r 0 2 3: for k = 1 to IterMax do 4: ρ = r k 1 r k 1 5: if k > 1 then

70 70/115 6: µ k 1 = ρ/ρ 0 7: p k = r k 1 + µ k 1 p k 1 8: else 9: p k = r 0 10: end if 11: q k = Ap k 12: ξ k = ρ/(p k q k) 13: x (k) = x (k 1) + ξ k p k 14: r k = r k 1 ξ k q k 15: relres = r k 2 /β 16: if relres< ε then 17: break 18: end if 19: ρ 0 = ρ 20: end for

71 71/115 21: if relres< ε then 22: 输出近似解 x (k) 及相关信息 23: else 24: 输出算法失败信息 25: end if 注记 编程练习 :CG 方法求解二维 poisson 方程. 在优化领域, CG 方法中的向量 p k 称为搜索方向, 而 ξ k 称为步长. 因此, 搜索方法是 A 共轭的, 这就是方法名称的由来. 注记 在 CG 方法中, 需要额外存储四个向量 x, p, Ap 和 r, 每步迭代的主要运算为一个矩阵向量乘积和两个内积.

72 72/115 CG 方法也具有最优性质. 定理 10 设 A 对称正定, 则 x (k) = arg min x x A (4.32) x x (0) +K k 当且仅当 与 GMRES 方法不同的是, CG 方法极小化的是绝对误差的 A 范数. x (k) x (0) + K k 且 b Ax (k) K k. (4.33) ( 板书 ) 注记 如果 A 对称, 但不定, 此时 A 范数就没有定义, 因此最优性结论不再成立. 而从方法推导过程可知, 此时仍可以使用 CG 方法, 但由于 LDL 分解不一定存在, 因此方法可能会中断.

73 73/115 例 3 用 CG 方法求解线性方程组 Ax = b, 其中 2 1. A = B I+I B R N 2 N , B = 画出 n = 64 时的相对残量下降曲线 : Poisson_CG.m N N 1 1, b =. 1 N 2 1.

74 74/115

75 To be continued... 75/115

76 76/ 极小残量法 (MINRES) 这里假定 A 对称, 但不定. 此时 CG 方法不再适用. 与 GMRES 方法类似, 我们取约束空间为 L k x (0) + K k 中寻找近似解 x (k), 满足 = AK k, 即在仿射空间 x (k) x (0) + K k 且 b Ax (k) AK k. 这就是极小残量法 (MINRES ). 它也可以看作是 GMRES 的对称情形, 即作用在对称方程组上的 GMRES 方法. 但由于 A 对称, 因此在计算 K k 的正交基时, 可以采用 Lanczos 方法, 即三项递推方法, 从而节省运算量. 与 GMRES 方法类似, MINRES 具有下面的最优性质.

77 77/115 定理 11 设 A R n n 是对称矩阵, x (k) x (0) + K k 是 MINRES 方法迭代 k 步后得到的近似解, 则 x (k) = arg min b Ax 2, x x (0) +K k 即 x (k) 是残量在仿射空间 x (0) + K k 中的唯一最小值点. 由 GMRES 方法的推导过程可知, 当 MINRES 方法迭代 k 步之后, 近似解可表示为 其中 x (k) = x (0) + V k y k (4.34) y k = arg min y R k βe 1 T k+1,k y 2 (4.35)

78 78/115 这里 T k+1,k 是由 Lanczos 方法生成的 (k + 1) k 的上 Hessenberg 矩阵 α 1 β 1. β T k+1,k = βk 1. β k 1 与 CG 方法类似, 我们无需存储所有的基向量 v i s. 事实上, x (k) 可以通过直接更新 x (k 1) 得到. 下面我们就给出推导过程. α k β k 最小二乘问题的求解 我们仍然用 QR 分解来求解最小二乘问题 (4.35). 设 T k+1,k = Q k+1 R k+1,k

79 79/115 是 T k+1,k 的 QR 分解, 其中 Q k+1 R (k+1) (k+1) 是正交矩阵, R k+1,k R (k+1) k 是一个上三角矩阵. 由于 T k+1,k 是三对角的, 所以 R k+1,k 也只有三条对角线非零, 即 R k+1,k = τ (1) 1 τ (2) 1 τ (3) 1 τ (1) 2 τ (2) 2 τ (3) τ (1) k 2 τ (2) k 2 τ (3) k 2 τ (1) k 1 τ (2) k 1 τ (1) k 0 0 (k+1) k [ R k 0 ],

80 80/115 其中 R k 是由 R k+1,k 的前 k 行组成的 k 阶矩阵. 于是, 我们有 βe 1 T k+1,k y 2 = βe 1 Q k+1 R k+1,k y) 2 ( [ ] = Q k+1 Q k+1 (βe R k 1) y) 0 [ ] 2 = [Q k+1,k, q k+1 ] R k y (βe 1 ), (4.36) 0 2 其中 Q k+1,k 是由 Q k+1 的前 k 列组成的子矩阵, q k+1 是 Q k+1 的最后一列. 如果 R k 非奇异, 则最小二乘问题 (4.35) 的解为 因此, y k = βr 1 k Q k+1,k e 1. x (k) = x (0) + βv k R 1 k Q k+1,k e 1. (4.37)

81 81/115 T k+1,k 的 QR 分解的具体实现 首先, 由于 T k+1,k 是上 Hessenberg 矩阵, 因此可以采用 Givens 变换来做. 其次, 与 GMRES 类似, 我们可以借助递推方法来减少运算量, 即在 T k,k 1 的 QR 分解的基础上, 通过一次 Givens 变换得到 T k+1,k 的 QR 分解. 假定 T k,k 1 的 QR 分解为 T k,k 1 = (G k 1 G k 2 G 1 ) R k,k 1 = Q k R k,k 1, 其中 G i 表示 Givens 变换, Q k = (G k 1 G k 2 G 1 ), R k,k 1 是上三角

82 82/115 阵 : R k,k 1 = τ (1) 1 τ (2) 1 τ (3) 1 τ (1) 2 τ (2) 2 τ (3) (3) τ. k 3... (2) τ k 2 τ (1) k 为了保证矩阵乘积的相容性, 我们假定 G i 的维数是自动增长的, 即 I i 1 c i s i G i = s i c i Rk k, i = 1, 2,... k 1. I k i 1

83 83/115 对于给定矩阵 B, 变换 B G i B 仅仅修改 B 的第 i 和第 (i + 1) 行. 因此, T k+1,k = 且 Q 1 k β k 1 α k T k,k β k 1 α k 0 β k [ ] Q k 0 = 0 1 = G k 1 G k 2 G β k 1 α k R k,k 1 Q 1 k 0 β k = G k 1 G k β k 1 α k 0. 0 β k 1 α k [ 0. 0 = τ (3). k 2 τ (2) k 1 α k Q k ] T k+1,k,

84 84/115 构造 Givens 变换 G k = k 1 I c k s k s k c k R (k+1) (k+1), 消去 Tk+1,k 最后一列的最后一个元素 β k, 即选取 c k = α k τ (1) k, s k = β k τ (1) k 其中 τ (1) k = α k 2 + β2 k. 由于 R k,k 1 的最后一行全为 0, 且左乘 G k 只会影响 Tk+1,k 的最后两行, 因此

85 G k Tk+1,k = 0. 0 τ (3) k 2 τ (2) k 1 τ (1) k 0 0 R k,k 1 = τ (1) 1 τ (2) 1 τ (3) (3) τ k 3... (2) τ k 2 τ (3) k 2 τ (1) k 1 τ (2) k 1 τ (1) k 0 0 R k+1,k. 于是我们就得到 T k+1,k 的 QR 分解 T k+1,k = Q k+1 R k+1,k, 其中 ( 这里我们假定 Givens 变换 G i 的维数是自动增长的 ) [ ] Q k 0 Q k+1 = G k 0 1 = (G kg k 1... G 1 ). 需要指出的是, R k,k 1 是 R k+1,k 的前 k 行和前 k 列. 因此, R k 1 是 R k 85/115

86 86/115 的 (k 1) 阶顺序主子矩阵. 下面我们给出从 x (k) 到 x (k+1) 的递推公式. 设 通过直接计算, 由 P k R k = V k 可知 p 1 = v 1 /τ (1) 1, ( ) p 2 = v 2 τ (2) 1 p 1 p i = P k = [p 1, p 2,..., p k ] V k R 1 k. /τ (1) ( v i τ (2) i 1 p i 1 τ (3) i 2 p i 2 2, (4.38) ) /τ (1) i, i = 3, 4,.... 又 V k = [V k 1, v k ], 且 R k 1 是 R k 的 (k 1) 阶顺序主子矩阵. 所以 P k = [P k 1, p k ].

87 87/115 记 βq k+1 e 1 的前 k 个元素组成的向量为 ξ (k), 即 [ ] βq k+1 e ξ (k) 1 =, ξ k+1 其中 ξk+1 表示 βq k+1 e 1 的最后一个元素. 易知 Q k+1,k e 1 和 Q k+1 e 1 的前 k 个分量是一样的, 即 ξ (k) = βq k+1,k e 1. 因此 y (k) = βr 1 k Q k+1,k e 1 = R 1 k ξ(k). (4.39) 又 [ ] Q k 0 Q k+1 = G k 0 1, 所以, Q k+1 (1, 1 : k 1) = Q k (1, 1 : k 1), 即 Q k+1 e 1 和 Q k e 1 的前 k 1 个分量是一样的. 于是我们可以得到 ξ (k) 和 ξ (k 1) 之间的关系

88 88/115 式 : ξ (k) = [ ξ (k 1) ξ k ], k = 2, 3,..., 其中 ξ k 表示 ξ (k) 的最后一个分量. 因此, 对所有正整数 k, ξ (k) 可以表示为 ξ (k) = [ξ 1, ξ 2,..., ξ k ], k = 1, 2,....

89 89/115 根据 (4.37), 我们有 x (k) = x (0) + V k R 1 k Q k+1,k (βe 1) = x (0) + P k ξ (k) = x (0) + [P k 1, p k ] [ ξ (k 1) ξ k = x (0) + P k 1 ξ (k 1) + ξ k p k ] = x (k 1) + ξ k p k, (4.40) 其中 p k 可通过递推公式 (4.38) 来计算, 而 ξ k 是 βq k+1 e 1 的第 k 个分量. 又 βq k+1 e 1 = Q k+1 (βe 1) = G k Gk 1 G 1 (βe 1 ), 即向量 βq k+1 e 1 可以通过将 Givens 变换依次作用在 βe 1 上得到.

90 下面考虑残量的计算. 由 (4.39) 可知 r k 2 = b Ax (k) 2 = r 0 AV k y k 2 = βv 1 V k+1 T k+1,k y k 2 = V k+1 (βe 1 Q k+1 R k+1,k y k ) 2 ( [ ] ) = Q k+1 Q k+1 (βe R k 1) y k 0 [ ] [ ] 2 ξ (k) R k y k = ξ k+1 0 = ξ k+1, 其中 ξk+1 是 Q k+1 (βe 1) = G k G k 1 G 1 (βe 1 ) 的最后一个分量. 综上所述, 根据 p k 的递推公式 (4.38) 和 x (k) 的递推公式 (4.40), 我们就可以构造下面的 MINRES 方法. 2 90/115

91 91/115 算法 4.3 Minimal Residual (MINRES) Method 1: 给定初值 x (0), ( 相对 ) 精度要求 ε > 0 和最大迭代步数 IterMax 2: 计算 r 0 = b Ax (0) 和 β = r 0 2 3: 令 β 0 = 0, v 0 = 0, p 1 = 0, p 0 = 0 4: v 1 = r 0 /β, ξ = βe 1 5: for k = 1, 2,..., n do 6: w k = Av k β k 1 v k 1 7: α k = (w k, v k ) 8: w k = w k α k v k 9: β k = w k 2 10: if k > 2 then [ ] % apply G k 1 G k 2 to the last column of T k+1,k [ ] [ ] 11: τ (3) k 2 c k 2 s k 2 0 = β k 1 s k 2 c k 2 β k 1 12: end if

92 13: if k [ > 1 then ] [ 14: τ (2) k 1 α k = 15: end if c k 1 s k 1 s k 1 c k 1 ] [ ] βk 1 α k 16: if α k > β k then % form the Givens rotation G k 17: set c k = 1+γ 1, s 2 k = c k γ where γ = β k / α k 18: else 19: set s k = 1 1+γ 2, c k = s k γ where γ = α k / β k 20: end if 21: τ (1) k = c k α k + s j β k % apply G k to last column of T k+1,k [ ] [ ] [ ] 22: ξ k c k s k = ξ k+1 s k c k ξ k 0 % apply G k to ξ ( ) 23: p k = p k 2 v k τ (2) k 1 p k 1 τ (3) k 2 /τ (1) k where p 0 = p 1 = 0 24: x (k) = x (k 1) + ξ k p k 92/115

93 93/115 25: relres = ξ k+1 /β 26: if relres< ε then % check convergence 27: break 28: end if 29: v k+1 = w k /β k 30: end for 31: if relres < ε then 32: 输出近似解 x (k) 及相关信息 33: else 34: 输出算法失败信息 35: end if 与 GMRES 相比, MINRES 方法的优点是

94 94/115 MINRES 充分利用了系数矩阵的对称性, 这使得每步迭代的运算量不会随着迭代步数的增加而增加 ; 在迭代过程中, 我们不需要存储所有的 v i 和 p i, 因此存储量也与迭代步数无关. 注记 CG 方法也可以用来求解对称不定线性方程组, 但可能会产生中断. 在 MINRES 方法中, 我们用了三项递推公式, 好处是可以有效节省运算量和存储量, 但缺点是, 随着迭代步数的增加, 舍入误差的影响会越来越大. Sleijpen, van der Vorst 和 Modersitzki 曾经指出, 在 MINRES 中, 舍入误差是按 κ(a) 2 的速度增长的. 而在 GMRES 中, 舍入误差只按 κ(a) 的

95 速度增长. 因此, 对于坏条件问题, 在使用 MINRES 需要加倍小心, 特别是迭代步数比较大的情况下. 此时, 我们也可以选择 SYMMLQ 方法. 95/115

96 96/ SYMMLQ 方法 在 SYMMLQ 方法中, 我们取 L k = K k, 即在仿射空间 x (0) + K k 中寻找近似解 x (k), 满足 x (k) x (0) + K k 且 b Ax (k) K k. SYMMLQ (Symmetric LQ) 方法的收敛速度可能不如 MINRES, 但对于坏条件问题, SYMMLQ 通常比 MINRES 更稳定. 注记 当 A 对称正定时, SYMMLQ 方法与 CG 方法是等价的.

97 设 x (k) 是 SYMMLQ 在 x (0) + K k 中找到的近似解, 则 x (k) 可表示为 x (k) = x (0) + V k y k (4.41) 其中 y k 是下面方程的解 : T k y = βe 1 (4.42) 如果 A 是不定的, 则 T k 也可能是不定的, 因此就不能用 LDL 分解来解方程组 (4.42). 此时, 我们采用 LQ 分解 ( 类似于 QR 分解 ), 即 T k = L k Q k (4.43) 其中 Lk R k k 是下三角矩阵, Q k R k k 是正交矩阵. 于是 x (k) = x (0) + V k T 1 k (βe 1) = x (0) + V k Q 1 k L k (βe 1) = x (0) + P k z (k), 97/115

98 其中 P k V k Q k Rn k, z (k) L 1 k (βe 1) R k. 由于 T k 是三对角矩阵, 因此 T k 的 LQ 分解可以通过一系列的 Givens 变换来说实现. 同样地, T k+1 的 LQ 分解可以在 T k 的 LQ 分解的基础上, 通过一次 Givens 变换得到. 下面给出具体推导过程. k = 1 时 此时 T k = α 1, 无需做任何分解, k = 2 时 [ ] 当 k = 2 时, T k = α 1 β 1 β 1 α 2. 只需一次 Givens 变化即可得到 T k 的 LQ 分解 : T k = L 2 G 1. 98/115

99 99/115 T k T k+1 设 T k 的 LQ 分解为 T k = L k (G 1 G 2 G k 1 ) L k Q k, 其中 Q k = (G 1 G 2 G k 1 ). 这里 G i 是 Givens 变换 : I i 1 c i s i G i = s i c i Rk k, i = 1, 2,..., k 1. I k i 1

100 100/115 由于 T k 是三对角, 因此 Lk = T k Q k 只有三条对角线非零, 即 L k = l (1) 1 l (2) 2 l (1) 2 l (3) 3 l (2) 3 l (1) l (3) k 1 l (2) k 1 l (1) k 1 l (3) k l (2) k l(1) k. 注记 这里最右下角的元素上面有个波浪号, 是用于区分 Lk 和 L k, 其中 L k 表示 Lk+1 的 k 阶顺序主子矩阵. 即 Lk 并不是 Lk+1 的子矩阵.

101 101/115 令 则 Q k+1 [ ] Q k [ ] T k+1 Q k+1 = T 0. k Q k 0 0 β 0 1 k 0 0 β k α k+1 其中 R (k+1) (k+1), = L k 0. 0 β k 0 0 l (3) k+1 β k α k+1 [ 0 0 l (3) k+1 β k ] = [ 0 0 β k ] Q k = [ 0 0 β k ] Gk 1. 即 l (3) k+1 = s k 1β k, β k = c k 1 β k. 构造 Givens 变换 G k, 消去 T k+1 Q k+1 倒数第二行的最后一个分量 β k,,

102 102/115 即 其中 c k = k 1 G k = I c k s k R (k+1) (k+1), s k c k (1) l k l (1) k, s k = β k l (1) k, l (1) k = ( l(1) ) 2 k + β 2 k. 令 Q k+1 = G k Q k+1. 由于右乘 G k 时只会影响到最后两列的元素, 所以 T k+1 Q k+1 = L 0. k 0 β k 0 0 l (3) β k+1 k α k+1 G k = l (1) 1 0 l (2) 2 l (1) 2 l (3) 3 l (2) 3 l (1) l (3) k l (2) k l (1) k l (3) k+1 l (2) k+1.. l(1) k+1 L k+1.

103 103/115 于是, 我们就得到 T k+1 的 LQ 分解 T k+1 = L k+1 Q k+1.

104 104/115 记 Lk 的 k 1 阶顺序主子矩阵为 L k 1, L k+1 的 k 阶顺序主子矩阵为 L k. 由于 L k 和 Lk 只有第 (k, k) 位置上的元素不同, 因此 L k 1 也是 L k 的 k 1 阶顺序主子矩阵. 于是我们就得到下面的递推关系 : L j 1 是 L j 的 j 1 阶顺序主子矩阵, j = 1, 2,... 令 z (k) = L 1 k (βe 1) R k. 由于 z (k) = [ ] z (k) = z (k 1) z k L 1 k, z (k) = (βe 1), 所以 [ ] z (k 1) 其中 z (k 1) = L 1 k 1 (βe 1), z k 和 z k 分别表示 z (k) 和 z (k) 的最后一个元素. 因此, 我们可以将 z (k) 和 z (k) 写为 z (k) = [z 1,..., z k 1, z k ], z (k) = [z 1,..., z k 1, z k ], k = 1, 2,..., z k,

105 105/115 其中 z k = l(1) k l(1) k 由 L k z (k) = βe 1 可知 z 1 = β/l (1) 1, z 2 = l (2) 2 z 1 /l (1) 2, ( ) z i = l (3) i z i 2 + l (2) i z i 1 /l (1) i, i = 3, 4,..., k. z k. (4.44) 分别记 Pk = V k Q k 和 Pk+1 = V k+1 Q k+1 的前 k 1 列和前 k 列组成的矩阵为 P k 1 和 P k, 则

106 106/115 [ ] P k+1 = V k+1 Q k+1 = [V Q k 0 k, v k+1 ] 0 1 [ ] = Pk, v k+1 G k G k k 1 = [P k 1, p k, v k+1 ] I c k s k (4.45) s k c k = [P k 1, p k, p k+1 ] = [P k, p k+1 ], 其中 p k 和 p k+1 分别表示 Pk 和 Pk+1 的最后一列. 因此, P k 和 P k 1 有下面的递推关系 P k = [P k 1, p k ],

107 107/115 其中 p k 表示 P k 的最后一列. 所以, 我们可以将 P k 统一写为 P k = [p 1, p 2,..., p k ], k = 1, 2,.... 易知 p 1 = P 1 = v 1. 由 (4.45) 可知 p 1 = v 1, p k = c k p k s k v k+1, p k+1 = s k p k + c k v k+1, k = 1, 2,.... (4.46) 令 x (k) = x (0) + P k z (k), k = 1, 2,..., 则 [ x (k) = x (0) +P k z (k) = x (0) +[P k 1, p k ] z (k 1) z k ] = x (k 1) +z k p k, k = 1, 2,.... (4.47)

108 108/115 于是有 x (k+1) = x (0) + P k+1 z (k+1) = x (0) + [P k, p k+1 ] [ z (k) z k+1 ] = x (0) + P k z (k) + z k+1 p k+1 = x (k) + z k+1 p k+1. (4.48) 有了这个关系式后, 我们在实际计算中只需计算 x (k). 由定理 5 可知 r k = b Ax (k) = β k (e k y k)v k+1. (4.49) 因此, 为了计算残量的值, 我们需要计算 y k, 即解方程 (4.42). 但事实上, 我们只需知道 y k 最后一个分量即可, 无需把整个 y k 都计算出来. 我们注意到 T k 是对称的, 因此 T k y k = T k y k = βe 1, 所以 L k y k = Q k (βe 1 ).

109 109/115 观察上述等式两边的最后一个元素可知 l(1) k e k y k = e k L k y k = e k Q k(βe 1 ) = βe k (G 1G 2 G k 1 ) e 1 = β(g 1 G 2 G k 1 e k ) e 1 = βs 1 s 2 s k 1. 由 Givens 变换 G k 的具体构造公式可知, s k l(1) k + c k β k = 0. 所以 ( ) r k = β k (e k y k)v k+1 = βs 1 s 2 s k 1 β k / l (1) k v k+1 = (βs 1 s 2 s k 1 s k /c k )v k+1, 即 r k 2 = βs 1 s 2 s k 1 s k /c k = c k 1 s k /c k r k 1 2.

110 110/115 算法 4.4 SYMMLQ 方法 1: 给定初值 x (0) 和 ( 相对 ) 精度要求 ε > 0 2: 计算 r 0 = b Ax (0) 和 β = r 0 2 3: v 1 = r 0 /β, p 1 = v 1, ξ 0 = β, x (0) = x (0) 4: for k = 1, 2,..., do 5: w k = Av k β k 1 v k 1 where β 0 = 0 and v 0 = 0 6: α k = (w k, v k ) 7: w k = w k α k v k 8: β k = w k 2 9: if k = 1 then 10: l(1) k 11: end if = α k 12: if k = 2 then 13: βk 1 = β k 1 14: end if

111 111/115 15: if k > 2 then % apply G k 2 to the last row of T k [ ] [ ] [ ] 16: l (3) c k 2 s k 2 k β k 1 = 0 β k 1 s k 2 c k 2 17: end if 18: if k > 1 then % form the Givens rotation G k 1 19: if l (1) k 1 > β k 1 then 20: set c k 1 = 1, s 1+γ 2 k 1 = c k 1 γ where γ = β k 1 / l (1) 21: else 22: set s k 1 = 1 1+γ 2 k 1 = s k 1 γ where γ = 23: end if 24: end if k 1 l (1) k 1 /β k 1 25: if k > 1 then % apply G k 1 to the last two columns of T k Qk 1 ( l(1) ) 2 26: l (1) k 1 = k 1 + β 2 k 1

112 112/115 27: 28: end if [ l (2) k l(1) k ] [ ] [ ] c k 1 s k 1 = βk 1 α k s k 1 c k 1 29: % compute z k 30: if k = 2 then 31: z 1 = β/l (1) 1 32: end if 33: if k = 3 then 34: z 2 = l (2) 2 z 1 /l (1) 2 35: end if 36: if k > 3 then ( ) 37: z k 1 = l (3) k 1 z k 3 + l (2) k 1 z k 2 /l (1) k 1 38: end if 39: % compute p k 40: if k = 1 then

113 113/115 41: p 1 = v 1 42: end if 43: if k > 1 then 44: p k 1 = c k 1 p k 1 s k 1 v k 45: p k = s k 1 p k 1 + c k 1 v k 46: end if 47: % update x (k) 48: if k > 1 then 49: x (k 1) = x (k 2) + z k 1 p k 1 50: end if 51: % check convergence 52: if k > 1 then 53: ξ k 1 = (c k 2 s k 1 /c k 1 ) ξ k 2 54: if ξ k 1 < ε then ( ) 55: x (k) = x k 1 + l (1) k z k/ l (1) k p k

114 114/115 56: stop 57: end if 58: end if 59: v k+1 = w k /β k 60: end for 注记 由 (4.49) 可知, 在 SYMMLQ 方法中, 残量是相互正交的. 注记 如果 A 是对称正定的, 则 SYMMLQ 与 CG 等价, 此时 SYMMLQ 极小化 x (k) x A. 如果 A 是不定的, 则没有这个最优性质. 但可以证明, SYMMLQ 在仿射空间 x (0) + AK k (A, r 0 ) 中极小化 x (k) x 2, 参见 [?? ].

115 115/115