多元总体最小二乘问题的牛顿解法

王乐洋，赵英文，陈晓勇，臧德彦

1. 东华理工大学测绘工程学院，江西 南昌 330013; 2. 流域生态与地理环境监测国家测绘地理信息局重点实验室，江西 南昌 330013; 3. 江西省数字国土重点实验室，江西 南昌 330013



多元总体最小二乘问题的牛顿解法

王乐洋1,2,3，赵英文1，陈晓勇1,2,3，臧德彦1,2,3

1. 东华理工大学测绘工程学院，江西 南昌 330013; 2. 流域生态与地理环境监测国家测绘地理信息局重点实验室，江西 南昌 330013; 3. 江西省数字国土重点实验室，江西 南昌 330013

transformation

Foundation support： The National Natural Science Foundation of China (Nos. 41204003; 41161069; 41304020; 41464001); The National Department Public Benefit Research Foundation (Surveying，Mapping and Geoinformation)(No. 201512026); The Natural Science Foundation of Jiangxi Province (No. 20151BAB203042); The Science and Technology Project of the Education Department of Jiangxi Province (Nos. GJJ150595; KJLD12077; KJLD14049); The Found of Key Laboratory of Watershed Ecology and Geographical Environment Monitoring (No. WE2015005); The Scientific Research Foundation of ECIT(No. DHBK201113); Innovation Fund Designated for Graduate Students of Jiangxi Province (Nos. YC2015-S266; YC2015-S267); Innovation Fund Designated for Graduate Students of ECIT(No. DHYC-2015005); The Found of the Key Laboratory of Mapping from Space, NASG (No.K201502)

摘要：为提高多元总体最小二乘问题参数估值的解算效率，推导了基于牛顿法的多元加权总体最小二乘算法；分析比较了基于牛顿法的多元加权总体最小二乘解和基于拉格朗日乘数法多元加权总体最小二乘解之间的关系，根据协因数传播律给出了多元总体最小二乘平差的16种协因数阵的近似计算公式。新算法能够解决观测矩阵和系数矩阵元素具有相关性的问题，并且可以把观测矩阵和系数矩阵的随机元素和常数元素纳入到一个协因数阵中进行处理。算例结果表明，本文提出的多元总体最小二乘问题的牛顿解法可行且收敛速度更快。

关键词：多元总体最小二乘；牛顿法；协因数传播律；仿射变换

总体最小二乘(total least squares，TLS)方法是近30多年来发展起来的一种能同时顾及观测值误差和模型系数矩阵误差的数学方法，其理论及应用研究是目前国内外研究的热点问题[1]。文献[2]首次提出总体最小二乘概念，随着变量误差模型(errors-in-variables，EIV)由单变量模型扩展到多元变量模型[3]，总体最小二乘也相应扩展到多元总体最小二乘(multivariate TLS，MTLS)[4-5]。某些模型的系数矩阵会因未知参数写成向量形式含有重复出现的随机元素和大量的零元素；多元总体最小二乘将参数向量和观测向量改写成矩阵形式，避免了这种现象，进而削弱了系数矩阵中行列之间的相关性，而且使得系数矩阵的构造更加简单。

