非線性LTS估計的截斷凝聚光滑化方法

肖瑜

（華東交通大學理學院，江西南昌330013）

非線性LTS估計的截斷凝聚光滑化方法

肖瑜

（華東交通大學理學院，江西南昌330013）

LTS估計是一類穩(wěn)健估計，其模型可以轉(zhuǎn)化為min-min-sum型非凸非光滑優(yōu)化問題，其目標函數(shù)是從m個光滑函數(shù)中取m?個的組合求和，即使m不是很大，組合數(shù)也會非常大，導致min函數(shù)的組成函數(shù)個數(shù)非常多，求解非常困難。針對這個問題的特點，并結(jié)合解一般minmax問題的截斷凝聚光滑化方法，給出了解LTS估計問題的截斷凝聚光滑化算法。數(shù)值結(jié)果表明了算法的有效性和高效率。

LTS估計；凝聚函數(shù)；截斷凝聚光滑化

在數(shù)據(jù)擬合中，最小二乘估計（the least squares estimator，簡記為LS估計），l1估計及l(fā)∞估計等方法都會受到異常點的影響。1985年，Rousseeuw提出一種穩(wěn)健估計——LTS估計［1-3］（The least trimmed squares esti?mator）。設(shè)觀測數(shù)據(jù)(ui,vi)，i=1,…,m。

其中：x為待估計的參數(shù)；ηi為對應(yīng)的殘量。那么LTS估計為

LTS估計模型雖然目標函數(shù)連續(xù)，但是非凸非光滑，求解比較困難。目前，解LTS估計大概有3類算法。一種是精確求解［1］。由于LTS估計實際上是對某個子集（包含m?個數(shù)據(jù)點的子集，特別的，對線性擬合問題取n個數(shù)據(jù)點）的最小二乘估計，可以通過對所有含m?個數(shù)據(jù)點的子集做最小二乘估計，然后取殘量平方和最小的。這種方法需要做次最小二乘估計，若m和m?稍大，總計算量就非常大，如m=40,m?=32時，需要計算7千多萬次最小二乘估計。還有一類隨機近似算法［1-3，7］，其中使用最廣泛的是Rousseeuw和Leroy于1987年提出的的PROGRESS算法，其具體步驟如下［1］：從m個數(shù)據(jù)點中隨機取m?個，記為?；對?做最小二乘估計，得到最小二乘解x*；計算x*對應(yīng)的目標函數(shù)值，即；重復以上步驟若干次，取使目標函數(shù)最小的x*為近似解。PROGRESS算法中，設(shè)重復次數(shù)為p，那么得到原問題解的概率為還有啟發(fā)式算法，如文獻［8］中運用微分進化算法來計算LTS估計。這些算法都有各自的優(yōu)缺點，第1類算法雖然能精確求解，但是計算量非常大。第2、3類算法可以通過控制重復次數(shù)來控制計算量，但是這些算法是依概率得到最優(yōu)解的，重復次數(shù)越多，概率越大，但計算量也越大，并且還是不能保證得到最優(yōu)解，甚至局部最優(yōu)解。

LTS問題式（1）可以等價變形為min-min-sum規(guī)劃問題。對于max型的非光滑函數(shù)，基于Jaynes的最大熵原理，李興斯提出了一類光滑函數(shù)，稱為凝聚函數(shù)，并給出了解min-max等問題的凝聚函數(shù)法［9］。在計算上，當max函數(shù)的組成函數(shù)很多時，凝聚函數(shù)的梯度和Hessen計算量很大，影響了凝聚函數(shù)法計算效率。在文［10-11］中，作者提出了解min-max問題的截斷凝聚光滑化算法。在每個迭代點處，在不影響算法的收斂速度的前提下，均自適應(yīng)的選取一部分函數(shù)凝聚。這樣，參與梯度、Hessen計算的函數(shù)也只有很少一部分，從而大大的減少了計算量。

