基于GEMS算法的潛變量高斯圖模型結構學習

鄭倩貞，徐平峰，曹 蕾

(長春工業大學數學與統計學院，吉林 長春 130012)

高斯圖模型能夠清晰直觀地反應變量間的相互關系，被廣泛應用于高維情形.在對實際問題進行圖模型結構學習時，僅考慮觀測變量有時并不能正確反應變量間的相互關系，因此需考慮潛變量對可觀測變量的影響，在給定潛變量時探討可觀測變量間的條件關系.Chandrasekaran等[1]將可觀測變量的邊緣協方差陣的逆陣分解為一個稀疏陣和一個低秩陣，提出了懲罰似然的方法，對稀疏陣和低秩陣分別施加1范數及核范數懲罰，并結合凸優化和代數幾何對潛變量圖模型選擇問題進行了研究.Yuan[2]基于Chandrasekaran等人的研究，將懲罰似然的核范數懲罰項替換為對低秩陣的秩的約束條件，提出了潛變量GLasso(LVglasso)方法，并結合EM算法對高維情形下的潛變量圖模型選擇問題進行了模擬研究.Lauritzen等[3]對懲罰似然做了與Yuan相似的處理，采用插補的方法，結合EM算法和GLasso算法進行模型選擇.

但上述與EM相結合的方法需要先給定1組正則化參數，然后對每個正則化參數利用EM算法求懲罰似然的最小值點.如果正則化參數選取不當，將會導致每次迭代的模型離真模型越來越遠，而且增加計算時間.本文基于期望模型選擇(EMS)算法[4]的思想，在每次迭代時從候選模型中選取期望信息準則最小的模型作為下一步的當前模型，下一次迭代時在當前模型下求候選模型的期望信息準則的值.但由于可能的模型太多，在模型選擇時遍歷全部模型不可行，因此只選出部分模型作為候選模型.這里的候選模型也可以通過1組正則化參數來確定，但每次的正則化參數不一定相同.稱這種方法為廣義期望模型選擇(GEMS)算法.模擬實驗顯示，基于GEMS的LVglasso方法收斂速度快，計算時間短.

1 高斯圖模型

Yuan等[5]提出通過最小化負1懲罰對數似然的方法去估計高斯圖模型的協方差逆陣Ω，懲罰似然為

2 LVglasso方法

其中：S-L?0表示S-L為正定矩陣；L0表示L為非負定矩陣；為ΩO的估計，為的估計；為可觀測樣本的對數似然函數，即

基于Chandrasekaran等[1]提出的懲罰似然，Yuan[2]提出了計算更加方便的LVglasso方法：

其中：0≤r≤p，S?=S-diag(S).限制條件rank(L)≤r相當于假設存在r個潛變量.

考慮完全數據x=(x1，…，xn)T=(xO，xH)，xi=(xO，i，xH，i)T，其中xO，i為第i個樣本的可觀測數據，xH，i為第i個樣本的不可觀測數據，i=1，…，n.Ω的LVglasso估計為

3 GEMS算法

EMS算法[4]是一種迭代算法，用于處理缺失數據情形下的模型選擇問題.該算法的每次迭代都需給定當前模型Mc和當前模型下的參數θc∈ΘMc，并依次進行期望步(E步)和模型選擇步(MS步)，直至滿足停止準則得到最優的模型估計M*和參數估計θ*∈ΘM*.本文的GEMS算法與EMS算法類似，不同之處在于：GEMS算法的MS步不遍歷全部模型，而是通過GLasso算法找出候選模型，候選模型可由1組正則化參數來確定，且每次迭代的正則化參數可能不同.在這些候選模型中選擇期望BIC最小的模型，然后將該模型及其對應的參數作為下一次迭代的當前模型和當前參數.從部分而非全部模型中選擇最優模型可大大減小計算成本，提高計算效率，尤其是在高維問題中.考慮潛變量高斯圖模型的結構學習問題，現有當前模型G(t)及當前參數Ω(t)，則第(t+1)次迭代如下：

