04-Krylov 子空间方法

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

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

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

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

5 .如何计算近似解 设 {v1 , v2 , . . . , vm } 和 {w1 , w2 , . . . , wm } 分别是 K 和 L 一组基. 令 V = [v1 , v2 , . . . , vm ], W = [w1 , w2 , . . . , wm ]. W ∈ Rn×m , V ∈ Rn×m ˜ ∈ x(0) + K, 因此存在向量 y ∈ Rm 使得 x 由于 x ˜ = x(0) + V y ˆ≜x x ˜ − x(0) = V y 由正交性条件 (4.4) 可知 r0 − AV y ⊥ wi , i = 1, 2, . . . , m , 即 ⊺ ⊺ W AV y = W r0 . ⊺ 在实际计算中, W AV 通常 若 W ⊺ AV 非奇异, 则解得 y = (W ⊺ AV )−1 W ⊺ r0 . 因此 可以直接形成, 无需另外计 算. ⊺ ⊺ ˜ = x(0) + V (W AV )−1 W r0 x . 5/106 

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

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

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

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

10 .Arnoldi 过程: 计算 Km 的一组正交基 算法 2.1 基于 Gram-Schmidt 正交化的 Arnoldi 过程 1: 给定非零向量 r, 计算 v1 = r/∥r∥2 2: for j = 1, 2, . . . , m − 1 do 3: wj = Avj 4: for i = 1, 2, . . . , j do 5: hij = (wj , vi ) 6: end for ∑j 7: wj = wj − hij vi 做计算时一定要注意任何可 i=1 8: hj+1,j = ∥wj ∥2 能存在的 “中断”. 9: if hj+1,j = 0 then 10: break 11: end if 12: vj+1 = wj /hj+1,j 13: end for 10/106 

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

12 .Arnoldi 过程基本性质 根据上面的引理, 我们可以立即得到下面的结论. 定理 3 如果算法 2.1 不提前终止, 则向量 v1 , v2 , . . . , vm 构成 Km 的 一组标准正交基, 其中 Km (A, r) = span{r, Ar, . . . , Am−1 r}. 12/106 

13 .Arnoldi 过程重要性质 ∑ j 由 Arnoldi 过程可知 hj+1,j vj+1 = Avj − hij vi , 因此 ∑ j i=1 wj = wj − hij vi i=1 ∑ j vj+1 = wj /hj+1,j Avj = hj+1,j vj+1 + hij vi i=1   h1j    .  h1j  .   .       h2j  hj+1,j     = [v1 , v2 , . . . , vj+1 ]  .  = [v1 , . . . , vj+1 , vj+2 , . . . , vm+1 ]   .   .   0   .  hj+1,j  ..    0 = Vm+1 Hm+1,m (:, j), 13/106 

14 .其中 Vm+1 = [v1 , v2 , . . . , vm+1 ],   h11 h12 h13 ··· h1,m−1 h1,m h ···  wj = Avj  21 h22 h23 h2,m−1 h2,m    for i = 1, 2, . . . , j do  0 h32 h33 ··· h3,m−1 h2,m    hij = (wj , vi )  0 ···  Hm+1,m =  0 h43 h4,m−1 h4,m  ∈ R(m+1)×m end for ∑  . . .. .. ..  wj = wj − j  .. ..  hij vi  . . .  i=1   hj+1,j = ∥wj ∥2  0 0 0 · · · hm,m−1 hm,m  0 0 0 ··· 0 hm+1,m 这里 hij 是由 Arnoldi 过程所定义的. 14/106 

15 .Arnoldi 过程重要性质 (cont.) 定理 4 设 Vm = [v1 , v2 , . . . , vm ], 则 这两个公式在后面的算法推 ⊺ 导过程中非常重要 AVm = Vm+1 Hm+1,m = Vm Hm + hm+1,m vm+1 em (4.5) ⊺ Vm AVm = Hm (4.6) 其中 em = [0, . . . , 0, 1]⊺ ∈ Rm , Hm = Hm+1,m (1 : m, 1 : m) ∈ Rm×m , 即 Hm 是由 Hm+1,m 的前 m 行组成的 Hessenberg 矩阵. Hm A Vm = Vm+1 = Vm + Hm+1,m 15/106 

16 .Arnoldi 过程重要性质 (cont.) 如果 Arnoldi 方法提前终止, 则我们可以得到一个不变子空间. 定理 5 如果算法 2.1 在第 k (k ≤ m) 步终止, 即 hk+1,k = 0, 则有 一旦出现不变子空间, Arnoldi AVk = Vk Hk , 过程就会被中断, 这是好事 还是坏事? 即 Kk 是 A 的一个不变子空间. 16/106 

17 .考虑到数值稳定性, 通常使用 修正的 Gram-Schmidt 过程 (MGS) 算法 2.2 基于 MGS 的 Arnoldi 过程 1: 给定非零向量 r 2: v1 = r/∥r∥2 3: for j = 1, 2, . . . , m − 1 do 4: wj = Avj 标准 Gram-Schmidt 过程: 5: for i = 1, 2, . . . , j do wj = Avj 6: hij = (wj , vi ) for i = 1, 2, . . . , j do 7: wj = wj − hij vi hij = (wj , vi ) 8: end for end for ∑j 9: hj+1,j = ∥wj ∥2 wj = wj − hij vi i=1 10: if hj+1,j = 0 then 11: break 12: end if 13: vj+1 = wj /hj+1,j 14: end for 17/106 

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

