04_Krylov 子空间方法 (II)

1 .第四讲 Krylov 子空间方法 (II) 5 收敛性分析 6 基于双正交化过程的迭代方法 7 免转置迭代方法 8 正规方程的迭代方法 

2 .5 收敛性分析 5.1 CG 方法的收敛性 5.2 GMRES 方法的收敛性 理论上, 若不考虑舍入误差, 则 Krylov 子空间方法都将在 n 步内终止. 但实际应用中, 我们通常希望收敛所需的迭代步数远远小于 n. 另一方面, 由于舍入误差的存在, 使得方法不一定都能在 n 步内收敛. 一般来讲, 如果 A 是可以酉 对角化的, 则 CG 方法和 GM- 收敛性分析不仅仅给出理论结果, 还有助于寻求改善迭代的方法 RES 方法的收敛性主要取决 于 A 的特征值分布. 但对于一般矩阵, 则比较复 本节主要讨论 CG 方法和 GMRES 方法的收敛性. 杂, 目前还没有定论. 2/76 

3 .5.1 CG 方法的收敛性 设 A 对称正定 . 令 x(k) 是 CG 方法在仿射空间 x(0) + Kk (A, r0 ) 中找到的近似解. 则由 CG 方法的最优性质可知 ∥x(k) − x∗ ∥A = min ∥x − x∗ ∥A (4.1) x∈x(0) +Kk (A,r0 ) 记 Pk 为所有次数不超过 k 的实系数多项式集合. 则 Kk (A, r0 ) 中任何 一个向量都可以写成 p(A)r0 形式, 其中 p ∈ Pk−1 . 设 x ∈ x(0) + Kk (A, r0 ), 则存在多项式 p(t) ∈ Pk−1 使得 x = x(0) + p(A)r0 . 3/76 

4 .记 ϵ0 ≜ x(0) − x∗ , 则 x − x∗ = x(0) + qk (A)r0 − x∗ = ϵ0 + p(A)(b − Ax(0) ) = ϵ0 + p(A)(Ax∗ − Ax(0) ) = (I − Ap(A))ϵ0 . 于是 ∥x − x∗ ∥2A = ∥(I − Ap(A))ϵ0 ∥2A ⊺ ⊺ = ϵ0 (I − Ap(A)) A(I − Ap(A))ϵ0 ⊺ ⊺ ≜ ϵ0 q(A) Aq(A)ϵ0 , (4.2) 其中 q(t) = 1 − t p(t) ∈ Pk 满足 q(0) = 1. 4/76 

5 .令 A = QΛQ⊺ 是 A 的特征值分解, 其中 Q 是正交矩阵, Λ = diag(λ1 , λ2 , . . . , λn ), λi > 0. 记 y = [y1 , y2 , . . . , yn ]⊺ ≜ Q⊺ ϵ0 . 由 (4.1) 和 (4.2) 可知 ⊺ ⊺ ∥x(k) − x∗ ∥2A = min ∥x − x∗ ∥2A = min ϵ0 q(A) Aq(A)ϵ0 x∈x(0) +Kk (A,r0 ) q∈Pk , q(0)=1 ⊺ ⊺ ⊺ = min ϵ0 Qq(Λ) Λq(Λ)Q ϵ0 q∈Pk , q(0)=1 ⊺ ⊺ = min y q(Λ) Λq(Λ)y q∈Pk , q(0)=1 ∑ n = min yi2 λi q(λi )2 q∈Pk , q(0)=1 i=1 ∑ n ≤ min max {q(λi )2 } yi2 λi q∈Pk , q(0)=1 1≤i≤n i=1 ⊺ = min max {q(λi ) } y Λy2 q∈Pk , q(0)=1 1≤i≤n 5/76 

