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

多元总体最小二乘问题的牛顿解法王乐洋;赵英文;陈晓勇;臧德彦【摘要】为提高多元总体最小二乘问题参数估值的解算效率，推导了基于牛顿法的多元加权总体最小二乘算法;分析比较了基于牛顿法的多元加权总体最小二乘解和基于拉格朗日乘数法多元加权总体最小二乘解之间的关系，根据协因数传播律给出了多元总体最小二乘平差的16种协因数阵的近似计算公式.新算法能够解决观测矩阵和系数矩阵元素具有相关性的问题，并且可以把观测矩阵和系数矩阵的随机元素和常数元素纳入到一个协因数阵中进行处理.算例结果表明，本文提出的多元总体最小二乘问题的牛顿解法可行且收敛速度更快.【期刊名称】《测绘学报》【年(卷)，期】2016(045)004【总页数】8页(P411-417,424)【关键词】多元总体最小二乘;牛顿法;协因数传播律;仿射变换【作者】王乐洋;赵英文;陈晓勇;臧德彦【作者单位】东华理工大学测绘工程学院，江西南昌330013;流域生态与地理环境监测国家测绘地理信息局重点实验室，江西南昌330013;江西省数字国土重点实验室江西南昌330013;东华理工大学测绘工程学院，江西南昌330013;东华理工大学测绘工程学院，江西南昌330013;流域生态与地理环境监测国家测绘地理信息局重点实验室，江西南昌330013;江西省数字国土重点实验室，江西南昌330013;东华理工大学测绘工程学院，江西南昌330013;流域生态与地理环境监测国家测绘地理信息局重点实验室，江西南昌330013;江西省数字国土重点实验室，江西南昌330013【正文语种】中文【中图分类】P207transformationFoundationsupport:TheNationalNaturalScienceFoundationofChina(Nos.41204003;41161069;41304020;41464001);TheNationalDepartmentPublicBenefitResearchFoundation(Surveying，MappingandGeoinformation)(No.201512026);TheNaturalScienceFoundationofJiangxiProvince(No.20151BAB203042);TheScienceandTechnologyProjectoftheEducationDepartmentofJiangxiProvince(Nos.GJJ150595;KJLD12077;KJLD14049);TheFoundofKeyLaboratoryofWatershedEcologyandGeographicalEnvironmentMonitoring(No.WE2015005);TheScientificResearchFoundationofECIT(No.DHBK201113);InnovationFundDesignatedforGraduateStudentsofJiangxiProvince(Nos.YC2015-S266;YC2015-S267);InnovationFundDesignatedforGraduateStudentsofECIT(No.DHYC2015005);TheFoundoftheKeyLaboratoryofMappingfromSpace,NASG(No.K201502)总体最小二乘(totalleastsquares，TLS)方法是近30多年来发展起来的一种能同时顾及观测值误差和模型系数矩阵误差的数学方法，其理论及应用研究是目前国内外研究的热点问题［1］。文献［2］首次提出总体最小二乘概念，随着变量误差模型(errors-in-variables，EIV)由单变量模型扩展到多元变量模型［3］，总体最小二乘也相应扩展到多元总体最小二乘(multivariateTLS，MTLS)［4-5］。某些模型的系数矩阵会因未知参数写成向量形式含有重复出现的随机元素和大量的零元素；多元总体最小二乘将参数向量和观测向量改写成矩阵形式，避免了这种现象，进而削弱了系数矩阵中行列之间的相关性，而且使得系数矩阵的构造更加简单。在国夕卜相关研究中，文献［5—7］认为观测矩阵和系数矩阵中的随机元素是等精度的，给出了基于奇异值分解法和拉格朗日乘数法的多元总体最小二乘算法；文献［8］讨论了基于拉格朗日乘数法的多元加权总体最小二乘算法，将随机元素的方差阵分解成两个矩阵取直积的形式，但是当观测数据含有两个以上的方差分量或是观测数据含有多类观测值时，随机变量的方差协方差阵就很难写成或是不能写成两个矩阵取直积的形式；文献［9］推导了基于奇异值分解法的多元加权总体最小二乘算法，但是这种算法只使用了行尺度权阵，其权阵不具有一般性。以上文献的所有算法都不适用于系数矩阵与观测矩阵中随机元素具有相关性的情况。文献［10］给出了基于拉格朗日乘数法的多元加权总体最小二乘算法，其权阵具有一般性，但没有具体给出系数矩阵含有常数元素的解法以及没有解决算法的精度评定问题，这种算法的解算效率也还需要验证。在国内相关研究中，文献［11］将基于奇异值分解法的多元总体最小二乘算法应用在三维坐标转换中；文献［12］给出的多元总体最小二乘算法受到正交Procrustes算法制约仅适用于观测值误差对称独立的情况，随机元素相关性问题也没有解决；文献［13］只给出了观测值等精度同方差情况下的多元总体最小二乘算法。对于系数矩阵含有常量的情况，文献［10,14—16］将系数矩阵拆分成分别包含随机元素和常数元素的系数矩阵；文献［11—13］采用消元法，分两步计算随机元素对应的未知参数和常数元素对应的未知参数；文献［8,17］使常数元素的方差和协方差为零，然后进行加权解算；文献［18］首次将EIV模型改写成更具一般性的PEIV(partialEIV)模型，解决了系数矩阵任意位置含有常数元素的平差问题。但以上方法都没有把系数矩阵和观测向量(或矩阵)中的所有元素纳入到一个协因数阵中进行处理。基于以上问题，本文推导了顾及系数矩阵含有常数元素情况下的多元总体最小二乘的牛顿解法，并讨论了多元总体最小二乘的牛顿法解和拉格朗日乘数法解之间的关系以及给出了多元总体最小二乘的精度评定公式，最后通过模拟数据和实测数据验证了本文算法的可行性和有效性。1.1多元变量EIV模型多元变量EIV模型［6-8］为式中，YeRnxm为观测矩阵；EYuRnxm为Y的误差矩阵；AuRnxt为系数矩阵；EAuRnxt为A的误差矩阵；MuRtxm为参数矩阵；V为系数矩阵和观测矩阵中随机元素的误差向量，这与文献［10］的表达一致；为未知单位权方差；QAA、QYY、QAY和QYA为系数矩阵随机元素和观测矩阵随机元素的协因数阵以及它们之间的互协因数阵；P为系数矩阵随机元素和观测矩阵随机元素的权阵；vec(・)为矩阵按列拉直运算。顾及系数矩阵A中某列含有常数元素，在文献［8,17—18］的基础上，系数矩阵和相应的协因数阵定义为式中，As是系数矩阵A的随机元素部分；Ad是系数矩阵A的常数元素部分；Es为随机元素的误差矩阵；Qi(i=12...,t)uRnxn(t+m)。此时，协因数阵Q为奇异阵，但在计算时并没有直接求取Q的逆，因此并不影响解算过程。1.2多元总体最小二乘问题的牛顿解法本文在文献［10］推导牛顿法思路的基础上，进一步推导了一般情形下多元总体最小二乘问题的牛顿解法。根据总体最小二乘问题平差准则VTPV=min，由求条件极值的拉格朗日乘数法构造如下目标函数「(V,三，入)=VTPV+2入T(vec(Y+EY)-(三TIn)vec(A+EA))式中，入uRnmxl为拉格朗日乘数向量；区为Kronecker-Zehfuss积。对式(5)中的元素分别求导得施=0式中，矩阵，这与文献[10]的表达一致，符号〃人”表示最优估值。由式(6)和式(7)，得到改正数最优估值为[10]进而得到牛顿法求解无约束问题的目标函数式中式(9)对vec(三)求一阶导数，得经推导可得式(10)第1部分由逆矩阵求导原理，得式(10)第2部分式中式中，ei=[0...010...0]T为m维列向量，除了第i个元素为1夕卜，其余的元素全部为0;i=[(j-1)/t]+1;[・]为取整运算；。由式(11)和式(12)得式中，区(A+EA)。式(9)对vec(三)求二阶导数，得并令式中，i,j=12...,tm;;MuRtmxtm;sqrt(・)为开方运算。则依据牛顿法原理得到未知参数的改正数以及式中，)为未知参数最佳估值；vec(三0)为未知参数近似值。多元总体最小二乘问题的牛顿解法的迭代步骤为(k为迭代次数):第1步，根据给定的A、Y、Q,利用式(14)和式(17)构造和，计算最小二乘解((Im®区A))-L(Im®作为迭代初始值。第2步，由上一步获得的数值，分别计算由式(15)和式(18)，进一步计算Gk+1、Hk+1和第3步，根据式(20)更新未知参数估值。第4步，重复第2和3步，直到当＜£何为很小正值，本文取10-12)时，并通过逆拉直运算计算作为未知参数的最佳估值。1.3牛顿法解和拉格朗日乘数法解之间的关系拉格朗日乘数法的多元总体最小二乘解为［10］(Im®式(21)在近似值处展开，写成式(20)的形式(Im®vec(H0)+vec(A)L式中，拉格朗日乘数法未知参数改正数为令利用矩阵反演公式［10,19］，得又因为联合式(23)、式(25)和式(26)，得(TX1)=vec(A)L-Rvec(A)L式中，ItmuRtmxtm为单位矩阵；。由式(27)可知，牛顿法未知参数改正数与拉格朗日乘数法未知参数改正数的区别在于R矩阵。拉格朗日乘数法把平差准则的条件极值转化成式(5)的自由极值，对目标函数式(5)中的变量分别求一阶偏导数，并令其为零来求取目标函数的极小值；牛顿法判断极小值不仅需要目标函数对变量的一阶导数为零，还需要二阶导数大于零，即Hessian矩阵正定，这使得牛顿法对全局的判断能力更强，寻找极小值的速度也更快［20］；是一阶导数意义下的未知参数改正数，)N是二阶导数意义下的未知参数改正数，目标函数的二阶导数部分组成R矩阵。一般情形下，当目标函数的二阶导数存在时，R矩阵也是存在的，)L可以理解为拉格朗日乘数法未知参数改正数的修正量，这也正是牛顿法比拉格朗日乘数法迭代效率高的原因。当R或)L为零时，牛顿法就等价于拉格朗日乘数法。2.1单位权方差估值公式根据文献［6—7］，多元总体最小二乘解法的单位权方差估值公式为式中，m(n-t)为自由度，mn为观测矩阵Y元素的个数，mt为参数矩阵三元素的个数。当系数矩阵含有常数元素时，-1；式(28)也与文献［6］的式(3.7)和文献［7］的式(37)相一致，但本式的表达具有一般性。2.2协因数阵计算公式多元总体最小二乘法中，数据多以矩阵形式存在。为计算协因数阵，把矩阵写成向量形式，得到式(29)将向量和近似写成L的函数，根据协因数传播律得到各个向量间的一阶近似协因数阵［21］。本文求出的多元总体最小二乘问题的未知参数的协因数阵也与文献［19,22—23］中给出的一般总体最小二乘问题的未知参数的协因数阵在表达形式上相一致。所有协因数阵如表1所示，其中，；；Q。3.1模拟算例根据文献［17］中给出的源坐标系的真值数据，本文通过给定的二维仿射变换参数真值计算得到目标坐标系的真值数据。在得到两组坐标系的真值数据后，依据文献［17］，本文在模拟数据中源坐标系所加误差的协因数阵为Qxy=I2®(0.01diag(［1,2,3,L5,4,2,7,2,L&3,6］))模拟数据中目标坐标系所加误差的协因数阵为QXY=I2®(0.01diag(［136,1,1,843,65452］))共模拟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法。由表2可知，N法、L-N法、Fang(2011)法和Schaffrin(2009)法这4种多元加权总体最小二乘法获得的参数结果一致，并且参数估值的差值范数要小于WLS法。从平均迭代次数来看，Schaffrin(2009)法的迭代次数最多，Fang(2011)法次之，L-N法的迭代次数要稍小于N法，为最小。本文提出的多元总体最小二乘牛顿解法的迭代效率要高于Fang(2011)法和Schaffrin(2009)法这两种解法。由文献［20］可知，因为牛顿法是一种用H(Hessian)矩阵导出的最速下降法，牛顿步径也即)N是好的搜索方向，而且牛顿法也具有二次收敛速度，所以，牛顿法的迭代次数最少。本文提出的牛顿法比Fang(2011)法(拉格朗日乘数法)有着更高的收敛速度，它们的区别主要是由)L矩阵决定的，而且当三0靠近时，牛顿法有着更高的收敛速度，这也解释了为何L-N法比N法有着更少的迭代次数。以上算例没有考虑系数矩阵元素与观测矩阵元素的相关性。为了验证本文的牛顿法能够解决更一般的多元总体最小二乘问题，在文献［19］的基础上，使用Matlab7.11.0软件，利用函数randn(100,13*4)Txrandn(100,13*4)生成协因数阵Qr,利用函数mvnrnd生成均值均为0，方差-协方差阵为0.001Qr的随机误差矩阵，并将其施加到真值数据中以作为计算的模拟数据。由于模拟数据使得两套坐标系具有了相关性，此时Schaffrin(2009)法不再适用，选用本文提出的牛顿法和拉格朗日-牛顿法，以及Fang(2011)法模拟计算10000次，计算时取Q=Qr,迭代法的终止条件均为10-12，并由获得的10000次的参数估值绘制出的参数估值频率直方图及其相应的拟合出的正态分布曲线图如图1所示。计算结果表明，牛顿法的平均迭代次数为6.9726，拉格朗日-牛顿法的平均迭代次数为6.8287,Fang(2011)法的平均迭代次数为12.1882。牛顿法和拉格朗日-牛顿法使用的协因数阵Q为奇异矩阵，但这并不影响未知参数的解算。Fang(2011)法的迭代次数最多，拉格朗日-牛顿法迭代次数要少于牛顿法，这也与上文的结论相符合。参数估值的频率直方图近似于正态分布，相应的正态分布拟合曲线大约也以参数的真值为对称轴。以上的结论证明了多元总体最小二乘牛顿解法的可行性和有效性。3.2实测算例选取文献[19]中二维坐标转换数据作为计算数据，共有5组在源坐标系和目标坐标系下的公共点。源坐标系和目标坐标系权阵为[19]Pxy=diag([18.7817,6.3774,12.6489，17.4769,22.2726,23.9823,13.6804，3.4656,3.7324,6.4377])PXY=diag([9.8361,5.5357,12.7369,12.0099,10.181,11.3661,11.147,5.8834,9.8322,7.5678])使用本文提出的牛顿法和拉格朗日-牛顿法，以及Fang(2011)法计算二维仿射变换参数，迭代法的终止条件均为10-12。3种方法的未知参数估值列于表3中。利用表3中的未知参数估值结果和表1中的公式计算未知参数估值协因数阵，并将其列于表4中。表3为3种方法的共同结果，其精度一致，且为0.147294126638。结果表明，牛顿法迭代了2次，拉格朗日-牛顿法迭代了1次，Fang(2011)法迭代了3次，这也说明了本文算法的合理性。本文推导了多元总体最小二乘问题的牛顿解法，并比较了基于牛顿法的多元加权总体最小二乘解和基于拉格朗日乘数法的多元加权总体最小二乘解之间的关系，给出了多元总体最小二乘算法精度评定的16种协因数阵的近似计算公式。通过算例验证，并与现有的解法作了比较，得到以下几点结论：牛顿法有着更少的迭代次数，更高的解算效率。在观测矩阵和系数矩阵元素对应的协因数阵中，本文取常数元素的协因数为0。这虽然导致协因数阵Q奇异，但并不影响参数的估计结果；本文的算法也能解决观测矩阵与系数矩阵元素具有相关性的问题。⑶牛顿法的解算效率会受到迭代初始值的影响，本文采用的拉格朗日-牛顿法因其迭代初始值更接近参数最优估值，而使得解算效率比起牛顿法有所提高。本文给出的多元总体最小二乘解法精度评定的公式只是近似公式，如何推导更为严密的总体最小二乘精度评定公式还需进一步研究。修回日期：2015-10-12Firstauthor:WANGLeyang(1983—),male,PhD,associateprofessor,majorsingeodeticinversionandgeodeticdataprocessing.E-mail:【相关文献】王乐洋.基于总体最小二乘的大地测量反演理论及应用研究[J].测绘学报,2012,41(4):629.WANGLeyang.ResearchonTheoryandApplicationofTotalLeastSquaresinGeodeticInversion[J].ActaGeodaeticaetCartographicaSinica,2012,41(4):629.[2]GOLUBGH,VANLOANCF.AnAnalysisoftheTotalLeastSquaresProblem[J].SIAMJournalonNumericalAnalysis,1980,17(6):883-893.[3]SPRENTP.ModelsinRegressionandRelatedTopics[M].London:Methuen&CoLtd,1969.王乐洋，许才军.总体最小二乘研究进展[J].武汉大学学报(信息科学版),2013,38(7):850856,878.WANGLeyang,XUCaijun.ProgressinTotalLeastSquares[J].GeomaticsandInformationScienceofWuhanUniversity,2013,38(7):850-856,878.[5]VANHUFFELS,VANDEWALLEJ.TheTotalLeastSquaresProblem:ComputationalAspectsandAnalysis[D].Philadelphia,Pennsylvania:SIAM,1991.SCHAFFRINB,FELUSYA.MultivariateTotalLeastsquaresAdjustmentforEmpiricalAffineTransformations[C]llVIHotine-MarussiSymposiumonTheoreticalandComputationalGeodesy:ChallengeandRoleofModernGeodesy.Berlin:Springer,2008,132:238-242.SCHAFFRINB,FELUSYA.OntheMultivariateTotalLeastsquaresApproachtoEmpiricalCoordinateTransformations.ThreeAlgorithms[J].JournalofGeodesy,2008,82(6):373-383.[8]SCHAFFRINB,WIESERA.EmpiricalAffineReferenceFrameTransformationsbyWeightedMultivariateTLSAdjustment[C]llGeodeticReferenceFrames.Berlin:Springer,2009,134:213-218.FELUSYA,BURTCHRC.OnSymmetricalThreedimensionalDatumConversion[J].GPSSolutions,2009,13(1):65-74.[10]FANGX.WeightedTotalLeastSquaresSolutionsforApplicationsinGeodesy[D].Hanover:LeibnizUniversityofHanover,2011.孙大双，张友阳,黄令勇,等.多元总体最小二乘在大旋转角三维坐标转换中的应用[J].测绘科学技术学报,2014,31(5):481-485.SUNDashuang,ZHANGYouyang,HUANGLingyong,etal.ApplicationofMultivariateTotalLeastSquaresMethodinThree-dimensionCoordinateConversionwithLargeRotationAngle[J].JournalofGeomaticsScienceandTechnology,2014,31(5):481-485.黄令勇，吕志平，任雅奇,等.多元总体最小二乘在三维坐标转换中的应用[J].武汉大学学报(信息科学版),2014,39(7):793-798.HUANGLingyong,LUZhiping,RENYaqi,etal.ApplicationofMultivariateTotalLeastSquareinThreedimensionalCoordinateTransformation[J].GeomaticsandInformationScienceofWuhanUniversity,2014,39(7):793-798.钱承军，陈义.多变量总体最小二乘在点云拼接中的应用[J].测绘与空间地理信息,2015,38(1):67-69,76.QIANChengjun,CHENYi.ApplicationofMultivariableTotalLeastSquaresintheRegistrationofPointClouds[J].Geomatics&SpatialInformationTechnology,2015,38(1):67-69,76.王乐洋，许才军,鲁铁定.边长变化反演应变参数的总体最小二乘方法[J].武汉大学学报(信息科学版),2010,35(2):181-184.WANGLeyang,XUCaijun,LUTieding.InversionofStrainParameterUsingDistanceChangesBasedonTotalLeastSquares[J].GeomaticsandInformationScienceofWuhanUniversity,2010,35(2):181-184.XUCaijun,WANGLeyang,WENYangmao,etal.StrainRatesintheSichuan-YunnanRegionBasedupontheTotalLeastSquaresHeterogeneousStrainModelfromGPSData[J].TerrestrialAtmosphericandOceanicSciences,2011,22(2):133-147.王乐洋,于冬冬,吕开云.复数域总体最小二乘平差[J].测绘学报,2015,44(8):866876.DOI:10.11947/j.AGCS.2015.20130701.WANGLeyang,YUDongdong,LUKaiyun.ComplexTotalLeastSquaresAdjustment[J].ActaGeodaeticaetCartographicaSinica,2015,44(8):866876.DOI:10.11947/j.AGCS.2