劉彬，范瑞星，劉浩然，張力悅，王海羽，張春蘭

（1. 燕山大學信息科學與工程學院，河北 秦皇島 066004；2. 河北省特種光纖與光纖傳感重點實驗室，河北 秦皇島 066004）

1 引言

貝葉斯網絡（BN, Bayesian network）是結合圖論和概率論來表示因果知識的概率圖模型，是用于不確定領域中推理和預測的最佳方式之一[1]。貝葉斯網絡可以用圖論的語言直觀地揭示問題的結構，并利用該結構降低概率推理的計算復雜度。由于貝葉斯網絡直觀易懂，在風險分析、機器學習、信息學等研究領域[2-3]都有應用。

貝葉斯網絡的構建包含結構學習、參數學習和推理學習。結構學習是基礎與核心，完備數據下的結構學習方法主要有3種：基于依賴性測試的方法[4]、基于評分搜索的方法[5]和混合方法[6]，其中常見的結構學習方法是基于評分搜索的方法，即在所有節(jié)點的結構空間內按照一定的搜索策略及評分準則構建貝葉斯網絡結構。

基于評分搜索的方法學習貝葉斯網絡結構是一種 NP問題[7]，國內外學者通常利用啟發(fā)式算法來解決此類問題。Tsamardinos等[8]提出了一種基于依賴性測試和爬山算法的最大最小爬山（MMHC,max-min hill-climbing）算法，該算法雖然改善了檢索策略，降低了搜索空間復雜度，但由于搜索空間的縮小易導致算法陷入局部最優(yōu)。劉浩然等[9]提出了基于最大支撐樹（MWST, most weight supported tree）和蟻群算法（ACO, ant colony optimization）的混合搜索節(jié)點序算法（MAK, MWST-ACO-K2），該算法在處理小型網絡時可取得較理想的結果，但是與其他基于節(jié)點序搜索算法類似，需要對種群中所有個體運行K2算法得到對應的適應度值，在大網絡中存在時間復雜度較高、結果較差等問題。Wang等[1]提出了基于離散水循環(huán)算法的貝葉斯結構學習算法（BEWCA-BN, binary encoding water cycle algorithm for BN structures learning），該算法根據邏輯算子提出了改進的水循環(huán)算法更新個體，其中蒸發(fā)策略雖然提高了算法跳出局部最優(yōu)的能力，但在小網絡中通常需要花費更多的時間尋找最優(yōu)解。Contaldi等[10]提出了基于精英遺傳算法的貝葉斯結構學習算法（AESL-GA, adaptive elite-based structure learner using genetic algorithm），該算法采用自適應的控制參數來避免參數設置對結果的影響，在小網絡中學習到了較優(yōu)的網絡結構，但在大網絡中由于縮小搜索空間導致學習到的結果不太理想。

上述啟發(fā)式算法應用于貝葉斯網絡結構學習時由于參數設置對搜索過程的影響，存在搜索效率較低、易陷入局部最優(yōu)等問題，而Seyedall等[11]提出的樽海鞘算法（SSA, salp swarm algorithm）具有參數較少、操作簡單、易于實現等優(yōu)點。SSA將種群劃分為引領者（leader）與跟隨者（follower），通過引領者領導跟隨者形成 slaps鏈進行種群尋優(yōu)。在處理復雜問題時由于引領者對跟隨者的引領作用，容易使整個種群過早地聚集于局部最優(yōu)解[12]。

本文將SSA應用到貝葉斯網絡結構學習，同時為了避免SSA陷入局部最優(yōu)，將全局尋優(yōu)能力強的差分進化（DE, differential evolution）算法[13]與SSA融合，構建了一種混合樽海鞘-差分進化（HBSS-DE,hybrid binary salp swarm-differential evolution）算法。該算法首先設置規(guī)模因子將種群劃分[14]為較好種群和較差種群。然后構建樽海鞘搜索策略更新較好種群，提高算法的收斂速度；構建差分搜索策略更新較差種群，避免陷入局部最優(yōu)，并在合并子種群時利用變異算子增大搜索范圍。最后通過種群的迭代搜索最佳結構。

2 HBSS-DE算法的構建

