分布式變量選擇—MCP 正則化

王格華, 王璞玉, 張 海

(西北大學數學學院，西安 710069)

1 引言

考慮線性回歸Y=Xβ+ε，其中X ∈RN×p是輸入變量，Y ∈RN是輸出變量，β ∈Rp是該模型的系數向量，ε是隨機誤差.當模型稀疏時，即真實β的大部分分量等于零，我們考慮用稀疏正則化方法來解決此問題.通常，正則化方法具有如下形式

其中Y= [Y1,Y2,··· ,YN]T ∈RN, X= [X1,X2,··· ,XN]T ∈RN×p，λ ＞0 是調控參數，用來控制模型的復雜度.當k= 0 時，上述正則化框架對應于AIC 或BIC[1,2]，此時它被稱為L0正則化.雖然L0正則化可以得到最稀疏的解，但求解它對應于一個NP 難問題.因此，眾多學者提出了不同的正則化方法來逼近L0正則化，其中包括Lasso 方法(當k= 1 時)[3-7].為了得到更稀疏的解，許多學者提出了非凸正則化方法，包括Adaptive Lasso[8]、Smoothly Clipped Absolute Deviation(SCAD)[9]、Minimax Concave Penalty(MCP)[10]、L1/2正則化[11]等.特別地，MCP 方法的提出作為一種經典的非凸正則化方法，具有許多理想的性質，并且已經被廣泛研究[12-15].它具有以下形式

其中λ ≥0, γ ＞1.

上述正則化方法的提出和發展為我們處理和分析海量高維數據提供了有效工具，但其僅適合于單機數據處理，即數據存于同一個機器中.隨著信息的發展，大量的生物科學和醫學數據由于數據量巨大而采用分布式存儲方式.因此，設計和研究適合于處理分布式數據的機器學習算法是當前數據分析關注點之一.針對分布式數據存儲的新特點，Mateos 等[16]提出了將分布式計算與正則化方法相結合的分布式Lasso 方法，Wang 等[17]提出了分布式L1/2正則化方法.鑒于非凸正則化對變量選擇和特征提取的優越性，我們關注基于分布式計算的非凸MCP 正則化方法.

本文研究分布式MCP 方法.首先，依數據所屬計算機的不同，提出分布式MCP 模型；其次，基于ADMM 算法，提出分布式MCP 算法并證明它的收斂性.最后，通過模擬實驗和真實數據實驗，證明所提方法在處理海量分布式存儲的數據的有效性.

2 分布式MCP 模型

模型(3)將所有數據拆為J個部分，其與原模型(1)等價.

為了使每臺計算機能夠不依賴其他計算機獨立的進行計算，我們考慮用局部變量{βj}替換全局變量β，其中βj為每臺計算機的局部估計.那么模型(3)可改寫為下述分布式MCP 模型

3 分布式MCP 算法

本節基于ADMM 算法求解分布式MCP，并給出算法的收斂性分析.

3.1 ADMM 算法

ADMM 算法于20 世紀70 年代早期被提出[18,19]，隨后被眾多學者推廣[20-24].近年來，ADMM 算法在計算機視覺，信號處理等諸多領域被廣泛應用.本節借助ADMM 算法求解分布式模型(4).首先，考慮添加輔助局部變量，則模型(4)可以改寫為

模型(5)要求網絡中所有局部變量彼此相等，注意到我們假設圖G= (J,E)是連通的，那么(3)和(5)是等價的.

為了方便表述，我們將模型(5)寫成向量的形式

其中

其中A*=[Ip×p,··· ,Ip×p]T ∈R(J×p)×(p),0∈R(J×p)×p是零矩陣，A1=[A*,0,··· ,0],A2=[0,A*,··· ,0],··· ,AJ=[0,··· ,0,A*].

模型(6)的二次增廣拉格朗日函數為

其中V ∈R(J×J×p)×1為拉格朗日乘子，c ＞0 為預選參數.

ADMM 算法分為以下三個步驟：

步驟1 更新Z

步驟2 更新β

步驟3 更新拉格朗日乘子

在分析ADMM 算法的收斂性之前，我們給出一些記號.令?f表示函數f的次梯度，M表示矩陣ATA的最小正特征值，M′表示矩陣BTB的最小正特征值.如果一個函數可微且梯度是李普希茲連續的，我們稱這個函數李普希茲可微.顯然，模型(6)中的f(β)是李普希茲可微的，且李普希茲常數為Lf.對于任何Z1, Z2，都有下式成立，存在ω ≥1，

對于[Z2]的每個維度[Z2]i，若[Z2]i= 0，則[Z2]i的廣義次梯度在[-1,1]之間，如果ω足夠大，那么對于任何Z1，我們都有

若[Z2]i/=0，則[Z2]的二階梯度為0 或-1/γ，那么

故(11)成立.

下面，我們給出幾個引理和主要的收斂定理.

