陳少東 李志強(qiáng)

(北京化工大學(xué) 數(shù)理學(xué)院， 北京 100029)

引 言

廣義線性模型由Nelder等[1]在研究指數(shù)族分布時引入，該模型不僅能刻畫連續(xù)型數(shù)據(jù)，還可刻畫計數(shù)數(shù)據(jù)、屬性數(shù)據(jù)等離散型數(shù)據(jù)，因此在醫(yī)學(xué)、經(jīng)濟(jì)學(xué)以及社會科學(xué)等領(lǐng)域得到廣泛應(yīng)用。

隨著數(shù)據(jù)時代的到來，數(shù)據(jù)呈現(xiàn)海量化趨勢。數(shù)據(jù)海量化主要體現(xiàn)在兩個方面：一是能獲取的變量維度越來越高；二是可獲取的觀測樣本越來越多。在數(shù)據(jù)量不是很大時，面對高維的變量如何進(jìn)行變量選擇的問題已得到了廣泛的研究。在以往文獻(xiàn)中，有許多研究人員利用懲罰函數(shù)進(jìn)行變量選取，例如Tibshirani[2]提出的least absolute shrinkage and selection operator (Lasso)懲罰、Fan等[3]提出的smoothly clipped absolute deviation (SCAD)懲罰以及Zhang[4]提出的minimax concave penalty (MCP)懲罰等。

為了解決基于懲罰函數(shù)的廣義線性模型變量選擇問題，Breheny和Huang[5]給出一種基于坐標(biāo)下降法的計算方法，該算法能有效求出模型局部最優(yōu)解；然而，在海量數(shù)據(jù)環(huán)境下，由于該方法使用了全部海量數(shù)據(jù)進(jìn)行迭代計算，導(dǎo)致出現(xiàn)計算效率較低甚至計算機(jī)存儲空間不足等問題。為了克服這些瓶頸，需要探究新的估計算法來解決海量數(shù)據(jù)下廣義線性模型的變量選擇問題。

分治算法是在海量數(shù)據(jù)下估計統(tǒng)計模型參數(shù)的一類主流方法。目前已有部分文獻(xiàn)利用分治算法研究海量數(shù)據(jù)下的統(tǒng)計模型。其中，Lin等[6]提出了一種解決海量數(shù)據(jù)非線性估計方程的估計方法；在變量選取方面，方方等[7]研究了海量數(shù)據(jù)下Logistic模型的平均算法；Chen等[8]給出了一種非凸懲罰下具有典則連接的廣義線性模型的變量選擇估計方法。基于分治算法處理海量數(shù)據(jù)的詳細(xì)介紹可參考文獻(xiàn)[9]。

然而，以上研究的相關(guān)工作還有待進(jìn)一步深入，例如采用文獻(xiàn)[7]中的模型平均估計方法難以進(jìn)行高效的運(yùn)算；文獻(xiàn)[6]、文獻(xiàn)[8]只研究了海量數(shù)據(jù)下具有典則連接的廣義線性模型。但在實際建模應(yīng)用中，大部分連接函數(shù)都是非典則連接，因此研究一般廣義線性模型的懲罰估計具有更重要的現(xiàn)實意義。

本文基于文獻(xiàn)[5]在Logistic模型下的估計方法，推導(dǎo)出一般廣義線性模型在SCAD懲罰和MCP懲罰下的迭代公式，并將該方法推廣到了海量數(shù)據(jù)下一般廣義線性模型中。該方法避免了在迭代過程中一次性加載全部數(shù)據(jù)，從而有效地克服了基于原始海量數(shù)據(jù)的估計方法帶來的計算機(jī)存儲空間不足等問題。與此同時，考慮到當(dāng)前在解決海量數(shù)據(jù)問題時，大部分處理工具采用分布式并行框架，本文結(jié)合分治算法與分布式計算的共同特征，給出了算法在Spark集群環(huán)境的計算步驟。模擬結(jié)果表明，估計算法在克服內(nèi)存瓶頸的同時提高了計算效率，其估計精度要優(yōu)于隨機(jī)抽樣得到的估計，而且能有效應(yīng)用于集群環(huán)境。

1 廣義線性模型及懲罰方法

1.1 廣義線性模型

一般廣義線性模型滿足

(1)

式中，f(·)為概率密度函數(shù),yi∈為響應(yīng)變量，xi∈p為觀測變量，β∈p為未知參數(shù)，θ為典則參數(shù)，a(·)、b(·)、c(·)為已知函數(shù)，μi為均值函數(shù)，V(μi)為方差函數(shù)，g為連接函數(shù)，為常數(shù)，i=1,…,N。

為了估計模型的未知參數(shù)β，通常采用極大似然估計方法，其似然函數(shù)可表示為

(2)

1.2 基于懲罰的變量選擇

