基于混合粒子群優化的貝葉斯網絡結構學習方法

尉永清，陳小雪，伊 靜，孟媛媛

1(山東師范大學 信息科學與工程學院，濟南 250358)2(山東警察學院 公共基礎部，濟南 250014)3(山東省分布式計算機軟件新技術重點實驗室，濟南 250358)4(山東建筑大學 信息科學與工程學院，濟南 250101)

1 引 言

貝葉斯網絡以概率論和圖論為理論基礎，它將概率統計與圖論相結合，主要應用于解決復雜系統不確定性數據推理和分析，已經成為處理不確定性問題的有效工具[1，2].目前貝葉斯網絡的研究是人工智能[3，4]、基因分析[5]、金融投資與分析[6]、故障診斷[7]、信息融合[8]等領域不確定知識表達和推理技術的主流方法.

貝葉斯網絡研究主要集中在貝葉斯網絡推理、貝葉斯網絡學習和基于貝葉斯網絡的應用三個方面.其中，貝葉斯網絡學習作為貝葉斯網絡的基礎，主要是通過分析數據來獲得貝葉斯網絡的過程.貝葉斯網絡學習主要包括參數學習和結構學習兩種.參數學習是指已知網絡結構，通過不斷學習來確定網絡參數.貝葉斯網絡結構學習指的是通過結合包含專家知識在內的先驗信息，以尋找出與樣本數據集擬合得最好的網絡結構[9].貝葉斯網絡結構學習主要包括基于獨立性測試的方法[10，11]和基于評分搜索[12-14]的方法.基于評分搜索的方法主要取決于評分函數的確定和搜索策略的選擇.與參數學習相比結構學習復雜的多，同時，貝葉斯網絡結構學習問題已經證明是一個NP難問題[15，16]，一般采用啟發式搜索算法解決評分搜索問題.貝葉斯結構網絡結構往往是由專家給出，相對費時、費力.當面對特別復雜的系統時，僅僅依靠專家經驗來構建貝葉斯網絡基本不可能實現，并且構造正確的貝葉斯網絡結構是進行參數學習的前提，因此本文主要研究網絡結構學習.

粒子群優化算法[17](Particle Swarm Optimization，PSO)是Kennedy和Eberhart等人通過對鳥群飛行、聚集行為的研究，在1995年提出的一種群體智能優化算法.該算法結構簡單、參數較少，在高維空間函數尋優方面具有收斂速度快，解的質量高、魯棒性好、極易實現等優點.因此在函數優化、神經網絡訓練 、模式分類、控制工程等領域受到了廣泛的關注和學習[18-21].將智能算法與貝葉斯網絡結合成為學者的研究熱點.2011年冀俊忠等人提出了一種基于蟻群優化的貝葉斯網絡結構學習算法.2013年汪春峰等人將無約束優化和遺傳算法結合，提出了一種學習貝葉斯網絡結構的限制型遺傳算法.2014年張平等人將人工蜂群算法應用到貝葉斯網絡結構學習中.2014年劉揚等人提出一種信息論結合粒子群優化的貝葉斯網絡結構學習算法.2014年高曉光等人提出一種互信息限制的貝葉斯網絡結構學習算法.

鑒于此，本文將貝葉斯網絡結構學習和粒子群優化算法看成一個組合優化的問題，將粒子群優化算法應用于貝葉斯網絡結構學習中，在此基礎上，并將遺傳算法良好的并行計算能力進行有效的結合，利用混合優化方法尋找最優的貝葉斯網絡.同時，結合貝葉斯網絡結構特點，設計基于遺傳算子的位置更新策略，主要是將粒子群算法中的位置更新策略與遺傳算子中的變異算子和交叉算子相結合，以增加學習的精度和效率.最后以標準的Alarm 和Aisa 網絡為實例，并與其他算法相比較體現了本文算法在貝葉斯網絡結構學習中較強的學習能力.

2 相關介紹

2.1 貝葉斯網絡

