蘇溫慶, 郭 驍, 張 海

(西北大學(xué)數(shù)學(xué)學(xué)院，西安 710127)

0 引言

復(fù)雜網(wǎng)絡(luò)是近年來醫(yī)學(xué)、計(jì)算機(jī)科學(xué)和統(tǒng)計(jì)學(xué)等領(lǐng)域研究的熱點(diǎn)問題之一。Albert 和Barab′asi[1]研究發(fā)現(xiàn)，不同于隨機(jī)網(wǎng)絡(luò)，在基因網(wǎng)絡(luò)、社交網(wǎng)絡(luò)和社會(huì)網(wǎng)絡(luò)等大多數(shù)真實(shí)網(wǎng)絡(luò)中，網(wǎng)絡(luò)節(jié)點(diǎn)的度服從冪律分布，即網(wǎng)絡(luò)中節(jié)點(diǎn)的度為d的概率正比于d?γ，這里γ> 0。這類網(wǎng)絡(luò)稱為無標(biāo)度網(wǎng)絡(luò)。無標(biāo)度網(wǎng)絡(luò)的典型特征是網(wǎng)絡(luò)中存在少數(shù)中心點(diǎn)(hubs)擁有極其多的邊連接，而大多數(shù)節(jié)點(diǎn)只有很少量的邊連接。在實(shí)際中，少數(shù)中心點(diǎn)往往對無標(biāo)度網(wǎng)絡(luò)起著主導(dǎo)作用。因此，在無標(biāo)度的假設(shè)下研究網(wǎng)絡(luò)的結(jié)構(gòu)，從而理解真實(shí)網(wǎng)絡(luò)的功能和特征不僅具有理論意義，更具有實(shí)際價(jià)值。

圖模型[2]作為一種分析網(wǎng)絡(luò)結(jié)構(gòu)的方法被廣泛使用。圖模型的基本思想在于通過圖表示隨機(jī)變量聯(lián)合分布，并基于分布研究變量之間相互關(guān)系。圖中的節(jié)點(diǎn)表示隨機(jī)變量，邊度量隨機(jī)變量間的關(guān)系。根據(jù)邊是否有方向，圖模型可分為有向圖模型與無向圖模型。本文基于有向圖模型開展真實(shí)網(wǎng)絡(luò)結(jié)構(gòu)分析。有向圖模型可度量隨機(jī)變量間的因果關(guān)系，兩變量間的箭頭Xi →Xj表示Xi是Xj的直接原因。因果關(guān)系在生物科學(xué)、醫(yī)藥學(xué)和人工智能等領(lǐng)域廣泛存在。因此，研究有向圖模型具有重要意義。有向無環(huán)圖(Directed Acyclic Graphs, DAGs)是指邊有向，但不存在有向環(huán)的圖，作為有向圖模型的一個(gè)特類，近年來受到廣泛關(guān)注。相比于無向圖模型，DAGs 的估計(jì)是一個(gè)困難的問題，主要難點(diǎn)在于DAGs 的結(jié)構(gòu)學(xué)習(xí)往往對應(yīng)于一個(gè)NP 難問題，且由于觀測等價(jià)性，可能無法估計(jì)邊的方向。大多數(shù)現(xiàn)有估計(jì)DAGs 的方法可大致歸為兩類。第一類是基于約束的方法。Spirtes 等[3]提出的PC 算法和Tsamardinos 等[4]提出的MMHC 算法均為此類方法。這類方法通過一系列條件獨(dú)立性的檢測來判斷節(jié)點(diǎn)之間是否存在邊，主要缺點(diǎn)在于實(shí)際問題往往難以滿足上述方法的假設(shè)。第二類是基于準(zhǔn)則的方法，即通過最優(yōu)化給定的評分函數(shù)來估計(jì)DAGs。Lam 和Bacchus[5]通過最小描述長度構(gòu)造評分函數(shù)，Heckerman 等[6]通過貝葉斯方法來估計(jì)DAGs。網(wǎng)絡(luò)的一個(gè)主要特征在于實(shí)際中存在大量稀疏網(wǎng)絡(luò)，所謂稀疏網(wǎng)絡(luò)是指網(wǎng)絡(luò)中存在的邊是稀疏的。因此，目前發(fā)展起來的稀疏正則化方法越來越多的被應(yīng)用于網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)建模。Shaojaie 和Michailidis[7]在假定節(jié)點(diǎn)序已知的條件下，提出通過罰對數(shù)似然估計(jì)方法直接估計(jì)高斯DAGs 的鄰接矩陣。Fu 和Zhou[8]在假定節(jié)點(diǎn)序未知的條件下，通過l1正則化方法研究稀疏DAGs 的結(jié)構(gòu)學(xué)習(xí)。隨后，Aragam 和Zhou[9]將該方法進(jìn)一步推廣到非凸懲罰下。文獻(xiàn)[10]分析了高斯分布下高維稀疏DAGs 的l0懲罰估計(jì)的理論性質(zhì)。

