用共轭梯度分解求解最小二乘问题
作者：蒲小丽
来源：《新校园·中旬刊》2011年第11期
        摘 要：本文先讨论了求解对称正定线性方程组的共轭梯度法.然后对系数矩阵列满秩的线性方程组运用正则化方法将其转化为对称正定线性方程组后再运用实用共轭梯度法进行求解，并举实例证明。
        关键词：共轭梯度法；正则化方法；最小二乘问题；线性方程组
       
        一、引言
        在实际的科学与工程问题中，常常将问题归结为一个线性方程组的求解问题，而求解线性方程组的数值解法大体上可分为直接法和迭代法两大类.直接法是指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法.因此，直接法又称为精确法.迭代法则是采取逐次逼近的方法，亦即从一个初始向量出发，按照一定的计算格式，构造一个向量的无穷序列,
其极限才是方程组的精确解，只经过有限次运算得不到精确解.当线性方程组的系数矩阵为对称正定矩阵时，我们常用共轭梯度法求解，目前有关的方法与理论已经相当成熟，并且已成为求解大型稀疏线性方程组最受欢迎的一类方法.
        二、最小二乘问题
        定义1[1] 给定矩阵A∈Rm×n，A列满秩及向量b∈Rm，确定x∈Rn使得||b-Ax||2=||r(x)||2=■||r(y)||■=■||Ay-b||■.
        该为问题称为最小二乘问题，简记为LS（Least Squares）问题，其中r(x)称为残向量.
        三、共轭梯度法
        共轭梯度法简称CG(conjugate gradient)方法,又称共轭斜量法,是一种求解大型稀疏对称正定方程组十分有效的方法.它是一种变分方法,对应于求一个二次函数的极值.
        考虑线性方程组Ax=b的求解问题，其中A是给定的n阶对称正定矩阵，b是给定的n维向量， x是待求解的n维向量.为此,定义二次泛函ψ(x)=xTAx-2bTx.
        定理1设A对称正定,求方程组Ax=b的解，等价于求二次泛函ψ(x)的极小点.
        下面给出共轭梯度法全局最优性定理：
        定理2用共轭梯度法计算得到的近似解xk满足
        ψ(xk)=min{ψ(x):x∈x0+K(A,r0,k)}
        或
        ||xk-x*||A=min{||x-x*||A:x∈x0+K(A,r0,k)}，
        其中||x||A=■，x* 是方程组Ax=b的解，K(A,r0,k)是由所定义的Krylov子空间.
        然而实际用时，由于误差出现，使rk之间的正交性很快损失，以致于其有限步终止性已不再成立.此外，在实际应用共轭梯度法时，由于一般n很大，以至于迭代O(n)次所耗费的计算时间就已经使用户无法接受了.因此，实际上将共轭梯度法作为一种迭代法使用，而且通常是||rk||是否已经很小及迭代次数是否已经达到最大允许的迭代次数kmax来终止迭代.从而得到解对称正定线性方程组的实用共轭梯度法，其算法如下：
        x=初值k =0；r=b-Ax;ρ=rTr while(■>ε||b||2)and(k
        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又十分方便的应用问题有益.
        四、共轭梯度法求解最小二乘问题的正则化方程组（法方程组）
        定理3当A列满秩时，求最小二乘问题的解等价于解ATAx=ATb.
        应用共轭梯度法于对称正定方程组ATAx=ATb来求方程组Ax=b，A∈Rm×n(m≥n)且A列满秩的最小二乘解，即为Krylov子空间法中的正则化方法.由A列满秩有ATA对称正定，则方程组ATAx=ATb，A∈Rm×n(m≥n)存在唯一解.
        下面给出其实用共轭梯度法的详细算法且算法中不出现计算ATA情形:
        x=初值k=0;b=ATb;m=Ax;r=b-ATm;
        while(■>ε||b||2)and(k
        p=r else β=ρ/■;p=r+βp end
        n=Ap;ω=ATn;α=ρ/pTω;x=x+αp r=r-αω;■=ρ;ρ=rTr end
        注：算法中采用了两次矩阵向量积来避免出现计算ATA情形.
        算例
        编写实用共轭梯度法的Matlab程序求解方程组ATAx=ATb，其中
        A=[1 -1;-1 1;2 -2;-3 1], b=[1 ２ ３ ４]Ｔ.
        解 先建立conjgrad.m文件，内容如下:
        function [x]=conjgrad(A,b,x)
        A=[1 -1;-1 1;2 -2;-3 1];
        b=[1;2;3;4];x=rand(2,1);k=0;b=A'*b;m=A*x;r=b-A'*m;T=r'*r;
        While norm(r)>1e-10*norm(b) k=k+1;
        if k==1 p=r;else B=T/t;p=r+B*p;end
        n=A*p;w=A'*n;a=T/(p'*w);x=x+a*p;r=r-a*w;t=T;T=r'*r;end
        运行后得
        ans =-2.4167 -3.2500
        即有方程组的数值解x=-2.4167-3.2500.
        可验证通过Matlab程序求得的数值解满足精度要求.
        五、总结
        本文首先给出最小二乘问题的定义，随后从盲人下山法开始讨论了共轭梯度法的具体推导过程及其相关性质与算法.继而重点给出正则化方法的实用共轭梯度算法进行检验.最后，
需要说明虽然正则化方法是求一般线性方程组Ax=b，A∈Rm×n(m≥n)正则化最小二乘问题且A列满秩的最小二乘解的一种方法且简单易行，但是也有许多不足之处，如m>n时一般无解；ATA形成时运算量大，A中某些信息会丢失；当A病态时其收敛性速度由于K2(ATA)=K22(A)很大变得非常之慢等，故为了避免正则化方法的缺点，可运用残量极小化方法或残量正交化方法等更好的方法来解决此类问题.
       
        【参考文献】
        [1]徐树方，高立，张平文.《数值线性代数》[M].北京：北京大学出版社.2000，（139-151）.
        [2]施光燕，董加礼.《最优化方法》[M].北京：高等教育出版社.1999，（47-52）.
        [3]杜延松，沈艳军，覃太贵.《数值分析及实验》[M].北京：科学出版社.2004.
        [4]胡守信，李伯年.《基于MATLAB的数学实验》[M].北京：科学出版社.2004，（22-27）.
        [5]谷同祥，刘兴平，迟学斌.《对称正定问题多维搜索方向共轭梯度法的收敛性理论》[J].计算数学.2002，26（1）:117-128.

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。