Logistic組稀疏回歸模型的Bayes建模及變分推斷

沈圓圓, 曹文飛, 韓國棟

(陜西師范大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，西安 710119)

1 引言

Logistic 回歸是統(tǒng)計分析中基本且十分重要的預(yù)測模型，在風(fēng)險預(yù)測、流行病學(xué)及生物分析等領(lǐng)域具有廣泛應(yīng)用[1-3].在Logistic 回歸中，模型的變量選擇與參數(shù)估計是解決統(tǒng)計問題的關(guān)鍵.考慮一般的Logistic 回歸模型

這里

表示 sigmoid 函數(shù)，w∈RN表示回歸系數(shù)向量，xn= (xn1,xn2,···,xnN)T是第n個輸入向量，yn={-1,1}是對應(yīng)的二元響應(yīng)變量，M表示觀測樣本的數(shù)目.所有的輸入向量放到一起構(gòu)成M×N設(shè)計矩陣所有的響應(yīng)變量放到一起構(gòu)成響應(yīng)向量 y=(y1,y2,···,yM)T.

Logistic 回歸的一般目標(biāo)是基于一組觀測樣本(xn,yn),n=1,2,···,M，設(shè)計某種準(zhǔn)則估計回歸系數(shù)w.最大似然準(zhǔn)則是參數(shù)估計中經(jīng)典準(zhǔn)則之一.Logistic 回歸模型(1)經(jīng)過最大似然準(zhǔn)則得到如下估計

在實際問題中，該似然估計一般是可行的，但對于高維問題，該估計往往會失效，而且該估計不具有變量選擇的功能.在應(yīng)對高維問題并嵌入變量選擇功能方面，Lasso[4]是一個十分成功的準(zhǔn)則.Logistic 回歸模型(1)經(jīng)過Lasso 準(zhǔn)則得到如下估計

上述Logistic Lasso 估計通過施加?1懲罰誘導(dǎo)出稀疏性，從而實現(xiàn)了變量選擇.與最佳子集選擇、逐步回歸等方法相比，Lasso 估計有巨大的計算優(yōu)勢.但是，Lasso 做變量選擇時僅依賴于單個輸入變量而沒有注意到各個變量之間的關(guān)系，因此常常會選擇出一些不太重要的變量.此外，Lasso 解依賴于因子間的正交方式，即對一個分類預(yù)測子選擇不同的正交方式將產(chǎn)生不同的解，這是不合理的，因為因子選擇和估計問題的解不應(yīng)取決于因子是怎樣編碼的.

Group Lasso[5,6]適當(dāng)擴(kuò)展了Lasso 懲罰，在一定程度上克服了上述問題.Kim 等人[7]首次將Group Lasso 應(yīng)用于Logistic 回歸模型(1)，并設(shè)計出了一個有效的推斷算法.在此基礎(chǔ)上，Meier 等人[8]提出來一個改進(jìn)的方法，得到了Logistic Group Lasso 估計

其中Ii是屬于第i個變量組的指標(biāo)集，i=1,2,···,G，G是組的數(shù)目，為第i組的數(shù)目.顯然，(4)中的懲罰項可視為介于?1和?2之間的懲罰，并且包含?1范數(shù)作為特殊情形.注意到在(4)中，由于稀疏性被施加在組上而不是單個變量上，因而該估計考慮到了信號元素之間的相互關(guān)系，并且該估計在群組正交重新參數(shù)化下不變，從而一定程度上克服了Lasso 的缺點.Logistic Group Lasso 估計(4)在高維統(tǒng)計、壓縮感知[9,10]、機(jī)器學(xué)習(xí)[6,11]以及圖像與網(wǎng)絡(luò)分析[12-14]等領(lǐng)域均顯示出良好性能.

