基于GSA算法改进的K均值聚类

娄 奥，姚敏立，袁 丁

(火箭军工程大学 作战保障学院，陕西 西安 710025)

0 引 言

传统K均值聚类存在对初始聚心敏感、全局搜索能力弱和人工设定聚类数等缺点。为改善K均值聚类的性能，欧阳浩等[1]引入小生境和禁忌算法的思想，田诗宵等[2]则提出基于密度峰值优化的算法改进方案，这些方法显著提升了正确率，但未考虑算法稳定性的问题。Goel等[3]和Pizzuti等[4]引入智能算法思想，把进化算子应用到迭代寻优的过程中，提高了算法鲁棒性，但仍存在参数较多、运算繁琐等问题。

万有引力搜索算法(gravitational search algorithms，GSA)是伊朗科学家Esmat Rashed等提出的一种启发式群智能算法。有研究结果表明，当优化基准测试函数时，GSA算法的寻优精度和速度等都明显优于遗传算法、模拟退火算法等其它智能算法[5]。

本文针对传统K均值聚类的缺点，提出一种基于GSA算法的改进K均值聚类。该算法以均方误差作为GSA的适应度函数，全局搜索聚类质量最优的初始聚类中心；引入种群成熟度因子避免GSA陷入局部最优；设置戴维森堡丁指数为聚类质量评价指标，确定最佳的聚类数。本文将改进后K均值聚类与工程实践中常用的最大最小距离聚类、传统K均值聚类、K-means++聚类等算法作实验对比。结果验证，该算法具有更好的聚类效果和稳定性。

1 GSA算法

牛顿第二定律指出，自然界任何两个物体之间都是相互吸引的，它们之间存在的力称作引力。引力的大小与两物体质量的乘积成正比，与距离的平方成反比。

GSA算法作为一种通用型优化算法，吸收牛顿第二定律的特点，在求解优化问题时，不仅搜索粒子位置，还考虑粒子质量。粒子质量用来评价粒子位置的优劣，粒子位置越好，质量越大。粒子之间相互吸引并向质量较大的粒子位置方向移动，通过迭代运算，质量最大的粒子位置成为最优解。

(1)

(2)

其中，fitnessi(t) 表示算法第t次迭代时粒子i的适应度值。对于最小值优化问题，最优适应度best_fit(t) 和最差适应度worst_fit(t) 分别为

(3)

(4)

用于最大值优化问题时则与之相反。

第t次迭代时，定义粒子i和粒子j在第d维的相互吸引力大小为

(5)

其中，Mi(t) 和Mj(t) 分别表示粒子i和粒子j的惯性质量；ε为常量；G(t) 表示第t次迭代时的引力系数；Rij(t) 表示粒子i和粒子j之间的欧式距离；计算公式分别为

(6)

(7)

其中，G0是引力系数初值；a是引力系数的衰减因子，一般为定值；T表示最大迭代次数。

第t次迭代时，粒子i在第d维所受的总作用力为

(8)

其中，rand1表示一个[0,1]区间内的随机数；Kbest初始值为N， 随迭代次数增加线性减小至1，定义为

(9)

其中，final_per表示对其它粒子产生作用力的粒子百分比。

第t次迭代时，粒子i在第d维上的加速度为

(10)

每次迭代进化时，粒子i的速度和位置由以下公式进行更新

(11)

(12)

其中，rand2表示一个[0,1]区间内的随机数。

2 K均值聚类算法

K均值聚类是无监督的硬聚类算法，把样本点到中心点的某种距离作为优化目标，目的是使簇内各个对象的相似度尽可能高，且各簇之间相似度尽可能小。

K均值聚类通常把欧式距离作为相似度测度，取误差的平方和作为聚类准则函数。在算法运行时，K均值聚类通常有两种迭代终止条件供选择：①聚类准则函数值尽可能小；②聚类中心点位置未更新。

K均值聚类的具体流程如下：

input: 聚类中心个数k；N个样本点组成的数据集

步骤1 随机选择k个样本点作为初始聚类中心。

步骤2 根据最小距离准则，把剩余样本点赋给k个簇。

步骤3 计算每个簇内所有点的平均值，得到新的聚类中心集合。

步骤4 循环步骤2和步骤3，直至聚类中心不变或误差平方和减小至阈值。

output:k个两两之间相似度很低的簇

3 基于GSA算法改进的K均值聚类

经典K均值聚类随机选取初始聚心和人工选择聚类数，这严重影响算法的搜索精度和稳定性。考虑GSA算法的全局寻优能力强，而K均值聚类又有局部精确解的搜索能力，将两者有机结合，优势互补，是本文算法改进的研究方向。

