超定方程最小二乘解
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
x1 + y2 = 的向量的线性组合问题。由于两个向量1，2 不构成三维空间的一组基，所以一般情况下 这一问题无解。而由向量1，2 张成的子空间 span{1，2}是一张平面，记为。则超定方 程组的最小二乘解实际上是求 X*，使 GX* 恰好等于 在平面上的投影。而最小二乘解所 对应的残差向量则垂直于向量 GX*。事实上，由正规方程组
X=( R T QTQR)-1(QR)Tb = R-1QTb 由此可知，用分解方法求超定方程组的最小二乘解只需求解上三角方程组
RX=QTb 四、GPS 定位解算原理
如果不考虑 GPS 接收机钟差，也不考虑信号传播过程的电离层延迟、对流层延迟和多 径延迟误差，则 GPS 定位模型为
j (x x j )2 ( y y j )2 (z z j )2 ，( j = 1，2，…，N )
GTGX GT L
当为矩阵 GT G 非奇异时，其解为
X (GTG)1GT L
得定位解
x1 x0 x ， y1 y0 y ， z1 z0 z
如果概略坐标（x0，y0，z0）误差较大，则需要做迭代计算，以（x1，y1，z1）代替（x0，y0，
z0）重复上面计算过程。
X=
Columns 1 through 3
9.613338291042342e+005 -5.674076369901314e+006 2.740537661301808e+006
Column 4
-3.600000000248986e+004
数值实验问题 1．一个测绘员要测量出在某个基准点上三个山头 A、B、C 的高度，首先从基准点处
下面数据来自于 Applied mathematics and computation 119(2001)21—34
文章题目：Alternative algorithms for the GPS static positioning solution
作者：John B. Lundberg
表 1 GPS 卫星数据
矩阵的每一行代表一个方程，m 行代表 m 个线性联立方程。 n 列代表 n 个变量。如果 m 是独立方程数，根据 m<n、m=n、m>n 确定方程是 ‘欠定’、‘适定’ 还是 ‘超定’。
超定方程组：方程个数大于未知量个数的方程组。 对于方程组 Ra=y，R 为 n×m 矩阵，如果 R 列满秩，且 n>m 超定方程一般是不存在解的矛盾方程。 例如，如果给定的三点不在一条直线上， 我们将无法得到这样一条直线，使得这条直
MATLAB 仿真程序（gps0.m）
%卫星的坐标，伪距
xyz=[16414028.668, 660383.618, 20932036.907
16896800.648, -18784061.365, -7418318.856
9339639.616, -14514964.658, 20305107.161
GTGX=GTb
得 GT (b –GX * ) = 0
rmin =b – GX*
上式的几何意义可解释为：最小二乘解的残
差向量与超定方程组的系数矩阵 G 的所有
列向量正交。从而
GX*
所以
(X*)T GT (b –GX * ) = 0 (GX *，b –GX * ) = 0
三、最小二乘解的两种方法 超定方程组的最小二乘解是指正规方程组
9682727.508, -24060519.485, 3985404.530];
ro=[24658975.31743
22964286.41228
21338550.64536
23606547.29359
24263298.50401
20758264.10823
21847468.81689
20352077.19349];
end
b=G\L;
%解几何方程
error=abs(b(1:3)); P=P+b(1:3);k=k+1;
%计算位值修正量绝对值 %修正定位数据
end
toc
%程序运行计时结束
format long e
X=[P;b(4)]'
%显示满足精度的定位数据
Format
计算结果
elapsed_time = 0.0160
其中， x x x0 ，y y y0 ， z z z0 ，
lj
x
j R~ j
x0
，mj
y
j R~ j
y
0
，nj
z
j R~ j
z0
Lj j R~ j ， R~ j (x j x0 )2 ( y j y0 )2 (z j z0 )2
||
b
GY
||22
||
(b
GX
*
)
G( X
*
Y
)
||
2 2
上式右端利用内积，得
||
b
GY
||
2 2
||
b
GX *
||
2 2
|| G( X *
Y)
||22 ||
b
GX *
||
2 2
从而有 || b – GY ||2 ≥ || b – GX*||2
等式仅当 Y=X*时成立。所以 X*是超定方程组 GX=b 的最小二乘解。 命题 2：如果 X*是超定方程组 GX=b 的最小二乘解，则 X*满足正规方程组 GTGX=GTb
-18335582.591, -11640868.305, 15028599.071
-2077142.705, -20987755.987, -15879741.196
-4957166.885, -23306741.039, 12039027.096
17977519.820, -13089823.312, 14331151.065
欠定方程组: 方程个数小于未知量个数的方程组。 对于方程组 Ra=y，R 为 n×m 矩阵，且 n<m。则方程组有无穷多组解，此时称方程组
为欠定方程组。 内点法和梯度投影法是目前解欠定方程组的常用方法。
超定方程组最小二乘解课程设计
最小二乘法广泛地应用于工程计算中，用最小二乘法消除（平滑）误差，用最小二乘法 从有噪声的数据中提取信号，从海量数据中找出数据变化的趋势，……。甚至利用简单函数 计算复杂函数的近似值，我们并不期望它的近似值多么精确（事实上很多时候也不用很精 确），尽管如此还是希望计算出的近似数据与原始数据之间有相似之处。如果从线性代数角 度来理解最小二乘法，实际上是将一个高维空间的向量投影到低维子空间所涉及的工作。
GPS-7
17977519.820 -13089823.312 14331151.065 21847468.81689
GPS-8
9682727.508
-24060519.485 3985404.530
20352077.19349
x
பைடு நூலகம்
y
z
Clock offset
96133.382910 -5674076.3699 2740537.6613 -36000.0000
其中，(x，y，z) 为接收机坐标，(x j，y j，z j)为卫星 Pj 的坐标； j 是伪距（实际上可得到
的距离观测量）。当已知接收机的概略位置 r0 = (x0，y0，z0)时，可以用牛顿迭代法结合最 小二乘原理实现精确定位。对模型右边应用 Taylor 级数展开，并略去高次项，得到
l jx m jy n jz b Lj
设视界内的卫星数为 N。将上述方程组写为矩阵形式，得到 GX=L
其中， X [x y z b]T ，L=[L1，L1，……，LN]T 而
l1 m1 n1 1
G
l
2
m2
n2 1
l
N
mN
nN
1
用最小二乘法求解超定方程组，得正规方程
GPS-4
-18335582.591 -11640868.305 15028599.071 23606547.29359
GPS-5
-2077142.705 -20987755.987 -15879741.196 24263298.50401
GPS-6
-4957166.885 -23306741.039 12039027.096 20758264.10823
证
由题设， || b GX *
||2
min
XRm
|| b GX
||2 ，利用 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)]+ t2 (Gv，Gv)
一、超定方程组的最小二乘解
当方程组 GX=b 的方程数多于未知数个数时，对应的系数矩阵 G 的行数大于列数，此 时方程组被称为是超定方程组。设 G=(giu)m×n，当 m>n 时即所谓的高矩阵，绝大多数情况下， 超定方程组没有古典意义下的解。超定方程组的最小二乘解是一种广义解，是指使残差 r = b – GX 的 2-范数达取极小值的解，即
由题设 t=0 是 g(t)的极小值点。由极值必要条件，得 g(0) 0 。即
(GX*，Gv) – (Gv，b)=0 将左端整理化简，便得
(Gv，GX* – b ) =0 利用内积性质，得
( v，GT(GX* – b ) )=0 由 v 的任意性，得 GT(GX* – b ) =0
二、最小二乘解的几何意义 首先考虑一个简单的超定方程组
2 4
11
3
1
5 2
x
y
3
6
该方程组的右端向量是三维向量，系数矩阵的每一列也是三维向量，但待求的未知向量却是
二维向量。将系数矩阵按列分块，G =[1，2]，记右端向量为。则方程组求解问题可表示 为求组合系数 x 和 y 使
||
b
GX
*
||2
min
XRm
||
b
GX
||2
该问题是一个优化问题。 命题 1：如果 X*是正规方程组 GTGX=GTb 的解，则 X*是超定方程组 GX=b 的最小二乘解 证 由题设可得，GT (b – GX*)=0。对任意 n 维向量 Y，显然有
(X* – Y)TGT (b – GX*)=0 考虑残差 2-范数平方，由
ID
X
Y
Z
GPS-1
16414028.668 660383.618
20932036.907 24658975.31743
GPS-2
16896800.648 -18784061.365 -7418318.856 22964286.41228
GPS-3
9339639.616
-14514964.658 20305107.161 21338550.64536
线同时经过给定这三个点。 也就是说给定的条件（限制）过于严格， 导致解不存在。在实 验数据处理和曲线拟合问题中，求解超定方程组非常普遍。比较常用的方法是最小二乘法。 形象的说，就是在无法完全满足给定的这些条件的情况下，求一个最接近的解。
曲线拟合的最小二乘法要解决的问题，实际上就是求以上超定方程组的最小二乘解的 问题。
GTGX=GTb 的解。如果系数矩阵(GTG)可逆，则正规方程组有唯一解。此时，最小二乘解可以形式地写 为如下形式
X=(GTG)-1GTb 两种常用的方法如下 1．用对称矩阵的三角分解法解正规方程组 GTGX = GTb；记 A=GTG，则 A 是对称矩阵，由 三角分解 A = L D LT，其中 L 是下三角矩阵，D 是对角矩阵。将这一算法写过三个过程： ①解下三角方程组：LY1 = GTb； ②解对角方程组：DY2 = Y1 ； ③解上三角方程组：LTY3 = Y2 2．用矩阵的 QR 分解直接求解超定方程组 由 QR 分解（正交三角分解）G=QR，其中 Q 是正交矩阵，R 是上三角矩阵。将 QR 分解代 入最小二乘解表达式中，得
P=[0;0;0];
%设定位数据初始值即概略坐标
error=10;k=0;
tic while error>1e-8
%程序运行计时开始 %设置定位精度
for j=1:8
Dj=xyz(j,:)-P'; dj=norm(Dj);
%计算接收机到第 j 颗卫星向量及其长度
G(j,1:3)=Dj/dj;G(j,4)=-1; L(j,:)=dj-ro(j); %计算几何矩阵第 j 行及方程右端值
X=( R T QTQR)-1(QR)Tb = R-1QTb 由此可知，用分解方法求超定方程组的最小二乘解只需求解上三角方程组
RX=QTb 四、GPS 定位解算原理
如果不考虑 GPS 接收机钟差，也不考虑信号传播过程的电离层延迟、对流层延迟和多 径延迟误差，则 GPS 定位模型为
j (x x j )2 ( y y j )2 (z z j )2 ，( j = 1，2，…，N )
GTGX GT L
当为矩阵 GT G 非奇异时，其解为
X (GTG)1GT L
得定位解
x1 x0 x ， y1 y0 y ， z1 z0 z
如果概略坐标（x0，y0，z0）误差较大，则需要做迭代计算，以（x1，y1，z1）代替（x0，y0，
z0）重复上面计算过程。
X=
Columns 1 through 3
9.613338291042342e+005 -5.674076369901314e+006 2.740537661301808e+006
Column 4
-3.600000000248986e+004
数值实验问题 1．一个测绘员要测量出在某个基准点上三个山头 A、B、C 的高度，首先从基准点处
下面数据来自于 Applied mathematics and computation 119(2001)21—34
文章题目：Alternative algorithms for the GPS static positioning solution
作者：John B. Lundberg
表 1 GPS 卫星数据
矩阵的每一行代表一个方程，m 行代表 m 个线性联立方程。 n 列代表 n 个变量。如果 m 是独立方程数，根据 m<n、m=n、m>n 确定方程是 ‘欠定’、‘适定’ 还是 ‘超定’。
超定方程组：方程个数大于未知量个数的方程组。 对于方程组 Ra=y，R 为 n×m 矩阵，如果 R 列满秩，且 n>m 超定方程一般是不存在解的矛盾方程。 例如，如果给定的三点不在一条直线上， 我们将无法得到这样一条直线，使得这条直
MATLAB 仿真程序（gps0.m）
%卫星的坐标，伪距
xyz=[16414028.668, 660383.618, 20932036.907
16896800.648, -18784061.365, -7418318.856
9339639.616, -14514964.658, 20305107.161
GTGX=GTb
得 GT (b –GX * ) = 0
rmin =b – GX*
上式的几何意义可解释为：最小二乘解的残
差向量与超定方程组的系数矩阵 G 的所有
列向量正交。从而
GX*
所以
(X*)T GT (b –GX * ) = 0 (GX *，b –GX * ) = 0
三、最小二乘解的两种方法 超定方程组的最小二乘解是指正规方程组
9682727.508, -24060519.485, 3985404.530];
ro=[24658975.31743
22964286.41228
21338550.64536
23606547.29359
24263298.50401
20758264.10823
21847468.81689
20352077.19349];
end
b=G\L;
%解几何方程
error=abs(b(1:3)); P=P+b(1:3);k=k+1;
%计算位值修正量绝对值 %修正定位数据
end
toc
%程序运行计时结束
format long e
X=[P;b(4)]'
%显示满足精度的定位数据
Format
计算结果
elapsed_time = 0.0160
其中， x x x0 ，y y y0 ， z z z0 ，
lj
x
j R~ j
x0
，mj
y
j R~ j
y
0
，nj
z
j R~ j
z0
Lj j R~ j ， R~ j (x j x0 )2 ( y j y0 )2 (z j z0 )2
||
b
GY
||22
||
(b
GX
*
)
G( X
*
Y
)
||
2 2
上式右端利用内积，得
||
b
GY
||
2 2
||
b
GX *
||
2 2
|| G( X *
Y)
||22 ||
b
GX *
||
2 2
从而有 || b – GY ||2 ≥ || b – GX*||2
等式仅当 Y=X*时成立。所以 X*是超定方程组 GX=b 的最小二乘解。 命题 2：如果 X*是超定方程组 GX=b 的最小二乘解，则 X*满足正规方程组 GTGX=GTb
-18335582.591, -11640868.305, 15028599.071
-2077142.705, -20987755.987, -15879741.196
-4957166.885, -23306741.039, 12039027.096
17977519.820, -13089823.312, 14331151.065
欠定方程组: 方程个数小于未知量个数的方程组。 对于方程组 Ra=y，R 为 n×m 矩阵，且 n<m。则方程组有无穷多组解，此时称方程组
为欠定方程组。 内点法和梯度投影法是目前解欠定方程组的常用方法。
超定方程组最小二乘解课程设计
最小二乘法广泛地应用于工程计算中，用最小二乘法消除（平滑）误差，用最小二乘法 从有噪声的数据中提取信号，从海量数据中找出数据变化的趋势，……。甚至利用简单函数 计算复杂函数的近似值，我们并不期望它的近似值多么精确（事实上很多时候也不用很精 确），尽管如此还是希望计算出的近似数据与原始数据之间有相似之处。如果从线性代数角 度来理解最小二乘法，实际上是将一个高维空间的向量投影到低维子空间所涉及的工作。
GPS-7
17977519.820 -13089823.312 14331151.065 21847468.81689
GPS-8
9682727.508
-24060519.485 3985404.530
20352077.19349
x
பைடு நூலகம்
y
z
Clock offset
96133.382910 -5674076.3699 2740537.6613 -36000.0000
其中，(x，y，z) 为接收机坐标，(x j，y j，z j)为卫星 Pj 的坐标； j 是伪距（实际上可得到
的距离观测量）。当已知接收机的概略位置 r0 = (x0，y0，z0)时，可以用牛顿迭代法结合最 小二乘原理实现精确定位。对模型右边应用 Taylor 级数展开，并略去高次项，得到
l jx m jy n jz b Lj
设视界内的卫星数为 N。将上述方程组写为矩阵形式，得到 GX=L
其中， X [x y z b]T ，L=[L1，L1，……，LN]T 而
l1 m1 n1 1
G
l
2
m2
n2 1
l
N
mN
nN
1
用最小二乘法求解超定方程组，得正规方程
GPS-4
-18335582.591 -11640868.305 15028599.071 23606547.29359
GPS-5
-2077142.705 -20987755.987 -15879741.196 24263298.50401
GPS-6
-4957166.885 -23306741.039 12039027.096 20758264.10823
证
由题设， || b GX *
||2
min
XRm
|| b GX
||2 ，利用 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)]+ t2 (Gv，Gv)
一、超定方程组的最小二乘解
当方程组 GX=b 的方程数多于未知数个数时，对应的系数矩阵 G 的行数大于列数，此 时方程组被称为是超定方程组。设 G=(giu)m×n，当 m>n 时即所谓的高矩阵，绝大多数情况下， 超定方程组没有古典意义下的解。超定方程组的最小二乘解是一种广义解，是指使残差 r = b – GX 的 2-范数达取极小值的解，即
由题设 t=0 是 g(t)的极小值点。由极值必要条件，得 g(0) 0 。即
(GX*，Gv) – (Gv，b)=0 将左端整理化简，便得
(Gv，GX* – b ) =0 利用内积性质，得
( v，GT(GX* – b ) )=0 由 v 的任意性，得 GT(GX* – b ) =0
二、最小二乘解的几何意义 首先考虑一个简单的超定方程组
2 4
11
3
1
5 2
x
y
3
6
该方程组的右端向量是三维向量，系数矩阵的每一列也是三维向量，但待求的未知向量却是
二维向量。将系数矩阵按列分块，G =[1，2]，记右端向量为。则方程组求解问题可表示 为求组合系数 x 和 y 使
||
b
GX
*
||2
min
XRm
||
b
GX
||2
该问题是一个优化问题。 命题 1：如果 X*是正规方程组 GTGX=GTb 的解，则 X*是超定方程组 GX=b 的最小二乘解 证 由题设可得，GT (b – GX*)=0。对任意 n 维向量 Y，显然有
(X* – Y)TGT (b – GX*)=0 考虑残差 2-范数平方，由
ID
X
Y
Z
GPS-1
16414028.668 660383.618
20932036.907 24658975.31743
GPS-2
16896800.648 -18784061.365 -7418318.856 22964286.41228
GPS-3
9339639.616
-14514964.658 20305107.161 21338550.64536
线同时经过给定这三个点。 也就是说给定的条件（限制）过于严格， 导致解不存在。在实 验数据处理和曲线拟合问题中，求解超定方程组非常普遍。比较常用的方法是最小二乘法。 形象的说，就是在无法完全满足给定的这些条件的情况下，求一个最接近的解。
曲线拟合的最小二乘法要解决的问题，实际上就是求以上超定方程组的最小二乘解的 问题。
GTGX=GTb 的解。如果系数矩阵(GTG)可逆，则正规方程组有唯一解。此时，最小二乘解可以形式地写 为如下形式
X=(GTG)-1GTb 两种常用的方法如下 1．用对称矩阵的三角分解法解正规方程组 GTGX = GTb；记 A=GTG，则 A 是对称矩阵，由 三角分解 A = L D LT，其中 L 是下三角矩阵，D 是对角矩阵。将这一算法写过三个过程： ①解下三角方程组：LY1 = GTb； ②解对角方程组：DY2 = Y1 ； ③解上三角方程组：LTY3 = Y2 2．用矩阵的 QR 分解直接求解超定方程组 由 QR 分解（正交三角分解）G=QR，其中 Q 是正交矩阵，R 是上三角矩阵。将 QR 分解代 入最小二乘解表达式中，得
P=[0;0;0];
%设定位数据初始值即概略坐标
error=10;k=0;
tic while error>1e-8
%程序运行计时开始 %设置定位精度
for j=1:8
Dj=xyz(j,:)-P'; dj=norm(Dj);
%计算接收机到第 j 颗卫星向量及其长度
G(j,1:3)=Dj/dj;G(j,4)=-1; L(j,:)=dj-ro(j); %计算几何矩阵第 j 行及方程右端值