HBSS-DE算法利用最大支撐樹與爬山算法建立初始種群P0，然后將P0降序排列更新為種群P，并設置規(guī)模因子將種群P劃分為較好種群P1與較差種群P2。構建樽海鞘搜索策略更新P1，建立差分搜索策略更新P2。合并為種群P時利用兩點變異算子增加種群的多樣性，并根據規(guī)模因子重新劃分P進入下次迭代。迭代結束輸出種群P中評分最高的樽海鞘個體，即最佳的貝葉斯結構。

2.1 初始種群的建立及劃分

根據數據樣本計算目標網絡各節(jié)點間的互信息，利用2個節(jié)點間的互信息可以得出這2個節(jié)點是否相關[15]，不存在相關關系的節(jié)點必然不存在因果關系。以變量X、Y為例，互信息I(X,Y)為

其中，p(X,Y)為變量X和變量Y的聯合概率，p(X)為變量X的概率，p(Y)為變量Y的概率。

根據互信息計算各節(jié)點之間的權重，利用最大支撐樹原則生成一個候補結構[16]，可以有效地縮小搜索的空間。該候補結構中，除樹形結構之外的節(jié)點利用爬山算法中的加邊、減邊、轉邊算子得到樽海鞘種群。其中每個樽海鞘代表的貝葉斯網絡結構可用如圖1所示的鄰接矩陣表示。

圖1 樽海鞘個體

采用貝葉斯信息準則（BIC, Bayesian information criterion）評分函數對樽海鞘個體評價。根據BIC評分函數的可分解性，當BN中的局部結構改變時，為了減少重復計算的次數，只需利用式(2)計算改變的局部結構G1的評分fnew(G1,D)，然后代入式(3)即可得到樽海鞘個體的評分f(G,D)。

其中，Xi表示節(jié)點，π(Xi)表示Xi的父節(jié)點；mijk表示數據中滿足π(Xi)組合，即取值為j且Xi=k的樣本數；qi表示π(Xi)取值共有qi種組合；ri表示Xi共有ri種取值；D表示數據的樣本；G表示鄰接矩陣，G2表示G中未改變的局部結構，fold(G2,D) 表示G2的評分。

選取樽海鞘種群中高于平均評分的個體作為初始種群P0。文獻[14]在種群劃分階段將初始種群隨機分成2個相等的子種群更新個體，由于個體分布隨著種群迭代動態(tài)變化，導致迭代后期無法提高算法的局部搜索能力。為了保證種群全局搜索與局部搜索的平衡，本文根據種群的進化情況建立自適應劃分種群的規(guī)模因子q，具體的形式為

其中，NP為P0的種群規(guī)模，Tm為最大迭代次數，t為迭代次數，■*■為向下取整函數，K為自適應調整規(guī)模因子q的參數，是根據個體評分與該種群平均評分的相對值建立的，具體的形式為

其中，fmax為P的最高評分值，f為P的平均評分值，fmin為P的最低評分值。

由于在迭代初期，種群的個體分布較分散，即K的值較小，因此q的值較小保證算法的全局搜索。隨著迭代的進行，種群的個體分布趨于集中，即K的值也逐漸增大，因而q的值較大保證算法的局部搜索。

根據樽海鞘個體的評分將P0降序排列更新為種群P，利用規(guī)模因子q將種群劃分為較好種群P1與較差種群P2，具體形式為

其中，P1的種群規(guī)模NP1=q，P2的種群規(guī)模NP2= N P- N P1。

2.2 子種群更新

為了提高HBSS-DE算法的收斂速度，利用自適應的變異算子與交叉算子改進 SSA的引領者與跟隨者，建立樽海鞘搜索策略更新較好種群P1。同時為了提高算法的收斂精度，改進DE算法中的變異算子與交叉算子，建立差分搜索策略更新較差種群P2。

其中，rand為[0,1]之間的隨機數，“1”為加邊操作，“0”為減邊操作為變異位置的值。

其中，Crmax為交叉概率上限；Crmin為交叉概率下限；M為控制交叉方式的參數，它是根據個體的差異，即當代種群中個體最佳評分與平均評分的相對值建立的，具體的形式如式(13)所示。