懲罰函數(shù)方法是一種常用的變量選擇方法，在進(jìn)行變量選擇時，通過在式(2)上添加懲罰項，即可得到懲罰目標(biāo)函數(shù)，現(xiàn)將其轉(zhuǎn)化為極小化問題

(3)

式中，βj為β的第j維分量。

SCAD懲罰和MCP懲罰是兩種常用的非凸懲罰函數(shù)，與Lasso方法相比，這兩種懲罰能保證Oracle性質(zhì)，因此本文將圍繞這兩種懲罰方法進(jìn)行研究。

1.2.1SCAD懲罰函數(shù)

SCAD懲罰函數(shù)在βj∈[0,∞)上定義為

(4)

式中，λ≥0，γ>2。

當(dāng)βj較小時，SCAD和Lasso具有一樣的懲罰力度；隨著βj增大至λ，SCAD懲罰力度開始逐漸下降；當(dāng)βj>γλ時，懲罰力度降為0。SCAD懲罰保留了Lasso懲罰的稀疏性等優(yōu)點，并修正了Lasso懲罰過度壓縮系數(shù)的問題。

1.2.2MCP懲罰函數(shù)

MCP懲罰函數(shù)在βj∈[0,∞)上定義為

(5)

式中，λ≥0，γ>1。

MCP懲罰起初與Lasso懲罰有相等的懲罰力度，然后隨著βj的增大，懲罰力度逐漸減小，當(dāng)βj>γλ時，懲罰力度降到0。在高維線性回歸中，MCP懲罰方法是一種幾乎無偏且準(zhǔn)確的變量選擇方法。

1.2.3Logistic模型懲罰估計方法

為了解決Logistic模型在上述兩種懲罰下的估計問題，Breheny和Huang[5]結(jié)合了文獻(xiàn)[10]給出的用于求解廣義線性模型的迭代加權(quán)最小二乘法，將優(yōu)化問題轉(zhuǎn)化為線性模型下加權(quán)最小二乘的優(yōu)化問題，然后采用以下方法進(jìn)行估計。

設(shè)對數(shù)似然函數(shù)的第m次迭代的估計值為m，Logistic模型在β=m處的近似懲罰似然函數(shù)為

(6)

記Xj為第j個自變量所對應(yīng)的整體觀測值，mj為m的第j維分量，為了求解式(6)的極小值，文獻(xiàn)[5]給出具體估計步驟如下：

(4)基于坐標(biāo)下降法重復(fù)步驟(1)～(3)直到收斂。

當(dāng)懲罰函數(shù)為SCAD時，f(zj,λ,γ)滿足

當(dāng)懲罰函數(shù)為MCP時，f(zj,λ,γ)滿足

這里，

2 海量數(shù)據(jù)下一般廣義線性模型的懲罰估計方法

2.1 一般廣義線性模型下的懲罰估計

Logistic模型僅僅是廣義線性模型中的一個特例，在實際問題中，因變量除了服從兩點分布外，還可能服從多項分布、泊松分布等，且在建模時，大部分連接函數(shù)都是非典則連接。因此有必要進(jìn)一步研究一般情形下的廣義線性模型的估計方法。

(7)

因此結(jié)合式(6)的懲罰函數(shù)和式(7)，可以利用與式(6)類似的計算步驟得到一般廣義線性模型的懲罰估計。

上述方法在求解過程中需要把所有的觀測數(shù)據(jù)一次性加載到計算機(jī)內(nèi)存當(dāng)中，當(dāng)觀測數(shù)據(jù)量N很大時，可能會帶來計算機(jī)內(nèi)存不足等問題。因此，需要對該算法進(jìn)行改進(jìn)，以適應(yīng)海量數(shù)據(jù)情形。

2.2 海量數(shù)據(jù)下一般廣義線性模型的懲罰估計

為了解決海量數(shù)據(jù)下具有典則連接的廣義線性模型的變量選擇問題，Chen等[8]采用分治思想，給出了一種具有典則連接的廣義線性模型的變量選擇方法,其基本思想是：首先將全部數(shù)據(jù)劃分為K個互不相交的子數(shù)據(jù)集；接著分別求解每塊子數(shù)據(jù)集下的變量選擇及參數(shù)估計結(jié)果；然后通過投票的方式解決不同子數(shù)據(jù)集下選取的變量可能不同的問題，確定最終選取的變量；最后通過加權(quán)方式對K個估計結(jié)果進(jìn)行聚合，將其作為最終估計結(jié)果。本文延續(xù)上述思路，給出了具有一般連接函數(shù)的廣義線性模型變量選擇的估計方法。

為了求解第k塊數(shù)據(jù)集Dk下的壓縮統(tǒng)計量k，記Xkj為Dk下第j個自變量所對應(yīng)的觀測值，mk為第m次迭代的估計值，mkj為mk的第j維分量，通過1.2.3節(jié)方法進(jìn)行估計：

