共轭梯度法实验报告

1、共轭梯度法实验报告 数值代数实验报告 一、 实验名称： 用共轭梯度法解线性方程组。 二、实验目的： 进一步熟悉理解掌握共轭梯度法解法思路，提高 matlab 编程能力 。 三、 实验要求： 已知线性方程矩阵，应用共轭梯度法在相关软件编程求解线性方程 组的解。 四、实验原理： 1共轭梯度法： 考虑线性方程组 ax b 的求解问题，其中 a 是给定的 n 阶对称正定矩阵， b 是给定的 n 维向量， x 是待求解 的 n 维向量. 为此, 定义二次泛函 t t (x) x ax 2b x . 定理 1 设 a 对称正定 , 求方程组 ax b 的解，等价于求二次泛函 (x) 的极小值点 . 定理

2、1 表明，求解线性方程组问题就转化为求二次泛函 (x) 的极小值点问题 . 求解二次函数极小值问题，通常好像盲人下山那样，先给定一个初始向量 x 0 ，确定 一个下山方向 p 0 ，沿着经过点 x 0 而方向为 p 0 的直线 x x 0 p 0 找一个点 x x p ， 1 0 0 0 使得对所有实数 有 x p x p ， 0 0 0 0 0 即在这条直线上 x 1 使 (x) 达到极小. 然后从 x 1 出发，再确定一个下山的方向 p 1 ，沿着 直 线 x x p 再 跨 出 一 步 ， 即 找 到 1 使 得 x 在 x 2 x 1 1 p 1 达 到 极 小 ： 1 1 x p x

3、 p . 1 1 1 1 1 重复此步骤，得到一串 0 , 1 , 2 ,l 和 p 0 , p 1 , p 2 ,l ， 称 p k 为 搜索方向 ， k 为步长. 一般情况下，先在 x k 点找下山方向 p k ，再在直线 x x p 上确定步长 k 使 k k x p x p k k k k k , 最后求出 x x p . 然而对不同的搜索方向和步长，得到各种不同的算法 . k 1 k k k 1 由此，先考虑如何确定 k . 设从 x k 出发，已经选定下山方向 p k . 令 f x p k k t t x p a x p 2b x p k k k k k k 2 t 2 t p

4、ap r p x ， k k k k k 其中 r k b ap k . 由一元函数极值存在的必要条件有 t t f 2 p ap 2r p 0 k k k k 所确定的 即为所求步长 ，即 k 步长确定后，即可算出 t r p k k k t p ap k k . x x p . k 1 k k k t 此时，只要 r p 0 ，就有 k k 即 x x . k 1 k x x x p x 2 k 1 k k k k k t r p k k 2 t 2 t 0 p ap r p k k k k k k t p ap k k 再考虑如何确定下山方向 p k . 易知负梯度方向是 (x) 减小最

5、快的方向， 但简单分 析就会发现负梯度方向只是局部最佳的下山方向， 而从整体来看并非最佳 . 故采用新 的方法寻求更好的下山方向 共轭梯度法 . 下面给出共轭梯度法的具体计算过程 : 给定初始向量 x 0 ，第一步仍选用负梯度方向为下山方向，即 p 0 r 0 ，于是有 t r r 0 0 ,x x p ,r b ax . 0 t 1 0 0 0 1 0 p ap 0 0 对以后各步， 例如第 k+1 步（k 1）, 下山方向不再取 r k ，而是在过点由向量 r k 和 p k 1 所张成的二维平面 2 x | x x k r k p k 1 , , r 内找出使函数 下降最快的方向作为新的

6、下山方向 p k . 考虑 在 2 上的限制： , ( x k r k p k ) 1 t (x r p ) a(x r p ) k k k 1 k k k 1 t 2b (x r p ) . k k k 1 计算 关于 , 的偏导得： t t t 2 r ar r ap r r , k k k k 1 k k t t 2 r ap p ap , t 其中最后一式用到了 r p 1 0 ，这可由 r k 的定义直接验证 . 令 k k 1 k 1 k 1 k k 0 ， 即知 在 2 内有唯一的极小值点 x% x r p ， k 0 k 0 k 1 2 其中 0 和 0 满足 t t t r

7、ar r ap r r 0 k k 0 k k 1 k k , t t r ap p ap 0 k k 1 0 k 1 k 1 0. 由于 r k 0 必有 0 0 ，所以可取 1 0 p x% x r p k k k k 1 0 0 作为新的下山方向 . 显然，这是在平面 2 内可得的最佳下山方向 . 令 0 k ，则可 1 得 t r ap k k 1 . k 1 t p ap k 1 k 1 t 注：这样确定的 p k 满足 p ap 1 0 ，即 p k 与 p k 1 是相互共轭的 . k k 0 总结上面的讨论，可得如下的计算公式： t r p k k ， x k 1 x k k

8、p k ， k t p ap k k r b ax ， k 1 k 1 t r ap k 1 k ， p k 1 r k 1 k p k . k t p ap k k 在实际计算中，常将上述公式进一步简化，从而得到一个形式上更为简单而且对称 的计算公式 . 首先来简化 r k 1 的计算公式： r 1 b ax 1 b a(x p ) r ap . k k k k k k k k 因为 ap k 在计算 k 是已经求出，所以计算 r k 1 时可以不必将 x k 1 代入方程计算，而是 从递推关系 r k 1 b k ap k 得到. 再来简化 k 和 k 的计算公式 . 此处需要用到关系式

