共轭梯度法实验报告

数值代数实验报告

一、实验名称：用共轭梯度法解线性方程组。

二、实验目的：进一步熟悉理解掌握共轭梯度法解法思路，提高matlab编程能力。 三、实验要求：已知线性方程矩阵，应用共轭梯度法在相关软件编程求解线性方程组的解。 四、实验原理：

1．共轭梯度法：

考虑线性方程组

Axb

的求解问题，其中A是给定的n阶对称正定矩阵，b是给定的n维向量，x是待求解的n维向量.为此,定义二次泛函

(x)xTAx2bTx.

定理1 设A对称正定,求方程组Axb的解，等价于求二次泛函(x)的极小值点. 定理1表明，求解线性方程组问题就转化为求二次泛函(x)的极小值点问题.求解二次函数极小值问题，通常好像盲人下山那样，先给定一个初始向量x0，确定一个下山方向p0，沿着经过点x0而方向为p0的直线xx0p0找一个点

x1x00p0，

使得对所有实数有

x00p0x0p0，

即在这条直线上x1使(x)达到极小.然后从x1出发，再确定一个下山的方向p1，沿着直

线xx1p1再跨出一步，即找到1使得x在x2x11p1达到极小：

x11p1x1p1.

重复此步骤，得到一串

0,1,2,xxkpk上确定步长k使

和 p0,p1,p2,，

称pk为搜索方向，k为步长.一般情况下，先在xk点找下山方向pk，再在直线

xkkpkxkpk,

最后求出xk1xkkpk.然而对不同的搜索方向和步长，得到各种不同的算法.

-可编辑修改-

。

由此，先考虑如何确定k.设从xk出发，已经选定下山方向pk.令 fxkpk

xkpkAxkpk2bTxkpk

T 2pkApk2rkTpkxk，

T其中rkbApk.由一元函数极值存在的必要条件有

Tf2pkApk2rkTpk0

所确定的即为所求步长k，即 步长确定后，即可算出

此时，只要rkTpk0，就有

rkTpk. kTpkApkxk1xkkpk.

xk1xkxkkpkxk

TApk2krkTpk k2pkrkTpk2即xk1xk.

TpkApk0

再考虑如何确定下山方向pk.易知负梯度方向是(x)减小最快的方向，但简单分析就会发现负梯度方向只是局部最佳的下山方向，而从整体来看并非最佳.故采用新的方法寻求更好的下山方向——共轭梯度法. 下面给出共轭梯度法的具体计算过程:

给定初始向量x0，第一步仍选用负梯度方向为下山方向，即p0r0，于是有

r0Tr00T,x1x00p0,r1bAx0.

p0Ap0对以后各步，例如第k+1步（k1）,下山方向不再取rk，而是在过点由向量rk和pk1所张成的二维平面

2{x|xxkrkpk1,,R}

内找出使函数下降最快的方向作为新的下山方向pk.考虑在2上的：

,(xkrkpk1)

(xkrkpk1)TA(xkrkpk1)

2bT(xkrkpk1). 计算关于,的偏导得： 2rTArrTAprTr,kkkk1kk

TTT2rAppkk1k1Apk1,0r其中最后一式用到了rkpk1，这可由的定义直接验证.令 k0， 即知在2内有唯一的极小值点

xxk0rk0pk1，

-可编辑修改-

。

其中0和0满足

0rkTArk0rkTApk1rkTrk, TT0rkApk10pk1Apk10.1由于rk0必有00，所以可取

pk0p

00k1作为新的下山方向.显然，这是在平面2内可得的最佳下山方向.令k10，则可

0得

rkTApk1k1T.

pApk1k1TApk10，即pk与pk1是相互共轭的. 注：这样确定的pk满足pkxxkrk总结上面的讨论，可得如下的计算公式：

rkTpk ， xk1xkkpk， kTpkApk rk1bAxk1，

rkT1Apk， pk1rk1kpk. kTpkApk在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称的计算公式.首先来简化rk1的计算公式：

rk1bAxk1bA(xkkpk)rkkApk.

因为Apk在计算k是已经求出，所以计算rk1时可以不必将xk1代入方程计算，而是从递推关系rk1bkApk得到.

再来简化k和k的计算公式.此处需要用到关系式

rkTrk1rkTpk1rkT1pk0, k1,2,.

从而可导出

1rkT1rkT1rk1,， 1Tk1TTpkApkpkrkrk1pkrk

k1Tk1Trrpr kkk1k1krk.

由此可得

kkrkTrkrkT1rk1kT,， kT..

pkApkrkrk从而有求解对称正定方程组的共轭梯度法算法如下：

x0初值

r0bAx0；k0

while rk0

kk1 if k1 p0r0

-可编辑修改-

。

else

k2rkT1rk1rkT2rk2 pk1rk1k2pk2 end

Tk1rkT1rk1pk1Apk1

xkxk1k1pk1

rkrk1k1Apk1

end

xxk

注：该算法每迭代一次仅需要使用系数矩阵A做一次矩阵向量积运算. 定理2 由共轭梯度法得到的向量组ri和pi具有如下基本性质： （1）piTrj0， 0ijk; （2）riTrj0， ij，0i,jk; （3）piTApj0， ij，0i,jk; （4）span{r0,其中

,rk}span{p0,,pk}(A,r0,k1)，

,Akr0}，

