陳麗敏，楊靜，張健沛

(1.哈爾濱工程大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，黑龍江哈爾濱150001;2.牡丹江師范學(xué)院計(jì)算機(jī)科學(xué)與技術(shù)系，黑龍江牡丹江157011)

信息網(wǎng)絡(luò)普遍存在，如社會(huì)信息網(wǎng)絡(luò)、DBLP書目網(wǎng)絡(luò)，這些網(wǎng)絡(luò)由多種類型數(shù)據(jù)構(gòu)成，不同類型數(shù)據(jù)之間彼此關(guān)聯(lián)，稱之為異構(gòu)信息網(wǎng)絡(luò)。對(duì)異構(gòu)信息網(wǎng)絡(luò)聚類分析能更好地理解網(wǎng)絡(luò)的隱藏結(jié)構(gòu)以及每個(gè)類的數(shù)據(jù)所代表的角色［1］，而基于概率模型的異構(gòu)信息網(wǎng)絡(luò)聚類算法［2-3］只針對(duì)具體的應(yīng)用領(lǐng)域設(shè)計(jì)函數(shù)，不具有普遍性，且收斂性也不穩(wěn)定。異構(gòu)網(wǎng)絡(luò)的鏈接推理［4］時(shí)，往往需要最初的聚類劃分更精確一些，而傳統(tǒng)的高階異構(gòu)聚類算法［5-6］復(fù)雜度太高不適合異構(gòu)信息網(wǎng)絡(luò)。異構(gòu)信息網(wǎng)絡(luò)經(jīng)常是動(dòng)態(tài)的而非靜態(tài)的，動(dòng)態(tài)同構(gòu)數(shù)據(jù)的演化聚類［7］的研究已經(jīng)做了很多，而動(dòng)態(tài)的異構(gòu)信息網(wǎng)絡(luò)的演化聚類分析，目前只有ENetClus［8］算法，該算法沒有關(guān)注聚類質(zhì)量，更側(cè)重跟蹤類的變化及分析類的推進(jìn)、形成與消失。Khoa［9］使用近似 commute time嵌入聚類同構(gòu)數(shù)據(jù)集，取得了很好的效果。受該思想啟發(fā)，本文從相容二部圖的角度提出一種基于嵌入技術(shù)的動(dòng)態(tài)異構(gòu)信息網(wǎng)絡(luò)的演化聚類算法，具有較高的聚類質(zhì)量，且計(jì)算速度也比較快。

1 時(shí)間平滑二部圖

定義 1:給定G=＜V，E，W＞，V=X0∪X1，其中X0與X1為2個(gè)不同類型的數(shù)據(jù)集。若?〈xi，xj〉∈E，則xi∈X0且xj∈X1，稱G為二部圖。

r(xi，xj)表示二部圖G的結(jié)點(diǎn)xi與xj的關(guān)系，若結(jié)點(diǎn)xi與xj有關(guān)系，則結(jié)點(diǎn)xi與xj之間有邊存在，否則無邊存在。

給定t時(shí)刻的二部圖Gt，?xi，xj∈Gt。使用代價(jià)函數(shù)時(shí)間平滑Gt兩結(jié)點(diǎn)的關(guān)系:

式中:rt(xi，xj)表示t時(shí)刻結(jié)點(diǎn)xi與xj的時(shí)間平滑的關(guān)系，rO(xi，xj)表示t時(shí)刻二部圖兩結(jié)點(diǎn)的原始關(guān)系，即沒有時(shí)間平滑的關(guān)系。rt－1(xi，xj)表示先前時(shí)間結(jié)點(diǎn)xi與xj之間的關(guān)系。SC(·)為快照代價(jià)函數(shù)，表示t時(shí)刻rt(xi，xj)與rO(xi，xj)的相似度。TC()為時(shí)間代價(jià)函數(shù)，表示t時(shí)刻rt(xi，xj)與rt－1(xi，xj)的相似度。

函數(shù)SC()和TC()可選擇多種衡量指標(biāo)。取SC()=，使cost最小的rt(xi，xj)就是最佳關(guān)系，則兩結(jié)點(diǎn)最佳的關(guān)系為

