金 焱 胡云安 張瑾 黃雋

(1海軍航空工程學(xué)院控制工程系，煙臺(tái) 264001)(2中國(guó)人民解放軍91213部隊(duì)裝備部，煙臺(tái) 264007)(3海軍航空工程學(xué)院指揮系，煙臺(tái) 264001)

貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)是目前貝葉斯網(wǎng)絡(luò)研究中的熱點(diǎn)和難點(diǎn)問(wèn)題，其目的是要通過(guò)對(duì)數(shù)據(jù)集的分析，發(fā)現(xiàn)節(jié)點(diǎn)之間的依賴關(guān)系，從而獲得與數(shù)據(jù)集吻合最好的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu).貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)是一個(gè)NP難題，在實(shí)際中，人們主要通過(guò)啟發(fā)式方法來(lái)獲取近似最優(yōu)的網(wǎng)絡(luò)結(jié)構(gòu).貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)的學(xué)習(xí)方法一般分為基于獨(dú)立性測(cè)試的方法和基于評(píng)分搜索的方法.前者通過(guò)計(jì)算節(jié)點(diǎn)之間的互信息或條件互信息來(lái)挖掘節(jié)點(diǎn)之間的因果依賴關(guān)系，從而構(gòu)建出網(wǎng)絡(luò)結(jié)構(gòu)，其優(yōu)點(diǎn)是比較直觀，貼近貝葉斯網(wǎng)絡(luò)的語(yǔ)義，把條件獨(dú)立性測(cè)試和網(wǎng)絡(luò)結(jié)構(gòu)的搜索分離開，不足之處是對(duì)條件獨(dú)立性測(cè)試產(chǎn)生的誤差非常敏感［1］；后者利用一定的評(píng)分準(zhǔn)則，通過(guò)在候選網(wǎng)絡(luò)空間中進(jìn)行評(píng)分搜索，得到評(píng)分最高的網(wǎng)絡(luò)結(jié)構(gòu)作為最優(yōu)網(wǎng)絡(luò)結(jié)構(gòu)，其優(yōu)點(diǎn)是過(guò)程簡(jiǎn)單規(guī)范，在稠密的網(wǎng)絡(luò)結(jié)構(gòu)下學(xué)習(xí)效率較高，但該類算法往往由于采用爬山法的搜索策略而不能找到最優(yōu)的網(wǎng)絡(luò)結(jié)構(gòu).基于獨(dú)立性測(cè)試的方法的代表是Cheng［2］等提出的三段式學(xué)習(xí)算法；基于評(píng)分搜索的方法的代表是Cooper和Herskovits提出的K2算法［3］.后來(lái)研究學(xué)者又在此基礎(chǔ)上通過(guò)融合先驗(yàn)知識(shí)、改進(jìn)評(píng)分函數(shù)、改進(jìn)搜索策略、結(jié)合其他組合優(yōu)化算法，以及混合兩類結(jié)構(gòu)學(xué)習(xí)方法等途徑提出了很多改進(jìn)算法.

