肖瑜

（华东交通大学理学院，江西南昌330013）

非线性LTS估计的截断凝聚光滑化方法

肖瑜

（华东交通大学理学院，江西南昌330013）

LTS估计是一类稳健估计，其模型可以转化为min-min-sum型非凸非光滑优化问题，其目标函数是从m个光滑函数中取m͂个的组合求和，即使m不是很大，组合数也会非常大，导致min函数的组成函数个数非常多，求解非常困难。针对这个问题的特点，并结合解一般minmax问题的截断凝聚光滑化方法，给出了解LTS估计问题的截断凝聚光滑化算法。数值结果表明了算法的有效性和高效率。

LTS估计；凝聚函数；截断凝聚光滑化

在数据拟合中，最小二乘估计（the least squares estimator，简记为LS估计），l1估计及l∞估计等方法都会受到异常点的影响。1985年，Rousseeuw提出一种稳健估计——LTS估计［1-3］（The least trimmed squares esti⁃mator）。设观测数据(ui,vi)，i=1,…,m。

其中：x为待估计的参数；ηi为对应的残量。那么LTS估计为

LTS估计模型虽然目标函数连续，但是非凸非光滑，求解比较困难。目前，解LTS估计大概有3类算法。一种是精确求解［1］。由于LTS估计实际上是对某个子集（包含m͂个数据点的子集，特别的，对线性拟合问题取n个数据点）的最小二乘估计，可以通过对所有含m͂个数据点的子集做最小二乘估计，然后取残量平方和最小的。这种方法需要做次最小二乘估计，若m和m͂稍大，总计算量就非常大，如m=40,m͂=32时，需要计算7千多万次最小二乘估计。还有一类随机近似算法［1-3，7］，其中使用最广泛的是Rousseeuw和Leroy于1987年提出的的PROGRESS算法，其具体步骤如下［1］：从m个数据点中随机取m͂个，记为ℵ；对ℵ做最小二乘估计，得到最小二乘解x*；计算x*对应的目标函数值，即；重复以上步骤若干次，取使目标函数最小的x*为近似解。PROGRESS算法中，设重复次数为p，那么得到原问题解的概率为还有启发式算法，如文献［8］中运用微分进化算法来计算LTS估计。这些算法都有各自的优缺点，第1类算法虽然能精确求解，但是计算量非常大。第2、3类算法可以通过控制重复次数来控制计算量，但是这些算法是依概率得到最优解的，重复次数越多，概率越大，但计算量也越大，并且还是不能保证得到最优解，甚至局部最优解。

LTS问题式（1）可以等价变形为min-min-sum规划问题。对于max型的非光滑函数，基于Jaynes的最大熵原理，李兴斯提出了一类光滑函数，称为凝聚函数，并给出了解min-max等问题的凝聚函数法［9］。在计算上，当max函数的组成函数很多时，凝聚函数的梯度和Hessen计算量很大，影响了凝聚函数法计算效率。在文［10-11］中，作者提出了解min-max问题的截断凝聚光滑化算法。在每个迭代点处，在不影响算法的收敛速度的前提下，均自适应的选取一部分函数凝聚。这样，参与梯度、Hessen计算的函数也只有很少一部分，从而大大的减少了计算量。

