黃飛虎，陳松燦

南京航空航天大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，南京 211106

1 引言

無向概率圖模型是一類用于刻畫一組隨機(jī)變量之間條件相關(guān)性的強(qiáng)大統(tǒng)計(jì)工具，目前已被廣泛應(yīng)用于機(jī)器學(xué)習(xí)、計(jì)算機(jī)視覺、生物信息學(xué)與社會學(xué)等領(lǐng)域[1-4]。高斯圖模型（Guassian graphical model，GGM）為一類流行的無向概率圖模型，能很好地刻畫一組正態(tài)分布隨機(jī)變量的條件相關(guān)性。具體地，假設(shè)隨機(jī)向量x=(x1,x2,…,xp)T∈Rp服從多元正態(tài)分布N(μ,Σ)，與之對應(yīng)的無向圖為G(V,E)，其中V={x1,x2,…,xp}為頂點(diǎn)集，E=V×V代表邊集，那么對于任意(i,j)?E，xi⊥xj|x(i,j)表示隨機(jī)變量xi與xj條件獨(dú)立。對于(i,j)?E當(dāng)且僅當(dāng)(Σ-1)ij=0，即協(xié)方差矩陣逆（也稱為精度矩陣）的(i,j)元素為0。因此，精度矩陣的稀疏模式能刻畫圖模型的結(jié)構(gòu)。由此可知，概率圖模型結(jié)構(gòu)的估計(jì)可等價(jià)于稀疏精度矩陣的估計(jì)。目前已存在大量對于圖模型與精度矩陣估計(jì)的工作[5-9]，它們大致可分為三類：第一類通過利用其他變量來稀疏擬合每個(gè)變量而得到每個(gè)點(diǎn)的近鄰。例如，文獻(xiàn)[5]通過利用套索模型（Lasso[10]）擬合每個(gè)變量而提出了近鄰選擇估計(jì)器，該方法可視為一種偽似然估計(jì)方法。第二類通過直接最小化?1范數(shù)懲罰的負(fù)對數(shù)似然。例如，文獻(xiàn)[3,6]通過直接求解?1范數(shù)懲罰的負(fù)對數(shù)似然估計(jì)高斯圖模型。文獻(xiàn)[7]利用有效的塊坐標(biāo)下降方法求解該?1范數(shù)懲罰的對數(shù)似然問題，提出了著名的圖套索（graphical Lasso）。第三類通過利用樣本協(xié)方差直接估計(jì)稀疏精度矩陣。例如，文獻(xiàn)[8]通過求解一系列稀疏線性規(guī)劃問題估計(jì)稀疏精度矩陣。文獻(xiàn)[9]提出了一個(gè)帶約束的?1范數(shù)最小估計(jì)器（constrained?1-minimization for inverse matrix estimation，CLIME）估計(jì)稀疏精度矩陣。

盡管GGM能很好地刻畫正態(tài)分布的數(shù)據(jù)，但其要求正態(tài)分布假設(shè)過于苛刻。事實(shí)上，人們所采集到的數(shù)據(jù)往往面臨兩類問題：（1）數(shù)據(jù)很少嚴(yán)格服從正態(tài)分布；（2）數(shù)據(jù)通常含有少量噪聲。為了處理問題（1），文獻(xiàn)[11]將正態(tài)分布推廣到非參數(shù)正態(tài)分布（non-parameter normal distribution，nonparanormal），進(jìn)而提出了一類半?yún)?shù)概率圖模型。具體地，如果存在一些單變量的單調(diào)可微函數(shù){fi}p i=1，有f(x)=(f1(x1),f2(x2),…,fp(xp))T服從多元正態(tài)分布N(0,Σ)，那么x=(x1,x2,…,xp)T服從非參數(shù)正態(tài)分布NPN(0,Σ,f)。同時(shí)，由單變量函數(shù){fi}p i=1的單調(diào)可微性，稀疏精度矩陣Σ-1同樣刻畫了隨機(jī)變量(x1,x2,…,xp)的相關(guān)性，即給定其他變量xi與xj條件獨(dú)立當(dāng)且僅當(dāng)(Σ-1)ij=0。為了同時(shí)解決問題（1）與（2），文獻(xiàn)[12-13]采用基于非參排序的統(tǒng)計(jì)量（Spearman’s rho 或Kendall’s tau）估計(jì)相關(guān)矩陣，提出了魯棒的估計(jì)方法用于學(xué)習(xí)半?yún)?shù)概率圖模型。總之，這些半?yún)D模型的估計(jì)方法的基本流程為：首先利用基于截?cái)嗟恼龖B(tài)計(jì)分（normal scoring[11]）或基于非參排序的統(tǒng)計(jì)量[12-13]估計(jì)出相關(guān)矩陣，然后把它代入現(xiàn)有圖模型估計(jì)器，學(xué)習(xí)出稀疏的精度矩陣，即得到相應(yīng)的圖結(jié)構(gòu)。

