基于無標度先驗的有向無環圖結構學習

蘇溫慶, 郭 驍, 張 海

(西北大學數學學院，西安 710127)

0 引言

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

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

上述正則化方法大多通過l1型懲罰函數來引入稀疏先驗。但l1型罰項并不能反映網絡的無標度特征。近來，部分學者通過引入無標度先驗，對無向圖的結構學習開展研究。Liu 和Ihler[11]直接利用網絡節點的度服從冪律分布這一先驗信息，提出了罰項為Log 型與l1型懲罰函數復合的模型，用l1范數近似節點的度，并提出用重賦權迭代算法來求解這一模型。進一步，Guo 等[12]用lq(0

本文在假定節點序已知的條件下，研究具有無標度特征DAGs 的結構學習問題。通過引入網絡中節點度的信息和邊的稀疏先驗，構造罰項為Log 型與lq(0

1 有向無環圖結構學習模型

本節討論具有無標度特征的稀疏有向無環圖結構學習模型。首先我們給出有向無環圖模型的數學表示，然后基于無標度信息，構建新的稀疏有向無環圖結構學習模型。

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

我們稱(B0,?0)是對應于Σ0的DAG。對于一個DAG 對應的鄰接矩陣，由于圖中不存在有向環，因此，可通過置換把鄰接矩陣變為下三角矩陣。在本文中，我們關注節點序已知的情形，因此，同樣可通過提前置換把B0變為主對角線元素為0 的下三角矩陣。不失一般性，我們以下假定B0為下三角矩陣。

其中λ> 0 是調控參數，該參數中包含了尺度參數γ。顯然，該模型可劃分為p個獨立的子模型，其中第j個子模型對應B的第j列。

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

2 算法及其收斂性分析

1: 令k =0，輸入τ >0, ?>0 及初始值βt(0)j =βtj。2: 對于k =0,1,2,···，重復下式直至收斂或達到要求的迭代次數β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 輸出結果。

3 實驗

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

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

我們通過ROC 曲線來度量各方法的優劣。該曲線反應了調控參數λ變化時各方法估計真實網絡的整體表現。ROC 曲線的橫軸是FPR(False Positive Rate)，縱軸是TPR(True Positive Rate)，這里FPR=FP/(FP+TN), TPR=TP/(TP+FN)，其中TP(True Positive)為真實網絡有邊，模型估計為有邊的個數；FP(False Positive)為真實網絡無邊，模型估計為有邊的個數；FN(False Negative)為真實網絡有邊，模型估計為無邊的個數；TN(True Negative)為真實網絡無邊，模型估計為無邊的個數。

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

圖1 實驗生成有向網絡，維度分別為50、100、200

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

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

圖3 不同尺度參數γ 下的ROC 曲線

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

圖4 實驗網絡和對應的ROC 曲線

實驗3 我們通過真實數據驗證所提方法的有效性。我們研究ecoli 網絡，實驗中使用Kao 等[22]提供的23 個時間點上7 個轉錄因子和40 個受調控基因的基因表達數據，即樣本個數為23。圖5(a)為已知網絡(包含基因名稱)。圖5(b)為實驗結果。可以看出本文方法的表現優于其它兩種方法。這與數值實驗中結果一致。

圖5 ecoli 調節網絡及三種方法的ROC 曲線

以上實驗中我們固定?和τ均為0.1。現在我們討論?和τ的取值對本文方法結果的影響。對于?，給定實驗網絡為圖1(a)，樣本大小n=100，信號強度β=1, τ=0.1，分別給定?= 0.1,0.01,0.001。對于τ，給定實驗網絡為圖1(a)，樣本大小n= 100，信號強度β= 1, ?= 0.1，分別給定τ= 0.1,0.01,0.001。圖6 為實驗結果。從圖6 可以看出我們的方法對?和τ的取值不敏感，因此預先固定?和τ是合理的。

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

注4 考慮到運算效律，在實驗中均采用一步重賦權迭代。實驗表明一步重賦權迭代求解效果明顯。

4 結論

本文研究在無標度先驗下，節點序已知的有向無環圖結構學習問題。通過引入網絡中節點度的信息和邊的稀疏先驗，構造罰項為Log 型與lq(0