第八讲预处理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 relres0和最大选代步数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 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
.242.第八讲预处理16:set m= j and breakendif17:18:Uj+1 = wj/hj+1,j19:if [hj.5l>[hi+1,|then%Form theGivens rotationG,hi+1.20:Vi+Sj=CjTwhereT =Ci =-hj.s21:elsehi.s22:Sj= V Cj= SjT where T hi+1.end if23:24:hjj=cjhjj+sjhj+1,j,hi+1j=0% Apply G, to last column of Hj+1.j5]-[g5% Apply G, to the right-hand side25:CE+1sici26:ifSi+il/β<then% Check convergence: I/,ll2 = [6j+1l27:set m = j and breakendif28:29:endfor30: 计算y(m)= R=1s(m),其中 Rm =H(1 : m,1:m),ε(m) = E(1 : m)31:计算近似解(m)=z(0)+Vmy(m)需要指出的是，在左预处理GMRES方法中，停机准则中的ISi+1I并不是原方程组残量r的模，而是预处理残量=P-1(b－Az()）的模，如果要计算真实残量，除了显式计算外，没有其它更好的办法.所以，如果停机准则是基于真实残量的，则需要额外去计算，这就会增加很多额外的运算量.幸运的是，右预处理方式没有这个问题右预处理方式如果用P做右预处理，则对应的预处理方程组为AP-lu=b,r=P-lu(8.9)给定选代初值(0)，则u(0)=P(0).相应的预处理初始残量为ro=b-AP-1u(0)=b- Az(0) = r(0)即预处理残量与原方程组残量是一样的.因此无需计算u(0).设u（m）=u(0）+Vmm是迭代m步后的近似解，则对应的原方程组近似解(m）为r(m) = P-1u(m) = P-1(u(0) + Vmym) = r(0) + P-1Vmym:原方程组的残量为rm=b- Ar(m)= ro- AP-1Vmym=βui -Vm+1Hm+1,mym=Vm+1(Bei -Hm+1,mJm).由于m是最小二乘问题min Ilβei-Hm+1,myll2VER
· 242 · 第八讲 预处理 16: set m = j and break 17: end if 18: vj+1 = wj/hj+1,j 19: if |hj,j | > |hj+1,j | then % Form the Givens rotation Gj 20: cj = √ 1 1+τ 2 , sj = cj τ where τ = hj+1,j hj,j 21: else 22: sj = √ 1 1+τ 2 , cj = sj τ where τ = hj,j hj+1,j 23: end if 24: hj,j = cjhj,j + sjhj+1,j , hj+1,j = 0 % Apply Gj to last column of Hj+1,j 25: [ ξj ξj+1] = [ cj sj −sj cj ] [ξj 0 ] % Apply Gj to the right-hand side 26: if |ξj+1|/β < ε then % Check convergence: ∥r˜j∥2 = |ξj+1| 27: set m = j and break 28: end if 29: end for 30: 计算 y (m) = R−1 m ξ (m) , 其中 Rm = H(1 : m, 1 : m), ξ (m) = ξ(1 : m) 31: 计算近似解 x (m) = x (0) + Vmy (m) 需要指出的是, 在左预处理 GMRES 方法中, 停机准则中的 |ξj+1| 并不是原方程组残量 rk 的模, 而 是预处理残量 r˜k = P −1 (b − Ax(k) ) 的模. 如果要计算真实残量, 除了显式计算外, 没有其它更好的办 法. 所以, 如果停机准则是基于真实残量的, 则需要额外去计算 rk, 这就会增加很多额外的运算量. 幸运 的是, 右预处理方式没有这个问题. 右预处理方式 如果用 P 做右预处理, 则对应的预处理方程组为 AP −1u = b, x = P −1u. (8.9) 给定迭代初值 x (0) , 则 u (0) = P x(0) . 相应的预处理初始残量为 r˜0 = b − AP −1u (0) = b − Ax(0) = r (0) , 即预处理残量与原方程组残量是一样的. 因此无需计算 u (0) . 设 u (m) = u (0) + Vmym 是迭代 m 步后的 近似解, 则对应的原方程组近似解 x (m) 为 x (m) = P −1u (m) = P −1 (u (0) + Vmym) = x (0) + P −1Vmym, 原方程组的残量为 rm = b − Ax(m) = r0 − AP −1Vmym = βv1 − Vm+1Hm+1,mym = Vm+1(βe1 − Hm+1,mym). 由于 ym 是最小二乘问题 min y∈Rm ∥βe1 − Hm+1,my∥2
8.1预处理选代方法.243-的解、因此IIrml2=lBei-Hm+1,mymll2=β-lqi(m+1)l这与不带预处理的GMRES方法是一样的.因此在实际求解过程中我们可以直接得到原方程组残量的模，而无需计算u（m）和（m）.这是右预处理方式与左预处理方式的主要区别右预处理GMRES方法具体描述如下：算法8.5.右预处理GMRES方法1：给定初值(0）和（相对）精度要求e>02:计算ro=b-Az(0)和β=llroll23: 令 Ui = Fo/B, E = βe14:for j =1,2....do5:i,= P-1uj% Apply preconditioning6:wj=Awjfori= 1, 2.....j do7:8:hij =(wj,ui)9:Wj = wj - hijuiend for10:11:hj+1,j = Iwjll212:fori=1,2,....j-1do% ApplyG-1....,Gi to thelast column ofHi+1.hijhi.j13:hi+1.jhi+1end for14:if hij+1.j=0 then15:set m = j and break16:endif17:18:Uj+1=Wj/hj+1.j19:if [hj.jl>|hi+1.| then%Form theGivens rotation G,h$+l.120:Cj = V+ j = CT where T =hj.selse21:his22:I+,Cj=SjTwhereTsi=hi+1.j23:endif24:%ApplyG,to lastcolumn ofHi+1.5hj.j=cjhj,j+ sjhi+1.j,hi+1.j=0SiSCiF25:% Apply G, to the right-hand side0Si+1Si26:% Check convergence: Ir,ll2 = 6i+1if Sj+1l/β<then27:set m = j and break28:endif29:endfor
8.1 预处理迭代方法 · 243 · 的解, 因此 ∥rm∥2 = ∥βe1 − Hm+1,mym∥2 = β · |q1(m + 1)|. 这与不带预处理的 GMRES 方法是一样的. 因此在实际求解过程中我们可以直接得到原方程组残量的 模, 而无需计算 u (m) 和 x (m) . 这是右预处理方式与左预处理方式的主要区别. 右预处理 GMRES 方法具体描述如下: 算法 8.5. 右预处理 GMRES 方法 1: 给定初值 x (0) 和 (相对) 精度要求 ε > 0 2: 计算 r0 = b − Ax(0) 和 β = ∥r0∥2 3: 令 v1 = ˜r0/β, ξ = βe1 4: for j = 1, 2, . . . do 5: w˜j = P −1vj % Apply preconditioning 6: wj = Aw˜j 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 16: set m = j and break 17: end if 18: vj+1 = wj/hj+1,j 19: if |hj,j | > |hj+1,j | then % Form the Givens rotation Gj 20: cj = √ 1 1+τ 2 , sj = cj τ where τ = hj+1,j hj,j 21: else 22: sj = √ 1 1+τ 2 , cj = sj τ where τ = hj,j hj+1,j 23: end if 24: hj,j = cjhj,j + sjhj+1,j , hj+1,j = 0 % Apply Gj to last column of Hj+1,j 25: [ ξj ξj+1] = [ cj sj −sj cj ] [ξj 0 ] % Apply Gj to the right-hand side 26: if |ξj+1|/β < ε then % Check convergence: ∥rj∥2 = |ξj+1| 27: set m = j and break 28: end if 29: end for
.244.第八讲预处理30:计算y(m)=R=1(m),其中Rm=H(1:m,1:m),(m)=(1:m)31:计算近似解r(m)=r(0)+P-1Vmy(m)两边预处理方式如果预处理子P是由乘积形式给出的，即P=PLPR,则可以构造下面的两边预处理方程组P-"APr"u= Pr"b, = Pr"u.将GMRES用于求解上述线性方程组，则可得到两边预处理GMRES方法.具体推导过程与前面类似，留作练习需要指出的是，对于两边预处理GMRES方法，其预处理后的残量为=L-1(6-Ar(k))十由于左预处理方式和两边预处理方式所得到的残量并不是原方程组的残量，因此它们都可能会导致选代提前终止或延迟终止，特别是当预处理子P的条件数较大时.因此，通常情况下建议采用右预处理方式8.1.3左、右预处理GMRES的最优性质设(0）是选代初始值，记（）和）分别是左预处理GMRES方法和右预处理GMRES方法迭代步后得到的近似解。根据GMRES方法的最优性质可知，对于左预处理GMRES方法，有2(2) =P-1(6 - Ar)12(8.10)argminEr(0)+K(P-1A,P--1ro而对于右预处理GMRES,我们有(k))=P-1u(A),其中u(3) =[6 - AP-1ul2argminuEu(0)+K(AP-1,ro)通过变量代换=P-1u,可得)=6- Arl2.(8.11)argminEr(0)+P-1K(AP=1,ro)又P-1Kk(AP-1,ro)=P-1 span(ro,AP-1ro, (AP-1)?ro,...,(AP-1)k-1ro)=span[P-1ro, P-1AP-1ro,P-1(AP-1)?ro,.**,P-1(AP-1)k-1ro)= span[P-1ro, (P-1A)P-1ro, (P-1A)?p-1ro,.... (P-1A)k-1 p-1ro)= K(P-1A, P-1ro),由（8.10)和（8.11)可知，(k）和（(）属于同一个子空间中.不同之处在于，(k）极小化的是预处理后残量的范数,即 IIP-1(b-Az)l2,而()极小化的是原方程组残量的范数，即Ib-Aazl2
· 244 · 第八讲 预处理 30: 计算 y (m) = R−1 m ξ (m) , 其中 Rm = H(1 : m, 1 : m), ξ (m) = ξ(1 : m) 31: 计算近似解 x (m) = x (0) + P −1Vmy (m) 两边预处理方式 如果预处理子 P 是由乘积形式给出的, 即 P = PLPR, 则可以构造下面的两边预处理方程组 P −1 L AP −1 R u = P −1 L b, x = P −1 R u. 将 GMRES 用于求解上述线性方程组, 则可得到两边预处理 GMRES 方法. 具体推导过程与前面类似, 留 作练习. 需要指出的是, 对于两边预处理 GMRES 方法, 其预处理后的残量为 r˜k = L −1 (b − Ax(k) ). † 由于左预处理方式和两边预处理方式所得到的残量并不是原方程组的残量, 因此它们都可能会 导致迭代提前终止或延迟终止., 特别是当预处理子 P 的条件数较大时. 因此, 通常情况下建议采 用右预处理方式. 8.1.3 左、右预处理 GMRES 的最优性质 设 x (0) 是迭代初始值, 记 x (k) L 和 x (k) R 分别是左预处理 GMRES 方法和右预处理 GMRES 方法迭代 k 步后得到的近似解. 根据 GMRES 方法的最优性质可知, 对于左预处理 GMRES 方法, 有 x (k) L = argmin x∈x(0)+Kk(P −1A,P −1r0) ∥P −1 (b − Ax)∥2. (8.10) 而对于右预处理 GMRES, 我们有 x (k) R = P −1u (k) , 其中 u (k) = argmin u∈u(0)+Kk(AP −1,r0) ∥b − AP −1u∥2. 通过变量代换 x = P −1u, 可得 x (k) R = argmin x∈x(0)+P −1Kk(AP −1,r0) ∥b − Ax∥2. (8.11) 又 P −1Kk(AP −1 , r0) = P −1 span{r0, AP −1 r0,(AP −1 ) 2 r0, . . . ,(AP −1 ) k−1 r0} = span{P −1 r0, P −1AP −1 r0, P −1 (AP −1 ) 2 r0, . . . , P −1 (AP −1 ) k−1 r0} = span{P −1 r0,(P −1A)P −1 r0,(P −1A) 2P −1 r0, . . . ,(P −1A) k−1P −1 r0} = Kk(P −1A, P −1 r0), 由 (8.10) 和 (8.11) 可知, x (k) L 和 x (k) R 属于同一个子空间中. 不同之处在于, x (k) L 极小化的是预处理后残 量的范数, 即 ∥P −1 (b − Ax)∥2, 而 x (k) R 极小化的是原方程组残量的范数, 即 ∥b − Ax∥2
8.2预处理技术.245 定理8.1设(）和）分别是左预处理GMRES方法和右预处理GMRES方法选代k步后得到的近似解,且选代初始值均为(0),则r(k)是P-1(b-Az)l2在(0)+K(P-1A,P-1ro)的极小值点,而r(k）则是6-Arl2在同一个子空间中的极小值点R8.2预处理技术Finding a good preconditioner to solve a given sparse linear system is often viewed as a combinationof art and science.Theoretical results are rare and some methods work surprisingly well, often despiteexpectations.—Saad[64],2003预处理能否取得成功的关键是能否找到一个好的预处理子。预处理子的构造与问题本身是密切相关的，通用的“好预处理子”是不存在的，关于预处理技术的理论分析很少，大多数情况下只能根据经验来构造，尽管如此，在实际应用中，这些根据经验构造出来的预处理子往往能取得好好的效果，有时甚至会大大出乎人们的意外一般来说，预处理子可以分为两大类(a）代数预处理子（AlgebraicPreconditioner)，即仅仅根据所给的矩阵来构造预处理子(b）专用预处理子(Problem-SpecificPreconditioner)，即根据问题的物理背景所构造的预处理子显然，由于专用预处理子充分利用了问题的物理背景知识，所以它们往往具有很好的数值表现，如多重网格，区域分解，快速变换等等，但它们严重依赖于原问题的物理背景，因此不具有通用性我们这里只介绍代数预处理子，即仅仅根据所给的系数矩阵来构造预处理方法.这种预处理方法具有一定的通用性预处理子的选取通常与具体问题本身是密切有关的，一般的选取方法有·设有矩阵分裂A=M-N则M可作为一个预处理子；·不完全LU分解；·近似逆，即选取矩阵P，使得P-1～A-1·对角矩阵，块对角矩阵，三角矩阵，块三角矩阵等等8.2.1矩阵分裂预处理子考虑线性方程组Ar = b,AERnxn对A做如下的矩阵分裂A-M-N(8.12)其中M非奇异，则可以得到下面的选代方法(k+1) = M-1N(k) + M-1b = 2(k) + (M-1b - M-1Ar(k) . k = 0,1,2
8.2 预处理技术 · 245 · 定理 8.1 设 x (k) L 和 x (k) R 分别是左预处理 GMRES 方法和右预处理 GMRES 方法迭代 k 步后得到的近 似解, 且迭代初始值均为 x (0) . 则 x (k) L 是 ∥P −1 (b − Ax)∥2 在 x (0) + Kk(P −1A, P −1 r0) 的极小值点, 而 x (k) R 则是 ∥b − Ax∥2 在同一个子空间中的极小值点. 8.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 [64], 2003. 预处理能否取得成功的关键是能否找到一个好的预处理子. 预处理子的构造与问题本身是密切相 关的, 通用的 “好预处理子” 是不存在的. 关于预处理技术的理论分析很少, 大多数情况下只能根据经验来构造. 尽管如此, 在实际应用中, 这 些根据经验构造出来的预处理子往往能取得好好的效果, 有时甚至会大大出乎人们的意外. 一般来说, 预处理子可以分为两大类 (a) 代数预处理子 (Algebraic Preconditioner), 即仅仅根据所给的矩阵来构造预处理子. (b) 专用预处理子 (Problem-Specific Preconditioner), 即根据问题的物理背景所构造的预处理子. 显然, 由于专用预处理子充分利用了问题的物理背景知识, 所以它们往往具有很好的数值表现, 如多重 网格, 区域分解, 快速变换等等. 但它们严重依赖于原问题的物理背景, 因此不具有通用性. 我们这里只介绍代数预处理子, 即仅仅根据所给的系数矩阵来构造预处理方法. 这种预处理方法具 有一定的通用性. 预处理子的选取通常与具体问题本身是密切有关的. 一般的选取方法有 • 设有矩阵分裂 A = M − N, 则 M 可作为一个预处理子; • 不完全 LU 分解; • 近似逆, 即选取矩阵 P, 使得 P −1 ≈ A−1 ; • 对角矩阵, 块对角矩阵, 三角矩阵, 块三角矩阵等等. 8.2.1 矩阵分裂预处理子 考虑线性方程组 Ax = b, A ∈ R n×n , 对 A 做如下的矩阵分裂: A = M − N (8.12) 其中 M 非奇异, 则可以得到下面的迭代方法 x (k+1) = M−1Nx(k) + M−1 b = x (k) + ( M−1 b − M−1Ax(k) ) . k = 0, 1, 2, . . .
:246.第八讲预处理这等价于求解下面的方程组M-1Ar = M-1b.(8.13)这就是与矩阵分裂（8.12）相对应的预处理线性方程组.将Krylov子空间方法用于求解方程组（8.13)，就得到预处理Krvlov子空间方法.矩阵M就是由矩阵分裂（8.12）所定义的预处理子理论上讲，任何一个矩阵分裂都可以定义一个预处理子.但为了使得预处理子能有很好的预处理效果，往往需要其在一定意义下与A充分接近设A=D-L-U,其中D,-L,-U分别是A的对角部分，严格下三角部分和严格上三角部分,并假定D非奇异则由我们之前讨论的定常迭代法，可以立即得到下面的预处理子·Jacobi预处理子，即取A的对角部分作为预处理子：Pj = D.·G-S预处理子，即取A的下三角部分作为预处理子：PGs =D-L.·SOR预处理子，即(wD - L).PsOR=·SSOR预处理子，即Pso = u(2- (1-) + wl] -[1 - u) + u] .由于SSOR对参数w的取值不是很敏感，因此我们通常令w=1，即PsSOR = (D - L)D-1(D - U).8.2.2不完全LU分解One of the most popular preconditioners is the class of incomplete factorization preconditioners which areusually effective for large and sparse matrices. The incomplete factorization methods were first introduced by Buleev[13, 14] in the late 1950s, and independently by Varga [77]; see also [45]. However, a major breakthrough tookplace in 1977, by Meijerink and van der Vorst [54] where the incomplete Cholesky-conjugate gradient (ICCG)algorithm was proposed. Since then, a number of improvements and extensions have been made, including level offll and drop tolerance-based incomplete factorizations, and so on.When the LU decomposition of a sparse matrix A is carried out through Gaussian elimination, fill-in usuallytakes place.This means that the triangular factors L and U are considerably less sparse than A.However, bydropping part of the fill-in in the processing of the factorization, we can obtain a simple but useful preconditionerof theform P=LU,where L and U arethe incompleteLUfactors.In general, the incomplete LU factorization (ILU) is a process that computes a sparse unit lower-triangularmatrix L and a sparse upper-triangular matrix U such thatA=LU-R
· 246 · 第八讲 预处理 这等价于求解下面的方程组 M−1Ax = M−1 b. (8.13) 这就是与矩阵分裂 (8.12) 相对应的预处理线性方程组. 将 Krylov 子空间方法用于求解方程组 (8.13), 就 得到预处理 Krylov 子空间方法. 矩阵 M 就是由矩阵分裂 (8.12) 所定义的预处理子. 理论上讲, 任何一个矩阵分裂都可以定义一个预处理子. 但为了使得预处理子能有很好的预处理效 果, 往往需要其在一定意义下与 A 充分接近. 设 A = D − L − U, 其中 D, −L, −U 分别是 A 的对角部分, 严格下三角部分和严格上三角部分, 并 假定 D 非奇异. 则由我们之前讨论的定常迭代法, 可以立即得到下面的预处理子: • Jacobi 预处理子, 即取 A 的对角部分作为预处理子: PJ = D. • G-S 预处理子, 即取 A 的下三角部分作为预处理子: PGS = D − L. • SOR 预处理子, 即 PSOR = 1 ω (ωD − L). • SSOR 预处理子, 即 PSSOR = 1 ω(2 − ω) [ (1 − ω)D + ωL] D−1 [ (1 − ω)D + ωU] . 由于 SSOR 对参数 ω 的取值不是很敏感, 因此我们通常令 ω = 1, 即 PSSOR = (D − L)D−1 (D − U). 8.2.2 不完全 LU 分解 One of the most popular preconditioners is the class of incomplete factorization preconditioners which are usually effective for large and sparse matrices. The incomplete factorization methods were first introduced by Buleev [13, 14] in the late 1950s, and independently by Varga [77]; see also [45]. However, a major breakthrough took place in 1977, by Meijerink and van der Vorst [54] where the incomplete Cholesky-conjugate gradient (ICCG) algorithm was proposed. Since then, a number of improvements and extensions have been made, including level of fill and drop tolerance-based incomplete factorizations, and so on. When the LU decomposition of a sparse matrix A is carried out through Gaussian elimination, fill-in usually takes place. This means that the triangular factors L and U are considerably less sparse than A. However, by dropping part of the fill-in in the processing of the factorization, we can obtain a simple but useful preconditioner of the form P = LU, where L and U are the incomplete LU factors. In general, the incomplete LU factorization (ILU ) is a process that computes a sparse unit lower-triangular matrix L and a sparse upper-triangular matrix U such that A = LU − R