文章快速检索  
  高级检索
附有不等式约束的加权整体最小二乘算法
曾文宪1, 方兴1, 刘经南1,2, 姚宜斌1    
1. 武汉大学 测绘学院,湖北 武汉 430079;
2. 武汉大学 卫星导航定位技术研究中心,湖北 武汉 430079
摘要:针对现有附有不等式约束的整体最小二乘算法的缺陷,以partial EIV(errors-in-variables)模型为基础,在整体最小二乘准则下,通过将附有不等式约束的EIV模型的求解转换为标准的附有不等式约束的最优化问题,并采用惩罚函数法等方法得到了附有不等式约束的加权整体最小二乘新算法。新算法将现有算法的特殊权阵限制条件扩展到了一般性权矩阵,将要求系数矩阵元素全部随机的限定条件扩展到了可同时包含随机和非随机元素的一般情况,并且新算法解决了现有算法计算量受制于约束方程数量的缺陷。实例计算表明,本文提出的算法简单、有效,具有普遍适用性。
关键词整体最小二乘估计     EIV模型     不等式约束     非线性算法    
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
First author: ZENG Wenxian(1975—), female, PhD, majors in the theory and method of surveying data processing. E-mail: wxzeng@sgg.whu.edu.cn
Corresponding author: FANG Xing E-mail: xfang@sgg.whu.edu.cn
Abstract: 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.
Key words: total least squares     errors-in-variables model     inequality constraints     nonlinear program    
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表示yn×1随机误差向量,且ey期望为0、协因数阵Qy=w-1σ2(wσ2分别表示观测向量的权阵以及单位权方差);β表示t×1的参数向量;A表示n×t的系数矩阵;EA表示An×t随机误差矩阵,EA的期望为0且协因数阵为QA=ω-1σ2(ω表示系数矩阵观测元素的权阵),通常假定EAey不相关;G表示不等式约束方程的k×t系数矩阵;z表示k×1的常数向量。

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

式中,仅将A中的随机量提取出来构成m×1估计量a参与平差计算,a表示相应于a的观测向量,ea表示am×1随机误差向量;(h+Ba)=vec(AEA),vec为向量化符号,表示将n×t矩阵(AEA)按列顺序排列为nt×1的列向量;hnt×1列向量,对应于vec(AEA)中的非随机项,其中随机项部分取0,Ba对应于vec(AEA)中的随机项,B是由常数构成的nt×m矩阵。

式(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)算法

非线性ICPEIV模型(2)的待估参数为βa,根据整体最小二乘准则

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

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

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

罚函数法的基本思想是通过对目标函数式(4)中的f(βa,)增加一个二次型惩罚项,将式(4)转化为无约束条件的目标函数形式

式中,μ表示惩罚参数(μ>0);ci(β)=max[0,-Giβ+zi](i=1,2,…,k);GiG的第i行;ziz的第i个元素。式(5)存在一阶导数且一阶导数连续。

当式(5)中的二次型惩罚项为不等于0的正数时,参数解在可行解范围之外;当惩罚项等于0时,表明参数解必定落在满足式(4)的可行解区域范围内。当惩罚参数μ增大并趋于无穷时,表示在目标函数中加大了约束条件控制程度。因此,从较小的μk值(k=1,2,…,n)开始并逐步增大其数值,将式(5)对待估参数βa分别求一阶偏导并令其等于0,从而构成方程组或者采用牛顿法等无约束最优化方法均可求得模型的估计值,其中对应于最小的fICWTLS的参数解即为ICPEIV模型的最优解。

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

(1) 给定初始值:μ0>0及τ>0,β00(可采用不考虑不等式约束方程情况下的式(1)的TLS解作为初始值)。

(2) 固定μ0值,根据式(5)估计,直到参数前后两次估计值小于设定的τ

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

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

(5) 系列μ值下,最小fICWTLS值对应的参数解即为最优估计值

对于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的全部元素为随机量,Ay中所有观测元素不相关且中误差均为0.1,其中参数真值表 2。为了改进模型的估计结果,利用参数的先验信息设计了18个不等式约束方程(见表 1),其中,0≤βi≤0.8(i=1,2,3,4)以及-0.5≤βj≤0(j=5,6,7)可分别表示为

(I4⊗[-1, 1]T)β≥([1,1,1,1]T⊗[0,0.8]T)
(I3⊗[-1, 1]T)β≥([1,1,1]T⊗[-0.5,0]T)
式中,I4I3分别表示4阶和3阶的单位阵。

