基于信念傳播的分布式最大權(quán)匹配算法

張源

（東南大學(xué) 移動通信國家重點實驗室，江蘇 南京 210096）

1 引言

考慮一個一般的無向圖G=(V, E)，其中V代表所有節(jié)點的集合，E代表所有邊的集合?？紤]E的某個子集，如果其中任意2條邊都沒有交點，則稱該子集為圖G的一個匹配（match）。假設(shè)每條邊e∈E都被賦予一個正權(quán)重 we，則本文研究最大權(quán)匹配（MWM）問題，即尋找具有最大權(quán)重和的匹配。在形式上，MWM問題可以描述為下列二值線性整數(shù)規(guī)劃問題：

其中，xe=1或0分別代表邊e是或不是匹配的邊，N(i)是所有與節(jié)點i相連的邊構(gòu)成的集合。

最大權(quán)匹配問題是一個非?；镜臄?shù)學(xué)問題[1]，并且能被用于解決大量實際問題。在無線網(wǎng)絡(luò)領(lǐng)域中，有各種無線資源分配與調(diào)度問題，其核心就是尋找給定圖的最大權(quán)匹配如文獻[2,3]。因此，高效的MWM算法，對于解決無線網(wǎng)絡(luò)中的資源調(diào)配問題是非常重要的。事實上，在圖論、組合優(yōu)化等數(shù)學(xué)領(lǐng)域中針對該問題的研究已經(jīng)相當(dāng)充分與完整，并已提出了大量算法[4,5]。然而，這些算法都是中心式的，即需要有某個中心節(jié)點知道整個圖的全局信息。在無線網(wǎng)絡(luò)中，一般來說是不存在這種能掌握全局信息的中心節(jié)點的，而必須采用能夠分布式執(zhí)行的最大權(quán)匹配算法。因此，這些算法是無法直接應(yīng)用于無線網(wǎng)絡(luò)中。

在文獻[6,7]中，已經(jīng)提出了幾種分布式MWM算法。這些算法都是基于貪婪（greedy）原則，非常簡單，其找出的匹配質(zhì)量一般較低（本文第4節(jié)中的仿真結(jié)果也說明了這一點）。近年來，又出現(xiàn)了一類全新的基于信念傳播（belief propagation）的方法。該方法來自人工智能與機器學(xué)習(xí)領(lǐng)域，并且已被成功應(yīng)用于信道編碼問題[8]。文獻[9,10]將該方法應(yīng)用于MWM問題，并提出了一類基于信念傳播的分布式MWM算法。該算法內(nèi)容如下：首先，該算法以迭代方式運行。在每次迭代過程中，每條邊f(xié)=(k,i)向鄰邊e=(i,j)發(fā)送消息mf→e，其計算公式為

其中，符號(x)+代表取x與0中的大者。根據(jù)所有接收到的消息，每條邊e計算度量。

根據(jù)該度量，每條邊e進行如下判斷

其中，問號代表不確定情形。下文將稱該算法為MWM問題的BP算法。文獻[9,10]已證明該算法當(dāng)圖G為二分圖（bipartite）時等價于拍賣（auction）算法[11]，并且當(dāng)MWM唯一時保證收斂至最優(yōu)解。然而，當(dāng)圖G是一般含圈圖時，則只證明了當(dāng)式(1)的松弛線性規(guī)劃沒有分?jǐn)?shù)解時才能收斂至最優(yōu)解。值得一提的是，針對最大權(quán)獨立子集（MWIS）問題的類似算法在文獻[12]中也得到了研究。

2 BP算法

本文主要分析BP算法的不足，揭示其產(chǎn)生的原因，然后針對這些不足對原有BP算法進行改進。本節(jié)提及的各命題及證明統(tǒng)一在第3節(jié)中給出。事實上，BP算法主要存在兩方面的不足。首先，BP算法非常容易出現(xiàn)振蕩現(xiàn)象，從而導(dǎo)致不收斂。為此，引入了信念傳播圈的概念，并且在引理1中給予證明，式(2)更新消息時，如果圖中存在信念傳播圈，那么該圈中相鄰邊間的消息有可能在二值間振蕩取值。這就是BP算法中的振蕩現(xiàn)象。第4節(jié)中的仿真結(jié)果也展現(xiàn)了這一點。為了消除振蕩，建議修改式(2)中的消息計算公式為