19 .3 一般线性方程组的 Krylov 子空间方法 给定一个迭代初始值 x(0) ∈ Rn , 令 Km = span{r0 , Ar0 , A2 r0 , . . . , Am−1 r0 }, 其中 r0 = b − Ax(0) 为初始残量. 下面介绍分别基于 L = Km 和 L = AKm 的投影方法: • 完全正交方法: L = Km • 广义极小残量法: L = AKm 19/106 

20 .3.1 完全正交方法 (FOM): L = Km 寻找 ˜ ∈ x(0) + Km x 满足 b − A˜ x ⊥ Km 由 Arnoldi 过程可知 ⊺ ⊺ ⊺ Vm AVm = Hm 且 Vm r0 = Vm (βv1 ) = βe1 , 其中 β = ∥r0 ∥2 , e1 = [1, 0, . . . , 0]⊺ ∈ Rm . ˜ ∈ x(0) + Km 可知, 存在向量 y˜ ∈ Rm 使得 由x 矩阵 Vm 的列向量构成 Km (0) 的一组基 x ˜=x + Vm y˜. 由正交性条件 b − A˜ x ⊥ Km 可得 与 Km 正交等价于与 Km 的 ⊺ Vm (b − A˜ x) = 0, 一组基正交 20/106 

21 .即 ⊺ ⊺ ⊺ 0 = Vm (b − Ax(0) − AVm y˜) = Vm r0 − Vm AVm y˜ = βe1 − Hm y˜. −1 如果矩阵 Hm 非奇异, 则 y˜ = βHm e1 , 所以 −1 ˜ = x(0) + βVm Hm x e1 21/106 

22 .算法 3.1 完全正交方法 (FOM) 1: 给定初值 x(0) , 计算 r0 = b − Ax(0) , β = ∥r0 ∥2 和 v1 = r0 /β 2: for j = 1, 2, . . . , m − 1 do 3: wj = Avj 4: for i = 1, 2, . . . , j do 5: hij = (wj , vi ) 6: wj = wj − hij vi 7: end for 8: hj+1,j = ∥wj ∥2 9: if hj+1,j = 0 then 10: m = j, break 11: end if 12: vj+1 = wj /hj+1,j 13: end for 14: 求解线性方程组 Hm y˜ = βe1 并计算近似解 x ˜ = x(0) + Vm y˜ 22/106 

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

24 .FOM (cont.) 残量估计 ˜ ∈ x(0) + Km 由算法 3.1 得到的近似解, 则残量表达式为 定理 6 设 x ⊺ r˜ ≜ b − A˜ x = −hm+1,m (em y˜) vm+1 −1 其中 y˜ = βHm e1 , em = [0, . . . , 0, 1]⊺ ∈ Rm . 因此, ⊺ ∥˜ r∥2 = hm+1,m |em y˜| = hm+1,m |˜ y (m)| (板书) AVm = ⊺ Vm Hm + hm+1,m vm+1 em 24/106 

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

26 .算法 3.2 实用 FOM 方法 1: 给定初值 x(0) 和 (相对) 精度要求 ε > 0 2: 计算 r0 = b − Ax(0) , β = ∥r0 ∥2 和 v1 = r0 /β 3: for j = 1, 2, . . . do 4: wj = Avj 5: for i = 1, 2, . . . , j do 6: hij = (wj , vi ) 7: wj = wj − hij vi 8: end for 9: hj+1,j = ∥wj ∥2 10: 求解方程 Hj y˜ = βe1 11: if hj+1,j |˜ y (j)|/β < ε then % relative residual 12: break 13: end if 14: vj+1 = wj /hj+1,j 15: end for 16: 计算近似解 x ˜ = x(0) + Vj y˜ 26/106 

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

28 .GMRES 的推导过程 任意向量 x ∈ x(0) + Km 都可表示为 x = x(0) + Vm y, 其中 y ∈ Rm , 故 b − Ax = b − A(x(0) + Vm y) = r0 − AVm y = βv1 − Vm+1 Hm+1,m y = Vm+1 (βe1 − Hm+1,m y), (4.9) 其中 e1 = [1, 0, . . . , 0]⊺ ∈ Rm+1 . 又 Vm+1 的列向量标准正交的, 所以 ∥b − Ax∥2 = ∥Vm+1 (βe1 − Hm+1,m y)∥2 = ∥βe1 − Hm+1,m y∥2 . (4.10) 28/106 

29 .GMRES 的推导过程 (cont.) 因此, 极小化问题 (4.8) 的解可表示为 ˜ = x(0) + Vm y˜ x , (4.11) 其中 y˜ 是下面最小二乘问题的解 min ∥βe1 − Hm+1,m y∥2 . (4.12) y∈Rm 注记 与 FOM 方法中求解一个线性方程组不同的是, 在 GMRES 方法中, 我们是求解一个最小二乘问题. 29/106