基于差分進化策略的貝葉斯網絡結構學習方法

牛艷飛，馬 潔

(北京信息科技大學自動化學院，北京 100192)

1 引言

貝葉斯網絡(Bayesian Network，BN)經常用于處理涉及不確定知識的問題，是一種綜合了圖論和概率論的圖模型，既可以定性地表示節點之間的因果關系，又可以定量地描述變量間的聯合概率分布。在數據挖掘領域，貝葉斯網絡一直是眾多學者的研究重點，并且在各領域均有重要的應用，如故障溯源[1]、醫療診斷[2]、金融分析[3]等。

貝葉斯網絡的構建主要由結構學習和參數學習兩部分組成，其中結構學習因其對先驗知識的依賴而成為貝葉斯網絡構建的一大難點[4]。本文針對貝葉斯網絡的結構學習，提出一種混合部分互信息量和改進差分進化算法的PMI-DE算法，該算法首先基于變量之間的部分互信息描述網絡節點的關聯度并構建初始種群，再通過引入動態因子改進差分進化算法的交叉策略，平衡算法的全局收斂和局部搜索能力，從而得到最優的貝葉斯網絡結構。最后在仿真中將該算法與遺傳算法和爬山算法進行對比，分析該算法的優劣。

2 貝葉斯網絡

2.1 貝葉斯網絡的概念

貝葉斯網絡[5-6]由變量之間的有向邊和節點概率表構成，其數學描述可以用BN=來表示。

其中G=是一個二進制矩陣，用來表示貝葉斯網絡的結構。X={x1，x2，…，xn}決定矩陣的大小，代表貝葉斯網絡中的各個節點變量，E={a11…aij…ann}(aij=0or1)決定矩陣中的元素，反映了網絡節點xi和xj之間的因果關系。圖1為一個貝葉斯網絡結構的矩陣表達。

圖1 貝葉斯網絡結構示意圖

P是貝葉斯網絡中各個節點xi∈X的聯合概率分布P(xi|Pa(xi))，用來定量描述各個節點之間的依賴關系，其中Pa(xi)為xi在網絡結構G中的父節點集合。貝葉斯網絡全節點之間的聯合概率分布如式(1)

(1)

2.2 貝葉斯網絡的結構學習

從數據集D={d1，d2，d3…dn}中尋找一個與變量關系匹配度最高的貝葉斯網絡結構是結構學習的主要目標，但是該過程的計算復雜度會隨著節點的個數增多而成為一個無法進行精確求解的NP-hard問題[7]。式(2)展示了變量節點個數與貝葉斯網絡結構個數之間的關系

(2)

可見想要通過遍歷的方式對貝葉斯網絡結構進行精確查找是不可行的，對此，國內外已經有很多貝葉斯網絡結構尋優方法，比較主流的方法是通過合理的結構評價函數和啟發式搜索算法對網絡結構進行尋優，最終得到準確的網絡結構。Cooper等人[8]提出的K2算法利用節點序縮小了搜索空間并且可以避免產生非法結構，但是K2算法嚴重依賴于節點序的正確性；Tsamardinos等人[9]提出了最大最小爬山算法(Max-Min Hill-Climbing，MMHC)，先通過獨立性測試縮小搜索空間，再利用爬山算法進行遍歷確定貝葉斯網絡結構，該方法得到了比較廣泛的應用；Saoussen[10]等提出PSO-K2算法，通過粒子群算法對節點的順序進行尋優，再將最優結果與K2算法結合得到貝葉斯網絡的結構；劉浩然等[11]提出MAK(MWST-ACO-K2)算法，該算法融合最大支撐樹和蟻群算法搜索節點序，在處理小型網絡時可取得較理想的結果；魏中強等[12]提出了CMI-PK2算法，其先通過條件互信息得到節點序，再將K2算法的搜索策略進行改進來提高學習的速度，得到了較好的學習效果。

本文提出了一種基于部分互信息和差分進化算法的混合算法(PMI-DE算法)，該算法首先利用部分互信息度量網絡節點之間的相關性，并在此基礎上產生初始種群，然后利用差分進化算法對貝葉斯網絡結構進行尋優，通過引入動態交叉因子來改進交叉策略，使得算法能夠自適應地選擇交叉對象，最終得到最優貝葉斯網絡結構。

3 PMI-DE算法研究

3.1 部分互信息

互信息是信息論中的一個重要概念，經常被用來度量隨機變量之間的相關性，若已知兩個隨機變量X={x1，x2，x3…xn}Y={y1，y2，y3…ym}，則X和Y之間的互信息如式(3)

(3)

式中，p(x)表示X的先驗概率，p(x，y)為X和Y的聯合概率。

條件互信息用來衡量隨機變量X和Y在條件Z下的相關關系。定義如式(4)所示

(4)

