王格華, 王璞玉, 張 海

(西北大學(xué)數(shù)學(xué)學(xué)院，西安 710069)

1 引言

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

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

其中λ ≥0, γ ＞1.

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

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

2 分布式MCP 模型

模型(3)將所有數(shù)據(jù)拆為J個(gè)部分，其與原模型(1)等價(jià).

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

3 分布式MCP 算法

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

3.1 ADMM 算法

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

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

為了方便表述，我們將模型(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)的二次增廣拉格朗日函數(shù)為

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

ADMM 算法分為以下三個(gè)步驟：

步驟1 更新Z

步驟2 更新β

步驟3 更新拉格朗日乘子

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

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

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

故(11)成立.

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

引理1 對(duì)于所有的k，下面的結(jié)論成立：

且

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

(A4)可由下式得到

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

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

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

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

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

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

下面，我們證明增廣拉格朗日函數(shù)對(duì)任意k都有界.

對(duì)于任意一個(gè)Zk+1，存在β′，使得AZk+1+Bβ′=0，那么

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

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

證明 首先，有

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

成立.同時(shí)，我們有

又因?yàn)?/p>

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

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

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

3.2 參數(shù)估計(jì)

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

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

根據(jù)數(shù)據(jù)所屬計(jì)算機(jī)的不同，我們將(8)分解為J個(gè)子問題

下面，我們給出坐標(biāo)下降法的收斂性分析.

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

證明 假設(shè)預(yù)選參數(shù)c ＞1/(γ*J)，則對(duì)于每個(gè)計(jì)算機(jī)j，和所有維度m, O是關(guān)于[Zj]m的嚴(yán)格凸函數(shù).O的嚴(yán)格凸性和連續(xù)性滿足文獻(xiàn)[26]中定理4.1 和定理5.1 的條件，進(jìn)而坐標(biāo)下降算法收斂到坐標(biāo)最小值.由于O的所有方向?qū)?shù)都存在，因此每個(gè)坐標(biāo)最小值也是局部最小值.

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

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

同樣，通過使用坐標(biāo)下降法，給出局部βj的第m個(gè)坐標(biāo)的顯式解

注1 如果計(jì)算機(jī)i與計(jì)算機(jī)j不連接，則Eij=0.這意味著我們?cè)诟耑j和βj時(shí)，只需要計(jì)算機(jī)j的鄰居的信息.

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

算法1 分布式MCP 算法

對(duì) ?j ∈J, 令β =(1,··· ,1)T, Z =(1,··· ,1)T, Vij =(1,··· ,1)T,局部運(yùn)行for k =1,2,···第j 臺(tái)計(jì)算機(jī)在Nj 與鄰居傳遞信息對(duì)于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 折交叉驗(yàn)證

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

重復(fù)上述步驟，我們可以計(jì)算每個(gè)λi的均方誤差MSE，使得

最終，我們選擇預(yù)測(cè)誤差最小的λdcv=arg minλi{e(λi)}作為調(diào)控參數(shù).

4 實(shí)驗(yàn)

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

4.1 實(shí)驗(yàn)1(模擬實(shí)驗(yàn))

在這個(gè)實(shí)驗(yàn)里，對(duì)于我們考慮的線性模型

我們?cè)O(shè)置變量數(shù)量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.當(dāng)應(yīng)用分布式MCP 方法時(shí)，我們將所有數(shù)據(jù)劃分為J=5 組，并選擇鄰接矩陣

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

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

表1 Lasso 和分布式MCP 在實(shí)驗(yàn)1 中的結(jié)果

4.2 實(shí)驗(yàn)2

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

圖1 當(dāng)n=100000, p=100 時(shí)，分布式MCP 方法變量選擇的結(jié)果

上述實(shí)驗(yàn)結(jié)果表明，對(duì)于非分布式算法，單機(jī)很難處理的數(shù)據(jù)量較大的問題，分布式MCP 正則化可以對(duì)其進(jìn)行變量選擇.進(jìn)一步，所提出的方法適合當(dāng)今分布式處理數(shù)據(jù)的需求.

4.3 實(shí)驗(yàn)3

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

表2 Lasso 和分布式MCP 在7 個(gè)數(shù)據(jù)集上的運(yùn)行時(shí)間結(jié)果

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

4.4 實(shí)驗(yàn)4(前列腺數(shù)據(jù))

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

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

圖2 Lasso 選擇結(jié)果

圖3 分布式MCP 選擇結(jié)果

5 結(jié)論

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