例如，t時(shí)刻DBLP書目網(wǎng)絡(luò)的papers與authors構(gòu)成一個(gè)二部圖G，G中只有不同類型的結(jié)點(diǎn)之間存在關(guān)系。t時(shí)刻的G中，先前時(shí)間的papers已經(jīng)全部更換，但先前時(shí)間的authors會(huì)保留在t時(shí)刻的二部圖中，因此t時(shí)刻的二部圖G，僅使用異構(gòu)數(shù)據(jù)的關(guān)系無法體現(xiàn)先前時(shí)間authors之間的聯(lián)系。而先前時(shí)間authors之間的聯(lián)系影響著t時(shí)刻papers的劃分。所以需要表達(dá)先前時(shí)間所有結(jié)點(diǎn)間的關(guān)系。cij表示二部圖結(jié)點(diǎn)i，j的commute time距離，cij能夠表達(dá)二部圖的所有結(jié)點(diǎn)間的關(guān)系。cij是2個(gè)結(jié)點(diǎn)間所有路徑的平均值，故cij能夠表達(dá)一段時(shí)間兩結(jié)點(diǎn)的關(guān)系。令ct(xi，xj)表示t時(shí)刻及先前時(shí)間結(jié)點(diǎn)xi與xj的 commute time 距離。令rt－1(xi，xj)=ct－1(xi，xj)，rO(xi，xj)是t時(shí)刻二部圖Gt中兩結(jié)點(diǎn)的原始關(guān)系。由式(1)獲得t時(shí)刻時(shí)間平滑二部圖，若rt'(xi，xj)≠0，則兩結(jié)點(diǎn)之間存在邊，否則無邊。rt'(xi，xj)≠0的數(shù)目就是的邊的數(shù)目，中2個(gè)相同類型的結(jié)點(diǎn)之間也可能存在邊。充分體現(xiàn)了t時(shí)刻及先前時(shí)間所有結(jié)點(diǎn)間的關(guān)系。

式(1)只計(jì)算t時(shí)刻屬于的結(jié)點(diǎn)間的關(guān)系。先前時(shí)間結(jié)點(diǎn)的關(guān)系會(huì)導(dǎo)致不夠稀疏，通過計(jì)算Gt中也屬于的結(jié)點(diǎn)的k近鄰，可構(gòu)造稀疏的時(shí)間平滑二部圖。由2.2節(jié)可快速計(jì)算的近似commute time嵌入，而由嵌入可計(jì)算任意兩結(jié)點(diǎn)的ct(xi，xj)，因此每次只要存儲(chǔ)t時(shí)刻近似commute time嵌入即可。

2 時(shí)間平滑二部圖的近似commute time嵌入

2.1 時(shí)間平滑二部圖的commute time嵌入

設(shè)G'是加權(quán)無向圖，L是G'的Laplacian矩陣，L+是L的偽逆矩陣，則

性質(zhì) 1［10］:cij=Gvol'(ei－ej)TL+(ei－ej)

其中Gvol'是G'的權(quán)重總和，即Gvol'=∑wij;wij為G'的結(jié)點(diǎn)i、j構(gòu)成的邊的權(quán)重;ei是第i個(gè)元素為1的單位列向量，即

設(shè)Λ與Φ分別是時(shí)間平滑二部圖的Laplacian矩陣L的特征值構(gòu)成的對(duì)角矩陣及特征值對(duì)應(yīng)的特征向量矩陣，特征值λ1≤λ2≤…≤λn。L+是的矩陣L的偽逆矩陣，則由性質(zhì)1，的任意兩結(jié)點(diǎn)i，j的cij為

則cij是空間第i列向量與第j列向量的歐式距離的平方，稱為時(shí)間平滑二部圖的commute time嵌入［11］。

直接計(jì)算ψ需花費(fèi)O(n3)時(shí)間分解特征矩陣。設(shè)有n個(gè)結(jié)點(diǎn)s條邊，定向的邊，令

則 Bs×n是一個(gè)有向邊－點(diǎn)入射矩陣。令是由邊的權(quán)值構(gòu)成的對(duì)角矩陣，則的Laplacian矩陣 L=BTB［11］。

2.2 時(shí)間平滑二部圖的近似commute time嵌入

其中，Qkr×s是行向量獨(dú)立同分布的隨機(jī)矩陣，

