基于自適應懲罰的潛變量高斯圖模型結構學習

鄭倩貞, 徐平峰

(1. 長春工業大學 數學與統計學院, 長春 130012; 2. 東北師范大學 前沿交叉研究院, 長春 130024)

圖模型的結構學習問題, 是指從給定的數據集中估計出反映隨機變量間獨立性結構的圖, 也稱為模型選擇問題[1]. 高斯圖模型的參數估計和模型選擇問題等價于協方差逆陣的估計問題及協方差逆陣中零元素的識別問題, 協方差逆陣的元素為零表示對應兩個隨機變量在給定其余隨機變量時具有條件獨立性[2]. 隨著高維復雜數據的不斷增多, 研究稀疏圖模型的估計問題尤為重要. 為估計稀疏協方差逆陣, 常用的方法為懲罰似然方法, 例如: 文獻[3-5]基于懲罰似然方法研究了變量完全觀測時的高斯圖模型選擇問題. 但在實際生活中, 觀測變量通常與一些潛在的不可觀測變量(即潛變量)相關. 此時, 觀測變量間的圖模型結構不一定具有稀疏性, 需考慮潛變量對觀測變量的影響, 在給定潛變量時探討觀測變量間的條件獨立性. 通常情況下, 潛變量會對圖模型選擇問題造成很大困難, 因為潛變量的個數以及潛變量與觀測變量之間的關系可能都是未知的.

Chandrasekaran等[6]首先對含潛變量的圖模型選擇問題進行了研究, 提出了當潛變量和觀測變量服從聯合高斯分布時, 觀測變量邊緣協方差陣的逆陣可分解為稀疏陣和低秩陣之和, 其中稀疏陣對應給定潛變量時觀測變量間的條件獨立性. 基于稀疏低秩分解, 文獻[6]提出了正則化極大似然分解框架進行潛變量圖模型的結構學習, 其中正則項包含對稀疏陣施加的LASSO(least absolute shrinkage and selection operator)懲罰以及對低秩陣施加的核范數懲罰. 由于這兩種懲罰均不含自適應權重, 因此稱這兩種懲罰為非自適應懲罰, 由此得到的稀疏陣和低秩陣估計稱為非自適應懲罰似然估計. 文獻[6]將優化問題視為對數-行列式半正定規劃問題, 采用Wang等[7]提出的對數-行列式近端點算法(log-determinant proximal point algorithm, LogdetPPA)進行求解. 針對文獻[6]中的優化問題, Ma等[8]提出了兩種交替方向的求解方法: 一種是將目標問題視為一致性優化問題, 利用傳統的交替方向乘子法(alternating direction method of multipliers, ADMM)[9]進行求解; 另一種是基于近端梯度的交替方向法(proximal gradient-based alternating direction method, PGADM). 這兩種方法都是把包含3個變量的原始問題轉換成包含2個變量塊的優化問題. PGADM算法速度比LogdetPPA更快, 并具有全局收斂性. 此外, Ma[10]提出了一種交替近端梯度方法求解文獻[6]中的優化問題, 并證明了算法的全局收斂性; Meng等[11]對文獻[6]的方法做了進一步的理論研究; 文獻[12-15]等也對含潛變量圖模型結構學習進行了相關研究.

研究表明, 非自適應懲罰得到的估計, 包括LASSO懲罰得到的稀疏協方差逆陣估計[16]以及核范數懲罰在降秩回歸問題中得到的低秩系數陣估計[17], 均存在偏差較大的問題, 而自適應懲罰通常可降低估計偏差. 因此, 本文采用自適應懲罰似然方法處理含潛變量的高斯圖模型結構學習問題, 對稀疏陣和低秩陣進行估計, 以得到給定潛變量時觀測變量間的條件獨立關系. 本文對稀疏陣部分施加自適應LASSO懲罰[18], 對低秩陣部分施加自適應核范數懲罰[17]. 與非自適應懲罰的優化問題類似, 采用ADMM算法優化求解自適應懲罰似然的最小化問題, 并且在求解過程中仍具有顯式表達式, 以確保算法的計算效率. 本文模擬比較了自適應懲罰與非自適應懲罰在潛變量高斯圖模型結構學習和參數估計上的性能. 結果表明, 自適應懲罰顯著優于非自適應懲罰, 有效降低了稀疏陣和低秩陣的估計偏差, 能更準確地學習觀測變量間的條件獨立關系.

