任東平 郭建喜 郝小禮 蔣濤

1.海軍勤務學院基礎(chǔ)部；2.海軍裝備部裝備保障大隊

變量消元法（Variable Elimination， VE）是貝葉斯網(wǎng)絡眾多推理算法中最基本的一個，其推理的快慢和復雜度主要取決于消元的順序。尋找最優(yōu)消元順序是一個非確定性多項式難解算法（Nondeterminism Polynomial Hard， NP-Hard）問題，在實際中常采用啟發(fā)式搜索來求解。為了提高變量消元法的推理速度，在此對最小度、最大勢、最小缺邊和最小增加復雜度搜索方法進行了研究，以亞洲網(wǎng)絡為例，分析計算了上述搜索方法的復雜度和消元順序，通過MATLAB R2018a對上述不同搜索方法分別進行網(wǎng)絡構(gòu)建和推理，最后通過推理時間分析比較了4種搜索方法的性能。實驗結(jié)果表明最小增加復雜度搜索方法優(yōu)于其他搜索方法，其平均耗時最少為0.012s，可加快貝葉斯網(wǎng)絡的推理過程。

貝葉斯網(wǎng)絡源于Pearl[1]提出的不確定知識表示模型，因其有著概率論的嚴謹理論基礎(chǔ)，以圖模型清晰地表示隨機事件之間的相互依賴關(guān)系，成為人工智能、模式識別和故障診斷等領(lǐng)域的研究熱點[2]。貝葉斯網(wǎng)絡推理主要有兩種方法：精確推理和近似推理。精確推理主要有變量消元法（Variable Elimination）、團樹傳播算法（Clique Tree Propagation）、多樹傳播（Polytree Propagation）和圖約簡算法（Graph Reduction）等。其中，變量消元法是最基本的推理算法，也較容易理解，其優(yōu)勢在于它的通用性和簡易性，能夠很好地解決復雜貝葉斯網(wǎng)絡的推理，尤其是復雜的故障診斷貝葉斯網(wǎng)。

變量消元法的復雜性主要取決于其消元的先后順序，尋找最優(yōu)消元順序是一個NP-難解問題[3]，目前普遍采用近似算法求解，主要有最小度搜索[4]、最大勢搜索[2，5]、最小缺邊搜索[2]，向光軍通過對變量消元法的分析，提出了最小缺邊搜索（Minimum Deficiency Search， MDS）和變量消元并行運行的算法[6]，但該算法沒有找到最優(yōu)的消元順序，仍存在一定的缺陷；高文宇提出了一種新的搜索算法——最小增加復雜度搜索[7]，但只是通過仿真實驗進行說明，而沒有用到實例中，未給出其搜索方法的具體實現(xiàn)過程，也沒有和上述搜索方法進行對比。針對以上問題，本文以貝葉斯網(wǎng)絡經(jīng)典模型亞洲網(wǎng)絡來驗證上述搜索算法在變量消元法中的優(yōu)劣。

1 貝葉斯網(wǎng)絡

貝葉斯網(wǎng)絡（Bayesian Network）是一個有向無環(huán)圖，其中的節(jié)點表示事件中的隨機變量，有向邊則表示節(jié)點之間的相互依賴關(guān)系。每個節(jié)點都有概率分布，根節(jié)點的概率分布屬于邊緣分布，也可稱為先驗概率；非根節(jié)點則是條件概率。

在概率論中，可以用條件概率鏈[8]的形式表示其聯(lián)合概率，其形式如式（1）所示：

公式（1）也可以稱為鏈規(guī)則，其中Y為隨機變量，P（Y1,Y2,… ,Yk）是Y1,Y2,… ,Yk的聯(lián)合 概率，P（Yi|Yi-1, … ,Y1）是Yi,Yi-1, … ,Y1的條件概率。

將貝葉斯網(wǎng)絡條件獨立性假設(shè)用于鏈規(guī)則公式（1）可得到如式（2）所示：

其中P（Yi|Parent（Yi））為貝葉斯網(wǎng)絡節(jié)點Yi的條件概率，Parent（Yi）表示Yi的直連父節(jié)點。