(4)基于坐標(biāo)下降法重復(fù)步驟(1)～(3)直到收斂。這里，k、mki、kj均為β=mk時代入相應(yīng)公式計算得到的值。

因為不同的子數(shù)據(jù)集選取的變量可能不同，所以可先采用投票方式確定選取的變量

(8)

當(dāng)vj=1表明第j個變量被選中，否則剔除該變量,I(·)為示性函數(shù)，w∈(0,K]為閾值。設(shè)A={j|vj=1,j=1,…,p}為被選中變量集。令A(yù)k表示由k的第j維分量組成的子向量，其中j∈A。

為了對K個數(shù)據(jù)集上的估計結(jié)果Ak進(jìn)行聚合，將似然函數(shù)的二階近似導(dǎo)數(shù)作為權(quán)重矩陣進(jìn)行加權(quán)，得到最終的聚合估計結(jié)果

(9)

式中SAk=ATSkA。令v=diag(v1,…,vp)，則A為由v的第{j|j∈A}列元素組成的子矩陣。

由式(9)可知，在變量選擇及參數(shù)估計的實現(xiàn)過程中，僅需保留每個子數(shù)據(jù)集Dk上的壓縮變量Ak、SAk即可得到整塊數(shù)據(jù)集的近似估計。因為每個小數(shù)據(jù)集上的估計是獨(dú)立的，所以對計算機(jī)性能的依賴大大降低，只需要能保證每個小數(shù)據(jù)集的估計能正常進(jìn)行即可。

2.3 基于Spark集群的并行實現(xiàn)

當(dāng)前處理海量數(shù)據(jù)的主流工具是基于分布式搭建的。分布式計算通過將一個大型計算任務(wù)拆分為多個子任務(wù)，然后將這些子計算任務(wù)分發(fā)給集群中的多個計算機(jī)節(jié)點進(jìn)行計算，以此來完成并行化計算。

在本文所提算法的估計過程中，數(shù)據(jù)壓縮是耗時最久、占用計算資源最多的步驟，然而基于分治算法，每個子數(shù)據(jù)集在進(jìn)行數(shù)據(jù)壓縮過程中是獨(dú)立的，該特性使其適合以分布式方式實現(xiàn)。

Spark和Hadoop是目前處理海量數(shù)據(jù)的主流工具，借助Hadoop的存儲框架HDFS對數(shù)據(jù)按塊進(jìn)行分布式存儲，并采用適合處理迭代式計算的Spark計算框架，便能將該算法以并行的方式實現(xiàn)，從而提高計算效率。

并行計算流程如圖1所示，每個計算節(jié)點從HDFS中讀取不同的數(shù)據(jù)塊，然后以并行的方式對每一個子數(shù)據(jù)塊上數(shù)據(jù)進(jìn)行壓縮，接著將所有數(shù)據(jù)塊上的結(jié)果廣播給主節(jié)點，讓主節(jié)點完成聚合步驟，最后輸出結(jié)果。

3 數(shù)值模擬

本節(jié)對算法進(jìn)行數(shù)值模擬，以驗證算法的有效性。有效性通過估計精度和時間消耗兩個角度來衡量，模擬通過普通單機(jī)和Spark集群兩種方式完成。選擇常用的兩點分布，即yi|xi～B(1,μi)，其中μi通過以下兩種方式來建模。

(1)在典則連接下，即為常用的Logistic模型

(2)在非典則連接下，選擇Probit模型

式中，Φ(·)為標(biāo)準(zhǔn)正態(tài)分布的累積函數(shù)。

模擬時選取的數(shù)據(jù)量N=106，β0的維數(shù)為200，其中非零個數(shù)為24，自變量獨(dú)立同分布，服從標(biāo)準(zhǔn)正態(tài)分布。

3.1 單機(jī)計算

本次實驗在Ubuntu 16.0.4操作系統(tǒng)下，使用Python 3.7.0完成，計算機(jī)配置為CPU 3.40 GHz、內(nèi)存8 G。所有實驗重復(fù)50次取平均值作為最終結(jié)果。

為了檢驗估計效果，重新隨機(jī)生成了100萬個樣本點，并將這些樣本下的模型分類準(zhǔn)確率作為估計效果的評價指標(biāo)，結(jié)果如表1所示。在所有的實驗中，變量選取結(jié)果均是正確的非零變量。由表1可以看出，Logistic模型的分類準(zhǔn)確率可達(dá)0.928以上，Probit模型的分類準(zhǔn)確率可達(dá)0.958以上，并且經(jīng)過分塊后，兩個模型的分類準(zhǔn)確率下降基本控制在0.1%內(nèi)，這說明該估計的精確度較高且穩(wěn)定。