上述正則化方法大多通過l1型懲罰函數(shù)來引入稀疏先驗(yàn)。但l1型罰項(xiàng)并不能反映網(wǎng)絡(luò)的無標(biāo)度特征。近來，部分學(xué)者通過引入無標(biāo)度先驗(yàn)，對無向圖的結(jié)構(gòu)學(xué)習(xí)開展研究。Liu 和Ihler[11]直接利用網(wǎng)絡(luò)節(jié)點(diǎn)的度服從冪律分布這一先驗(yàn)信息，提出了罰項(xiàng)為Log 型與l1型懲罰函數(shù)復(fù)合的模型，用l1范數(shù)近似節(jié)點(diǎn)的度，并提出用重賦權(quán)迭代算法來求解這一模型。進(jìn)一步，Guo 等[12]用lq(0

本文在假定節(jié)點(diǎn)序已知的條件下，研究具有無標(biāo)度特征DAGs 的結(jié)構(gòu)學(xué)習(xí)問題。通過引入網(wǎng)絡(luò)中節(jié)點(diǎn)度的信息和邊的稀疏先驗(yàn)，構(gòu)造罰項(xiàng)為Log 型與lq(0

1 有向無環(huán)圖結(jié)構(gòu)學(xué)習(xí)模型

本節(jié)討論具有無標(biāo)度特征的稀疏有向無環(huán)圖結(jié)構(gòu)學(xué)習(xí)模型。首先我們給出有向無環(huán)圖模型的數(shù)學(xué)表示，然后基于無標(biāo)度信息，構(gòu)建新的稀疏有向無環(huán)圖結(jié)構(gòu)學(xué)習(xí)模型。

有向無環(huán)圖是圖中只包含有向邊，不包含有向環(huán)的一類圖。考慮一個(gè)DAG，G=(V,E)，這里V= 1,2,··· ,p對應(yīng)于圖的p個(gè)節(jié)點(diǎn)，且節(jié)點(diǎn)j(j= 1,2,··· ,p)表示隨機(jī)變量Xj，E ?V×V代表有向邊。邊(i,j)∈E可表示為i →j。此時(shí)節(jié)點(diǎn)i稱為節(jié)點(diǎn)j的父代節(jié)點(diǎn)，記為i ∈pa(j)。在本文中，我們主要關(guān)注高斯有向無環(huán)圖模型。DAGs 的p個(gè)節(jié)點(diǎn)對應(yīng)的隨機(jī)變量Xi,··· ,Xp，服從均值為0，方差為Σ0的p維高斯分布，即Xi,··· ,Xp ～N(0,Σ0)。

我們稱(B0,?0)是對應(yīng)于Σ0的DAG。對于一個(gè)DAG 對應(yīng)的鄰接矩陣，由于圖中不存在有向環(huán)，因此，可通過置換把鄰接矩陣變?yōu)橄氯蔷仃?。在本文中，我們關(guān)注節(jié)點(diǎn)序已知的情形，因此，同樣可通過提前置換把B0變?yōu)橹鲗蔷€元素為0 的下三角矩陣。不失一般性，我們以下假定B0為下三角矩陣。

