基于殘差向量l1范數(shù)最小化與基追蹤的多元線性模型參數(shù)估計方法

馮志強(qiáng)，張鴻燕

（1.海南師范大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，海南 ?？?571158；2.海南師范大學(xué) 信息科學(xué)技術(shù)學(xué)院，海南 海口 571158）

多元線性模型在科學(xué)、技術(shù)、經(jīng)濟(jì)、管理、軍事等領(lǐng)域中有著廣泛的應(yīng)用，估計多元線性模型的參數(shù)是一個在理論上與實踐上均有重大價值的實際問題。設(shè)X1,X2,…,Xn是n個自變量（解釋變量），p1,p2,…,pn是模型參數(shù)，Y是因變量，則多元模型可以表示為

當(dāng)d大于臨界值dt時，認(rèn)為點V對于該超平面即為外點（outliers）（見圖1）。需要注意的是：dt通常是由實際的物理觀測過程與問題的精度要求所決定的，某些情況下則難以明確dt的經(jīng)驗數(shù)值。

圖1 外點示例Figure 1 Demonstration of outliers

各種最小二乘方法適用的前提是假定噪聲服從正態(tài)分布（Gauss-Markov 假設(shè)）而且樣本數(shù)據(jù)中不存在外點（outliers）。如果噪聲不服從正態(tài)分布，采用以上方法估計會存在較大的誤差；如果觀測數(shù)據(jù)中存在外點，則參數(shù)估計的性能會很差甚至完全不能用。隨機(jī)抽樣一致性（RANdom SAmple Consensus，RANSAC）方法很好地解決了外點帶來的問題[7-8]，它在外點的比例小于50%的條件下都可以使用。RANSAC方法是一種隨機(jī)算法，在攝影測量與計算機(jī)視覺領(lǐng)域得到了廣泛使用，但是在其他領(lǐng)域受關(guān)注的程度不高，而且這種迭代型的算法的時間復(fù)雜度比較高，不適合實時計算。

如果同時考慮噪聲的統(tǒng)計分布偏離正態(tài)分布的情形以及存在外點的情形，線性模型參數(shù)估計的優(yōu)化準(zhǔn)則應(yīng)當(dāng)選擇s= 1對應(yīng)的l1范數(shù)作為誤差的代價函數(shù)。采用l1范數(shù)最小化方法估計線性模型參數(shù)的難點在于代價函數(shù)不再是可微的，求解比較困難，但是估計的精度與穩(wěn)健性很好。

在l1范數(shù)最小化的意義下估計線性模型（1）的參數(shù)向量p歸結(jié)為求解超定線性方程組Ap=b，相應(yīng)的最優(yōu)化問題是

該問題的解稱為方程組Ap=b在l1范數(shù)最小化意義下的解，簡稱為方程組的l1范數(shù)解。相應(yīng)地，問題（P0）簡稱為l1范數(shù)最小化問題。

求解問題（P0）的方法已有不少，主要有線性規(guī)劃方法[9]、投影下降法[10]、有效集法[11]、區(qū)間方法[12]等。這些方法各有優(yōu)缺點，如線性規(guī)劃方法使問題規(guī)模擴(kuò)大，投影下降算法結(jié)構(gòu)復(fù)雜性高，不易被工程技術(shù)人員掌握，有效集法算法簡單直觀，但是無法處理矩陣A行虧秩的情況，區(qū)間方法過于復(fù)雜并且無最優(yōu)解區(qū)間的檢驗，比較麻煩。文獻(xiàn)[13]針對問題（P0）提出了基于最小l1殘差向量的數(shù)值解法，它將求解過程轉(zhuǎn)化為一個約束不可微最優(yōu)化問題的求解與相容線性方程組的求解。遺憾的是，該方法只適用于系數(shù)矩陣A是滿秩的情況，即rank(A)=n。

本文進(jìn)一步完善了文獻(xiàn)[13]中的方法，得到基于基追蹤與Moore-Penrose 廣義逆的新型算法，該算法具有如下優(yōu)點：增添了對系數(shù)矩陣A缺秩情況的處理，解除了對于滿秩條件的限制；算法結(jié)構(gòu)簡單、易于操作、直觀性強(qiáng)并且適合處理大規(guī)模問題。

為了后續(xù)討論的方便，這里給出將會用到的一些數(shù)學(xué)記號及其解釋（見表1）。

表1 部分?jǐn)?shù)學(xué)符號說明Table1 Description of some mathematical notations

1 最小殘差向量的性質(zhì)及等價形式

最小l1殘差向量的定義與性質(zhì)如下[14]：

定義1（最小l1殘差向量） 設(shè)popt為問題（P0）的解，稱ropt=Apopt-b為問題（P0）的最小l1殘差向量。如果能求得最小l1殘差向量，那么相容線性方程組Ap=b+ropt的解即為問題（P0）的解。

