第七讲子空间迭代方法 1 Krylov 子空间 2 GMRES 算法 3 共轭梯度法 (CG) 4 收敛性分析 5 其它 Krylov 子空间迭代算法

Size: px

Start display at page:

Download "第七讲子空间迭代方法 1 Krylov 子空间 2 GMRES 算法 3 共轭梯度法 (CG) 4 收敛性分析 5 其它 Krylov 子空间迭代算法"

猛吴
5 years ago
Views:

1 第七讲子空间迭代方法 1 Krylov 子空间 2 GMRES 算法 3 共轭梯度法 (CG) 4 收敛性分析 5 其它 Krylov 子空间迭代算法

2 2/55 基本思想在一个维数较低的子空间中寻找解析解的一个最佳近似. 子空间迭代算法的主要过程可以分解为下面三步 : (1) 寻找合适的子空间 ; (2) 在该子空间中求最佳近似 ; (3) 若这个近似解满足精度要求, 则停止计算 ; 否则, 重新构造一个新的子空间, 并返回第 (2) 步. 这里主要涉及到的两个关键问题是 : (1) 如果选择和更新子空间 ; (2) 如何在给定的子空间中寻找最佳近似. 关于第一个问题, 目前较成功的解决方案就是使用 Krylov 子空间.

3 3/55 1 Krylov 子空间设 A R n n, r R n, 则由 A 和 r 生成的 m 维 Krylov 子空间定义为 K m = K m (A, r) span{r, Ar, A 2 r,..., A m 1 r}, m n. 设 dim K m = m, 令 v 1, v 2,..., v m 是 K m 的一组基, 则 x K m 可表示为 x = y 1 v 1 + y 2 v y m v m V m y 寻找最佳近似 x (m) 转化为 (1) 寻找一组合适的基 v 1, v 2,..., v m ; (2) 求出 x (m) 在这组基下面的表出系数 y (m).

4 4/55 基的选取 : Arnoldi 过程最简单的基 : {r, Ar, A 2 r,..., A m 1 r} 非正交, 稳定性得不到保证. Arnoldi 过程 : 将 {r, Ar, A 2 r,..., A m 1 r} 单位正交化 1: v 1 = r/ r 2 2: for j = 1 to m do 3: z = Av j 4: for i = 1 to j do % MGS 正交化过程 5: h i,j = (v i, z), z = z h i,j v i 6: end for 7: h j+1,j = z 2 % if h j+1,j = 0, break, endif 8: v j+1 = z/h j+1,j 9: end for

5 5/55 Arnoldi 过程的矩阵表示记 V m = [v 1, v 2,..., v m ], h 1,1 h 1,2 h 1,3 h 1,m h 2,1 h 2,2 h 2,3 h 2,m h 3,2 h 3,3 h 3,m H m+1,m = h m,m 1 h m,m R (m+1) m, h m+1,m 则由 Arnoldi 过程可知 Av j = h 1,j v 1 + h 2,j v h j,j v j + h j+1,j v j+1. 所以有 AV m = V m+1 H m+1,m = V m H m + h m+1,m v m+1 e T m, (7.1)

6 6/55 其中 H m = H m+1,m (1 : m, 1 : m), e m = [0,..., 0, 1] T R m. 由于 V m 是列正交矩阵, 上式两边同乘 V T m 可得 V T m AV m = H m. (7.2) 等式 (7.1) 和 (7.2) 是 Arnoldi 过程的两个重要性质.

7 7/55 Lanczos 过程若 A 对称, 则 H m 为对称三对角, 记为 T m, 即 α 1 β 1. β T m =. (7.3) βm 1 β m 1 α m Lanczos 过程的性质与三项递推公式 ( 令 v 0 = 0 和 β 0 = 0) AV m = V m T m + β m v m+1 e T m, V T m AV m = T m. (7.4) (7.5) β j v j+1 = Av j α j v j β j 1 v j 1, j = 1, 2,...,

