附有不等式约束的加权整体最小二乘算法

曾文宪,方 兴,刘经南,2,姚宜斌

1.武汉大学测绘学院,湖北武汉 430079;2.武汉大学卫星导航定位技术研究中心,湖北武汉 430079

附有不等式约束的加权整体最小二乘算法

曾文宪1,方 兴1,刘经南1,2,姚宜斌1

1.武汉大学测绘学院,湖北武汉 430079;2.武汉大学卫星导航定位技术研究中心,湖北武汉 430079

针对现有附有不等式约束的整体最小二乘算法的缺陷,以partial EIV(errors-in-variables)模型为基础,在整体最小二乘准则下,通过将附有不等式约束的EIV模型的求解转换为标准的附有不等式约束的最优化问题,并采用惩罚函数法等方法得到了附有不等式约束的加权整体最小二乘新算法。新算法将现有算法的特殊权阵限制条件扩展到了一般性权矩阵,将要求系数矩阵元素全部随机的限定条件扩展到了可同时包含随机和非随机元素的一般情况,并且新算法解决了现有算法计算量受制于约束方程数量的缺陷。实例计算表明,本文提出的算法简单、有效,具有普遍适用性。

整体最小二乘估计;EIV模型;不等式约束;非线性算法

1 引 言

整体最小二乘估计(total least squares, TLS)作为EIV[1](errors-in-variables)模型的严密估计方法,目前已广泛应用于大地测量等众多科学研究和工程应用领域。文献[2]提出了整体最小二乘准则。文献[3]基于正交回归原理推导了TLS数值算法。假定观测值不相关情况下,文献[4]首次提出了真正统计意义上的TLS算法。文献[5]提出了著名的奇异值分解(singular value decomposition,SVD)算法。文献[6]研究了稳健整体最小二乘算法(robust TLS)等。文献[7]证明了当观测数趋于无穷大时,TLS估计具有渐进无偏性。大地测量领域针对普遍存在的观测值不等精度、相关的情况,文献[8—11]提出了加权TLS算法(weighted TLS,WTLS),其中,文献[11]研究了最一般性权矩阵条件下的WTLS算法。文献[1]提出了基于partial EIV模型(PEIV)的WTLS算法,能够将结构性等各种形式的系数矩阵纳入统一的模型形式求解。其他TLS算法的发展情况见文献[12]。

当参数估计存在先验信息时,EIV模型扩展为附有约束的EIV模型。相当多的文献研究了附有等式约束的EIV模型(equality-constrained EIV,ECEIV),如文献[13]、文献[14]以及文献[15]分别提出了3种不同的附有等式约束的TLS算法(equality-constrained TLS,ECTLS)。

某些先验信息要求对模型附加不等式约束条件,就形成了附有不等式约束的EIV模型(inequality-constrained EIV,ICEIV),如变形监测中某因素引起的变形系数理论上要大于某一数值,应对该估计值进行最小值约束;坐标转换模型的TLS算法中要求强制附合到某中心点,即该点坐标的改正数不能超出一定数值范围等。关于附有不等式约束的EIV模型目前仅检索到两篇研究文献,文献[16]运用拉格朗日乘数法将ICTLS问题转化为广义线性互补问题(linear complementarity problem,LCP),再通过排列组合的方法进行求解。文献[17]提出了基于穷举法的附有不等式约束的整体最小二乘算法(inequality-constrained TLS,ICTLS)。上述文献开启了ICTLS算法研究的先河,但存在以下问题有待解决:①算法的计算量随约束方程个数的增长呈指数增长,如当约束方程个数为50时,组合或搜索次数达到约250≈1015,因此,约束方程较多时,算法的计算量急剧增长甚至导致无法计算;②算法仅适用于系数矩阵元素全部为随机元素的EIV模型以及等精度、不相关观测值。

本文以partial EIV模型为基础,研究了附有不等式约束的加权整体最小二乘(inequality-constrained weighted TLS,ICWTLS)算法。本文提出的ICWTLS算法适用于随机和非随机元素并存、或呈现结构性特征的一般性系数矩阵以及不等精度、相关观测值,并且算法的计算效率远高于基于穷举法的ICTLS算法,计算量不受约束方程个数的制约。

2 附有不等式约束的partial-EIV模型

ICEIV模型的数学表达式为[7,17])