6 . ⊺ = min max {q(λi )2 } ϵ0 Aϵ0 q∈Pk , q(0)=1 1≤i≤n = min max {q(λi )2 } ∥ϵ0 ∥2A , q∈Pk , q(0)=1 1≤i≤n 即 ∥x(k) − x∗ ∥2A ≤ min max {q(λi )2 }. ∥x(0) − x∗ ∥2A q∈Pk , q(0)=1 1≤i≤n 定理 1 设 A 对称正定, 则由 CG 方法得到的近似解 x(k) 满足 ∥x(k) − x∗ ∥A ≤ min max |q(λi )|. (4.3) ∥x(0) − x∗ ∥A q∈Pk , q(0)=1 1≤i≤n 需要指出的是, (4.3) 中上界是紧凑的, 即对任意 k, 总存在某个右端项 b 或某个初始值 x(0) (与 k 和 A 有关), 使得 (4.3) 中的等号成立. 也就是 说, (4.3) 中的上界描述了 CG 方法在最坏情况下的收敛情况. 6/76 

7 .由于 (4.3) 中的上界依赖于 A 的所有特征值. 在通常情况下, A 的特征 值是很难计算的. 因此该结论的实用性不大. 但是, 如果我们知道 A 的最大和最小特征值的话, 则可进一步缩放为 ∥x(k) − x∗ ∥A ≤ min max |q(λi )| ≤ min max |q(λ)|. ∥x(0) − x∗ ∥A q∈Pk , q(0)=1 1≤i≤n q∈Pk , q(0)=1 λmin ≤λ≤λmax 然后利用 Chebyshev 多项式的性质, 我们可以得到下面的结论. 定理 2 设 A 对称正定, 最大与最小特征值分别为 λmax 和 λmin , 则 (√ )k ∥x(k) − x∗ ∥A κ(A) − 1 ≤2 √ , (4.4) ∥x(0) − x∗ ∥A κ(A) + 1 其中 κ(A) = λmax /λmin 是 A 的 (谱) 条件数. (板书) 7/76 

8 .注记 上界 (4.3) 和 (4.4) 有着本质性的差异. 它们的值可能会相差很大. 上 界 (4.4) 描述的是区间 [λmin , λmax ] 的所有特征值分布的最坏情况. 所 以 (4.4) 并不能完全描述 CG 方法的收敛性质. 注记 A 的条件数越小, 则收敛越快. 但是, 如果 A 的条件数很大, 却并不 一定表明 CG 会收敛很慢. 注记 事实上, 在很多实际应用中, 人们观察到 CG 方法往往具有超收敛 性, 即方法的平均收敛速度随着迭代步数的增加而不断加快. 8/76 

9 .定理 3 (CG 方法的超收敛性) 设 A 对称正定, 且其特征值为 0 < λn ≤ · · · ≤ λn+1−i ≤ b1 ≤ λn−i ≤ · · · ≤ λj+1 ≤ b2 ≤ λj ≤ · · · ≤ λ1 . 则当 k ≥ i + j 时有 ( )k−i−j { ( )∏ j ( )} ∥x(k) − x∗ ∥A b−1 ∏ n λ − λℓ λℓ − λ ≤2 max , ∥x(0) − x∗ ∥A b+1 λ∈[b1 ,b2 ] ℓ=n+1−i λℓ ℓ=1 λℓ 其中 ( ) 12 b2 b= ≥ 1. b1 由定理结论可知, 当 b1 与 b2 非常接近时, i + j 步后, CG 收敛会非常快. 9/76 

10 .如果 A 的特征值聚集在 1 附近, 即除了有限多个特征值以外, A 的特 征值都在区间 [1 − ε, 1 + ε] 中, 其中 ε 很小, 则有下面的结论. 推论 4 设 A 对称正定, 特征值为 0 < δ ≤ λn ≤ · · · ≤ λn+1−i ≤ 1 − ε ≤ λn−i ≤ · · · ≤ λj+1 ≤ 1 + ε ≤ λj ≤ · · · ≤ λ1 . 则当 k ≥ i + j 时有 ( )i ∥x(k) − x∗ ∥A 1+ε ≤2 εk−i−j . (4.5) ∥x(0) − x∗ ∥A δ 10/76 