到目前為止，上述圖模型的建立均基于同一分布數(shù)據(jù)，因此不適合刻畫異構(gòu)性或動態(tài)性的數(shù)據(jù)。例如，采集了包含正常與病狀的腦影像數(shù)據(jù)[14]，如果利用上述圖模型分別構(gòu)建正常大腦與病狀大腦的各自腦網(wǎng)絡(luò)，則會忽略它們的共性結(jié)構(gòu)；如果利用上述圖模型總體估計(jì)單個(gè)腦網(wǎng)絡(luò)，則會忽略它們之間的差異結(jié)構(gòu)。因此，為了能更好地挖掘這些異構(gòu)數(shù)據(jù)的結(jié)構(gòu)信息，聯(lián)合學(xué)習(xí)多個(gè)圖模型已成為一個(gè)研究主題，典型的工作有文獻(xiàn)[14-20]。例如，文獻(xiàn)[15]利用層次稀疏結(jié)構(gòu)懲罰能很好學(xué)習(xí)出多個(gè)圖模型的共性結(jié)構(gòu)。文獻(xiàn)[16]通過利用組套索（group Lasso）[21]與兩兩融合套索（fused Lasso）[22]的結(jié)構(gòu)懲罰學(xué)習(xí)多個(gè)圖模型的共性結(jié)構(gòu)，提出了聯(lián)合圖套索（joint graphical Lasso）。同時(shí)，文獻(xiàn)[14]利用有序融合套索聯(lián)合學(xué)習(xí)多個(gè)有序的概率圖模型。為了使這些聯(lián)合圖模型能更好地勝任矩陣變量的數(shù)據(jù)，如腦功能性磁共振成像（functional magnetic resonance imaging，fMRI）數(shù)據(jù)及股票交易數(shù)據(jù)等，文獻(xiàn)[20]提出了聯(lián)合矩陣變量的高斯圖模型。另外，針對隨著時(shí)間光滑變化的異構(gòu)數(shù)據(jù)，目前工作[23-25]提出了相應(yīng)的動態(tài)GGM學(xué)習(xí)動態(tài)的條件相關(guān)性?？傮w上，這些工作均利用核光滑方法估計(jì)出相應(yīng)的協(xié)方差矩陣，再把已估計(jì)的協(xié)方差矩陣代入已有的圖模型估計(jì)器得到相應(yīng)的動態(tài)圖結(jié)構(gòu)。

同樣地，盡管上述聯(lián)合或動態(tài)的圖模型能較好地分析異構(gòu)數(shù)據(jù)的條件相關(guān)性，但是它們均建立在嚴(yán)格的正態(tài)分布假設(shè)下。由于當(dāng)前高維的異構(gòu)數(shù)據(jù)通常很難嚴(yán)格服從正態(tài)分布且常含噪聲，上述聯(lián)合的與動態(tài)的圖模型仍然很難勝任這些異構(gòu)數(shù)據(jù)。例如，對于采集不同病狀下的腦影像數(shù)據(jù)，由于疾病的易變性通常使其服從一些尾部較重的分布。另外，在采集數(shù)據(jù)過程中由于儀器不穩(wěn)定，再加上志愿者頭部的運(yùn)動及呼吸心跳的影響，往往采集到的數(shù)據(jù)都帶有一定噪聲。為了處理上述問題，本文提出聯(lián)合半?yún)?shù)圖模型學(xué)習(xí)這些異構(gòu)數(shù)據(jù)的條件相關(guān)性。同時(shí)，針對光滑變化的異構(gòu)數(shù)據(jù)（如時(shí)序的fMRI數(shù)據(jù)），提出聯(lián)合的動態(tài)半?yún)?shù)圖模型。在建模上，將基于非參排序的相關(guān)矩陣估計(jì)方法與結(jié)構(gòu)融合圖套索方法相結(jié)合，提出了半?yún)?shù)融合圖套索估計(jì)器。特別針對動態(tài)圖模型，提出了一種新的核光滑Kendall’s tau相關(guān)矩陣?？傊?，本文主要貢獻(xiàn)如下：

（1）提出了聯(lián)合的半?yún)?shù)圖模型用于學(xué)習(xí)非正態(tài)分布異構(gòu)數(shù)據(jù)的條件相關(guān)性，且其較目前已有的聯(lián)合圖模型更靈活、魯棒。

（2）進(jìn)一步針對光滑變化的動態(tài)異構(gòu)數(shù)據(jù)，提出了聯(lián)合動態(tài)半?yún)?shù)圖模型。

（3）采用了有效的ADMM（alternating direction method of multipliers）方法對提出的模型進(jìn)行求解。

（4）利用一些人工數(shù)據(jù)與真實(shí)數(shù)據(jù)（如腦影像、股票交易數(shù)據(jù)）同時(shí)驗(yàn)證了模型的有效性。

