共轭梯度法的预处理
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
这是最常用的预处理方阵
SSOR分裂
将A分裂为 A D CL CLT，其中D diag(a11, a22,ann) , CL为严格
下三角矩阵，则预处理矩阵C为：
C
1
(D CL)D1/ 2
(2 )
其中(0 1) 是参数，而且有
A CCT
1
(D CL)D1(D CL)T
(2 )
一、知识回顾
（3）预处理共轭梯度方法
由于预处理矩阵的构造方式很多，所以就形成了多种预处理 方法，总的来说，大致可归为如下几种途径：
1、取预处理矩阵（也称预优矩阵）为A的一个小带宽部分 （如三对角或对角线部分）；
2、通过A的分裂，尤其是线性稳定迭代法中的矩阵 A的分裂 构造预处理矩阵（如矩阵分裂法）；
二、算法
（1）共轭梯度法算法
令 x0 0, r0 b Ax0 , p0 r0
对 k 0,1, 2,L 直到收敛完成
k (rk , rk ) / (pk , Apk ),
xk1 xk k pk ,
rk1 rk k Apk ,
k 1
(rk , (rk 1 ,
rk ) rk 1)
,
结束。
一、知识回顾
（2）预处理共轭梯度法原理 将Ax=b变形为C-1Ax=C-1b。
令 A~ C1AC T , ~x CT x,b C1b 则
~~ ~
Ax b A x b
~
要求
当A
其中
是c1o是对nA称d的2正( A最定)大矩特c阵o征n时d值2，(，Aco)nnd是2 (AA的) 最小1n 特K征(值A)，K(A)是A的条件数
3、通过A的各种近似分解得到预处理矩阵（如不完全分解）； 4、通过矩阵A的多项式构造预处理矩阵； 5、子结构、区域分裂、EBE预处理途径。
一、知识回顾
对角预优矩阵法 如果对称正定线性方程组( 1) 的系数矩阵 A 的对角元素相差 较大, 则可取预优矩阵C-1为
C 1 D1/ 2 , D diag (a11, a22,, ann ), A (aij )nn
A=[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;
C=chuli(A,0.009); A1=inv(C)*A*inv(C'); b1=inv(C)*b;
%SSOR此处程序改为
x=x+alpha*inv(C')*p
对角预优矩阵
x=
7.85971307544586 0.422926408295008 -0.0735922390240463 -0.540643016894626 0.0106261628540363
m= 4
tol= 4.731753657562764e-04
7.85971307544587 0.422926408294999 -0.0735922390236594 -0.540643016894632 0.0106261628540375
n= 5
tol =
6.745229690695623e-07
（2）预处理共轭梯度法代码
以对角预优矩阵为例
clc;clear
ep=10^(-2); [x,n,tol]=cg1(A1,b1,ep,C);
function C=chuli1(A) D=diag(diag(A)); C=D^(-1/2);
function C=chuli(A,w) D=diag(diag(A));L=tril(A,-1); C=(D-w*L)*D^(-1/2);
共轭梯度法的预处理
设 A 为 n 阶对称正定矩阵, 考虑线性方程组
Ax b
(1)
其中x 是未知向量,b是右端已知向量。
一、知识回顾
1、共轭梯度法 如果我们对于固定的k ，在 Rn 空间中取 k 个线性无关的向
量
p1, p2 ,L , pk
使 xk 满足
(xk ) min(x)
Span( p1 , p2 ,L , pk )
当 kn 发点。
时，
xn
一定是方程 (1)的精确解，这是共轭梯度法的出
一、知识回顾
2、Biblioteka Baidu处理共轭梯度法
（1）预处理共轭梯度法的引进
共轭梯度法理论上是一种直接方法，即按照这个算法得 到的 xn 应当是方程 Ax b 的精确解。但由于其数值的稳定 性不好，实际上由于舍入误差的影响 xn 不可能满足原方程。 如果把它看成一种迭代法，当A条件数较大时，它收敛的很 慢。因此引进了改进共轭梯度法的手段——预处理共轭梯度 法。
function [x,n,tol]=cg1(A,b,ep,C) r=b;p=r;n=0;x=zeros(length(b),1); while n<10
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)/(r1'*r1); p=r+beta*p; n=n+1; end
zk 1 M 1rk 1,
k 1 (zk 1,zk 1) / (rk , zk ),
结束。
pk 1 zk 1 k 1 pk
三、代码
（1）共轭梯度法代码
• function [x,n,tol]=cg(A,b,ep)
• r=b;p=r;n=0;x=zeros(length(b),1);
• while n<10
• 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;
• end
x=
pk1 rk k 1 pk ,
二、算法
（2）预处理共轭梯度算法
令
x0 0,
r0 b Ax0 , z0 M 1r0 , p0 z0 ,
对 k 0,1,L 直到收敛完成。
k (rk , rk ) / (pk , Apk ),
xk 1 xk k pk ,
rk1 rk k Apk ,
SSOR分裂
将A分裂为 A D CL CLT，其中D diag(a11, a22,ann) , CL为严格
下三角矩阵，则预处理矩阵C为：
C
1
(D CL)D1/ 2
(2 )
其中(0 1) 是参数，而且有
A CCT
1
(D CL)D1(D CL)T
(2 )
一、知识回顾
（3）预处理共轭梯度方法
由于预处理矩阵的构造方式很多，所以就形成了多种预处理 方法，总的来说，大致可归为如下几种途径：
1、取预处理矩阵（也称预优矩阵）为A的一个小带宽部分 （如三对角或对角线部分）；
2、通过A的分裂，尤其是线性稳定迭代法中的矩阵 A的分裂 构造预处理矩阵（如矩阵分裂法）；
二、算法
（1）共轭梯度法算法
令 x0 0, r0 b Ax0 , p0 r0
对 k 0,1, 2,L 直到收敛完成
k (rk , rk ) / (pk , Apk ),
xk1 xk k pk ,
rk1 rk k Apk ,
k 1
(rk , (rk 1 ,
rk ) rk 1)
,
结束。
一、知识回顾
（2）预处理共轭梯度法原理 将Ax=b变形为C-1Ax=C-1b。
令 A~ C1AC T , ~x CT x,b C1b 则
~~ ~
Ax b A x b
~
要求
当A
其中
是c1o是对nA称d的2正( A最定)大矩特c阵o征n时d值2，(，Aco)nnd是2 (AA的) 最小1n 特K征(值A)，K(A)是A的条件数
3、通过A的各种近似分解得到预处理矩阵（如不完全分解）； 4、通过矩阵A的多项式构造预处理矩阵； 5、子结构、区域分裂、EBE预处理途径。
一、知识回顾
对角预优矩阵法 如果对称正定线性方程组( 1) 的系数矩阵 A 的对角元素相差 较大, 则可取预优矩阵C-1为
C 1 D1/ 2 , D diag (a11, a22,, ann ), A (aij )nn
A=[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;
C=chuli(A,0.009); A1=inv(C)*A*inv(C'); b1=inv(C)*b;
%SSOR此处程序改为
x=x+alpha*inv(C')*p
对角预优矩阵
x=
7.85971307544586 0.422926408295008 -0.0735922390240463 -0.540643016894626 0.0106261628540363
m= 4
tol= 4.731753657562764e-04
7.85971307544587 0.422926408294999 -0.0735922390236594 -0.540643016894632 0.0106261628540375
n= 5
tol =
6.745229690695623e-07
（2）预处理共轭梯度法代码
以对角预优矩阵为例
clc;clear
ep=10^(-2); [x,n,tol]=cg1(A1,b1,ep,C);
function C=chuli1(A) D=diag(diag(A)); C=D^(-1/2);
function C=chuli(A,w) D=diag(diag(A));L=tril(A,-1); C=(D-w*L)*D^(-1/2);
共轭梯度法的预处理
设 A 为 n 阶对称正定矩阵, 考虑线性方程组
Ax b
(1)
其中x 是未知向量,b是右端已知向量。
一、知识回顾
1、共轭梯度法 如果我们对于固定的k ，在 Rn 空间中取 k 个线性无关的向
量
p1, p2 ,L , pk
使 xk 满足
(xk ) min(x)
Span( p1 , p2 ,L , pk )
当 kn 发点。
时，
xn
一定是方程 (1)的精确解，这是共轭梯度法的出
一、知识回顾
2、Biblioteka Baidu处理共轭梯度法
（1）预处理共轭梯度法的引进
共轭梯度法理论上是一种直接方法，即按照这个算法得 到的 xn 应当是方程 Ax b 的精确解。但由于其数值的稳定 性不好，实际上由于舍入误差的影响 xn 不可能满足原方程。 如果把它看成一种迭代法，当A条件数较大时，它收敛的很 慢。因此引进了改进共轭梯度法的手段——预处理共轭梯度 法。
function [x,n,tol]=cg1(A,b,ep,C) r=b;p=r;n=0;x=zeros(length(b),1); while n<10
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)/(r1'*r1); p=r+beta*p; n=n+1; end
zk 1 M 1rk 1,
k 1 (zk 1,zk 1) / (rk , zk ),
结束。
pk 1 zk 1 k 1 pk
三、代码
（1）共轭梯度法代码
• function [x,n,tol]=cg(A,b,ep)
• r=b;p=r;n=0;x=zeros(length(b),1);
• while n<10
• 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;
• end
x=
pk1 rk k 1 pk ,
二、算法
（2）预处理共轭梯度算法
令
x0 0,
r0 b Ax0 , z0 M 1r0 , p0 z0 ,
对 k 0,1,L 直到收敛完成。
k (rk , rk ) / (pk , Apk ),
xk 1 xk k pk ,
rk1 rk k Apk ,