眾所周知，K2算法有2個(gè)前提條件，就是要事先給出節(jié)點(diǎn)的順序和每個(gè)節(jié)點(diǎn)的父節(jié)點(diǎn)個(gè)數(shù)的上限，在此基礎(chǔ)上K2算法具有非常優(yōu)秀的學(xué)習(xí)性能.但是，這2個(gè)前提條件在很多情況下是無(wú)法滿足的，這就給K2算法的應(yīng)用造成了困難.模擬退火(simulated annealing，SA)［4］是一種基于蒙特卡洛迭代求解的全局組合尋優(yōu)算法，其全局尋優(yōu)能力是以嚴(yán)密的退火計(jì)劃為保證的，即足夠高的初始溫度、緩慢的退火速度、大量的迭代次數(shù)及同一溫度下足夠的擾動(dòng)次數(shù)［5］.因此效率低下一直是SA存在的主要問(wèn)題，針對(duì)這一問(wèn)題學(xué)者們也提出了很多改進(jìn)算法［6-8］，使SA逐步從理論走向?qū)嶋H應(yīng)用.在改進(jìn)算法中，高溫階段模型擾動(dòng)遍歷解空間的能力是確保SA最終能夠找到全局最優(yōu)解的關(guān)鍵因素之一.在貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)中，學(xué)者們也提出了許多基于SA的學(xué)習(xí)算法［9-13］，但是算法在高溫階段模型的全局?jǐn)_動(dòng)能力還有待進(jìn)一步驗(yàn)證.鑒于這一問(wèn)題，本文提出將K2算法和SA算法相結(jié)合進(jìn)行貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí).本算法(K2SA)的主體為模擬退火，只是對(duì)貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)模型的擾動(dòng)方式進(jìn)行了修改，具體為隨機(jī)交換當(dāng)前節(jié)點(diǎn)序中2個(gè)節(jié)點(diǎn)的位置以產(chǎn)生新的節(jié)點(diǎn)序，然后利用K2算法得到網(wǎng)絡(luò)結(jié)構(gòu)，再對(duì)此網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行評(píng)分，作為SA的能量函數(shù)值繼續(xù)進(jìn)行SA的過(guò)程.

1 K2算法與模擬退火算法分析

1.1 K2 算法分析

K2算法是由Cooper和Herskovits于1992年提出的基于評(píng)分搜索的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)算法［3］，是貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)領(lǐng)域的經(jīng)典算法，具有非常優(yōu)異的學(xué)習(xí)性能.

我們知道貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)的目的就是要通過(guò)對(duì)數(shù)據(jù)集D的分析，找到與D吻合最好的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)BS，即后驗(yàn)概率P(BS最大的網(wǎng)絡(luò)結(jié)構(gòu).因?yàn)镻(D)為常數(shù)，所以，由P(BS=P(BS，D)/P(D)知，就是要找到使得概率 P(BS，D)最大的網(wǎng)絡(luò)結(jié)構(gòu)BS.即

K2算法就是以式(1)作為評(píng)分函數(shù)來(lái)學(xué)習(xí)貝葉斯網(wǎng)絡(luò)結(jié)構(gòu).從式(1)可以看出，評(píng)分函數(shù)是可分解的，即它可以看作是每個(gè)節(jié)點(diǎn)Xi，i=1，2，…，n與其父節(jié)點(diǎn)集πi構(gòu)成的n個(gè)局部結(jié)構(gòu)的評(píng)分的乘積.令

則只要最大化每個(gè)局部結(jié)構(gòu)的評(píng)分g(Xi，πi)即可最大化G(BS，D)，也必然最大化整個(gè)貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)的評(píng)分(式(1)).依據(jù)這一思想，在給出節(jié)點(diǎn)序ρ以及每個(gè)節(jié)點(diǎn)父節(jié)點(diǎn)個(gè)數(shù)上限u的情況下，K2利用貪婪搜索依次為每個(gè)節(jié)點(diǎn)找到其父節(jié)點(diǎn)集，從而最終構(gòu)建出完整的貝葉斯網(wǎng)絡(luò).具體方法為，對(duì)每一個(gè)節(jié)點(diǎn) Xi，i=1，2，…，n，從節(jié)點(diǎn)序中排在它之前的節(jié)點(diǎn)集中不斷地找出節(jié)點(diǎn)加入Xi的父節(jié)點(diǎn)集πi，使得Xi與πi構(gòu)成的局部結(jié)構(gòu)的評(píng)分g(Xi，πi)不斷增加，直到πi加入節(jié)點(diǎn)不再增加局部結(jié)構(gòu)的評(píng)分為止，在此過(guò)程中，限制每個(gè)節(jié)點(diǎn)的父節(jié)點(diǎn)的數(shù)量不超過(guò)u.

1.2 模擬退火算法分析

