基于編輯距離的序列聚類算法的優化

孫啟航，楊鶴標

(江蘇大學 計算機科學與通信工程學院，江蘇 鎮江 212013)

1 概 述

序列聚類的基礎問題就是序列的相似性度量[1]。直觀看來，當序列中對應位置的對象的值相似時，才被認為是相似的。在現有的研究中，根據度量范圍的不同，序列相似性度量算法可分為如下兩類：

(1)局部序列相似性度量。在生物信息學研究中，DNA序列或蛋白質序列中往往存在著最能體現序列整體特征的關鍵片段，這些關鍵片段在刻畫序列特點時占有相當高的比重，因此可以用來度量序列的相似性[2]。基于局部序列的相似性度量算法，核心就是要提取出能表征不同序列的局部序列，從而在局部序列的基礎上度量全序列的整體相似性。

(2)全序列相似性度量。序列的特征通過整個序列來體現，這時度量序列的相似性必須從序列的全局考慮。例如有5個患者P1、P2、P3、P4、P5在不同時間點的臨床行為序列如圖1所示。

PatientP1’sclinicalsequence:14 21 30 34 26 67 90 45 70 29PatientP2’sclinicalsequence:33 21 62 92 17 76 19 43 70 29PatientP3’sclinicalsequence:33 95 62 34 17 67 19 45 57 56PatientP4’sclinicalsequence:72 72 54 54 46 68 53 70 57 56PatientP5’sclinicalsequence:11 38 59 80 22 22 16 65 57 56

圖1 患者的臨床行為序列

盡管患者P3與患者P4、P5的臨床行為序列有一致的頻繁子模式57、56，但是患者P3與P2、P3與P1、P1與P2在對應時刻的相同醫療行為數量分別為4,3,3，因此P1、P2、P3整體上的臨床行為更相似。編輯距離的相似性度量算法能夠準確反映序列全局相似度[3]，但其時間復雜度相對較高。

2 相關技術

K均值是基于劃分的聚類方法[4-6]，該算法首先指定K值，并在數據集中選擇K個點作為簇的初始質心。然后將每條序列都分配到距離它最近的質心，當完成一遍劃分后，更新每個簇的質心。重復分配序列的過程，直到質心不再變化，聚類結果才是穩定的，最終所有數據都被劃分到K個簇中。

相比于K均值聚類，二分K均值[7]避免了兩個序列之間距離的直接計算，時間復雜度是線性增長的。序列的二分K均值聚類算法可作如下描述[8]：把所有序列初始化為一個簇加入簇表；從簇表中選取一個總體相似度水平最低的簇C；利用K-means方法將簇C二分為C1和C2；將C1和C2加入簇表中；重復以上步驟，直到產生K個簇時停止。

編輯距離[9]是計算兩個序列之間相似度的算法。求解的常用思路是動態規劃法[10]，即對于兩個序列A={a1,a2,…,am}和B={b1,b2,…,bn},用dis(i,j)表示序列A的前i(0≤i≤m)個項通過編輯轉化為序列B的前j(0≤j≤n)個項所需的最小代價。序列A和B的相似度計算方法具有子結構和子問題重疊性質，可以通過不斷地遞歸求解dis(i,j)最優解得到。具體的相似度計算步驟如下：

(1)參數初始化。

distance(A，B)=dis(m,n)；dis(0,0)=0；

dis(0,j)=j；dis(i,0)=i

(2)遞歸計算序列相似度。

1≤i≤m,1≤j≤n

(1)

其中，當ai=bi時，k(i,j)=0,否則k(i,j)=1。

通過在動態規劃矩陣中對式(1)的計算，矩陣最右下角的dis(m,n)可以通過不斷得到子問題的最優解來迭代解決，最終解得兩個序列的編輯距離。

3 帶剪枝策略的序列聚類算法

3.1 基于全序列比對的編輯距離上下界

定義1：序列是由項表中的項目組成的有序項目集合，記為S={a1,a2,…,an},其中ak∈E(1≤k≤n)，ak稱為序列中的項。T為臨床行為項目表中項目的個數，T=|E|。序列S的長度為n，n=|S|。

(2)

標簽距離和編輯距離的關系如定理1所示，這給后文聚類算法的剪枝策略提供了理論基礎。

定理1：對于給定的兩個臨床行為序列S1和S2，標簽距離TD(S1,S2)與兩條序列的長度之和|S1|+|S2|分別是編輯距離ED(S1,S2)的下限和上限，即TD(S1,S2)≤ED(S1,S2)≤|S1|+|S2|。