2 相關(guān)工作

本文首先介紹非參數(shù)正態(tài)分布與半?yún)?shù)概率圖模型。存在一系列單值單調(diào)且可微函數(shù){fi}p i=1與對稱正定矩陣Σ且diag(Σ)=I，那么稱隨機(jī)向量x=(x1,x2,…,xp)T服從非參數(shù)正態(tài)分布NPN(0,Σ,f)，當(dāng)且僅當(dāng)f(x)=(f1(x1),f2(x2),…,fp(xp))服從多元正態(tài)分布N(0,Σ)。文獻(xiàn)[12-14]證明矩陣Ω=Σ-1的稀疏模式能刻畫x=(x1,x2,…,xp)T的條件相關(guān)性（即Ωij=0?xi⊥xj|x{}i,j），且基于該非參數(shù)正態(tài)分布提出了半?yún)?shù)圖模型。

下面介紹半?yún)?shù)圖模型的估計(jì)方法。文獻(xiàn)[11]提出了基于正態(tài)計(jì)分的半?yún)?shù)圖模型估計(jì)方法，而文獻(xiàn)[12-13]提出了一類基于非參排序方法估計(jì)該半?yún)?shù)模型，其不僅比基于正態(tài)計(jì)分的方法具有更優(yōu)的收斂率，且更加魯棒。具體地，首先利用基于非參排序的統(tǒng)計(jì)量（Spearman’s rho 或Kendall’s tau）估計(jì)相關(guān)矩陣Σ，然后將它代入已有圖模型估計(jì)精度矩陣Ω=Σ-1，即半?yún)?shù)圖稀疏結(jié)構(gòu)。例如，基于非參排序的Kendall’s tau相關(guān)系數(shù)τkl估計(jì)如下：

然后相關(guān)矩陣Σ=(Σkl)通過Kendall’s tau相關(guān)系數(shù)矩陣Γ?=(τ?kl)估計(jì)得到[26-27]，其中：

3 聯(lián)合半?yún)?shù)圖模型

下面提出聯(lián)合半?yún)?shù)圖模型用于學(xué)習(xí)非正態(tài)分布異構(gòu)數(shù)據(jù)的條件相關(guān)性。該問題等價(jià)于學(xué)習(xí)多個(gè)具有一些共性結(jié)構(gòu)的半?yún)?shù)圖模型。事實(shí)上，本文研究聯(lián)合半?yún)?shù)圖模型的動機(jī)源于一些重要的應(yīng)用。例如，利用一些來自同一種病多種亞型的腦影像數(shù)據(jù)[14]，通過聯(lián)合學(xué)習(xí)不同病狀的腦網(wǎng)絡(luò)可挖掘出疾病的發(fā)展情況。

具體地，給定K類獨(dú)立同分布樣本服從非參數(shù)正態(tài)分布NPN(0,Σk,fk)，[K]={1,2,…,K}。通常為了獲得稀疏的圖結(jié)構(gòu)，求解下面的?1范數(shù)懲罰的負(fù)對數(shù)似然問題：

其中為的樣本協(xié)方差矩陣。函數(shù)為存在且未知的隱函數(shù)，因此不能直接求得相關(guān)矩陣那么類似于文獻(xiàn)[12-13]采用基于非參排序方法直接估計(jì)它。具體地，可通過上述Kendall’s tau統(tǒng)計(jì)量估計(jì)每類的相關(guān)矩陣

考慮到多個(gè)半?yún)?shù)圖模型共享一些結(jié)構(gòu)，即精度矩陣共享一些稀疏結(jié)構(gòu)，因此提出了半?yún)?shù)融合圖套索方法聯(lián)合估計(jì)這些圖模型。具體地，求解如下的結(jié)構(gòu)正則化的負(fù)對數(shù)似然問題：

其中，為負(fù)對數(shù)似然項(xiàng)；為稀疏懲罰項(xiàng)，使得每個(gè)圖模型稀疏；P(Ω)=為有序融合套索懲罰項(xiàng)，使得相鄰的圖模型更相似。這里λ1與λ2為非負(fù)的正則化參數(shù)，其中λ1控制每個(gè)圖的稀疏率，而λ2控制相鄰的圖相似程度。當(dāng)λ2=0時(shí)，問題（1）可解耦為K個(gè)稀疏正則化的負(fù)對數(shù)似然問題，那么該聯(lián)合模型退化為半?yún)?shù)圖模型[12-13]。

4 聯(lián)合動態(tài)半?yún)?shù)圖模型

下面提出聯(lián)合動態(tài)半?yún)?shù)圖模型用于學(xué)習(xí)光滑變化的非正態(tài)分布異構(gòu)數(shù)據(jù)的條件相關(guān)性，其動機(jī)源于一些有意義的應(yīng)用。例如，利用時(shí)序的fMRI數(shù)據(jù)學(xué)習(xí)人類不同年齡段的腦網(wǎng)絡(luò)[25]，以了解大腦發(fā)育情況。首先，定義一類新的動態(tài)半?yún)?shù)概率圖模型。

