共轭梯度法的预处理

1、共轭梯度法的预处理共轭梯度法的预处理设 A 为 n 阶对称正定矩阵, 考虑线性方程组(1)其中x 是未知向量,b是右端已知向量。Axb一、知识回顾一、知识回顾1、共轭梯度法 如果我们对于固定的 ，在 空间中取 个线性无关的向量 使 满足 当 时， 一定是方程 (1)的精确解，这是共轭梯度法的出发点。knRk12,kp ppkx12(,)(x )min (x)kkSpan pppknnx一、知识回顾一、知识回顾2、预处理共轭梯度法（1）预处理共轭梯度法的引进 共轭梯度法理论上是一种直接方法，即按照这个算法得到的 应当是方程 的精确解。但由于其数值的稳定性不好，实际上由于舍入误差的影响 不可能满足

2、原方程。如果把它看成一种迭代法，当A条件数较大时，它收敛的很慢。因此引进了改进共轭梯度法的手段预处理共轭梯度法。nxAxbnx一、知识回顾一、知识回顾（2）预处理共轭梯度法原理 将Ax=b变形为C-1Ax=C-1b。 令 则 要求 当 是对称正定矩阵时， 其中 是A的最大特征值， 是A的最小特征值，K(A)是A的条件数 AxbAxb22( )( )condAcondAA12( )( )ncondAK A1nbCbxCxACCATT11,一、知识回顾一、知识回顾（3）预处理共轭梯度方法由于预处理矩阵的构造方式很多，所以就形成了多种预处理方法，总的来说，大致可归为如下几种途径：1、取预处理矩阵（也

3、称预优矩阵）为A的一个小带宽部分（如三对角或对角线部分）；2、通过A的分裂，尤其是线性稳定迭代法中的矩阵 A的分裂构造预处理矩阵（如矩阵分裂法）；3、通过A的各种近似分解得到预处理矩阵（如不完全分解）；4、通过矩阵A的多项式构造预处理矩阵；5、子结构、区域分裂、EBE预处理途径。一、知识回顾一、知识回顾对角预优矩阵法如果对称正定线性方程组( 1) 的系数矩阵 A 的对角元素相差较大, 则可取预优矩阵C-1为这是最常用的预处理方阵nnijnnaAaaadiagDDC)(),(,22112/11SSOR分裂将A分裂为 ，其中 , 为严格下三角矩阵，则预处理矩阵C为：其中 是参数，而且有TLLCCD

4、A1122(,)nnDdiag aaa)10(TLLTCDDCDCCA)()()2(11LC2/1)()2(1DCDCL二、算法二、算法（1）共轭梯度法算法 令 对 直到收敛完成 结束。 000000,xrbAxpr0,1,2,k 1111111(r ,r )/(p ,Ap ),(r ,r ),(r,r),kkkkkkkkkkkkkkkkkkkkkkxxprrApprp二、算法二、算法（2）预处理共轭梯度算法 令 对 直到收敛完成。 结束。010000000,xrbAx zM r pz0,1,k 1111111,1111(r ,r )/(p ,Ap ),(zz)/(r ,z ),kkkkkkk

5、kkkkkkkkkkkkkkkkkxxprrApzM rpzp三、代码三、代码（1）共轭梯度法代码function x,n,tol=cg(A,b,ep)r=b;p=r;n=0;x=zeros(length(b),1);while n10 alpha=(r*r)/(p*A*p);r1=r; x=x+alpha*p; tol=norm(r); if norm(r)ep break; end r=r-alpha*A*p; beta=(r*r)/(r1*r1); p=r+beta*p; n=n+1;endx = 7.859713075445870.422926408294999-0.073592239

6、0236594-0.5406430168946320.0106261628540375n = 5tol =6.745229690695623e-07（2）预处理共轭梯度法代码以对角预优矩阵为例clc;clearA=0.2 0.1 1 1 0;0.1 4 -1 1 -1;1 -1 60 0 -2;1 1 0 8 4;0 -1 -2 4 700;b=1 2 3 4 5;%SSOR此处程序改为C=chuli1(A); A1=C*A*C; b1=C*b;ep=10(-2);x,n,tol=cg1(A1,b1,ep,C);function C=chuli1(A) D=diag(diag(A); C=D(

7、-1/2);function C=chuli(A,w) D=diag(diag(A);L=tril(A,-1); C=(D-w*L)*D(-1/2);C=chuli(A,0.009);A1=inv(C)*A*inv(C);b1=inv(C)*b;function x,n,tol=cg1(A,b,ep,C)r=b;p=r;n=0;x=zeros(length(b),1);while n10 alpha=(r*r)/(p*A*p);r1=r; x=x+alpha*C*p; tol=norm(r); if norm(r)ep break; end r=r-alpha*A*p; beta=(r*r)/

8、(r1*r1); p=r+beta*p; n=n+1;end%SSOR此处程序改为x=x+alpha*inv(C)*pSSOR分裂x = 7.859713075445870.422926408295008-0.0735922390240462-0.5406430168946260.0106261628540363m= 4tol = 6.567520686375984e-05对角预优矩阵x= 7.859713075445860.422926408295008-0.0735922390240463-0.5406430168946260.0106261628540363m= 4tol=4.73175

9、3657562764e-04方法方法迭代迭代次数次数xk误差误差共轭梯度法5 7.85971307544587,0.422926408294999-0.0735922390236594,-0.5406430168946320.01062616285403756.745229690695623e-07对角预优矩阵47.859713075445865;0.422926408295008;-0.073592239024046;-0.540643016894626;0.0106261628540366.567520686375984e-05SSOR分裂47.859713075445860;0.4229

10、26408295008;-0.073592239024046;-0.540643016894626;0.0106261628540364.731753657562764e-04此矩阵并没有显示出预处理的优势，所以我们又取矩阵A=diag(1 10 100 1000 10000)，输出结果如下方法方法迭代迭代次数次数xk误差误差共轭梯度法51;0.200000000000000;0.030000000000001;0.003999999999562;5.000000000000054e-043.895966164939117e-06对角预优矩阵11;0.200000000000000;0.030000000000000;0.004000000000000;5.000000000000000e-042.775557561562891e-17SSOR分裂11;0.200000000000000;0.030000000000000;0.004000000000000;5.000000000000001e-041.2688