其中，fmax為P中最高評分值，為P的平均評分值，fmin為P中最小評分值。由式(13)可知，在算法的迭代初期，M的值較小，因而Cr的值較小，算法能以較大的概率選擇單點交叉法更新個體。隨著迭代的進行，M的值逐漸增大，Cr的值也逐漸增大，樽海鞘個體以較大的概率選擇兩點交叉法更新個體。

將自適應因子代入式(14)的交叉算子中更新P1中的跟隨者個體即

其中，F1是單點交叉，即隨機選擇一列的位置，然后交換 2個鄰接矩陣對應列的位置；F2是兩點交叉，即交換2個鄰接矩陣對應兩列的位置。以圖2所示的單點交叉為例，交換矩陣A和矩陣C的最后一列，得到矩陣E。

圖2 單點交叉

1) 差分搜索策略

改進的變異算子先從P2中隨機選取3個樽海鞘個體（其中），然后將代入式(15)產生差分矩陣根據式(8)得到的變異位置集合R。利用式(16)將中的變異位置值轉化為變異概率并代入式(17)得到的變異位置值

其中，k2? [ 1,NP2]；r1、r2、r3為 [1 ,NP2]之間的隨機整數，并且r1≠r2≠r3；rand為 0～1的隨機數；c1(i,j)為變異概率；b=20為帶寬因子[17]；F=0.5為

改進的交叉算子是通過式(12)得到P2的自適應因子Cr，代入式(18)將變異個體與原個體交叉得到更新個體即

2) 選擇操作

2.3 合并子種群

為了能夠用更少的時間搜索到最佳網絡結構，文獻[14]將更新之后的子種群P1、P2合并為種群P，使種群中的個體在搜索空間中彼此共享位置信息。但是種群個體的分布情況是動態(tài)變化的，即在算法的迭代初期，個體分布比較分散，種群P1、P2中相同的個體較少；隨著迭代的進行，個體趨于集中，種群P1、P2中相同的個體也逐漸增多，種群P的多樣性易喪失。

由于迭代后期種群的多樣性遭到破壞，算法易陷入局部最優(yōu)。為保證算法迭代后期有能力跳出局部最優(yōu)，本文將種群P1、P2中相同的個體通過兩點變異增加種群的多樣性。根據式(20)將P2中與P1評分相同的個體保存到種群P3，評分不同的個體保存到種群P4。

其中，k1? [ 1,NP1]，k2? [ 1,NP2]，f(xk

t1+1)、f(xk

t+21)分別是種群P1、P2中每個個體的評分。

將P3中的個體通過兩點變異增加種群多樣性，其中兩點變異是對鄰接矩陣中任意一列的2個元素進行取反，如圖3所示對矩陣Q的右上角2個元素進行取反，得到矩陣W。

圖3 兩點變異

將更新種群3P與1P、4P合并成新的種群0P，根據個體評分降序排列更新種群P，利用規(guī)模因子重新劃分種群進入下一次迭代。迭代結束后保留種群P中的評分最佳的樽海鞘個體，即最優(yōu)的貝葉斯結構。

2.4 HBSS-DE算法的整體流程

本文2.1節(jié)～2.3節(jié)介紹了算法的原理，接下來，對算法的具體步驟進行簡要描述。

步驟1 輸入數據樣本，通過互信息建立最大支撐樹，初始化參數Tm，t=1。

步驟2 利用爬山算法對最大支撐樹進行初步搜索得到樽海鞘種群，并選擇高于平均值的樽海鞘個體建立初始種群0P。

步驟3 將0P通過評分降序排序更新為種群P，并通過自適應因子劃分種群P為較好種群1P與較差種群2P。

步驟4 種群1P中采用樽海鞘搜索策略進行更新，提高算法的收斂速度。

步驟5 種群P2中采用差分搜索策略進行更新，提高算法的收斂精度。

步驟6 將2個子種群中相同個體利用兩點變異增加種群的多樣性，變異操作之后合并為種群P0。

步驟7 若滿足t＜Tm，則t=t+1，并轉至步驟3；否則，輸出最高評分值對應的貝葉斯結構。

2.5 HBSS-DE算法的收斂性分析

根據 Solis等[18]提出的概率測度法分析隨機搜索算法HBSS-DE的收斂性。

