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

肖瑜

（華東交通大學理學院，江西南昌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