采用凝聚函數(shù)光滑化min-sum型目標函數(shù)，然后用截斷凝聚光滑化牛頓法求解。記q={} 1,2,…,m，S={I∈q|#(I)=m?}，#(I)表示集合I所含元素數(shù)目。由于問題的特殊性（min-sum型函數(shù)的下標集合包含從1,…,m中取m?個數(shù)的所有可能的組合），如果直接采用［10］中的截斷方法，在每個迭代點x處，需要對所有的組合I∈S，求出殘差平方和（計算量是），再進行次函數(shù)值大小比較計算。實際上，我們可以利用問題的特殊性，只要對m個函數(shù)值求出第m?小的，即求出（計算量為O(m)），然后根據(jù)的值來設(shè)計截斷準則。

1 LTS估計的截斷凝聚光滑化牛頓法

顯然，LTS估計（1）可以寫成如下形式

問題（2）的一階最優(yōu)性條件［13-14］：

命題1設(shè)（2）在x*處取得的最小值，那么

我們采用如下凝聚函數(shù)光滑化（2）的目標函數(shù)）

其中：t＞0是凝聚參數(shù)，當t→0+時，F(xiàn)t(x)→F(x)。因為LTS估計含有大量的下標組合，文獻［10］中mini?max問題的截斷方式并不是很適合。下面，我們嘗試尋求合適的截斷方式。對給定的xˉ，取參數(shù)μ＞0，令

定義截斷凝聚函數(shù)為

接下來，我們給出FSˉt(x)與精確凝聚函數(shù)Ft(x)的函數(shù)值，梯度及Hesse陣的一些估計式。首先定義函數(shù)

上述命題的證明過程繁瑣，此處不贅述，具體過程可參閱文獻［12］。

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

推論1對任意x∈Rn,t＞0,ε1＞0,ε2＞0，若（4）中的μ按照如下方式選擇

其中：γ,ω滿足γ≥γ(x),ω≥ω(x)，則有

由推論1知，在計算過程中，可以按照式（7）選取截斷參數(shù)μ，來控制截斷凝聚誤差。接下來，我們采用截斷凝聚方式（4）～（7），結(jié)合光滑化牛頓法，得到解LTS估計問題的截斷凝聚光滑化牛頓法。具體的算法如下：

算法1（截斷凝聚光滑化牛頓法）

初始值：x0∈Rn。

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

步驟1 設(shè)i=0,k=0,s=1,xk,i=0。

征值：然后計算Cholesky因子R，使得B(xk,i)=RRT，以及計算R的倒條件數(shù)c(R)。如果c(R)≥κ1，轉(zhuǎn)步驟4，否則轉(zhuǎn)步驟5。

步驟6 計算步長λk,i=βl，其中l(wèi)≥0為滿足下面條件的最小正整數(shù)～Ftk(xk,i+λk,ihk,i)-Ftk(xk,i)≤αλk,i＜

步驟7 設(shè)xk,i+1=xk,i+λk,ihk,i,i=i+1。按照（7）和（4）計算μ和。如果，轉(zhuǎn)步驟8，否則轉(zhuǎn)步驟3。

步驟8 如果s=1，計算t*使得，轉(zhuǎn)步驟9，否則設(shè)tk+1=1/(s(k+2))，k=k+1，i=0，轉(zhuǎn)步驟2。

算法1 有如下收斂結(jié)果（具體證明過程與［10］中算法2的收斂性證明類似，此處不贅述，可參閱文獻［10，12］）：

定理1設(shè)二次連續(xù)可微，水平集Ω有界，是算法1產(chǎn)生的序列。那么，存在無窮子序列N′?N，使得當k→N′∞時，有,并且0∈?F(x*)。

2 數(shù)值結(jié)果

我們用matlab語言實現(xiàn)了截斷光滑化牛頓法（簡記為TSN）。因為［1］中的精確求解方法計算時間太長，如對例1求解時間超過24×7小時，因此我們只將TSN算法與PROGRESS算法，以及凝聚光滑化牛頓法（即在TSN算法中略去截斷技巧）進行了比較。

例1（Circle Fitting問題［15］）LTS估計不僅可以估計v=f(u,x)這類型模型，也可以估計形如f(u,x)=0的模型，如Circle Fitting問題——找一個合適的球，使得所有的數(shù)據(jù)點到球面的距離的范數(shù)最小。我們用matlab產(chǎn)生40個數(shù)據(jù)點ui∈R3，使得其分布在某個球面上，然后對所有的數(shù)據(jù)點做擾動（其中有8個數(shù)據(jù)點的擾動程度大于其他點），最后得到的數(shù)據(jù)點參閱［12］表6.3。設(shè)待求的球面的中心為o∈R3，半徑為r，則ui對應(yīng)的誤差項

例2（圓錐擬合問題［12］）這個例子的數(shù)據(jù)由matlab程序生成，只是對其中一部分數(shù)據(jù)點增加了比較大的擾動作為異常點。首先產(chǎn)生標準圓錐面上的數(shù)據(jù)點

其中：ri和γi分別是[0.1,10]和[0,2π]上均勻分布的偽隨機數(shù)。然后在z?i上加一個服從N(0,0.3)分布的擾動項，并對數(shù)據(jù)進行旋轉(zhuǎn)和平移

其中：T(?)為變換矩陣

具體的數(shù)據(jù)參閱［12］表6.4（總共30個數(shù)據(jù)點，含6個異常點）。

表1 例1的數(shù)值結(jié)果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的數(shù)值結(jié)果Tab.2 The numerical results of Example 2

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

3 結(jié)論

LTS估計是一類穩(wěn)健估計，一般采用隨機算法計算其近似解。從數(shù)值優(yōu)化的角度考慮此問題，將其模型可以轉(zhuǎn)化為min-min-sum型非凸非光滑優(yōu)化問題，并給出截斷凝聚光滑化牛頓法（TSN）。不但理論上可行，數(shù)值結(jié)果也說明了算法具有很高的計算效率。與傳統(tǒng)的算法PROGRESS相比，截斷凝聚光滑化牛頓法能在更短的計算時間內(nèi)達到更優(yōu)的目標函數(shù)值。與沒有加截斷策略的凝聚光滑化牛頓法（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]曾明華,肖瑜,黃細燕.多層次交通網(wǎng)絡(luò)的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—），女，講師，博士，主要研究方向為數(shù)值優(yōu)化。

1005-0523（2014）04-0059-06