(A,r0,k1)span{r0,Ar0,定理3 用共轭梯度法计算得到的近似解xk满足

通常称之为Krylov子空间.下面给出共轭梯度法全局最优性定理：

xkminx:xx0(A,r0,k)

或

xkx*其中xAAminxx*,rk和p0,A:xx0(A,r0,k)，

(A,r0,k)是由所定义的Krylov子空间. xTAx，x*是方程组Axb的解，

,pk分别是Krylov子空间(A,r0,k1)的正

定理2表明，向量组r0,交基和共轭正交基.由此可知，共轭梯度法最多n步便可得到方程组的解x*.因此，理论上来讲，共轭梯度法是直接法.

然而实际使用时，由于误差的出现，使rk之间的正交性很快损失，以致于其有限步终止性已不再成立.此外，在实际应用共轭梯度法时，由于一般n很大，以至于迭代On次所耗费的计算时间就已经使用户无法接受了.因此，实际上将共轭梯度法作为一种迭代法使用，而且通常是rk是否已经很小及迭代次数是否已经达到最大允许的迭代次数kmax来终止迭代.从而得到解对称正定线性方程组的实用共轭梯度法，其算法如下：

x初值

-可编辑修改-

。

k0;rbAx;rTr

while

b2andkkmax

kk1

if k1 pr else

;prp end

Ap;pT;xxp rr;;rTr end

算法中，系数矩阵A的作用仅仅是用来由已知向量p产生向量Ap，这不仅可以充分利用A的稀疏性，而且对某些提供矩阵A较为困难而由已知向量p产生向量

Ap又十分方便的应用问题是十分有益的。

2．算例 1210 1 2 3 4x1 1 9 -1 2 -3x27运用共轭梯度法求解下述方程，并解释你所观察到的结果 22 -1 7 3 -5x3 14 3 2 3 12 -117x解：已知共轭梯度法的MATLAB4程序代码如下所示：

 124 -3 -5 -1 15x5function[x,n]=conjgrad(A,b,x0) %采用共轭梯度法求线性方程组Ax=b的解 %线性方程组的系数矩阵：A %线性方程组中的常数向量：b %迭代初始向量：x0 %线性方程组的解：x

%求出所需精度的解实际的迭代步数：n r1=b-A*x0; p=r1; n=0;

for i=1:rank(A)

if(dot(p,A*p)<1.0e-50) break; end

alpha=dot(r1,r1)/dot(p,A*p);

-可编辑修改-

。

x=x0+alpha*p; r2=r1-alpha*A*p; if(r2<1.0e-50) break; end

belta=dot(r2,r2)/dot(r1,r1); p=r2+belta*p; n=n+1; end

用共轭梯度法求解，在MATLAB命令窗口中输入求解程序：

>> A=[10 1 2 3 4 2 4 2 1 1;1 9 -1 2 -3 3 -3 -1 9 2;2 -1 7 3 -5 4 -5 7 -1 3;3 2 3 12 -1 5 -1 3 2 4;4 -3 -5 -1 15 1 3 2 3 2;2 3 4 5 1 10 4 5 3 4;4 -3 -5 -1 3 4 9 1 -3 2;2 -1 7 3 2 5 1 7 -1 1;1 9 -1 2 -3 3 -3 -1 12 10;1 2 3 4 2 4 2 1 10 15];

>> b=[12;-27;14;-17;12;12;-27;14;-17;12]; >> x0=[0;0;0;0;0;0;0;0;0;0]; >> [x,n]=conjgrad(A,b,x0) 输出的计算结果为： x = 8.5669 -7.1981 -7.3479 -13.9388 1.1303 11.5188 -26.66 -2.2388 0.0829 13.2786

输出的迭代次数为 n = 10

共轭梯度法理论上只要迭代n步，就能得出方程组的解，但是由于存在计算误差，即分数向小数转化时存在舍入误差，很难保证在第n步时得到准确解。

-可编辑修改-

。

3． 总结

本文首先给出最小二乘问题的定义，随后从盲人下山法开始讨论了共轭梯度法的具体推导过程及其相关性质与算法.继而重点给出正则化方法的实用共轭梯度算法并举例进行检验.最后，需要说明虽然正则化方法是求一般线性方程组Axb，

ARmnmn且A列满秩的最小二乘解的一种方法且简单易行，但是也有许多不足之处，如mn时一般无解；ATA形成时运算量大，A中某些信息会丢失；当A病

2(A)很大变得非常之慢等，故为了避免正则化方态时其收敛性速度由于2(ATA)2法的缺点，还可运用残量极小化方法或残量正交化方法等更好的方法来解决此类问题.

在实际的科学与工程问题中，常常将问题归结为一个线性方程组的求解问题，而求解线性方程组的数值解法大体上可分为直接法和迭代法两大类.直接法是指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法.因此，直接法又称为精确法.迭代法则是采取逐次逼近的方法，亦即从一个初始向量出发，按照一定的计算格式，构造一个向量的无穷序列,其极限才是方程组的精确解，只经过有限次运算得不到精确解.当线性方程组的系数矩阵为对称正定矩阵时，我们常用共轭梯度法（或简称CG法）求解，目前有关的方法与理论已经相当成熟，并且已成为求解大型稀疏线性方程组最受欢迎的一类方法.

-可编辑修改-

。

欢迎您的下载， 资料仅供参考！

致力为企业和个人提供合同协议，策划案计划书，学习课件等等

打造全网一站式需求

-可编辑修改-