在国外相关研究中，文献[5—7]认为观测矩阵和系数矩阵中的随机元素是等精度的，给出了基于奇异值分解法和拉格朗日乘数法的多元总体最小二乘算法；文献[8]讨论了基于拉格朗日乘数法的多元加权总体最小二乘算法，将随机元素的方差阵分解成两个矩阵取直积的形式，但是当观测数据含有两个以上的方差分量或是观测数据含有多类观测值时，随机变量的方差协方差阵就很难写成或是不能写成两个矩阵取直积的形式；文献[9]推导了基于奇异值分解法的多元加权总体最小二乘算法，但是这种算法只使用了行尺度权阵，其权阵不具有一般性。以上文献的所有算法都不适用于系数矩阵与观测矩阵中随机元素具有相关性的情况。文献[10]给出了基于拉格朗日乘数法的多元加权总体最小二乘算法，其权阵具有一般性，但没有具体给出系数矩阵含有常数元素的解法以及没有解决算法的精度评定问题，这种算法的解算效率也还需要验证。在国内相关研究中，文献[11]将基于奇异值分解法的多元总体最小二乘算法应用在三维坐标转换中；文献[12]给出的多元总体最小二乘算法受到正交Procrustes算法制约仅适用于观测值误差对称独立的情况，随机元素相关性问题也没有解决；文献[13]只给出了观测值等精度同方差情况下的多元总体最小二乘算法。对于系数矩阵含有常量的情况，文献[10,14—16]将系数矩阵拆分成分别包含随机元素和常数元素的系数矩阵；文献[11—13]采用消元法，分两步计算随机元素对应的未知参数和常数元素对应的未知参数；文献[8,17]使常数元素的方差和协方差为零，然后进行加权解算；文献[18]首次将EIV模型改写成更具一般性的PEIV(partial EIV)模型，解决了系数矩阵任意位置含有常数元素的平差问题。但以上方法都没有把系数矩阵和观测向量(或矩阵)中的所有元素纳入到一个协因数阵中进行处理。

基于以上问题，本文推导了顾及系数矩阵含有常数元素情况下的多元总体最小二乘的牛顿解法，并讨论了多元总体最小二乘的牛顿法解和拉格朗日乘数法解之间的关系以及给出了多元总体最小二乘的精度评定公式，最后通过模拟数据和实测数据验证了本文算法的可行性和有效性。

1多元总体最小二乘解法

1.1多元变量EIV模型

多元变量EIV模型[6-8]为

Y=(A+EA)Ξ-EY

(1)

(2)

顾及系数矩阵A中某列含有常数元素，在文献[8,17—18]的基础上，系数矩阵和相应的协因数阵定义为

(3)

(4)

式中，As是系数矩阵A的随机元素部分；Ad是系数矩阵A的常数元素部分；Es为随机元素的误差矩阵；Qi(i=1,2,…，t)∈Rn×n(t+m)。

此时，协因数阵Q为奇异阵，但在计算时并没有直接求取Q的逆，因此并不影响解算过程。

1.2多元总体最小二乘问题的牛顿解法

本文在文献[10]推导牛顿法思路的基础上，进一步推导了一般情形下多元总体最小二乘问题的牛顿解法。根据总体最小二乘问题平差准则VTPV=min，由求条件极值的拉格朗日乘数法构造如下目标函数

Γ(V,Ξ,λ)=VTPV+2λT(vec(Y+EY)-

(ΞT⊗In)vec(A+EA))

(5)

式中，λ∈Rnm×1为拉格朗日乘数向量；⊗为Kronecker-Zehfuss积。

对式(5)中的元素分别求导得

(6)

(7)

由式(6)和式(7)，得到改正数最优估值为[10]

(8)

进而得到牛顿法求解无约束问题的目标函数

(9)

式中

式(9)对vec(Ξ)求一阶导数，得

(10)

经推导可得式(10)第1部分

(11)

由逆矩阵求导原理，得式(10)第2部分

(12)

式中

(13)

(14)

由式(11)和式(12)得

(15)

式(9)对vec(Ξ)求二阶导数，得

(16)

并令

(17)

则

(18)

依据牛顿法原理得到未知参数的改正数

(19)

以及

(20)

多元总体最小二乘问题的牛顿解法的迭代步骤为(k为迭代次数)：

作为迭代初始值。

第2步，由上一步获得的数值，分别计算

由式(15)和式(18)，进一步计算Gk+1、Hk+1和

第3步，根据式(20)更新未知参数估值。

作为未知参数的最佳估值。

1.3牛顿法解和拉格朗日乘数法解之间的关系

拉格朗日乘数法的多元总体最小二乘解为[10]

(21)

式(21)在近似值处展开，写成式(20)的形式

(22)

式中，拉格朗日乘数法未知参数改正数为

(23)