定義1（動態(tài)半?yún)?shù)圖模型）如果隨機(jī)變量對(X,T)服從動態(tài)半?yún)?shù)概率圖模型，其相應(yīng)的動態(tài)圖為G(t)=(V,E(t))，那么其滿足如下條件：

（1）X|T=t～NPN(0,Σ(t),f)，其中T～g(t)為定義在[0,1]上的連續(xù)函數(shù)；

（2）動態(tài)圖G(t)=(V,E(t))包括固定點(diǎn)集合V，動態(tài)邊集合E(t)，其中邊的權(quán)重隨著時(shí)間變量t∈[0,1]變化，且其圖結(jié)構(gòu)也可以隨之改變，即精度矩陣Ω(t)隨著時(shí)間t變化；

（3）xi⊥xj|{x{i,j},T=t}當(dāng)且僅當(dāng) (i,j)?E(t)。

人們感興趣的時(shí)間變量T屬于有界區(qū)間，因此其可以轉(zhuǎn)化到區(qū)間[0,1]。不失一般性，本文均假設(shè)t∈[0,1]。接下來，為了估計(jì)該動態(tài)半?yún)?shù)圖模型，利用一種新的核光滑Kendall’s tau相關(guān)系數(shù)矩陣Γ(t)=(τkl(t))。具體地，當(dāng)每個(gè)時(shí)間點(diǎn)t∈[0,1]（即每個(gè)分布）只采一個(gè)樣本時(shí)，核光滑Kendall’s tau相關(guān)系數(shù)τkl(t)估計(jì)如下：

其中，ω(t,ti,tj)=Kh(t-ti)Kh(t-tj)。

當(dāng)每個(gè)時(shí)間點(diǎn)t∈[0,1]采m≥2個(gè)獨(dú)立同分布樣本時(shí)，核光滑Kendall’s tau相關(guān)系數(shù)τkl(t)估計(jì)如下：

這里，Kh(·)=K(·/h)為對稱核函數(shù)，其中h>0 為帶寬參數(shù)。例如，高斯核，其中帶寬參數(shù)h控制圍繞時(shí)間點(diǎn)ti的窗口。具體地，較小的h表明估計(jì)的圖模型隨時(shí)間變化的頻率較高，而較大的h表明估計(jì)的圖模型隨時(shí)間變化的頻率較低。然后，相關(guān)矩陣Σ(t)由核光滑Kendall’s tau相關(guān)系數(shù)矩陣Γ?(t)=(τ?kl(t))估計(jì)可得，具體為：

最后，把已估計(jì)出的相關(guān)矩陣代入已有的圖模型估計(jì)器（如graphical lasso[7]或CLIME[9]）可以得到稀疏精度矩陣Ω(t)，即動態(tài)圖結(jié)構(gòu)。

考慮到動態(tài)圖模型隨著時(shí)間變化依然保持一定的共性結(jié)構(gòu)，本文采用上文的聯(lián)合學(xué)習(xí)思想，提出聯(lián)合的動態(tài)半?yún)?shù)圖模型。具體為，把已估計(jì)出的相關(guān)矩陣 {Σ?(tk)}K k=1代入上文提出的半?yún)?shù)融合圖套索估計(jì)器（1），可以聯(lián)合估計(jì)多個(gè)時(shí)間點(diǎn)的精度矩陣即稀疏圖結(jié)構(gòu)。

5 模型優(yōu)化

本文利用交替方向乘子方法（ADMM[28]）求解問題（1）。ADMM是一類非常適用于求解帶等式約束問題的優(yōu)化方法，可表示如下：

其中，λ為拉格朗日乘子；ρ為懲罰參數(shù)。首先給出上述問題（2）的增廣拉格朗日函數(shù)：

那么ADMM采用Gauss-Seidel迭代求解問題（2），在第t+1步迭代表示如下：

下面應(yīng)用ADMM具體求解問題（1）。首先把問題（1）改寫為如下等式約束問題：

問題（3）的增廣拉格朗日函數(shù)可表示如下：

然后利用ADMM求解問題（3），在第t+1步迭代表示如下：

接下來，將分別介紹問題（4a）與（4b）的具體求解。首先對于問題（4a），其可以分解為K個(gè)獨(dú)立問題。對于k∈[K]：

然后對其目標(biāo)函數(shù)微分得到：

易知Ωk與矩陣Ak=Σ?k-Λk-ρZk共享特征向量，且其特征值滿足如下關(guān)系：