11 .5.2 GMRES 方法的收敛性 正规矩阵情形 设 A 是正规矩阵, 即 A = U ΛU ∗ , 其中 Λ = diag(λ1 , λ2 , . . . , λn ) 的对角线元素 λi ∈ C 为 A 的特征值. 设 x ∈ x(0) + Kk (A, r0 ), 则存在多项式 p(t) ∈ Pk−1 使得 x = x(0) + p(A)r0 . 于是 b − Ax = b − Ax(0) − Ap(A)r0 = (I − Ap(A))r0 ≜ q(A)r0 , (4.6) 其中 q(t) = 1 − t p(t) ∈ Pk 满足 q(0) = 1. 直接计算可知 ∥b − Ax∥2 = ∥q(A)r0 ∥2 = ∥U q(Λ)U ∗ r0 ∥2 ≤ ∥U ∥2 ∥U ∗ ∥2 ∥q(Λ)∥2 ∥r0 ∥2 = ∥q(Λ)∥2 ∥r0 ∥2 . 11/76 

12 .又 Λ 是对角矩阵, 所以 ∥q(Λ)∥2 = max |q(λi )|. 1≤i≤n 设 x(k) 是由 GMRES 方法得到的近似解. 由 GMRES 方法的最优性可 知, x(k) 极小化残量的 2 范数. 因此, ∥b − Ax(k) ∥2 = min ∥b − Ax∥2 x∈x(0) +Kk (A,r0 ) = min ∥q(A)r0 ∥2 q∈Pk , q(0)=1 ≤ min ∥q(Λ)∥2 ∥r0 ∥2 q∈Pk , q(0)=1 = ∥r0 ∥2 min max |q(λi )|. q∈Pk , q(0)=1 1≤i≤n 于是, 我们有下面的结论. 12/76 

13 .定理 5 设 A ∈ Rn×n 是正规矩阵, x(k) 是 GMRES 方法得到的近似 解, 则 ∥b − Ax(k) ∥2 ≤ min max |q(λi )|. (4.7) ∥r0 ∥2 q∈Pk , q(0)=1 1≤i≤n 注记 需要指出的是, 上界 (4.10) 是紧凑的. 与 CG 方法类似, 我们也可以选取一个包含 A 的所有特征值的区域 Ω ⊂ C, 然后将定理 5 中的上界缩放为 min max |q(λ)|. (4.8) q∈Pk , q(0)=1 λ∈Ω 通常 Ω 必须是连通的, 否则 (4.8) 的求解非常困难. 即使是两个区间的 并都没法求解. 另外, Ω 不能包含原点. 13/76 

14 .圆盘区域 首先考虑一种最简单的情形, 即 Ω 是复平面上的一个圆盘. 由于 A 是实的, 其特征值分布是关于实轴对称的, 所以可设 Ω 的圆心 为 C(c, 0), 半径为 r > 0. 这里假定 r < |c|, 即原点不在 Ω 内. ( )k c−z 定义多项式 qk (z) = , 则 qk (z) ∈ Pk 且 qk (0) = 1. 所以 c k c − λi rk min max |q(λi )| ≤ max |qk (λi )| = max ≤ . q∈Pk , q(0)=1 1≤i≤n 1≤i≤n 1≤i≤n c |c|k 由此可见, 当 r 越小或 c 越大时, 右式趋向于 0 的速度就越快. 注记 设 A 是正规矩阵, 其特征值包含在一个圆盘内, 则圆盘的半径越小 或者圆心离原点越远, 则 GMRES 收敛越快. 14/76 

15 .如果用椭圆盘来代替圆盘, 则可得出更紧凑的上界. 设 A 的特征值全部包含在椭圆盘 E(c, d, a) 内, 其中 d > 0 为焦距, a > 0 为半长轴长, C(c, 0) 为椭圆中心. 并且假定原点不在椭圆盘内. Im(z) Im(z) a a d o c Re(z) o c Re(z) d 15/76 

16 .令 Tk (z) 为 k 次复 Chebyshev 多项式, 则多项式 Tk ( c−z ) T˜k (z) = d c Tk ( d ) 就是极小极大问题 min max |q(λ)| (4.9) q∈Pk , q(0)=1 λ∈E(c,d,a) 的渐进最优解. 于是 min max |q(λi )| ≤ min max |q(λ)| ≤ max |T˜k (λ)|. q∈Pk , q(0)=1 1≤i≤n q∈Pk , q(0)=1 λ∈E(c,d,a) λ∈E(c,d,a) 假定椭圆盘 E(c, d, a) 的长轴在实轴上, 则由复 Chebyshev 多项式的性 质可知 Tk ( ad ) max |T˜k (λ)| = . λ∈E(c,d,a) |Tk ( dc )| 16/76 

