5-1-1 预处理CG方法 设A∈Rn×n对称正定，要求预处理子P也对称正定. 考虑使用两边预处理方式：设 P=LLT L-1A(L)-1u=L-1b,x=(D)-1u 巴用CG方法求解，迭代k步后的近似解记为网 →因=(L)-1因 图记kL-1b-L-1A(LT)-1u(周 →预处理后的残量 http://sath.ecnu.edu.cn/-jypan 16/72
511 预处理 CG 方法 设 A ∈ R n×n 对称正定, 要求预处理子 P 也对称正定. 考虑使用 两边预处理方式: 设 P = LL⊺ L −1A(L ⊺ ) −1 u = L −1 b, x = (L ⊺ ) −1 u 用 CG 方法求解, 迭代 k 步后的近似解记为 u (k) → x (k) = (L ⊺ ) −1u (k) 记 ˜rk ≜ L −1 b − L −1A(L ⊺ ) −1u (k) → 预处理后的残量 http://math.ecnu.edu.cn/~jypan 16/72
CG基本框架：Ax=b PCG:L-1AL-Tu=L-1b 1:给定初值xo) 1:给定初值o) 2:计算0=b-Axo) 2:计算0=b-Ax0) 3:p1=T0 3:0=L-10,1=70 4:fork=1,2,. do 4:for k =1,2,...do (Tk-1,Tk-1)】 5: 5k= (Tk-1,Tk-1) 5: (Apk,Pk) (亿-1AL-Tk,) 6: x因=k-1)+kpk 6: 因=k-1)+kPk 7: Tk Tk-1 EkApk 7 Tk =Tk-1-EkL-1AL-Tp (,i) 8: (Tk;Tk) Uk= 8: k= (Tk-1,Tk-1 (k-1,ik-1） 9: Pk+1=Tk十HkPk 9: k+1=ik十kpk 10:end for 10:end for 文算法得到的是因和，因此需要考虑原方程的近似解和残量 http://math.ecnu.edu.cn/-jypan 17/72
CG 基本框架: Ax = b 1: 给定初值 x (0) 2: 计算 r0 = b − Ax(0) 3: p1 = r0 4: for k = 1, 2, . . . do 5: ξk = (rk−1, rk−1) (Apk, pk) 6: x (k) = x (k−1) + ξkpk 7: rk = rk−1 − ξkApk 8: µk = (rk, rk) (rk−1, rk−1) 9: pk+1 = rk + µkpk 10: end for PCG: L −1AL−⊺ u = L −1 b 1: 给定初值 x (0) 2: 计算 r0 = b − Ax(0) 3: ˜r0 = L −1 r0, p˜1 = ˜r0 4: for k = 1, 2, . . . do 5: ξk = (˜rk−1, ˜rk−1) (L−1AL−⊺ p˜k, p˜k) 6: u (k) = u (k−1) + ξkp˜k 7: ˜rk = ˜rk−1 − ξkL −1AL−⊺ p˜k 8: µk = (˜rk, ˜rk) (˜rk−1, ˜rk−1) 9: p˜k+1 = ˜rk + µkp˜k 10: end for 算法得到的是 u (k) 和 ˜rk, 因此需要考虑原方程的近似解和残量. http://math.ecnu.edu.cn/~jypan 17/72
原方程的解x因=Lt因和残量=L7 留对迭代公式进行适当的改写，引入一个辅助变量 p合L-Tk →pk+1=LTpk+1=L-TTk十kLTpR=L-Tk十kPk Tk=Lik=L(Tk-1-EkLAL TPR)=Tk-1-EkAPk Pkt1=L-Tie+Hkpk L-TLIrk+ukpk=P-lrk+Hkpk: x因=L-T因=L-T(k-1)+k)=xk-1)+kP%, i-i-i)=(-1nk-1,I-1nk-=（k-1,-T-1k-=（k-1,P-1nk- A=户A-,=(L-AL-Tk,)=(AL-,(亿-T】 (Apk:Pk) (Tk:i)(Lrk:Lir)(rk:P-ire) k三 (k-1,k-1(L-1k-1,L-1k-1)(k-1,P-1-1 (记kPr1r) http://aath.ecnu.edu.cn/-jypan 18/72
原方程的解 x (k) = L −⊺ u (k) 和残量 rk = L˜rk 对迭代公式进行适当的改写, 引入一个辅助变量 pk ≜ L −⊺ p˜k =⇒ pk+1 = L −⊺ pk+1 = L −⊺ ˜rk + µkL −⊺ p˜k = L −⊺ ˜rk + µkpk rk = L˜rk = L(˜rk−1 − ξkL −1AL−⊺ p˜k) = rk−1 − ξkApk, pk+1 = L −⊺ ˜rk + µkpk = L −⊺ L −1 rk + µkpk = P −1 rk + µkpk, x (k) = L −⊺ u (k) = L −⊺ (u (k−1) + ξkp˜k) = x (k−1) + ξkpk, ξk = (˜rk−1, ˜rk−1) (L−1AL−⊺ p˜k, p˜k) = (L −1 rk−1, L −1 rk−1) (L−1AL−⊺ p˜k, p˜k) = (rk−1, L −⊺L −1 rk−1) (AL−⊺ p˜k,(L−⊺ p˜k) = (rk−1,P −1 rk−1) (Apk, pk) µk = (˜rk, ˜rk) (˜rk−1, ˜rk−1) = (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) http://math.ecnu.edu.cn/~jypan 18/72
算法5.1预处理CG(PCG)算法 1:给定初值to),(相对)精度要求e>0和最大迭代步数IterMax 2:计算m=b-Ao)和B=川ml2 3:20=P10,p1=20 4:计算p=0 5:for k=1 to IterMax do 6: qk=Apk,Ek=p/(P限qk) 7: x因=k-1)+EkPk 8: Tk Tk-1 -Ekqk 9: relres =llrall2/B 10: If relres<e then break end if%收敛性检测 11: Po=P 12: 2k P-1Tk; p=以 13: Hk=P/Po 14: Pk+1=2十4kPk 15:end for http://math.ecnu.edu.cn/-jypan 19/72
算法 5.1 预处理 CG (PCG) 算法 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, ξk = ρ/(p ⊺ k qk) 7: x (k) = x (k−1) + ξkpk 8: rk = rk−1 − ξkqk 9: relres = ∥rk∥2/β 10: If relres< ε then break end if % 收敛性检测 11: ρ0 = ρ 12: zk = P −1 rk, ρ = r ⊺ k zk 13: µk = ρ/ρ0 14: pk+1 = zk + µkpk 15: end for http://math.ecnu.edu.cn/~jypan 19/72
几点说明 ·矩阵L并没有出现在算法中，这意味着并不需要对P做Cholesky分解 在PCG算法中，每步迭代都需要计算4=P1,一般情况下，都是通过 求解方程组Pk=Tk来实现. http://aath.ecnu.edu.cn/-jypan 20/72
几点说明 ▶ 矩阵 L 并没有出现在算法中, 这意味着并不需要对 P 做 Cholesky 分解 ▶ 在 PCG 算法中, 每步迭代都需要计算 zk = P −1 rk, 一般情况下, 都是通过 求解方程组 Pzk = rk 来实现. http://math.ecnu.edu.cn/~jypan 20/72