其中，{αi}ip=1為矩陣Ωk的特征值；{βi}ip=1為矩陣Ak的特征值。因此，對矩陣Ak進(jìn)行特征值分解為Ak=UkBkUk，那么可得Ωk=UkDkUk，其中Dk為特征值{αi}ip=1組成的對角矩陣。

同樣，問題（4b）可以分解p2個(gè)獨(dú)立的融合套索問題。

對于1≤i,j≤p：

且子問題（5）可用標(biāo)準(zhǔn)融合套索的近似算子求解[29]。由于{Zk}為對稱矩陣，只要求解個(gè)子問題（5）。由于問題（4a）與（4b）均可分解為獨(dú)立的子問題，可以考慮利用并行框架來加速本文算法。

6 人工數(shù)據(jù)實(shí)驗(yàn)

下面利用一些人工數(shù)據(jù)驗(yàn)證本文模型的有效性。具體地，對于學(xué)習(xí)異構(gòu)數(shù)據(jù)的條件相關(guān)性，即學(xué)習(xí)多個(gè)半?yún)?shù)概率圖模型（semi-parameter probability graphical model，SPGM），本文的聯(lián)合半?yún)?shù)圖模型（joint semi-parameter graphical Lasso，JSPGL）將與標(biāo)準(zhǔn)的半?yún)?shù)圖模型[12-13]（semi-parameter graphical Lasso，SPGL）及聯(lián)合的GGM[14,16]（joint graphical Lasso，JGL）比較。對于動態(tài)的異構(gòu)數(shù)據(jù)條件相關(guān)性，即學(xué)習(xí)多個(gè)動態(tài)半?yún)?shù)概率圖模型（dynamic semi-parameter probability graphical model，DSPGM），本文的聯(lián)合動態(tài)半?yún)?shù)圖模型（joint dynamic semi-parameter graphical Lasso，JDSPGL）將與動態(tài)半?yún)?shù)圖模型（dynamic semiparameter graphical Lasso，DSPGL）及動態(tài)的GGM[21-23]（dynamical graphical Lasso，DGL）比較。

在實(shí)驗(yàn)中，為了突出本文模型的有效性，讓DGL與融合圖套索框架結(jié)合來參與比較。文中所有模型參數(shù)通過十重交叉驗(yàn)證得到。同時(shí)，所有實(shí)驗(yàn)均重復(fù)50次，下面報(bào)告的實(shí)驗(yàn)結(jié)果為其平均值。另外，上述所有動態(tài)圖模型，均選擇帶寬參數(shù)h=1。最后，所有算法均在Matlab軟件平臺上運(yùn)行，且在英特爾i5-3470處理器、16 GB內(nèi)存的計(jì)算機(jī)上執(zhí)行。

6.1 人工數(shù)據(jù)的生成

本節(jié)介紹一些人工數(shù)據(jù)的生成。不失一般性，本文只關(guān)注學(xué)習(xí)Erd?s-Rényi（ER）網(wǎng)絡(luò)。具體地，首先生成一個(gè)稀疏率92%的ER網(wǎng)絡(luò)，然后由生成的ER網(wǎng)絡(luò)復(fù)制K份，再對每個(gè)網(wǎng)絡(luò)隨機(jī)減少p/4個(gè)邊，最后得到多個(gè)具有一定相似結(jié)構(gòu)的網(wǎng)絡(luò)，其相應(yīng)的鏈接矩陣為為了使得這些鏈接矩陣符合精度矩陣，進(jìn)行如下賦值：

其中，Ek,k∈[K]表示圖邊集合。最后，在矩陣的對角元素加上相應(yīng)的正數(shù)以保證它們對稱正定。為了方便，令n=n1=n2=…=nK。接下來，讓每個(gè)正態(tài)分布生成n個(gè)數(shù)據(jù)點(diǎn)為了驗(yàn)證本文模型對正態(tài)分布假設(shè)的放松，與文獻(xiàn)[12]類似，再對數(shù)據(jù)進(jìn)行高斯累積分布函數(shù)轉(zhuǎn)化，如下：

其中為標(biāo)準(zhǔn)的高斯累積分布函數(shù)。因此，得到轉(zhuǎn)化數(shù)據(jù)服從非參數(shù)正態(tài)分布f)，k∈[K]。

下面介紹產(chǎn)生動態(tài)半?yún)?shù)圖模型的過程。同樣地，首先生成一個(gè)稀疏率92%的ER網(wǎng)絡(luò)，然后由生成的ER網(wǎng)絡(luò)復(fù)制n個(gè)，再對每個(gè)網(wǎng)絡(luò)隨機(jī)減少p/4個(gè)邊，最后得到多個(gè)具有一定相似結(jié)構(gòu)的網(wǎng)絡(luò)，其相應(yīng)的鏈接矩陣為為了使得這些鏈接矩陣符合動態(tài)結(jié)構(gòu)，對其進(jìn)行如下賦值：