引理1 對于所有的k，下面的結論成立：

且

由于M和M′分別為ATA和BTB的最小正特征值，則(A2)和(A3)顯然成立.

(A4)可由下式得到

其中，利用引理1(A1)，(a)成立；注意到f是李普希茲可微的，(b)成立.

引理2 若c ＞max{Lf/2M′,LfM′,ω/JM}，則當k →∞時，L(Zk+1,βk+1,V k+1)收斂.

證明 首先，我們證明在選取合理的預選參數c的情況下，L(Z,β,V)在每次ADMM迭代期間遞減.

我們將增廣拉格朗日函數之差拆分為

考慮第一項其中，利用引理1(A1)和f為李普希茲可微的，(f)成立；利用引理1(A3)和引理1(A4)，(g)成立.由上述兩項得到的不等式，ADMM 算法中的迭代滿足下式

如果c ＞max{Lf/2M′,ω/JM}，則在ADMM 迭代期間L(Z,β,V)單調遞減.

下面，我們證明增廣拉格朗日函數對任意k都有界.

對于任意一個Zk+1，存在β′，使得AZk+1+Bβ′=0，那么

其中，在(h)中我們假設c ＞LfM′，則L(Zk+1,βk+1,V k+1)有下界，進而當k →∞時，L(Zk+1,βk+1,V k+1)收斂.

引理3 存在dk+1∈?L(Zk+1,βk+1,V k+1)，使得當k →∞時，→0.

證明 首先，有

因為?V L=AZk+1+Bβk+1=(1/c)(V k+1-V k)，且由引理1(A4)，不等式

成立.同時，我們有

又因為

其中，利用Zk+1的最優性條件，(i)成立，則有

基于上述三個引理，我們給出ADMM 算法收斂性定理.

定理2 若c ＞max{2Lf/M′,LfM′,ω/M}，則由ADMM 算法(8)—(10)產生的序列收斂到L(Z,β,V)的駐點.

3.2 參數估計

為了求解ADMM 算法的三個步驟，我們改寫增廣拉格朗日函數

首先，考慮{Zj}的更新，令

根據數據所屬計算機的不同，我們將(8)分解為J個子問題

下面，我們給出坐標下降法的收斂性分析.

引理4 對于m= 1,2,··· ,p，令Oj,m為目標函數O關于變量Zj的坐標m的函數，Zj其它坐標是固定的.如果預選參數c ＞1/(γ*J)，則對于所有坐標m，Oj,m是關于[Zj]m的一個凸函數.

證明 假設預選參數c ＞1/(γ*J)，則對于每個計算機j，和所有維度m, O是關于[Zj]m的嚴格凸函數.O的嚴格凸性和連續性滿足文獻[26]中定理4.1 和定理5.1 的條件，進而坐標下降算法收斂到坐標最小值.由于O的所有方向導數都存在，因此每個坐標最小值也是局部最小值.

從引理4 中我們可以看出，當預選參數c ＞1/(γ*J)時，對于任何一個j ∈J，盡管目標函數O包含非凸的MCP 罰項，它仍是一個關于[Zj]m的凸函數.在XTX和c滿足某些條件時，目標函數L(Z,β,V)為凸函數.在這種特殊情況下，(15)收斂于目標函數(13)的全局最小值.

下面，我們考慮β的更新，它也可以分解為J個子問題

同樣，通過使用坐標下降法，給出局部βj的第m個坐標的顯式解

注1 如果計算機i與計算機j不連接，則Eij=0.這意味著我們在更新Zj和βj時，只需要計算機j的鄰居的信息.

結合(15),(17)和(9)，我們提出如下分布式MCP 算法.

算法1 分布式MCP 算法

對 ?j ∈J, 令β =(1,··· ,1)T, Z =(1,··· ,1)T, Vij =(1,··· ,1)T,局部運行for k =1,2,···第j 臺計算機在Nj 與鄰居傳遞信息對于m=1,2,··· ,p利用(15)更新[Zk+1 j ].end for l=1,2,··· ,p利用(17)更新[βk+1 j ].end通過(10)更新拉格朗日乘子V.end

3.3 分布式k 折交叉驗證

下面，我們給出分布式k折交叉驗證算法，以此來選擇懲罰參數λ.

重復上述步驟，我們可以計算每個λi的均方誤差MSE，使得

最終，我們選擇預測誤差最小的λdcv=arg minλi{e(λi)}作為調控參數.

4 實驗

本小節，我們通過4 個實驗來說明算法的有效性.其中，實驗1 通過模擬數據說明算法與非分布式算法有相同的變量選擇結果.實驗2、實驗3 通過設計大數據量數據說明算法適合于大數據的分析處理.實驗4 通過實際數據說明算法的實用性.

4.1 實驗1(模擬實驗)

在這個實驗里，對于我們考慮的線性模型