8 8/55 Lanczos 过程 1: Set v 0 = 0 and β 0 = 0 2: v 1 = r/ r 2 3: for j = 1 to m do 4: z = Av j 5: α j = (v j, z) 6: z = z α j v j β j 1 v j 1 7: β j = z 2 8: if β j = 0 then break, end if 9: v j+1 = z/β j 10: end for

9 9/55 Krylov 子空间算法的一般过程 (1) 令 m = 1 (2) 定义 Krylov 子空间 K m (A, r 0 ); (3) 找出仿射空间 x (0) + K m 中的最佳近似解 ; (4) 如果这个近似解满足精度要求, 则迭代结束 ; 否则令 m m + 1, 返回第 (2) 步.

10 10/55 Krylov 子空间迭代算法基本框架 1: 选取初始向量 x (0) 2: 计算 r 0 = b Ax (0), v 1 = r 0 / r 0 2 3: 寻找最佳近似解 : x (1) x (0) + K 1 = x (0) + span{v 1 } 4: if x (1) 满足精度要求 then 5: 终止迭代 6: end if 7: for m = 2 to n do 8: 调用 Arnoldi 或 Lanczos 过程计算向量 v m 9: 寻找最佳近似解 : x (m) x (0) +K m = x (0) +span{v 1,..., v m } 10: if x (m) 满足精度要求 then 11: 终止迭代 12: end if 13: end for

11 11/55 如何计算 x (0) + K m 中的最佳近似 x (m) 首先, 我们必须给出最佳的定义, 不同的定义会导致不同的算法. 最直接的方式 : x (m) x 2 达到最小. 但由于 x 不知道, 因此不实用. 什么是最佳 (1) r m 2 = b Ax (m) 2 达到最小 A 对称 MINRES, A 非对称 GMRES (2) A 对称正定, 极小化 x x (m) A CG ( 共轭梯度法 ) 本讲主要介绍 GMRES 算法和 CG 算法.

12 2 GMRES 算法 GMRES 算法是目前求解非对称线性方程组的最常用算法之一. 最佳近似解的判别方法为使得 r m 2 = b Ax (m) 2 最小对任意向量 x x (0) + K m, 可设 x = x (0) + V m y, 其中 y R m. 于是 r = b Ax = r 0 AV m y = V m+1 (βe 1 H m+1,m y), 这里 β = r 0 2. 由于 V m+1 列正交, 所以 r 2 = V m+1 (βe 1 H m+1,m y) 2 = βe 1 H m+1,m y 2. 于是最优性条件就转化为 y (m) = arg min y R m βe 1 H m+1,m y 2. (7.6) 用基于 Givens 变换的 QR 分解来求解即可. 12/55

13 13/55 GMRES 算法的基本框架算法 2.1 GMRES 迭代算法基本框架 1: 选取初值 x (0), 停机标准 ε > 0, 以及最大迭代步数 IterMax 2: r 0 = b Ax (0), β = r 0 2 3: v 1 = r 0 /β 4: for j = 1 to IterMax do 5: w = Av j 6: for i = 1 to j do % Arnoldi 过程 7: h i,j = (v i, w) 8: w = w h i,j v i 9: end for 10: h j+1,j = w 2 11: if h j+1,j = 0 then 12: m = j, break

14 14/55 13: end if 14: v j+1 = w/h j+1,j 15: relres = r j 2 /β % 相对残量 16: if relres< ε then % 检测是否收敛 17: m = j, break 18: end if 19: end for 20: 解最小二乘问题 (7.6), 得到 y (m) 21: x (m) = x (0) + V m y (m)

15 15/55 实施细节需要解决下面两个问题 : (1) 如何计算残量 r m b Ax (m) 的范数? (2) 如何求解最小二乘问题 (7.6)? 这两个问题可以同时处理.

16 16/55 最小二乘问题的求解设 H m+1,m 的 QR 分解为 H m+1,m = Q T m+1r m+1,m, 其中 Q m+1 是正交矩阵, 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 R m 1 y, 0 2 其中 R m R m m 非奇异 ( 假定 H m+1,m 不可约 ). 所以 y (m) = βr 1 m q 1 (1:m), r m 2 = b Ax (m) 2 = βe 1 H m+1,m y (m) 2 = β q 1 (m + 1), 其中 q 1 (m + 1) 表示 q 1 的第 m + 1 个分量.