模擬退火的思想最早是由Metropolis等人于1953年提出的，1983年Kirkpatrick等將其用于組合優(yōu)化［4，14］.SA 的出發(fā)點(diǎn)是基于統(tǒng)計(jì)物理學(xué)中對(duì)固體物質(zhì)的退火過(guò)程與組合優(yōu)化問(wèn)題之間的相似性.SA從某一個(gè)初始溫度開始，隨著溫度參數(shù)的不斷降低，結(jié)合Metropolis接受概率的突跳特性，在解空間中隨機(jī)搜索目標(biāo)函數(shù)的全局最優(yōu)解.SA的求解過(guò)程如下:

1)設(shè)定迭代初值k=0，設(shè)定初始溫度t=T0，隨機(jī)產(chǎn)生初始狀態(tài)s=s0，計(jì)算狀態(tài)s的目標(biāo)函數(shù)值E=E(s)；

2)由當(dāng)前狀態(tài)s隨機(jī)產(chǎn)生新狀態(tài)sn，計(jì)算狀態(tài)sn的目標(biāo)函數(shù)值En=E(sn)；

3) 產(chǎn)生［0，1］之間的隨機(jī)數(shù) r，若 r＜min{1，exp［－(En－E)/t］}，則 s=sn；

4)若不滿足等溫搜索終止條件，則轉(zhuǎn)2)；

5)若不滿足算法終止條件，則降溫，令k=k+1，t=Tk，轉(zhuǎn) 2)；

6)輸出算法搜索結(jié)果s.

SA算法中，退火溫度t控制著搜索過(guò)程向著目標(biāo)函數(shù)值減小的優(yōu)化方向前進(jìn)，同時(shí)又以接受概率exp［－(En－E)/t］來(lái)接受劣化狀態(tài)，SA 也因此可以跳出局部最優(yōu).從理論上說(shuō)，只要初始溫度足夠高，退火速度足夠慢，終止溫度較低以及各等溫溫度下足夠多次的模型擾動(dòng)，SA算法最終就能收斂到全局最優(yōu).這也是造成SA算法效率低下的原因.

2 K2SA算法

從上文的分析可以看出，SA算法最終能夠收斂到全局最優(yōu)是以高溫時(shí)模型擾動(dòng)遍歷解空間的能力，嚴(yán)密的退火計(jì)劃的相互配合為前提的.理論上，SA算法要求模型擾動(dòng)能夠幾乎遍歷整個(gè)解空間，因此若采用增加邊、刪除邊和轉(zhuǎn)向邊等方法對(duì)當(dāng)前貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行局部修改產(chǎn)生新的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)，則模型擾動(dòng)的遍歷能力值得懷疑.為了解決這個(gè)問(wèn)題，本文將K2與SA算法相結(jié)合，提出了K2SA貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)算法.算法流程如圖1所示.

1)模型擾動(dòng)方式

圖1 K2SA算法流程

K2SA采用如下方法產(chǎn)生新的網(wǎng)絡(luò)結(jié)構(gòu):在當(dāng)前節(jié)點(diǎn)序中，隨機(jī)選取2個(gè)節(jié)點(diǎn)，交換它們的位置，其他節(jié)點(diǎn)位置不變，產(chǎn)生網(wǎng)絡(luò)節(jié)點(diǎn)的一個(gè)排列ρ，并將其作為新的節(jié)點(diǎn)序，然后利用K2算法進(jìn)行結(jié)構(gòu)學(xué)習(xí)生成新的網(wǎng)絡(luò)結(jié)構(gòu).每個(gè)節(jié)點(diǎn)的父節(jié)點(diǎn)個(gè)數(shù)無(wú)限制，即在節(jié)點(diǎn)序ρ中，排在一個(gè)節(jié)點(diǎn)之前的所有節(jié)點(diǎn)都是該節(jié)點(diǎn)的候選父節(jié)點(diǎn).這里將網(wǎng)絡(luò)結(jié)構(gòu)的貝葉斯評(píng)分一并輸出，并對(duì)其取負(fù)值作為目標(biāo)函數(shù)值(因?yàn)榕c樣本數(shù)據(jù)吻合越好的網(wǎng)絡(luò)結(jié)構(gòu)的貝葉斯評(píng)分越高).