1 含潛變量的無向高斯圖模型及其自適應懲罰似然

無向高斯圖模型是指與無向圖G=(V,E)相關的多元正態分布模型, 其中頂點集V中的每個頂點v∈V表示一個高斯隨機變量Xv, 邊集E中的每條邊(u,v)∈E蘊含了給定其余隨機變量Xv′(v′∈V{u,v})時隨機變量Xu與Xv之間的條件相關性.當采用無向高斯圖模型對實際問題進行分析時, 可能會出現隨機變量為潛變量的情況, 即變量不存在觀測值.因此, 本文考慮含潛變量的無向高斯圖模型結構學習問題, 考察給定潛變量時觀測變量間的條件獨立性關系.

基于懲罰似然方法, 考慮對稀疏陣S施加自適應LASSO懲罰[18], 對低秩陣L施加自適應核范數懲罰[17].自適應LASSO懲罰和自適應核范數懲罰分別為加權版的LASSO懲罰和核范數懲罰. 最初, 自適應LASSO懲罰用于線性回歸的變量選擇問題, 自適應核范數懲罰用于高維多元降秩回歸問題. 對于潛變量無向高斯圖模型問題, 自適應懲罰似然估計為

其中:S-L?0表示矩陣S-L為正定陣,L0表示矩陣L為半正定陣;λ>0和β>0均為調整參數, 參數λ控制懲罰的強度, 參數β用于權衡稀疏陣和低秩陣兩項懲罰;σ1(L)≥…≥σp(L)≥0為矩陣L的奇異值;為自適應LASSO懲罰,為自適應核范數懲罰,Sij為矩陣S中第(i,j)位置的元素.當自適應懲罰項中的權重W=(Wij)p×p及w=(w1,w2,…,wp)T滿足Wij=1,wi=1(i,j=1,2,…,p)時, 式(1)中的估計退化為非自適應懲罰似然估計.

自適應懲罰要求權重Wij(或wi)隨著|Sij|(或σi(L))的增大而減小, 以降低估計偏差.為使權重具有上下界, 本文采用類似于文獻[19]的權重形式, 即對任意的i,j=1,2,…,p,

(2)

2 優化方法

下面采用ADMM算法求解式(1)中的最小化問題, 并給出調整參數λ和β的選取方法.

2.1 ADMM算法

ADMM算法是一種迭代算法, 廣泛應用于求解線性約束下的優化問題中, 如稀疏低秩分解問題[20]. 該算法通過分解-協調過程將一個大的全局問題分解成多個小的局部子問題, 從而通過整合各子問題的解得到原問題的最優解[9]. 考慮約束條件M=S-L, 將式(1)中的目標問題轉化成ADMM形式下的優化問題:

于是, 增廣Lagrange函數為

其中Y∈p×p為Lagrange乘子,μ>0為懲罰參數, 〈·,·〉表示矩陣內積, ‖·‖F表示矩陣的Frobenius范數.給定第t次ADMM迭代的估計值(Mt,St,Lt,Yt,μt), 則第(t+1)次迭代的更新步驟為

(3)

其中ρ>1為放大因子.下面給出關于矩陣M,S,L子問題的解.

1) 更新M.

式(3)中矩陣M的優化子問題可表示為

A=Udiag(e1(A),…,ep(A))UT,

2) 更新S.

上述優化問題完全可分, 即對任意i,j=1,2,…,p, 均有

其中sign(·)為符號函數, (·)+=max{·,0}.

3) 更新L.

令B=St+1-Mt+1-Yt/μt, 則式(3)中稀疏陣L的優化子問題可寫成

對矩陣B做特征值分解

B=Qdiag(e1(B),…,ep(B))QT,

其中e1(B)≥…≥ep(B)為矩陣B的特征值,Q為正交陣,QQT=QTQ=I.從而可推得

Lt+1=Qdiag((e1(B)-λw1/μt)+,…,(ep(B)-λwp/μt)+)QT.

2.2 調整參數的選取

調整參數λ和β控制模型的復雜度, 不同的參數設定可得到不同稀疏度的S及不同秩的L.對于優化問題(1), 可采用K折交叉驗證的方式選擇最優參數組合(λ*,β*).將樣本分成K折互不相交的子集, 記為Tk(k=1,2,…,K).定義K折交叉驗證的得分函數為