式中,y表示n×1的观测向量;ey表示y的n×1随机误差向量,且ey期望为0、协因数阵Qy＝w－1σ2(w和σ2分别表示观测向量的权阵以及单位权方差);β表示t×1的参数向量;A表示n×t的系数矩阵;EA表示A的n×t随机误差矩阵, EA的期望为0且协因数阵为QA＝ω－1σ2(ω表示系数矩阵观测元素的权阵),通常假定EA与ey不相关;G表示不等式约束方程的k×t系数矩阵;z表示k×1的常数向量。

现有ICTLS算法[16-17]假定式(1)中观测数据的权阵ω和w均为单位阵,并且A中元素全部为随机量。当A中存在非随机的固定元素或者呈现结构性特征时,必须对系数矩阵的协因数阵[18]或者对系数矩阵[15]进行特殊处理后求解。为了得到一般性系数矩阵下的统一算法,以partial EIV模型[1]为基础,将传统的ICEIV模型(1)改写为如下形式

式(2)构成了附有不等式约束的partial EIV模型(inequality-constvained partial EIV,ICPEIV),与传统的ICEIV模型[16-17]比较,其优势主要体现为:①将ICEIV模型中系数矩阵元素全部为随机量的限定扩展到了同时包含随机和非随机元素的一般情况;②ICPEIV模型可解析结构性系数矩阵,保证了算法的统一性;③ICPEIV模型A中参与平差计算的待估量个数m小于等于ICEIV中A中元素个数nt,尤其当系数矩阵中的随机量个数较少时,ICPEIV模型形式可以大大减少待估计量。因此,ICPEIV模型更具一般性,以下在不限定观测数据权阵的等精度和相关性的一般情况下,即ICPEIV模型(2)中权矩阵ω和w为任意正定对称阵,笔者提出了普遍适用的附有不等式约束的加权整体最小二乘新算法。

3 (ICWTLS)算法

可以求出其最优估计值。由式(2)可得

将ey和ea代入式(3),则ICPEIV模型(2)的ICWTLS算法转化为以下附有不等式约束的最优化问题

通过将ICPEIV模型的估计转换为附有不等式约束的最优化问题,避免了现有算法采用穷举法引起的算法受限于不等式方程个数的缺陷。根据最优化理论,提出了基于罚函数法[19]的ICWTLS新算法。

式中,μ表示惩罚参数(μ＞0);ci(β)＝max[0,－Giβ＋zi](i＝1,2,…,k);Gi为G的第i行;zi为z的第i个元素。式(5)存在一阶导数且一阶导数连续。

基于罚函数法的ICWTLS算法计算过程如下。

(3)如果参数解满足收敛条件,结束计算,转至步骤(4)。

(4)按给定规则增大μ0,转至步骤(2)进行迭代计算。

对于ICPEIV模型估计转换后的最优化式(4),除上述基于罚函数的ICWTLS算法外,同样可以采用最优估计理论的有效集法、序列二次规划法、内点算法等[19-20]得到相应的ICWTLS算法。基于上述不同最优估计方法的ICWTLS算法的特点是下一步要讨论的问题。

4 实例分析

为了说明本文提出的附有不等式约束的整体最小二乘算法的应用,笔者共模拟两个实例进行了计算。实例1的主要目的是验证和比较ICWTLS新算法与现有基于穷举法的ICTLS算法[17]的优缺点及计算效率,因此,笔者设计了一组观测数据等精度并且系数矩阵全部为随机元素的ICEIV模型数据。实例2选择了附有不等式约束的平面拟合模型,该模型的系数矩阵同时存在随机和非随机元素,现有算法无法解算,通过模拟一组不等精度的试验数据,采用ICWTLS算法求出了模型的参数解。

4.1 实例1(等精度、系数矩阵为随机元素)

实例1的模拟数据见表1,模型包含7×1待估参数向量β以及70×1待估系数向量¯a,系数矩阵A的全部元素为随机量,A和y中所有观测元素不相关且中误差均为0.1,其中参数真值β～见表2。为了改进模型的估计结果,利用参数的先验信息设计了18个不等式约束方程(见表1),其中,0≤βi≤0.8(i＝1,2,3,4)以及－0.5≤βj≤0 (j＝5,6,7)可分别表示为

式中,I4和I3分别表示4阶和3阶的单位阵。

表1 附有不等式约束的平差模型数据Tab.1 Data set of the inequality-constrained adjustment model