2)初始溫度的確定

初始溫度足夠高，可以使幾乎所有產(chǎn)生的候選結(jié)構(gòu)都能被接受，從而保證算法最終能夠收斂到全局最優(yōu)，但是這樣花費(fèi)的計(jì)算時(shí)間又將會(huì)增加.因此，初始溫度的確定應(yīng)折中考慮優(yōu)化質(zhì)量和優(yōu)化效率.K2SA采用如下方法確定初始溫度:隨機(jī)產(chǎn)生一組狀態(tài)，確定兩兩狀態(tài)間目標(biāo)函數(shù)之差的最大值，再由函數(shù)T0=－/lnPi確定初始溫度，其中Pi為初始接受概率.

3)溫度更新方法

單純加快溫度下降速度并不能保證算法以較快的速度收斂到全局最優(yōu)，溫度下降的速度必須與狀態(tài)產(chǎn)生函數(shù)相匹配.通常，各溫度下模型擾動(dòng)的次數(shù)越多，溫度可以下降越快.為了提高算法的效率，采用指數(shù)形式的退火計(jì)劃tk=T0αk.

4)等溫搜索終止條件

等溫搜索終止條件也稱為Metropolis抽樣穩(wěn)定準(zhǔn)則.SA的收斂性要求在每個(gè)溫度下產(chǎn)生候選解的個(gè)數(shù)趨于無(wú)窮大，以使相應(yīng)的馬氏鏈達(dá)到平穩(wěn)概率分布，在實(shí)際中這顯然不可行.K2SA采用以下等溫搜索終止條件:連續(xù)若干步等溫?cái)_動(dòng)目標(biāo)函數(shù)值變化小于一定的閾值；擾動(dòng)一定的次數(shù).以上2條之一滿足即終止該溫度下的等溫搜索.

5)SA算法終止條件

K2SA中的SA過(guò)程的終止由以下條件控制:達(dá)到最大迭代次數(shù)；連續(xù)若干步降溫得到的最優(yōu)目標(biāo)函數(shù)值變化小于一定的閾值.以上條件之一滿足即終止SA搜索過(guò)程.

另外，為了避免模擬退火的搜索過(guò)程中由于執(zhí)行概率接受環(huán)節(jié)而遺失當(dāng)前遇到的最優(yōu)解，K2SA加入了記憶功能，可將搜索過(guò)程中到目前為止遇到的最優(yōu)結(jié)構(gòu)存儲(chǔ)下來(lái).在SA的搜索過(guò)程結(jié)束后，再利用爬山法在SA過(guò)程獲得的最優(yōu)結(jié)構(gòu)附近進(jìn)行局部搜索，以期進(jìn)一步提高解的質(zhì)量.

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

為了驗(yàn)證本文提出的K2SA算法的性能，采用Asia網(wǎng)［15］對(duì)K2，SA和K2SA進(jìn)行了對(duì)比仿真實(shí)驗(yàn).Asia網(wǎng)又稱為Chest-clinic網(wǎng)，是一個(gè)根據(jù)專家知識(shí)構(gòu)建的醫(yī)療診斷網(wǎng)絡(luò)，由8個(gè)節(jié)點(diǎn)組成，每個(gè)節(jié)點(diǎn)都只有是與否2種狀態(tài)，網(wǎng)絡(luò)結(jié)構(gòu)如圖2所示.

圖2 Asia網(wǎng)結(jié)構(gòu)