但是，計算Logistic Group Lasso 估計(4)時需要通過交叉驗證方法確定模型參數(shù)λ的最優(yōu)值，該過程無論在數(shù)據(jù)量還是計算上都需要很大的代價.而且基于該準(zhǔn)則的點估計導(dǎo)致對新樣本x0的預(yù)測是一個硬的二元決策，而在理想情況下，我們期望在該新樣本上的預(yù)測是一個條件分布，并通過此條件分布來捕獲預(yù)測中的不確定性.

為了克服上述優(yōu)化方法的弊病，本文提出一個新的Bayes 概率方法來解決Logistic 組稀疏建模與推斷問題.具體來說，首先使用高斯-方差混合公式建模組稀疏參數(shù)向量先驗的分層結(jié)構(gòu)；其次，基于變分Bayes 方法導(dǎo)出所有參數(shù)的后驗推斷公式；最后，基于參數(shù)的后驗公式導(dǎo)出新樣本的概率預(yù)測公式.值得指出的是，用高斯-方差混合公式可以帶來諸多優(yōu)點，如參數(shù)向量的高斯分布建模使得積分計算容易實現(xiàn)、共軛先驗機(jī)制可以給出參數(shù)的閉式估計、以及參數(shù)的自動化估計將省略掉調(diào)制超參數(shù)的復(fù)雜過程等.

本文內(nèi)容組織如下：第2 節(jié)，建立Logistic 組稀疏回歸的概率模型；在第3 節(jié)，基于變分Bayes 方法，推導(dǎo)模型中相關(guān)參數(shù)的后驗推斷公式，并給出有效算法框架；基于回歸參數(shù)的后驗概率公式，在第4 節(jié)中導(dǎo)出新樣本的預(yù)測密度；第5 節(jié)中的多組模擬數(shù)值實驗將驗證所提出方法具有良好的分類性能；最后一節(jié)總結(jié)本文工作并展望后續(xù)的研究問題.

2 Logistic組稀疏回歸的Bayes建模

本節(jié)我們將提出Logistic 組稀疏回歸模型的Bayes 概率版本.首先，給出數(shù)據(jù)似然的分布形式；其次，利用高斯-方差混合分布給出回歸參數(shù)w 組稀疏的概率刻畫；最后，構(gòu)建出所有分布參數(shù)的超先驗分布，從而為整個生成過程建立一個全Bayes 概率模型.

下文中diag(v)為對角陣，v 為對角線元素組成的向量，〈·〉是關(guān)于相應(yīng)分布的期望.

2.1 數(shù)據(jù)似然

數(shù)據(jù)yn=1 或yn=-1 依賴于N維輸入xn.假設(shè)對數(shù)似然比

關(guān)于xn是線性的，yn=1 的條件似然由sigmoid 函數(shù)給出，即

易見

此時，條件分布密度函數(shù)可建模為

2.2 參數(shù)w建模

假定參數(shù) w 可劃分為G組，wi表示第i組 (i=1,2,···,G)，該組有di個參數(shù).特別地，當(dāng)G=N時，每組都只有1 個參數(shù)，即退化為不分組情形.

各組之間的先驗一般假設(shè)是獨立的，這種假設(shè)在很多實際問題中經(jīng)常被采用，例如多頻帶信號重建[15]、微陣列分析[16]等.在此假設(shè)下，回歸參數(shù)w 條件分布的乘積形式為

于是，回歸參數(shù)w 的整體概率建模問題轉(zhuǎn)化為各分組回歸參數(shù)wi的概率建模問題.我們將采用高斯-方差混合公式[17]的分層概率建模機(jī)制來實現(xiàn)回歸參數(shù)w 的組稀疏性質(zhì).

對zi積分，可得wi的邊緣分布

其中p(zi)稱為混合分布.選擇一個恰當(dāng)?shù)幕旌戏植伎梢员ＷC(8)式中積分能給出顯式解，從而有利于后續(xù)參數(shù)推斷.因此，混合分布的選擇對參數(shù)推斷是至關(guān)重要的.