17 17/55 H m+1,m 的 QR 分解的递推计算方法由于 H m+1,m 是上 Hessenberg 矩阵, 因此我们采用 Givens 变换. [ ] (1) 当 m = 1 时, H 21 = h 11 h 21, 构造 Givens 变换 G 1 使得 [ ] G 1 H 21 = = R 21, 即 H 21 = G T 1 R (2) 假定存在 G 1, G 2,..., G 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 T mr m,m 1.

18 18/55 为了书写方便, 这里假定 G i 的维数自动扩张, 以满足矩阵乘积的需要. (3) 考虑 H m+1,m 的 QR 分解. 易知 [ ] H m,m 1 h m H m+1,m =, 其中 h m = [h 1m, h 2m,..., h mm ] T. 0 h m+1,m 所以有 [ ] Q m 0 H m+1,m = 0 1 [ R m,m 1 Q m h m 0 h m+1,m ] = R m 1 hm 1 0 ĥ mm, 0 h m+1,m 其中 hm 1 是 Q m h m 的前 m 1 个元素组成的向量, ĥmm 是 Q m h m 的最后一个元素.

19 构造 Givens 变换 G m : I m G m = 0 c m s m R (m+1) (m+1), 0 s m c m 其中 c m = ĥm,m, s m = h m+1,m, h m,m h h m,m = ĥ2 m,m + h 2 m+1,m. 令 m,m 则 Q m+1 H m+1,m =G m [ ] Q m 0 Q m+1 = G m, 0 1 R m 1 hm 1 0 ĥ j,j = 0 h m+1,m R m 1 hm 1 0 hj,j 0 0 R m+1,m 所以可得 H m+1,m 的 QR 分解 H m+1,m = Q T m+1r m+1,m. 19/55

20 20/55 由 H m,m 1 的 QR 分解到 H m+1,m 的 QR 分解, 我们需要 (1) 计算 Q m h m, 即将之前的 m 1 个 Givens 变换作用到 H m+1,m 的最后一列的前 m 个元素上, 所以我们需要保留所有的 Givens 变换 ; (2) 残量计算 : r m 2 = βq 1 (m + 1) = βq m+1 (m + 1, 1), 即 G m G m 1 G 2 G 1 (βe 1 ) 的最后一个分量的绝对值. 由于在计算 r m 1 时就已经计算出 G m 1 G 2 G 1 (βe 1 ), 因此这里只需做一次 Givens 变换即可 ; (3) y (m) 的计算 : 当相对残量满足精度要求时, 需要计算 y (m) = R 1 m q 1 (1 : m), 而 q 1 即为 G m G m 1 G 2 G 1 (βe 1 ).

21 21/55 实用 GMRES 算法算法 2.2 实用 GMRES 算法 1: 给定初值 x (0), 停机标准 ε > 0, 最大迭代步数 IterMax 2: r 0 = b Ax (0), β = r 0 2 3: if β < ε then 4: 停止计算, 输出近似解 x (0) 5: end if 6: v 1 = r 0/β 7: ξ = βe 1 % 记录 q 1 8: for j = 1 to IterMax do 9: w = Av j 10: for i = 1 to j do % Arnoldi 过程 11: h i,j = (v i, w) 12: w = w h i,jv i 13: end for

22 22/55 14: h j+1,j = w 2 15: if h j+1,j = 0 then % 迭代中断 16: m = j, break 17: end if 18: v j+1 = w/h j+1,j 19: for i = 1 to j 1 do % 计算 G j 1 G 2G 1H j+1,j(1:j, j) [ ] [ ] [ ] 20: 21: end for h ij h i+1,j = c i s i s i c i h ij h i+1,j 22: if h jj > h j+1,j then % 构造 Givens 变换 G j 23: τ = h j+1,j/h jj, c j = 1/ 1 + τ 2, s j = c jτ 24: else 25: τ = h jj/h j+1,j, s j = 1/ 1 + τ 2, c j = s jτ 26: end if 27: h jj = c jh jj + s jh j+1,j % 计算 G jh j+1,j(1:j, j) 28: h j+1,j = 0