表 1 附有不等式约束的平差模型数据 Tab. 1 Data set of the inequality-constrained adjustment model
Ay
-0.841 7-5.006 60.874 6-1.010 70.954 84.982 7-4.972 4-2.967 4
-2.024 9-3.025 61.682 7-2.571 32.078 36.761 6-3.754 0-4.084 9
-2.870 6-0.959 52.453 9-3.846 03.123 78.486 4-2.621 8-4.897 2
-3.963 81.035 23.275 6-5.090 94.107 79.958 3-1.320 4-5.617 9
-5.026 33.080 84.184 5-6.616 44.996 511.820 0-0.351 8-6.405 0
-5.891 84.974 44.838 5-8.006 75.949 613.519 60.892 5-7.291 5
-6.901 87.071 56.024 6-9.500 26.877 615.184 91.892 8-7.999 7
-8.011 19.104 96.673 9-10.855 68.011 916.989 33.452 1-8.952 5
-8.870 411.017 87.193 2-12.335 29.095 918.563 64.500 9-9.895 1
-9.821 913.015 38.084 5-13.563 69.965 920.217 45.774 7-10.776 2
Gz
1.253-0.37802.7840-4.2010.5413.510 0
03.251-0.23601.526-0.69300.970 0
-2.54301.48700.325-0.98701.130 0
3.4580-0.1452.1240-1.2563.1421.620 0
 0≤βi≤0.8 (i=1,2,3,4)    -0.5≤βj≤0 (j=5,6,7)

ICWTLS算法估计值ICWTLS以及基于穷举法的ICTLS算法估计值ICTLS表 2,表中同时计算了不考虑约束方程时的TLS解TLS,可以看到:

表 2 ICWTLS算法估计结果 Tab. 2 Estimation of ICWTLS
参数ICWTLSTLSICTLS
β10.200 00.152 6-0.616 60.152 6
β20.300 00.260 00.664 60.260 0
β30.800 00.696 1-0.451 30.696 1
β40.500 00.504 20.453 00.504 2
β5-0.100 0-0.032 8-1.452 6-0.032 8
β6-0.500 0-0.500 00.047 1-0.500 0
β7-0.200 0-0.160 9-0.360 1-0.160 9

(1) 由于模型数据等精度、系数矩阵全部为随机量,因此,采用本文算法求得的参数解ICWTLS和基于穷举法的参数解ICTLS相同,均满足全部不等式约束条件
GICWTLSz=[0.088 4 0.911 6 0.404 8 0.595 2 0.975 3 0.024 7 1.0 0 0.936 5 0.063 5 0.661 4 0.338 6 0.772 9 0.227 1 0 0 0.137 3 0]T≥0

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

(3) 若不顾及约束条件,采用TLS算法得到的参数解TLS代入不等式约束方程的结果向量如下,其中第1、3、14、15、16、18个元素均小于0,不满足不等式约束条件式
GTLSz=[-0.385 8 1.385 8 -0.039 0 1.039 0 0.970 7 0.029 3 0.714 1 0.285 9 0.733 1 0.266 9 0.491 4 0.508 6 1.103 7
-0.103 7 -0.329 3 -1.634 5 1.437 9 -0.993 4]T

(4) 由于不等式约束条件利用了参数的先验信息,因此本实例的ICWTLS显然比TLS更接近于参数真值,即,当平差模型存在可靠的参数的先验信息时,如参数在确定范围内取值、数字地形拟合模型边界范围的设定等,平差算法若能将这些先验信息作为约束条件加入计算过程中,可以有效改善模型的估计结果或者使得估计结果能够满足模型设定的条件。

表 2中,表示参数真值;ICWTLS
表示采用本文ICWTLS算法的参数估计值;TLS表示不考虑不等式约束方程的TLS算法的参数估计值;ICTLS表示采用基于穷举法的ICTLS算法的参数估计值。

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 =[b1 b2bn]Tcn×1 =[c1 c2cn]Tebn×1=[eb1 eb2ebn]Tecn×1 =[ec1 ec2ecn]T

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

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

表 3 附有不等式约束的平面拟合模型数据 Tab. 3 Data set of the inequality-constrained plane fitting model
点号AyGz
111.407(0.3)6.700(0.3)-9.930(0.3)1001.750 0
212.404(0.3)7.058(0.3)-10.769(0.3)-100-1.850 0
312.610(0.3)7.947(0.3)-11.487(0.3)0101.450 0
414.817(0.5)8.993(0.5)-11.190(0.5)0-10-1.550 0
515.334(0.5)9.691(0.5)-12.312(0.5)
616.115(0.5)12.046(0.5)-13.356(0.5)
716.886(0.5)12.322(0.5)-13.239(0.5)
817.317(0.7)13.053(0.7)-14.783(0.7)
919.067(0.7)14.794(0.7)-15.295(0.7)
1019.818(0.7)14.977(0.7)-17.204(0.7)

表 4 附有不等式约束的平面拟合模型的ICWTLS估计结果 Tab. 4 Estimation of inequality-constrained WTLS and TLS of the adjustment model
参数β1β2β3
ICWTLS1.7501.550-2.170
TLS4.4182.168-2.708
5 结 论

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

参考文献
[1] XU P L, LIU J N, SHI C. Total Least Squares Adjustment in Partial 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 Least-squares 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] ZHANG Songlin, TONG Xiaohua, ZHANG Kun. A Solution to EIV Model with Inequality Constraints and Its Geodetic Applications[J]. Journal of Geodesy, 2013, 87(1): 89-99.
[18] MAHBOUB 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.
http://dx.doi.org/10.13485/j.cnki.11-2089.2014.0173
中国科学技术协会主管、中国测绘地理信息学会主办。
0

文章信息

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

文章历史

收稿日期:2013-12-17
修回日期:2014-08-29

相关文章

工作空间