共轭梯度法实验报告.doc

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

的求解问题，其中A是给定的n阶对称正定矩阵，b是给定的n维向量，是待求解的n维向量.为此,定义二次泛函
.
定理1 设A对称正定,求方程组的解，等价于求二次泛函的极小值点.
定理1表明，求解线性方程组问题就转化为求二次泛函的极小值点问题.求解二次函数极小值问题，通常好像盲人下山那样，先给定一个初始向量，确定一个下山方向，沿着经过点而方向为的直线找一个点
，
使得对所有实数有
，
即在这条直线上使达到极小.然后从出发，再确定一个下山的方向，沿着直
线再跨出一步，即找到使得在达到极小： .
重复此步骤，得到一串
和 ，
称为搜索方向，为步长.一般情况下，先在点找下山方向，再在直线上确定步长使

最后求出.然而对不同的搜索方向和步长，得到各种不同的算法.
由此，先考虑如何确定.设从出发，已经选定下山方向.令

，
其中.由一元函数极值存在的必要条件有

所确定的即为所求步长，即
.
步长确定后，即可算出
.
此时，只要，就有

即.
再考虑如何确定下山方向.易知负梯度方向是减小最快的方向，但简单分析就会发现负梯度方向只是局部最佳的下山方向，而从整体来看并非最佳.故采用新的方法寻求更好的下山方向——共轭梯度法.
下面给出共轭梯度法的具体计算过程:
给定初始向量，第一步仍选用负梯度方向为下山方向，即，于是有
.
对以后各步，例如第k+1步（k1）,下山方向不再取，而是在过点由向量和所张成的二维平面

内找出使函数下降最快的方向作为新的下山方向.考虑在上的限制：

.
计算关于的偏导得：

其中最后一式用到了，这可由的定义直接验证.令
，
即知在内有唯一的极小值点
，
其中和满足

由于必有，所以可取

作为新的下山方向.显然，这是在平面内可得的最佳下山方向.令，则可得

注：这样确定的满足，即与是相互共轭的.
总结上面的讨论，可得如下的计算公式：
， ，
，
， .
在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称的计算公式.首先来简化的计算公式：
.
因为在计算是已经求出，所以计算时可以不必将代入方程计算，而是从递推关系得到.
再来简化和的计算公式.此处需要用到关系式
.
从而可导出
，

.
由此可得
， .
从而有求解对称正定方程组的共轭梯度法算法如下：
初值
；
while

if

else

end

end

注：该算法每迭代一次仅需要使用系数矩阵做一次矩阵向量积运算.
定理2 由共轭梯度法得到的向量组和具有如下基本性质：
（1），
（2）， ，
（3）， ，
（4），
其中
，
通常称之为Krylov子空间.下面给出共轭梯度法全局最优性定理：
定理3 用共轭梯度法计算得到的近似解满足

或
，
其中，是方程组的解，是由所定义的Krylov子空间.
定理2表明，向量组和分别是Krylov子空间的正交基和共轭正交基.由此可知，共轭梯度法最多n步便可得到方程组的解.因此，理论上来讲，共轭梯度法是直接法.
然而实际使用时，由于误差的出现，使之间的正交性很快损失，以致于其有限步终止性已不再成立.此外，在实际应用共轭梯度法时，由于一般很大，以至于迭代次所耗费的计算时间就已经使用户无法接受了.因此，实际上将共轭梯度法作为一种迭代法使用，而且通常是是否已经很小及迭代次数是否已经达到最大允许的迭代次数来终止迭代.从而得到解对称正定线性方程组的实用共轭梯度法，其算法如下：
初值

while

if

else

end

end
算法中，系数矩阵的作用仅仅是用来由已知向量产生向量，这不仅可以充分利用的稀疏性，而且对某些提供矩阵较为困难而由已知向量产生向量又十分方便的应用问题是十分有益的。
2． 算例
运用共轭梯度法求解下述方程，并解释你所观察到的结果

解：已知共轭梯度法的MATLAB程序代码如下所示：
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.8966
-2.2388
0.0829
13.2786

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

3． 总结
本文首先给出最小二乘问题的定义，随后从盲人下山法开始讨论了共轭梯度法的具体推导过程及其相关性质与算法.继而重点给出正则化方法的实用共轭梯度算法并举例进行检验.最后，需要说明虽然正则化方法是求一般线性方程组，且列满秩的最小二乘解的一种方法且简单易行，但是也有许多不足之处，如时一般无解；形成时运算量大，中某些信息会丢失；当病态时其收敛性速度由于很大变得非常之慢等，故为了避免正则化方法的缺点，还可运用残量极小化方法或残量正交化方法等更好的方法来解决此类问题.
在实际的科学与工程问题中，常常将问题归结为一个线性方程组的求解问题，而求解线性方程组的数值解法大体上可分为直接法和迭代法两大类.直接法是指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法.因此，直接法又称为精确法.迭代法则是采取逐次逼近的方法，亦即从一个初始向量出发，按照一定的计算格式，构造一个向量的无穷序列,其极限才是方程组的精确解，只经过有限次运算得不到精确解.当线性方程组的系数矩阵为对称正定矩阵时，我们常用共轭梯度法（或简称CG法）求解，目前有关的方法与理论已经相当成熟，并且已成为求解大型稀疏线性方程组最受欢迎的一类方法.