3.1 聚类质量评价指标

本文引入戴维森堡丁指数(Davies-Bouldin index，DBI)作为聚类质量评价指标获取最佳聚类数。有研究结果表明，DBI对聚类结果中数据成员变动的敏感性很高，可作为评价聚类数优劣的依据[6]。

DBI指数的实质是度量每个簇最大相似度的均值，其计算公式如下所示

(13)

式中：Si指第i个簇内数据到中心点的平均距离，代表各粒子的位置分散程度，如式(14)所示。N是簇的总数

(14)

其中，Ai代表第i个簇的中心点，Xij代表第i个簇内第j个样本点，Ti表示第i个簇内样本点总数。

式(13)中Rij定义为第i个和第j个簇的中心点之间的欧式距离

(15)

DBI指数越小，说明不同簇的相似度越小，即聚类质量越好。

3.2 种群成熟度因子

GSA算法虽然在迭代初期有较高的搜索速度，但临近收敛时极易陷入局部最优，效率也会大大降低。由文献[7]可知，可以依据种群多样性的变化寻找算法收敛的临界点。常见的判断种群多样性大小的方式有平均粒距和粒子之间适应度的差异[8]。但这两种方法都未考虑粒子未来的运动趋势，所以并不完善。

根据模糊理论，粒子之间越相似则种群多样性越低[9]。粒子之间的相似程度可以用它目前位置和历史最优位置的相似程度来表示。当种群内粒子的相似程度大于阈值时，说明该区域粒子分布较密，种群可能陷入局部最优或早熟收敛。基于此，本文设置种群成熟度因子判断GSA算法是否早熟收敛：

Y(t)=(Y1(t)Y2(t)Y3(t)…Yi(t)…YN(t))T

(16)

(17)

对矩阵Y(t)作归一化处理，得到矩阵Y′(t)

(18)

Y′(t) 矩阵中每个行向量可看作一个模糊集合，表示粒子当前位置和历史最优位置在不同维度上的隶属度。Y′(t) 中任意两个模糊集合Y′a(t) 和Y′c(t) 的相似度可以用贴进度σac(t) 表示，记为

(19)

设δ(t) 为种群成熟度因子，它反映的是种群中粒子的平均相似程度，如式(20)所示

(20)

δ(t) 的取值区间为[0,1]。δ(t) 越接近1且大小稳定时，说明在第t次迭代时种群的多样性越差，算法可能已收敛陷入局部最优。

3.3 适应度函数

在智能优化算法中，适应度函数体现算法寻优的目标和方向[10]。本文使用均方误差(mean square error，MSE)作为GSA算法的适应度函数，即改进后聚类的目标函数。均方误差值越小，说明簇内相似度越高，聚类效果越好。本文中均方误差的定义如下所示

(21)

其中，Ai表示第i个聚类中心；Xj是簇Ci内的任意样本点；k是簇的个数。N是数据集样本点的总数。

3.4 改进后K均值聚类的算法流程

为改善传统K均值聚类对初值敏感而易落入局部最优的问题，本文采用具有较强全局搜索能力的GSA算法优化聚类中心。

GSA是单目标优化算法，求解的目标为全局极值，因此本文视一个聚类中心的集合为GSA的优化对象单元。为方便运算，本文建立粒子编码模型表示一个聚类中心集合。具体方法如下：设种群中粒子X由聚类中心ai(i=1,2,3…,n) 组成，用一维数组表示，长度为n×m，m是样本向量的维数。则种群中粒子X的编码结构为

(22)

改进后K均值聚类主要包括3部分。第一部分，设置初始聚类数n0为2，应用GSA全局搜索适应度函数值最小的初始聚心，再引入标准K均值聚类迭代生成簇和更新聚心；第二部分，递增聚类数直至其上限nmax， 对不同的n重复GSA和K均值聚类的工作，最终计算nmax-1个聚类结果的DBI值；第三部分，取最小的DBI值对应的n为最佳聚类数，对应的聚类结果为最佳聚类。

改进后K均值聚类的流程如下所示：

input:N个样本点组成的数据集

步骤2 从数据集中随机选择n个样本作为初始聚类中心，构成一个粒子，粒子形式见式(22)。重复L次，得到L个粒子。对每个粒子使用K均值求聚类，得到L个簇；

步骤3 初始化每个粒子的参数初值；

步骤4 计算所有粒子的惯性质量、所受的作用力和加速度；