引理1 隨機搜索算法 HBSS-DE滿足f(N(z,ξ) ) ≥f(z)，若ξ∈Sgbest，則f(N(z,ξ) ) ≥f(ξ)。其中，f為HBSS-DE算法解決最大化問題的適應度函數；N為產生較優(yōu)最新解的算子；z為搜索空間S的最優(yōu)解空間Sgbest中的樽海鞘個體，也是在Sgbest上產生可接受函數值的上確界；ξ為算法在S上隨機生成的樽海鞘個體。

證明根據 2.2節(jié)中選擇操作的描述，可將HBSS-DE算法中對當前最優(yōu)解的選擇算子N定義為

引理2 HBSS-DE算法的最優(yōu)解空間Sgbest的概率測度大于0，即

證明假設HBSS-DE算法的貝葉斯結構群為X={x1,x2,…,xn} ,x∈S，其中搜索空間S是可列集或有限集，顯然它的 Lebesgue[19]總是大于 0，即L[S] ＞0 ，HBSS-DE算法的最優(yōu)解空間Sgbest屬于S的一個Borel[19]子集，由HBSS-DE算法的最優(yōu)解空間Sgbest的定義可知，證畢。

引理3 HBSS-DE算法中，當滿足L[Sgbest]＞0時，式(22)成立。

其中，μt(?)為第t次迭代結果的概率測度。

證明當滿足L[Sgbest] ＞0時，有0＜μi,t(Sgbest)＜1，可知由μt產生的對Sgbest的概率測度為

將μt(Sgbest)代入式(22)，可得

至此，引理3得證。根據Solis等[18]提出的概率測度法可知，當HBSS-DE算法滿足引理1～引理3時，式(25)成立。其中表示在第t次的結果xt落到Sgbest里的概率值為1，即算法經過有限次迭代后，HBSS-DE算法中一定有存在于Sgbest的個體。尋找到最優(yōu)解空間Sgbest之后，根據HBSS-DE算法的迭代性原則，在以后的所有迭代中，種群中所有個體都向Sgbest靠攏，最終收斂于Sgbest。

3 實驗仿真

為了驗證HBSS-DE算法的性能，在操作系統為Windows7，處理器為Intel i3 3.40 GHz CPU，內存為4 GB，軟件環(huán)境為Matlab 2010a，基于貝葉斯網絡工具箱 FullBNT-1.0.7中的 ASIA網絡、CAR網絡、ALARM網絡進行仿真實驗。其中，標準的ASIA網絡具有8個節(jié)點、8條邊，標準的CAR網絡具有12個節(jié)點、9條邊，標準的ALARM網絡具有37個節(jié)點、46條邊。3種標準網絡是通過訓練實際數據建立的貝葉斯網絡，其中節(jié)點之間的邊表示因果關系。在3種網絡中隨機生成數據量為500、1 000、3 000、5 000個的數據樣本，因為數據是隨機產生的，為降低數據的隨機性對實驗結果準確性的影響，每種樣本容量分別產生 10個，且每個數據單獨運行10次，即每組數據運行100次后取平均值作為實驗結果。

隨著節(jié)點的增加，采用評分搜索的方法學習全局最優(yōu)解是一種 NP問題[7]，而全局最優(yōu)解與較其稍差的局部最優(yōu)解的實際效益相差不大[18]，因而采用足夠好的解Sgbest來代替全局最優(yōu)解，將NP問題轉化為N-問題。HBSS-DE算法在3種網絡中的實驗結果分別如表1～表3所示，其中fbest是100次實驗中的最佳BIC評分，fav是100次實驗的平均BIC評分，fsn是標準網絡的BIC評分。利用驗證算法的局部最優(yōu)解Sgbest是否接近全局最優(yōu)解。

由表1～表3可知，在小網絡（如ASIA網絡、CAR網絡）中，當數據量＞500個時HBSS-DE算法搜索到的最佳結構評分fbest與標準網絡的評分fsn幾乎相同，即在小網絡中可以搜索標準網絡，這是因為小網絡中節(jié)點較少、搜索空間較小。而在大網絡中（如ALARM網絡），HBSS-DE算法的fbest比較接近fsn，這是因為隨著節(jié)點的增多，搜索空間變大，算法收斂于接近全局最優(yōu)解的局部最優(yōu)解。在ASIA網絡中σ的平均值為0.004 91，在CAR網絡中σ的平均值為0.010 4，在ALARM網絡中σ的平均值為0.023 9，總體而言HBSS-DE算法的局部最優(yōu)解Sgbest接近全局最優(yōu)解。