其中，t∈[0,1]。同時(shí)在矩陣{Ω(tk)}n k=1的對角元素加上相應(yīng)的正數(shù)以保證它們對稱正定。接下來，讓每個(gè)正態(tài)分布N(0,Ω(tk)-1)，k∈[n]生成1個(gè)數(shù)據(jù)點(diǎn)xk。因此，得到一些獨(dú)立非同分布的樣本{xk}n k=1，即每個(gè)樣本服從各自的分布。同樣地，與上述類似把它們轉(zhuǎn)化為獨(dú)立非同分布的樣本{yk}n k=1，即它們服從NPN(0,(Ωk)-1,f)。

6.2 評價(jià)標(biāo)準(zhǔn)

本節(jié)給出對圖模型結(jié)構(gòu)恢復(fù)的真陽性率（TPR）與假陽性率（FPR）來評價(jià)所有模型的性能。假定為已估計(jì)出的稀疏精度矩陣為真實(shí)的精度矩陣，給出指標(biāo)TPR與FPR的定義如下：

其中，為指標(biāo)函數(shù)。同時(shí)，為了驗(yàn)證本文模型的魯棒性，考慮對這些人工數(shù)據(jù)加一些噪聲。具體地，在每個(gè)樣本矩陣隨機(jī)選取[nr]個(gè)元素用5或-5代替，其中0≤r≤1為噪聲率。

6.3 實(shí)驗(yàn)結(jié)果

在實(shí)驗(yàn)中，利用4個(gè)半?yún)?shù)概率圖模型的學(xué)習(xí)作為評估模型效果，即K=4。同時(shí)，對于動態(tài)半?yún)?shù)概率圖模型的學(xué)習(xí)在n個(gè)時(shí)間點(diǎn)隨機(jī)選取4個(gè)時(shí)間點(diǎn)聯(lián)合估計(jì)作為評估模型效果。

由圖1可知，在學(xué)習(xí)非正態(tài)分布異構(gòu)數(shù)據(jù)的相關(guān)性時(shí)，本文的JSPGL優(yōu)于JGL與SPGL，也更加魯棒。由圖2可知，本文聯(lián)合模型在小樣本情況下依然優(yōu)于JGL與SPGL。由圖3可知，在學(xué)習(xí)動態(tài)的非正態(tài)分布異構(gòu)數(shù)據(jù)的相關(guān)性時(shí)，本文的JDSPG優(yōu)于DSPGL與DGL，也更加魯棒。同樣地，由圖4可知，本文聯(lián)合動態(tài)圖模型在小樣本情況下依然優(yōu)于其他方法。從圖3、圖4可知，JDSPGL并非很顯著地優(yōu)于DSPGL。由于這兩種方法估計(jì)相關(guān)性矩陣Σ(t)均用核光滑方法，它們在估計(jì)相關(guān)性矩陣時(shí)已經(jīng)把每個(gè)時(shí)間點(diǎn)的信息考慮進(jìn)去了，即已經(jīng)用了聯(lián)合學(xué)習(xí)思想。

Fig.1 ROC curves of estimating multiple SPGMs at different noise contamination levels(n=200 and p=200)圖1 多個(gè)半?yún)?shù)圖模型在不同程度噪聲污染下估計(jì)的ROC曲線（n=200與p=200）

Fig.2 ROC curves of estimating multiple SPGMs at different noise contamination levels(n=100 and p=200)圖2 多個(gè)半?yún)?shù)圖模型在不同程度噪聲污染下估計(jì)的ROC曲線（n=100與p=200）

Fig.3 ROC curves of estimating DSPGMs at different noise contamination levels(n=200 and p=200)圖3 動態(tài)半?yún)?shù)圖模型在不同程度噪聲污染下估計(jì)的ROC曲線（n=200與p=200）

Fig.4 ROC curves of estimating DSPGMs at different noise contamination levels(n=100 and p=200)圖4 動態(tài)半?yún)?shù)圖模型在不同程度噪聲污染下估計(jì)的ROC曲線（n=100與p=200）

7 真實(shí)數(shù)據(jù)實(shí)驗(yàn)

本文利用真實(shí)的腦影像數(shù)據(jù)與股票交易數(shù)據(jù)分別驗(yàn)證提出的聯(lián)合半?yún)?shù)圖模型（JSPGL）與聯(lián)合動態(tài)半?yún)?shù)圖模型（JDSPGL）的有效性。

腦影像數(shù)據(jù)（http://adni.loni.ucla.edu/）采集于32個(gè)老年癡呆（Alzheimer’s disease，AD）大腦、71個(gè)認(rèn)知障礙（mild cognitive impairment，MCI）大腦與62個(gè)正常（normal control，NC）大腦，且所有數(shù)據(jù)包括116個(gè)特征，每個(gè)特征代表每個(gè)解剖興趣區(qū)域。對于腦影像數(shù)據(jù)，將利用JSPGL聯(lián)合構(gòu)建三類大腦網(wǎng)絡(luò)，其為AD腦網(wǎng)絡(luò)、MCI腦網(wǎng)絡(luò)與NC腦網(wǎng)絡(luò)。通過估計(jì)這些腦網(wǎng)絡(luò)找到它們的共性與差異（見圖5）。