TD(S1,S2)≤ED(S1,S2)的證明見文獻[11]。下面證明ED(S1,S2)≤|S1|+|S2|。

在S1和S2的求解矩陣中，通過動態規劃將S1編輯為S2時，從矩陣最右下角的(m,n)位置開始回溯，直到(0,0)時終止。每一次回溯或者是縱坐標減1，或者是橫坐標減1，或者是橫坐標和縱坐標均減1，最多經過m+n步可以走到(0,0)坐標，所以ED(S1,S2)≤|S1|+|S2|。故可推出定理1,TD(S1,S2)≤ED(S1,S2)≤|S1|+|S2|。

由定理1可以推出下列結論：

推論1：給定三個臨床行為序列S1、S2、S3，若TD(S1,S2)≥|S1|+|S3|,則有ED(S1,S2)≥ED(S1,S3)。由定理1得：ED(S1,S2)≥TD(S1,S2)≥|S1|+|S3|≥ED(S1,S3)，故可得推論1。

定理2：給定三條臨床序列S1、S2、S3，S1與S2的編輯距離ED(S1,S2)已知，若ED(S1,S2)≥2×(|S1|+|S3|)，則ED(S1,S3)≤ED(S2,S3)。

證明：由于|S1|+|S3|≥ED(S1,S3)，所以ED(S1,S2)≥2(|S1|+|S3|)≥2ED(S1,S3)，通過移項得到ED(S1,S2)-ED(S1,S3)≥ED(S1,S3)。由編輯距離的定義可知其滿足三角不等式兩邊之差小于第三邊[12]，因而有ED(S2,S3)≥ED(S1,S2)-ED(S1,S3)，所以ED(S2,S3)≥ED(S1,S3)。

計算序列編輯距離的時間復雜度為O(m×n)[13]，計算標簽距離的時間復雜度同為O(m×n)。對于待劃分的序列和兩個簇的質心，如果符合推論1和定理2中的條件，則可以僅通過計算標簽距離來劃分數據點,從而有效降低序列之間編輯距離的計算量。

3.2 基于前綴子序列的相似度比較

若推論1和定理2中的兩個前提條件都不能滿足，也有可能不需要進行兩條序列的全序列比較，轉而通過計算序列的部分編輯距離來衡量它們的臨近度。對于給定的兩個序列A={a1,a2,…,am}和B={b1,b2,…,bn}，假設其動態規劃矩陣為T(T為一個m+1行、n+1列的矩陣)，則可以推導出下面的性質：

性質1：|dis(i,j)-dis(i-Δi,j-Δj)|≤1。

其中，dis(i,j)表示動態規劃矩陣中第i行第j列的值，0≤i≤m,0≤j≤n，Δi,Δj∈(0,1),并且i-Δi≥0,j-Δj≥0。

證明：計算兩條序列編輯距離的動態規劃矩陣是由m+n-1個初始狀態dis(i,j)=0、dis(0,j)=j,dis(i,0)=i(0≤i≤m,0≤j≤n)和式(1)生成的，若ai=bj，則k(i,j)=0，否則k(i,j)=1。由歸納法證明dis(i,j)-dis(i-Δi,j-Δj)不會超過1，性質1得證。

由性質1表明，在反映編輯距離的動態規劃矩陣中，任何一個坐標對應的值減去緊鄰其左邊、上邊或者左上角的值，只會有0和1兩種結果，不會超過1。圖2所示的動態規劃矩陣揭示了性質1，且由性質1可以證明定理3。

圖2 序列animals與bicycle的動態規劃矩陣

定理3：對于給定的兩條序列A={a1,a2,…,am}和B={b1,b2,…,bn}(m>n)，其最小編輯代價的動態規劃矩陣為T，則當i≤j≤n時，存在dis(i,i)≤dis(j,j)。

證明定理3只要證明dis(i,i)≤dis(i+1,i+1)即可。根據動態規劃矩陣T的生成過程，可以對以下三種情況進行分析：

(1)若dis(i+1,i+1)=dis(i,i)+k(i+1,i+1)，因為k(i+1,i+1)的取值為0或1，所以dis(i+1,i+1)≥dis(i,i)。

(2)若dis(i+1,i+1)=dis(i+1,i)+1，根據性質1可得|dis(i+1,i)-dis(i,i)|≤1。分三種情況討論：