23 23/55 29: [ ξ j ξ j+1 ] [ c j = s j s j c j ] [ ] ξ j 0 % 计算 G j(βg j 1 G 2G 1e 1) 30: relres = ξ j+1 /β % 相对残量 31: if relres < ε then 32: m = j, break 33: end if 34: end for 35: m = j 36: y (m) = H(1 : m, 1 : m)\ξ(1 : m) % 最小二乘问题, 回代求解 37: x (m) = x (0) + V my (m) 38: if relres < ε then 39: 输出近似解 x 及相关信息 40: else 41: 输出算法失败信息 42: end if

24 24/55 GMRES 算法的中断在上面的 GMRES 算法中, 当执行到某一步时有 h j+1,j = 0, 则算法会中断 (breakdown). 如果出现这种中断, 则我们就可以找到精确解. 定理设 A R n n 非奇异且 r 0 0. 若 h i+1,i 0, i = 1, 2,..., k 1, 则 h k+1,k = 0 当且仅当 x (k) 是方程组的精确解. ( 不考虑舍入误差 ) ( 板书 )

25 25/55 带重启的 GMRES 算法由于随着迭代步数的增加, GMRES 算法的每一步所需的运算量和存储量都会越来越大. 因此当迭代步数很大时, GMRES 算法就不太实用. 重启技术事先设定一个重启迭代步数 k, 如果 GMRES 达到这个迭代步数时仍不收敛, 则计算出 x (0) + K k 中的最佳近似解 x (k), 然后令 x (0) = x (k), 重新开始新的 GMRES 迭代.

26 算法 2.3 带重启的 GMRES 算法 1: 设定重启步数 k ( n) 2: 给定初值 x (0), 停机标准 ε > 0, 最大迭代步数 IterMax 3: r 0 = b Ax (0), β = r 0 2 4: if β < ε then 5: 停止计算, 输出近似解 x = x (0) 6: end if 7: for iter=1 to ceil(itermax/k) do % 外循环 8: v 1 = r 0 /β 9: ξ = βe 1 10: for j = 1 to k do 11: 调用 GMRES 循环 12: end for 13: m = j 14: y (m) = H(1 : m, 1 : m)\ξ(1 : m) 26/55

27 27/55 15: x (m) = x (0) + V m y (m) 16: if relres < ε then % 收敛, 退出循环 17: break 18: end if 19: x (0) = x (m) % 重启 GMRES 20: r 0 = b Ax (0), β = r : end for 22: if relres < ε then 23: 输出近似解 x (m) 及相关信息 24: else 25: 输出算法失败信息 26: end if

28 28/55 带重启的 GMRES 算法需要注意的问题 (1) 如何选取合适的重启步数 k? 一般只能依靠经验来选取, 如 k = 20, 50. (2) 不带重启的 GMRES 算法能保证算法的收敛性, 但带重启的 GM- RES 算法却无法保证, 有时可能出现停滞现象 (stagnation).

29 29/55 3 共轭梯度法 (CG) 最佳近似 : x x (m) A 最小首先给出最佳近似解 x (m) 的一个性质. 定理设 A 对称正定, 则 x (m) = arg min x x A (7.7) x x (0) +K m 当且仅当 x (m) x (0) + K m 且 b Ax (m) K m. (7.8) ( 板书 )

30 30/55 Lanczos 过程 Lanczos 过程的三项递推公式 : AV m = V m+1 T m+1,m = V m T m + β m v m+1 e T m, Vm T AV m = T m, 其中 T m = tridiag(β i, α i+1, β i+1 ). 由前面的结论可知, 此时我们需要在 x (0) + K m 寻找最优解 x (m), 满足 b Ax (m) K m. (7.9) 下面就根据这个性质推导 CG 算法的迭代公式.

