总体最小二乘的迭代解法

第35卷第11期2010年11月武汉大学学报 信息科学版
G eomatics and Infor matio n Science o f Wuhan U niv ersity V ol.35N o.11
N ov.2010
收稿日期:2010 09 15。

项目来源:国家自然科学基金资助项目(40874010);江西省自然科学基金资助项目(2007GZC0474,2008GQC 0001,2008GZS0041);
武汉大学地球空间环境与大地测量教育部重点实验室开放研究基金资助项目(080104);现代工程测量国家测绘局重点实验室开放研究基金资助项目(T JES0802);数字国土江西省重点实验室开放研究基金资助项目(DLLJ200506)。

文章编号:1671 8860(2010)11 1351 04文献标志码:A
总体最小二乘的迭代解法
鲁铁定
1,2
周世健
1,3
(1 东华理工大学地球科学与测绘工程学院,抚州市学府路56号,344000)
(2 武汉大学测绘学院,武汉市珞喻路129号,430079)
(3 江西省科学院,南昌市上访路,330029)
摘 要:针对总体最小二乘解算问题,应用测量平差中的间接平差原理推导了总体最小二乘的迭代逼近解算公式,通过与奇异值分解法进行比较,得出两种解算方法具有等价性。

实验数据分析验证了算法的有效性。

关键词:总体最小二乘;奇异值分解;迭代算法;测量平差中图法分类号:P207.2
在经典的高斯 马尔可夫模型中,通常情况下是假设偶然误差仅存在于观测向量中,而系数矩阵不含误差。

然而由于观测条件的限制,观测向量、系数矩阵都有可能存在误差。

文献[1]讨论了共线方程解算外方位元素的解算问题[1]和坐标转换问题[2],认为得到的计算结果精度更高。

在应用总体最小二乘法解决顾及系数矩阵误差的平差解算中,基本上是用Golub 和Lan Loan 奇异值分解解法[3]。

在解算算法研究方面,文献[4 7]探讨了总体最小二乘的解算以及求解问题,并用于解决坐标转换等问题。

文献[8]探讨了线性回归的总体最小二乘解算算法,得到系数矩阵误差的总体最小二乘法可以解决最小二乘解算在不同方向拟合的不一致问题。

因此,研究总体最小二乘的解算及精度评定对于解决测量实践中遇到的系数矩阵含有误差的平差问题具有现实意义。

本文基于测量平差中常用的间接平差原理,简化了过程,得到了与文献[7]相同的解算公式。

1 总体最小二乘的奇异值分解
在总体最小二乘中,函数模型为:
(A + A )X =L + L
(1)
误差向量的均值向量和协方差矩阵为[7]:
L vec ( A )~
00, L 00
A I n
=
2
0I
n 00
I m I n
(2)
式中,v ec ( )表示矩阵的拉直变换。

将式(1)表示为误差方程式的形式为:
(A +E A )X =L +e
(3
)
可以改写为:
A +E A L +e X -1=
A L +E X -1
=0
(4)式中,E =E A
e 。

用其表示的限制约束条件
为[4]:
tr (E E T )=tr (E A E T A +ee T
)=
vec (E A )T vec (E A )+e T e =min
(5)
对增广矩阵A L 进行奇异值分解[4]
:A L =
U 1
U 2
1
n -(m+)
V T =U 1 V T (6)
其中,
U 1=
U 11U 12m
1
(7)
武汉大学学报 信息科学版2010年11月
V T
=
1
00
2
m
1
V 11V 12V 21
V 22
m
1
T
m 1
(8)
由Eckart Yo ung Mirsky 矩阵逼近理论[9]
知,未知参数X 的估值为[10]:
X ^=-V 12 V -1
22=
-1v m+1,m +1
V 12
(9)
由于v i 为矩阵A L
T
A L 的特征向量,所以
奇异值与特征向量满足关系
[10]
:
A T
A A
T
L
L T
A L T
L
X ^-1= 2
m +1
X ^-1
进一步可以计算得到单位权方差以及参数的协方差矩阵[7]:
^e T ^e +(vec E ^A )T v ec E ^A =tr (^e ^e T +
E ^A E ^T A )=tr (E ^ E ^T )= 2
m +1
(10) 20(TLS )= 2
m+1/(n-m)
(11)
D(X ^) 20
(N -
2m+1
I m )
-1
N (N -
2m +1
I m )
-1
(12)
其中,N =A T A 。