為了下文論述，在這里引入貝葉斯網(wǎng)絡的一個經(jīng)典模型——亞洲網(wǎng)絡[9]。抽煙可能會引起肺癌或者支氣管炎，出訪亞洲可能會導致肺結(jié)核，這3種疾病的任何一種都有可能會導致呼吸困難。若有肺結(jié)核或者肺癌，X光胸透結(jié)果可能呈陽性，把這些因果關(guān)系結(jié)合起來，便得出了經(jīng)典的亞洲網(wǎng)絡模型，如圖1所示，其中V表示到過亞洲，T表示感染肺結(jié)核，S表示抽煙，L表示患有肺癌，B表示患有支氣管炎，E表示感染肺結(jié)核或者患有肺癌，X表示拍攝X光，D表示呼吸困難。用公式（2）表達如圖1所示中的各個節(jié)點的聯(lián)合概率分布可得到如式（3）所示：

圖1 亞洲網(wǎng)絡模型Fig.1 Asian network model

如表1所示為亞洲網(wǎng)絡節(jié)點概率，其中P（i）表示節(jié)點i的先驗概率，P（i|j）表示在j為真時節(jié)點i的條件概率，）表示在j為假時節(jié)點i的條件概率。

表1 亞洲網(wǎng)絡節(jié)點概率Tab.1 Probability of asian network nodes

2 變量消元法

假設(shè)F是關(guān)于變量集X={X1，X2，…，Xn}的一個函數(shù)，而f={f1，f2，…，fm}也是一組函數(shù)，其中f所涉及的變量都是變量集X的一個子集。如果則f為F的一個分解，f1，f2，…，fm稱為F的分解因子。消元方式可以分為兩種。一種是從F={X1，X2，…，Xn}出發(fā)，通過的方式完成X1的消元；另一種從f={f1，f2，…，fm}出發(fā)，將與X1的所有相關(guān)函數(shù){f1，f2，…，fk}，利用公式行消元，并將返回到f中，就實現(xiàn)了X1的消元。

通過以上敘述，若采用第一種消元方式，其計算的復雜度相對于變量的個數(shù)n成指數(shù)增長；若采用第二種消元方法，其計算的復雜度僅僅依賴于涉及X1的分解因子f1,f2,… ,fk的變量數(shù)量，其計算復雜度會大大減小。

假設(shè)用X來表示一個貝葉斯網(wǎng)絡N中所有變量節(jié)點的集合，用f表示N中所有節(jié)點概率分布的集合，依據(jù)貝葉斯網(wǎng)絡的定義，則f可視為N表示的聯(lián)合概率分布P(X)的一個分解。假設(shè)觀測到了證據(jù)集E={E1,E2, …,Em}，它們的取值記為e。在集合f的因子中，將所有表示證據(jù)的變量設(shè)置為實際觀測值，就可以得到另一組函數(shù)，記之為f ′。這一步驟稱為設(shè)置證據(jù)，f ′是函數(shù)P(Y，E=e)的一個分解，這里Y=X∩。

設(shè)Y的一個子集為Q，f ′從中依次消去所有在Y中，但不在Q中的變量，就可以得到另一個函數(shù)集合，記為f ′′。f ′′是P(Q，E=e)的一個分解，將f ′′中的所有因子相乘便得到P(Q，E=e)。按照條件概率的定義可得到如式(4)所示：

VE算法需要提供5個輸入量：（1）貝葉斯網(wǎng)絡N，同時也是聯(lián)合分布P（X）的一個分解；（2）證據(jù)變量的集合E；（3）證據(jù)變量集合E的取值e；（4）查詢變量的集合Q；（5）所有不在Q∪E中的變量的排序ρ，即為消元的順序。對于聯(lián)合分布P（X）的獲取采用公式（2），只需按照網(wǎng)絡拓撲圖寫出聯(lián)合分布即可；證據(jù)集變量集合E以及E的取值e均由人員觀測得到；查詢變量Q即為所求事件；消元順序ρ的獲取需要通過搜索方法來獲取，不同的搜索方法獲取不同的消元順序，以上即為變量消元法的輸入量。VE推理算法的代碼實現(xiàn)如表2所示。