其中λ> 0 是調(diào)控參數(shù)，該參數(shù)中包含了尺度參數(shù)γ。顯然，該模型可劃分為p個(gè)獨(dú)立的子模型，其中第j個(gè)子模型對應(yīng)B的第j列。

注1 在模型中添加大于零的常數(shù)?以保證log 函數(shù)有意義。求解過程中可將?視為調(diào)整參數(shù)，以得到模型的全局最優(yōu)解或局部最優(yōu)解，但通過實(shí)驗(yàn)發(fā)現(xiàn)該方法對?不敏感。因此，本文預(yù)先設(shè)定?為一個(gè)較小的常數(shù)。

2 算法及其收斂性分析

1: 令k =0，輸入τ >0, ?>0 及初始值βt(0)j =βtj。2: 對于k =0,1,2,···，重復(fù)下式直至收斂或達(dá)到要求的迭代次數(shù)βt(k+1)j =arg min βj Rτ(βj;βt j,βt(k)j )=arg min p∑1//Xj ?Xβj//2 2+λ2ωt j2 q βj i=1//βtj//qq+?(βt(k)ij +τ)1?q|βij|.

3: 令βt+1 j 等于步2 輸出結(jié)果。

3 實(shí)驗(yàn)

本節(jié)通過數(shù)值模擬和實(shí)際數(shù)據(jù)來說明我們方法有效性。在模型(6)中，考慮到lq(0

數(shù)值實(shí)驗(yàn)中數(shù)據(jù)的生成方式如下。首先，我們用Albert 和Barab′asi[1]提出的擇優(yōu)連接機(jī)制來生成具有無標(biāo)度特征的有向無環(huán)網(wǎng)絡(luò)。具體地，給定一個(gè)初始節(jié)點(diǎn)，對下一個(gè)新的節(jié)點(diǎn)，以正比于的概率選擇某個(gè)已有節(jié)點(diǎn)，并生成由新節(jié)點(diǎn)指向舊節(jié)點(diǎn)的邊，這里di為已有節(jié)點(diǎn)i的入度，γ為尺度參數(shù)。其次，對于鄰接矩陣B，若節(jié)點(diǎn)i是節(jié)點(diǎn)j的父代節(jié)點(diǎn)，則令βij為大于0 的常數(shù)β，否則令βij為0。對于協(xié)方差矩陣?，數(shù)值實(shí)驗(yàn)中令ωj為1。最后，通過線性結(jié)構(gòu)方程模型(2)來生產(chǎn)實(shí)驗(yàn)數(shù)據(jù)X。

我們通過ROC 曲線來度量各方法的優(yōu)劣。該曲線反應(yīng)了調(diào)控參數(shù)λ變化時(shí)各方法估計(jì)真實(shí)網(wǎng)絡(luò)的整體表現(xiàn)。ROC 曲線的橫軸是FPR(False Positive Rate)，縱軸是TPR(True Positive Rate)，這里FPR=FP/(FP+TN), TPR=TP/(TP+FN)，其中TP(True Positive)為真實(shí)網(wǎng)絡(luò)有邊，模型估計(jì)為有邊的個(gè)數(shù)；FP(False Positive)為真實(shí)網(wǎng)絡(luò)無邊，模型估計(jì)為有邊的個(gè)數(shù)；FN(False Negative)為真實(shí)網(wǎng)絡(luò)有邊，模型估計(jì)為無邊的個(gè)數(shù)；TN(True Negative)為真實(shí)網(wǎng)絡(luò)無邊，模型估計(jì)為無邊的個(gè)數(shù)。

實(shí)驗(yàn)1 我們分別比較網(wǎng)絡(luò)維度p，信號(hào)強(qiáng)度β，尺度參數(shù)γ不同時(shí)三種方法的表現(xiàn)。首先，給定尺度參數(shù)γ= 1，生成維度分別為p= 50,100,200 的網(wǎng)絡(luò)，圖1 為具體網(wǎng)絡(luò)。在不同網(wǎng)絡(luò)維度下，分別比較信號(hào)強(qiáng)度為β=2,1,0.3 時(shí)各方法表現(xiàn)。這里樣本大小固定為n=100，?和τ均固定為0.1。圖2 為實(shí)驗(yàn)結(jié)果。其次，討論當(dāng)網(wǎng)絡(luò)中hubs 的影響增大，即γ變大時(shí)三種方法的表現(xiàn)。給定網(wǎng)絡(luò)維度p= 50，樣本大小n= 100，信號(hào)強(qiáng)度β=1，?和τ均固定為0.1，分別給定γ=1,1.5,2，圖3 為實(shí)驗(yàn)結(jié)果。每條ROC 曲線均為重復(fù)20 次實(shí)驗(yàn)結(jié)果(以下數(shù)值實(shí)驗(yàn)均為重復(fù)20 次結(jié)果)。其中，RWq=0.5 表示本文提出的模型，我們選擇q= 0.5 為本文模型的代表。RWq= 1 表示罰項(xiàng)為Log 型與l1型懲罰函數(shù)的復(fù)合模型，通過重賦權(quán)迭代算法求解。L1-CD 表示Shojaie 和Michailidis[7]提出的l1懲罰對數(shù)似然估計(jì)方法，通過坐標(biāo)下降算法求解。

圖1 實(shí)驗(yàn)生成有向網(wǎng)絡(luò)，維度分別為50、100、200

從圖2 可以看出，幾乎在所有情況下，本文方法RWq= 0.5 的表現(xiàn)好于其它兩種方法。當(dāng)信號(hào)強(qiáng)度β固定而網(wǎng)絡(luò)維度p變化時(shí)，隨著維度增大，三種方法的表現(xiàn)基本上均逐漸變差。當(dāng)p= 50 時(shí)，RWq= 1 的表現(xiàn)稍遜于RWq= 0.5，但隨著維度增大，除β= 0.3 外，RWq= 1 的表現(xiàn)與RWq= 0.5 的差距逐漸變大。當(dāng)信號(hào)強(qiáng)度低(β= 0.3)時(shí)，三種方法的表現(xiàn)均變差，此時(shí)本文方法的表現(xiàn)仍好于其它兩種方法但優(yōu)勢不明顯。從圖3 可以看出，在不同尺度參數(shù)γ下，本文方法的表現(xiàn)均好于其它兩種方法。隨著γ增大，RWq= 0.5 和RWq= 1 的表現(xiàn)均變好，L1-CD 的表現(xiàn)變差。這說明本文方法對重建由中心節(jié)點(diǎn)主導(dǎo)的網(wǎng)絡(luò)上是有效的。

圖2 不同(p,β)下ROC 曲線

圖3 不同尺度參數(shù)γ 下的ROC 曲線

實(shí)驗(yàn)2 實(shí)驗(yàn)1 中均通過擇優(yōu)連接機(jī)制來生成具有無標(biāo)度特征的有向無環(huán)網(wǎng)絡(luò)，且均為節(jié)點(diǎn)入度服從冪律分布的網(wǎng)絡(luò)。為說明我們所提的模型對重構(gòu)節(jié)點(diǎn)出度服從冪律分布的網(wǎng)絡(luò)上仍有效，我們選取兩個(gè)節(jié)點(diǎn)出度服從冪律的真實(shí)網(wǎng)絡(luò)作為實(shí)驗(yàn)網(wǎng)絡(luò)，并通過線性結(jié)構(gòu)方程模型(2)生成實(shí)驗(yàn)數(shù)據(jù)X，即使用真實(shí)網(wǎng)絡(luò)但不使用真實(shí)數(shù)據(jù)。第一個(gè)網(wǎng)絡(luò)為大腸桿菌(escherichia coli, ecoli)轉(zhuǎn)錄調(diào)節(jié)網(wǎng)絡(luò)，Kao 等[22]提出利用網(wǎng)絡(luò)結(jié)構(gòu)分析來推斷大腸桿菌的轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)。他們提供了7 個(gè)轉(zhuǎn)錄因子和40 個(gè)受調(diào)控基因的大腸桿菌調(diào)節(jié)網(wǎng)絡(luò)，即網(wǎng)絡(luò)維度為p= 47，見圖4(a)。第二個(gè)網(wǎng)絡(luò)為DREAM3(Dialogue on Reverse Engineering Assessment and Methods)數(shù)據(jù)集[23]中的ecoli1 網(wǎng)絡(luò)，網(wǎng)絡(luò)維度為p=50，見圖4(b)。顯然ecoli 網(wǎng)絡(luò)中節(jié)點(diǎn)出度服從冪律分布，而ecoli1 網(wǎng)絡(luò)中同時(shí)存在高入度的中心點(diǎn)和高出度的中心點(diǎn)。實(shí)驗(yàn)中給定樣本大小n= 30，?和τ均固定為0.1。圖4(c)和圖4(d)為對應(yīng)ROC 曲線。本文方法的表現(xiàn)均好于其它兩種方法。這說明本文方法在網(wǎng)絡(luò)節(jié)點(diǎn)出度服從冪律分布時(shí)同樣是有效的。

