第八讲预处理Nothing will be more central to computational science in the next century than the art of transforminga problem that appears intractable into another whose solution can beapproximated rapidly.For Krylovsubspace matrix iterations, this is preconditioning._Trefethen&Bau[70],1997为了改善krylov选代方法的收敛性与可靠性，我们通常需要运用预处理技术（preconditioning).通俗地说，预处理就是将原来难以求解的问题转化成一个等价的但比较容易求解的新问题.预处理技术的研究是目前科学计算领域中的重要研究课题之一，8.1预处理送代方法对于线性方程组而言，预处理就是对系数矩阵进行适当的线性转换，转换为一个新矩阵考虑线性方程组Ar=b,AeRnxn非奇异，beRn.(8.1)如果在方程组的两边同时左乘一个非奇异矩阵PERnxn的逆，则可得P-1Ar = P-1b.(8.2)这个新方程组就是预处理后的方程组，P就称为预处理子（preconditioner）.当Krylov子空间方法被用于求解（8.2）时，就称为预处理Krylov子空间方法理论上讲，任何一个非奇异矩阵都可以作为预处理子，但一个好的预处理子P通常需满足下面两个要求：（1）P-1A具有更好的特征值分布及/或更小的条件数(2)Pz=r容易求解其中第一个要求是为了确保预处理后的线性方程组更容易求解，也就是说，选取的预处理方法是有效的.第二要求是因为在用Krylov子空间方法求解预处理后的方程组（8.2）时，每步选代都需要求解一个以P为系数矩阵的线性方程组，为了不增加太多额外的运算量，必须很容易求解预处理方式（8.2）称为左预处理：相应地，我们可以在方程组两边同时右乘P-1，这就是右预处理即AP-1u= b,r= p-lu(8.3)237
第八讲 预处理 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[70], 1997. 为了改善 krylov 迭代方法的收敛性与可靠性, 我们通常需要运用预处理技术 (preconditioning). 通俗 地说, 预处理就是将原来难以求解的问题转化成一个等价的但比较容易求解的新问题. 预处理技术的研 究是目前科学计算领域中的重要研究课题之一. 8.1 预处理迭代方法 对于线性方程组而言, 预处理就是对系数矩阵进行适当的线性转换, 转换为一个新矩阵. 考虑线性方程组 Ax = b, A ∈ R n×n 非奇异, b ∈ R n . (8.1) 如果在方程组的两边同时左乘一个非奇异矩阵 P ∈ R n×n 的逆, 则可得 P −1Ax = P −1 b. (8.2) 这个新方程组就是预处理后的方程组, P 就称为预处理子 (preconditioner ). 当 Krylov 子空间方法被用于 求解 (8.2) 时, 就称为预处理 Krylov 子空间方法. 理论上讲, 任何一个非奇异矩阵都可以作为预处理子. 但一个好的预处理子 P 通常需满足下面两 个要求: (1) P −1A 具有更好的特征值分布及/或更小的条件数; (2) P z = r 容易求解. 其中第一个要求是为了确保预处理后的线性方程组更容易求解, 也就是说, 选取的预处理方法是有效 的. 第二要求是因为在用 Krylov 子空间方法求解预处理后的方程组 (8.2) 时, 每步迭代都需要求解一个 以 P 为系数矩阵的线性方程组, 为了不增加太多额外的运算量, 必须很容易求解. 预处理方式 (8.2) 称为左预处理. 相应地, 我们可以在方程组两边同时右乘 P −1 , 这就是右预处理, 即 AP −1u = b, x = P −1u. (8.3) 237
.238.第八讲预处理另外，我们也可以将P分解成两个矩阵的乘积，即P=LR于是我们可以用下面的方式对原方程组(8.1)进行预处理L-1 AR-1u= b,r=R-lu(8.4)这就是两边预处理以上是三种常用的预处理方式.这三种方式预处理后的系数矩阵分别为P-1A，AP-1和L-1AR-1由于它们是相似的，所以具有相同的特征值分布，如果A是对称正定的，则使用共轭梯度法求解时，这三种方式的预处理效果基本上是一样的.但对于非对称（特别是非正规）情形，效果可能会相差很大在实际使用中，该选取哪种预处理方式，需要根据问题本身和所用的方法来确定，如对于对称正定线性方程组的CG方法，三种方式都可以，而对于GMRES方法，则选取右预处理比较合适.一方面是实际使用时，得到的残量范数与原方程组的残量范数是一样的，另一方面是，右预处理极小化的是原始残量范数，而左预处理极小化的是预处理后的残量这里需要指出的是，在实际求解预处理后的方程组时，我们并不会显式地计算P-1（除非P-1非常容易计算），更不会显式地计算P-1A8.1.1预处理CG方法设AERnxn对称正定，并假定预处理子P也是对称正定的.为了保证预处理后的系数矩阵仍然是对称正定的，我们考虑使用两边预处理方式.设P的Cholesky分解为P= LLT.于是我们得到下面的预处理方程组L-A(LT)-1u=-b,=(LT)-u(8.5)用CG方法求解上述方程组，选代步后，得到的近似解记为u()，预处理残量记为≤L-1b-L-1A(LT)-1u().于是，求解预处理方程组（8.5)的CG方法可描述如下：算法8.1.两边预处理CG方法1：给定初值(0)2:计算ro=b-Az(0)3:令To=L-1ro,P1=Fo4: for k =1,2,...do(Tk-1, Tk-1)5:Ek= (L-IA(LT)-1pk,Pk)u(k) = u(k-1) + Expk6:7:TK = TK-1 -SAL-1A(LT)-1pA(Fk,TA)8:μk=(Fk-1, Tk-1)9:Pk+1=Tk+μkPk
· 238 · 第八讲 预处理 另外, 我们也可以将 P 分解成两个矩阵的乘积, 即 P = LR. 于是我们可以用下面的方式对原方程组 (8.1) 进行预处理 L −1AR−1u = b, x = R −1u. (8.4) 这就是两边预处理. 以上是三种常用的预处理方式. 这三种方式预处理后的系数矩阵分别为 P −1A, AP −1 和L −1AR−1 . 由于它们是相似的, 所以具有相同的特征值分布. 如果 A 是对称正定的, 则使用共轭梯度法求解时, 这 三种方式的预处理效果基本上是一样的. 但对于非对称 (特别是非正规) 情形, 效果可能会相差很大. 在实际使用中, 该选取哪种预处理方式, 需要根据问题本身和所用的方法来确定. 如对于对称正定 线性方程组的 CG 方法, 三种方式都可以, 而对于 GMRES 方法, 则选取右预处理比较合适. 一方面是实 际使用时, 得到的残量范数与原方程组的残量范数是一样的, 另一方面是, 右预处理极小化的是原始残 量范数, 而左预处理极小化的是预处理后的残量. 这里需要指出的是, 在实际求解预处理后的方程组时, 我们并不会显式地计算 P −1 (除非 P −1 非常 容易计算), 更不会显式地计算 P −1A. 8.1.1 预处理 CG 方法 设 A ∈ R n×n 对称正定, 并假定预处理子 P 也是对称正定的. 为了保证预处理后的系数矩阵仍然是 对称正定的, 我们考虑使用两边预处理方式. 设 P 的 Cholesky 分解为 P = LL⊺ . 于是我们得到下面的预处理方程组 L −1A(L ⊺ ) −1u = L −1 b, x = (L ⊺ ) −1u. (8.5) 用 CG 方法求解上述方程组, 迭代 k 步后, 得到的近似解记为 u (k) , 预处理残量 记为 r˜k ≜ L −1 b − L −1A(L ⊺ ) −1u (k) . 于是, 求解预处理方程组 (8.5) 的 CG 方法可描述如下: 算法 8.1. 两边预处理 CG 方法 1: 给定初值 x (0) 2: 计算 r0 = b − Ax(0) 3: 令 r˜0 = L −1 r0, p˜1 = ˜r0 4: for k = 1, 2, . . . do 5: ξk = (˜rk−1, r˜k−1) (L−1A(L⊺)−1p˜k, p˜k) 6: u (k) = u (k−1) + ξkp˜k 7: r˜k = ˜rk−1 − ξkL −1A(L ⊺ ) −1p˜k 8: µk = (˜rk, r˜k) (˜rk−1, r˜k−1) 9: p˜k+1 = ˜rk + µkp˜k
8.1预处理选代方法.239.10:endfor由于算法8.1中得到的是预处理后的近似解u(）和预处理残量，因此我们还需要考虑原方程的近似解和残量的计算我们对上述算法中的选代公式进行适当的改写.首先引入一个辅助变量Pk：Pk ≤ (LT)-1pk,k =1,2,....于是有Pk+1 = (LT)-pk+1 = (LT)-1+(LT)-p=(LT)-1+Pk又由的表达式可知，原方程组的残量工三因此可得到下面的递推公式Tk=Li= L(Fk-1 -SAL-1A(LT)-1pk) = Tk-1 -EkApk,Pk+1 = (LT)-1k + μkpk =(LT)-L-1Tk + μkpk = P-1rk + μkpk2() = (LT)-1u(k) = (LT)-1(u(k-1) + EkPk) = (k-1) + EkPk,其中(L-1rk-1, L-rk-1)(rk-1,(LT)-1L-1rk-1) _ (rk-1, P-1rk-1)=(L-1A(LT)-1pk,Pk)(A(LT)-1pk, (LT)-1pk)(Apk,Pk)(rk, P-1rk)(L-Irk, L-1rk)k=(L-1rk-1, L-1rk-1)(rk-1, P-1rk-1)记zk会P-1rk，则可得下面的预处理CG方法.算法8.2.PCG：预处理CG方法1：给定初值(0)，(相对)精度要求ε>0和最大选代步数IterMax2:计算ro=b-Az(0)和β=roll23:令20=P1r0,P1=204计算p=z05: for k = 1 to IterMax do6:qk=Apk7:S=p/(pLq)2(k) = r(k-1) + EkPk8:9:rk=rk-1-Ekqk10:relres = Irkl/2 /β11:if relres< e thenbreak12:end if13:14:Po= p15:zk = P-1rk16:p=TIzk17:μk= p/po
8.1 预处理迭代方法 · 239 · 10: end for 由于算法 8.1 中得到的是预处理后的近似解 u (k) 和预处理残量 r˜k, 因此我们还需要考虑原方程的 近似解和残量的计算. 我们对上述算法中的迭代公式进行适当的改写. 首先引入一个辅助变量 pk : pk ≜ (L ⊺ ) −1 p˜k, k = 1, 2, . . . . 于是有 pk+1 = (L ⊺ ) −1 pk+1 = (L ⊺ ) −1 r˜k + µk(L ⊺ ) −1 p˜k = (L ⊺ ) −1 r˜k + µkpk. 又由 r˜k 的表达式可知, 原方程组的残量 rk = Lr˜k. 因此可得到下面的递推公式: rk = Lr˜k = L(˜rk−1 − ξkL −1A(L ⊺ ) −1 p˜k) = rk−1 − ξkApk, pk+1 = (L ⊺ ) −1 r˜k + µkpk = (L ⊺ ) −1L −1 rk + µkpk = P −1 rk + µkpk, x (k) = (L ⊺ ) −1u (k) = (L ⊺ ) −1 (u (k−1) + ξkp˜k) = x (k−1) + ξkpk, 其中 ξk = (L −1 rk−1, L−1 rk−1) (L−1A(L⊺)−1p˜k, p˜k) = (rk−1,(L ⊺ ) −1L −1 rk−1) (A(L⊺)−1p˜k,(L⊺)−1p˜k) = (rk−1, P −1 rk−1) (Apk, pk) , µk = (L −1 rk, L−1 rk) (L−1rk−1, L−1rk−1) = (rk, P −1 rk) (rk−1, P −1rk−1) . 记 zk ≜ P −1 rk, 则可得下面的预处理 CG 方法. 算法 8.2. PCG: 预处理 CG 方法 1: 给定初值 x (0), (相对) 精度要求 ε > 0 和最大迭代步数 IterMax 2: 计算 r0 = b − Ax(0) 和 β = ∥r0∥2 3: 令 z0 = P −1 r0, p1 = z0 4: 计算 ρ = r ⊺ 0 z0 5: for k = 1 to IterMax do 6: qk = Apk 7: ξk = ρ/(p ⊺ k qk) 8: x (k) = x (k−1) + ξkpk 9: rk = rk−1 − ξkqk 10: relres = ∥rk∥2/β 11: if relres< ε then 12: break 13: end if 14: ρ0 = ρ 15: zk = P −1 rk 16: ρ = r ⊺ k zk 17: µk = ρ/ρ0
.240.第八讲预处理18:Pk+1=Zk+μkpk19:end for20: if relres< then输出近似解（）及相关信息21:22:else23:输出算法失败信息24:endif我们注意到，矩阵L并没有出现在算法中，这意味着我们并不需要对P做Cholesky分解十在算法8.2中，每步选代都需要计算zk=P-1rk:一般情况下，都是通过求解方程组Pzk=r来实现，除非P-1非常容易得到，或者是直接给出的事实上，PCG方法8.2也可以从左预处理方式导出考虑左预处理方程组P-1Ar =P-1b.(8.6)易知P-1A是正定的，但通常不是对称的，因此我们不能直接对（8.6)实施CG方法.此时我们需要定义一个新的内积：P-内积，即(8.7)(r,y)p≤ (Pr,y)= (r, Py).由于P是对称正定的，所以这个定义是有效的.在该内积下，有(P-1Ar,y)p =(Ar,y) = (t, Ay) = (r, P(P-1Ay)) = (r, P-1Ay)p,即P-1A关于P-内积是自伴随的.也就是说，在P-内积意义下，P-1A是“对称”的.这时我们就可以用CG方法来求解预处理方程组（8.6),但需要将欧拉内积改为P-内积.引人辅助变量2k会P-1rk，则求解预处理方程组（8.6)的CG方法可描述如下：算法8.3.基于P-内积的CG方法1：给定初值(0),(相对)精度要求e>0和最大选代步数IterMax2:计算ro=b- Ar(0)和β=Iroll23:令20=P-1ro,P1=204:fork=1toIterMaxdo(rk-1, zk-1)(zk-1, zk-1)p5:Sk =(P-1 Apk, Pk)p(Apk,Pk)r(k) = r(k-1) + SxPk6:7:Tk = Tk-1 - EkApkrelres = rkl2/β8:9:if relres< then10:breakend if11:
· 240 · 第八讲 预处理 18: pk+1 = zk + µkpk 19: end for 20: if relres< ε then 21: 输出近似解 x (k) 及相关信息 22: else 23: 输出算法失败信息 24: end if † 我们注意到, 矩阵 L 并没有出现在算法中, 这意味着我们并不需要对 P 做 Cholesky 分解. † 在算法 8.2 中, 每步迭代都需要计算 zk = P −1 rk. 一般情况下, 都是通过求解方程组 P zk = rk 来 实现, 除非 P −1 非常容易得到, 或者是直接给出的. 事实上, PCG 方法 8.2 也可以从左预处理方式导出. 考虑左预处理方程组 P −1Ax = P −1 b. (8.6) 易知 P −1A 是正定的, 但通常不是对称的, 因此我们不能直接对 (8.6) 实施 CG 方法. 此时我们需要定义 一个新的内积: P-内积, 即 (x, y)P ≜ (P x, y) = (x, P y). (8.7) 由于 P 是对称正定的, 所以这个定义是有效的. 在该内积下, 有 (P −1Ax, y)P = (Ax, y) = (x, Ay) = (x, P(P −1Ay)) = (x, P −1Ay)P , 即 P −1A 关于 P-内积是自伴随的. 也就是说, 在 P-内积意义下, P −1A 是 “对称” 的. 这时我们就可以用 CG 方法来求解预处理方程组 (8.6), 但需要将欧拉内积改为 P-内积. 引入辅助变量 zk ≜ P −1 rk, 则求解预处理方程组 (8.6) 的 CG 方法可描述如下: 算法 8.3. 基于 P-内积的 CG 方法 1: 给定初值 x (0), (相对) 精度要求 ε > 0 和最大迭代步数 IterMax 2: 计算 r0 = b − Ax(0) 和 β = ∥r0∥2 3: 令 z0 = P −1 r0, p1 = z0 4: for k = 1 to IterMax do 5: ξk = (zk−1, zk−1)P (P −1Apk, pk)P = (rk−1, zk−1) (Apk, pk) 6: x (k) = x (k−1) + ξkpk 7: rk = rk−1 − ξkApk 8: relres = ∥rk∥2/β 9: if relres< ε then 10: break 11: end if
8.1预处理选代方法.241 12:zk= zk-1 - QkP-1Apk = P-1rk(zk, zk)p(rk, zk)13:μk :(Zk=1, Zk-1)p(rk-1, Zk-1)14:Pk+1 = Zk + μkPk15:end for十不难看出，算法8.3和算法8.2是完全一样的类似地，我们也可以从右预处理方式来推导PCG方法，具体过程留作练习8.1.2预处理GMRES方法对于非对称Krylov子空间方法，也有三种预处理方式.但与对称正定情形不同的是，这三种方式并不等价，而且有时效果会相差很大，左预处理方式首先考虑左预处理.设PERnxn是预处理子（非奇异），则左预处理方程组为P-1 Ar = P-1b.(8.8)记rk为原线性方程组的残量，即rk=b-A(k),则预处理方程组（8.8)的残量为=P-1rk将GMRES方法用于求解方程组（8.8），即可得下面的左预处理GMRES方法算法8.4.左预处理GMRES方法1:给定初值(0)和（相对)精度要求>02:计算r=P-1(6-Az(0))和β=/oll23:令=To/B,=ei4:forj=1,2....do5:i, = Auj6:w, = P-1ij% Apply preconditioning7:fori=1,2,....jdo8:hij = (wj,u)9:Wi=wj-hijviendfor10:11:hj+1, = wjll212:for i =1, 2,...,j -1 do% Apply Gj-1,...,Gi to thelast column of Hi+1.jhi.jCiSihij13:hi+1ihi+lend for14:15:if hj+1.j =0 then
8.1 预处理迭代方法 · 241 · 12: zk = zk−1 − αkP −1Apk = P −1 rk 13: µk = (zk, zk)P (zk−1, zk−1)P = (rk, zk) (rk−1, zk−1) 14: pk+1 = zk + µkpk 15: end for † 不难看出, 算法 8.3 和算法 8.2 是完全一样的. † 类似地, 我们也可以从右预处理方式来推导 PCG 方法, 具体过程留作练习. 8.1.2 预处理 GMRES 方法 对于非对称 Krylov 子空间方法, 也有三种预处理方式. 但与对称正定情形不同的是, 这三种方式并 不等价, 而且有时效果会相差很大. 左预处理方式 首先考虑左预处理. 设 P ∈ R n×n 是预处理子 (非奇异), 则左预处理方程组为 P −1Ax = P −1 b. (8.8) 记 rk 为原线性方程组的残量, 即 rk = b − Ax(k) , 则预处理方程组 (8.8) 的残量为 r˜k = P −1 rk. 将 GMRES 方法用于求解方程组 (8.8), 即可得下面的左预处理 GMRES 方法. 算法 8.4. 左预处理 GMRES 方法 1: 给定初值 x (0) 和 (相对) 精度要求 ε > 0 2: 计算 r˜0 = P −1 (b − Ax(0)) 和 β = ∥r˜0∥2 3: 令 v1 = ˜r0/β, ξ = βe1 4: for j = 1, 2, . . . do 5: w˜j = Avj 6: wj = P −1w˜j % Apply preconditioning 7: for i = 1, 2, . . . , j do 8: hij = (wj , vi) 9: wj = wj − hijvi 10: end for 11: hj+1,j = ∥wj∥2 12: for i = 1, 2, . . . , j − 1 do % Apply Gj−1, . . . , G1 to the last column of Hj+1,j 13: [ hi,j hi+1,j] = [ ci si −si ci ] [ hi,j hi+1,j] 14: end for 15: if hj+1,j = 0 then