但計(jì)算Y，涉及L+，直接計(jì)算 L+復(fù)雜度過高。根據(jù)文獻(xiàn)［12］方法，可分解計(jì)算Y。令()，則 Y=θL+，等價(jià)于計(jì)算 YL=θ。通過矩陣θ的每個(gè)行向量θi計(jì)算方程組yiL=θi，其中yi是矩陣Y的行向量。使用STSolve求解程序［13］能夠線性時(shí)間計(jì)算出yiL=θi每個(gè)近似的，由于，則

設(shè)Gt對(duì)應(yīng)的鄰接矩陣為 Wn0×n1，t時(shí)刻的近似commute time嵌入的算法如下:

算法1 ApCte(approximate commute time embedding)

輸入:t時(shí)刻二部圖Gt的Wn0×n1;

輸出:的近似commute time嵌入;

步驟:

1)查找Gt中屬于的結(jié)點(diǎn)，由指示數(shù)據(jù)計(jì)算這些結(jié)點(diǎn)的k近鄰，由式(1)計(jì)算的鄰接矩陣;

4)采用 STSolve 方法［13］計(jì)算 YL=θ 的每個(gè);

5)輸出的嵌入。

數(shù)據(jù)集X0與X1的樣本映射到了一個(gè)共同的子空間。的前n0個(gè)列向量指示數(shù)據(jù)集X0，后n1個(gè)列向量指示數(shù)據(jù)集X0。設(shè)Gt中屬于的結(jié)點(diǎn)數(shù)為nt，nt＜n，算法 1 的 1)步采用kd樹構(gòu)造Gt中屬于的結(jié)點(diǎn)的k近鄰需要O(ntlnnt)時(shí)間。是稀疏圖，鄰接矩陣有s個(gè)非零元素，則2)步計(jì)算B與及L的時(shí)間為O(2s)+O(s)+O(n)。因?yàn)橄∈杈仃嘊有2s個(gè)非零元素，對(duì)角矩陣有s個(gè)非零元素，故3)步計(jì)算 θ的時(shí)間為O(2skr+s)。4)步用STSolve 方法［13］計(jì)算的時(shí)間為(skr)。則算法 1的時(shí)間復(fù)雜度僅為(ntlnnt+4s+n+3skr)。

3 基于近似commute time嵌入的異構(gòu)信息網(wǎng)絡(luò)的演化聚類

3.1 模型的組成

定義2 給定一個(gè)由M+1種類型的數(shù)據(jù)集χ=構(gòu)成的信息網(wǎng)絡(luò)G=＜V，E，W＞，如果?e=〈xi，xj〉∈E，那么xi∈X0且xj∈Xmm≠0，則G稱為星型模式的異構(gòu)信息網(wǎng)絡(luò)，X0稱為目標(biāo)類型，Xm(m≠0)稱為屬性類型。

給定一個(gè)由M+1種類型的數(shù)據(jù)集構(gòu)成的星型模式的異構(gòu)信息網(wǎng)絡(luò)，其中是Xm的對(duì)象數(shù)目，X0為目標(biāo)數(shù)據(jù)集為屬性數(shù)據(jù)集。W(0m)∈Rn0×nm表示X0與Xm之間的關(guān)系，其中，元素表示X0的樣本與Xm的樣本的關(guān)系。如果與存在關(guān)系，則與有邊存在，邊的權(quán)重為，否則無邊該信息網(wǎng)絡(luò)包含M個(gè)關(guān)系矩陣

t時(shí)刻目標(biāo)數(shù)據(jù)集X0與屬性數(shù)據(jù)集Xm構(gòu)成一個(gè)二部圖，由式(1)可計(jì)算t時(shí)刻時(shí)間平滑二部圖，設(shè)的鄰接矩陣為。由 2.2節(jié)，可計(jì)算的近似commute time嵌入Y(0m)=，其中前n0個(gè)列向量指示目標(biāo)數(shù)據(jù)集X0，表示為，稱之為指示子集，后nm個(gè)列向量指示屬性數(shù)據(jù)集Xm，表示為Y(m)。稱指示X0第i個(gè)對(duì)象的數(shù)據(jù)為指示數(shù)據(jù)，1≤i≤n0.的指示數(shù)據(jù)與X0的對(duì)象存在一一對(duì)應(yīng)的關(guān)系。M個(gè)二部圖對(duì)應(yīng)M個(gè)近似commute time嵌入，則目標(biāo)數(shù)據(jù)集X0被M個(gè)指示子集所指示，X0的每個(gè)對(duì)象被M個(gè)指示數(shù)據(jù)所指示。