我們設置變量數量p=8，具體地β=(3,1.5,0,0,0.5,0,0,0), X=(X1,X2,··· ,X8).Xj～N(0,Σ),Σ = (σij), σij= 0.5|i-j|,1≤i,j ≤n.當應用分布式MCP 方法時，我們將所有數據劃分為J=5 組，并選擇鄰接矩陣

我們模擬了100 個數據集，每個數據集包含100 個觀測值.我們分別應用Lasso 和分布式MCP 處理這100 個數據集.其中，我們用LARS 求解L1正則化.當應用分布式MCP 方法時，用分布式5 折交叉驗證來選取調節參數λ，預選參數c= 0.1.我們記錄100 組數據集上的正確識別零系數的平均數(C.A.N)，以及錯誤識別的零系數的平均數(IC.A.N)，結果如表1 所示.

從表1 中可以看出，非分布式MCP 與分布式MCP 有相同的變量選擇結果，其C.A.N是4.45，多于Lasso 的C.A.N 的2.26，這表明，分布式MCP 的變量選擇結果比Lasso 稀疏.C.A.N 為在100 個數據集上，結果模型和真實模型中值均為0 的系數個數的平均值；IC.A.N 為在100 個數據集上，結果模型中值為零但在真實模型中值為非零的系數個數的平均值.

表1 Lasso 和分布式MCP 在實驗1 中的結果

4.2 實驗2

在這個實驗中，我們模擬數據集n= 100000, p= 100，其中(β1,β2,··· ,β20) =(50,49,··· ,30),(β21,β22,··· ,βp)=(0,0,··· ,0)，預選參數c=0.1, X=(X1,X2,··· ,Xp).Xj～N(0,Σ),Σ = (σij), σij= 0.5|i-j|,1≤i,j ≤n.當應用分布式MCP 方法時，我們將所有數據劃分為J= 5 組，并選擇鄰接矩陣E.分布式MCP 方法變量選擇的結果呈現為圖1.

圖1 當n=100000, p=100 時，分布式MCP 方法變量選擇的結果

上述實驗結果表明，對于非分布式算法，單機很難處理的數據量較大的問題，分布式MCP 正則化可以對其進行變量選擇.進一步，所提出的方法適合當今分布式處理數據的需求.

4.3 實驗3

在這個實驗中，我們比較了Lasso 和分布式MCP 方法的計算效率.我們分別令n=105,106,107,3×107,5×107和p=10,100，并生成7 個數據集.在模擬實驗中，我們使用LARS 來求解L1正則化.然后我們運行了Lasso 和分布式MCP 方法，并比較了兩種算法的計算時間.結果如表2 所示.當n=3×107,5×107, p=10 時，Lasso 失效，而分布式MCP 依然可以有效運行.

表2 Lasso 和分布式MCP 在7 個數據集上的運行時間結果

實驗表明，當n= 105,106,107, p= 10,100 時，Lasso 和分布式MCP 均可以有效地選擇變量；當n ＞3×107時，Lasso 失效，但是分布式MCP 仍然可以有效運行.

4.4 實驗4(前列腺數據)

在這個實驗中，前列腺病人數據來源于UCI 標準數據庫中前列腺病人的相關數據，數據集大小為97.我們研究病人的抗原水平與腫瘤體積記錄(lcavol)，前列腺重量記錄(lweight)，年齡(age)，良性前列腺增生量(lbph)，精囊浸潤(svi)，包膜穿透記錄(lcp)，Gleason 積分(gleason)以及Gleason4 或5 所占的百分比(pgg45)等8 個指標之間的關系.將8 個指標看做輸入變量，將病人的抗原水平看做輸出變量，在對數據進行預處理后，我們用Lasso 和分布式MCP 算法對數據集進行變量選擇，并觀察實驗結果.

在運行分布式MCP 方法時，本文將所有的樣本數據分配給5 臺計算機，調控參數λ由分布式5 折交叉驗證來選擇，鄰接矩陣為E，預選參數c= 0.1.Lasso 選擇結果如圖2 所示，共有5 個變量選進模型，均方誤差為0.7531.分布式MCP 選擇結果如圖3 所示，共有3 個變量選進模型，均方誤差0.5029.分布式MCP 選擇結果比Lasso 少選擇了兩個變量，這說明分布式MCP 方法可以更有效地解決該類問題，相比L1可以得到更稀疏的解.

圖2 Lasso 選擇結果

圖3 分布式MCP 選擇結果

5 結論

本文研究分布式MCP 方法，依據數據所屬計算機的不同，將MCP 模型轉化為等價的分布式MCP 模型.基于ADMM 算法，我們提出了分布式MCP 算法并證明了它的收斂性.分布式的變量選擇結果與非分布式方法變量選擇結果相同.實驗表明，分布式MCP 方法能夠有效的處理分布式存儲的數據.進一步，本文提出的方法可以推廣到其他分布式非凸方法.所有這些問題都在我們目前的研究之下.