超定方程组最小二乘解课程设计最小二乘法广泛地应用于工程计算中，用最小二乘法消除（平滑）误差，用最小二乘法从有噪声的数据中提取信号，从海量数据中找出数据变化的趋势，……。

甚至利用简单函数计算复杂函数的近似值，我们并不期望它的近似值多么精确（事实上很多时候也不用很精确），尽管如此还是希望计算出的近似数据与原始数据之间有相似之处。

如果从线性代数角度来理解最小二乘法，实际上是将一个高维空间的向量投影到低维子空间所涉及的工作。

一、超定方程组的最小二乘解当方程组GX=b 的方程数多于未知数个数时，对应的系数矩阵G 的行数大于列数，此时方程组被称为是超定方程组。

设G=(g iu )m ×n ，当m>n 时即所谓的高矩阵，绝大多数情况下，超定方程组没有古典意义下的解。

超定方程组的最小二乘解是一种广义解，是指使残差r = b – GX 的2-范数达取极小值的解，即22*||||min ||||GX b GX b mRX -=-∈ 该问题是一个优化问题。

命题1：如果X *是正规方程组G T GX=G Tb 的解，则X *是超定方程组GX=b 的最小二乘解证 由题设可得，G T(b – GX *)=0。

对任意n 维向量Y ，显然有(X * – Y )T G T (b – GX *)=0考虑残差2-范数平方，由22**22||)()(||||||Y X G GX b GY b -+-=-上式右端利用内积，得22*22*22*22||||||)(||||||||||GX b Y X G GX b GY b -≥-+-=-从而有|| b – GY ||2 ≥ || b – GX *||2等式仅当Y =X *时成立。

所以X *是超定方程组GX=b 的最小二乘解。

命题2：如果X *是超定方程组GX=b 的最小二乘解，则X *满足正规方程组G T GX=G Tb证 由题设，22*||||min ||||GX b GX b mRX -=-∈，利用2-范数与内积关系，知X *是下面二次函数的极小值点ϕ(X ) = (GX ，GX ) – 2(GX ，b ) + (b ，b )取任意n 维向量v ，对任意实数t ，构造一元函数g (t ) = ϕ(X * + t v )显然， g (t ) 是关于变量t 的二次函数g (t ) = (G (X * + t v )，G (X * + t v )) – 2(G (X * + t v )，b ) + (b ，b )= g (0) + 2t [(GX *，Gv ) – (Gv ，b )]+ t 2 (Gv ，Gv )由题设t =0是g (t )的极小值点。

由极值必要条件，得0)0(='g 。

即(GX *，Gv ) – (Gv ，b )=0将左端整理化简，便得(Gv ，GX * – b ) =0利用内积性质，得( v ，G T (GX * – b ) )=0由v 的任意性，得G T(GX * – b ) =0二、最小二乘解的几何意义 首先考虑一个简单的超定方程组⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-6311215342y x该方程组的右端向量是三维向量，系数矩阵的每一列也是三维向量，但待求的未知向量却是二维向量。

将系数矩阵按列分块，G =[α1，α2]，记右端向量为β。

则方程组求解问题可表示为求组合系数x 和y 使x α1 + y α2 = β的向量的线性组合问题。

由于两个向量α1，α2不构成三维空间的一组基，所以一般情况下这一问题无解。

而由向量α1，α2张成的子空间span{α1，α2}是一张平面，记为π。

则超定方程组的最小二乘解实际上是求X *，使GX * 恰好等于β 在平面π上的投影。

而最小二乘解所对应的残差向量则垂直于向量GX *。

事实上，由正规方程组G T GX=G T b得G T(b –GX *) = 0上式的几何意义可解释为：最小二乘解的残差向量与超定方程组的系数矩阵G 的所有列向量正交。

从而(X *)T G T (b –GX * ) = 0所以(GX *，b –GX * ) = 0三、最小二乘解的两种方法超定方程组的最小二乘解是指正规方程组G T GX=G T b的解。

如果系数矩阵(G TG )可逆，则正规方程组有唯一解。

此时，最小二乘解可以形式地写为如下形式X=(G T G )-1G T b两种常用的方法如下1．用对称矩阵的三角分解法解正规方程组G TGX = G Tb ；记A=G TG ，则A 是对称矩阵，由三角分解A = L D L T，其中L 是下三角矩阵，D 是对角矩阵。

将这一算法写过三个过程：①解下三角方程组：LY 1 = G Tb ； ②解对角方程组：DY 2 = Y 1 ；③解上三角方程组：L TY 3 = Y 22．用矩阵的QR 分解直接求解超定方程组由QR 分解（正交三角分解）G=QR ，其中Q 是正交矩阵，R 是上三角矩阵。

将QR 分解代入最小二乘解表达式中，得X=( R T Q T QR )-1(QR )T b = R -1Q T b由此可知，用分解方法求超定方程组的最小二乘解只需求解上三角方程组RX=Q T b四、GPS 定位解算原理如果不考虑GPS 接收机钟差，也不考虑信号传播过程的电离层延迟、对流层延迟和多径延迟误差，则GPS 定位模型为222)()()(j j j j z z y y x x -+-+-=ρ，( j = 1，2，…，N )其中，(x ，y ，z ) 为接收机坐标，(x j，y j，z j)为卫星P j 的坐标；j ρ是伪距（实际上可得到的距离观测量）。

当已知接收机的概略位置r 0 = (x 0，y 0，z 0)时，可以用牛顿迭代法结合最小二乘原理实现精确定位。

对模型右边应用Taylor 级数展开，并略去高次项，得到j j j j L b z n y m x l =+++δδδ其中，0x x x -=δ，0y y y -=δ，0z z z -=δ，j j j R x x l ~0-=，j j j R y y m ~0-=，j jj Rz z n ~0-=j j j R L ~-=ρ，202020)()()(~z z y y x x R j j j j -+-+-=设视界内的卫星数为N 。

将上述方程组写为矩阵形式，得到GX=L其中，T b z y xX ][δδδ=，L=[L 1，L 1，……，L N ]T 而⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡---=111222111NN Nn m l n m l n m l G 用最小二乘法求解超定方程组，得正规方程T T G GX G L =当为矩阵TG G 非奇异时，其解为1()T T X G G G L -=得定位解x x x δ+=01，y y y δ+=01，z z z δ+=01如果概略坐标（x 0，y 0，z 0）误差较大，则需要做迭代计算，以（x 1，y 1，z 1）代替（x 0，y 0，z 0）重复上面计算过程。

下面数据来自于Applied mathematics and computation 119(2001)21—34 文章题目：Alternative algorithms for the GPS static positioning solution作者：John B. LundbergMATLAB仿真程序（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 %程序运行计时开始while error>1e-8 %设置定位精度for j=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 %程序运行计时结束format long eX=[P;b(4)]' %显示满足精度的定位数据Format计算结果elapsed_time = 0.0160X =Columns 1 through 39.613338291042342e+005 -5.674076369901314e+006 2.740537661301808e+006 Column 4-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 + 2a2 xy +a3y2 + a4x + a5y + 1=0需要由五个参数确定。

原则上只要对行星的位置作5次观测就足以确定它的整个轨迹方程。

但由于存在测量误差，由5次观测所确定的轨迹极不可靠，需要进行多次观测，用最小二乘法来消除误差，得到有关轨迹参数的更精确的值。