共轭梯度法实验报告

1 数值代数实验报告数值代数实验报告 1 实验名称 实验名称 用共轭梯度法解线性方程组 2 2 实验目的 实验目的 进一步熟悉理解掌握共轭梯度法解法思路 提高 matlab 编程能力 3 3 实验要求 实验要求 已知线性方程矩阵 应用共轭梯度法在相关软件编程求解线性方程 组的解 4 4 实验原理 实验原理 1 1 共轭梯度法 共轭梯度法 考虑线性方程组 Axb 的求解问题 其中 A 是给定的 n 阶对称正定矩阵 b 是给定的 n 维向量 是待求x 解的 n 维向量 为此 定义二次泛函 2 TT xx Axb x 定理定理 1 1 设 A 对称正定 求方程组的解 等价于求二次泛函的极小值点 Axb x 定理 1 表明 求解线性方程组问题就转化为求二次泛函的极小值点问题 x 求解二次函数极小值问题 通常好像盲人下山那样 先给定一个初始向量 确定 0 x 一个下山方向 沿着经过点而方向为的直线找一个点 0 p 0 x 0 p 00 xxp 1000 xxp 使得对所有实数有 00000 xpxp 即在这条直线上使达到极小 然后从出发 再确定一个下山的方向 沿 1 x x 1 x 1 p 着直 线再跨出一步 即找到使得在达到极小 11 xxp 1 x 2111 xxp 11111 xpxp 重复此步骤 得到一串 和 012 012 pp p 称为搜索方向搜索方向 为步长步长 一般情况下 先在点找下山方向 再在直线 k p k k x k p 上确定步长使 kk xxp k kkkkk xpxp 2 最后求出 然而对不同的搜索方向和步长 得到各种不同的算法 1kkkk xxp 由此 先考虑如何确定 设从出发 已经选定下山方向 令 k k x k p kk fxp 2 T T kkkkkk xpA xpbxp 2 2 TT kkkkk p Apr px 其中 由一元函数极值存在的必要条件有 kk rbAp 220 TT kkkk fp Apr p 所确定的即为所求步长 即 k T kk k T kk r p p Ap 步长确定后 即可算出 1kkkk xxp 此时 只要 就有0 T kk r p 1kkkkkk xxxpx 2 2 20 T kk TT kkkk kk T kk r p p Apr p p Ap 即 1kk xx 再考虑如何确定下山方向 易知负梯度方向是减小最快的方向 但简单 k p x 分析就会发现负梯度方向只是局部最佳的下山方向 而从整体来看并非最佳 故采 用新的方法寻求更好的下山方向 共轭梯度法共轭梯度法 下面给出共轭梯度法的具体计算过程 给定初始向量 第一步仍选用负梯度方向为下山方向 即 于是有 0 x 00 pr 00 0100010 00 T T r r xxp rbAx p Ap 对以后各步 例如第 k 1 步 k1 下山方向不再取 而是在过点由向量和 k r k r 所张成的二维平面 1k p 21 kkk x xxrpR 内找出使函数下降最快的方向作为新的下山方向 考虑在上的限制 k p 2 1 kkk xrp 11 T kkkkkk xrpA xrp 1 2 T kkk bxrp 计算关于的偏导得 1 111 2 2 TTT kkkkkk TT kkkk r Arr Apr r r AppAp 其中最后一式用到了 这可由的定义直接验证 令 1 0 T kk r p k r 0 即知在内有唯一的极小值点 2 3 001kkk xxrp 其中和满足 0 0 001 01011 0 TTT kkkkkk TT kkkk r Arr Apr r r AppAp 由于必有 所以可取0 k r 0 0 0 1 00 1 kkkk pxxrp 作为新的下山方向 显然 这是在平面内可得的最佳下山方向 令 则可 2 0 1 0 k 得 1 1 11 T kk k T kk r Ap pAp 注 这样确定的满足 即与是相互共轭的 k p 1 0 T kk p Ap k p 1k p 总结上面的讨论 可得如下的计算公式 T kk k T kk r p p Ap 1kkkk xxp 11kk rbAx 1 T kk k T kk rAp p Ap 11kkkk prp 在实际计算中 常将上述公式进一步简化 从而得到一个形式上更为简单而且对称 的计算公式 首先来简化的计算公式 1k r 11 kkkkkkkk rbAxbA xprAp 因为在计算是已经求出 所以计算时可以不必将代入方程计算 而是 k Ap k 1k r 1k x 从递推关系得到 1kkk rbAp 再来简化和的计算公式 此处需要用到关系式 k k 111 0 TTT kkkkkk r rr prp 1 2 k 从而可导出 111 1 TT kkk k rrr 1 11 TTT kkkkkkk kk p Apprrp r 11 11 TT kkkkkk kk rrpr r 由此可得 T kk k T kk r r p Ap 11 T kk k T kk rr r r 从而有求解对称正定方程组的共轭梯度法算法如下 初值 0 x 00 rbAx 0k whilewhile 0 k r 1kk ifif 1k 4 00 pr elseelse 21122 TT kkkkk rrrr 1122kkkk prp endend 11111 TT kkkkk rrpAp 111kkkk xxp 111kkkk rrAp endend k xx 注 该算法每迭代一次仅需要使用系数矩阵做一次矩阵向量积运算 A 定理定理 2 2 由共轭梯度法得到的向量组和具有如下基本性质 i r i p 1 0 T ij p r 0 ijk 2 0 T ij r r ij 0 i jk 3 0 T ij p Ap ij 0 i jk 4 000 1 kk span rrspan ppA r k 其中 0000 1 k A r kspan r ArA r 通常称之为 KrylovKrylov 子空间子空间 下面给出共轭梯度法全局最优性定理 定理定理 3 3 用共轭梯度法计算得到的近似解满足 k x 00 min k xxxxA r k 或 00 min k AA xxxxxxA r k 其中 是方程组的解 是由所定义的 KrylovKrylov 子空 T A xx Ax xAxb 0 A r k 间 定理 2 表明 向量组和分别是 Krylov 子空间的正 0 k rr 0 k pp 0 1 A r k 交基和共轭正交基 由此可知 共轭梯度法最多 n 步便可得到方程组的解 因此 x 理论上来讲 共轭梯度法是直接法 然而实际使用时 由于误差的出现 使之间的正交性很快损失 以致于其有 k r 限步终止性已不再成立 此外 在实际应用共轭梯度法时 由于一般很大 以至于n 迭代次所耗费的计算时间就已经使用户无法接受了 因此 实际上将共轭梯度 O n 法作为一种迭代法使用 而且通常是是否已经很小及迭代次数是否已经达到最 k r 大允许的迭代次数来终止迭代 从而得到解对称正定线性方程组的实用共轭梯度 max k 5 法 其算法如下 初值x 0 k rbAx T r r whilewhile max 2 band kk 1kk ifif 1k pr elseelse prp endend T Appxxp T rrr r endend 算法中 系数矩阵的作用仅仅是用来由已知向量产生向量 这不仅可以ApAp 充分利用的稀疏性 而且对某些提供矩阵较为困难而由已知向量产生向量AAp 又十分方便的应用问题是十分有益的 Ap 2 2 算例算例 运用共轭梯度法求解下述方程 并解释你所观察到的结果 12 17 14 27 12 15 1 5 3 4 1 12 3 2 3 5 3 7 1 2 3 2 1 9 1 4 3 2 1 10 5 4 3 2 1 x x x x x 解 已知共轭梯度法的 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 6 end alpha dot r1 r1 dot p A p x x0 alpha p r2 r1 alpha A p if r2 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 7 共轭梯度法理论上只要迭代 n 步 就能得出方程组的解 但是由于存在计算误差 即分数向小数转化时存在舍入误差 很难保证在第 n 步时得到准确解 3 3 总结总结 本文首先给出最小二乘问题的定义 随后从盲人下山法开始讨论了共轭梯度法 的具体推导过程及其相关性质与算法 继而重点给出正则化方法的实用共轭梯度算 法并举例进行检验 最后 需要说明虽然正则化方法是求一般线性方程组 Axb 且列满秩的最小二乘解的一种方法且简单易行 但是也有许多不 m n AR mn A 足之处 如时一般无解 形成时运算量大 中某些信息会丢失 当mn T A AA 病态时其收敛性速度由于很大变得非常之慢等 故为了避免正则A 2 22 T A AA 化方法的缺点 还可运用残量极小化方法或残量正交化方法等更好的方法来解决此 类问题 在实际的科学与工程问题中 常常将问题归结为一个线性方程组的求解问题 而求解线性方程组的数值解法大体上可分