①dis(i+1,i)-dis(i,i)=1，此時dis(i+1,i+1)=dis(i,i)+2，dis(i,i)

②dis(i+1,i)-dis(i,i)=0，此時dis(i+1,i+1)=dis(i,i)+1，dis(i,i)

③dis(i+1,i)-dis(i,i)=-1，此時dis(i+1,i+1)=dis(i,i)。

(3)若dis(i+1,i+1)=dis(i+1,i)+1，參考情況(2)中的證明過程，仍可證得dis(i+1,i+1)≥dis(i,i)。

綜上所述，有dis(i,i)≤dis(i+1,i+1)。故定理3得證。

定理3表明，在對兩個序列進行全序列編輯距離的計算時，其長度相等的前綴子序列具有隨著子序列中項的增加而編輯距離非遞減的特性。例如圖2中兩種不同灰度的陰影部分所顯示的，有ED(an,bi)≤ED(ani,bic)。

在臨床行為異常檢測中，需要將待檢序列S和兩個簇的質心C1、C2進行編輯距離的比較，若已知S和C1的編輯距離ED(S,C1)，又知道S和C2某個等長前綴子序列的編輯距離比ED(S,C1)大，那么不用計算S和C2的全序列編輯距離就可以將S劃入C1所代表的簇中。對于長度不同的臨床序列，通過定理4，當在聚類過程中的編輯距離比較時，其中一個編輯距離只需求出等長前綴子序列的部分即可。

定理4：設有三條序列S、P、Q，其中S={s1,s2,…,sm}，P={p1,p2,…,pn}(m>n)，S作為待分配序列需要比較它與P、Q的編輯距離的大小。若S、Q的編輯距離ED(S,Q)已知,且存在一個k≤n使得不等式ED(s1,s2,…,sk,p1,p2,…,pk)-(m-n)≥ED(S,Q)成立，則不等式ED(S,P)≥ED(S,Q)成立。

證明：設S和P的最優動態規劃矩陣為T，將性質1中的絕對值符號去掉可以得到-1≤dis(i,j)-dis(i-Δi,j-Δj)≤1，故可以推得dis(m,n)+1≥dis(m-1,n)，dis(m-1,n)+1≥dis(m-2,n)，…，dis(n+1,n)+1≥dis(n,n)，將這m-n個不等式合并可得dis(m,n)≥dis(n,n)-(m-n)。根據定理3，若存在一個k≤n，則有dis(k,k)≤dis(n,n)。因此，若k≤n且ED(s1,s2,…,sk,p1,p2,…,pk)-(m-n)≥ED(S,Q)成立，可得ED(S,P)=dis(m,n)≥dis(n,n)-(m-n)≥dis(k,k)-(m-n)≥ED(S,Q)。

定理4的推導表明，在比較待檢序列到簇的距離時，有時并不需要計算S和Q的全序列編輯距離就可以與ED(S,Q)進行比較，這大大降低了二分K均值算法在計算臨床序列在計算編輯距離上的開銷。

3.3 序列聚類算法的實現過程

對于數據集中的t條序列，其平均長度為L，在以編輯距離為臨近度度量的前提下，如果采用兩兩比較的方式來計算序列的編輯距離，算法的時間復雜度為O(t2L2)，這嚴重影響了算法的聚類效率。為降低算法的復雜度，減少聚類時間，從以下幾個方面來考慮：

(1)在聚類算法的選取上，選取二分K均值的方法。相比于K均值聚類，二分K均值避免了兩個序列之間距離的直接計算，時間復雜度與數據集大小不是呈指數增長，而是線性增長的。

(2)在處理簇間數據點的移動時，采取計算序列和簇的標簽距離的方式避免不必要的計算。假設在應用二分K均值算法時，當前的兩個質心是C1和C2，S表示序列數據集中的任意一條序列，且兩個質心之間的編輯距離ED(C1，C2)已知。在尋找離S點最近的簇時，需要計算ED(S,C1)和ED(S,C2)并比較這兩者的大小。由推論1和定理2可知，當Max{TD(S,C1)，ED(C1,C2)/2}≥|C2|+|S|時，ED(S,C2)≤ED(S,C1)，序列S放入C2所代表的簇；當Max{TD(S,C2)，ED(C1,C2)/2}≥|C1|+|S|時，有ED(S,C1)≤ED(S,C2)，序列S放入C1所代表的簇。