貝葉斯網絡(Bayesian Networks，BNs)作為概率信息的載體，主要由聯合概率分布的圖形來表示.通常一個貝葉斯網絡由有向無環圖(簡稱DAG)和條件概率表(CPT)兩部分組成：DAG中的每個節點代表一個隨機變量，根據節點之間的概率依賴關系由有向邊連接；CPT主要表示節點間的依賴強度.

對于一個n元有限隨機變量V={v1，v2，v3…vn}，Θ為變量X的聯合概率分布.則貝葉斯網絡可由二元組BN=(G，Θ)來表示，其中G=(V，E)是一個有向無環圖，V={v1，v2，v3…vn}是一個節點集合，V中的變量與圖中的節點相對應.E為有向邊集合，表示隨機變量間的依賴關系.Θ是一個網絡參數向量，為節點的條件概率分布集合，表示節點之間的依賴程度向量Θ由若干分量組成，用Θ={P[vi|π(vi)]，vi∈V}表示貝葉斯再給定節點vi的父母節點時的條件概率，且每個分量表示網絡結構中節點vi的條件概率，為一個網絡參數.π(vi)表示節點vi的父母節點，根據概率論的鏈式規則，變量vi的聯合概率分布為：

(1)

若沒有父節點，則π(vi)=Φ.

2.2 評分標準

貝葉斯網絡結構學習主要思想是首先結合先驗知識找到一個和樣本數據集合擬合最好的網絡結構，然后通過評分函數對網絡結構與樣本集合匹配程度進行測試選取最優的網絡.對于基于評分搜索的方法主要解決的是評分函數的選取和搜索方法的選擇兩個問題.其中評分函數是衡量生成的貝葉斯網絡結構與數據樣本匹配程度的度量標準.本文選用使用廣泛的貝葉斯信息標準(Bayesian information criterion，BIC)作為適應度函數，適應度值越大則表明匹配程度越高.

對于一組隨機變量V={v1，v2，v3…vn}，D={D1，D2，D3…Dn}是與這組變量對應的獨立分布的樣本集.G是v1，v2，v3…vn為節點的貝葉斯網絡.貝葉斯網絡結構學習的主要目的是首先結合先驗知識找到一個和樣本數據集合D擬合最好的網絡結構G.然后通過評分函數對網絡結構與樣本集合匹配程度進行測試選取最優的網絡.BIC評分函數是在樣本滿足獨立同分布的假設前提下，用對數似然度量結構與數據的擬合程度.BIC評分函數為[24，25]：

(2)

評分函數確定之后，接下來的任務就是選擇一個合適的搜索算法應用于貝葉斯網絡結構中以尋找出分數最高的網絡結構.對于一個包含n個節點的有向無環圖，f(n)是組成無環圖所有可能的個數，Robinson給出的f(n)的計算公式如下[26]：

f(1)=1

(3)

從上式中可以看出，隨著n的增大貝葉斯網絡的空間不斷增大.貝葉斯網絡結構學習問題已經證明是一個NP-Hard問題，故采一般地確定性的精確算法尋找最優的網絡結構，通常很難達到目標.

3 貝葉斯網絡結構的學習過程

3.1 粒子群優化算法

粒子群優化算法(Particle Swarm Optimization，PSO)起源于對鳥群捕食過程中的遷徙和聚集的模擬，是一種基于全局搜索的群體智能算法[27].在粒子群算法中，粒子可表示成一個二元組(X，V)，其中，X、V分別表示粒子的位置、速度.每個粒子都代表著搜索空間中的一個解，所有的粒子都對應一個由優化函數所決定的適應度值，解的質量通過適應度值的大小來衡量.粒子群優化算法首先初始化一群隨機粒子，然后在每一次迭代中，每個粒子根據自身的最優值(Pbest)和群體的最優值(Gbest)不斷地更新自己的位置找到最優解.

假設在d維空間中，粒子的位置定義為Xi=(Xi1，Xi2…Xid)，Vi=(Vi1，Vi2…Vid)，粒子的速度和位置更新公式如下[28]：

(4)

(5)

文獻[29]中證明了粒子的進化過程與速度無關，從而提出了一種簡化的粒子群算法(Simple Particle Swarm Optimization，SPSO)如下：

(6)