本文中，我們選取伽馬分布

作為混合分布p(zi)的具體形式，其密度函數(shù)為

將此混合分布代入等式(8)，可以得到關(guān)于回歸參數(shù)wi的邊緣分布

該邊緣分布被稱作Mckay’s Bessel 函數(shù)分布[18]，其中Kλi是二階修正Bessel 函數(shù).

至此，我們得到了wi的邊緣分布(9).該分布具有尖峰厚尾形式，因而可以誘導(dǎo)出稀疏解.從(9)式中可以看到，該分布具有超參數(shù)ai與λ.為避免繁瑣的超參數(shù)調(diào)制，接下來將在第3 節(jié)給出ai的先驗建模，本文中對參數(shù)λ我們不對它有封閉形式的更新，只需要它有數(shù)值解，把它看作一個自由參數(shù)，具體值將在第3 節(jié)說明.

3 變分推斷

在本節(jié)中，我們將利用變分Bayes 方法[19,20]，給出后驗概率p(w,z,a|y)的近似推斷.為方便敘述，記為?:={w,z,a}，近似后驗概率分布記為q(?).q(?)與p(?|y)之間的Kullback-Leibler 散度(以下簡稱KL 散度)為

以下，我們將通過最小化KL 散度來求得q(?).

將p(?,y)=p(?|y)P(y)及聯(lián)合分布

代入上式，可以得到

其中L(q)稱為變分下界.對于固定的觀測y，易知上式左邊為常數(shù).因此，KL 散度最小等價于變分下界L(q)最大.進(jìn)一步，假設(shè)后驗參數(shù)變量相互獨立，可得q(w,z,a) =q(w)q(z)q(a).此時，變分下界為

由于數(shù)據(jù)似然(6)在指數(shù)族上不存在共軛先驗，導(dǎo)致上式中變分下界L(q)關(guān)于后驗概率分布q(w)、q(z)及q(a)的極大值點沒有顯式解.因此，需要對數(shù)據(jù)似然做進(jìn)一步的處理.

考慮sigmoid 函數(shù)的一個下界函數(shù)

值得注意的是，每個數(shù)據(jù)均有一個局部的變化參數(shù)ξn.用h(w,ξ)代替數(shù)據(jù)似然P(y|X,w)，代入(11)，可得新的變分下界

其中??k表示在?中移除掉?k.注意到，上式左邊參數(shù)的后驗分布是通過右邊參數(shù)的分布表達(dá)出來，而右邊參數(shù)的分布又是未知的，因此需要通過一個迭代過程才能給出最終解.以下，我們將根據(jù)(14)式，計算出當(dāng)?k分別等于w、z 及a 時的具體分布；最后，通過迭代給出詳細(xì)的算法框架.

3.1 估計參數(shù)w

根據(jù)(14)式，當(dāng)?k等于w 時，可得

于是

3.2 估計參數(shù)z

此時僅需對每個q(zi)進(jìn)行估計.根據(jù)(14)式，可得zi的變分后驗估計為

其中

由此分布形式可知q(zi)服從廣義逆高斯(GIG)分布，其期望為

其中Σwi表示Σw中對應(yīng)于第i組的子矩陣.的后驗估計可通過該分布的矩[21]給出

3.3 估計超參數(shù)ai與變分參數(shù)ξn

于是，ai的期望為

其中〈zi〉用(18)式計算.

對于局部參數(shù){ξn}，由Bishop[19]可知，可以通過最大化邊緣似然下界來確定，具體的更新表達(dá)式為

3.4 VBLGL算法

綜上，可得

根據(jù)上述概率分布形式，通過迭代更新相關(guān)參數(shù)的矩直到滿足收斂條件，可以得到有關(guān)參數(shù)的估計.在迭代過程中，全局參數(shù)ka設(shè)置為10-6，θa設(shè)置為10-6，每個λi= 1，最大迭代步數(shù)為10000 步，收斂閾值tol 為10-7.為敘述方便，簡稱我們算法為VBLGL(Variational Bayesian Logistic Group Lasso)，整個算法迭代過程見算法1.