令

(24)

利用矩阵反演公式[10,19]，得

(25)

又因为

(26)

联合式(23)、式(25)和式(26)，得

(27)

2多元总体最小二乘解法精度评定

2.1单位权方差估值公式

根据文献[6—7]，多元总体最小二乘解法的单位权方差估值公式为

(28)

式中，m(n-t)为自由度，mn为观测矩阵Y元素的个数，mt为参数矩阵Ξ元素的个数。

2.2协因数阵计算公式

多元总体最小二乘法中，数据多以矩阵形式存在。为计算协因数阵，把矩阵写成向量形式，得到

(29)

表1多元总体最小二乘法协因数阵

Tab.1Cofactor matrice of multivariate total least squares adjustment

L^V^L^ξLQ-^X5Q-^X5-Q^ΣT^XT4^V-^X5^X50Q^ΣT^XT4^LQ-^X50Q-^X50^ξ-^X4^ΣQ^X4^ΣQ0^X3

3算例分析

3.1模拟算例

根据文献[17]中给出的源坐标系的真值数据，本文通过给定的二维仿射变换参数真值计算得到目标坐标系的真值数据。在得到两组坐标系的真值数据后，依据文献[17]，本文在模拟数据中源坐标系所加误差的协因数阵为

Qxy=I2⊗(0.01diag([1,2,3,1,5,4,2,7,2,1,8,3,6]))

模拟数据中目标坐标系所加误差的协因数阵为

QXY=I2⊗(0.01diag([1,3,6,1,1,8,4,3,6,5,4,5,2]))

共模拟500次，分别使用本文提出的牛顿法(N法)和拉格朗日-牛顿法(L-N法)，以及Fang(2011)法[10]、Schaffrin(2009)法[8]和加权最小二乘法(WLS)计算二维仿射变换参数，迭代法的终止条件均为10-12，将获得的二维仿射变换参数的平均值，以及二维仿射变换参数平均值与其真值的差值2范数和平均迭代次数列于表2中。

参数为矩阵时，文献[10]并没有给出在系数矩阵含有常数列时的具体公式算法。依据其提出的在参数为向量时系数矩阵含有常数列的加权总体最小二乘算法，本文编程实现了这种将系数矩阵分成两部分的多元加权总体最小二乘算法，并简称为Fang(2011)法。拉格朗日-牛顿法(L-N法)是指把Fang(2011)法(拉格朗日乘数法)获得的第一次多元加权总体最小二乘参数解作为本文牛顿法迭代的初始值；在选取迭代初始值时，如果使用多元最小二乘参数解作为初始值计算式(20)，则为牛顿法(N法)，如果使用多元加权总体最小二乘参数解作为初始值计算式(20)，则为拉格朗日-牛顿法(L-N法)；在能够获得多元加权总体最小二乘参数解的情况下，建议使用收敛速度更快的L-N法。

以上算例没有考虑系数矩阵元素与观测矩阵元素的相关性。为了验证本文的牛顿法能够解决更一般的多元总体最小二乘问题，在文献[19]的基础上，使用Matlab 7.11.0软件，利用函数randn(100,13*4)T×randn(100,13*4)生成协因数阵Qr，利用函数mvnrnd生成均值均为0，方差-协方差阵为0.001 Qr的随机误差矩阵，并将其施加到真值数据中以作为计算的模拟数据。

由于模拟数据使得两套坐标系具有了相关性，此时Schaffrin(2009)法不再适用，选用本文提出的牛顿法和拉格朗日-牛顿法，以及Fang(2011)法模拟计算10 000次，计算时取Q=Qr，迭代法的终止条件均为10-12，并由获得的10 000次的参数估值绘制出的参数估值频率直方图及其相应的拟合出的正态分布曲线图如图1所示。