股票交易數(shù)據(jù)（http://finance.yahoo.com/）收集于標(biāo)準(zhǔn)普爾500指數(shù)公司從2003年1月到2008年1月每天股票交易數(shù)據(jù)。該數(shù)據(jù)包括452家公司的1 258條收盤價(jià)格?？紤]到該股票交易數(shù)據(jù)隨著時(shí)間較光滑變化，本文利用JDSPGL學(xué)習(xí)這452家公司在股票交易中動態(tài)的條件相關(guān)性。

這些真實(shí)數(shù)據(jù)沒有已知的結(jié)構(gòu)信息，因此本文類似于文獻(xiàn)[7]利用Kullback-Leible（KL）損失定量地驗(yàn)證模型估計(jì)的性能。對于多類數(shù)據(jù)如腦影響數(shù)據(jù)，首先把每類數(shù)據(jù)[nk]劃分為M份{D1,D2,…,DM}，然后定義KL-loss如下：

其中，是在訓(xùn)練樣本（[nk]減去Dm）上估計(jì)得到的；Sm為測試樣本Dm的樣本協(xié)方差矩陣。對于動態(tài)數(shù)據(jù)如股票交易數(shù)據(jù)，首先把所有數(shù)據(jù)[n]劃分為{D1,D2,…,DM}，然后定義KL-loss如下：

其中是在訓(xùn)練樣本（[n]減去Dm）上估計(jì)得到的。

由表1可知，在腦影像數(shù)據(jù)實(shí)驗(yàn)上，本文JSPGL的性能優(yōu)于SPGL與JGL。同時(shí)，由圖5可知，NC腦網(wǎng)絡(luò)與MCI腦網(wǎng)絡(luò)的差異要小于NC腦網(wǎng)絡(luò)與AD腦網(wǎng)絡(luò)，因此JSPGL學(xué)習(xí)得到的腦網(wǎng)絡(luò)同時(shí)具有較好的解釋性。由表2可知，在股票交易數(shù)據(jù)上，本文JDSPGL的性能優(yōu)于DSPGL與DGL。

Table 1 5-flod KL-loss on brain imaging dataset表1 圖模型在腦影像數(shù)據(jù)上的5重KL-loss

Table 2 5-flod KL-loss on stock trading dataset表2 圖模型在股票數(shù)據(jù)上的5重KL-loss

Fig.5 Brain networks estimated by joint semi-parameter graphical model圖5 聯(lián)合半?yún)?shù)圖模型估計(jì)的腦網(wǎng)絡(luò)

8 總結(jié)

本文提出了聯(lián)合半?yún)?shù)概率圖模型用于學(xué)習(xí)非正態(tài)分布異構(gòu)數(shù)據(jù)的條件相關(guān)性。同時(shí)，針對光滑變化的異構(gòu)數(shù)據(jù)，提出了聯(lián)合動態(tài)半?yún)?shù)圖模型。將基于非參排序的相關(guān)矩陣估計(jì)方法與結(jié)構(gòu)融合圖套索方法相結(jié)合，提出了一類半?yún)?shù)融合圖套索方法來估計(jì)提出的模型。特別針對動態(tài)半?yún)?shù)圖模型，提出了一種新的核光滑Kendall’s tau相關(guān)矩陣。由于放寬了正態(tài)分布的假設(shè)，使得本文模型比當(dāng)前聯(lián)合高斯圖模型更靈活。由于采用了基于非參排序的相關(guān)矩陣估計(jì)方法，使得本文模型更魯棒。在未來工作中，將提出的聯(lián)合動態(tài)圖模型推廣到混合變量的半?yún)?shù)圖模型[30]。

：

[1]Lauritzen S L.Graphical models[M].Oxford:Oxford University Press,1996.

[2]Liu Jianwei,Cui Lipeng,Luo Xionglin.Survey on the sparse learning of probabilistic graphical model[J].Chinese Journal of Computers,2016,39(8):1597-1611.

[3]Banerjee O,Ghaoui L E,d'Aspremont A.Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[J].Journal of Machine Learning Research,2008,9:485-516.

[4]Huang Shuai,Li Jing,Sun Liang,et al.Learning brain connectivity of Alzheimer’s disease from neuroimaging data[C]//Proceedings of the 23rd Annual Conference on Neural Information Processing Systems,Vancouver,Dec 7-10,2009.Red Hook:CurranAssociates,2009:808-816.

[5]Meinshausen N,Bühlmann P.High-dimensional graphs and variable selection with the Lasso[J].The Annals of Statistics,2006,34(3):1436-1462.