2 总体最小二乘迭代解算
2.1 解算原理
将误差方程式(3)可以表示为[7]:
L =AX -e +E A X =AX +[-I ,X T
I ]
e e A
T
(13)
e
vec (E A )
~
00,
20I n 00
I m I n
(14)
式中,e A =vec (E A )。

式(13)表示为矩阵形式为:
e =AX -L (15)
e =[I n -X T
I n ]
e e A
=e -E A X (16)
由协因数传播律容易得到协因数阵的关系为
[11]
:
Q e =[I n -X T
I n ]
I n 00
I m I n
I n
-X I n
=
I n +X T X I n =(1+X T X )I n
(17) 总体最小二乘的目标函数式(5)等价于: e T
Q -1
e e =m in (18)
求解目标函数的自由极值得:
!e ^ T Q -1 e e ^ !X ^=2e ^ T Q -1 e A -2e ^ T e ^ X ^T
(1+X ^T X ^)2
=0
(19)
从式(17)可以看出,矩阵Q
e 为X 的函数,所以与最小二乘不同的是其导数增加了第二项,式(19)
转置整理得:
A T
Q -1
e e
^ -X ^e ^ T e ^ (1+X ^T
X ^)2
=0
(20)
将式(15)和式(16)代入式(20)得:
A T AX ^-A T L =X ^e ^ T e ^ /(1+X ^T
X ^)=
X ^(AX ^-L )T (AX ^-L )/(1+X ^T
X ^)=X ^^!
(21)
式中,^!!=(L -AX ^)T (L -AX ^)/(1+X ^T X ^)。

可以看出,式(21)与文献[7]中的结果相同。

上式可进一步表示为:
X ^=(A T A )-1(A T L +X ^^v )
(22)
如果在对式(18)求自由极值时未考虑Q -1
e 为
X 的函数,将其作为常数进行推导,可得:
A T Q -1 e AX ^-A T Q -1
e L =0
(23)
将式(15)代入式(23),整理得到参数的估值为:X ^=(A T Q -1 e A )-1A T Q -1 e L =[A T
(1+
X T X )I n A ]-1A T (1+X T X )L =(A T A )-1A T L
(24)
分析可得,式(24)的结果和最小二乘结果完全相同。

在线性回归和顾及自变量误差的参数估计中,文献[11 13]将自变量误差直接引入,并舍去二次项进行近似,然后按照最小二乘法进行计算,这样计算出来的参数估值和未顾及自变量误差的结果完全一致,这和文献[8]的结论也一致。

所以如果在解算中将Q e 近似作为常数矩阵,会导致结果与最小二乘结果一致,在公式推导中,应同时将Q e 作为参数的函数矩阵。

单位权方差以及由式(22)可以得到参数的协方差矩阵为:
20(T LS )=^v /(n -m)
(25)
D(X ^) 2
0(N -^!I m )
-1
N (N -^v I m )-1
(26)
2.2 迭代解算步骤
式(22)可以通过迭代解算得到参数估值,计算步骤如下:
1)^!=0,X ^(1)=N -1c ,N =A T A ,c =A T L
2)^!(i)
=(L -A X ^(i)
)T
(L -AX ^(i)
)/
(1+(X ^(i))T X ^(i)
)3)X ^(i +1)
=N -1
(c +X ^(i)
^!(i)
)4)当X ^(i +1)-X ^(i)<∀时,计算结束。