9、t t t r r 1 r p 1 r 1 p 0, k 1,2,k . k k k k k k 从而可导出 1 t t r r r , ， k 1 k 1 k 1 1 1 k t t t p ap p r r p r k k k k k 1 k k 1 t 1 t k k r r p r r . k k k 1 k 1 k k k k 由此可得 t t r r r r k k 1 1 . k k , . ， k t k t p ap r r k k k k 从而有求解对称正定方程组的共轭梯度法算法如下： x 初值 0 r b ax ； k 0 0 0 while r k 0 k k 1 if

10、 k 1 p r 0 0 3 else t t k r k r k r k r k 2 1 1 2 2 p r p k 1 k 1 k 2 k 2 end t t k r k r k p k ap k 1 1 1 1 1 x x p k k 1 k 1 k 1 r r ap k k 1 k 1 k 1 end x x k 注：该算法每迭代一次仅需要使用系数矩阵 a 做一次矩阵向量积运算 . 定理 2 由共轭梯度法得到的向量组 r 和 p i 具有如下基本性质： i t （1） p r 0， 0 i j k; i j t （2） r r 0 ， i j ， 0 i, j k; i j t （3）

11、 p ap 0， i j ， 0 i, j k; i j （4） span r ,k , r k span p ,k , p k (a, r ,k 1) ， 0 0 0 其中 k ( a,r , k 1) span r , ar ,k , a r ， 0 0 0 0 通常称之为 krylov 子空间. 下面给出共轭梯度法全局最优性定理： 定理 3 用共轭梯度法计算得到的近似解 x k 满足 x min x : x x (a,r ,k) k 0 0 或 x x x x x x a r k ， * min * : 0 ( , 0 , ) k a a 其中 t x x ax ， x * 是方程组 a

12、x b 的解， ( a, r 0 ,k) 是由所定义的 krylov 子空间. a 定理 2 表明，向量组 r 0 ,k ,r k 和 p 0 ,k , p k 分别是 krylov 子空间 ( a, r 0 ,k 1) 的正 交基和共轭正交基 . 由此可知，共轭梯度法最多 n 步便可得到方程组的解 x * . 因此， 理论上来讲，共轭梯度法是直接法 . 然而实际使用时，由于误差的出现，使 r k 之间的正交性很快损失，以致于其有 限步终止性已不再成立 . 此外，在实际应用共轭梯度法时，由于一般 n 很大，以至于 迭代 o n 次所耗费的计算时间就已经使用户无法接受了 . 因此，实际上将共轭梯

13、度 法作为一种迭代法使用， 而且通常是 r 是否已经很小及迭代次数是否已经达到最大 k 允许的迭代次数 k max 来终止迭代 . 从而得到解对称正定线性方程组的实用共轭梯度 法，其算法如下： x 初值 4 k 0; r b ax; t r r while b and k k max 2 k k 1 if k 1 p r else %; p r p end t ; ; ap p x x p t r r ; % ; r r end 算法中，系数矩阵 a 的作用仅仅是用来由已知向量 p 产生向量 ap ，这不仅可以 充分利用 a 的稀疏性，而且对某些提供矩阵 a 较为困难而由已知向量 p 产生向量

14、 ap 又十分方便的应用问题是十分有益的。 2算例 10 1 2 3 4 x 1 12 运用共轭梯度法求解下述方程，并解释你所观察到的结果 1 9 -1 2 - 3 x 27 2 2 -1 7 3 - 5 x 3 14 3 2 3 12 -1 x 17 解：已知共轭梯度法的 matlab 程序代码如下所示： 4 4 - 3 - 5 -1 15 x 12 functionx,n=conjgrad(a,b,x0) 5 %采用共轭梯度法求线性方程组 ax=b 的解 %线性方程组的系数矩阵： a %线性方程组中的常数向量： b %迭代初始向量： x0 %线性方程组的解： x %求出所需精度的解实际的迭

15、代步数： 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); 5 x=x0+alpha*p; r2=r1-alpha*a*p; if(r21.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 -

16、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.34

17、79 -13.9388 1.1303 11.5188 -26.8966 -2.2388 0.0829 13.2786 输出的迭代次数为 n = 10 共轭梯度法理论上只要迭代 n 步，就能得出方程组的解，但是由于存在计算误差， 即分数向小数转化时存在舍入误差，很难保证在第 n 步时得到准确解。 6 3 总结 本文首先给出最小二乘问题的定义，随后从盲人下山法开始讨论了共轭梯度法 的具体推导过程及其相关性质与算法 . 继而重点给出正则化方法的实用共轭梯度算 法并举例进行检验 . 最后，需要说明虽然正则化方法是求一般线性方程组 ax b， m n a r m n 且 a 列满秩的最小二乘解的一种方法且简单易行，但是也有许多不 足之处，如 m n 时一般无解； t a a 形成时运算量大， a 中某些信息会丢失；当 a病 态时其收敛性速度由于 t 2 2 (a a) 2 (a) 很大变得非常之慢等，故为了避免正则化方 法的缺点，还可运用残量极小化方法或残量正交化方法等更好的方法来解决此类问 题. 在实际的科学与工程问题中，常常将问题归结为