31 31/55 CG 算法的推导首先, 设 x (m) = x (0) + V m z (m), 其中 z (m) R m. 由 (7.9) 可知因此, 0 = V T m (b Ax (m) ) = V T m (r 0 AV m z (m) ) = βe 1 T m z (m), z (m) = T 1 m (βe 1 ). 设 T m 的 LDL T 分解为 T m = L m D m L T m. 于是 x (m) = x (0) +V m z (m) = x (0) +V m Tm 1 (βe 1 ) = x (0) +(V m L T m )(βdm 1 L 1 m e 1 ). 如果 x (m) 满足精度要求, 则计算结束. 否则我们需要计算 x (m+1) = x (0) + V m+1 T 1 m+1 (βe 1) = x (0) + (V m+1 L T m+1 )(βd 1 m+1 L 1 m+1 e 1). 这里 T m+1 = L m+1 D m+1 L T m+1.

32 32/55 记 P m V m L T m = [ p 1, p 2,..., p m ] R n m, y m βd 1 m L 1 m e 1 = [η 1,..., η m ] T R m. P m 和 y m 的递推关系式 ( 由 T m+1 的 LDL T 分解可得 ) P m+1 = V m+1 L T m+1 = [ P m, p m+1 ] y m+1 = βdm+1 1 L 1 m+1 e 1 = [ym, T η m+1 ] T, m = 1, 2,.... p m+1 的递推关系式 p m+1 = l m p m + v m+1.

33 x (m+1) 的递推关系式 [ x (m+1) = P m+1 y m+1 = [ P m, p m+1 ] y m η m+1 ] = x (m) + η m+1 p m+1 r m+1 的递推关系式 ( 收敛性判断 ) r m+1 = b Ax (m+1) = b A(x (m) + η m+1 p m+1 ) = r m η m+1 A p m+1 另一方面, 我们有 r m = b Ax (m) = r 0 AV m z (m) = β m (e T mz (m) )v m+1, 即 r m 与 v m+1 平行. 记 r m = τ m v m+1, 其中 τ 0 = β = r 0 2, τ m = β m (e T mz (m) ), m = 1, 2, /55

34 34/55 p m+1 的递推关系式 ( 定义 p m = τ m 1 p m ) p m+1 = τ m p m+1 = τ m (v m+1 l m p m ) = r m + µ m p m, (7.10) 其中 µ m = l m τ m /τ m 1, m = 1, 2,.... x (m+1) 和 r m+1 的新递推关系式 x (m+1) = x (m) + η m+1 p m+1 = x (m) + ξ m+1 p m+1, r m+1 = r m η m+1 A p m+1 = r m ξ m+1 Ap m+1, (7.11) (7.12) 其中 ξ m+1 = η m+1 /τ m, m = 1, 2,....

35 35/55 系数 ξ m+1 和 µ m 的计算方法引理下面的结论成立 : (1) r 1, r 2,..., r m 相互正交 ; (2) p 1, p 2,..., p m 相互 A- 共轭 (A- 正交 ), 即当 i j 时有 p T i Ap j = 0. ( 板书 ) 在等式 (7.10) 两边同时左乘 p T m+1a 可得 p T m+1ap m+1 = p T m+1ar m + µ m p T m+1ap m = rmap T m+1. 再用 rm T 左乘方程 (7.12) 可得 0 = rmr T m+1 = rmr T m ξ m+1 rmap T m+1.

36 36/55 于是 ξ m+1 = rt mr m r T map m+1 = rmr T m p T m+1 Ap. (7.13) m+1 等式 (7.10) 两边同时左乘 p T ma 可得 0 = p T map m+1 = p T mar m + µ m p T map m = µ m = rt map m p T map m. 为了进一步减少运算量, 将上式简化. 用 r T m+1 左乘方程 (7.12) 可得 r T m+1r m+1 = r T m+1r m ξ m+1 r T m+1ap m+1 = ξ m+1 r T m+1ap m+1. 于是 ξ m+1 = rt m+1r m+1 r T m+1 Ap m+1 = ξ m = rt mr m r T map m.