仿真所用計(jì)算機(jī)為HP Compaq 8000 Elite CMT Business PC，其處理器為四核 Intel Core 2 Quad Q95002.83GHz，4G 內(nèi)存，1T 硬盤，ATI Radeon HD 4650顯卡.K2算法所用的節(jié)點(diǎn)序?yàn)?Visit to Asia，Tuberculosis，Smoking，Lung cancer，Tuberculosis or cancer，X ray result，Bronchitis，Dyspnoea)，父節(jié)點(diǎn)上限為Asia網(wǎng)節(jié)點(diǎn)個(gè)數(shù)8；SA算法采用指數(shù)形式的退火計(jì)劃，初始溫度為10，衰減率為0.999，迭代次數(shù)為10000，等溫搜索次數(shù)為20，算法從一個(gè)隨機(jī)產(chǎn)生的網(wǎng)絡(luò)開始，通過(guò)在對(duì)當(dāng)前網(wǎng)絡(luò)的加邊、減邊和轉(zhuǎn)向邊操作產(chǎn)生的所有網(wǎng)絡(luò)中隨機(jī)選取一個(gè)作為新狀態(tài)；K2SA算法也采用指數(shù)形式的退火計(jì)劃，初始溫度由第2節(jié)所述方法確定，初始接受概率設(shè)為 0.95，溫度衰減率為0.999，等溫搜索終止條件為達(dá)到最大等溫搜索次數(shù)20，或者連續(xù)5次未接受新狀態(tài)，算法終止條件為達(dá)到最大迭代次數(shù)10000，或者連續(xù)10次降溫都未接受新狀態(tài).學(xué)習(xí)所用數(shù)據(jù)樣本根據(jù)Asia網(wǎng)參數(shù)由Monte Carlo法通過(guò)邏輯抽樣生成，在樣本量為500，1000，5000，10000 的情況下分別進(jìn)行了10次仿真實(shí)驗(yàn)，表1所示為10次仿真結(jié)果的平均值.表中，IE、ME和RE對(duì)應(yīng)列分別為學(xué)習(xí)得到的網(wǎng)絡(luò)結(jié)構(gòu)與真實(shí)的Asia網(wǎng)相比增加的邊、缺失的邊和反向的邊的數(shù)量.從表1可以看出，隨著樣本量的增加，K2和K2SA算法得到的網(wǎng)絡(luò)越來(lái)越接近于真實(shí)的網(wǎng)絡(luò)，而SA算法的效果卻沒有明顯的改善.在3種算法中，K2的學(xué)習(xí)效果最好，K2SA次之，SA最差.究其原因，因?yàn)镵2融入了先驗(yàn)知識(shí)(節(jié)點(diǎn)序和父節(jié)點(diǎn)個(gè)數(shù)的上限)，并且采用了較好的啟發(fā)式搜索算法，而且由于已知了節(jié)點(diǎn)序，所以K2的學(xué)習(xí)結(jié)果中沒有出現(xiàn)反向的邊；SA由于采用對(duì)當(dāng)前模型進(jìn)行局部隨機(jī)修改的模型擾動(dòng)方式，候選模型遍歷解空間的能力不足，造成了SA的學(xué)習(xí)效果不理想；K2SA改進(jìn)了SA的模型擾動(dòng)方式，使得候選模型遍歷解空間的能力得到了提高，最終獲得了令人滿意的學(xué)習(xí)效果.

表1 K2，SA和K2SA算法仿真結(jié)果

為了進(jìn)一步說(shuō)明3種算法的效率，表2給出了4種樣本量情況下，10次仿真運(yùn)行時(shí)間的平均值.從表中可以看出，K2算法的效率相當(dāng)高；SA算法為了獲得較好的學(xué)習(xí)結(jié)果，需要進(jìn)行大量的迭代運(yùn)算，效率極低；K2SA算法改進(jìn)了模型擾動(dòng)方式，以及等溫搜索和迭代終止條件，效率要大大高于SA算法，但遠(yuǎn)低于K2算法.

表2 K2，SA和K2SA算法運(yùn)行時(shí)間 s

4 結(jié)語(yǔ)