設(shè)X0劃分為H個(gè)類，β(m)是矩陣的權(quán)重，其中數(shù)據(jù)屬于M個(gè)類，這M個(gè)類分別位于不同的指示子集，這M個(gè)類設(shè)置相同的類標(biāo)號(hào)。令

從相容的角度，式(2)的目標(biāo)函數(shù)F取得最小值，則目標(biāo)數(shù)據(jù)集X0聚類達(dá)到最佳。顯然，式(2)的全局最優(yōu)解是NP難問題。

3.2 快速算法的推理

3.2.1 類標(biāo)號(hào)設(shè)置

在目標(biāo)數(shù)據(jù)集X0中隨機(jī)選擇H個(gè)對(duì)象，指示這H個(gè)對(duì)象的指示數(shù)據(jù)在各自的指示子集中作為H個(gè)類的初始中心點(diǎn)，指示同一個(gè)對(duì)象的中心點(diǎn)，令其所在類的標(biāo)號(hào)一致，則其他指示同一個(gè)對(duì)象的指示數(shù)據(jù)或者都屬于第j個(gè)類，或者都不屬于第j個(gè)類，1≤j≤H。

3.2.2 加權(quán)距離總和

X0的一個(gè)對(duì)象被M個(gè)指示數(shù)據(jù)所指示，這M個(gè)指示數(shù)據(jù)到各自指示子集的類的中心點(diǎn)的距離都影響著這個(gè)對(duì)象所屬類的分配。設(shè)qi∈X0，指示qi，則權(quán)重距離總和決定了qi所屬的類。即

式中:j是qi所屬類的標(biāo)號(hào)，也是所屬類的標(biāo)號(hào)。

3.2.3F的極小值

式(2)的F也可以表示為指示同一個(gè)對(duì)象的指示

給定M個(gè)指示子集的類的初始中心點(diǎn)，首先由式(3)劃分指示子集的類，此時(shí)令F=F0;的類的中心點(diǎn)不變，然后計(jì)算的每個(gè)類的新中心點(diǎn)，新中心點(diǎn)取值所在類的所有指示數(shù)據(jù)的平均值，令

故F1≤F0。

而當(dāng)?shù)念愄鎿Q了新的中心點(diǎn)的中心點(diǎn)不變，由式(3)重新劃分類，此時(shí)令F=F2，則有F2≤F1。

算法2 EClu-pACte(evolutionary clustering algorithm based on approximate commute time embedding for heterogeneous information network)

輸 入:t時(shí) 刻的，聚類數(shù)H;

輸出:t時(shí)刻目標(biāo)數(shù)據(jù)集X0的類;

步驟:

1)form=1∶Mdo

{①由算法1計(jì)算的嵌入Y(0m);

②確定指示X0的指示子集;}

3)do

{form=1∶Mdo

{①計(jì)算式(3)劃分的H個(gè)類;

②重新確定每個(gè)類的新的中心點(diǎn)，并建立類標(biāo)號(hào);

}while式(4)收斂;

4)輸出目標(biāo)數(shù)據(jù)集X0的類。

4 實(shí)驗(yàn)

4.1 實(shí)驗(yàn)數(shù)據(jù)

從DBLP信息網(wǎng)絡(luò)選取真實(shí)數(shù)據(jù)建立實(shí)驗(yàn)數(shù)據(jù)集，取 venues、authors、papers和 terms建立書目網(wǎng)絡(luò)。選取4個(gè)學(xué)術(shù)區(qū)域建立小數(shù)據(jù)集Ssmall，這4個(gè)區(qū)域包括 database、data mining、information retrieval和machine learning。每個(gè)區(qū)域取5個(gè)有代表性的conference，共20個(gè)會(huì)議，20個(gè)會(huì)議的所有authors、papers及出現(xiàn)在論文題目中的所有terms。papers為目標(biāo)數(shù)據(jù)集，venues、authors和terms為屬性數(shù)據(jù)集，建立一個(gè)星型模式的異構(gòu)信息網(wǎng)絡(luò)。本文使用Ssmall分析kr對(duì)聚類準(zhǔn)確率的影響。

