超定方程最小二乘解

矩阵的每一行代表一个方程，m行代表m个线性联立方程。n列代表n个变量。如果m是独立方程数，根据mn、m=n、mn确定方程是‘欠定’、‘适定’还是‘超定’。超定方程组：方程个数大于未知量个数的方程组。对于方程组Ra=y，R为n×m矩阵，如果R列满秩，且nm超定方程一般是不存在解的矛盾方程。例如，如果给定的三点不在一条直线上，我们将无法得到这样一条直线，使得这条直线同时经过给定这三个点。也就是说给定的条件（限制）过于严格，导致解不存在。在实验数据处理和曲线拟合问题中，求解超定方程组非常普遍。比较常用的方法是最小二乘法。形象的说，就是在无法完全满足给定的这些条件的情况下，求一个最接近的解。曲线拟合的最小二乘法要解决的问题，实际上就是求以上超定方程组的最小二乘解的问题。欠定方程组:方程个数小于未知量个数的方程组。对于方程组Ra=y，R为n×m矩阵，且nm。则方程组有无穷多组解，此时称方程组为欠定方程组。内点法和梯度投影法是目前解欠定方程组的常用方法。超定方程组最小二乘解课程设计最小二乘法广泛地应用于工程计算中，用最小二乘法消除（平滑）误差，用最小二乘法从有噪声的数据中提取信号，从海量数据中找出数据变化的趋势，……。甚至利用简单函数计算复杂函数的近似值，我们并不期望它的近似值多么精确（事实上很多时候也不用很精确），尽管如此还是希望计算出的近似数据与原始数据之间有相似之处。如果从线性代数角度来理解最小二乘法，实际上是将一个高维空间的向量投影到低维子空间所涉及的工作。一、超定方程组的最小二乘解当方程组GX=b的方程数多于未知数个数时，对应的系数矩阵G的行数大于列数，此时方程组被称为是超定方程组。设G=(giu)m×n，当mn时即所谓的高矩阵，绝大多数情况下，超定方程组没有古典意义下的解。超定方程组的最小二乘解是一种广义解，是指使残差r=b–GX的2-范数达取极小值的解，即22*||||min||||GXbGXbmRX该问题是一个优化问题。命题1：如果X*是正规方程组GTGX=GTb的解，则X*是超定方程组GX=b的最小二乘解证由题设可得，GT(b–GX*)=0。对任意n维向量Y，显然有(X*–Y)TGT(b–GX*)=0考虑残差2-范数平方，由22**22||)()(||||||YXGGXbGYb上式右端利用内积，得22*22*22*22||||||)(||||||||||GXbYXGGXbGYb从而有||b–GY||2≥||b–GX*||2等式仅当Y=X*时成立。所以X*是超定方程组GX=b的最小二乘解。命题2：如果X*是超定方程组GX=b的最小二乘解，则X*满足正规方程组GTGX=GTb证由题设，22*||||min||||GXbGXbmRX，利用2-范数与内积关系，知X*是下面二次函数的极小值点(X)=(GX，GX)–2(GX，b)+(b，b)取任意n维向量v，对任意实数t，构造一元函数g(t)=(X*+tv)显然，g(t)是关于变量t的二次函数g(t)=(G(X*+tv)，G(X*+tv))–2(G(X*+tv)，b)+(b，b)=g(0)+2t[(GX*，Gv)–(Gv，b)]+t2(Gv，Gv)由题设t=0是g(t)的极小值点。由极值必要条件，得0)0(g。即(GX*，Gv)–(Gv，b)=0将左端整理化简，便得(Gv，GX*–b)=0利用内积性质，得(v，GT(GX*–b))=0由v的任意性，得GT(GX*–b)=0二、最小二乘解的几何意义首先考虑一个简单的超定方程组6311215342yx该方程组的右端向量是三维向量，系数矩阵的每一列也是三维向量，但待求的未知向量却是二维向量。将系数矩阵按列分块，G=[1，2]，记右端向量为。则方程组求解问题可表示为求组合系数x和y使x1+y2=的向量的线性组合问题。由于两个向量1，2不构成三维空间的一组基，所以一般情况下这一问题无解。而由向量1，2张成的子空间span{1，2}是一张平面，记为。则超定方程组的最小二乘解实际上是求X*，使GX*恰好等于在平面上的投影。而最小二乘解所对应的残差向量则垂直于向量GX*。事实上，由正规方程组GTGX=GTb得GT(b–GX*)=0上式的几何意义可解释为：最小二乘解的残差向量与超定方程组的系数矩阵G的所有列向量正交。从而GX*rmin=b–GX*(X*)TGT(b–GX*)=0所以(GX*，b–GX*)=0三、最小二乘解的两种方法超定方程组的最小二乘解是指正规方程组GTGX=GTb的解。如果系数矩阵(GTG)可逆，则正规方程组有唯一解。此时，最小二乘解可以形式地写为如下形式X=(GTG)-1GTb两种常用的方法如下1．用对称矩阵的三角分解法解正规方程组GTGX=GTb；记A=GTG，则A是对称矩阵，由三角分解A=LDLT，其中L是下三角矩阵，D是对角矩阵。将这一算法写过三个过程：①解下三角方程组：LY1=GTb；②解对角方程组：DY2=Y1；③解上三角方程组：LTY3=Y22．用矩阵的QR分解直接求解超定方程组由QR分解（正交三角分解）G=QR，其中Q是正交矩阵，R是上三角矩阵。将QR分解代入最小二乘解表达式中，得X=(RTQTQR)-1(QR)Tb=R-1QTb由此可知，用分解方法求超定方程组的最小二乘解只需求解上三角方程组RX=QTb四、GPS定位解算原理如果不考虑GPS接收机钟差，也不考虑信号传播过程的电离层延迟、对流层延迟和多径延迟误差，则GPS定位模型为222)()()(jjjjzzyyxx，(j=1，2，…，N)其中，(x，y，z)为接收机坐标，(xj，yj，zj)为卫星Pj的坐标；j是伪距（实际上可得到的距离观测量）。当已知接收机的概略位置r0=(x0，y0，z0)时，可以用牛顿迭代法结合最小二乘原理实现精确定位。对模型右边应用Taylor级数展开，并略去高次项，得到jjjjLbznymxl其中，0xxx，0yyy，0zzz，jjjRxxl~0，jjjRyym~0，jjjRzzn~0jjjRL~，202020)()()(~zzyyxxRjjjj设视界内的卫星数为N。将上述方程组写为矩阵形式，得到GX=L其中，TbzyxX][，L=[L1，L1，……，LN]T而111222111NNNnmlnmlnmlG用最小二乘法求解超定方程组，得正规方程TTGGXGL当为矩阵TGG非奇异时，其解为1()TTXGGGL得定位解xxx01，yyy01，zzz01如果概略坐标（x0，y0，z0）误差较大，则需要做迭代计算，以（x1，y1，z1）代替（x0，y0，z0）重复上面计算过程。下面数据来自于Appliedmathematicsandcomputation119(2001)21—34文章题目：AlternativealgorithmsfortheGPSstaticpositioningsolution作者：JohnB.Lundberg表1GPS卫星数据IDXYZGPS-116414028.668660383.61820932036.90724658975.31743GPS-216896800.648-18784061.365-7418318.85622964286.41228GPS-39339639.616-14514964.65820305107.16121338550.64536GPS-4-18335582.591-11640868.30515028599.07123606547.29359GPS-5-2077142.705-20987755.987-15879741.19624263298.50401GPS-6-4957166.885-23306741.03912039027.09620758264.10823GPS-717977519.820-13089823.31214331151.06521847468.81689GPS-89682727.508-24060519.4853985404.53020352077.19349xyzClockoffset96133.382910-5674076.36992740537.6613-36000.0000MATLAB仿真程序（gps0.m）%卫星的坐标，伪距xyz=[16414028.668,660383.618,20932036.90716896800.648,-18784061.365,-7418318.8569339639.616,-14514964.658,20305107.161-18335582.591,-11640868.305,15028599.071-2077142.705,-20987755.987,-15879741.196-4957166.885,-23306741.039,12039027.09617977519.820,-13089823.312,14331151.0659682727.508,-24060519.485,3985404.530];ro=[24658975.3174322964286.4122821338550.6453623606547.2935924263298.5040120758264.1082321847468.8168920352077.19349];P=[0;0;0];%设定位数据初始值即概略坐标error=10;k=0;tic%程序运行计时开始whileerror1e-8%设置定位精度forj=1:8Dj=xyz(j,:)-P';dj=norm(Dj);%计算接收机到第j颗卫星向量及其长度G(j,1:3)=Dj/dj;G(j,4)=-1;L(j,:)=dj-ro(j);%计算几何矩阵第j行及方程右端值endb=G\L;%解几何方程error=abs(b(1:3));%计算位值修正量绝对值P=P+b(1:3);k=k+1;%修正定位数据endtoc%程序运行计时结束formatlongeX=[P;b(4)]'%显示满足精度的定位数据Format计算结果elapsed_time=0.0160X=Columns1through39.613338291042342e+005-5.674076369901314e+0062.740537661301808e+006Column4-3.600000000248986e+004数值实验问题1．一个测绘员要测量出在某个基准点上三个山头A、B、C的高度，首先从基准点处观测，测得山头高度分别为x1=1237(ft)、x2=1914(ft)、x3=2417(ft)为了进一步确认测量数据。测绘员爬上第一个山头，测得第二个山头相对于第一个山头的高度为711(ft)，第三个山头相对于第一个山头的高度为1177(ft)。最后它再爬上第二个山头测得第三个山头相对于第二个山头的高度为475(ft)。用最小二乘法处理六个测量数据。2．数据拟合的最小二乘方法源于天文学中对行星或慧星这类天体的轨道计算。1795年，高斯在计算行星的椭圆轨道时提出并使用了这种方法，这一方法由勒让德于1805年首次公布。由开普列的研究成果，行星在其轨道平面上的运行轨迹是一个椭圆，而椭圆方程a1x2+2a2xy+a3y2+a4x+a5y+1=0需要由五个参数确定。原则上只要对行星的位置作5次观测就足以确定它的整个轨迹方程。但由于存在测量误差，由5次观测所确定的轨迹极不可靠，需要进行多次观测，用最小二乘法来消除误差，得到有关轨迹参数的更精确的值。最小二乘近似将几十次甚至上百次的观察所产生的高维空间问题降维到五参数的椭圆轨迹模型的五维空间处理。行星位置的10个观测点数据如下x1.020.950.870.770.670.560.440.300.1