3.2 粒子編碼

由于PSO算法的搜索空間是由貝葉斯網絡結構組成.假設一個貝葉斯網絡結構有n個節點組成，則可以用一個n×n維的鄰接矩陣來表示.本文采用文獻[28，30]矩陣編碼方式進行編碼，如果把“粒子的位置”也定義為一個有向無環圖，則粒子的當前位置也可以用矩陣G表示，G={gij}，其中gij定義如下：

(7)

本文算法中的貝葉斯網絡結構的個體可表示為：g11g21…gn1g12g22…gn2…g1ng2n…gnn

圖1 粒子的位置表示Fig.1 Particle position

根據上述編碼方式，圖1的貝葉斯網絡結構編碼為0011001000010000

3.3 基于遺傳算子的位置更新策略

通過對粒子的變異操作，提高了全局搜索能力，避免陷入局部最優解；將粒子分別與局部最優值和個體最優值交叉，保護了粒子的最優狀態值也同時增強了種群的多樣性.在提高搜索效率的同時減小了粒子更新的隨機性，通過對貝葉斯網絡邏輯關系的推理，獲得準確的貝葉斯網絡結構.采用標準的Alarm[22]和Asia[23]網絡驗證了本文算法用于貝葉斯網絡結構學習的可行性及優越性，與其他算法相比取得了較好的實驗效果.

3.3.1 變異算子操作

(8)

式中，r1是[0，1)的隨機數，A代表變異操作，ω代表變異概率.

本文算法對貝葉斯網絡結構設計了3種變異操作如圖2所示.

圖2 變異操作后的撲拓圖Fig.2 A graph of variation operations

圖2中(a)表示當前模型，圖2(b)表示從1→4減邊操作，圖2(c)表示從2→3反轉邊操作，圖2(d)表示從2→4加邊操作，圖2(e)表示從2→4多邊操作，因為導致回路，記為不合法的多邊操作.通過加邊、刪邊或反轉邊三種操作之后，再計算網絡得分.由圖2可以看出，粒子在交叉變異之后可能產生非法的拓撲圖，如圖2(e)所示.對此本文采用文獻[34]所提出的修復算子對非法的貝葉斯網絡結構進行去環操作.具體修復步驟如下：

1)求出與網絡拓撲圖相對應矩陣的傳遞閉包；

2)通過判斷閉包矩陣對角線上的元素是否全為0，判斷是否為合法的網絡結構.如果全為0，代表合法的網絡結構；否則，保留主對角線上非0元素對應的節點(這些節點均位于環內).

3)任取環內的一點，求出它位于閉環內的所有父節點；

4)對于這些父節點指向該節點的任一邊，進行刪除或者反轉操作，使得網絡結構中不存在有向環.

3.3.2 交叉算子操作

交叉是產生新個體的核心，因為交叉后的子代是父代的繼承和重組，所以種群適應度較高，交叉的結果直接影響后代的質量和收斂的結果.為了豐富種群的多樣性本文采用兩點交叉的方式.在粒子與個體最優粒子交叉的過程中，將會產生新粒子，主要是由父母粒子的公共部分和隨機部分組成，公式如下：

(9)

(10)

式中，Cp代表新粒子與個體最優粒子的交叉操作，Cg表示個體最優粒子與全局最優粒子的交叉過程.c1，c2是學習因子，c1代表個體最優粒子的交叉概率，c2全局最優粒子的交叉概率，r2，r3是[0，1)的隨機數.

3.3.3 粒子位置更新

粒子群算法通過遺傳算子的操作，粒子的位置更新公式為：

(11)

慣性權重可以衡量算法的全局搜索能力和局部搜索能力，當ω過大會時有利于全局搜素，當ω過小時有利于局部搜索.學習因子是調節粒子自身經驗和群體經驗的權重，主要反應粒子之間的交流，影響粒子的運動軌跡.c1，c2過大過小都會對結果造成不利影響.通過以上分析，通過公式(12～14)線性的調節ω，c1，c2.

(12)

(13)

(14)

式中，t表示當前迭代次數，Itmax代表最大的迭代次數.ωs，c1s，c2s代表初始值，ωe，c1e，c2e代表最終值.