17 .因此 Tk ( ad ) min max |q(λi )| ≤ . q∈Pk , q(0)=1 1≤i≤n |Tk ( dc )| 所以我们可得下面的结论. 定理 6 设 A 是正规矩阵, x(k) 是由 GMRES 得到的近似解, 则 虽然 T˜k (z) 通常不是极小极 ∥b − Ax(k) ∥2 Tk ( ad ) 大问题 (4.9) 的解, 但上界 (4.10) ≤ . (4.10) ∥r0 ∥2 |Tk ( dc )| 却往往能给出 GMRES 方法 收敛速度的一个很好的估计. 由 Chebyshev 多项式的性质可知, 上界 (4.10) 可以用下面的值来近似 ( √ )k a + a2 − d2 √ . c + c2 − d2 17/76 

18 .非正规情形 当 A 不是正规矩阵时, 在一般情况下, GMRES 方法的收敛性很难分析. 如果 A ∈ Rn×n 是可对角化的, 即 A = XΛX −1 , 其中 X ∈ Cn×n 非奇异, Λ = diag(λ1 , λ2 , . . . , λn ) 的对角线元素 λi ∈ C 为 A 的特征值, 则 ∥b − Ax(k) ∥2 = min ∥b − Ax∥2 = min ∥q(A)r0 ∥2 . x∈x(0) +Kk (A,r0 ) q∈Pk , q(0)=1 (4.11) 相类似地, 我们可以得到下面的结论. 18/76 

19 .定理 7 设 A = XΛX −1 , 其中 X ∈ Cn×n 非奇异, Λ 是对角矩阵, x(k) 是 GMRES 方法得到的近似解, 则 ∥b − Ax(k) ∥2 ≤ ∥X∥2 ∥X −1 ∥2 min max |q(λi )| ∥r0 ∥2 q∈Pk , q(0)=1 1≤i≤n = κ(X) min max |q(λi )|, (4.12) q∈Pk , q(0)=1 1≤i≤n 其中 κ(X) 是 X 的谱条件数. 如果 A 接近正规, 则 κ(X) ≈ 1. 此时上界 (4.12) 在一定程度上能描述 GMRES 的收敛速度. 但如果 X 远非正交, 则 κ(X) 会很大, 此时结论 中的上界就失去实际意义了. 需要指出的是, 上面的分析并不意味着非正规矩阵就一定比正规矩阵 收敛慢. 事实上, 对任意一个非正规矩阵, 总存在一个相应的正规矩阵, 使得 GMRES 方法的收敛速度是一样的 (假定初始残量相同). 19/76 

20 .虽然 GMRES 方法的收敛性与系数矩阵的特征值有关, 但显然并不仅 仅取决于特征值的分布. 事实上, 我们有下面的结论. 定理 8 对于任意给定的特征值分布和一条不增的收敛曲线, 则总存 在一个矩阵 A 和一个右端项 b, 使得 A 具有指定的特征值分布, 且 GMRES 方法的收敛曲线与给定的收敛曲线相同. 20/76 

21 .例 1 考虑线性方程组 Ax = b 其中   0 1    0 1     A= . . . . , b = e1 . . .     0 1  a0 a1 · · · an−1 当 a0 ̸= 0 时, A 非奇异. 易知, A 的特征值多项式为 p(x) = λn − an−1 λn−1 − an−2 λn−2 − · · · − a1 λ − a0 . 方程组的精确解为 x = [−a1 /a0 , 1, 0, . . . , 0]⊺ . 以零向量为迭代初值, 则 GMRES 迭代到第 n 步时才收敛. (前 n − 1 步残量范数不变). (GMRES_example01.m) 21/76 

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

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