表1 分類準(zhǔn)確率與分塊數(shù)的關(guān)系

a—Logistic；b—Probit。

表2給出了單機(jī)環(huán)境下參數(shù)估計消耗的時間與分塊數(shù)的關(guān)系。由表2可知，隨著分塊數(shù)的增加，兩個模型的參數(shù)估計過程所消耗的時間逐漸減少，這說明了分塊方法在一定程度上可以提高模型參數(shù)的估計效率。

表2 單機(jī)環(huán)境下計算時間與分塊數(shù)關(guān)系

a—Logistic；b—Probit。

結(jié)合表1和表2的實驗結(jié)果，分塊方法能夠在保證分類準(zhǔn)確率微小下降的情況下提升模型的參數(shù)估計效率。

表3和表4給出了本文的聚合估計方法與隨機(jī)抽樣方法的估計結(jié)果對比。選取的評價指標(biāo)為估計偏差，其計算公式為：設(shè)為估計結(jié)果，β0為真實值，則估計偏差定義為這里‖·‖1為一范數(shù)。因模擬所用數(shù)據(jù)是隨機(jī)產(chǎn)生的，所以可將分塊后的每一塊數(shù)據(jù)集下的估計結(jié)果視為一次以的比例進(jìn)行隨機(jī)抽樣得到的估計結(jié)果。由表中結(jié)果可知，兩個模型的結(jié)果是一致的：固定分塊數(shù)K的情況下，因為聚合方法使用了全部數(shù)據(jù)進(jìn)行估計，所以其估計偏差要比基于隨機(jī)抽樣得到的估計偏差小。該結(jié)論表明：在計算機(jī)內(nèi)存有限的情況下，采用聚合的估計方法總要優(yōu)于隨機(jī)抽樣方法。

表3 Logisitc下聚合方法與抽樣方法偏差比較

a—SCAD；b—MCP。

表4 Probit下聚合方法與抽樣方法偏差比較

a—SCAD；b—MCP。

3.2 Spark集群計算

在海量數(shù)據(jù)環(huán)境下，實際上數(shù)據(jù)多數(shù)是以分布式的方式存儲在多臺計算機(jī)上。Spark是目前用于海量數(shù)據(jù)計算的主流分布式計算框架之一，本文通過Spark計算框架來完成模擬。實驗在Ubuntu 16.0.4系統(tǒng)下，通過4臺普通計算機(jī)組成的Spark集群來完成，計算環(huán)境為Python 3.7.0及Spark 2.3.2。其中，計算節(jié)點的配置為CPU 3.40 GHz、內(nèi)存4 G。

表5給出的是Spark集群下，模型參數(shù)估計消耗的時間與分塊數(shù)的關(guān)系。由表5可知，相比于單機(jī)環(huán)境，采用集群方法可以提高模型參數(shù)估計的效率，這也說明了分塊方法在集群環(huán)境中的可行性。

4 實證分析

選取搜狗實驗室(https:∥www.sogou.com/labs)公開的2012年6月—7月科技板塊和股票板塊的新聞數(shù)據(jù)集建立Probit模型，以此來檢驗本文算法在實際應(yīng)用中的可行性。其中隨機(jī)抽取40 138條科技新聞和39 971條股票新聞用于估計模型參數(shù)，將剩余的9 844條科技新聞和4 437條股票新聞用于測試，將其分類準(zhǔn)確度作為模型分類效果評價指標(biāo)。

表5 Spark集群下計算時間與分塊數(shù)的關(guān)系

a—Logistic；b—Probit。

整個建模流程為：首先對新聞文本進(jìn)行分詞并過濾掉無用的停用詞，然后通過term frequency- inverse document frequency (TF- IDF)方法將每個新聞文本轉(zhuǎn)化為1 000維的詞向量，最后通過本文算法完成Probit模型的變量選擇。結(jié)果表明，采用MCP懲罰在測試集上得到的準(zhǔn)確率為0.955 8，采用SCAD懲罰在測試集上得到的準(zhǔn)確率為0.951 4。由以上結(jié)果可知，兩種懲罰方法在實證數(shù)據(jù)集上都有較好的表現(xiàn)，對測試集的新聞分類準(zhǔn)確度可達(dá)0.95以上，這也進(jìn)一步證明了本文給出的方法對于解決實際問題是可行的。

5 結(jié)束語

本文研究了海量數(shù)據(jù)下一般廣義線性模型的變量選擇估計問題，結(jié)合分治思想與坐標(biāo)下降法，給出了基于MCP懲罰和SCAD懲罰的估計方法，該方法有效降低了估計時對計算機(jī)內(nèi)存的依賴，克服了計算時內(nèi)存不足的瓶頸。數(shù)值模擬表明該方法能進(jìn)一步提高計算效率，并通過實證分析證明了本文估計方法在解決實際問題時的可行性。