鄭倩貞, 徐平峰

(1. 長(zhǎng)春工業(yè)大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 長(zhǎng)春 130012; 2. 東北師范大學(xué) 前沿交叉研究院, 長(zhǎng)春 130024)

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

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

研究表明, 非自適應(yīng)懲罰得到的估計(jì), 包括LASSO懲罰得到的稀疏協(xié)方差逆陣估計(jì)[16]以及核范數(shù)懲罰在降秩回歸問題中得到的低秩系數(shù)陣估計(jì)[17], 均存在偏差較大的問題, 而自適應(yīng)懲罰通?？山档凸烙?jì)偏差. 因此, 本文采用自適應(yīng)懲罰似然方法處理含潛變量的高斯圖模型結(jié)構(gòu)學(xué)習(xí)問題, 對(duì)稀疏陣和低秩陣進(jìn)行估計(jì), 以得到給定潛變量時(shí)觀測(cè)變量間的條件獨(dú)立關(guān)系. 本文對(duì)稀疏陣部分施加自適應(yīng)LASSO懲罰[18], 對(duì)低秩陣部分施加自適應(yīng)核范數(shù)懲罰[17]. 與非自適應(yīng)懲罰的優(yōu)化問題類似, 采用ADMM算法優(yōu)化求解自適應(yīng)懲罰似然的最小化問題, 并且在求解過程中仍具有顯式表達(dá)式, 以確保算法的計(jì)算效率. 本文模擬比較了自適應(yīng)懲罰與非自適應(yīng)懲罰在潛變量高斯圖模型結(jié)構(gòu)學(xué)習(xí)和參數(shù)估計(jì)上的性能. 結(jié)果表明, 自適應(yīng)懲罰顯著優(yōu)于非自適應(yīng)懲罰, 有效降低了稀疏陣和低秩陣的估計(jì)偏差, 能更準(zhǔn)確地學(xué)習(xí)觀測(cè)變量間的條件獨(dú)立關(guān)系.

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

無向高斯圖模型是指與無向圖G=(V,E)相關(guān)的多元正態(tài)分布模型, 其中頂點(diǎn)集V中的每個(gè)頂點(diǎn)v∈V表示一個(gè)高斯隨機(jī)變量Xv, 邊集E中的每條邊(u,v)∈E蘊(yùn)含了給定其余隨機(jī)變量Xv′(v′∈V{u,v})時(shí)隨機(jī)變量Xu與Xv之間的條件相關(guān)性.當(dāng)采用無向高斯圖模型對(duì)實(shí)際問題進(jìn)行分析時(shí), 可能會(huì)出現(xiàn)隨機(jī)變量為潛變量的情況, 即變量不存在觀測(cè)值.因此, 本文考慮含潛變量的無向高斯圖模型結(jié)構(gòu)學(xué)習(xí)問題, 考察給定潛變量時(shí)觀測(cè)變量間的條件獨(dú)立性關(guān)系.

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

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

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

(2)

2 優(yōu)化方法

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

2.1 ADMM算法

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

于是, 增廣Lagrange函數(shù)為

其中Y∈p×p為L(zhǎng)agrange乘子,μ>0為懲罰參數(shù), 〈·,·〉表示矩陣內(nèi)積, ‖·‖F(xiàn)表示矩陣的Frobenius范數(shù).給定第t次ADMM迭代的估計(jì)值(Mt,St,Lt,Yt,μt), 則第(t+1)次迭代的更新步驟為

(3)

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

1) 更新M.

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

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

2) 更新S.

上述優(yōu)化問題完全可分, 即對(duì)任意i,j=1,2,…,p, 均有

其中sign(·)為符號(hào)函數(shù), (·)+=max{·,0}.

3) 更新L.

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

對(duì)矩陣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 調(diào)整參數(shù)的選取

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

3 數(shù)值模擬

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

為評(píng)價(jià)模型選擇性能, 比較真陽率(TPR)、 陽性預(yù)測(cè)率(PPV)及馬修斯相關(guān)系數(shù)(MCC):

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

表1 不同情形下TPR,PPV,MCC的均值及標(biāo)準(zhǔn)差Table 1 Means and standard deviations of TPR,PPV,MCC in different situations

為比較參數(shù)估計(jì)結(jié)果, 考慮矩陣M,S,L的估計(jì)誤差, 分別為

表2列出了不同情形下參數(shù)估計(jì)誤差的均值及標(biāo)準(zhǔn)差.由表2可見, 所有情形下自適應(yīng)懲罰的參數(shù)估計(jì)誤差均小于非自適應(yīng)懲罰.由自適應(yīng)懲罰得到的參數(shù)估計(jì)更接近于真實(shí)值.

表2 不同情形下矩陣M,S,L的估計(jì)誤差均值及標(biāo)準(zhǔn)差Table 2 Means and standard deviations of estimation errors of matrices M,S,L in different situations

4 實(shí)例分析

下面采用自適應(yīng)懲罰似然方法對(duì)枯草芽孢桿菌核黃素(維生素B2)生產(chǎn)數(shù)據(jù)集(該數(shù)據(jù)集可在文獻(xiàn)[21]的補(bǔ)充材料中下載)進(jìn)行圖模型結(jié)構(gòu)學(xué)習(xí). 該數(shù)據(jù)集的樣本量為71, 變量個(gè)數(shù)為4 089, 其中有4 088個(gè)變量表示不同基因表達(dá)水平的對(duì)數(shù), 1個(gè)變量表示核黃素生產(chǎn)率的對(duì)數(shù). 本文僅對(duì)文獻(xiàn)[21]中簡(jiǎn)化后的數(shù)據(jù)集riboflavinv100.csv做圖模型推斷, 其中包含經(jīng)驗(yàn)方差最大的100個(gè)基因表達(dá)水平變量和1個(gè)測(cè)量核黃素生產(chǎn)率的變量.

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

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