對(duì)1993－2008年上述4個(gè)學(xué)術(shù)區(qū)域的20個(gè)會(huì)議的所有authors、papers及出現(xiàn)在論文題目中的所有terms，分析目標(biāo)數(shù)據(jù)集papers的演化情況。

本文算法EClu-pACte的所有實(shí)驗(yàn)均采用文獻(xiàn)［13］的一種近乎線性時(shí)間的求解程序計(jì)算的嵌入數(shù)據(jù)集，該方法用于對(duì)角占優(yōu)矩陣。所有算法均在MATLAB環(huán)境中實(shí)現(xiàn)。

4.2 關(guān)系矩陣的確定

X0表示目標(biāo)數(shù)據(jù)集 papers，X1、X2與X3分別表示屬性數(shù)據(jù)集 authors，venues與 terms。則X0與的原始關(guān)系為的元素為

4.3 參數(shù)kr取值分析

選取小數(shù)據(jù)集Ssmall，對(duì)于給定的kr，取u=50，α=1，參數(shù)kr的變化對(duì)papers聚類質(zhì)量的影響，如圖1所示。實(shí)驗(yàn)說明當(dāng)kr＞50時(shí)準(zhǔn)確率曲線已經(jīng)趨于平滑，取kr=60很適合，則其他實(shí)驗(yàn)也取kr=60。

圖1 kr對(duì)聚類準(zhǔn)確率的影響Fig.1 The influence of kron clustering accuracy

4.4 參數(shù)α取值分析

參數(shù)α是用來平衡快照質(zhì)量與歷史質(zhì)量的。本實(shí)驗(yàn)的目標(biāo)數(shù)據(jù)集papers在不同的時(shí)間點(diǎn)(年份)是不同的，沒有重復(fù)的，papers不存在連續(xù)性，但屬性數(shù)據(jù)集存在連續(xù)性。屬性數(shù)據(jù)集的連續(xù)性影響t時(shí)刻目標(biāo)數(shù)據(jù)集papers的劃分，使得t時(shí)刻的papers與先前時(shí)間的papers存在著關(guān)聯(lián)。本次實(shí)驗(yàn)取k=7，α的取值對(duì)聚類準(zhǔn)確率的影響如圖2所示，其中圖2的聚類準(zhǔn)確率取每個(gè)時(shí)刻(1994－2008)共計(jì)15個(gè)聚類準(zhǔn)確率的平均值。說明α=0.8聚類質(zhì)量最好。

圖2 α對(duì)聚類準(zhǔn)確率的影響Fig.2 The influence of α on clustering accuracy

4.5 計(jì)算速度分析

劃分1994年的papers，比較α取不同值時(shí)不同算法的計(jì)算速度。如表1所示，當(dāng)α=1時(shí)，t時(shí)刻的二部圖沒有時(shí)間平滑，異構(gòu)信息網(wǎng)絡(luò)是稀疏的，故計(jì)算速度很快，與ENetClus算法計(jì)算速度幾乎一致。當(dāng)0＜α＜1時(shí)，計(jì)算速度略有下降，原因是需要計(jì)算中屬于的結(jié)點(diǎn)的k近鄰，以構(gòu)造稀疏的時(shí)間平滑的二部圖。但采用kd樹構(gòu)造Gt中屬于的結(jié)點(diǎn)的k近鄰僅需要O(ntlnnt)時(shí)間。實(shí)驗(yàn)的目標(biāo)數(shù)據(jù)集是papers，中的papers肯定不會(huì)出現(xiàn)在t時(shí)刻的中，故每次只需存儲(chǔ)指示屬性數(shù)據(jù)集authors、venues和terms指示子集即可，因此計(jì)算速度相對(duì)還是比較快的。

表1 計(jì)算速度比較Table 1 Comparison of computing speed s

4.6 準(zhǔn)確率分析