從而得到一個改進的BP算法。進一步地，在定理1中證明了，如果根據(jù)式(5)更新消息，則可以消除振蕩問題。

除了振蕩現(xiàn)象，BP算法還可能出現(xiàn)所謂不確定現(xiàn)象。該現(xiàn)象的起因是，當(dāng)邊e根據(jù)式(4)進行判決時，如果be=0，應(yīng)該如何處理？如果判決xe=1，則最終解可能不正確，即出現(xiàn)匹配邊相鄰的情形；如果判決為xe=0，則最終解可能不是最優(yōu)的。這就是BP算法中的不確定現(xiàn)象。文獻[9,10]中并沒有研究這一問題。本文首先在引理4中證明，如果be＞0，所有與e相鄰的邊f(xié)都有bf＜0，不會出現(xiàn)bf=0的情形。接下來，在引理5中證明，如果be=0，則所有與e相鄰的邊f(xié)都有bf≤0，并且邊e每側(cè)有且僅有一條邊的度量為0，即不確定。綜合這2條引理，可以說明，在圖中所有的不確定邊一定構(gòu)成一個圈。在定義3中稱這樣的圈為不確定圈。在引理7中證明，在不確定圈上執(zhí)行信念傳播算法，仍然得到不確定圈。因此，BP算法是不適用于不確定圈的，必須采用新方法處理不確定情形。根據(jù)引理8，建議采用下列簡單的方法。如果邊e發(fā)現(xiàn)自己是不確定的，即be=0，則邊e要沿著不確定圈發(fā)送探針（probe），分別將所有奇數(shù)邊與偶數(shù)邊的權(quán)重和保存在不同的變量中，即

然后返回邊e，其中，邊1為邊e，邊2為邊e順時針一側(cè)的鄰邊，依次類推。假設(shè)不確定圈包括n條邊。如果n=2k，則根據(jù)下式進行判決

如果n=2k+1，則需要進一步計算

其中，w2k+1是邊e逆時針一側(cè)的鄰邊的權(quán)重，然后根據(jù)下式進行判決

綜上，BP算法主要在振蕩與不確定2個方面存在不足。針對這些不足，改進BP算法將主要包括如下2個階段。在第1階段，相鄰邊之間根據(jù)式(5)交換消息直至收斂。在第2階段，每條邊根據(jù)式(3)計算度量，如果 be＞0則判決 xe=1，如果be＜0則判決xe=0，如果be=0，則根據(jù)式(7)或式(9)進行判決。最后，在定理2中證明了該改進算法是正確的。計算機仿真實驗結(jié)果表明，該改進算法可以顯著提高信念傳播類算法在一般圖中尋找最大權(quán)匹配的能力，具體內(nèi)容將在第4節(jié)中給出。

3 算法分析

定義1 考慮任意邊f(xié)=(k,i)，如果

則稱h*為f在節(jié)點k側(cè)的下游邊。

定義2 考慮圖1中所示的一組邊1，2，…，n，n+1，并且邊k+1是邊k的下游邊（1≤k≤n），如果n+1=1，則稱這n條邊構(gòu)成一個信念傳播圈；如果n+1=?，則邊n一定位于圖的邊緣，并稱這n條邊構(gòu)成一條信念傳播路徑。

圖1 信念傳播圈與路徑

引理1 給定信念傳播圈包括n條邊，如果n為奇數(shù)，則其中任意相鄰邊間消息的取值有可能不收斂，并且是在2個數(shù)值之間振蕩取值。

證明 考慮信念傳播圈中的任意一條消息m1→n。根據(jù)定義2，則該消息的計算公式為

定義函數(shù)

可以通過歸納法證明，當(dāng)n為奇數(shù)時，函數(shù)yn(x)的形狀必為圖2(a)～圖2(d)、圖2(e)和圖2(j)中之一；當(dāng)n為偶數(shù)時，函數(shù)yn(x)的形狀必為圖2(f)～圖2(i)、圖2(e)和圖2(j)中之一。