采用凝聚函数光滑化min-sum型目标函数，然后用截断凝聚光滑化牛顿法求解。记q={} 1,2,…,m，S={I∈q|#(I)=m͂}，#(I)表示集合I所含元素数目。由于问题的特殊性（min-sum型函数的下标集合包含从1,…,m中取m͂个数的所有可能的组合），如果直接采用［10］中的截断方法，在每个迭代点x处，需要对所有的组合I∈S，求出残差平方和（计算量是），再进行次函数值大小比较计算。实际上，我们可以利用问题的特殊性，只要对m个函数值求出第m͂小的，即求出（计算量为O(m)），然后根据的值来设计截断准则。

1 LTS估计的截断凝聚光滑化牛顿法

显然，LTS估计（1）可以写成如下形式

问题（2）的一阶最优性条件［13-14］：

命题1设（2）在x*处取得的最小值，那么

我们采用如下凝聚函数光滑化（2）的目标函数）

其中：t＞0是凝聚参数，当t→0+时，Ft(x)→F(x)。因为LTS估计含有大量的下标组合，文献［10］中mini⁃max问题的截断方式并不是很适合。下面，我们尝试寻求合适的截断方式。对给定的xˉ，取参数μ＞0，令

定义截断凝聚函数为

接下来，我们给出FSˉt(x)与精确凝聚函数Ft(x)的函数值，梯度及Hesse阵的一些估计式。首先定义函数

上述命题的证明过程繁琐，此处不赘述，具体过程可参阅文献［12］。

对任意x0∈Rn,t0＞0,记Ω={x|F(x)≤Ft0(x0)}。由上面的命题，得到如下推论：

推论1对任意x∈Rn,t＞0,ε1＞0,ε2＞0，若（4）中的μ按照如下方式选择

其中：γ,ω满足γ≥γ(x),ω≥ω(x)，则有

由推论1知，在计算过程中，可以按照式（7）选取截断参数μ，来控制截断凝聚误差。接下来，我们采用截断凝聚方式（4）～（7），结合光滑化牛顿法，得到解LTS估计问题的截断凝聚光滑化牛顿法。具体的算法如下：

算法1（截断凝聚光滑化牛顿法）

初始值：x0∈Rn。

参数：t0＞0,t̂∈(0,1);α,β,κ1∈(0,1);θ∈(0,(1-α)/32),δ＞0,γ,ω充分大使得γ≥mx∈aRxnγ(x),ω≥mx∈aRxnω(x)；函数εa(t),εb(t)，τ(t),σ(t):(0,∞)→(0,∞)满足εb(t)≥εa(t)＞τ(t),lt→im0+τ(t)=0,ε1(t)=θτ(t),ε2(t)＞0。

步骤1 设i=0,k=0,s=1,xk,i=0。

征值：然后计算Cholesky因子R，使得B(xk,i)=RRT，以及计算R的倒条件数c(R)。如果c(R)≥κ1，转步骤4，否则转步骤5。

步骤6 计算步长λk,i=βl，其中l≥0为满足下面条件的最小正整数～Ftk(xk,i+λk,ihk,i)-Ftk(xk,i)≤αλk,i＜

步骤7 设xk,i+1=xk,i+λk,ihk,i,i=i+1。按照（7）和（4）计算μ和。如果，转步骤8，否则转步骤3。

步骤8 如果s=1，计算t*使得，转步骤9，否则设tk+1=1/(s(k+2))，k=k+1，i=0，转步骤2。

算法1 有如下收敛结果（具体证明过程与［10］中算法2的收敛性证明类似，此处不赘述，可参阅文献［10，12］）：

定理1设二次连续可微，水平集Ω有界，是算法1产生的序列。那么，存在无穷子序列N′⊂N，使得当k→N′∞时，有,并且0∈∂F(x*)。

2 数值结果

我们用matlab语言实现了截断光滑化牛顿法（简记为TSN）。因为［1］中的精确求解方法计算时间太长，如对例1求解时间超过24×7小时，因此我们只将TSN算法与PROGRESS算法，以及凝聚光滑化牛顿法（即在TSN算法中略去截断技巧）进行了比较。

例1（Circle Fitting问题［15］）LTS估计不仅可以估计v=f(u,x)这类型模型，也可以估计形如f(u,x)=0的模型，如Circle Fitting问题——找一个合适的球，使得所有的数据点到球面的距离的范数最小。我们用matlab产生40个数据点ui∈R3，使得其分布在某个球面上，然后对所有的数据点做扰动（其中有8个数据点的扰动程度大于其他点），最后得到的数据点参阅［12］表6.3。设待求的球面的中心为o∈R3，半径为r，则ui对应的误差项

例2（圆锥拟合问题［12］）这个例子的数据由matlab程序生成，只是对其中一部分数据点增加了比较大的扰动作为异常点。首先产生标准圆锥面上的数据点

其中：ri和γi分别是[0.1,10]和[0,2π]上均匀分布的伪随机数。然后在z͂i上加一个服从N(0,0.3)分布的扰动项，并对数据进行旋转和平移

其中：T(∙)为变换矩阵

具体的数据参阅［12］表6.4（总共30个数据点，含6个异常点）。

表1 例1的数值结果Tab.1 The numerical results of Example 1

其中：o0=(-1.2,1.2,-1.2,1.2,-1.2,1.2,-1.2,1.2)T,r0=1。

表2 例2的数值结果Tab.2 The numerical results of Example 2

其中：x0=(-0.1,-0.1,-2,1.5,-1.5,0.6)T。

3 结论

LTS估计是一类稳健估计，一般采用随机算法计算其近似解。从数值优化的角度考虑此问题，将其模型可以转化为min-min-sum型非凸非光滑优化问题，并给出截断凝聚光滑化牛顿法（TSN）。不但理论上可行，数值结果也说明了算法具有很高的计算效率。与传统的算法PROGRESS相比，截断凝聚光滑化牛顿法能在更短的计算时间内达到更优的目标函数值。与没有加截断策略的凝聚光滑化牛顿法（SN）相比，两种算法得到的解完全一样，但是我们的算法在计算时间远远少于SN算法。

[1]ROUSSEEUW P J.Leastmedian of squares regression[J].J Amer Statist Assoc,1984,79:871-880.

[2]ROUSSEEUW P J.Multivariate estimation with high breakdown point[C]//Mathematical Statistics and Applications,Dordrecht: Reidel,1985:283-297.

[3]ROUSSEEUW P J,LEROY A M.Robust regression and outlier detection[M].New York:John Wiley&Sons Inc,1987:87.

[4]ZAMAN A,ROUSSEEUW P J,ORHAN M.Econometric applications of high-breakdown robust regression techniques[J].Econo⁃metrics Letters,2001,71(1):1-8.

[5]YE M,HARALICK R M.Optical flow from a least-trimmed squares based adaptive approach[C]//International Conference on Pat⁃tern Recognition ICPR 2000,Barcelona:IEEE,2000:1052-1055.

[6]WANG H,SUTER D.Using symmetry in robustmodel fitting[J].Pattern Recognition Letters,2003,24(16):2953-2966.

[7]ROUSSEEUW P J,HUBERT M.Recent developments in PROGRESS[C]//L1-Statistical Procedures and Related Topics,CA:In⁃stitute of Mathematical Statistics,1997:201-214.

[8]CIZEK P.Robust estimation in nonlinear regression and limited dependent variablemodels[R].Berlin:Humboldt University of Berlin,2002:189.

[9]LI X S.An aggregate functionmethod for nonlinear programming[J].Science in China,1991,34(2):1467-1473.

[10]XIAO Y,YU B.A truncated aggregate smoothing newtonmethod forminimax problems[J].Appl Math Comput,2010,216:1868-1879.

[11]XIAO Y,YU B,XIONG H J.Truncated aggregate homotopymethod for nonconvex nonli-near programming[J].Optimization Methods&Software,2014,29(1):160-176.

[12]肖瑜.截断凝聚光滑化算法[D].大连:大连理工大学,2010.

[13]ROLEFF K.A stablemultiple exchange algorithm for linear SIP[C]//Lecture Notes in Control and Information Science,Berlin: Springer,1979:83-96.

[14]CLARKE F H.Optimization and nonsmooth analysis[M].New York:John Wiley&Sons Inc,1983:56.

[15]GANDER W,GOLUB G H.TREBEL R.Least-squares Fitting of Circles and Ellipses[J].BIT Numerical Mathematics,1994,34 (4):558-578.

[16]曾明华,肖瑜,黄细燕.多层次交通网络的UE与SO混合均衡与效率损失[J].华东交通大学学报,2012,29(2):57-62.

Truncated Aggregate Smoothing Method for Nonlinear LTS Estimator

Xiao Yu
(School of Basic Science,East China Jiaotong University,Nanchang 330013,China)

The computing of the nonlinear least trimmed squares(LTS)estimator is considered.LTS is a robust esti⁃mator and can be converted to amin-min non-convex and non-smooth programming problem.For the data set with sizem,the objective function is theminimum of all them͂-subsets'residual sum of squares.Even ifmis not big,the number of the subsetsmay be very large whichmakes computing LTS estimator difficult.For such a special kind of problem,an appropriate truncated criteria standard is given and then an efficient truncated smooth⁃ing Newtonmethod is proposed.The numerical results show the efficiency.

LTS estimator;aggregate function;truncation smoothing

O221.2

A

2014-03-07

国家自然科学基金资助项目（11171051，11261019）

肖瑜（1981—），女，讲师，博士，主要研究方向为数值优化。

1005-0523（2014）04-0059-06