聚类分析是数据挖掘、模式识别等研究方向的重要研究内容之一, 主要是将数据集中相似的样本尽可能划分为相同的簇, 而把相异的样本尽可能划归为不同的簇.经过几十年的发展, 已经形成众多经典聚类方法^{[1, 2]}.在广受欢迎的新聚类算法方面, Frey和Dueck在2007年基于因子图的信念传播和最大化算法提出近邻传播聚类(affinity propagation, 简称AP)^{[3, 4]}, Rodriguez和Laio在2014年基于聚类中心的刻画提出快速搜索密度峰值聚类算法(clustering by fast search and find of density peaks)^[5], 为聚类算法的设计提供了一种全新思路.

K-means与K-medoids是基于划分的经典聚类算法.K-means以其低计算复杂度受到广泛欢迎, K-medoids因强鲁棒性、对噪声数据及异常点处理能力强等优势, 同样被广泛应用^{[6, 7]}.一直以来, 很多学者都致力于改进K-medoids聚类算法, 以期在效率上追赶K-means算法, 同时又维持其优点.近年来, 对K-medoids聚类效率的研究可以概括为如下4个方面.

1) 以抽样为基础, 减少聚类样本数及(或)中心点交换次数.

Kaufman提出的CLARA(clustering large applications)^[8]算法针对整个数据集进行多重抽样, 以抽样子集的最优中心点充当整个数据集的中心点进行PAM(partitioning around medoids)聚类.相对于聚类子集抽样, Ng等人提出的CLARANS(CLARA based on randomized search)^[9]算法是对中心点进行抽样, 随机地选择一个代表对象(中心点)、一个非代表对象(非中心点)计算交换代价, 当随机选择次数达到最大阈值后停止搜索, 返回当前最优结果.Barioni提出的PAM-SLIM(slim-tree applied to PAM)^[10]算法首先把数据集建成Slim-tree结构, 然后以Slim-tree中的一层作为抽样集进行PAM聚类, 得到的最优中心点作为整个数据集的最优中心点.

2) 以优化中心点为基础, 减少聚类过程中的迭代次数或者中心点交换次数.

Park等人^[11]提出的快速K-medoids算法(fast K-medoids, 简称FKMEDOIDS)以密度排序选择初始中心点, 并事先计算出所有样本的距离矩阵, 在聚类过程中, 凡涉及样本点距离计算时, 只需要调用该距离矩阵的值, 然而, 保存距离矩阵需要占据大量空间内存, 是一种典型的以空间换取时间的算法.Zadegan等人^[12]提出的排序K-medoids算法与FKMEDOIDS都事先计算相异度矩阵, 不同点在于, 基于排序更新中心点可以减少中心点交换次数, 获得聚类效率的提升.预先计算距离矩阵的思想也被二分K-medoids(bisecting K-medoids, 简称BPAM)^[13]聚类所采用, 并且在基因共表达数据上验证了该方法的实用性.Lai等人^[14]提出了方差增强K-medoids聚类, 该方法的中心点是逐个增加, 这与PAM直接给定类的个数不同; 其次, 该方法还引入了Polygon算子, 沿类中心点分割线计算方差, 以此选择合适的中心点.

3) 以并行计算为基础, 实现并行K-medoids聚类.

Jiang等人^[15]实现了基于Hapdoop分布式计算平台上的K-medoids聚类.Yue等人^[16]基于迭代MapReduce过程提出并行化K-medoids算法.Han提出了基于具备分布式、极大并行性和非确定等特点的P系统优化K-medoids算法^{[17, 18]}.

4) 以避免重复计算为基础, 减少诸如点之间距离等的重复计算.

CLATIN^[19]算法利用中心点之间的三角形不规则网络(triangular irregular network, 简称TIN)确定交换中心点后的受影响子集, 并只计算该子集的交换代价, 以此决定是否交换中心点, 有效提高了聚类效率.2002年, Chu等人^{[20, 21]}提出了部分距离(partial distance)、前一次中心点指标(previous medoid index)、三角不等式消除条件(triangular inequality elimination criteria)优化PAM的方法, 这些方法可以有效减少重复的距离计算.2007年, Chu等人^[22]在原有基础上推导出一些新的不等式, 包括涉及3个中心点之间的距离关系、新增样本点与坐标原点的距离关系并保存在内存中以供直接调用.2008年, Chu等人^[23]再次对上述优化方法进行总结, 舍弃一些优化效果不佳的不等式, 比如含3个中心点情形等, 并探讨了这些优化方法的组合情形及相应的实验结果.