證明如下：當(dāng) n=1時，y1(x)=(w1-x)+，其形狀為圖 2(a)中所示；當(dāng) n=2時，y2(x)=(w1-(w2-x)+)+=(w1-y1(x))+，當(dāng) w1＞w2時，其形狀為圖 2(g)中所示；當(dāng)w1=w2時，其形狀為圖2(f)中所示；當(dāng)w1＜w2時，其形狀為圖2(h)中所示。假設(shè)n=k時命題成立。當(dāng)n=k+1時，yk+1(x)=(w1-yk(x))+。如果k為奇數(shù)，如果yk(x)的形狀為圖2(a)中所示，依據(jù)w1與yk(x)大小關(guān)系的不同，yk+1(x)的形狀為圖2(f)～圖2(h)中所示之一；如果yk(x)的形狀為圖2(b)中所示，則yk+1(x)的形狀為圖2(e)～圖2(h)中所示之一；如果yk(x)的形狀為圖2(c)中所示，則yk+1(x)的形狀為圖2(e)～圖2(h)和圖2(i)中所示之一；如果yk(x)的形狀為圖2(d)中所示，則yk+1(x)的形狀為圖2(e)～圖2(h)和圖2(i)中所示之一；如果yk(x)的形狀為圖2(e)中所示，則yk+1(x)的形狀為圖2(j)中所示；如果yk(x)的形狀為圖2(j)中所示，則yk+1(x)的形狀為圖2(e)和圖2(j)中所示之一。如果k為偶數(shù)，推導(dǎo)過程是類似的。從而命題得證。

消息 m1→n的取值就是方程 x=yn(x)的解。具體來說，就是根據(jù)xt+1=yn(xt)的迭代方程確定m1→n取值，其中上標(biāo)代表迭代次數(shù)。分3種情形分析如下。

1) 當(dāng)n為偶數(shù)時，yn(x)的形狀為圖2(e)～圖2(j)中之一，該迭代過程一定收斂。

2) 當(dāng)n為奇數(shù)時，如果yn(x)的形狀為圖2(e)和圖2(j)中之一，或者yn(x)的形狀為圖2(a)～圖2(d)中之一，即其中包含有一段斜率為負(fù) 1的斜線段，但方程 x=yn(x)的解不在該斜線段上，則該迭代過程一定收斂。

圖2 消息計算模式

3) 當(dāng)n為奇數(shù)時，如果yn(x)的形狀為圖2(a)～圖2(d)中之一，且方程x=yn(x)的解就落在其中的斜線段上，則可以證明該迭代過程不收斂，且 m1→n在2個數(shù)值之間振蕩取值。記斜線段的自變量范圍為a≤x≤b，方程x=yn(x)的解為x*。假設(shè)迭代從某x0（a≤x0＜x*≤b）開始，則 x1=yn(x0)=2x*-x0，x2=yn(x1)=x0等。這樣，m1→n將在 x0與 x1之間振蕩取值。類似地，如果迭代從某x0（a≤x*＜x0≤b）開始，m1→n也將在2個數(shù)值之間振蕩取值。

綜上，當(dāng)n為奇數(shù)時，信念傳播圈中消息計算過程可能不收斂。如果不收斂，則消息將在2個數(shù)之間振蕩取值。證畢。

引理 2 給定信念傳播路徑，其中任意相鄰邊間消息的取值一定收斂。

證明 假設(shè)信念傳播路徑包括n條邊，考慮其中任意一條消息 me→e-1。根據(jù)定義 2，該消息的計算公式為

該計算過程是確定的，因此一定是收斂的。所以，信念傳播路徑中相鄰邊間消息計算的過程都是收斂的。證畢。

定理1 如果根據(jù)式(5)計算消息，給定信念傳播圈，其中任意相鄰邊間消息的取值一定是收斂的。

證明 假設(shè)信念傳播圈包括n條邊，并考慮其中消息 m1→n。根據(jù)式(5)更新消息，其取值將通過明，分3種情形，分析如下。

1) 當(dāng)n為偶數(shù)時，yn(x)的形狀為圖2(e)～圖2(j)

2) 當(dāng)n為奇數(shù)時，如果yn(x)的形狀為圖2(e)～圖 2(j)中之一，或者 yn(x)的形狀為圖 2(a)～圖 2(d)中之一，且方程x=yn(x)的解不在斜線段上，則類似可知迭代過程收斂。