步骤5 更新所有粒子的速度和编码值；

步骤6 计算种群成熟度因子δ(t)。 若超过阈值δ′或与δ(t-1) 相差小于0.001，则进入步骤7，否则进入步骤8；

步骤7 选出适应度值最小的粒子，对其用K均值算法求聚类。进入步骤9；

步骤8 判断迭代次数k是否达到最大值kmax， 若是则进入步骤7，否则回到步骤4；

步骤9 计算DBI指数，令n=n+1， 返回步骤2。若n>nmax， 则进入步骤10；

步骤10 选出DBI值最小时的聚类数nbest。 当n=nbest时得到的聚类即算法的输出结果。

output: 两两之间相似度很低的簇。

改进后K均值聚类的伪代码如下所示：

The improved k-means clustering based on GSA

(1) Set the initial value and upper limit of cluster number.

(2) Letnequals 2,selectnsamples randomly to form particles, repeatLtimes.

(3) Identifier particles according to Eq.(22)

(4) Initialize parameters of GSA.

(5) Whilen

(6) Forkfrom 1 tokmaxstep 1

(9) Calculateδ(t) according to Eq.(16)～Eq.(20)

(10) Ifδ(t)>δ′ or |δ(t)-δ(t-1)|≤0.001 then

(11) Break

(12) End if

(13) End for

(14) Cluster the particle withbest_fit(t) according to Eq.(3) and Eq.(21)

(15) Calculate DBI according to Eq.(13)～Eq.(15)

(16)n=n+1

(17) End while

(18) Electnbestwith minimun DBI

(19) Return the cluster whenn=nbest

4 实验验证

本文在Intel(R)Core(TM)i3-2100CPU@3.10GHz，内存12 G，Windows7系统和Matlab2014b的环境下，对最大最小距离聚类、K均值聚类、K-means++聚类和本文改进后K均值聚类进行仿真实验验证。

实验数据采用UCI公开标准数据库中的Iris、Glass、Wine和Balance-scale等4种数据集。其中，Glass属于多类别数据集，Wine属于高维数据集，Balance-scale属于多样本数据集。每种数据集的基本情况见表1。

表1 UCI数据集基本情况

设本文算法的种群粒子总数为50。统计聚类后各簇的个体样本到中心点的距离之和的均值，用来检验算法的寻优能力和稳定性。记录4种算法20次独立运算的最小值(Min)，最大值(Max)，平均值(Ave)和标准差(Std)，结果见表2～表5。

表2 Iris数据集实验结果对比

表3 Glass数据集实验结果对比

表4 Wine数据集实验结果对比

表5 Balance-scale数据集实验结果对比

对于Iris和Balance-scale数据集，本文算法的最小值、最大值、平均值等指标都最小，说明各簇内个体的间距很小，聚类中心的选取合适。本文算法的标准差接近0，仅次于最大最小距离聚类，说明其性能稳定，对多样本数据集鲁棒性较好。

对于Glass数据集，k-means++聚类的最小值最小，为6.4729。但本文算法的最大值、平均值等都远远小于最大最小距离、传统K均值和K-means++聚类，说明本文算法对于多类别数据集亦善于求聚类，且聚类质量较高。

对于Wine数据集，本文算法的平均值最小，标准差也优于K均值和K-means++聚类，说明其对多特征数据集也能取得不错的聚类效果。

为进一步检验本文改进后K均值聚类的有效性，记录4种算法20次独立运算的平均正确率，见表6。

表6 算法平均正确率对比/%

由表6可知，与传统K均值聚类对比，本文算法对4个数据集测试的平均正确率分别提高0.91、1.40、2.25和2.88个百分点。这充分验证了本文的K均值聚类改进策略的有效性。

综上所述，本文改进后K均值聚类对于多样本、多类别和高维数据集等都有良好的搜索性能和稳定性，与最大最小距离、传统K均值和K-means++等聚类对比，其聚类结果更理想。

5 结束语

本文提出一种基于GSA算法的改进K均值聚类。采用粒子编码策略，把聚类中心集合作为GSA算法的寻优对象，克服标准K均值聚类对初始聚心敏感的问题。引入种群成熟度因子，避免算法陷入局部最优。引入聚类质量评价指标，确定最佳聚类数。同时算法把均方误差作为适应度函数，指引全局寻优方向。实验结果表明，本文算法具有更高的正确率和更好的稳定性，对不同类型的数据集都有理想的聚类效果。进一步提高算法的聚类正确率及拓展在一些领域的应用，如改进径向基神经网络辨识算法，将是下一步的研究方向。