3.4 算法實現步驟

本文提出一種混合粒子群優化的貝葉斯網絡結構學習算法，簡稱MMPSO算法.在該算法中，粒子所處的位置代表一種網絡結構圖，用鄰接圖來表示每一個候選網絡結構，如何產生初始解以及對每一個解的鄰域進行搜索是解決問題的關鍵.

由于初始種群的選擇對PSO算法的尋優性能影響較大，對此，本文通過采用文獻[35]提出的最大權生成樹(MWST)算法，產生網絡結構的初始邊集，然后生成與貝葉斯網絡結構擬合度最好的樹結構.以此作為初始結構圖模型，通過對其任意添加邊、刪除邊或反轉邊的操作以產生該圖的所有鄰接圖.在相應鄰接矩陣中選用一定數量的矩陣作為初始種群，然后在利用MMPSO搜索出最優的網絡結構.具體的步驟如下：

實驗步驟：

Step1. 設置初始參數：種群規模N，交叉概率，最大迭代次數等.

Step2. 產生初始種群：通過MWST算法得到初始結構圖，并由此生成該圖的所有鄰接圖，選取N個結構圖作為初始粒子群；

圖3 本文算法的流程圖Fig.3 Flow chart of algorithm

Step3. 根據公式(2)計算每個粒子的BIC評分值作為適應度值；并且根據個體極值找出全局最大的評分值和對應的初始全局粒子，同時更新每個粒子的最優解；

Step4. 根據公式(8)，調節權重系數對粒子進行變異操作，產生新粒子，以保證良好的全局搜索能力和局部搜索能力.

Step5. 根據公式(9)，將新粒子與個體最優粒子進行交叉操作，同時更新個體極值與全局極值；

Step6. 根據公式(10)，將粒子的個體極值與全局極值進行交叉操作，更新新的個體極值與全局極值.

Step7. 通過公式(2)重新計算每個粒子的BIC值，并更新種群的全局最優解.

Step8. 判斷粒子是否早熟收斂，如果早熟收斂則執行混沌優化搜索策略，否則繼續執行粒子群算法.

Step9. 檢查終止條件(達到足夠好的位置或達到最大迭代次數)；如果達到條件，運行終止并輸出全局最優解.否則，返回Step 4.

4 實驗結果與分析

4.1 數據集

本文仿真實驗的硬件環境是內存4G，CPU為Intel(R)CoreTM，2 Duo 3.30GHz.實驗平臺Windows7，程序實現采用Matlab2014b，以及軟件工具包FullBNT-1.0.4.

表1 標準測試網絡數據組成Table 1 Standard test network data

為了測試MMPSO的性能，本文采用通用的Benchmark數據集Alarm網絡、Asia網絡、網絡來完成實驗.具體數據組成如表1所示.

4.2 參數設置及評價指標

根據標準的Alarm、Asia網絡，首先利用BNT工具依照標準概率表生成樣本量為500、1000、2000、3000、5000的數據集，然后使用本文算法分別對產生的多組數據進行結構學習，并將運行的結果與MMHC算法[36]、GA算法[31]、PSO算法[28]、MMACO算法[37]進行對比，以此來驗證本文算法的效果.在仿真實驗中，設置初始粒子種群為20，算法的參數設置如下：(注：初始化種群的規模為N，螞蟻數量為m，信息揮發因子為ρ，權重參數α，β，交叉概率記為pc，變異概率記為pm)最大迭代次數Tmax=100.

表2 參數設置Table 2 Parameter setting

實驗采用的評價指標是BIC評分值和結構漢明距離(SHD)，評價各種算法所學到的最優結構與標準網絡結構之

間的差異.其中

1)標準BIC評分準則，表示算法學習到最優結構時的BIC得分值，分值越高結果越好.

圖4 標準的Aisa網Fig.4 Aisa network圖5 學習后的Aisa網Fig.5 Learned Aisa network

2)結構漢明距離(SHD)：用來度量學習所得網絡結構與真實網絡結構之間的差異程度，漢明距離用網絡結構中多余邊、反轉邊和丟失邊三者之和來表示.取值越小代表結果越好.

