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

沈圓圓, 曹文飛, 韓國棟

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

1 引言

Logistic 回歸是統計分析中基本且十分重要的預測模型，在風險預測、流行病學及生物分析等領域具有廣泛應用[1-3].在Logistic 回歸中，模型的變量選擇與參數估計是解決統計問題的關鍵.考慮一般的Logistic 回歸模型

這里

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

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

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

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

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

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

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

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

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

2 Logistic組稀疏回歸的Bayes建模

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

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

2.1 數據似然

數據yn=1 或yn=-1 依賴于N維輸入xn.假設對數似然比

關于xn是線性的，yn=1 的條件似然由sigmoid 函數給出，即

易見

此時，條件分布密度函數可建模為

2.2 參數w建模

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

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

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

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

其中p(zi)稱為混合分布.選擇一個恰當的混合分布可以保證(8)式中積分能給出顯式解，從而有利于后續參數推斷.因此，混合分布的選擇對參數推斷是至關重要的.

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

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

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

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

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

3 變分推斷

在本節中，我們將利用變分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)及聯合分布

代入上式，可以得到

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

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

考慮sigmoid 函數的一個下界函數

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

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

3.1 估計參數w

根據(14)式，當?k等于w 時，可得

于是

3.2 估計參數z

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

其中

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

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

3.3 估計超參數ai與變分參數ξn

于是，ai的期望為

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

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

3.4 VBLGL算法

綜上，可得

根據上述概率分布形式，通過迭代更新相關參數的矩直到滿足收斂條件，可以得到有關參數的估計.在迭代過程中，全局參數ka設置為10-6，θa設置為10-6，每個λi= 1，最大迭代步數為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)更新每個重復di次，依次更新i=1,2,···,G，直到得到一個

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

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

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

步驟7輸出參數〈w〉以及其他參數的結果.

4 預測密度

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

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

兩邊取對數，得

其中

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

5 模擬實驗

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

5.1 協變量不分組情形

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

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

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

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

5.2 協變量分組情形

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

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

6 總結與展望

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

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