算法1VBLGL 算法

步驟1初始化：ξ ∈RM,ka=10-6,θa=10-6, tol=10-7；

步驟2更新〈w〉：將Z-1和ξ值代入(16)更新〈w〉；

步驟3更新zi：用(17)更新每個重復(fù)di次，依次更新i=1,2,···,G，直到得到一個

步驟4更新ai：用(18)更新〈zi〉，并代入(19)更新〈ai〉.同上，依次更新G次得到一個a；

步驟5更新ξ：用(20)更新ξn，依次對n= 1,2,···,M進(jìn)行更新，直到得到一個ξ；

步驟6重復(fù)步驟2 到步驟5，直到前后兩次〈w〉差的?2范數(shù)小于tol 時，迭代終止；

步驟7輸出參數(shù)〈w〉以及其他參數(shù)的結(jié)果.

4 預(yù)測密度

對于新樣本xn，用sigmoid 函數(shù)的下界代替數(shù)據(jù)似然P(y|X,w)，同時用回歸參數(shù)w 的近似后驗分布q(w)代替p(w|X,y)，可得預(yù)測密度為

將q(w)代入上式并注意到高斯密度函數(shù)可以完全積分，可得

兩邊取對數(shù)，得

其中

在(21)中，最大化參數(shù)ξ，即可得該參數(shù)的更新公式

5 模擬實驗

本節(jié)中，將通過模擬實驗驗證本文所提出的方法是可行的，并具有良好的預(yù)測效果.模擬實驗分成兩組：第一組對應(yīng)協(xié)變量不分組情形，此時所有的方法都退化為解決Logistic 稀疏回歸問題；第二組對應(yīng)協(xié)變量分組情形，該組實驗考慮Logistic 組稀疏回歸問題.在第一組實驗里，我們比較了Meier 等人[8]針對協(xié)變量分組的Logistic 組稀疏回歸問題提出的Logistic group Lasso 方法(簡記為LGL 方法)，以及Drugowitsch[22]針對協(xié)變量不分組的Logistic 稀疏回歸問題提出的自動關(guān)聯(lián)決策的變分Bayes 方法(簡記為VBARD 方法).在第二組實驗里，我們比較了LGL 方法.在這兩組實驗里，我們的方法都簡記為VBLGL 方法.

5.1 協(xié)變量不分組情形

該部分的實驗在兩組不同大小的數(shù)據(jù)集上進(jìn)行.第一組數(shù)據(jù)集中訓(xùn)練集的產(chǎn)生過程如下：首先，通過Matlab 中的命令rand(300,500)產(chǎn)生300×500 的訓(xùn)練集設(shè)計矩陣Xtr，該矩陣中的元素服從(0,1)區(qū)間上的均勻分布；其次，產(chǎn)生大小為500×1 的回歸參數(shù)向量wtr，其非零元素有50 個，每個元素均服從標(biāo)準(zhǔn)正態(tài)分布；最后，通過訓(xùn)練集設(shè)計矩陣Xtr與回歸參數(shù)向量wtr生成反應(yīng)向量

不難看出，反應(yīng)向量由取值為-1 或者1 的300 個輸出樣本構(gòu)成.通過類似的方式，我們可生成樣本大小為1000 的測試集，具體數(shù)據(jù)記為Xte，wte以及yte.第二組數(shù)據(jù)集的產(chǎn)生過程類似第一組，但訓(xùn)練集和測試集中樣本的大小都為1000，維數(shù)仍為500.