3 數值模擬

在模擬實驗中, 對自適應懲罰(γ1≠0,γ2≠0)和非自適應懲罰(γ1=0,γ2=0)在含潛變量的無向高斯圖模型上的模型選擇和參數估計性能進行比較.對于自適應懲罰情況, 僅考慮γ1=γ2=1.

為評價模型選擇性能, 比較真陽率(TPR)、 陽性預測率(PPV)及馬修斯相關系數(MCC):

其中TP為真陽類個數, TN為真陰類個數, FP為假陽類個數, FN為假陰類個數. 表1列出了不同情形下TPR,PPV和MCC的均值及標準差. 由表1可見, 在所有情形下由自適應懲罰得到的PPV和MCC均顯著優于非自適應懲罰. 對于TPR, 當n=500時非自適應懲罰性能更好, 但當樣本量增大(即n=1 000)時, 兩種懲罰的TPR幾乎一樣好, 均接近于1. 因此, 基于自適應懲罰似然的潛變量圖模型結構學習性能更好.

表1 不同情形下TPR,PPV,MCC的均值及標準差Table 1 Means and standard deviations of TPR,PPV,MCC in different situations

為比較參數估計結果, 考慮矩陣M,S,L的估計誤差, 分別為

表2列出了不同情形下參數估計誤差的均值及標準差.由表2可見, 所有情形下自適應懲罰的參數估計誤差均小于非自適應懲罰.由自適應懲罰得到的參數估計更接近于真實值.

表2 不同情形下矩陣M,S,L的估計誤差均值及標準差Table 2 Means and standard deviations of estimation errors of matrices M,S,L in different situations

4 實例分析

下面采用自適應懲罰似然方法對枯草芽孢桿菌核黃素(維生素B2)生產數據集(該數據集可在文獻[21]的補充材料中下載)進行圖模型結構學習. 該數據集的樣本量為71, 變量個數為4 089, 其中有4 088個變量表示不同基因表達水平的對數, 1個變量表示核黃素生產率的對數. 本文僅對文獻[21]中簡化后的數據集riboflavinv100.csv做圖模型推斷, 其中包含經驗方差最大的100個基因表達水平變量和1個測量核黃素生產率的變量.

為比較含潛變量和不含潛變量(即固定低秩陣L=0)時觀測變量間的圖結構, 將樣本分為訓練集和測試集, 其中訓練集包含57個樣本, 測試集包含14個樣本.首先, 基于訓練集樣本得到稀疏陣和低秩陣的估計; 然后, 基于測試集樣本計算負的對數似然.圖1為由自適應懲罰(γ1=γ2=1)得到的圖模型選擇結果, 其中(A)為給定潛變量時觀測變量間的條件獨立性關系, (B)為不考慮潛變量時觀測變量間的條件獨立性關系. 圖1中橫軸從左至右(縱軸從上至下)依次對應第1～101個觀測變量, 黑色表示對應觀測變量間條件相關, 白色表示對應變量間條件獨立. 當考慮潛變量時, 估計出的圖模型共包含531條邊(占觀測變量對總數的10.5%)以及13個潛變量(對應低秩陣的秩), 基于測試集的負對數似然為49.54. 當不考慮潛變量時, 估計出的圖模型共包含1104條邊(占觀測變量對總數的21.9%), 基于測試集的負對數似然為50.48. 顯然, 考慮潛變量時能得到更稀疏的圖模型, 且負對數似然更小, 表明考慮潛變量對觀測變量的影響能更好地擬合枯草芽孢桿菌核黃素生產數據集.

綜上所述, 本文提出了一種基于自適應懲罰似然的潛變量高斯圖模型結構學習方法. 首先, 在觀測似然后面加上兩項自適應懲罰項, 分別為稀疏陣部分的自適應LASSO懲罰以及低秩陣部分的自適應核范數懲罰; 然后, 通過ADMM算法最小化懲罰似然以求解稀疏陣和低秩陣的參數估計, 從而得到給定潛變量時觀測變量間的條件獨立關系. 在模擬實驗中, 通過與非自適應懲罰的比較, 驗證了自適應懲罰似然方法在模型選擇和參數估計方面性能均更好.