基于编辑距离的序列聚类算法的优化

孙启航，杨鹤标

(江苏大学 计算机科学与通信工程学院，江苏 镇江 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.