计算结果表明，牛顿法的平均迭代次数为6.972 6，拉格朗日-牛顿法的平均迭代次数为6.828 7，Fang(2011)法的平均迭代次数为12.188 2。牛顿法和拉格朗日-牛顿法使用的协因数阵Q为奇异矩阵，但这并不影响未知参数的解算。Fang(2011)法的迭代次数最多，拉格朗日-牛顿法迭代次数要少于牛顿法，这也与上文的结论相符合。参数估值的频率直方图近似于正态分布，相应的正态分布拟合曲线大约也以参数的真值为对称轴。以上的结论证明了多元总体最小二乘牛顿解法的可行性和有效性。

表2　5种方法参数估计结果

图1　所估参数的直方图及其正态分布拟合曲线Fig.1　Histogram and fitting curve of normal distribution of estimated parameters

3.2实测算例

选取文献[19]中二维坐标转换数据作为计算数据，共有5组在源坐标系和目标坐标系下的公共点。

源坐标系和目标坐标系权阵为[19]

Pxy=diag([18.781 7,6.377 4,12.648 9，

17.476 9,22.272 6,23.982 3,13.680 4，

3.465 6,3.732 4,6.437 7])

PXY=diag([9.836 1,5.535 7,12.736 9,

12.009 9,10.181,11.366 1,11.147,

5.883 4,9.832 2,7.567 8])

使用本文提出的牛顿法和拉格朗日-牛顿法，以及Fang(2011)法计算二维仿射变换参数，迭代法的终止条件均为10-12。3种方法的未知参数估值列于表3中。利用表3中的未知参数估值结果和表1中的公式计算未知参数估值协因数阵，并将其列于表4中。

表3　3种方法参数估计结果

表4　参数估值协因数阵

4结论

本文推导了多元总体最小二乘问题的牛顿解法，并比较了基于牛顿法的多元加权总体最小二乘解和基于拉格朗日乘数法的多元加权总体最小二乘解之间的关系，给出了多元总体最小二乘算法精度评定的16种协因数阵的近似计算公式。通过算例验证，并与现有的解法作了比较，得到以下几点结论：

(1) 牛顿法有着更少的迭代次数，更高的解算效率。

(2) 在观测矩阵和系数矩阵元素对应的协因数阵中，本文取常数元素的协因数为0。这虽然导致协因数阵Q奇异，但并不影响参数的估计结果；本文的算法也能解决观测矩阵与系数矩阵元素具有相关性的问题。

(3) 牛顿法的解算效率会受到迭代初始值的影响，本文采用的拉格朗日-牛顿法因其迭代初始值更接近参数最优估值，而使得解算效率比起牛顿法有所提高。

本文给出的多元总体最小二乘解法精度评定的公式只是近似公式，如何推导更为严密的总体最小二乘精度评定公式还需进一步研究。

参考文献：

[1]王乐洋.基于总体最小二乘的大地测量反演理论及应用研究[J]. 测绘学报, 2012, 41(4): 629.