圖4 實(shí)驗(yàn)網(wǎng)絡(luò)和對應(yīng)的ROC 曲線

實(shí)驗(yàn)3 我們通過真實(shí)數(shù)據(jù)驗(yàn)證所提方法的有效性。我們研究ecoli 網(wǎng)絡(luò)，實(shí)驗(yàn)中使用Kao 等[22]提供的23 個(gè)時(shí)間點(diǎn)上7 個(gè)轉(zhuǎn)錄因子和40 個(gè)受調(diào)控基因的基因表達(dá)數(shù)據(jù)，即樣本個(gè)數(shù)為23。圖5(a)為已知網(wǎng)絡(luò)(包含基因名稱)。圖5(b)為實(shí)驗(yàn)結(jié)果?？梢钥闯霰疚姆椒ǖ谋憩F(xiàn)優(yōu)于其它兩種方法。這與數(shù)值實(shí)驗(yàn)中結(jié)果一致。

圖5 ecoli 調(diào)節(jié)網(wǎng)絡(luò)及三種方法的ROC 曲線

以上實(shí)驗(yàn)中我們固定?和τ均為0.1?，F(xiàn)在我們討論?和τ的取值對本文方法結(jié)果的影響。對于?，給定實(shí)驗(yàn)網(wǎng)絡(luò)為圖1(a)，樣本大小n=100，信號(hào)強(qiáng)度β=1, τ=0.1，分別給定?= 0.1,0.01,0.001。對于τ，給定實(shí)驗(yàn)網(wǎng)絡(luò)為圖1(a)，樣本大小n= 100，信號(hào)強(qiáng)度β= 1, ?= 0.1，分別給定τ= 0.1,0.01,0.001。圖6 為實(shí)驗(yàn)結(jié)果。從圖6 可以看出我們的方法對?和τ的取值不敏感，因此預(yù)先固定?和τ是合理的。

圖6 不同? 和τ 下的ROC 曲線

注4 考慮到運(yùn)算效律，在實(shí)驗(yàn)中均采用一步重賦權(quán)迭代。實(shí)驗(yàn)表明一步重賦權(quán)迭代求解效果明顯。

4 結(jié)論

本文研究在無標(biāo)度先驗(yàn)下，節(jié)點(diǎn)序已知的有向無環(huán)圖結(jié)構(gòu)學(xué)習(xí)問題。通過引入網(wǎng)絡(luò)中節(jié)點(diǎn)度的信息和邊的稀疏先驗(yàn)，構(gòu)造罰項(xiàng)為Log 型與lq(0