3) 當(dāng)n為奇數(shù)時，如果yn(x)的形狀為圖2(a)～圖2(d)中之一，且方程x=yn(x)的解就落在其中斜線段上，

引理 3 考慮與邊 f=(k,i)一側(cè)相鄰的任意 2條邊 e=(i,j)與 h=(i,l)，則 mf→e=mf→h。

證明 根據(jù)消息的定義可得。證畢。

引理4 如果邊e的度量be＞0，則所有與e相鄰的邊f(xié)都有 bf＜0。

證明 根據(jù)條件，be=we-mf*→e-mg*→e＞0，其中f*與 g*是分別從兩端與 e相鄰且消息取值最大的邊，如圖 3所示?？紤]邊 e的任一相鄰邊 f，則bf=wf-mh*→f-me*→f，其中，h*與 e*是分別與 f相鄰且消息取值最大的邊。根據(jù)消息計算公式we-mg*→e=me→f， 有 be=me→f-mf*→e＞0 ， 即 me→f＞mf*→e。再根據(jù) e*的含義，則有 me*→f≥me→f＞mf*→e。另一方面，考慮 bf，根據(jù)消息計算公式wf-mh*→f=mf→e，有 bf=mf→e-me*→f。將 me*→f＞mf*→e代入，則有 bf＜mf→e-mf*→e，再根據(jù) f*的含義，則有mf→e≤mf*→e，從而得到 bf＜0。證畢。

圖3 引理4的證明

引理5 如果邊e的度量be=0，則所有與e相鄰的邊f(xié)都有bf≤0；進一步地，如果假設(shè)每條邊的下游邊是唯一的，則邊e每側(cè)有且僅有一條邊的度量為0。

證明 如圖 4所示，be=we-mf*→e-mg*→e=0，其中f*與g*是分別從兩端與e相鄰且消息取值最大的邊。由于 we-mg*→e=me→f，則 be=me→f*-mf*→e=0，即me→f=mf*→e?？紤]邊 e的任一相鄰邊 f，則bf=wf-mh*→f-me*→f，其中，h*與 e*是分別與 f相鄰且消息取值最大的邊。由于wf-mh*→f=mf→e，再考慮到 e*與 f*的含義，則有 bf=mf→e-me*→f≤mf*→e-me*→f= me→f-me*→f≤0。接下來考慮該引理中的第2部分，以邊e的i節(jié)點一側(cè)為例?？紤]邊f(xié)*，其度量為bf*=wf*-mh*→f*-me*→f*，其中，h*與 e*是分別與 f*相鄰且消息取值最大的邊。由于wf*-mh*→f*=mf*→e，則有bf*=mf*→e-me*→f*?？梢宰C明 e*=e。根據(jù) f*的含義可知 mf*→e≥me*→e，根據(jù) be=0 可知 me→f*=mf*→e，因此有 me→f*=mf*→e≥me*→e=me*→f*。再根據(jù) e*的含義可知 me*→f*≥me→f*，從而有 me*→f*=me→f*。由于假設(shè)邊f(xié)*的下游邊是唯一的，則有e*=e。這樣，可以得到bf*=mf*→e-me→f*=0，并且是唯一的。證畢。

圖4 引理5的證明

定義 3 由所有度量為零的邊構(gòu)成的圈為不確定圈。

引理6 不確定圈一定是信念傳播圈。

證明 在不確定圈上，根據(jù)引理5的證明，其中任意2條相鄰邊都互為下游邊，從而符合信念傳播圈的定義。證畢。

引理 7 在不確定圈上執(zhí)行信念傳播算法，仍然得到不確定圈。

證明 在不確定圈上，根據(jù)引理5的證明，其中任意2條相鄰邊都互為下游邊。因此，不論是在原圖中還是單獨在不確定圈上，所執(zhí)行的性能傳播算法是完全相同的，從而仍然得到不確定圈。證畢。

引理8 假設(shè)不確定圈包括n條邊，記為1，2，…，n。如果n為偶數(shù)，即n=2k，則最大權(quán)匹配的權(quán)和為max(w1+w3+…+w2k-1, w2+w4+…+w2k)；如果n為奇數(shù)，即 n=2k+1，則最大權(quán)匹配的權(quán)和為 max(w1+w3+…+w2k-1, w3+w5+…+w2k+1, w2+w4+…+w2k)。