表2 VE算法的實現(xiàn)過程Tab.2 Implementation process of VE algorithm

3 消元復雜度分析

從變量消元法的實現(xiàn)過程可以看出，最耗費時間和內(nèi)存的步驟是對Elim(f，Z)這一函數(shù)的調(diào)用。這一函數(shù)的運算量遠超其他步驟的運算，所以說整個變量消元法的計算復雜度主要是由這一函數(shù)來決定的。

函數(shù)Elim(f，Z)表示從f中挑出所有包含Z的函數(shù){f1，…，fk}，將它們相乘，得到中間函數(shù)G，再從G中消去Z。設(shè)X1，…，Xk是G中所有除Z之外的變量，如果把函數(shù)表示成多維表，則G所需存儲的函數(shù)個數(shù)為綜上所述，是函數(shù)G復雜性的一個度量，從而也是Elim(f，Z)復雜性的一個度量，稱之為變量Z的消元成本，記為δ(Z)[2，6]，如式（5）所示。

用變量消元法進行推理時，需要消去多個變量，設(shè)它們依次為Z1, Z2, …,Zn。在這里用總消元成本來衡量變量消元法的總復雜度。以亞洲網(wǎng)絡為例，函數(shù)集合為f={P（V） ,P（S）,P（T|V） ,P（L|S）,P（B|S）,P（D|B,E）,P（X|E）,P（E|T,L）}。設(shè)證據(jù)節(jié)點為X=1，E=1，計算V的可信度，即求P(V|X=1，E=1)的概率，消元順序為ρ=＜D,B,L,S,T＞，網(wǎng)絡中所有變量均取二值。消去D時，需要計算G=P(D|B，E)，共涉及3個變量：D、B、E，所以D的消元成本δ(D)=2×2×2=8。依次類推，其余4個節(jié)點B，L，S，T的消元成本分別為 4，8，2，4。則按此消元順序進行消元的總成本為26。

所謂結(jié)構(gòu)圖就是考慮一個函數(shù)的集合F={f1，…，fm}，F(xiàn)的結(jié)構(gòu)圖可定義為從一個空圖出發(fā)，對F中的每一個變量，在圖中相應的添加一個節(jié)點，對于任意兩個變量X和Y，如果它們同時出現(xiàn)在同一個因子fi中，則在它們對應的節(jié)點之間添加一條邊。按此規(guī)則構(gòu)建的圖模型為結(jié)構(gòu)圖。如圖2所示為亞洲網(wǎng)絡結(jié)構(gòu)圖。

圖2 亞洲網(wǎng)絡結(jié)構(gòu)圖Fig.2 Asia network structure diagra m

從一組函數(shù)f中消去一個變量Z的成本，能從該函數(shù)的結(jié)構(gòu)圖中算得。如前所述，從f中消去Z的成本其中X1，…，Xl是除Z之外的變量，Z是函數(shù)所涉及到的所有變量，f1，…，fk是所有涉及Z的因子。根據(jù)上文對結(jié)構(gòu)圖的定義，在f的結(jié)構(gòu)圖中，{X1，…，Xl}剛好是所有與Z相鄰節(jié)點的集合，記為nb(Z)，因此有如式（6）所示：

在亞洲網(wǎng)絡中，從f中消去變量T需要計算G=P(T|V)P(E|T，L)，涉及到4個節(jié)點{T，V，E，L}，所以在結(jié)構(gòu)圖中，與T相鄰的節(jié)點使nb（T） ={V,L,E}，所以δ(T)也可以按公式（5）計算，即

4 消元順序構(gòu)造

4.1 最小度法