37 37/55 即 r T map m = r T mr m /ξ m, 于是 µ m = rt map m p T map m = rt mr m p T map m 1 ξ m = rt mr m r T m 1 r m 1 (7.14) 注意, 以上递推公式是从 m = 1 开始的. 因此 m = 0 时需要另外推导. 首先, 由 p 1 的定义可知 p 1 = P 1 = V 1 L T 1 = v 1 = p 1 = τ 0 p 1 = βv 1 = r 0. 其次, 由 Lanczos 过程可知 T 1 = α 1 = v T 1 Av 1. 注意到 β = r T 0 r 0, 于是 x (1) = x (0) + V 1 T1 1 (βe 1 ) = x (0) + β v1 T Av v 1 = x (0) + rt 0 r 0 1 p T 1 Ap p 1. 1

38 38/55 令 ξ 1 = rt 0 r 0 p T 1 Ap ( 注 : 之前的 ξ m+1 计算公式 (7.13) 只对 m 1 有定义 ), 则 1 当 m = 0 时关于 x (m+1) 的递推公式仍然成立. 最后考虑残量. 易知 r 1 = b Ax (1) = b Ax (0) rt 0 r 0 p T 1 Ap Ap 1 = r 0 ξ 1 Ap 1, 1 即当 m = 0 时关于 r m+1 的递推公式也成立.

39 39/55 共轭梯度法算法 3.1 共轭梯度法 1: 给定初值 x (0), 停机准则 ε > 0, 最大迭代步数 IterMax 2: r 0 = b Ax (0) 3: β = r 0 2 4: if β < ε then 5: 停止迭代, 输出近似解 x (0) 6: end if 7: for m = 1 to IterMax do 8: ρ = rm 1r T m 1 9: if m > 1 then 10: µ m 1 = ρ/ρ 0 11: p m = r m 1 + µ m 1 p m 1 12: else

40 13: p m = r 0 14: end if 15: q m = Ap m 16: ξ m = ρ/(p T mq m ) 17: x (m) = x (m 1) + ξ m p m 18: r m = r m 1 ξ m q m 19: relres = r m 2 /β 20: if relres < ε then 21: 停止迭代, 输出近似解 x (m) 22: end if 23: ρ 0 = ρ 24: end for 25: if relres < ε then 26: 输出近似解 x (m) 及相关信息 27: else 28: 输出算法失败信息 40/55

41 41/55 29: end if CG 算法的每个迭代步的主要运算为一个矩阵向量乘积和两个向量内积 ;

42 4 收敛性分析 CG 算法的收敛性设 x 是解析解, x (m) 是 CG 算法在 x (0) + K m 中找到的近似解, 即 x (m) = arg min x x A. x x (0) +K m 记 P k 为所有次数不超过 k 的多项式的集合. 对任意 x x (0) + K m, 存在 p(t) P m 1, 使得 x = x (0) + p(a)r 0. 于是有 x x = ε 0 + p(a)(b Ax (0) ) = ε 0 + p(a)(ax Ax (0) ) q(a)ε 0, 其中 ε 0 = x (0) x, 多项式 q(t) = 1 tp(t) P m 且 q(0) = 1. 所以 x x 2 A = ε T 0 q(a) T Aq(A)ε 0. 42/55

43 43/55 设 A = QΛQ T, Λ = diag(λ 1, λ 2,..., λ n ), 记 y = [y 1, y 2,..., y n ] T Q T ε 0. x (m) x 2 A = min x x (0) +K m x x 2 A = min q P m, q(0)=1 εt 0 Qq(Λ) T Λq(Λ)Q T ε 0 n = min yi 2 λ i q(λ i ) 2 q P m, q(0)=1 i=1 min max {q(λ i) 2 } q P m, q(0)=1 1 i n n yi 2 λ i i=1 = min q P m, q(0)=1 max 1 i n {q(λ i) 2 } y T Λy = min q P m, q(0)=1 max 1 i n {q(λ i) 2 } ε T 0 Aε 0 = min q P m, q(0)=1 max 1 i n {q(λ i) 2 } ε 0 2 A,