證明 在不確定圈中，可能的匹配方案就是依次選擇邊。具體分析如下。如果n是偶數(shù)，則有2種選擇，即 x2i=1（1≤i≤k），相應(yīng)的權(quán)和為 w1+w3+…+w2k-1；或者 x2i-1=1（1≤i≤k），相應(yīng)的權(quán)和為w2+w4+…+w2k。因此，最大權(quán)匹配的權(quán)和為max(w1+w3+…+w2k-1, w2+w4+…+w2k)。如果n是奇數(shù)，則有3種選擇，即x1=x3=…=x2k-1=1，相應(yīng)的權(quán)和為 w1+w3+…+w2k-1；或者 x2=x4=…=x2k=1，相應(yīng)的權(quán)和為 w2+w4+…+w2k；或者 x3=x5=…=x2k+1=1，相應(yīng)的權(quán)和為w3+w5+…+w2k+1。因此，最大權(quán)匹配的權(quán)和為 max(w1+w3+…+w2k-1, w3+w5+…+w2k+1,w2+w4+…+w2k)。證畢。

定理 2 在改進算法執(zhí)行過程中，所有得到的解都是滿足匹配要求的。

證明 需要證明的是，迭代算法的每次輸出一定滿足匹配要求，即所有被選中的邊都是不相鄰的。在改進算法中，一條邊被選中有2種可能。1）該邊不是不確定的。在這種情況下，根據(jù)引理4可知，所有與該邊相鄰的邊的度量都將嚴(yán)格的小于零，是一定不會被選中的，從而保證了正確性。2）該邊是不確定的。在這種情況下，根據(jù)引理5可知，所有與該邊相鄰的邊的度量或者嚴(yán)格小于零，或者也是不確定的?？紤]與該邊相鄰的不確定邊。由于改進算法中采用式(7)或式(9)處理不確定情形，并能達到引理8中給出的權(quán)和性能，位于同一不確定圈上的不確定邊一定是交替成為匹配邊的，因此被選中的邊一定都是不相鄰的，從而保證了正確性。證畢。

4 仿真

首先考慮比較小型的含圈圖，有4個節(jié)點，并且任意2個節(jié)點間都有邊相連，即一個4階完全圖（complete graph）。圖中每條邊都隨機產(chǎn)生一個0～1之間的正實數(shù)作為權(quán)重，然后分別運行BP算法與改進的BP算法，并記錄下每2條相鄰邊之間發(fā)送消息的迭代計算過程。圖 5中給出了一次典型的運行結(jié)果，從實驗結(jié)果可以看出，BP算法采用式(2)計算消息，消息的取值出現(xiàn)了明顯的振蕩現(xiàn)象，而改進的BP算法由于采用式(5)計算消息，其消息的計算過程是收斂的，從而解決了這一問題。進一步地，統(tǒng)計了每種仿真場景下所有 1 000次仿真的收斂情況，BP算法與改進的BP算法都有可能出現(xiàn)不收斂，其中BP算法有286次不收斂，而改進的BP算法只有42次不收斂，明顯優(yōu)于BP算法。綜合該組仿真結(jié)果可以知道，改進BP算法算法雖然未能完全解決一般圖情形下不收斂的問題，但與BP算法相比已經(jīng)有了非常明顯的改善。因此，該組仿真結(jié)果表明，改進BP算法的收斂性能明顯優(yōu)于BP算法。