表1 ASIA網絡HBSS-DE算法的BIC評分

表2 CAR網絡HBSS-DE算法的BIC評分

表3 ALARM網絡HBSS-DE算法的BIC評分

根據文獻[21]，列出分析算法性能的評價指標。

Ext（execution time）：找到最佳網絡所需要的時間（單位為s）。

TP（true positive）：真正例，即預測結構與標準結構中相同且方向為正的邊。

FP（false positive）：假正例，即預測結構中邊的方向為正，標準結構中為負。

FN（false negative）：假負例，即預測結構中邊的方向為負，標準結構中為正。

TPR（true positive rate）：正確率，即預測結構中TP與標準結構中邊的比值，

FPR（false positive rate）：精確率，即預測結構中TP與預測結構全部邊的比值，

F（F-measure）：綜合評價指標，即正確率和召回率的調和平均值，

W：100次實驗中HBSS-DE算法的實驗值ψ與其他算法的實驗值ε相比提升的百分比，即

將 HBSS-DE算法與 AESL-GA、MAK、BEWCA-BN、MMHC算法進行對比實驗。其中本文算法自適應因子的上下限分別為0.85、0.2。根據文獻[1]設BEWCA-BN尺度因子c的上下限分別為0.9、0.7，搜索強度參數d的上下限分別為0.2、0.01。根據文獻[7]設 MMHC獨立性測試的置信度α為0.01，統計閾值設置為0.05。根據文獻[8]設MAK信息素強度系數為 1，信息素蒸發(fā)系數為 0.7，啟發(fā)式因子權重為2。根據文獻[9]設AESL-GA獨立性測試的置信度為0.01，精英個體閾值α為0.09。實驗結果如下：表4～表6是不同算法在不同網絡的平均BIC評分，圖4～圖6是不同算法在不同網絡的結構對比，表7～表9是不同算法在不同網絡的執(zhí)行時間對比。

由表4～表6可知，當數據量為500時，在ASIA網絡與CAR網絡中，本文算法和BEWCA-BN算法均能找到相對較優(yōu)的評分，但隨著數據量的增大，本文算法相比于其他算法能學習到更優(yōu)的評分。根據W計算，當數據量為500、1 000、3 000、5 000時，本文算法與其他算法相比，BIC評分的平均提高值分別為然后利用得到整體 BIC提升百分比。本文算法在 ASIA網絡中的 BIC評分比MAK、MMHC、AESL-G、BEWCA-BN 算法平均提升了 3.71%、4.34%、3.85%、1.86%，在 CAR網絡中平均提升了4.41%、3.64%、2.49%、1.47%，在 ALARM 網絡中平均提升了 6.96%、9.08%、10.02%、3.84%?？傮w而言，HBSS-DE算法能夠找到評分更優(yōu)的貝葉斯網絡結構，在大型網絡中比其他算法的提升效果更佳。這是由于本文算法可以在較好個體周圍局部搜索，在較差個體周圍全局搜索，并通過規(guī)模因子保證了算法的局部搜索與全局搜索的平衡，即在迭代前期全局搜索，迭代后期局部搜索。同時對于2個子種群中相同個體利用兩點變異策略增加了種群的多樣性，避免陷入局部最優(yōu)。

表4 不同算法在ASIA網絡中不同數據量下的平均BIC評分

表5 不同算法在CAR網絡中不同數據量下的平均BIC評分

表6 不同算法在ALARM網絡中不同數據量下的平均BIC評分

圖4 ASIA網絡不同數據下的結構對比

圖5 CAR網絡不同數據下的結構對比

圖6 ALARM網絡不同數據下的結構對比

表7 不同算法在ASIA網絡中不同數據量下的執(zhí)行時間

表8 不同算法在CAR網絡中不同數據量下的執(zhí)行時間

表9 不同算法在ALARM網絡中不同數據量下的執(zhí)行時間