ENetClus算法取文獻(xiàn)［8］的最佳參數(shù)，本文算法取α=0.8，其他參數(shù)同上述實(shí)驗(yàn)。準(zhǔn)確率比較如圖3所示，本文算法的聚類質(zhì)量要明顯好于ENet-Clus算法。與ENetClus算法相比，本文算法能夠比較真實(shí)地反映在t時(shí)刻及先前時(shí)間數(shù)據(jù)對(duì)象的關(guān)系。若t時(shí)刻的數(shù)據(jù)對(duì)象與先前的數(shù)據(jù)對(duì)象不存在任何關(guān)系，本文算法能夠真實(shí)地反映出來。而ENetClus算法是從類的角度時(shí)間平滑，而不管前后兩個(gè)時(shí)刻數(shù)據(jù)是否確實(shí)存在關(guān)系，因此比較粗糙。若前后時(shí)刻數(shù)據(jù)存在關(guān)系，ENetClus能大致反映其關(guān)系，否則ENetClus反映的是錯(cuò)誤的關(guān)系。而且基于概率模型的ENetClus算法受應(yīng)用領(lǐng)域限制，通用性不強(qiáng)。

圖3 演化聚類準(zhǔn)確率比較Fig.3 Comparison of evolutionary clustering accuracy

5 結(jié)論

通過理論分析及實(shí)驗(yàn)驗(yàn)證說明:

1)本文采用時(shí)間平滑二部圖充分反映了某時(shí)刻及先前時(shí)間結(jié)點(diǎn)間的關(guān)系，聚類質(zhì)量高于以往的算法;

2)實(shí)驗(yàn)的運(yùn)行時(shí)間說明利用稀疏性，采用線性時(shí)間求解程序加快了計(jì)算速度;

3)本文算法通用性強(qiáng)，適合于任何的星型模式的異構(gòu)信息網(wǎng)絡(luò)，不受異構(gòu)信息網(wǎng)絡(luò)的應(yīng)用領(lǐng)域所限制。但本文算法參數(shù)較多，每個(gè)時(shí)間平滑二部圖的關(guān)系的權(quán)重都需要人為確定，如何自動(dòng)選取最佳的關(guān)系權(quán)重，還需進(jìn)一步研究。

［1］SUN Yizhou，HAN Jiawei.Mining heterogeneous information networks:principles and methodologies［J］.Synthesis Lectures on Data Mining and Knowledge Discovery，2012，3(2):1-159.

［2］SUN Yizhou，YU Yintao，HAN Jiawei.Rankclus:rankingbased clustering of heterogeneous information networks with star network schema［C］//Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining，ACM，2009:797-806.

［3］WANG R，Shi C ，PHILIP S Y，WU B.Integrating clustering and ranking on hybrid heterogeneous information network［M］.Berlin Advances in Knowledge Discovery and Data Mining.2013:583-594.

［4］AGGARWAL C，XIE Y，PHILIP S Y，On dynamic link inference in heterogeneous networks［C］//SDM.2013:415-426.

［5］GAO Bin，LIU Tieyan，MA Weiying.Star-structured highorder heterogeous data co-clustering based on cosistent information theory［C］//ICDM＇06.Hong Kong，China，2006:880-884.

［6］LONG Bo，ZHANG Zhongfei，WU Xiaoyun，et al.Spectral clustering for multi-type relational data［C］//Proceedings of the 23rd International Conference on Machine Learning，ACM.2006:585-592.

［7］AGGARWAL C，SUBBIAN K.Evolutionary network analysis:a survey［J］.ACM Computing Surveys(CSUR)，2014，47(1):10.

［8］GUPTA M，AGGARWAL C，HAN J，et al.Evolutionary clustering and analysis of bibliographic networks［C］//Advances in Social Networks Analysis and Mining(ASONAM).Kaohsiung，Taiwan，2011:63-70.

［9］KHOA N L D，CHAWLA S.Large scale spectral clustering using resistance distance and Spielman-Teng solvers［J］.Discovery Science，2012:7-21.

［10］QIU H，HANCOCK E R.Clustering and embedding using commute times［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2007，29(11):1873-1890.

［11］SPIELMAN D A，SRIVASTAVA N，Graph sparsi fi cation by effective resistances［J］.SIAM Journal on Computing，2011，40(6):1913-1926.

［12］SPIELMAN D A，TENG Shanghua.Nearly-linear time algorithms for preconditioning and solving symmetric，diagonally dominant linear systems［J］.SIAM Journal on Matrix Analysis and Applications，2014，35(3):835-885.

［13］KOUTIS I，MILLER G L，TOLLIVER D.Combinatorial preconditioners and multilevel solvers for problems in computer vision and image processing［J］.Computer Vision and Image Understanding，2011，115(12):1638-1646.