纵观上述基于抽样的优化算法, 基于中心点优化及并行化加速K-medoids聚类的算法均未针对PAM算法本身进行改进, 因此, 针对PAM优化的算法可以直接应用到上述3种优化方法中.而避免重复计算的方法直接优化PAM算法, 此优化最大的优势在于:可以保持PAM算法所有的优秀性质, 并且可以直接应用在基于抽样、中心点优化、并行化等K-medoids聚类中, 使得聚类效率可以叠加.基于以上综述及文献[8, 20-23]中的加速思想, 本文将以距离不等式为基础提出SPAM(speed up PAM)优化算法, 主要包括以下3个方面.

1.将文献[8]与文献[20-23]中能够加速K-medoids聚类且具有普适性的距离不等式分别整理成引理3.1、引理3.2, 应用到SPAM算法中;

2.提出新的距离不等式定理3.1及定理3.2加速K-medoids聚类过程, 并证明其正确性;

3.提出以O(n+K²)额外内存空间开销辅助引理3.1、引理3.2、定理3.1的结合, 并应用到SPAM算法中, 使得K-medoids聚类复杂度由O(K(n-K)²)降低至O((n-K)²).

本文所涉及的4种方法——PAM, TPAM, FKEMEOIDS, SPAM及引理、定理之间具有紧密联系, 如:SPAM算法是基于TPAM算法、定理3.1及空间加速, 而TPAM算法又是基于PAM算法及引理3.2, FKMEDOIDS与TPAM均是基于PAM算法的改进.图 1可以清晰简洁地反映它们之间这种层次关系, 同时也能直观反映方法之间的差异, 如FKMEDOIDS与TPAM方法都是基于PAM方法.然而FKMEDOIDS方法利用了距离矩阵, 同时对中心点进行了优化, 这些并没有应用在TPAM方法中.

1 K-medoids算法简介

设数据集D包含n个d维欧式空间中的对象{x₁, x₂, …, x_n}, x_i∈R^d, K-medoids聚类就是把D中的对象分配到K个中心点O={O₁, O₂, …, O_K}所代表的簇C₁, C₂, …, C_K中, 使得对于1≤i, j≤K, C_i⊂D且${C_i} \cap {C_j} = \mathit{\Phi }, \cup _{i = 1}^K = D$, 并且使得簇内对象相互相似, 而与其他簇中的对象相异.围绕中心点划分PAM(partitioning around medoids)算法是K-medoids聚类的一种主流实现, PAM算法流程如下^[24]:

输入:数据集D, 簇个数K;

输出:K个簇的集合.

1)  从D中随机选择K个对象作为初始的代表对象;

2)  Repeat

3)   将每个剩余的对象分配到最近的代表对象所代表的簇;

4)   随机地选择一个非代表对象O_random;

5)   计算用O_random代替代表对象O_j的总代价S;

6)   If S < 0, thenO_random替换O_j, 形成新的K个代表对象的集合;

7) Until S不发生变化

本文以上述PAM算法为基础进行改进研究.

2 K-medoids加速理论 2.1 K-medoids距离加速理论

在PAM算法步骤3)中, 需要反复计算非代表对象与每一个中心点之间的距离, 并通过比较, 将非代表对象重新分配到合适的簇中, 每更换一个中心点, 都需要重复上述步骤并计算交换总代价.文献[8]将初始聚类(聚类过程中, 初次对非代表对象进行划分)与后续聚类(聚类过程中, 除初次划分之外的所有划分)分开处理, 在后续聚类过程中, 一般分为4种情形对所有非代表对象进行重新分派, 并计算交换总代价, 我们将上述4种情形整理成如下引理.