由圖4～圖6可知，本文算法的TP值在ASIA網絡中比MAK、MMHC、AESL-GA、BEWCA-BN算法平均多了1.49、1.34、1.02、0.79，在CAR網絡中比其他算法平均多了2.21、1.71、1.23、0.56，在ALARM網絡中比其他算法平均多了7.41、4.11、4.81、3.65。本文算法在ASIA網絡、CAR網絡、ALARM 網絡中 FP的值均小于其他算法，且在 3個網絡中的TPR值均大于其他算法。根據W計算，當數據量為500、1 000、3 000、5 000時，本文算法與其他算法相比，F值平均提高的百分比分別為然后利用計算整體提高百分比。本文算法的F值在ASIA網絡中比MAK、MMHC、AESL-GA、BEWCA-BN算法平均提高了22.54%、17.89%、15.17%、8.57%，在CAR網絡中比其他算法平均提高了25.74%、19.29%、13.32%、9.64%，在 ALARM 網絡中比其他算法平均提高了24.71%、15.19%、10.91%、6.36%?？傮w而言，HBSS-DE算法能夠找到更準確的貝葉斯網絡結構，這是由于在較差種群中采用差分搜索策略，即在較差個體周圍進行全局搜索增加了種群的多樣性；通過自適應因子保證種群在迭代前期的全局搜索能力強；同時合并子種群時的兩點變異增大種群的搜索范圍，增加了種群跳出局部最優(yōu)的能力。

由表7～表9可知，當ASIA網絡的數據量為500和1 000時，AESL-GA運行時間最短，HBSS-DE算法運行時間比AESL-GA略長一些，這是由于在小數據量下HBSS-DE算法為了提高算法的全局搜索，通過規(guī)模因子將種群劃分為2個子種群更新個體，并將子種群中相同的個體利用兩點交叉變異擴大搜索的范圍。雖然避免了陷入局部最優(yōu)，但同時增加了運行時間。根據W計算，當數據量為500、1 000、3 000、5 000時，HBSS-DE算法與其他算法相比，執(zhí)行時間平均縮短的百分比分別為其整體縮短的百分比根據得到。與MAK算法、MMHC算法、BEWCA-BN算法相比，HBSS-DE算法在 ASIA網絡中的運行時間平均縮短了15.25%、21.77%、5.32%。在CAR網絡中，HBSS-DE算法比 MAK、MMHC、AESL-GA、BEWCA-BN算法平均縮短了18.57%、13.26%、15.18%、5.12%。當ALARM網絡數據量為5 000時，HBSS-DE算法比BEWCA-BN算法耗時略長，這是由于節(jié)點較多的網絡搜索的范圍較大，在迭代的過程中差分搜索策略雖然增大了搜索的空間，提高了算法的全局搜索能力；合并子種群時的兩點交叉變異雖然增加了種群的多樣性，避免了算法陷入局部最優(yōu)，但是差分搜索策略與兩點交叉變異也增加了算法的搜索時間。HBSS-DE算法在ALARM網絡中比MAK、MMHC、AESL-GA、BEWCA-BN算法平均縮短了28.23%、92.63%、19.04%、1.42%?？傮w而言，HBSS-DE算法在大型網絡中比其他算法的搜索效率的優(yōu)勢更明顯。這是因為樽海鞘搜索策略的局部搜索能力強，能夠快速地收斂到最優(yōu)解。而自適應因子保證了迭代后期的局部搜索，提高了搜索效率；同時在評價個體時僅計算個體改變的部分結構，減少在評價個體時消耗的無用時間。

4 結束語

本文提出了基于混合樽海鞘-差分算法的貝葉斯網絡結構學習算法HBSS-DE。該算法通過設置規(guī)模因子，將種群劃分為2個子種群，利用構建的樽海鞘搜索策略與差分搜索策略分別更新不同的子種群，在合并子種群時利用兩點交叉變異增加種群的多樣性。仿真實驗證明，HBSS-DE算法能夠找到最佳的貝葉斯網絡結構，其中規(guī)模因子平衡了局部搜索與全局搜索，樽海鞘搜索策略提高了算法的尋優(yōu)效率，差分搜索策略提高了算法的收斂精度。與其他算法相比，HBSS-DE算法的收斂精度與尋優(yōu)效率均有提升。