定理1 設(shè)有兩個超定線性方程組Ap=b與Bp=b，若存在可逆陣P使得B=Ap，則這兩個超定線性方程組的最小l1殘差向量相同。

引入最小l1殘差向量的概念與基本性質(zhì)的目的是確定超定線性方程的l1范數(shù)解的結(jié)構(gòu)與表示，達(dá)成這個目標(biāo)的關(guān)鍵是建立以下的定理。

于是定理2得證。

需要說明的是：為了提高本文提出的多元線性模型參數(shù)估計方法的適應(yīng)性，本文對文獻(xiàn)[13]給出的結(jié)果做了改進(jìn)，即在式（9）與式（11）中用矩陣的Moore-Penrose偽逆A?1代替了文獻(xiàn)[13]中所采用的矩陣逆A-11。

2 求解l1 范數(shù)最小化問題的稀疏優(yōu)化方法

2.1 稀疏優(yōu)化的等價轉(zhuǎn)換

對于最小l1殘差向量ropt，可以通過現(xiàn)有的優(yōu)化方法求得它的一個稀疏解，即使得r中分量取0值的數(shù)目盡可能多。換句話說，求解問題（P1）就轉(zhuǎn)換為l1范數(shù)最小化的稀疏優(yōu)化問題，這可以通過基追蹤方法來求解[15]。

2.2 算法描述

綜合上述推導(dǎo)過程，可以得到問題（P0）的求解算法，相應(yīng)的偽代碼描述如算法1所示：

算法1：超定線性方程組最小l1 范數(shù)求解算法輸入：A ∈Rm × n,b ∈Rm ×1,m ＞n ≥2輸出：p ∈Rn ×1 1.初始化m為矩陣A的行數(shù)，n為矩陣A的列數(shù)；2.如果m ＜n輸出錯誤信息：“矩陣A的行數(shù)必須大于列數(shù)”；3.取矩陣A的前n行構(gòu)成矩陣A1：A1 ←A(1:n)；4.取矩陣A的n + 1～m行構(gòu)成矩陣A2：A2 ←A(n + 1:m)；5.計算矩陣A1的Moore-Penrose逆：A?1 ←CalcMoorePenroseInverse（A1）；6.計算矩陣:C ←A2A?1；7.計算向量:w ←Cb(1:n) - b(n + 1:m)；8.設(shè)定矩陣:D ←[-C Im - n]；9.設(shè)定矩陣:Φ ←[D -D]；10.構(gòu)建線性規(guī)劃問題目標(biāo)函數(shù)中的系數(shù)向量:c ←[1,1,...,1]T ∈R2m ×1；11.求解標(biāo)準(zhǔn)線性規(guī)劃:αopt ←LinprogSolver(c,Aeq = Φ,beq = w,lb = 02m × 1)；12.求得最小模剩余向量:ropt ←αopt(1:m) - αopt(m + 1:2m)；13.計算矩陣A的Moore-Penrose逆：A? ←CalcMoorePenroseInverse（A）；14.計算最優(yōu)參數(shù)向量：popt ←A?1(b + ropt)；15.返回popt，計算結(jié)束。

偽代碼中的第11行用到了標(biāo)準(zhǔn)線性規(guī)劃問題的求解器，第5行與第13行用到了Moore-Penrose 逆的求解器。前者有大量成熟的標(biāo)準(zhǔn)工具與軟件包，后者也有備受歡迎的Greville列遞推算法等多種方法[1，17]。

2.3 誤差分析

若矩陣A行缺秩，由式（9）可知其并不影響算法對于問題（P0）的求解，故在這里只考慮矩陣A列缺秩的情況所帶來的誤差。

若矩陣A的列秩為rank(A)=k＞0，即列秩虧缺為n-k，則經(jīng)過一定次數(shù)的初等變換后，可以得到如下關(guān)系：

3 驗證與確認(rèn)

3.1 滿秩情況

問題：已知

構(gòu)造約束最優(yōu)化問題

即為超定線性方程組Ap=b的最小l1范數(shù)解。

3.2 列缺秩情況

在MATLAB 環(huán)境中分別構(gòu)造線性模型的系數(shù)矩陣A∈Rm×n與向量u∈Rn×1，其中A、u中的每個元素均服從標(biāo)準(zhǔn)正態(tài)分布[18]。接著對A中的元素做出如下轉(zhuǎn)換：

A(:,n-i) =A(:,1:k)fi,i= 0,1,…,n-k- 1， （30）

其中fi每次均隨機(jī)產(chǎn)生，構(gòu)造的矩陣A為列秩為k的矩陣（其中列缺秩為n-k）；最后令線性模型的右端向量b=Au，u即為符合式（24）的精確解。

令m= 10,n= 5，對列缺秩矩陣A構(gòu)造的超定線性方程組進(jìn)行求解（算法1）所得到的部分測試數(shù)據(jù)如表2所示。