44 44/55 引理设 x (m) 是 CG 算法迭代 m 步后得到的近似解. 则 x (m) x A x (0) x A min q P m, q(0)=1 max q(λ i). 1 i n 当 A 的特征值不知道时, 可用区间代替, 即 x (m) x A min max q(λ). x (0) x A q P m, q(0)=1 λ n λ λ 1 由 Chebyshev 多项式的最佳逼近性质可知, 上式的解为 ( ) T 2t (λ1+λ n) ( ) m m λ 1 λ n 1 κ(a) + 1 q(t) = ) = q(t). T m ( λ1+λn 2 κ(a) 1 λ 1 λ n

45 45/55 定理设 A R n n 对称正定, x (m) 是 CG 算法迭代 m 步后得到的近似解. 则 ( ) m x (m) x A κ(a) 1 2, x (0) x A κ(a) + 1 其中 κ(a) = λ 1 /λ n.

46 46/55 CG 算法的超收敛性如果我们能够获得 A 的更多的特征值信息, 则能得到更好的误差限. 定理设 A 对称正定, 特征值为 0 < λ n λ n+1 i b 1 λ n i λ j+1 b 2 λ j λ 1. 则当 m i + j 时有 x (m) x A x (0) x A 2 ( b 1 b + 1 其中 b = (b 2 /b 1 ) ) m i j max λ [b 1,b 2] { n k=n+1 i ( ) j ( ) } λ λk λk λ λ k k=1 λ k (7.15) 由此可知, 当 b 1 与 b 2 非常接近时, 迭代 i + j 步后, CG 收敛会非常快!

47 47/55 推论设 A 对称正定, 特征值为 0 < δ λ n λ n+1 i 则当 m i + j 时有 1 ε λ n i λ j ε λ j λ 1. x (m) ( ) i x A 1 + ε 2 ε m i j. (7.16) x (0) x A δ

48 通常 Ω 必须是连通的, 否则求解非常困难, 即使两个区间的并都没法求解. 48/55 GMRES 算法的收敛性正规矩阵情形 : A = UΛU 定理设 A R n n 是正规矩阵, x (m) 是 GMRES 得到的近似解, 则 b Ax (m) 2 min r 0 max q(λ i). (7.17) 2 q P m, q(0)=1 1 i n 需要指出的是, 上界 (7.17) 是紧凑的. 设 Ω C 是包含 A 的所有特征值的一个区域 ( 不能包含原点 ), 则 b Ax (m) 2 r 0 2 min q P k, q(0)=1 max λ Ω q(λ).

49 49/55 非正规情形设 A R n n 可对角化, 即 A = XΛX 1, 则 b Ax (k) 2 = min b Ax 2 = min q(a)r 0 2. (7.18) x x (0) +K k (A,r 0) q P k, q(0)=1 相类似地, 我们可以得到下面的结论. 定理设 A = XΛX 1, 其中 X C n n 非奇异, Λ 是对角矩阵, x (k) 是 GMRES 算法得到的近似解, 则 b Ax (k) 2 X 2 X 1 2 min r 0 max q(λ i) 2 q P k, q(0)=1 1 i n 其中 κ(x) 是 X 的谱条件数. = κ(x) min max q(λ i), (7.19) q P k, q(0)=1 1 i n

50 50/55 如果 A 接近正规, 则 κ(x) 1. 此时上界 (7.19) 在一定程度上能描述 GMRES 的收敛速度. 当如果 X 远非正交, 则 κ(x) 会很大, 此时该上界就失去实际意义了. 需要指出的是, 上面的分析并不意味着非正规矩阵就一定比正规矩阵收敛慢. 事实上, 对任意一个非正规矩阵, 总存在一个相应的正规矩阵, 使得 GMRES 算法的收敛速度是一样的. 虽然 GMRES 算法的收敛性与系数矩阵的特征值有关, 但显然并不仅仅取决于特征值的分布. 事实上, 我们有下面的结论. 定理对于任意给定的特征值分布和一条不增的收敛曲线, 则总存在一个矩阵 A 和一个右端项 b, 使得 A 具有指定的特征值分布, 且 GMRES 算法的收敛曲线与给定的收敛曲线相同.

51 51/55 例考虑线性方程组 Ax = b 其中 A = , b = e 1. a 0 a 1 a 2 a n 1 当 a 0 0 时, A 非奇异. 易知, A 的特征值多项式为 p(x) = λ n a n 1 λ n 1 a n 2 λ n 2 a 1 λ a 0. 方程组的精确解为 x = [ a 1 /a 0, 1, 0,..., 0].

52 52/55 以零向量为迭代初值, 则 GMRES 迭代到第 n 步时才收敛. ( 前 n 1 步残量范数不变 ). 程序参见 GMRES_example01.m

53 53/55 如果 A 不可以对角化我们在分析 GMRES 算法的收敛性时, 通常会想办法用一个新的极小化问题来近似原来的极小化问题 (7.18). 当然, 这个新的极小化问题应该是比较容易求解的. 事实上, 我们有 b Ax (k) min q(a)r q P k, q(0)=1 = r 0 2 r 0 2 max v 2=1 min q(a)v 2 (7.20) q P k, q(0)=1 min q(a) 2. (7.21) q P k, q(0)=1 不等式 (7.20) 右端代表的是在最坏情况下的 GMRES 收敛性, 而且是紧凑的, 即它是所能找到的不依赖于 r 0 的最好上界. 但我们仍然不清楚, 到底是 A 的那些性质决定着这个上界 [? ].

54 54/55 可以证明, 当 A 是正规矩阵时, 上界 (7.20) 和 (7.21) 是相等的 [?? ]. 但是, 对于大多数非正规矩阵而言, 这两者是否相等或者非常接近, 迄今仍不太清楚. 最后需要指出的是, 算法的收敛性也依赖于迭代初值和右端项. 所以上定理中的上界描述的都是最坏情况下的收敛速度. 也就是说, 在实际计算中, 算法的收敛速度可能会比预想的要快得多.

55 55/55 5 其它 Krylov 子空间迭代算法对称非对称正规方程 CG (1952) MINRES (1975) SYMMLQ (1975) SQMR (1994) FOM (1981) GMRES (1984) BiCG (1976) QMR (1991) CGS (1989) BiCGStab (1992) TFQMR (1993) FGMRES (1993) 对称正定, 正交投影法 (Galerkin) 对称不定, 斜投影法 (Petrov-Galerkin) 对称不定对称不定正交投影法, Arnoldi 斜投影法 (Petrov-Galerkin), Arnoldi 双正交 (biorthogonalization) 双正交 (biorthogonalization) Transpose free Transpose free, smoother convergence than CGS Transpose free, smoother convergence than CGS CGLS (1982) 最小二乘 ( 法方程 ) LSQR (1982) 最小二乘 ( 法方程 )

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

第四讲 Krylov 子空间方法投影方法 Krylov 子空间与 Arnoldi 过程一般线性方程组的 Krylov 子空间方法对称线性方程的 Krylov 子空间方法收敛性分析基于双正交化过程的迭代方法免转置迭代方法正规方程的迭代方法 2/115 大规模稀疏线性方程组 Ax = b, A R n n, b R n. Krylov 子空间方法是子空间方法的成功代表首选方法 Krylov

第七讲 子空间迭代方法 1 Krylov 子空间 2 GMRES 算法 3 共轭梯度法 (CG) 4 收敛性分析 5 其它 Krylov 子空间迭代算法

第七讲子空间迭代方法 1 Krylov 子空间 2 GMRES 算法 3 共轭梯度法 (CG) 4 收敛性分析 5 其它 Krylov 子空间迭代算法