在結(jié)構(gòu)圖中，最小度搜索法就是把度數(shù)最小的節(jié)點變量放到消元順序隊列的末尾，并在網(wǎng)絡中刪去該節(jié)點和連接該節(jié)點的有向邊。若有多個節(jié)點的度數(shù)相同，則任選其一，重復上述過程，直到消去結(jié)構(gòu)圖中所有的節(jié)點。此消元法稱為最小度方法。以亞洲網(wǎng)絡模型為例，設(shè)證據(jù)節(jié)點為X=1，E=1，其中1表示事件為真。分別求V，T，S，L，B，D的后驗概率。即求P(V|X=1，E=1)，P(T|X=1，E=1)，P(S|X=1，E=1)，P(L|X=1，E=1)，P(B|X=1，E=1)，P(D|X=1，E=1)的概率。采用最小度搜索方法求消元順序，總消元順序為ρ=＜E,D,B,L,S,T,X,V＞。求節(jié)點V時的消元順序為ρ=＜D,B,L,S,T＞，由第二節(jié)中的公式（3）可得其聯(lián)合分布為如式（7）所示：

消去元素D如式（8）所示：

消去元素B如式（9）所示：

消去元素L如式（10）所示：

消去元素S如式（11）所示：

消去元素T如式（12）所示：

計算V的后驗概率分布得到如式（13）所示：

計算V的最大假設(shè)檢驗得到如式（14）所示：

上述過程即為變量消元法的推理計算過程，在MATLAB R2018a中構(gòu)建亞洲網(wǎng)絡，求得P （V= 1 |X= 1 ,E= 1 ） =0.01，

耗時0.0750s（后續(xù)不同搜索方法的推理計算過程均如上，不再贅述）。同樣求T時的消元順序為ρ=＜ D,B,L,S,V ＞；求S時的消元順序為ρ=＜ D,B,L,T,V ＞；求L時的消元順序為ρ=＜ D,B,S,T,V ＞；求B時的消元順序為ρ=＜ D,L,S,T,V ＞ ；求 D 時的消元順序為ρ=＜ B,L,S,T,V ＞?？傻酶鱾€節(jié)點的推理時間分別為0.0239s，0.0127s，0.0164s，0.0146s，0.0264s。其平均推理耗時為0.0282s。

4.2 最大勢搜索

最大勢搜索是根據(jù)以下規(guī)則在無向圖中對所有節(jié)點進行編號：在有n個節(jié)點的無向圖中，在第i步中，選擇擁有最多已編號的相鄰節(jié)點，且這些節(jié)點未編號，將其編號為n-i+1。如果存在多個這樣的節(jié)點，則任選其一。當無向圖中所有節(jié)點都被編完號后，將編號按照由小到大的順序排序，其順序就是對應節(jié)點的消元順序，此為最大勢搜索方法的一般思想。運用最大勢搜索方法求亞洲網(wǎng)絡消元順序：第1步先對B節(jié)點編號為8，第2步可對節(jié)點S或者D編號，這里將D節(jié)點編號為7，第3步對節(jié)點S編號為6，第4步對節(jié)點L編號為5，第5步對E節(jié)點編號為4，第6步可對節(jié)點T或者X編號，這里選擇將T編號為3，第7步將節(jié)點V編號為2，第8步將最后一個節(jié)點X編號為1。至此所有節(jié)點均已編號，其消元順序為ρ=＜ X，V，T，E，L，S，D，B ＞。除去證據(jù)節(jié)點X和E，可的相應節(jié)點的消元順序，求V時的消元順序為ρ=＜ T，L，S，D，B ＞；求T時的消元順序為ρ=＜ V，L，S，D，B ＞；求S時的消元順序為ρ=＜ V，T，L，D，B ＞；求L時的消元順序為ρ=＜ V，T，S，D，B ＞；求B時的消元順序為ρ=＜ V，T，L，S，D ＞；求 D 時的消元順序為ρ=＜ V，T，L，S，B ＞。可得各個節(jié)點的推理時間分別為0.0258s，0.0126s，0.0118s，0.0126s，0.0118s，0.0146s。其平均推理耗時為0.0149s。

4.3 最小缺邊搜索