4.3 實驗結果分析

為了驗證本文算法用于貝葉斯網絡結構的效果，主要從三個方面進行實驗：算法的可行性、算法的準確性、算法的收斂性.

圖6 標準的Alarm網Fig.6 Alarm network圖7 學習之后的Alarm網Fig.7 Learned Alarm network

4.3.1 可行性分析

本文以經典的Alarm網和Asia網作為仿真模型來評價算法的性能.圖4，圖5分別標準和學習之后的Aisa網絡，圖6，圖7分別標準和學習之后的Alarm網絡.從上圖中可以看出，MMPSO在Aisa網的結構中，學習的最優效果是僅出現從long到smoking的反向邊，準確率高于其他算法.在復雜度更

表3 不同算法對Alarm網絡學習的結果Table 3 Results of different algorithms for learning Alarm networks

高的Alarm網的結構學習中，與標準的網絡對比，分別多3條邊，少一條邊，反向一條邊.(多11至15，27至37，33至15；少：17-26，27-36；反：30至29)最有效果是僅有一條反向邊，較其他算法精確率較高.

表4 不同算法對Aisa網絡學習的結果Table 4 Results of different algorithms for learning Aisa networks

4.3.2 準確性分析

為了驗證MMPSO算法學習網絡結構的準確性，將本文算法與與MMHC算法[36]、GA算法[31]、PSO算法[28]、MMACO算法[37]分別對產生的多組數據進行結構學習，并將運行結果進行比較.表3，表4給出了各種算法獨立運行10次得到的BIC評分的均值和SHD值，其中表示獨立運行10次BIC評分的均值和方差.

從表3至表4的數據中可以看出，在樣本容量相同的情況下，MMPSO算法學習的網絡結構更接近真實網絡，說明了本文算法的有效性.在Alarm數據集上，MMHC學習效果最差，本文算法所學習的效果最好，MMACO次之，雖然隨著樣本數量的逐漸增大，MMACO與MMPSO的得分值逐漸接近，但是本文算法在保證得分較高的情況下，SHD取值最小.在Alarm數據集上，與其他算法相比，MMPSO的SHD取值與MMABC取值相同，但MMPSO的得分值高于MMACO，得到了比較好的結果.

4.3.3 收斂性分析

為了驗證MMPSO算法的收斂性能，以Aisa為例，用2000組數據做實驗并將本文算法與PSO和MMACO算法的結果進行對比.用2000組數據做實驗并將MMPSO與PSO和MMACO算法的結果進行對比.

表5 3種算法學習Aisa網絡的結果比較Table 5 Results of learning Aisa networks

設置每種算法的最大迭代次數為50次，初始粒子為20，表5列出了三種算法分別運行10次，學習Aisa網絡的BIC得分和相應有效時間的均值.由表8和收斂結果圖可以看出，

圖8 學習Aisa網絡的結果收斂圖Fig.8 Convergence of Aisa network

本文算法雖然在時間上略高于PSO算法，但是MMPSO算法學習到的網絡結構得分最高，收斂效果最好，得到了最好的結構.由此證明了MMPSO算法的有效性以及可行性.

5 總結與展望

本文將粒子群優化法和遺傳算法的優點進行有效的結合，并將混沌優化策略引入混合算法中形成了一種混合優化算法，同時將這種算法應用到貝葉斯網絡結構學習方法中.通過對粒子的變異操作，提高了全局搜索能力，避免陷入局部最優解；將粒子分別與局部最優值和個體最優值交叉，保護了粒子的最優狀態值也同時增強了種群的多樣性.在提高搜索效率的同時減小了粒子更新的隨機性.通過實驗結果可以看出，本文在樣本容量相同的情況下，學習的網絡結構更接近真實網絡，并且在最短的有效時間內，學習到分值較高的網絡結構，收斂效果最好.但是，在實驗中發現，某些明顯錯誤邊，特別是在多余邊問題上的改善不是很明顯，這是算法的不足之處.算法仍需改進優化，將是下一步研究的重點.