式中，p(x|z)p(y|z)分別為X、Y在條件Z下的先驗概率，p(x，y|z)為X和Y在條件Z下的聯合概率，p(x，y，z)為X、Y、Z的聯合概率。

雖然互信息和條件互信息都能夠度量變量之間的相關性，但是互信息往往會過高的估計兩個隨機變量之間的相關性，以此為依據構建的貝葉斯網絡結構會產生過多的冗余邊，而條件互信息則會過低的估計兩個隨機變量之間的相關性，這將導致貝葉斯網絡出現嚴重的邊缺失現象，為了平衡互信息和條件互信息的這種缺陷，Zhao等人[13]提出了使用部分互信息(part mutual information，PMI)來度量隨機變量之間的相關性。

已知隨機變量X和Y在條件Z下獨立表示為式(5)

p(x|z)p(y|z)=p(x，y|z)

(5)

則X和Y在條件Z下的部分獨立性定義如式(6)

p*(x|z)p*(y|z)=p(x，y|z)

(6)

PMI(X，Y|Z)由式(7)所示的部分獨立性和KL距離(Kullback-Leibler divergence)[14]得出，定義如下

PMI(X，Y|Z)=

D(p(x，y，z)||p*(x|z)p*(y|z)p(z))

(7)

式中p(x，y，z)為X、Y、Z的聯合概率，D(p(x，y，z)||p*(x|z)p*(y|z)p(z))表示p(x，y，z)到p*(x|z)p*(y|z)p(z)的KL距離。PMI(X，Y|Z)也可以寫成式(8)的形式

(8)

3.2 差分進化算法改進

差分進化算法(Differential Evolution，DE)是由Storn等人[15]于1997年提出的，其主要包含變異、交叉、選擇三大部分。與遺傳算法的不同之處在于其變異操作是從種群中隨機挑選兩個個體，將差作為第三個個體的變化來源進行變異。該算法主要用于解決連續變量的全局尋優問題，但是貝葉斯網絡由離散的二進制矩陣表示，為了將差分進化算法應用于結構學習，將算法的變異和交叉算子進行了相應的修改。

3.2.1 適應度函數

適應度函數是演化算法的導向，用來衡量當前的迭代過程中最優的網絡結構是否滿足設定要求。評價貝葉斯網絡結構的函數根據原理可以劃分為兩類，一是基于貝葉斯統計的評分函數，如K2評分、BDeu評分；二是基于信息論的評分函數，如MIT評分[16]。本文采用可分解的BIC評分函數[17]度量網絡結構的優劣程度，其定義如式(9)

(9)

BIC((xi，pa(xi))|D)

(10)

則有

(11)

由式(10)和式(11)可見，BIC評分是可分解的，但是BIC評分一般為負值，為了便于比較和盡可能保留最終結果的有效數字，本文定義適應度函數如式(12)

(12)

3.2.2 初始種群的構建

初始種群的確定是算法的重要組成部分，為了挑選優秀的初始種群，本文首先通過對n個節點利用部分互信息進行打分，得到節點相關關系權重矩陣W={w11，w12，w13…wij…wnn}，其中wij的計算如式(13)

(13)

其中x∈xi，y∈xj表示節點xi和xj的類別，Z表示類別向量

Gi={g11，g12，g13…gij…gnn}

(14)

該過程首先利用節點之間的相關關系權重矩陣，限制了初始種群中冗余邊和缺失邊的數量，同時概率矩陣保證了初始種群的多樣性，為后續的變異和交叉操作提供了可靠的搜索基礎。

3.2.3 變異算子

與連續變量的尋優不同，在貝葉斯網絡結構的尋優中，個體是二進制矩陣，傳統的變異算子無法直接應用，故設計了針對于二進制矩陣的變異算子。

首先從初始種群中隨機選擇一個個體Gj={g11，g12，g13…gij…gnn}，j=1，…，N，將Gj變化成為當代種群的最優結構Gbest需要進行的操作用轉移矩陣Mov來表示

Mov=Gbest-Gj={m11，m12，m13…mij…mnn}

(15)

其中mij代表了轉移過程的三種操作

(16)

將待變異的個體按照轉移矩陣進行變異，并由變異概率cp進行控制：

(17)

變異過程中網絡結構中會產生無意義的gij=-1元素，其代表了該結構變異之前此處無邊，故不需要進行轉移矩陣中的刪除邊操作，將其置零，最終得到變異之后的個體Ui，i=1，…，N

3.2.4 交叉算子

交叉操作的目的是在個體的鄰域進行局部搜索，傳統差分進化算法的交叉操作全程由一個固定的交叉概率控制，在算法中期極易陷入局部最優，為了平衡算法的全局尋優和局部搜索能力，本文引入了自適應的動態交叉因子，其控制待交叉個體的交叉對象來源Hi。動態因子計算如下

a=t/Tmax