無向圖中某個節(jié)點Z的缺邊數(shù)就是在用函數(shù)Elim(f，Z)消去Z時需要添加的邊數(shù)。最小缺邊搜索一邊計算節(jié)點的缺邊數(shù)，一邊將節(jié)點從無向圖中刪除。每一步都需計算節(jié)點的缺邊數(shù)，在消去節(jié)點時選擇缺邊數(shù)最小的節(jié)點進行刪除，當同時有多個節(jié)點的缺邊數(shù)相同時，任選其一，如此循環(huán)，直到所有節(jié)點均被刪除。仍以亞洲網(wǎng)絡模型為例，在初始時，節(jié)點{V，S，T，L，B，E，X，D}的缺邊數(shù)分別為 {0，1，2，2，2，8，0，0}，任選其一，這里選擇 V。消去 V后節(jié)點{S，T，L，B，E，X，D}的缺邊數(shù)分別為 {1，0，2，2，8，0，0}，任選其一，這里選擇X。消去X后節(jié)點{S，T，L，B，E，D}的缺邊數(shù)分別為{1，0，2，2，6，0}，任選其一，這里選擇T。消去T后節(jié)點{S，L，B，E，D}的缺邊數(shù)分別為{1，1，3，3，0}，消去D。消去D后節(jié)點{S，L，B，E}的缺邊數(shù)分別為{1，1，1，1}，由于節(jié)點S，L，B，E構(gòu)成一個環(huán)，所以缺邊數(shù)相同，可消去任意一個，這里選擇消去E，其余3個節(jié)點可任意消去，這里選擇的順序為L，S，B。綜上所述，最小缺邊搜索法的消元順序為ρ=＜V，X，T，D，E，L，S，B ＞。除去證據(jù)節(jié)點X和E，可得其他節(jié)點的消元順序，求V時的消元順序為ρ=＜ T，D，L，S，B ＞；求T時的消元順序為ρ=＜ V，D，L，S，B ＞；求 S 時的消元順序為ρ=＜ V，T，D，L，B ＞；求L時的消元順序為ρ=＜V，T，D，S，B ＞；求B 時的消元順序為ρ=＜ V，T，D，L，S ＞；求D時的消元順序為ρ=＜ V，T，L，S，B ＞?？傻酶鱾€節(jié)點的推理時間分別為 0.0150s，0.0122s，0.0119s，0.0119s，0.0121s，0.0125s。其平均推理耗時為0.0126s。

4.4 最小增加復雜度搜索

其主要思想是在結(jié)構(gòu)圖中消元時，刪除某個節(jié)點，與該節(jié)點相連的所有邊均會被刪除，為了降低圖的復雜度，往往是希望刪除的越多越好，在這里將刪去邊的條數(shù)用ed表示；當節(jié)點刪除后，需要將刪除后的圖構(gòu)成一個完備圖，所以會在某些節(jié)點之間添加一些邊，使其構(gòu)成完備圖，但添加的邊又會增加圖的復雜度，所以希望添加的邊越少越好，將添加的邊數(shù)記為ea，用ea/ed可得一個新的衡量指標，即最小增加復雜度，用IC表示。選擇消元順序時，先計算出所有節(jié)點的IC值，選擇IC值最小的進行消元，若有多個相同值，則任選其一，接著計算余下節(jié)點的IC值，如此重復，直到所有節(jié)點被消去，可得出消元順序。最小增加復雜度搜索方法是一個動態(tài)的搜索過程，從整個網(wǎng)絡圖出發(fā)，對其不斷地循環(huán)判斷，以尋得最優(yōu)解，提高推理整體的運行時間，避免了局部消元過快而整體較慢的不足。此算法適用于故障診斷時對多查詢變量的消元，能夠提高整體的運行時間。以亞洲網(wǎng)絡為例，運用最小增加復雜度搜索方法，如表3所示為初始狀態(tài)節(jié)點的IC值。

表3 初始狀態(tài)節(jié)點IC值Tab.3 Initial state node IC value

由表3可知，節(jié)點V，X，D的IC值均為0，可任選其一進行消元，這里選擇節(jié)點V，消去節(jié)點V后可得其余節(jié)點的IC值，如表4所示。