本文通過(guò)對(duì)節(jié)點(diǎn)序的擾動(dòng)，利用K2算法改進(jìn)了基于模擬退火的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)算法的局部模型擾動(dòng)方式，提高了模型全局?jǐn)_動(dòng)能力，從而使得K2SA算法能夠獲得近似全局最優(yōu)貝葉斯網(wǎng)絡(luò)結(jié)構(gòu).同時(shí)，由于模擬退火算法本身存在效率不高以及K2SA算法模型結(jié)構(gòu)產(chǎn)生策略的問(wèn)題，K2SA算法的效率略顯不足.

References)

［1］范敏，黃席樾，石為人，等.一種改進(jìn)的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)算法［J］.系統(tǒng)仿真學(xué)報(bào)，2008，20(17):4613-4617.Fan Min，Huang Xiyue，Shi Weiren，et al.Improved Bayesian networks structure learning algorithm［J］.Journal of System Simulation，2008，20(17):4613-4617.(in Chinese)

［2］Cheng J，Greiner R，Kelly J，et al.Learning Bayesian networks from data:an information-theory based approach［J］.Artificial Intelligence，2002，137(1/2):43-90.

［3］Cooper F G，Herskovits E.A Bayesian method for the induction of probabilistic networks from data［J］.Machine Learning，1992，9(4):309-347.

［4］Kirkpatrick S，Gelatt C Jr，Vecchi M.Optimization by simulated annealing［J］.Science，1983，220:671-680.

［5］陳華根，吳健生，王家林，等.模擬退火算法機(jī)理研究［J］.同濟(jì)大學(xué)學(xué)報(bào):自然科學(xué)版，2004，32(6):802-805.Chen Huagen，Wang Jiansheng，Wang Jialin，et al.Mechanism study of simulated annealing algorithm［J］.Journal of Tongji University:Natural Science，2004，32(6):802-805.(in Chinese)

［6］Gelfand S B.Analysis of simulated annealing for optimization［D］.Cambridge:Massachusetts Institute of Technology，1987.

［7］Ingber L.Very fast simulated annealing［J］.Math Conput Modeling，1989，12:967-973.

［8］陳華根，李麗華，許惠平，等.改進(jìn)的非?？焖倌M退火算法［J］.同濟(jì)大學(xué)學(xué)報(bào):自然科學(xué)版，2006，34(8):1121-1125.Chen Huagen，Li Lihua，Xu Huiping，et al.Modified very fast simulated annealing algorithm［J］.Journal of Tongji University:Natural Science，2006，34(8):1121-1125.(in Chinese)

［9］Singh M，Valtorta M.Construction of Bayesian network structures from data:a brief survey and an efficient algorithm［J］.International Journal of Approximate Reasoning，1995，12(2):111-131.

［10］Elidan G，Ninio M，F(xiàn)riedman N，et al.Data perturbation for escaping local maxima in learning［C］//The Eighteenth National Conference on Artificial Intelligence.Edmonton，Alberta，Canada，2002:132-139.

［11］Hartemink A J，Gifford D K，Jaakkola T S，et al.Combining location and expression data for principled discovery of genetic regulatory network models［C］//Pacific Symposium on Biocomputing. Kauai， HI，USA，2002，7:437-449.

［12］Chickering D M，Geiger D，Heckerman D.Learning Bayesian networks:search methods and experimental results［C］//Proceedings of the Fifth International Workshop on Artificial Intelligence and Statistics.Fort Lauderdale，F(xiàn)L，USA，1995:112-128.

［13］張少中，王秀坤，丁華.基于模擬退火的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)算法［J］.計(jì)算機(jī)科學(xué)，2002，31(10):196-208.Zhang Shaozhong，Wang Xiukun，Ding Hua.An algorithm for Bayesian networks structure learning based on simulated annealing［J］.Computer Science，2002，31(10):196-208.(in Chinese)

［14］王凌.智能優(yōu)化算法及其應(yīng)用［M］.北京:清華大學(xué)出版社，2001:17-35.

［15］Lauritzen S L，Spiegelhalter D J.Local computations with probabilities on graphical structures and their application to expert systems［J］.Journal of the Royal Statistical Society B，1988，50(2):157-194.