表2 列缺秩矩陣A部分測試數(shù)據(jù)Table 2 Partial test data of column deficient rank matrix A

從數(shù)據(jù)中可以得到如下關(guān)系（考慮到計算機(jī)精度原因，其實際求解精度的數(shù)量級為10-15）：

將式（31）代入式（27）可知此時β= 0，且有

3.3 外點情況

不失方法上的一般性，取m= 100,n= 30，設(shè)定的理想?yún)?shù)向量p～1,p～2,p～3∈R30×1如表3所示：

表3 測試所用3種參數(shù)向量類型Table 3 Three types of parameter vector used in the test

在MATLAB 環(huán)境中構(gòu)造系數(shù)矩陣A∈R100×30，其中n= 30,m= 100,A中的每個元素均服從標(biāo)準(zhǔn)正態(tài)分布[18]。距離臨界值取為dt= 0.15。令bi=Ap～i,i= 1,2,3，在數(shù)據(jù)集Ω的y分量中添加強(qiáng)正態(tài)噪聲ξ～N(0,5)使得點集合

如圖2 所示，取ηoutlier=20%，對右端向量b中20%的元素進(jìn)行擾動，圖2 橫坐標(biāo)為數(shù)據(jù)點的分量序號i，縱坐標(biāo)的值為右端向量b中yi分量的值(1 ≤i≤100)，添加強(qiáng)噪聲擾動ξi后的分量坐標(biāo)yi+ξi用虛線表示與加粗的點表示。

圖2 在右端向量中添加強(qiáng)噪聲使得數(shù)據(jù)點集Ω中20%的點成為外點Figure 2 Adding strong noise to the component y in b for the data set Ω such that the percentagte of outlers is 20%

本算例中，用于參數(shù)估計的點對是100個，25%比例的外點意味著有25個外點。從圖3中誤差曲線的趨勢可以看出，算法是適用于各種數(shù)據(jù)類型的。從相對誤差的數(shù)量級可以看出：最小l1范數(shù)解在外點比例不超過25%時具有非常好的穩(wěn)健性，參數(shù)估計的相對誤差的數(shù)量級控制在10-10以下，在外點比例接近27%～28%時，相對誤差的數(shù)量級急劇增加到10-2；對于傳統(tǒng)的最小二乘法，在只有一個外點的情況下其相對誤差的數(shù)量級為10-1，此時參數(shù)估計的精度通常難以滿足應(yīng)用需求，隨著外點比例的增加，相對誤差接近于1，這意味著誤差向量與精確解的二范數(shù)或能量基本相同，得出的參數(shù)估計值已經(jīng)完全失效了。

圖3 不同參數(shù)類型下相對誤差與外點比例的關(guān)系Figure 3 Dependence of relative error on proportion of outliers for different types of parameters

若定義圖3中最小l1范數(shù)解的相對誤差增加到10-3時對應(yīng)的外點比例為臨界誤差外點比例，記為ηc，定義m n為數(shù)據(jù)的冗余程度，則在不同數(shù)據(jù)冗余度下的臨界誤差外點比例如圖4所示：

圖4 臨界誤差外點比例與數(shù)據(jù)冗余度的關(guān)系0.5Figure 4 Dependence of critical proportion of outliers on data redundancy

從圖4中可以看出：臨界外點比例ηc與數(shù)據(jù)冗余度m n并非單純的線性關(guān)系。隨著數(shù)據(jù)冗余度m n的增大，臨界外點比例ηc也對應(yīng)增大（算法的性能更好）。尤其是當(dāng)數(shù)據(jù)樣本大于100（m＞100）時，較低的數(shù)據(jù)冗余度（m n＜5）也能使得算法求解的臨界外點比例ηc＞50%，而基于最小二乘法的RANSAC方法對應(yīng)的臨界外點比例必須滿足要求ηc＜50%。這表明基于l1范數(shù)的參數(shù)估計方法的穩(wěn)健性優(yōu)于RANSAC方法。

4 結(jié)論

本文利用殘差向量l1范數(shù)最小化與基追蹤準(zhǔn)則給出了求解多元線性模型參數(shù)估計的最小l1范數(shù)解的一種穩(wěn)健的新方法，通過相容線性方程組的引入極大地簡化了原問題的求解步驟，并且針對參數(shù)矩陣A缺秩情況下算法的可行性進(jìn)行了驗證，證實了算法的有效性。數(shù)值仿真實驗表明：在包含稀疏噪聲的情況下，最小l1范數(shù)解優(yōu)于最小二乘解；當(dāng)噪聲導(dǎo)致出現(xiàn)一定比例的外點時，最小二乘法則失效不能使用，但是本文提出的新方法在普遍的情形下當(dāng)外點比例不超過25%時能有效地恢復(fù)真實的參數(shù)數(shù)值，而在數(shù)據(jù)冗余性比較高的場合，允許的臨界外點比例超過50%。