若上述條件都不滿足，才需要計算ED(S,C1)和ED(S,C2)的值。但是，也有可能只需計算S與其中一個質心的等長前綴子序列的編輯距離而不必計算全部就可以比較它們的大小。步驟(3)的剪枝策略描述了S到兩個質心的編輯距離比較方法。

(3)對于臨床行為序列數據集中的每一個序列S，創建一個維數為2T的概貌向量Vs={n1,n2,…,nT,d1,d2,…,dT}，其中ni(1≤i≤T)是S的標簽向量中的值，表示臨床行為項表中第i個項在S中出現的次數，di表示序列S中的項ei到序列第一個項的距離之和。概貌向量在序列標簽的基礎上增加了反映序列中項的位置信息的維度，如果兩條臨床序列差異很大，它們所對應的概貌向量的曼哈頓距離也會很大[14]。

分析定理4可以發現，若想通過該定理減少編輯距離的計算復雜度，需要事先知道待檢序列與哪一個簇的質心的編輯距離較小，通過計算三個序列概貌向量的曼哈頓距離來幫助判斷。序列概貌向量的曼哈頓距離在很大概率上能幫助確定編輯距離的大小，從而有效決定待檢臨床序列與簇的質心編輯距離的計算順序，最終達到簡化計算序列臨近度時間開銷的目的。

3.4 帶剪枝策略的二分K均值序列聚類算法

綜合上一節的三點陳述，帶剪枝策略的二分K均值臨床序列聚類算法PSClu(clustering algorithm based on similarity of prefix sequence)描述如下：

算法1：帶剪枝策略的臨床序列聚類算法PSClu。

輸入：含有t條臨床行為數據序列的數據集SequenceSet={S1,S2,…,St},參數k

輸出：序列的k個集合

/*起始狀態下，將原始序列數據集SequenceSet初始化為一個簇并放入簇表ClusterList中，此時簇的個數ClusterCount=1*/

ClusterList←SequenceSet；ClusterCount=1；

for(i=1;i≤t;i++)

將序列數據集SequenceSet掃描一趟，生成Si的標簽T(Si)以及概貌向量VSi;

end for;