接下來研究改進算法的性能。仿真場景如下。假設(shè)在 1km×1km的正方形區(qū)域內(nèi)隨機分布若干節(jié)點，其中，如果任意2個節(jié)點之間距離小于某一門限值，則令這2個節(jié)點間有邊相連，并且每條邊都隨機產(chǎn)生一個0～1之間的正實數(shù)作為權(quán)重，這樣就可以得到一個隨機的加權(quán)圖。在仿真中假設(shè)門限值為0.4km。針對給定圖，將分別運行3種不同的分布式MWM算法。第1種算法是文獻[6]中的基于貪婪的分布式MWM算法；第2種算法是文獻[9,10]中的BP算法；第3種算法就是本文提出的改進BP算法。另外，還將通過求解優(yōu)化模型(1)來計算給定圖中的最大權(quán)匹配，作為分布式算法性能的上限。值得說明的是，對傳統(tǒng)BP算法來說，由于傳統(tǒng)BP算法由于振蕩與不確定問題，很可能不收斂，最終輸出結(jié)果很可能是不匹配的，從而導(dǎo)致經(jīng)常得到超過最優(yōu)解的性能。因此本文是搜集了所有BP算法輸出正確匹配的結(jié)果（如圖6所示）。

在仿真中分別假設(shè)節(jié)點總數(shù)在10～30間取值。在不同的節(jié)點總數(shù)下，隨機生成100組不同的節(jié)點分布并得到相應(yīng)的圖，對每張圖分別運行前述 2種分布式MWM算法及最優(yōu)化方法，然后把100次仿真的結(jié)果平均起來作為該節(jié)點總數(shù)下的仿真結(jié)果。從圖6中可以看出，本文所提出的改進算法的性能與最優(yōu)化方法幾乎是相同的，與傳統(tǒng)的BP算法或貪婪算法相比性能改進是顯著的。這組仿真結(jié)果充分說明本文所提出的改進BP算法在性能上的優(yōu)越性。

圖5 消息的迭代計算過程

圖6 分布式匹配算法性能比較

5 結(jié)束語

本文研究了一般圖的最大權(quán)匹配的分布式算法。以基于信念傳播的分布式算法為基礎(chǔ)，分析了原有算法的振蕩問題并提出一種改進的相鄰節(jié)點間交換消息的方法，分析了原有算法中的不確定問題，并提出了一種新的處理方法以保證算法的正確性，從而得到一種新的改進算法。仿真結(jié)果表明，與已有算法相比，改進算法可以收斂至與最優(yōu)解非常接近的程度，這種算法具有更好的收斂與權(quán)重和性能。

[1] DASGUPTA S, PAPADIMITRIOU C, VAZIRANI U. Algorithms[M].New York, McGraw Hill, 2008.

[2] TASSIULAS L, SARKAR S. Maximin fair scheduling in wireless ad hoc networks[J]. IEEE Journal on Selected Areas in Communications,2005, 23(1): 163-173.

[3] CHEN L, LOW S H, CHIANG M, et al. Cross-layer congestion control, routing and scheduling design in ad hoc wireless networks[A].IEEE INFOCOM 2006[C]. Barcelona, Spain, 2006.1-13.

[4] COOK W J, CUNNINGHAM W E, PULLEYBLANK W R, et al. Combinatorial Optimization[M]. New Jersey: Dover Publications, 1998.

[5] SCHRIJVER A. Combinatorial Optimization-Polyhedra and Efficiency[M]. Berlin: Springer, 2004.

[6] HOEPMAN J. Simple distribute weighted matchings[EB/OL].http://arxiv.org /abs/cs/0410047, 2004.

[7] HOEPMAN J, KUTTEN S, LOTKER Z. Efficient distributed weighted matchings on trees[A]. Structural Information and Communication Complexity[C]. Berlin, Springer, 2006.

[8] KSCHISCHANG F R, FREY B J, LOELOGER E A. Factor graph and the sum-product algorithm[J]. IEEE Transactions on Information Theory, 2001, 47(2): 498-519.

[9] BAYATI M, SHAH D, SHARMA M. Max-product for maximum weight matching: convergence, correctness, and LP duality[J]. IEEE Transactions on Information Theory, 2008, 54(3): 1241-1251.

[10] SANGHAVI S, MALIOUTOV D, WILLSKY A S. Belief propagation and LP relaxation for weighted matching in general graphs[J]. IEEE Transactions on Information Theory, 2011, 57(4): 2203-2212.

[11] BERTSEKAS D P. Auction algorithm for network flow problems: a tutorial introduction[J]. Computational Optimization and Applications,1992, 1(1): 7-66.

[12] SANGHAVI S, SHAH D, WILLSKY A S. Message passing for maximum weight independent set[J]. IEEE Transactions on Information Theory, 2009, 55(11): 4822-4834.