24 .6 基于双正交化过程的迭代方法 6.1 双正交化过程 6.2 BiCG 方法 6.3 QMR 方法 GMRES 方法的一个典型缺点是: 随着迭代步数的增加, 每个迭代步 的运算量和存储量都会随之增加. 因此, 就有学者考虑, 能不能将 CG 方法的三项递推技巧运用于求解 非对称的线性方程组. 双正交化过程就是这种思想的一个具体实现. 24/76 

25 .6.1 双正交化过程 双正交化过程 (biorthogonalization) 也称为双边 Lanczos 过程 (two-sided Lanczos Process), 是指同时构造出子空间 K(A, r0 ) 的基 {v1 , v2 , . . .} 和 子空间 K(A⊺ , r˜0 ) 的基 {w1 , w2 , . . .}, 满足 这样做的优点是可以利用三 项递推公式, 从而节省运算 (vi , wj ) = 0, i ̸= j, 量. 但需要注意的是, 这两组 基本身并不保证正交. 即这两组基之间相互正交. 双正交化过程也可以看作是将对称矩阵的 Lanczos 过程推广到非对 称情形. 在 Lanczos 过程中, vk+1 是通过将 Avk 与 vk 和 vk−1 正交化后得到的. 类似地, 在双正交化过程中, vk+1 则是由 Avk 与 wk 和 wk−1 正交化后 所得到. 同时, wk+1 是由 A⊺ wk 与 vk 和 vk−1 正交化后得到. 25/76 

26 .算法 6.1 双正交化过程 1: 给定两个非零向量 r0 , r˜0 满足 (˜ r0 , r0 ) ̸= 0 2: 计算 β = ∥r0 ∥2 , 令 β0 = 0, γ0 = 0 3: v1 = r0 /β, w1 = β˜ r0 , r0 ) % make sure that w1⊺ v1 = 1 r0 /(˜ 4: for k = 1, 2, . . . , do 5: vˆk+1 = Avk − βk−1 vk−1 6: wˆk+1 = A⊺ wk − γk−1 wk−1 7: αk = (vk , w ˆk+1 ) % αk = wk⊺ Avk 8: v˜k+1 = vˆk+1 − αk vk 9: w˜k+1 = w ˆ − αk wk √ k+1 10: γk = |(˜ vk+1 , w ˜k+1 )| 11: if γk = 0 then 12: break % ghost break down 13: end if 26/76 

27 .14: βk = (˜ vk+1 , w ˜k+1 )/γk % so that βk γk = (˜ vk+1 , w ˜k+1 ) 15: vk+1 = v˜k+1 /γk 16: wk+1 = w ˜k+1 /βk 17: end for 两个关键耦合递推公式: γk vk+1 = Avk − αk vk − βk−1 vk−1 , (4.15) ⊺ βk wk+1 = A wk − αk wk − γk−1 wk−1 , k = 1, 2, . . . , m. (4.16) 27/76 

28 .向量组 {vk }m k=1 和 {wk }k=1 之间的双正交性 m 引理 9 设算法 6.1 的前 m−1 步不中断, 则向量组 {vk }m k=1 和 {wk }k=1 m 是双正交的, 即 ⊺ Wm Vm = I, (4.17) 其中 Vm = [v1 , v2 , . . . , vm ], Wm = [w1 , w2 , . . . , wm ]. (板书) 28/76 

29 .需要指出的是, 算法中 βk 和 γk 的选取是自由的, 但需要满足条件 βk γk = (˜ ˜k+1 ). 这个约束是为了确保 vk+1 , w (˜ vk+1 , w˜k+1 ) (vk+1 , wk+1 ) = = 1. β k γk 将耦合三项递推过程 (4.15) 和 (4.16) 写成矩阵形式即为 AVk = Vk+1 Tk+1,k = Vk Tk + γk vk+1 e⊺k , ⊺ k = 1, 2, . . . , m, A⊺ Wk = Wk+1 Tk+1,k = Wk Tk⊺ + βk wk+1 e⊺k ,   (4.18)   α1 β 1 0   .  γ1 . . . .. .   ..        其中 ek =   ∈ R , Tk+1,k =  k . .. . .. β  , 0  k−1     γk−1 αk  1 γk (k+1)×k Tk 是 Tk+1,k 的前 k 行. 29/76