WANG Leyang. Research on Theory and Application of Total Least Squares in Geodetic Inversion[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(4): 629.

[2]GOLUB G H, VAN LOAN C F. An Analysis of the Total Least Squares Problem[J]. SIAM Journal on Numerical Analysis, 1980, 17(6): 883-893.

[3]SPRENT P. Models in Regression and Related Topics[M]. London: Methuen & Co Ltd, 1969.

[4]王乐洋, 许才军. 总体最小二乘研究进展[J]. 武汉大学学报(信息科学版), 2013, 38(7): 850-856, 878.

WANG Leyang, XU Caijun. Progress in Total Least Squares[J]. Geomatics and Information Science of Wuhan University, 2013, 38(7): 850-856, 878.

[5]VAN HUFFEL S, VANDEWALLE J. The Total Least Squares Problem: Computational Aspects and Analysis[D]. Philadelphia, Pennsylvania: SIAM, 1991.

[6]SCHAFFRIN B, FELUS Y A. Multivariate Total Least-squares Adjustment for Empirical Affine Transformations[C]∥VI Hotine-Marussi Symposium on Theoretical and Computational Geodesy: Challenge and Role of Modern Geodesy. Berlin: Springer, 2008, 132: 238-242.

[7]SCHAFFRIN B, FELUS Y A. On the Multivariate Total Least-squares Approach to Empirical Coordinate Transformations. Three Algorithms[J]. Journal of Geodesy, 2008, 82(6): 373-383.

[8]SCHAFFRIN B, WIESER A. Empirical Affine Reference Frame Transformations by Weighted Multivariate TLS Adjustment[C]∥Geodetic Reference Frames. Berlin: Springer, 2009, 134: 213-218.

[9]FELUS Y A, BURTCH R C. On Symmetrical Three-dimensional Datum Conversion[J]. GPS Solutions, 2009, 13(1): 65-74.

[10]FANG X. Weighted Total Least Squares Solutions for Applications in Geodesy[D]. Hanover: Leibniz University of Hanover, 2011.

[11]孙大双, 张友阳, 黄令勇, 等. 多元总体最小二乘在大旋转角三维坐标转换中的应用[J]. 测绘科学技术学报, 2014, 31(5): 481-485.

SUN Dashuang, ZHANG Youyang, HUANG Lingyong, et al. Application of Multivariate Total Least Squares Method in Three-dimension Coordinate Conversion with Large Rotation Angle[J]. Journal of Geomatics Science and Technology, 2014, 31(5): 481-485.

[12]黄令勇, 吕志平, 任雅奇, 等. 多元总体最小二乘在三维坐标转换中的应用[J]. 武汉大学学报(信息科学版), 2014, 39(7): 793-798.

HUANG Lingyong, LÜ Zhiping, REN Yaqi, et al. Application of Multivariate Total Least Square in Three-dimensional Coordinate Transformation[J]. Geomatics and Information Science of Wuhan University, 2014, 39(7): 793-798.

[13]钱承军, 陈义. 多变量总体最小二乘在点云拼接中的应用[J]. 测绘与空间地理信息, 2015, 38(1): 67-69, 76.

QIAN Chengjun, CHEN Yi. Application of Multivariable Total Least Squares in the Registration of Point Clouds[J]. Geomatics & Spatial Information Technology, 2015, 38(1): 67-69, 76.

[14]王乐洋, 许才军, 鲁铁定. 边长变化反演应变参数的总体最小二乘方法[J]. 武汉大学学报(信息科学版), 2010, 35(2): 181-184.

WANG Leyang, XU Caijun, LU Tieding. Inversion of Strain Parameter Using Distance Changes Based on Total Least Squares[J]. Geomatics and Information Science of Wuhan University, 2010, 35(2): 181-184.

[15]XU Caijun, WANG Leyang, WEN Yangmao, et al. Strain Rates in the Sichuan-Yunnan Region Based upon the Total Least Squares Heterogeneous Strain Model from GPS Data

[J]. Terrestrial Atmospheric and Oceanic Sciences, 2011, 22(2): 133-147.

[16]王乐洋, 于冬冬, 吕开云. 复数域总体最小二乘平差[J]. 测绘学报, 2015, 44(8): 866-876. DOI: 10.11947/j.AGCS.2015.20130701.

WANG Leyang, YU Dongdong, LÜ Kaiyun. Complex Total Least Squares Adjustment[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(8): 866-876. DOI: 10.11947/j.AGCS.2015.20130701.

[17]MAHBOUB V. On Weighted Total Least-squares for Geodetic Transformations[J]. Journal of Geodesy, 2012, 86(5): 359-367.

[18]XU Peiliang,LIU Jingnan,SHI Chuang.Total Least Squares Adjustment in Partial Errors-in-variables Models: Algorithm and Statistical Analysis[J]. Journal of Geodesy, 2012, 86(8): 661-675.

[19]JAZAERI S, AMIRI-SIMKOOEI A R, SHARIFI M A. Iterative Algorithm for Weighted Total Least Squares Adjustment[J]. Survey Review, 2014, 46(334): 19-27.

[20]BOYD S, VANDENBERGHE L. Convex Optimization[M]. Cambridge: Cambridge University Press, 2004.

[21]王乐洋, 鲁铁定. 总体最小二乘平差法的误差传播定律[J]. 大地测量与地球动力学, 2014, 34(2): 55-59.

WANG Leyang, LU Tieding. Propagation Law of Errors in Total Least Squares Adjustment[J]. Journal of Geodesy and Geodynamics, 2014, 34(2): 55-59.

[22]ZHOU Yongjun, KOU Xinjian, ZHU Jianjun, et al. A Newton Algorithm for Weighted Total Least-squares Solution to a Specific Errors-in-variables Model with Correlated Measurements[J]. Studia Geophysica et Geodaetica, 2014, 58(3): 349-375.

[23]王乐洋, 许才军, 鲁铁定. 病态加权总体最小二乘平差的岭估计解法[J]. 武汉大学学报(信息科学版), 2010, 35(11): 1346-1350.

WANG Leyang, XU Caijun, LU Tieding. Ridge Estimation Method in Ill-posed Weighted Total Least Squares Adjustment[J]. Geomatics and Information Science of Wuhan University, 2010, 35(11): 1346-1350.

(责任编辑：丛树平)

修回日期： 2015-10-12

First author： WANG Leyang (1983—), male, PhD, associate professor, majors in geodetic inversion and geodetic data processing.

E-mail： wleyang@163.com

A Newton Algorithm for Multivariate Total Least Squares Problems

WANG Leyang1,2,3,ZHAO Yingwen1,CHEN Xiaoyong1,2,3,ZANG Deyan1,2,3

1. Faculty of Geomatics, East China University of Technology, Nanchang 330013, China; 2. Key Laboratory of Watershed Ecology and Geographical Environment Monitoring, NASG，Nanchang 330013, China; 3. Jiangxi Province Key Lab for Digital Land, Nanchang 330013, China

Abstract：In order to improve calculation efficiency of parameter estimation, an algorithm for multivariate weighted total least squares adjustment based on Newton method is derived. The relationship between the solution of this algorithm and that of multivariate weighted total least squares adjustment based on Lagrange multipliers method is analyzed. According to propagation of cofactor, 16 computational formulae of cofactor matrices of multivariate total least squares adjustment are also listed. The new algorithm could solve adjustment problems containing correlation between observation matrix and coefficient matrix. And it can also deal with their stochastic elements and deterministic elements with only one cofactor matrix. The results illustrate that the Newton algorithm for multivariate total least squares problems could be practiced and have higher convergence rate.

Key words：multivariate total least squares; Newton method; propagation of cofactor; affine

第一作者简介：王乐洋(1983—)，男，博士，副教授，研究方向为大地测量反演及大地测量数据处理。

收稿日期：2015-05-11

基金项目：国家自然科学基金(41204003;41161069;41304020;41464001)；测绘地理信息公益性行业科研专项(201512026)；江西省自然科学基金(20151BAB203042)；江西省教育厅科技项目(GJJ150595;KJLD12077;KJLD14049)；流域生态与地理环境监测国家测绘地理信息局重点实验室开放基金(WE2015005)；东华理工大学博士科研启动金(DHBK201113)；江西省研究生创新专项资金(YC2015-S266;YC2015-S267)；东华理工大学研究生创新专项资金(DHYC-2015005)；对地观测技术国家测绘地理信息局重点实验室开放基金(K201502)

中图分类号：P207

文献标识码：A

文章编号：1001-1595(2016)04-0411-07

引文格式：王乐洋，赵英文，陈晓勇,等.多元总体最小二乘问题的牛顿解法[J].测绘学报，2016,45(4)：411-417,424. DOI:10.11947/j.AGCS.2016.20150246.

WANG Leyang,ZHAO Yingwen,CHEN Xiaoyong,et al.A Newton Algorithm for Multivariate Total Least Squares Problems[J]. Acta Geodaetica et Cartographica Sinica,2016,45(4):411-417,424. DOI:10.11947/j.AGCS.2016.20150246.