第八讲预处理技术

5/20
预处理方式的选择
若 A 对称正定, 用 CG 求解, 这三种方式的预处理效果基本一致. 但若 A 非对称 (特别是非正规) 情形, 效果可能会相差很大. 实际使用中, 选取哪种预处理方式, 根据问题本身和所用的算法来确定. 对于 GMRES, 右预处理比较合适, 因为预处理后的残量就是真实残量. † 需要指出的是, 在实际求解预处理后的方程组时, 我们并不需要显式 地计算 P −1 (除非 P −1 非常容易计算), 更不需要显式地计算 P −1 A.
1: 5: 6: 7: 8: 9: 10:
x(k) = x(k−1) + ξk pk rk = rk−1 − ξk Apk zk = zk−1 − αk P −1 Apk = P −1 rk (rk , zk ) ( zk , z k ) P µk = = ( z k −1 , z k −1 ) P (rk−1 , zk−1 ) pk+1 = zk + µk pk end for
14/20
右预处理
AP −1 u = b, x = P −1 u.
给定迭代初值 x0 , 则 u(0) = P x(0) . 相应的预处理初始残量为 r ˜0 = b − AP −1 u(0) = b − Ax(0) = r(0) . 因此无需计算 u(0) . 设 u(m) = u(0) + Vm ym 是迭代 m 步后的近似解, 则 x(m) = P −1 u(m) = P −1 (u(0) + Vm ym ) = x(0) + P −1 Vm ym , 真实残量为 rm = b − Ax(m) = r0 − AP −1 Vm ym = Vm+1 (βe1 − Hm+1,m ym ). 因此 ∥rm ∥2 = ∥βe1 − Hm+1,m ym ∥2 = β · |q1 (m + 1)|.
6/20
1.1
预处理 CG 算法
设 A ∈ Rn×n 对称正定, 并设 预处理子 P 也对称正定. 为了保证预处理后的系数矩阵仍然是对称正定的, 我们考虑使用两边预 处理方式, 及设 P 的 Cholesky 分解为 P = LLT . 于是我们得到下面的预处理方程组 L−1 AL−T u = L−1 b, x = L−T u.
16/20
两边预处理
如果预处理子 P 是由乘积形式给出的, 即 P = PL PR , 则可以构造下面的两边预处理方程组
−1 −1 −1 PL APR u = PL b, −1 x = PR u.
将 GMRES 用于求解上述线性方程组, 则可得到两边预处理 GMRES 方法. 需要指出的是, 对于两边预处理 GMRES 方法, 其预处理后的残量为 r ˜k = L−1 (b − Ax(k) ).
3/20
预处理子的构造
理论上讲, 任何一个非奇异矩阵都可以作为预处理子. 但一个好的预处理 子 P 通常需满足下面两个要求: (1) P −1 A 具有更好的特征值分布及/或更小的条件数; (2) P z = r 容易求解. 其中第一个要求是为了确保预处理后的线性方程组更容易求解, 也就是 说, 选取的预处理方法是有效的. 第二要求是因为在用 Krylov 子空间方法求解预处理后的方程组 (9.2) 时, 每步迭代都需要求解一个以 P 为系数矩阵的线性方程组, 为了不增加太 多额外的运算量, 必须很容易求解.
原始解与真实残量
首先引入一个辅助变量: pk L−T p ˜k , 于是有
pk+1 = (L )−1 pk+1 = (L )−1 r ˜k + µk (L )−1 p ˜k = (L )−1 r ˜k + µk pk . 又 rk = Lr ˜k , 因此可得递推公式: rk = Lr ˜k = L(˜ rk−1 − ξk L−1 A(L )−1 p ˜k ) = rk−1 − ξk Apk , pk+1 = (L )−1 r ˜k + µk pk = (L )−1 L−1 rk + µk pk = P −1 rk + µk pk , x(k) = (L )−1 u(k) = (L )−1 (u(k−1) + ξk p ˜k ) = x(k−1) + ξk pk , 其中 (L−1 rk−1 , L−1 rk−1 ) (rk−1 , P −1 rk−1 ) = , (L−1 A(L )−1 p ˜k , p ˜k ) (Apk , pk ) (L−1 rk , L−1 rk ) (rk , P −1 rk ) µ k = −1 = . (L rk−1 , L−1 rk−1 ) (rk−1 , P −1 rk−1 ) ξk =
如果在方程组的两边同时左乘一个非奇异矩阵 P ∈ Rn×n 的逆, 则可得 P −1 Ax = P −1 b. 这就是预处理后的方程组, P 就称为预处理子 (preconditioner). 用 Krylov 子空间方法求解 (9.2) 时, 就称为 预处理 Krylov 子空间方法. (9.2)
9/20
记 zk
P −1 rk , 即预处理残量, 则可得下面的算法.
算法 1.2 预处理 CG 方法 (PCG) 给定初值 x(0) , 计算 r0 = b − Ax(0) 2: 令 z0 = P −1 r0 , p1 = z0 , 计算 ρ = r0 z0 3: for k = 1, 2, . . . do 4: ρ0 = ρ, qk = Apk , ξk = ρ/(pk qk )
1: 5: 6: 7: 8: 9:
u(k) = u(k−1) + ξk p ˜k r ˜k = r ˜k−1 − ξk L−1 A(L )−1 p ˜k µk = (˜ rk , r ˜k ) (˜ rk−1 , r ˜k−1 )
p ˜k+1 = r ˜k + µk p ˜k end for
8/20
为了改善 krylov 迭代方法的收敛性与可靠性, 我们通常需要运用预处理 技术 (preconditioning). 通俗地说, 预处理就是将原来难以求解的问题转化成一个等价的但比较 容易求解的新问题. 预处理技术的研究是目前科学计算领域中的重要研究课题之一.
2/20
1
预处理迭代方法
对于线性方程组, 预处理就是通过适当的线性转换, 转换为一个新方程组. 考虑线性方程组 Ax = b, A ∈ Rn×n 非奇异, b ∈ Rn . (9.1)
17/20
2
预处理技术
Finding a good preconditioner to solve a given sparse linear system is often viewed as a combination of art and science. Theoretical results are rare and some methods work surprisingly well, often despite expectations. — Saad, 2003.
11/20
△
算法 1.3 基于 P -内积的 CG 方法 给定初值 x(0) , 计算 r0 = b − Ax(0) 2: 令 z0 = P −1 r0 , p1 = z0 3: for k = 1, 2, . . . do (zk−1 , zk−1 )P (rk−1 , zk−1 ) 4: ξk = = (P −1 Apk , pk )P (Apk , pk )
第八讲 预处理技术
1 2 预处理迭代方法 预处理技术
Nothing will be more central to computational science in the next century than the art of transforming a problem that appears intractable into another whose solution can be approximated rapidly. For Krylov subspace matrix iterations, this is preconditioning. — Trefethen & Bau, 1997.
1: 5: 6: 7: 8: 9: 10: 11:
x(k) = x(k−1) + ξk pk rk = rk−1 − ξk qk zk = P −1 rk ρ = rk zk µk = ρ/ρ0 pk+1 = zk + µk pk end for
注: 矩阵 L 并没有出现在算法中, 因此无需对 P 做 Cholesky 分解.
15/20
这与不带预处理的 GMRES 算法是一样的. 因此在实际求解过程中我们可 以直接得到原方程组真实残量的模, 而无需计算 u(m) 和 x(m) . 这是右预 处理方式与左预处理方式的主要区别. 算法 1.4 右预处理 GMRES 方法
1: 2: 3: 4: 5: 6: 7: 8: 9:
给定初值 x(0) , 计算 r0 = b − Ax(0) 令 v1 = r ˜0 /β , ξ = βe1 for j = 1, 2, . . . do w ˜j = P −1 vj % Apply preconditioning wj = Aw ˜j ... ... end for −1 (m) 计算 y (m) = Rm ξ , 其中 Rm = H (1 : m, 1 : m), ξ (m) = ξ (1 : m) 计算近似解 x(m) = x(0) + P −1 Vm y (m)
4/20
预处理方式
• 左预处理: P −1 Ax = P −1 b • 右预处理: AP −1 y = b, x = P −1 y x = R−1 y , P = LR • 两边预处理: L−1 AR−1 y = b,
这三种方式预处理后的系数矩阵分别为 P −1 A, AP −1 和 L−1 AR−1 . 由于它们是相似的, 所以有相同的特征值分布.
13/20
† 需要指出的是, 在左预处理 GMRES 算法中, 停机准则中的 ξj +1 并不 是原方程组真实残量的模, 而是预处理残量 r ˜k = P −1 (b − Ax(k) ) 的 模. 如果要计算真实残量, 除了显式计算外, 没有其它更好的办法. 所 以, 如果停机准则是基于真实残量的, 则需要额外计算真实残量.
用 CG 方法求解上述方程组, 迭代 k 步, 得到的近似解记为 u(k) , 预处理残量记为 r ˜k L−1 b − L−1 AL−T u(k) .
7/20
预处理 CG 方法
算法 1.1 两边预处理 CG 方法 给定初值 x(0) , 计算 r0 = b − Ax(0) 2: 令 r ˜0 = L−1 r0 , p ˜1 = r ˜0 3: for k = 1, 2, . . . do (˜ rk−1 , r ˜k−1 ) 4: ξk = −1 (L A(L )−1 p ˜k , p ˜k )
不难看出, 该算法与前面的 PCG 算法是完全一样的.
12/20
1.2
预处理 GMRES 算法
对于非对称 Krylov 方法, 三种预处理方式并不等价, 有时效果相差很大. 左预处理 : P −1 Ax = P −1 b.
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:
给定初值 x(0) , 计算 r ˜0 = P −1 (b − Ax(0) ) 令 v1 = r ˜0 /β , ξ = βe1 for j = 1, 2, . . . do w ˜j = Avj wj = P −1 w ˜j % Apply preconditioning for i = 1, 2, . . . , j do hij = (wj , vi ), wj = wj − hij vi end for ... ... end for ... ...
10/20
左预处理 CG
事实上, PCG 算法 1.2 也可以从左预处理方式导出. P −1 Ax = P −1 b. 易知 P −1 A 正定, 但通常不对称, 因此不能直接实施 CG 算法. 定义 P -内积, 即 (x, y )P = (P x, y ) = (x, P y ). 在该内积下, 有 (P −1 Ax, y )P = (Ax, y ) = (x, Ay ) = (x, P (P −1 Ay )) = (x, P −1 Ay ), 即 P −1 A 关于 P -内积自伴随 (在 P -内积意义下, P −1 A 是 “对称” 的) 将 CG 算法中的内积改为 P -内积, 即可得左预处理 CG 算法.