(18)

式中t為當前迭代次數，Tmax為最大迭代次數，每個個體的交叉對象選擇由a來控制

(19)

交叉過程采用兩點交叉方式，即從交叉對象中隨機選取兩個元素，與待交叉個體進行元素互換，得到交叉后的個體Vi，i=1，…，N。

(20)

3.3 非法結構修復

算法的變異和交叉操作可能給網絡結構帶來部分非法結構，例如雙向邊以及環狀結構。為了保證迭代過程中個體的正確性，必須在變異和交叉操作后對個體進行非法結構判斷和修復[18]，修正的步驟如圖所示：

圖2 網絡結構修復流程圖

3.4 算法偽代碼

PMI-DE算法的偽代碼如下：

PMI-DE Algorithm

Input：Data set、N、cp、hp、Tmax

Output：Gbest

Step1：Initialization

1)constructs the partial mutual information weight matrix W by Formula (13)

2)construct the edge existence probability matrix P

3)for i=1：N

4) Generate individuals Gi

5) calculate F(Gi)

6) end

7) Find the best individual Gbest(0)，(F(Gbest)=Max(F(Gi)))

8) Pbest=G(0)

Step2：evolution

9) While t < Tmax

10) for i=1：N

11) if rand

12) Mov=Gbest(t)-Gi(t)

13) Ui(t)=Gi(t)+Mov

14) illegal structure restoration

15) end

16) Calculate the dynamic factoraand Select crossover Hi(t) by Formula (18) and (19)

17) Two points cross to produce Vi(t) by Formula (20)

18) illegal structure restoration

Step3：update

19) if F(Vi(t))

20) Gi(t+1)=Vi(t)

21) else

22) Gi(t+1)=Gi(t)

23) end

26) end

27)end

4 仿真研究

為了驗證PMI-DE算法的有效性，本節以小型Asia網絡和中型Car網絡為模型，分別生成樣本數為500、1000、2000的數據集，與遺傳算法(GS)、爬山算法(HC)進行比較，分析PMI-DE算法的優缺點。

為了評價不同的算法在相同的數據集中所得到的貝葉斯網絡結構的優劣程度，本文采用海明距離來度量算法所得的網絡結構與真實結構的差異程度，海明距離的計算如式(21)

H(G)=A(G)+D(G)+R(G)

(21)

式中，A(G)表示冗余邊，D(G)表示缺失邊，R(G)表示反向邊。網絡結構的海明距離越小，說明該結構與真實的網絡結構相似度越高，進而說明了算法能夠更加準確地從數據集中學習到貝葉斯網絡結構。

實驗采用Matlab軟件平臺，運行環境：操作系統Windows 10，處理器為Intel(R) Core(TM) i7-8750H，主頻2.20GHz，內存8GB。PMI-DE算法的相關參數初始化設置如下：種群規模：50，變異概率為0.5，交叉概率為0.5，最大迭代次數為30次。為了保證算法結果的有效性，分別將三種算法獨立運行10次，并取平均值記錄，結果如下。

表1 Asia網絡樣本錯誤邊對比

表2 Car網絡樣本錯誤邊對比

圖3 Asia網絡結構學習漢明距離對比

圖4 Car網絡結構學習漢明距離對比

可見，在訓練數據相同的情況下，PMI-DE算法所得的網絡結構漢明距離最小，即得到的網絡結構最優，說明了PMI-DE算法在對網絡結構進行尋優時具有比較高的精度。

為了測試PMI-DE算法的學習性能，選擇Asia網絡樣本數量為500、1000、2000的數據集，三種算法分別運行10次取平均值，將PMI-DE算法的迭代速度和最佳評分值與其它兩種算法進行對比。結果見表3和圖4。

表3 樣本數為1000的數據集中學習性能比較

圖5 三種算法運行效率對比

由表3和圖5可以看出，當樣本數較小時，三種算法的計算速度相當，但是當樣本數較大時，遺傳算法的計算速度表現出明顯的不足，這是因為遺傳算法在迭代的過程中容易陷入局部最優，導致在要求精度的前提下需要更多的迭代次數，PMI-DE算法和爬山算法的計算時間相當，但是由最優評分值可以看出在精度上PMI-DE算法要更加精確。

5 結論

從數據中學習貝葉斯網絡的結構是一個棘手的問題，本文提出一種結合部分互信息和差分進化算法的混合學習方法，該算法以節點之間的部分互信息為基礎構建初始種群，縮小了搜索空間，提升了算法的執行效率，再通過引入動態因子改進傳統差分進化算法的交叉算子，使算法能夠自適應地平衡全局尋優和局部搜索，防止算法陷入局部最優。在標準的貝葉斯網絡中進行仿真，結果表明該算法有良好的尋優效果和運行效率，為貝葉斯網絡的結構學習提供了新的混合方法。