引理3.1(加速引理Ⅰ).设O_i, i=1, 2, …, K代表第i个簇C_i的中心点, dist(P, E)代表样本点P, E∈D之间的距离.在K-medoids聚类中, 中心点组O={O₁, O₂, …, O_K}, j=1, 2, …, K由非代表对象${O'_j}$交换代表对象O_j成为簇C_j的中心点, 即, 交换后新中心点组为$O' = \{ {O_1}, ..., {O'_j}, ..., {O_K}\} $, 则对于任意非代表对象P∈D, 有:

(1) 当P∈C_i, i≠j, 若$dist(P, {O_i}) < dist(P, {O'_j})$, 则P∈C_i; 否则, P∈C_j;

(2) 当P∈C_j, 若$dist(P, {O_j}) > dist(P, {O'_j})$, 则P∈C_j; 否则, P∈C_i.其中, i满足dist(P, Q_i)=min_l{dist(P, Q_l)|Q_l∈Q'}.

引理3.1说明:非代表对象新簇标号的更新并不需要完全重新计算该非代表对象与新中心点组中所有中心点之间的距离, 这就可以节省无效的距离计算.这种基于距离不等式的方法被应用于Chu等人^[20-23]提出的TPAM算法中(注:原文中将不等式推广在随机抽样的CLARANS算法基础上, 为了在本文中参照, 将不等式推广在PAM算法上, 将其称为TPAM算法), 并提出一个新的距离不等式.尽管原文并未将其以严格的定理形式给出, 但是对其进行了严谨的数学阐述, 现将其整理为如下加速引理:

引理3.2(加速引理Ⅱ).设O_i, i=1, 2, …, K代表第i个簇C_i的中心点, dist(P, E)代表样本点P, E之间的距离.在K-medoids聚类中, 中心点组O={O₁, …, O_j, …, O_K}, j=1, 2, …, K由非代表对象${O'_j}$交换代表对象O_j成为簇C_j的中心点, 即, 交换后新中心点组为$O' = \{ {O_1}, ..., {O'_j}, ..., {O_K}\} $.则对于任意非代表对象P∈C_i, 若$dist({O_i}, {O'_j}) \ge 2dist(P, {O_i})$, 则P∈C_i.

引理3.2考虑了前一次聚类时的中心点O_i及新交换后的中心点${O'_j}$之间的关系, ${O'_j}$与簇C_j有着紧密联系, 我们发现, 簇C_j的前一次聚类中心点O_j与当前新中心点${O'_j}$也蕴含着与引理3.2类似的关系, 我们将其称为定理3.1.

定理3.1(加速定理Ⅰ).设O_i, i=1, 2, …, K代表第i个簇C_i的中心点, dist(P, E)代表样本点P, E之间的距离.在K-medoids聚类中, 中心点组O={O₁, …, O_j, …, O_K}, j=1, 2, …, K由非代表对象${O'_j}, {\kern 1pt} {\kern 1pt} {O'_j} \notin {C_i}$交换代表对象O_j成为簇C_j的新中心点, 即, 交换后新中心点组为$O' = \{ {O_1}, ..., {O'_j}, ..., {O_K}\} $.则对于任意非代表对象P∈C_i, i≠j, 若dist(P, Q_i)≤$dist({O_i}, {O_j})/2 + |\cos \theta |dist({O_j}, {O'_j})$, 则P∈C_i.其中, $\cos \theta = \frac{{dis{t^2}({O_i}, {O_j}) + dis{t^2}({O_j}, {{O'}_j}) - dis{t^2}({O_i}, {{O'}_j})}}{{2dist({O_i}, {O_j})dist({O_j}, {{O'}_j})}}.$

证明:分两种情况进行证明.

●情形1, P恰好在平面内${O_i}{O_j}{O'_j}$, 此时令P₀=P.

令${O'_j} \notin {C_j}$, 由于C_i, C_j不是同一个簇, 则C_i, C_j必定可以由O_iO_j的中垂线分开(在R^d空间中, 任意不在同一直线上的3个点必可以在同一平面内, 以下证明均在此平面${O_i}{O_j}{O'_j}$内), 如图 2所示.

令P'是P₀在直线O_iO_j上的投影, ${O''_j}$是${O'_j}$在直线O_iO_j上的投影, 易知:

$ dist(P', {O''_j}) \le dist({P_0}, {O'_j}). $

由于P₀∈C_i, 因此P'点不可能超过O_iO_j的中垂线(图 2垂直虚线), 故:

$ dist({O_i}, {O_j})/2-\cos \theta \cdot dist({O_j}, {O'_j}) \le dist(P', {O''_j}). $

根据余弦定理可知:

$ \cos \theta = \frac{{dis{t^2}({O_i}, {O_j}) + dis{t^2}({O_j}, {{O'}_j})-dis{t^2}({O_i}, {{O'}_j})}}{{2dist({O_i}, {O_j})dist({O_j}, {{O'}_j})}}. $

根据已知:

$ dist({P_0}, {O_i}) = dist(P, {O_i}) \le dist({O_i}, {O_j})/2 + |\cos \theta |dist({O_j}, {O'_j}) $

得到:

$ dist(P, {O_i}) \le dist(P', {O''_j}) \le dist({P_0}, {O'_j}) = dist(P, {O'_j}). $

因此, P∈C_i.注意:由图中所示, 去掉绝对值符号需要添加一个负号.

●情形2, 此时P不在平面${O_i}{O_j}{O'_j}$内, 令P₀为P在该平面内的投影.则:

$ dist\left( {{P_0}, {Q_i}} \right) \le dist\left( {P, {Q_i}} \right). $

由已知条件:

$ dist(P, {O_i}) \le dist({O_i}, {O_j})/2 + |\cos \theta |dist({O_j}, {O'_j}) $

可知:

$ dist({P_0}, {O_i}) \le dist({O_i}, {O_j})/2 + |\cos \theta |dist({O_j}, {O'_j}). $

此时与情形1一致, 有:

$ dist({P_0}, {O_i}) \le dist({P_0}, {O'_j}). $

又因为P₀是P在平面${O_i}{O_j}{O'_j}$内的投影, 故:

$ dist(P, {O_i}) \le dist(P, {O'_j}). $

因此, P∈C_i.

综上两种情形, 定理得证.

定理3.1与引理3.2都是以中心点之间的距离为基础, 但是选择了不同的距离.

$dist({O_i}, {O'_j})/2$与$dist({O_i}, {O_j})/2 + |\cos \theta |dist({O_j}, {O'_j})$的值有些相近, 特殊情形时也可以相等.

●当θ < π/2时, 绝对值前需要改成负号, 这样, 定理3.1要劣于引理3.2;

●但是θ > π/2时, 虽然$dist({O_i}, {O'_j})/2$是三角形${O_i}{O_j}{O'_j}$的最长边的一半, dist(O_i, Q_j)/2是较短边的一半, 然而定理3.1中增加了$\cos \theta |dist({O_j}, {O'_j})$, 使得$dist({O_i}, {O_j})/2 + |\cos \theta |dist({O_j}, {O'_j})$可能优于$dist({O_i}, {O'_j})/2$.

为了更严格估计θ的合理范围, 假设定理3.1优于引理3.2, 此时, 不等式(1)必须成立:

$ dist({O_i}, {O_j})/2 + |\cos \theta |dist({O_j}, {O'_j}) \ge dist({O_i}, {O'_j})/2 $

(1)

注意:限定θ > π/2, 由于式(1)两边均大于0, 故可以先乘以2, 然后直接平方, 且不改变不等式符号, 可得:

$ dis{t^2}({O_i}, {O_j}) + 4|\cos \theta |dist({O_j}, {O'_j})dist({O_i}, {O_j}) + 4|\cos \theta {|^2}dis{t^2}({O_j}, {O'_j}) \ge dis{t^2}({O_i}, {O'_j}) $

(2)

将cosθ用余弦定理代入式(2), 需要注意:由于θ > π/2, 加上绝对值符号后应添加一个负号, 并进行整理得:

$ 2|\cos \theta |dist({O_i}, {O_j}) \ge (1-4|\cos \theta {|^2})dist({O_j}, {O'_j}). $

上式同除以$2|\cos \theta |, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} dist({O_j}, {O'_j})$ > 0的数, 不等式符号不变, 则:

$ dist({O_i}, {O_j})/dist({O_j}, {O'_j}) \ge (1-4|\cos \theta {|^2})/(2|\cos \theta |) $

(3)

即:只要θ满足式(3), 就能保证式(1)成立.此时, 定理3.1优于引理3.2.图 3给出了θ与(1-4|cosθ|²)/(2|cosθ|)在区间[2π/3, π]上的图像.

当θ=2π/3时, (1-4|cosθ|²)/(2|cosθ|)=8.8818e-16~0, 结合图 3可知:θ > 2π/3时, 式(3)总是成立的, 亦即定理3.1要优于引理3.2.因此, 建议θ > 2π/3时, 应用定理3.1;反之, 则用引理3.2.

引理3.1、引理3.2与定理3.1都是针对交换单个中心点的算法, 如PAM, TPAM等, 但是诸如FKMEDOIDS聚类算法, 可能同时交换多个中心点, 此时, 上述引理与定理将不再适用.为此, 我们将引理3.1推广至同时交换多中心点情形, 有如下定理.

定理3.2(加速定理Ⅱ).设O_i, i=1, 2, …, K代表第i个簇C_i的中心点, dist(P, E)代表样本点P, E之间的距离.在K-medoids聚类中, 中心点组O={O₁, O₂, …, O_K}中任意m(1≤m≤K)个代表对象被非代表对象交换后成为新中心点组$O' = \{ {O'_1}, {O'_2}, ..., {O'_K}\} $, 则对于任意非代表对象P, 有:

(1) 当P∈C_i, O_i∈O∩O', 若$dist(P, {O_i}) < {\min _j}\{ dist(P, {O'_j})|{O'_j} \in O' - O\} $, 则P∈C_i, 否则重新分配P∈C_l.其中, l满足$dist(P, {O'_l}) = {\min _j}\{ dist(P, {O'_j})|{O'_j} \in O' - O\} ;$

(2) 当P∈C_i, O_i∉O', 若$dist(P, {O_i}) > {\min _j}\{ dist(P, {O'_j})|{O'_j} \in O' - O\} $, 则重新分配P∈C_l, 其中, l满足$dist(P, {O'_l}) = {\min _j}\{ dist(P, {O'_j})|{O'_j} \in O' - O\} $; 否则, 重新分配P∈C_l, 其中, l满足$dist(P, {O'_l}) = {\min _j}\{ dist(P, {O'_j})|{O'_j} \in O'\} .$

定理3.2是引理3.1的推广, 唯一不同点在于引理3.1考虑交换一个中心点, 而在定理3.2中考虑同时交换多个中心点, 此时只需要考虑中心点组的差集O'-O就可以完成相关证明, 简略的证明过程见附录.这个定理也表明:基于距离不等式思想的加速方法可以推广到多中心点交换情形, 并不局限于PAM, TPAM, SPAM等算法.

2.2 K-medoids空间加速

虽然加速引理Ⅰ(引理3.1)已经证明, 但想要起到加速效果, 就需要先计算dist(P, Q_i), $dist(P, {O'_j})$.在加速引理Ⅱ(引理3.2)及加速定理Ⅰ(定理3.1)中, 应用它们的前提也需要事先计算$dist({O_i}, {O'_j})$, dist(P, Q_i), dist(O_i, Q_j), $dist({O_j}, {O'_j})$.对于每一个非代表对象P来说, 我们可以注意到:dist(P, Q_i), $dist(P, {O'_j})$是需要重新计算的; 而$dist({O_i}, {O'_j}), dist({O_i}, {O_j}), dist({O_j}, {O'_j})$与P本身无关, 只跟其所属簇的中心点相关.对于每一次中心点交换之后, 如果先计算中心点之间的K²/2个距离, 那么在应用引理3.2、定理3.1时, 只需要再计算一次dist(P, Q_i), 而引理3.1还需要多计算$dist(P, {O'_j})$.而dist(P, Q_i)代表上一次聚类中, 非代表对象P与其所属簇中心点之间的距离, 实际上在前一次聚类中已经计算.我们可以增加相应的空间开销来保存其距离值, 而避免其再次计算.这也是本文空间加速的基础.

以牺牲空间开销来降低时间消耗, 是一种常见的做法.Chu等人^{[22, 23]}在提出引理3.2中距离不等式时, 也考虑以空间开销换取时间消耗的方法, 只不过, 他们是在聚类之前计算所有中心点与样本点距原点的距离$\sqrt {\sum\nolimits_{i = 1}^d {o_{3i}^2} } $及$\sqrt {\sum\nolimits_{i = 1}^d {x_i^2} } $, 对于一个数据集大小为n及簇个数为K的聚类来说, 其需要额外空间开销O(n+K).但是对于本文的引理3.1、引理3.2、定理3.1来说, 更多涉及非代表对象与其所属簇中心点之间的距离及中心点组之间的K²/2个距离, 因此, 本文需要额外增加O(n+K²)的空间开销, 以使引理3.1、引理3.2及定理3.1更充分地发挥其加速效应.而非代表对象与其所属簇中心点之间的距离在聚类中是必须计算的, 只需要将其保存即可.

3 SPAM算法流程

上文中已经引入了提高K-medoids聚类效率的引理3.1、引理3.2, 同时提出了定理3.1、定理3.2, 也阐述了内存空间与引理、定理之间的结合.在SPAM算法中, 由于只涉及单个中心点的交换, 因此只使用引理3.1、引理3.2、定理3.1, 它们的组织顺序如图 4所示, 其中, 条件判断中规定:如果满足定理(引理)条件, 则为Yes, 否则为No.

结合PAM算法流程, 本文所提出的SPAM算法流程如下.

输入:数据集D, 簇个数K;

输出:K个簇的集合D'.

1)  从D中随机选择K个样本点作为K个簇代表对象;

2)  Repeat

3)   将每一个非代表对象分配到最近的代表对象所代表的簇;

4)   随机选择一个非代表对象O_random及代表对象O_j进行交换;

5)   对于每一个非代表对象P, 依次根据定理3.1、引理3.2、引理3.1判别O_random交换O_j后P所属簇并计算该交换产生的代价;

6)   计算所有样本的交换总代价S;

7)   If S < 0, then O_random替换O_j, 形成新的K个代表对象(中心点组);

8)  Until代表对象(中心点组)不发生改变

9)  输出K个簇的集合D'

4 算法复杂度分析

从PAM与SPAM的算法流程上来看, 这两者相差很小, 都是在每一次聚类迭代(中心点组更新)中有K(n-K)个O_random, O_j交换对, 并且每一个交换对都需要重新分派n-K个非代表对象的所属簇标号, TPAM算法也是如此, 因此这3种方法的聚类总代价为O(K(n-K)²).相比之下, FKMEDOIDS聚类算法每一次聚类迭代的复杂度低得多, 仅为O(nK).我们进一步分析复杂度, 因为数值的比较、加减法运算相对于数值乘除法、开平方运算要快非常多, 因此在耗时方面, 前者相对于后者几乎可以忽略, 那么整个聚类过程的时间消耗主要集中于样本点距离之间的计算.在PAM算法中, 以引理3.1为基础, 该引理中的情形(1)、情形(2)至少需要计算两次距离, 分别为dist(P, Q_i), $dist(P, {O'_j})$, 按照距离计算次数且忽略比较运算, 复杂度仍为O(K(n-K)²).在SPAM方法中, 首先引入了空间加速方法, 以O(n+K²)的额外空间保存了K²/2个中心点之间的距离及n个样本点与其所属簇中心点之间的距离, 那么与PAM方法相比, SPAM方法可以节省计算dist(P, Q_i), 在引理3.2、定理3.1中, dist(P, Q_i), $dist({O_i}, {O'_j})$, dist(O_i, Q_j), $dist({O_j}, {O'_j})$都无需重新计算, 就可以确定非代表对象P在交换O_random, O_j后所属簇; 对于被交换中心点O_j所代表簇C_j来说, 簇内样本点很难满足定理3.1, 但是有一些样本点仍可以满足引理3.2的条件, 如果连引理3.2都无法满足, 那就需要转到引理3.1中计算, 仍只需计算$dist(P, {O'_j})$.如果按平均的方式认为n-K个非代表对象平均分配到K个簇中, 即每个簇的非代表对象样本点为(n-K)/K, 按照距离计算次数且忽略比较运算, SPAM算法复杂度近似为O(K(n-K)·(n-K)/K)=O((n-K)²).而TPAM复杂度应在O((n-K)²)与O(K(n-K)²)之间.虽然FKMEDOIDS每一次迭代的复杂度仅为O(nK), 然而却需要在聚类之前计算出所有样本之间的距离矩阵及优化初始中心点时的样本密度排序, 距离矩阵计算过程的时间复杂度为O(n²), 同时还需额外的O(n²)的空间开销以保存该距离矩阵, 这比SPAM的O(n+K²)空间开销要大一个量级; 同时, 排序的最快算法也需要O(nlogn)的时间复杂度.如果设聚类过程迭代了t次, 根据上述分析, 那么整个算法的复杂度及额外空间开销见表 1.

表 1中第3栏的额外空间开销指以PAM算法的空间开销为基础, 我们将比PAM算法多出来的内存空间开销称为额外空间开销.从整个算法的时间复杂度来看, SPAM要比PAM, TPAM低, 当迭代次数t较小时, FKMEDOIDS与SPAM差不多; 反之, SPAM复杂度要高于FKMEDOIDS.然而, 从额外空间开销来看, FKMEDOIDS是最大的, SPAM与TPAM是同等量级的.

5 实验结果及讨论

本节里将在实际及人工数据集上进行算法性能测试, 并与PAM, TPAM, FKMEDOIDS算法进行比较.程序代码采用C++实现, 在64G内存的服务器上运行.

5.1 数据集

本次实验中, 我们共使用6个数据集, 其中4个数据集来自UCI^[25], 并且人工生成两个随机数据集, 数据集的具体描述见表 2.其中, 6个数据集的样本量、属性个数、类别个数均不同, 对于人工数据集RandData, 我们以二维高斯分布为基础, 分别以均值向量[10, 10], [10, 30], [10, 50], [30, 10], [30, 30], [30, 50], [50, 10], [50, 30], [50, 50]及协方差矩阵[10, 0;0, 10]随机生成900个样本, 每个类别均为100个样本.对于该二维人工数据集, 我们可以直观地从图上看出其分布情况, 如图 5所示(图例C_i代表第i类数据).对于人工数据集RandData2, 与RandData的生成方式类似, 只是该数据集生成10维10个类别的1 000个数据, 每一个类别的数据, 都是以一个维度上均值10其余为0, 协方差矩阵为5(对角线上元素为5), 每个类别均为100个样本.

5.2 性能评价指标

虽然本文致力于提高K-medoids聚类效率, 但不以牺牲聚类准确性为代价, 为此, 首先考虑聚类总代价S.

$ S = \sum\limits_{i = 1}^K {\sum\limits_{j = 1}^{{n_i}} {dist({x_j}, {O_i})} } . $

其次就是考虑聚类总时间(该时间不包括从文本中读入数据及结果输出到文本中), 然而PAM, TPAM, SPAM都是随机初始化聚类中心, 可能导致算法迭代次数(中心点更换次数)不一样而影响总时间的比较, 文中将采用10次重复实验的平均值作为讨论结果.此外, 还将给出聚类总迭代次数与平均迭代的运行时间.本文SPAM算法采用引理3.1、引理3.2、定理3.1及空间加速来提高聚类效率, 本质上就是节省样本距离的计算次数.因此, 文中也同样给出平均迭代样本距离的计算次数.

5.3 结果与讨论

对PAM, TPAM, FKMEDOIDS, SPAM这4种方法在表 2中的6个数据集进行10次重复测试, 首先, 表 3给出了10次重复实验的最小值、均值、最大值.在3种方法中, 最小值、最大值几乎相等, libras数据的最大值有微小差异, 这说明SPAM方法与PAM, TPAM相比不会降低聚类精度.

将聚类总代价画成直观的直方图并注上相应数值, 如图 6所示, 除了FKMEDOIDS算法外, PAM, TPAM, SPAM在6个数据集上的聚类总代价没有明显差异, 均能达到聚类最优结果, 略微差异可以参考柱图顶端的Cost具体值; 同时还说明, 改进的SPAM算法保持了PAM算法的鲁棒性.FKMEDOIDS算法在数据集RandData及RandData2上的聚类总代价明显高于其他3个算法, 说明FKMEDOIDS算法较易陷入局部极值, 这跟它中心点的优化方式及中心点的更新方式紧密相关, 这种中心点的优化方式降低了K-medoids算法的鲁棒性, 牺牲了聚类质量.总体上来说, FKMEDOIDS算法性能可以接受, 但PAM, TPAM, SPAM这3种算法的聚类能力要更胜一筹.

表 4给出了4种方法的聚类时间, SPAM在所有数据集上都要比PAM, TPAM快.与FKMEDOIDS算法相比, 加速效率与数据集本身特性还是相关的, 这两个方法的差异也较大.尤其是在占用内存空间方面, FKMEDOIDS要比SPAM多一个量级.

虽然聚类结果相差无几, 但聚类过程中的迭代次数却有差异, 如图 7所示, 尤其是在Libras数据集上, PAM, TPAM, SPAM算法的迭代次数将近是FKMEDOIDS算法的3倍.而在Bupa, Iris数据集方面, SPAM的迭代次数均比FKMEDOIDS算法要少.这说明4种算法的迭代次数与数据集本身的特性及初始中心点相关.因此, 在聚类效率的比较方面, 仅仅比较聚类总时间(见表 4)不足以反映4种方法之间的优劣, 故需再考虑聚类平均迭代的运行时间.

4种算法的聚类平均迭代时间如图 8所示, 具体值参看柱图上端数字.从图中可以看出:SPAM算法比PAM算法几乎快一倍(Banknote:1.311倍, Bupa:1.580倍, Iris:0.828倍, Libras:1.172倍, RandData:1.516倍, RandData 2:1.502倍); 同时, 比TPAM算法的聚类速度也快很多; 维数较多时, 效果会更明显(如Libras, RandData2两个数据集).我们还可以注意到:在Iris, Libras, Randdata, RandData2数据集上, FKMEDOIDS算法比PAM, TPAM, SPAM都要快; 然而在其他两个数据集上, SPAM要快于FKMEDOIDS算法, 因此, SPAM与FKMEDOIDS算法的性能跟数据集特性相关, 两种算法的性能都有可能超越对方.但是, FKMEDOIDS算法额外开销了O(n²)的内存空间, 而SPAM仅仅额外开销了O(n+K²)内存空间, 比FKMEDOIDS少一个量级, 当SPAM聚类速度低于FKMEDOIDS时, 数值上也相差较小.综合起来, SPAM利用了更少资源提高聚类效率.

聚类平均迭代时间仅仅是从结果方面进行的评价, 本文提出采用引理3.1、引理3.2、定理3.1及空间加速方法加速K-medoids算法, 本质上是节省大量的重复样本间距离计算, 因此可以统计每次迭代过程中样本距离的计算次数, 来反映本文方法的提高效率的实质, 如图 9所示.

首先分析PAM, TPAM, SPAM这3种算法, 从图 9中可以明显反映出:SPAM算法计算样本距离次数要比PAM, TPAM算法少非常多, 这就可以解释为什么SPAM算法可以提高聚类效率.图中显示, FKMEDOIDS算法的样本距离计算次数几乎看不见.这是因为该算法在聚类之前进行了样本点的距离矩阵计算, 在后面的聚类过程中, 只需要调用距离矩阵的样本点距离值即可, 无需重新计算, 使得该算法需要额外的O(n²)内存空间, 这也合理解释了在算法时间复杂度上需要加上O(n²), 变成O(n²+tnK).虽然图 9中反映FKMEDOIDS计算距离次数最少, 但是表 4中的聚类总时间及图 8中的聚类平均迭代时间都反映出, SPAM算法在数据集Banknote, Bupa上的效率要高于FKMEDOIDS.

上述实验数据集均比较小(样本最多的Banknote仅1 372), 只能说明在较小样本集情形下, SPAM算法性能比PAM, TPAM优秀, 而与FKMEDOIDS性能难分伯仲; 对于不同的数据集, 性能上均有可能超越对方.为了更全面地验证SPAM算法性能, 我们随机生成一个10 000个样本的数据集, 分别聚成5, 10, 15, …, 100个簇, 比较聚类平均迭代时间, 见图 10.图中反映出PAM, TPAM, SPAM这3个算法随着聚类数目的增加, 耗时也随之增加, 但是3个算法的耗时增率却依次递减, 即:SPAM算法随着聚类数目的增多, 耗时的增加量要明显小于PAM, TPAM算法.对于FKMEDOIDS聚类算法来说, 总体趋势是随着聚类数目的增加, 聚类时间降低, 这跟其聚类之前需要计算距离矩阵与中心点优化的样本密度排序消耗时间相吻合.而与SPAM相比较而言, 小于50个簇时, SPAM效率要高于FKMEDOIDS; 大于50个簇时, 性能反转.总体而言, 对于大数据集来说, SPAM性能要优于PAM, TPAM; 当簇个数非常大时, 也无法保证SPAM性能超越FKMEDOIDS, 此时, 有待于进一步提高效率.

总体上来说, 在聚类的质量方面, PAM, TPAM, SPAM极其相近, 且略优于FKMEDOIDS算法; 在聚类效率方面, SPAM算法完全优于PAM, TPAM算法, 而与FKMEDOIDS算法互有优劣.这种结论上的差异, 可以从SPAM与FKMEDOIDS原理上得到印证.FKMEDOIDS算法事先计算并保存所有样本的距离矩阵, 而SPAM算法无需执行这一步骤; 在初始中心点选择方面, FKMEDOIDS算法需要计算样本密度值并排序, 选择前K个样本点作为初始中心点, 而SPAM随机选择K个样本点; 还有一个显著差异在中心点的更新方式上, 这也导致了定理3.1无法优化FKMEDOIDS算法, 而应该采用定理3.2的结论

6 结论

本文针对K-medoids算法效率较低的问题, 在PAM, TPAM算法的基础上, 利用距离不等式, 并结合O(n+K²)空间开销, 提出了SPAM算法, 进一步提升了K-medoids聚类的效率; 在时间复杂度方面, 在一次中心点更换过程中, 即t=1时, 从PAM, TPAM算法的O(K(n-K)²)降至SPAM算法的O((n-K)²); 同时, 本文的改进保持了K-medoids算法的一切优点, 并可以推广到现有的很多优化算法中, 如基于抽样、中心点优化、并行化等.在UCI实际数据集及人工数据集的实验均表明, SPAM算法的效率高于PAM, TPAM算法.与FKMEDOIDS算法的实验对比还表明:在比其少消耗内存空间的情况下, SPAM算法效率并不一定比FKMEDOIDS差; 与此同时, SPAM算法在诸如Libras属性个数特别多的数据集上效率提升较少, 仍需进一步加以改进.

附录

定理3.2的简要证明.

需要分4种情形.

1) P∈C_i,O_i∈O∩O',若$dist(P,{O_i}) < {\min _j}\{ dist(P,{O'_j})|{O'_j} \in O' - O\} $,此时说明在所有的新交换中心点中,点P与新中心点的距离均大于dist(P,Q_i),因此,P∈C_i;

2) P∈C_i,O_i∈O∩O',若$dist(P,{O_i}) < {\min _j}\{ dist(P,{O'_j})|{O'_j} \in O' - O\} $不成立,此时说明在所有的新交换中心点中, 存在新中心点,使得点P与该新中心点的距离小于dist(P,Q_i),因此,P需要重新分配到该新中心点(最小距离)所代表的簇中,即,满足$dist(P,{O'_l}) = {\min _j}\{ dist(P,{O'_j})|{O'_j} \in O' - O\} ;$

3) 当P∈C_i,O_i∉O',说明点P所在的簇i的中心点已经更换,若$dist(P,{O_i}) > {\min _j}\{ dist(P,{O'_j})|{O'_j} \in O' - O\} $,则说明存在新中心点,使得点P到该新中心点的距离要小于原始簇的距离dist(P,Q_i),此时需要P重新分配到该新中心点(最小距离)所代表的簇中,即,满足$dist(P,{O'_l}) = {\min _j}\{ dist(P,{O'_j})|{O'_j} \in O' - O\} ;$

4) 当P∈C_i,O_i∉O',说明点P所在的簇i的中心点已经更换,若$dist(P,{O_i}) > {\min _j}\{ dist(P,{O'_j})|{O'_j} \in O' - O\} $不成立,则说明所有新更换的中心点到P的距离均大于原始簇的距离dist(P,Q_i),此时需要比较所有中心点,以确定P所属的簇,即,满足$dist(P,{O'_l}) = {\min _j}\{ dist(P,{O'_j})|{O'_j} \in O'\} .$