(2)从算法的计算量比较,本文算法只需迭代计算23次,而穷举法理论上排列组合所需计算次数为218≈260 000。当模型约束方程个数较多的情况下,基于穷举法的ICTLS算法计算量要远大于本文提出的ICWTLS算法。当不等式约束方程过多时,基于穷举法的算法甚至无法在有效时间内求解。

表2 ICWTLS算法估计结果Tab.2 Estimation of ICWTLS

4.2 实例2(附有不等式约束的平面拟合模型)

线性回归模型是测绘领域常用的基本EIV模型之一,实例2选用平面拟合模型说明ICWTLS算法的应用,平面拟合模型式(1)中观测向量、系数矩阵、参数向量等形式如下

由上式可以看到平面拟合模型的系数矩阵同时包含了非随机和随机元素,即第1列元素为已知量,第2列和第3列是由观测数据构成的随机量,现有附有不等式约束的整体最小二乘算法无法处理这类情况。本文算法可以解算任意结构性系数矩阵形式的EIV模型,将上式表示为partial EIV模型式(2)的形式,式中的矩阵和向量可表示为

式中,1n×1表示元素均为1的n×1单位向量; 02n×1表示元素均为0的2n×1向量;0n×n表示n维的零方阵;I2n×2n表示2n维的单位阵;bn×1＝ [b1b2…bn]T;cn×1＝[c1c2…cn]T;ebn×1＝[eb1eb2…ebn]T;ecn×1＝[ec1ec2…ecn]T。

假定根据先验信息,要求模型的截距参数β1和斜率参数β2必须在如下数值范围内

式(7)相当于对平面拟合模型附加了4个参数β的不等式约束方程,相应的G和z见表3。笔者共模拟了10个拟合点数据,系数矩阵随机量和观测向量的中误差见对应观测数据括号内数值。估计结果见表4,可以看到,ICWTLS算法参数结果均满足所有不等式约束条件(式(7)),而TLS解并不满足不等式约束方程。即当平差模型存在可靠的先验信息时,将其作为附加约束条件,可得到满足设定条件的估计结果。

表3 附有不等式约束的平面拟合模型数据Tab.3 Data set of the inequality-constrained plane fitting model

表4 附有不等式约束的平面拟合模型的ICWTLS估计结果Tab.4 Estimation of inequality-constrained WTLS and TLS of the adjustment model

5 结 论

当EIV模型存在参数的先验信息时,如模型某些参数取值应在一定范围内或者地形拟合的边界条件等就构成了EIV模型的不等式约束条件。若平差计算顾及到约束条件,可以充分利用模型的先验信息改善估计结果或者使得估计结果满足设定的条件。目前仅有文献[16—17]对ICTLS进行了研究,但提出的算法只适用于系数矩阵全部元素随机、观测数据等权的特殊情况。此外,当约束方程数量较大时,建立在穷举法基础上的ICTLS算法的计算量随约束方程个数呈指数增长;当约束方程数量过大时,算法甚至可能无法在有效时间内进行计算。本文将ICEIV模型改写为更为一般化的ICPEIV模型,并且在整体最小二乘准则下,将模型的求解转化为标准的附有不等式约束的最优化问题,提出的ICWTLS新算法能采用统一的方法估计结构性系数矩阵、随机和非随机元素并存的各种情况,并且没有限制观测数据的精度和相关性。同时,算法的计算量不受约束方程个数的制约。论文通过两个实例对ICWTLS新算法进行了验证和说明,计算结果表明算法较好地解决了现有算法的限制问题,新算法在实际应用中简单、有效、具有普遍适用性。

[1] XU P L,LIU J N,SHI C.Total Least Squares Adjustment inPartial Errors-in-variables Models:Algorithm and Statistical Analysis[J].Journal of Geodesy,2012,86(8):661-675.

[2] ADCOCK R J.Note on the Method of Least Squares[J].The Analyst,1877,4(6):183-184.

[3] PEARSON K.On Lines and Planes of Closest Fit to Systems of Points in Space[J].Philosophical Magazine Series 6, 1901,2(11):559-572.

[4] DEMING W E.The Application of Least Squares[J].Philosophical Magazine Series 7,1931,11(68):146-158.

[5] 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.

[6] WASTON G A.Robust Counterparts of Errors-in-variables Problems[J].Computational Statistics&Data Analysis, 2007,52(2):1080-1089.

[7] VAN HUFFEL S,VANDEWALLE J.The Total Least Squares Problem:Computational Aspects and Analysis [M].Philadelphia:Society for Industrial and Applied Mathematics,1991.