表4 消去節(jié)點V后節(jié)點IC值Tab.4 Node IC value after eliminating node V

由表4可知，節(jié)點T，X，D的IC值均為0，可任選其一進行消元，這里選擇節(jié)點D，消去節(jié)點D后可得其余節(jié)點的IC值，如表5所示。

表5 消去節(jié)點D后節(jié)點IC值Tab.5 Node IC value after eliminating node D

由表5可知，節(jié)點T，X的IC值均為0，可任選其一進行消元，這里選擇節(jié)點X，消去節(jié)點X后可得其余節(jié)點的IC值，如表6所示。

表6 消去節(jié)點X后節(jié)點IC值Tab.6 Node IC value after eliminating node X

由表6可知，節(jié)點T的IC值均為0，消去節(jié)點T后可得其余節(jié)點的IC值，如表7所示。

表7 消去節(jié)點T后節(jié)點IC值Tab.7 Node IC value after eliminating node T

由表7可知，節(jié)點S，L，B，E的IC值均為0.5，可任選其一進行消元，這里選擇節(jié)點S，消去節(jié)點S后余下的節(jié)點L，B，E構(gòu)成一個三角圖，所有節(jié)點都是等效的，可任意組合。最終得出消元順序ρ=＜ V，D，X，T，X，L，B，E ＞。可得各個節(jié)點的推理時間分別為0.0119s，0.0130s，0.0124s，0.0116s，0.0115s，0.0116s。其平均推理耗時為0.012s。如圖3所示為4種推理算法平均耗時對比，可以明顯的看出最小增加復雜度搜索算法運行時間優(yōu)于其他3種搜索方法。

圖3 不同推理算法對比Fig.3 Comparison of different inference algorithms

5 結(jié)語

貝葉斯網(wǎng)絡有著廣泛的應用前景，可用于人工智能、故障診斷等領(lǐng)域，計算速度的快慢直接決定著該方法在上述領(lǐng)域能否較好的應用。故本文以亞洲網(wǎng)絡模型為例，在推理實例中論述了變量消元法的實現(xiàn)過程，通過MATLAB R2018a對亞洲網(wǎng)絡進行構(gòu)造和推理，通過實驗分別對比了變量消元法的4種消元順序構(gòu)造法，實驗表明最小增加復雜度搜索算法優(yōu)于其他搜索算法，可有效減少推理運算的時間。

引用

[1]PEARL J.Fusion,Propagation,and Structuring in Belief Networks[J].Artificial Intelligence,1986(3):241-288.

[2]張連文.貝葉斯網(wǎng)引論[M].北京:科學出版社,2006.

[3]COOPER G F.The Computational Complexity of Probabilistic Inference Using Bayesian Belief Networks[J].Artificial Intelligence,1990,42(2):393-405.

[4]AMESTOY P R,DAVIS T A,DUFF I S.Algorithm 837:AMD,an Approximate Minimum Degree Ordering Algorithm[J].ACM Transactions on Mathematical Software,2004,30(3):381-388.

[5]BERRY A,BLAIR J R S,HEGGERNES P,et al.Maximum Cardinality Search for Computing Minimal Triangulations of Graphs[J].Algorithmica (New York),2004,39(4):287-298.

[6]向光軍,孔兵,歐家欽.貝葉斯網(wǎng)絡VE推理算法的并行化研究[J].云南大學學報(自然科學版),2010,32(4):392-395.

[7]高文宇,張力.貝葉斯網(wǎng)最優(yōu)消元順序的近似構(gòu)造算法[J].計算機應用,2011,31(8):2072-2074.

[8]李儉川.貝葉斯網(wǎng)絡故障診斷與維修決策方法及應用研究[D].長沙:中國人民解放軍國防科學技術(shù)大學,2002.

[9]ZHANG N L.Hierarchical Latent Class Models for Cluster Analysis[C]//National Conference on Artificial Intelligence,2002:230-237.

[10][ZHANG N L,POOLE D.Exploiting Causal Independence in Bayesian Network Inference[J].The Journal of Artificial Intelligence Research,1996.