2.3 迭代解算与奇异值分解的关系
将解算结果代入式(15)得:
e ^ =[I n ,-X ^T I n ]
^e ^
e A =
^e -E ^A X ^=AX ^-L (27)
根据总体最小二乘准则,可以将式(27)看作条件平差问题,应用条件平差原理可得其法方程式为:
I n
-X ^T I
I n -X ^ I
K =e ^ =AX ^-L
1352
第35卷第11期鲁铁定等:总体最小二乘的迭代解法
解算法方程式可得联系数向量为:
K=(1+X^T X^)-1e=(1+X^T X^)-1(AX^-L)
进而可得系数矩阵和观测量的改正数向量为:
^e ^e A =
I n
-X^ I
K=
(1+X^T X^)-1(AX^-L)
-X^ I n(1+X^T X^)-1(AX^-L)
(28)将式(28)表示为:
^e=(1+X^T X^)-1(AX^-L)(29) ^e A=(-X^ I n)(1+X^T X^)-1(A X^-L)(30)由Kro necker积的性质易知,式(30)可表示为: E^A=-(1+X^T X^)-1(AX^-L)X^T(31)对^!进一步分析整理可得[7]:
^!=(L-A X^)T(L-A X^)/(1+X^T X^)=
[L T(L-AX^)+X^T(X^^!)]/(1+X^T X^)=
L T L-(A T L)T X^(32)联系式(22)和式(32),可表示为矩阵形式为[7]:
A T A A T L L T A L T L
X^
-1
=^!
X^
-1
,^!=!min∀0(33)
上式为特征值问题[10],参数的总体最小二乘解为最小特征值对应的特征向量[14],其最小特征值为^!=^e T^e+^e T A^e A=m in,这和Golub、Lan Lo an提出的奇异值分解方法完全一致。

由式(10)和式(32)可以看出,^!=^e T^e+^e T A^e A=2m+1,所以其精度估计公式一致。

3 算例及分析
以文献[15]例5 2的数据为样本观测值,共计25个点,假设回归直线为y=a+bx或x=c+ dy。

表1给出了所有的样本点数据,下面以两种拟合方案进行解算。

表1 观测样本值
T ab.1 V alue of the O bserv atio ns
编号Y X编号Y X
110.9835.3149.5739.1
211.1329.71510.9446.8
312.5130.8169.5848.5
48.4058.81710.0959.3
59.2761.4188.1170.0
68.7371.319 6.8370.0
7 6.3674.4208.8874.5
88.5076.6217.6872.1
97.8270.7228.4758.1
109.1457.5238.8644.6
118.2446.42410.3833.4
1212.1928.92511.0828.6
1311.8828.1 方案1 以x为自变量、y为因变量建立回归方程y=a+bx。

#最小二乘法:应用最小二乘平差解算得到回归方程为y=13.6284-0.0799x。

∃SVD分解方法:对增广矩阵进行分解得到正交矩阵V,根据前面公式可得:
^a
^b
=-1
v m+1,m+1
V12=
14.1952
-0.0897
得回归方程为y=14.1952-0.0897x。

%迭代逼近法:根据最小二乘计算得到参数的近似值!(1)=0.09737,进一步可得^a^b T= 14.1951-0.0897T,则得回归方程为y= 14.1952-0.0897x。

方案2 以y为自变量、x为因变量建立回归方程x=c+dy。

分别应用与方案1相同的方法,并将回归方程表示为和方案1相同的形式,则有:#最小二乘法:y=15.3022-0.1117x;∃SVD分解方法:y=14.1951-0.0897x;%迭代逼近法:y=14.1948-0.0897x。

从方案1和方案2的解算结果可以看出,用不同方法建立回归方程,因为最小二乘解算所建立的方程解算过程中没有顾及自变量误差的影响,因此导致了最小二乘解得到的回归方程式不一致,而用总体最小二乘解,由于同时顾及了自变量误差的影响,所以得到的回归方程具有一致性,解算结果检验了算法的正确性。

4 结 语
奇异值分解算法是一种常用的总体最小二乘计算方法,应用测量平差中的间接平差原理推导出了总体最小二乘的迭代逼近计算公式,通过与奇异值分解算法的比较阐述了两者在解算方法和精度估计方面的等价性。

从迭代计算的步骤可以看出,迭代解算公式类似于测量平差中最小二乘平差解算,方便使用。

推导过程和文献[7]中的推导过程比较,过程简洁,更容易理解。

本文考虑观测值向量和系数矩阵为等精度情形,如果不是等精度情况,协因数阵Q e的形式比较复杂,所以矩阵Q e对X的求导形式也比较复杂。

关于不等精度情况的解算公式、不等精度和等精度情况解算结果的影响大小是后续要进行的研究。