while(ClusterCount

從簇表ClusterList中選一個內部相似程度最差的簇C；

分別將簇C中的最長序列和最短序列作為兩個簇的質心CO1和CO2；

計算兩個質心的編輯距離ED(CO1,CO2)以及概貌向量VCO1和VCO2；

while(簇C中所有序列不會在CO1和CO2兩個質心代表的簇之間發生移動)

for(屬于簇C的每一條序列S’)

if(Max{TD(S’,CO1),ED(CO1,CO2)/2}≥|CO2|+|S’|)

將序列S’放到以CO2為質心的簇C2中；

else if(Max{TD(S’,CO2),ED(CO1,CO2)/2}≥|CO1|+|S’|)

將序列S’放到以CO1為質心的簇C1中；

else if(L1(VS',VCO1)≤L1(VS',VCO2))

計算ED(S’,CO1);

if(S’和CO2的等長前綴子序列滿足定理4前提條件中的不等式)

終止S’與CO2編輯距離的計算，并將序列S’放到以CO1為質心的簇C1中；

else

繼續利用動態規劃矩陣算完ED(S’,CO2),然后確定將S’放到哪個簇；

else

計算ED(S’,CO2);

if(S’和CO1的等長前綴子序列滿足定理4前提條件中的不等式)

終止S’與CO1編輯距離的計算，并將序列S’放到以CO2為質心的簇C2中；

else

繼續利用動態規劃矩陣算完ED(S’,CO1),然后確定將S’放到哪個簇；

end for;

更新兩個簇的質心CO1和CO2；

end while;

ClusterCount=ClusterCount+1；

將C1和C2添加到簇表ClusterList中；

end while;

4 實驗結果分析與驗證

為了驗證PSClu算法在對序列進行聚類分析時的性能，使用Java語言編碼，并在人工合成的序列數據集上與其他算法進行對比實驗。所有實驗均在主頻為2.8 GHz，操作系統為Windows7，可用內存為3.24 G的PC機上實現。

4.1 實驗數據集的生成

類似文獻[15]中實驗所使用的方法，該課題合成數據集的方法如下：給定項表E(表中項的個數設定為10，用a到j十個英文字母表示)，先隨機生成k條長度不同的根序列，并以這k條根序列代表k個簇，然后在每條根序列的基礎上，通過隨機變換序列長度和其中的項來合成其他全部序列。合成序列的特征通過如下的標識符進行描述：K表示合成的數據集中簇的個數，C表示數據集中序列的個數，L表示最短根序列的長度，△表示其他根序列長度在最短根序列長度的基礎上的遞增度，VL表示根序列所代表的簇中的其他序列長度相對于根序列變化的百分率約束，VP表示根序列所代表的簇中的各序列變化的元素相對于根序列變化百分率的約束。例如，對于數據集K6C6000 L50△50VL5VP10的解釋就是，該數據集中有6個簇，一共6 000條序列分布在這6個簇中，因為最短根序列的長度為50，所以其他根序列的長度在其基礎上以50位單位遞增，所以6條根序列的長度分別為50,100,150,200,250,300。每個根序列S所代表的簇中，其他序列都是通過變化S長度(隨機在任意位置上添加或刪除項)的至多5%以及改變項(隨機將序列中的任意項替換成其他項)的至多10%而合成的。

4.2 聚類效率的比較

對這兩種算法進行測試的四組合成數據集特征如表1所示。對于這四組數據集，控制住其他特征不變，將參數C由1 000增長到10 000。

表1 合成的四組序列數據集

PSClu與EDClu算法執行的時間對比見圖3。

圖3 EDCluster與PSClu聚類時間對比

從圖中可以看出，帶剪枝策略、使用序列局部相似性來代替全局相似性的PSClu算法的運行時間比不帶剪枝策略的EDClu要少。

5 結束語

針對序列聚類算法中序列相似性度量算法的不足，提出了PSClu序列聚類算法。通過編輯距離的上下界以及前綴子序列的相似性計算減少序列之間兩兩編輯距離的計算。實驗表明，PSClu算法在時間效率上明顯優于EDClu算法。由于在最差的情況下，PSClu的時間復雜度仍然是與序列均長L呈平方項時間復雜度的，今后的研究目標將采用隨機近似的方法進一步降低PSClu算法的時間復雜度。

[1] 戴東波,湯春蕾,熊 赟.基于整體和局部相似性的序列聚類算法[J].軟件學報,2010,21(4):702-717.

[2] 唐東明,朱清新,楊 凡,等.一種有效的蛋白質序列聚類分析方法[J].軟件學報,2011,22(8):1827-1837.

[3] WAGNER R A,FISCHER M J.The string-to-string correction problem[J].Journal of the ACM,1974,21(1):168-173.

[4] 袁 方,周志勇,宋 鑫.初始聚類中心優化的k-means算法[J].計算機工程,2007,33(3):65-66.

[5] 王 勇,唐 靖,饒勤菲,等.高效率的K-means最佳聚類數確定算法[J].計算機應用,2014,34(5):1331-1335.

[6] CELEBI M E,KINGRAVI H A,VELA P A.A comparative study of efficient initialization methods for the k-means clustering algorithm[J].Expert Systems with Applications,2013,40(1):200-210.

[7] HAMERLY G,ELKAN C.Alternatives to the k-means algorithm that find better clusterings[C]//Proceedings of the eleventh international conference on information and knowledge management.New York,NY,USA:ACM,2002:600-607.

[8] 戴 紅,常子冠,于 寧.數據挖掘導論[M].北京:清華大學出版社,2015.

[9] ILIOPOULOS C S,RAHMAN M S.New efficient algori-thms for the LCS and constrained LCS problems[J].Information Processing Letters,2008,106(1):13-18.

[10] SAKOE H,CHIBA S.Dynamic programming algorithm optimization for spoken word recognition[M].[s.l.]:Morgan Kaufmann Publishers Inc.,1990.

[11] BANG K S,LU H,SIAU K.An efficient index structure for spatial databases[J].Journal of Database Management,1996,7(3):3-16.

[12] RISTAD E S,YIANILOS P N.Learning string-edit distance[J].IEEE Transactions on Pattern Analysis & Machine Intelligence,1998,20(5):522-532.

[13] 姜 華,韓安琪,王美佳,等.基于改進編輯距離的字符串相似度求解算法[J].計算機工程,2014,40(1):222-227.

[14] 劉瑩霞.鏈碼技術和聚類分析在基因序列中的應用[D].廣州:華南理工大學,2012.

[15] AGRAWAL R, SRIKANT R. Mining sequential patterns[C]//Eleventh international conference on data engineering.Washington,DC,USA:IEEE Computer Society,1995:3-14.