(1) E步

進而，有Q函數

Q(G，Ω|G(t)，Ω(t))=Et(-2(G，Ω))+log(n)dfG= -nlogdetΩ+ntr[ΩEt(Σn)]+log(n)dfG.

(2) MS步

因此，可得到1組候選模型G={G(t)，Gλm，m=1，…，k}.針對每個候選模型G∈G，計算Q函數，從而得到Ω的估計Ω(t+1)=argminΩQ(G，Ω|G(t)，Ω(t))及其對應的圖模型G(t+1)=argminG∈GQ(G，Ω(t+1)|G(t)，Ω(t)).

將G(t+1)，Ω(t+1)作為下一次迭代的當前模型和當前參數，重復以上步驟直至滿足停止準則.

4 模擬實驗及實例分析

在模擬實驗中，考慮不同情形下潛變量圖模型的結構學習問題，對GEMS算法和EM算法在LVglasso估計求解問題上的模擬結果進行了比較.模擬實驗覆蓋p=48，98，148，198，h=2，r=2，5，n=500，1 000共16種情形，每種情形各模擬50次.真模型產生機制與Yuan[2]類似，兩者的不同之處體現在對潛變量的設定上.在本文真模型中，每個潛變量至少和2個可觀測變量、至多和(p-1)個可觀測變量有關.值得注意的是，當在進行n=1 000，p=198，r=2情形設定下的第38次EM算法模擬時，由R中lvglasso函數產生的Ω迭代初值為非對稱陣，所以該種情形只模擬了37次.16種情形的CPU平均運行時間如表1所示.

表1 不同情形下CPU平均運行時間 s

從表1中可看出，對于任意一種情況，EM算法的運行時長都要遠大于GEMS算法，達到5倍、10倍，甚至是15倍的差距.GEMS算法大大提升了潛變量圖模型選擇的速度.本文用于評價算法性能的指標為：

其中tp，tn，fp，fn分別為真陽類、真陰類、假陽類、假陰類的個數.圖1給出了所有情形下tpr，ppv和mcc的箱線圖.總體上看，GEMS較EM有更優的表現，但在極個別情況如n=500，p=198，r=2或5時EM的tpr值較大.同時可看出，樣本量越大，潛變量個數的假設越接近真實模型，模型推斷就越準確.

白色箱子代表EM算法，灰色箱子代表GEMS算法，橫軸為真模型可觀測變量的個數.

基于GEMS算法，對Wille等[7]論文中擬南芥植物類異戊二烯生物合成相關基因的數據進行了潛變量高斯圖模型結構學習，估計了各基因間的條件相關性.該數據的數據來源為https：∥static-content.springer.com/esm/art%3A10.1186%2Fgb-2004-5-11-r92/MediaObjects/13059_2004_896_MOESM1_ESM.txt，數據中共有118個樣本，每個樣本包含39個基因表達.若假設的潛變量個數不同，則推斷出的各基因間的條件相關性也不同.假設潛變量個數為r=1或r=3時的估計結果如圖2所示.當r=1時，共估計出174條邊，算法運行時間約為43 s；當r=3時，共估計出38條邊，算法運行時間約為10 s.

圖2 r=1和r=3時的基因圖模型

5 結語

本文簡要介紹了高斯圖模型及潛變量高斯圖模型下的LVglasso方法，給出了GEMS算法結合LVglasso下潛變量圖模型結構學習的算法步驟，并從模擬實驗的角度比較了GEMS算法和EM算法在潛變量圖模型選擇問題上的優劣.通過多種不同情形下的模擬實驗，可以發現，樣本量越大，潛變量個數的假設越接近真實模型，模型推斷就越準確.結合tpr，ppv，mcc以及CPU平均運行時間，無論在何種模擬情形，GEMS算法在模型選擇上的表現較EM算法優越.