[6]Yuan Ming,Lin Yi.Model selection and estimation in the Gaussian graphical model[J].Biometrika,2007,94(1):19-35.

[7]Friedman J,Hastie T,Tibshirani R.Sparse inverse covariance estimation with the graphical Lasso[J].Biostatistics,2008,9(3):432-441.

[8]Yuan Ming.High dimensional inverse covariance matrix estimation via linear programming[J].Journal of Machine Learning Research,2010,11:2261-2286.

[9]Cai T,Liu Weidong,Luo Xi.A constrained?1 minimization approach to sparse precision matrix estimation[J].Journal of the American Statistical Association,2011,106(494):594-607.

[10]Tibshirani R.Regression shrinkage and selection via the Lasso[J].Journal of the Royal Statistical Society:Series B Methodological,1996,58(1):267-288.

[11]Liu Han,Lafferty J D,Wasserman LA.The nonparanormal:semiparametric estimation of high dimensional undirected graphs[J].Journal of Machine Learning Research,2009,10(3):2295-2328.

[12]Liu Han,Han Fang,Yuan Ming,et al.High-dimensional semiparametric Gaussian copula graphical models[J].The Annals of Statistics,2012,40(4):2293-2326.

[13]Xue Lingzhou,Zou Hui.Regularized rank-based estimation of high-dimensional nonparanormal graphical models[J].The Annals of Statistics,2012,40(5):2541-2571.

[14]Yang Sen,Lu Zhaosong,Shen Xiaotong,et al.Fused multiple graphical Lasso[J].SIAM Journal on Optimization,2015,25(2):916-943.

[15]Guo Jian,Levina E,Michailidis G,et al.Joint estimation of multiple graphical models[J].Biometrika,2011,89(1):1-15.

[16]Danaher P,Wang Pei,Witten D M.The joint graphical Lasso for inverse covariance estimation across multiple classes[J].Journal of the Royal Statistical Society:Series B Statistical Methodology,2014,76(2):373-397.

[17]Zhu Yunzhang,Shen Xiaotong,Pan Wei.Structural pursuit over multiple undirected graphs[J].Journal of the American StatisticalAssociation,2014,109(508):1683-1696.

[18]Lee W,Liu Yufeng.Joint estimation of multiple precision matrices with common structures[J].Journal of Machine Learning Research,2015,16:1035-1062.

[19]Cai T,Li Hongzhe,Liu Weidong,et al.Joint estimation of multiple high-dimensional precision matrices[J].Statistica Sinica,2016,26(2):445-464.

[20]Huang Feihu,Chen Songcan.Joint learning of multiple sparse matrix Gaussian graphical models[J].IEEE Transactions on Neural Networks and Learning Systems,2015,26(11):2606-2620.

[21]Yuan Ming,Lin Yi.Model selection and estimation in regression with grouped variables[J].Journal of the Royal Statistical Society:Series B Statistical Methodology,2006,68(1):49-67.

[22]Tibshirani R,Saunders M,Rosset S,et al.Sparsity and smoothness via the fused Lasso[J].Journal of the Royal Statistical Society:Series B Statistical Methodology,2005,67(1):91-108.

[23]Zhou Shuheng,Lafferty J D,Wasserman LA.Time varying undirected graphs[J].Machine Learning,2010,80(2):295-319.

[24]Kolar M,Xing E P.On time varying undirected graphs[C]//Proceedings of the 14th International Conference on Artificial Intelligence and Statistics,Fort Lauderdale,Apr 11-13,2011:407-415.

[25]Qiu Huitong,Han Fang,Liu Han,et al.Joint estimation of multiple graphical models from high dimensional time series[J].Journal of the Royal Statistical Society:Series B Statistical Methodology,2016,78(2):487-504.

[26]Fang Hongbin,Fang Kaitai,Kotz S.The meta-elliptical distributions with given marginal[J].Journal of Multivariate Analysis,2002,82(1):1-16.

[27]Kruskal W H.Ordinal measures of association[J].Journal of the American Statistical Association,1958,53(284):814-861.

[28]Boyd S,Parikh N,Chu E,et al.Distributed optimization and statistical learning via the alternating direction method of multipliers[J].Foundations&Trends in Machine Learning,2011,3(1):1-122.

[29]Hoefling H.A path algorithm for the fused Lasso signal approximator[J].Journal of Computational and Graphical Statistics,2010,19(4):984-1006.

[30]Fan Jianping,Liu Han,Ning Yang,et al.High dimensional semiparametric latent graphical model for mixed data[J].Journal of the Royal Statistical Society:Series B Statistical Methodology,2017,79(2):405-421.

附中文參考文獻(xiàn):

[2]劉建偉,崔立鵬,羅雄麟.概率圖模型的稀疏化學(xué)習(xí)[J].計(jì)算機(jī)學(xué)報(bào),2016,39(8):1597-1611.