[8] SCHAFFRIN B,WIESER A.On Weighted Total Leastsquares Adjustment for Linear Regression[J].Journal of Geodesy,2008,82(7):415-421.

[9] SCHAFFRIN B,LEE I,CHOI Y,et al.Total Least-squares (TLS)for Geodetic Straight-line and Plane Adjustment[J].Bollettino di Geodesia e Scienze Affini,2006,65(3):141-168.

[10] SHEN Yunzhong,LI Bofeng,CHEN Yi.An Iterative Solution of Weighted Total Least Squares Adjustment[J].Journal of Geodesy,2010,85(4):229-238.

[11] FANG X.Weighted Total Least Squares:Necessary and Sufficient Conditions,Fixed and Random Parameters[J].Journal of Geodesy,2013,87(8):733-749.

[12] LIU Jingnan,ZENG Wenxian,XU Peiliang.Overview of Total Least Squares Methods[J].Geomatics and Information Science of Wuhan University,2013,38(5):505-512.(刘经南,曾文宪,徐培亮.整体最小二乘估计的研究进展[J].武汉大学学报:信息科学版,2013,38(5):505-512.)

[13] SCHAFFRIN B,FELUS Y A.On Total Least-squares Adjustment with Constraints[C]∥International Association of Geodesy Symposia.Sapporo:International Association of Geodesy,2005,128:417-421.

[14] MAHBOUB V,SHARIFI M A.On Weighted Total Least Squares with Linear and Quadratic Constraints[J].Journal of Geodesy,2013,87(3):607-608.

[15] FANG Xing.A Structured and Constrained Total Least-Squares Solution with Cross-covariances[J].Studia Geophysica et Geodaetica,2014,58(1):1-16.

[16] DE MOOR B.Total Linear Least Squares with Inequality Constraints[R].Delft:Department of Electrical Engineering, 1990.

[17] ZH ANG Songlin,TONG Xiaohua,ZH ANG Kun.A Solution to EIV Model with Inequality Constraints and Its Geodetic Applications[J].Journal of Geodesy,2013, 87(1):89-99.

[18] M A HBOUB V.On Weighted Total Least-squares for Geodetic Transformation[J].Journal of Geodesy,2012, 86(5):359-367.

[19] NOCEDAL J,WRIGHT S J.Numerical Optimization[M].New York:Springer,2006.

[20] FLETCHER R.Practical Methods of Optimization[M].New York:John Wiley&Sons,2000.

(责任编辑:陈品馨)

Weighted Total Least Squares Algorithm with Inequality Constraints

ZENG Wenxian1,FANG Xing1,LIU Jingnan1,2,YAO Yibin1
1.School of Geodesy and Geomatics,Wuhan University,Wuhan 430079,China;2.Research Center of GNSS, Wuhan University,Wuhan 430079,China

Since the inequality-constrained total least squares(ICTLS)is strongly limited due to the combinational difficulty,adjustment of the partial errors-in-variables(EIV)model which is equipped with inequality constraints is investigated.The original ICTLS problem to a standard optimization problem is reconfigured in this paper,which can be solved by existing methods such as penalty based methods.The novel ICWTLS(inequality-constrained weighted TLS)algorithm can deal with the ICTLS problem with a structure coefficient matrix and a general weight matrix,and successfully avoid the combinatorial difficulty.The examples illustrate that the new algorithm proposed in this paper is efficient and simple, which can be used in a general case in practice.

total least squares;errors-in-variables model;inequality constraints;nonlinear program

ZENG Wenxian(1975—),female,PhD, majors in the theory and method of surveying data processing.

FANG Xing

P207

A

1001-1595(2014)10-1013-06

国家自然科学基金(41474006;41404005;41231174;41174012);中央高校基本科研基金(2042014kf053)

2013-12-17

曾文宪(1975—),女,博士,主要从事测量数据处理理论与应用的研究。

E-mail:wxzeng＠sgg.whu.edu.cn

方兴

E-mail:xfang＠sgg.whu.edu.cn

ZENG Wenxian,FANG Xing,LIU Jingnan,et al.Weighted Total Least Squares Algorithm with Inequality Constraints[J].Acta Geodaetica et Cartographica Sinica,2014,43(10):1013-1018.(曾文宪,方兴,刘经南,等.附有不等式约束的加权整体最小二乘算法[J].测绘学报,2014,43(10):1013-1018.)

10.13485/j.cnki.11-2089.2014.0173

修回日期:2014-08-29