表1 展示出三種方法在這兩組數(shù)據(jù)集上10 次模擬的平均分類結(jié)果，其中mean(tr)表示在訓(xùn)練集上預(yù)測誤差均值，mean(te)表示在測試集上預(yù)測誤差均值，而std(tr)與std(te)分別表示在訓(xùn)練集和測試集上預(yù)測誤差的標(biāo)準(zhǔn)差.如表1 所示，通過對比指標(biāo)mean(tr)可以看到，所提出的方法VBLGL 取得了最小的平均訓(xùn)練誤差.這表明，與VBARD 及LGL方法相比，VBLGL 方法具有更好的擬合性能.指標(biāo)mean(te)的對比結(jié)果進(jìn)一步表明，在Logistic 非分組稀疏回歸問題上，VBLGL 方法相比于其他兩種方法具有競爭性的表現(xiàn).另外，從標(biāo)準(zhǔn)差std(tr)與std(te)指標(biāo)可以看出，三種方法取得的分類結(jié)果都是比較穩(wěn)定的.

表1: 協(xié)變量不分組情形下分類效果對比

5.2 協(xié)變量分組情形

在這部分實驗中，考慮協(xié)變量分組情形下的分類效果.為論述方便，僅考慮等分組情形.對于組大小不相等的情形，本文所提出的方法仍然適用.協(xié)變量的維數(shù)固定為500，每10 個相鄰元素為一組，總共分為50 組.非零元素組的數(shù)目記為p(亦稱為“組稀疏度”).在接下來的實驗中，固定組稀疏度參數(shù)p= 3.首先，產(chǎn)生回歸系數(shù)w∈R500×1，其中在w 中有p= 3 組值非零，其他組元素值全部為零；其次，固定回歸系數(shù)w 的值，按照上述5.1 節(jié)類似的方式分別產(chǎn)生多組不同樣本大小的訓(xùn)練集設(shè)計矩陣Xtr與反應(yīng)向量ytr，其樣本大小分別為200, 300, 400, 500, 750 和1000；最后，基于回歸系數(shù)w，我們生成測試樣本大小為1000 的測試集設(shè)計矩陣Xte∈R1000×500和反應(yīng)向量yte∈R1000×1.通過上述過程，我們得到6 組訓(xùn)練集和測試集.

表2 展示出兩種方法，在這6 組數(shù)據(jù)集上10 次模擬的平均分類結(jié)果，其中mean(tr)，mean(te)，std(tr)，std(te)的含義如5.1 部分一致.如表2 所示，從指標(biāo)mean(tr)與mean(te)上看，在數(shù)據(jù)集(1)-(6)上，VBLGL 方法在訓(xùn)練集、測試集上的預(yù)測誤差均值均小于LGL 方法.這一現(xiàn)象表明，VBLGL 方法不僅具有更好的擬合精度，而且有更好的分類預(yù)測效果.另外，從標(biāo)準(zhǔn)差指標(biāo)std(tr)、std(te)上來看，這兩種方法所取得的分類結(jié)果都比較穩(wěn)定.總體上來看，VBLGL 方法比LGL 方法更加穩(wěn)定.這一現(xiàn)象是合理的，因為VBLGL 方法的概率化建模和參數(shù)的自動學(xué)習(xí)可以提高算法的穩(wěn)定性.

6 總結(jié)與展望

針對Logistic 組稀疏回歸問題，本文在Bayes 概率框架下建立了參數(shù)估計的新模型，并基于變分Bayes 方法導(dǎo)出參數(shù)的后驗推斷公式和新樣本的概率預(yù)測公式.與建模Logistic 組稀疏回歸問題的優(yōu)化方法進(jìn)行比較，本文提出的方法具有模型可解釋性、參數(shù)的自動估計以及估計量擁有概率解等良好性質(zhì)，而且模擬實驗結(jié)果進(jìn)一步表明所提出的方法具有更好的分類預(yù)測效果.在未來的工作中，我們擬將本文所提出的方法推廣到多分類的Logistic 組稀疏回歸問題上.

表2: 協(xié)變量分組情形下分類效果對比