参 考 文 献
[1] 陈义,陆珏,郑波.总体最小二乘方法在空间后方交
会中的应用[J].武汉大学学报 信息科学版,
1353
武汉大学学报 信息科学版2010年11月
2008,33(12):1271 1274
[2] 陆珏,陈义,郑波.总体最小二乘方法在三维坐标转
换中的应用[J].大地测量与地球动力学,2008,l28
(5):77 81
[3] G olub G H,L an L oan F C.A n Analy sis o f t he
T otal L east Squar es Pro blem[J].SIA M Jour nal on
Numerical A nalysis,1980,17(6):883 893
[4] Schaffrin B,Felus Y A.On t he M ultiv ariat e T o tal
L east squares A ppr oach to Empir ical Co or dinate
T ransfor mations[J].T hr ee A lg or ithms J G eod,
2008,82:373 383
[5] Schaffr in B,Felus A Y.M ult ivar iate T ota l Least
squar es A djustment for Em pir ical Affine
T ransfor mations[C].T he6th H ot ine M ar ussi Sy mpo sium for T heor etical and Co mputational
G eodesy,Spring er,Ber lin,2007
[6] Schaffrin B,L ee I P,Felus Y A,et al.T o tal
L east squares(T LS)fo r Geo detic Straight line and
Plane A djustment[J].Boll Geo d Sci Affini,2006,
65(3):141 168
[7] Schaffrin B.A N ote on Co nstr ained T ot al Least
squar es Estimation[J].L inear A lg ebr a A ppl,2006,
417(1):245 258
[8] 鲁铁定,陶本藻,周世健.基于整体最小二乘法的线
性回归建模和解法[J].武汉大学学报 信息科学
版,2008,33(5):504 507
[9] Eckart C,Y oung G.T he A ppro ximation of O ne
M atrix by A no ther of L o wer Rank[J].
Psychometr ika,1936,1(3):211 218
[10]Van H uffel S,V andew alle J.T he T o tal Least
squar es P ro blem Computat ional A spects and
Analysis[M].Philadelphia:So ciety for Industrial
and A pplied M athematics,1991
[11]周世健,鲁铁定.双变量线性回归的解算[J].江西科
学,2008,26(1):109 111
[12]刘大杰,史文中,童小华,等.GIS空间数据的精度
分析与质量控制[M].上海:上海科学技术文献出
版社,1999
[13]王乐洋,许才军,鲁铁定.边长变化反演应变参数的
总体最小二乘方法[J].武汉大学学报 信息科学
版,2010,35(2):181 184
[14]Str ang G.Linear Algebra and Its A pplicat ions[M].
3rd ed.San D iego:Har court Br ace Jovanov ich,
1988
[15]腾素珍,冯敬海.数理统计学[M].大连:大连理工大
学出版社,2005
第一作者简介:鲁铁定,副教授,博士生,主要从事测绘数据处理理论和大地测量研究。

E mail:tdlu@
An Iterative Algorithm for Total Least Squares Estimation
L U T ieding1,2 ZH OU Shij ian1,3
(1 School of Geoscience and Surveying En gineering,East China Ins titute of Technology,56Xu efu Road,Fuzh ou344000,Chin a)
(2 Sch ool of Geod esy an d Geom atics,Wuhan University,129L uoyu Road,Wu han430079,Chin a)
(3 Jiangxi Acad emy of Sciences,Shangfan g Road,Nan chang,330029,Ch ina)
Abstract:Total least squares(TLS)appro ach aim s at estimating a matrix o f parameters from
a linear model w hen there are error s in both the o bserv ation vecto r L and the data m atrix A.
T he authors derived an iterative alg orithm to so lve the TLS pr oblem by using the principle o f indirect adjustm pared w ith the m ethod based on singular v alue decom positio n,the iter ative alg orithm coincides w ith the SVD algor ithm.T he calculated example has proved that the iterative alg orithm is validity and rationality.
Key w ords:TLS;sing ular value deco mposition;iter ative alg orithm;surveying adjustm ent
About the first author:LU Tieding,associate professor,P h.D candidate,m ajors in surveying data processing and geodesy.
E mail:tdlu@
1354。