0 引言

聚类是机器学习和人工智能的核心领域之一。在聚类算法中，K-means算法具有简单性、有效性和准确性等优点，成为聚类算法中最常用的算法^[1-4]。通过提前给定一个参数K，K-means算法把样例集划分为K个类，每个类的代表对象用这个类里面所有点的中心来表示。与其他聚类算法相比，K-means算法在保持良好的聚类性能的同时，有着更快的速度。在其他聚类算法中，常用K-means算法进行数据的预处理或者后处理，比如谱聚类^[5-7]和聚类集成算法^[8]。

在过去的几十年中，许多研究者致力于改善K-means算法的有效性和准确性^[9-12]，但是K-means算法的一个缺点却很少被关注：由K-means算法产生的聚类经常是不平衡的，即不同类中所含样本的数量的差异会很大。但是在许多聚类问题中，存在着问题域上的平衡先验以及数据收集过程中的平衡偏差，因此在实际应用中，经常期望聚类算法的输出结果是平衡的，即聚类结果中不同类中包含基本相同的样本点个数。例如在对手写数字图像进行聚类时，有理由相信10个数字出现的频率是相同的，所以希望聚类结果中每个类包含基本相同数量的样本。在图像分类领域，词袋模型(bag-of-words)^[13]使用K-means算法把图像特征聚类成不同的类别。文献^[14-16]指出，为了得到更高的图像分类准确率，需要平衡的聚类结果，但是使用K-means常常得不到这样的结果。

由上所述，需要一种能输出平衡聚类结果的聚类算法。本研究重点关注这样的问题：鼓励聚类结果是平衡的，但不严格要求平衡程度，即需要考虑定性的平衡要求，但不是严格的平衡限制(例如把1 000个样本分为10个类，每个类恰好有100个样本)。虽然已经存在一些关于定量限制的平衡聚类方法的研究^[17-18]，但试验证明本研究提出的定性的平衡聚类方法Charl较具有优势。这一算法是由两种经典的聚类算法，Hartigan-Wong算法和Lloyd算法组合而成。

与Lloyd算法相比，Hartigan- Wong算法更倾向于产生较平衡的聚类结果。本研究提出了一种简单有效的调整机制来动态调整聚类结果的平衡性，与Hartigan-Wong算法相比，新调整策略有着更好的灵活性。Charl算法把这种调整机制与Lloyd的批处理更新方式结合起来，确保了Charl算法的效率。试验结果表明:Charl算法不仅获得了更加平衡的聚类结果，也得到更低的代价函数值和更好的聚类性能。

1 Hartigan-Wong算法与Lloyd算法

给定一个拥有n个样本点的数据集，其中每个样本点的格式如下： {x₁,x₂,…,x_n}∈R^d。K-means聚类的任务是把这些点分成K个不同的类，使得相同类中的样本点相似度高，不同类间的样本点相似度低。

1.1 Lloyd算法

形式化来说，K-means算法应该找到所有样本点的分配系数 r_ik∈ {0,1}(1≤i≤n,1≤k≤K)，其中n为样本点的个数，K为聚成的类的个数。并且同时应该找到向量μ_k(1≤k≤K)，每个向量代表一个类的中心。如果x_i属于第k个类，则r_ik=1，否则r_ik=0。由上可知，第k个类的大小为n_k=$\sum\limits_{i=1}^{n}{{{r}_{ik}}}$。

在K-means算法当中，需要最小化代价函数

$J=\sum\limits_{i=1}^{n}{\sum\limits_{k=1}^{K}{{{r}_{ik}}\left\| {{x}_{i}}-{{\mu }_{k}} \right\|}}{{}^{2}}.$

(1)

式(1) 推导得出^[1-2]：

${{r}_{ij}}=\left\{ \begin{array}{*{35}{l}} 1, & \text{if}j=\underset{1\le k\le K}{\mathop{\arg \text{min}}}\,{{\left\| {{x}_{i}}-{{\mu }_{k}} \right\|}^{2}}; \\ 0, & \text{otherwise}. \\ \end{array} \right.$

(2)

${{\mu }_{k}}\frac{\sum\limits_{i=1}^{n}{{{r}_{ik}}{{x}_{i}}}}{\sum\limits_{i=1}^{n}{{{r}_{ik}}}}=\frac{\sum\limits_{i=1}^{n}{{{r}_{ik}}{{x}_{i}}}}{{{n}_{k}}}.$

(3)

Lloyd算法按照批处理的方式，通过两个阶段来找出r_ik和 μ_k的值。

首先固定所有类的中心 μ_k，使用式(2)更新r_ik；然后固定所有r_ik，通过式(3)更新 μ_k。 Lloyd算法反复进行这两个阶段直到结果收敛，保证J收敛到局部极小值^[1-2]。

Lloyd算法简单有效，但该算法的输出经常是不平衡的。对于不同的k(1≤k≤K)，n_k可能会在很大的范围内波动^[14-16]。

1.2 Hartigan-Wong算法

与Lloyd算法相比，Hartigan-Wong算法能得到更低的代价函数值J^[19-21]。与Lloyd的批处理更新不同，Hartigan-Wong是一种在线更新算法。每次迭代过程中，Hartigan-Wong算法先随机选择一个点x_i，然后把x_i从它所属的类中移除，并用式(3)来更新这个类的中心 μ_k。最后通过式(4)把x_i加入到第j个类：

${{r}_{ij}}=\left\{ \begin{array}{*{35}{l}} 1,\text{if}j=\underset{1\le k\le K}{\mathop{\arg \min }}\,\frac{{{n}_{k}}}{{{n}_{k}}+1}{{\left\| {{x}_{i}}-{{\mu }_{k}} \right\|}^{2}}; \\ 0,\text{otherwise}. \\ \end{array} \right.$

(4)

并再次使用式(3)来更新第j个类的中心 μ_j。

这两种算法主要有两个不同点。在线更新方式使Hartigan-Wong算法能获得更低的代价函数值J^[20]，但是Lloyd的批处理更新方式在计算上更高效。

更重要的是，它们在更新r_ik的方式上不同。Lloyd算法使用了平方距离${{\left\| {{x}_{i}}-{{u}_{k}} \right\|}^{2}}$，而Hartigan-Wong算法多了一个乘法项$\frac{{{n}_{k}}}{{{n}_{k}}+1}$。因为$\frac{x}{x+1}$(x>0)是严格递增函数，所以Hartigan-Wong算法更倾向于把样本点分配给包含样本点个数较小的类。举例来说：

${{\left\| {{x}_{i}}-{{\mu }_{1}} \right\|}^{2}}=10,{{n}_{1}}=5.$

(5)

${{\left\| {{x}_{i}}-{{\mu }_{2}} \right\|}^{2}}=9,\text{ }{{n}_{2}}=17.$

(6)

根据式(2)和(4)，Lloyd算法会把x_i分配给第2个类，而Hartigan-Wong算法会把x_i分给第1个类。这个差异的产生是因为n₁＜＜n₂。

需要指出的是，Hartigan-Wong算法中的$\frac{{{n}_{k}}}{{{n}_{k}}+1}$并不是特意为了把样本点分配给样本个数较小的类而设计的。存在这个乘法项，是因为把一个样本点x_i加入到一个合适的类产生了类大小的变化，例如某个类的样本个数的数目从n_k变为n_k+1。

2 2 Charl算法

通过有意放大式(4)的影响，可以便捷地得到一个符合平衡需求的聚类结果。

2.1 新的调整机制

以批处理的方式，使用式(7)把样本点x_i分配到新的类上：

${{r}_{ij}}=\left\{ \begin{array}{*{35}{l}} 1,\text{if}j=\underset{1\le k\le K}{\mathop{\arg \min }}\,{{\left( \frac{{{n}_{k}}}{{{n}_{k}}+1} \right)}^{p}}\left. {{\left\| {{x}_{i}}-{{\mu }_{k}} \right\|}^{2}} \right|; \\ 0,\text{otherwise}\text{.} \\ \end{array} \right.$

(7)

在式(7)中，要求p≥1。也就是说，希望保留并且扩大式(4)把样本点分配给包含样本点个数较少的类的倾向性。如果p＞＞1，那么这种倾向性将会变得更大。p≥1表明，跟其它类相比，即使某个类的中心到x_i的距离比较大，只要它的样本点的个数更少，样本点x_i还是可能被分配到该类。因此，在样本点的重分配过程中，样本点更容易被分配到具有较少样本点的类，因此这种算法会产生更加平衡的聚类结果。

新的调整策略和Lloyd、Hartigan-Wong算法的调整策略的重要区别是，Hartigan-Wong算法和Lloyd算法通过式(2)和式(4)可以保证代价函数J是非递增的，因此这两个算法最终都会收敛。在式(2)中，通过计算J的偏导数，可以很容易验证这个结论，即获得的r_ik保证代价函数J达到了这一阶段的全局最小值。在式(4)中，因为先把样本点x_i从所属的类移出，再将样本点x_i重新加入到某个类，所以也能保证代价函数J达到当前阶段的全局最小值^[21]。

但式(7)并不能保证聚类过程中代价函数随着迭代次数的增加而收敛。在下面的讨论中，将分析何时式(7)会增大代价函数J的值。

2.2 不当样本点的产生

假设样本点x_i在第t次迭代过程中被重分配到了第j个类，即r_ij=0。如果存在s(1≤s≤K)满足式(8)和(9)，则说明由于x_i的重分配，导致关于x_i的代价函数$\left( \sum\limits_{k=1}^{K}{{{r}_{ik}}{{\left\| {{x}_{i}}-{{u}_{k}} \right\|}^{2}}} \right)$的增加：

$J=\underset{1\le k\le K}{\mathop{\arg \min }}\,{{\left( \frac{{{n}_{k}}}{{{n}_{k}}+1} \right)}^{p}}{{\left\| {{x}_{i}}-{{\mu }_{k}} \right\|}^{2}},$

(8)

${{\left\| {{x}_{i}}-{{u}_{j}} \right\|}^{2}}>{{\left\| {{x}_{i}}-{{u}_{s}} \right\|}^{2}}.$

(9)

式(8)表明x_i会被重分配到第j组，式(9)表明代价函数值在重分配之后会增加。也就是说，当存在s满足式(9)和(10)时，产生了一个不当的样本点分配：

${{\left( \frac{{{n}_{j}}}{{{n}_{j}}+1} \right)}^{p}}{{\left\| {{x}_{i}}-{{\mu }_{j}} \right\|}^{2}}\le {{\left( \frac{{{n}_{j}}}{{{n}_{j}}+1} \right)}^{p}}{{\left\| {{x}_{i}}-{{u}_{s}} \right\|}^{2}}.$

(10)

由式(9)(10)解得：

$\frac{{{n}_{s}}\left( {{n}_{j}}+1 \right)}{\left( {{n}_{s}}+1 \right){{n}_{j}}}\ge {{\left\| \frac{{{x}_{i}}-{{\mu }_{j}}}{{{x}_{i}}-{{\mu }_{s}}} \right\|}^{\frac{2}{p}}}>1.$

(11)

定义容忍系数

$\rho =\left( \left\| \frac{{{x}_{i}}-{{u}_{j}}}{{{x}_{i}}-{{u}_{s}}} \right\| \right){{}^{\frac{2}{p}}},$

(12)

且在实际情况中n_j>0，则：

${{n}_{s}}\left( {{n}_{j}}-\rho {{n}_{j}}+1 \right)\ge \rho {{n}_{j}}>0.$

(13)

结合式(11)(12)，可以进一步得到：

$1＜\rho ＜1+\frac{1}{{{n}_{j}}},$

(14)

${{n}_{j}}＜\rho {{n}_{j}}＜{{n}_{j}}+1,$

(15)

$\frac{{{n}_{t}}}{{{n}_{j}}}\ge \frac{\rho }{{{n}_{j}}+1-\rho {{n}_{j}}}>1.$

(16)

容易证明，式(14)(15)(16)与式(9)(10)等价,且式(14)(15)(16)是发生不当分配的必要条件。

2.3 不当样本点分配的减少

如果在第t次迭代中，式(7)为x_i重分配到第j个类，那么在第t次迭代：

(1) 先根据式(12)计算ρ，此过程与t无关。

(2) 如果存在s(1≤s≤K)，同时满足式(14)和(16)，那么这次对x_i的重分配就会导致代价函数值J的增加。

通过分析式(8)~(16)，可以得到如下的推论：

(1) 根据式(8)和(9)，当p增加时，Charl算法更倾向于将样本点x_i分配到样本个数较少的类，但代价函数增大的风险也随之增大。当需要聚类结果更平衡时，可以令p为较大的值，反之，可令p为较小的值。

(2) 根据式(13)，当p增加时，ρ会减少，式(14)更容易被满足。因此较大的p更容易导致不当的样本点分配。

(3) 当分类结果接近均匀时，对于所有的s(1≤s≤K)，$\frac{{{n}_{s}}}{{{n}_{j}}}\approx 1$。此时并不需要聚类结果变得更加平衡，可将p的值设置为一个较小值。

(4) 如果式(14)成立，那么式(16)中的$\frac{\rho }{{{n}_{j}}+1-\rho {{n}_{j}}}$取决于ρn_j和n_j+1的接近程度。当ρn_j接近n_j+1时，产生不当的样本点分配的概率也随之增大。

(5) 最后需要强调的是，尽管在某样本点的分配上可能会使J增大，但在一次迭代过程中，所有样本点分配完成后，整体的代价函数值J仍然可能是下降的。

2.4 平衡程度评价指标和自适应调整

由前面分析可知，本研究需要一种评价聚类平衡性的指标，也需要一种动态调整p的方式。本研究使用聚类结果中不同类的标准差作为平衡程度的评价指标，计算式为：

$\sigma =\sqrt{\frac{1}{K}{{\sum\limits_{k=1}^{K}{\left( {{n}_{k}}-\frac{n}{K} \right)}}^{2}}}.$

(17)

如果聚类结果不平衡，则σ较大；如果聚类结果绝对平衡，则σ=0。由于不强制要求聚类结果绝对平衡，所以当聚类结果的σ≤$\frac{n}{5K}$时，认为这个聚类结果是平衡的，反之，则为不平衡。定义b=$\frac{\sigma }{n\cdot {{K}^{-1}}}$为平衡因子，则在本研究中定义当且仅当b≤0.2时，该聚类结果是平衡的。对于其他数据集，可以根据期望得到的聚类结果的平衡程度来确定平衡因子b。

动态调整p的方案如下:首先，初始化p=1，如果聚类结果不平衡(b≥0.2)，用p=1.25p来扩大p。如果聚类结果是平衡的，用p=p/1.25来缩小p。但是如果b≤0.2并且p＜1.25，则将p置成1。根据上述描述，Charl算法的具体过程如下。

1: 输入：x₁,x₂,…,x_n,K,Tmax。

2: 随机选择K个点作为初始中心μ_k(1≤k≤K)。

3: t→0，使用式(1)计算J_t。

4: p→1。

5: 若迭代次数没有达到T_max，进行以下迭代：

6: 使用式(7)更新r_ik(1≤i≤n,1≤k≤K)。

7: 使用式(3)更新聚类中心 μ_k。

8: 重新计算类的大小n_k，计算标准差

${{\sigma }_{t}}=\sqrt{\frac{1}{K}\sum\limits_{k=1}^{K}{{{\left( {{n}_{k}}-\frac{n}{K} \right)}^{2}}}}.$

10: 如果σ_t>$\frac{n}{5K}$，则p→1.25p。

11: 否则p→max(p/1.25,1)。

12: t→t+1,使用式(1)计算J_t。

13: 如果J_t≥J_t-1，终止迭代。

Charl算法采用批处理的更新方式，所以具有跟Lloyd一样的高效计算的优点。在每次迭代过程中，根据聚类结果是否均衡，自适应地调整p。通过这一策略，Charl算法倾向于得到更平衡的聚类结果。在满足以下两个条件之一时，终止Charl算法的迭代过程：

(1) Charl算法迭代数达到最大迭代次数。

(2) 某次迭代后，代价函数值J上升。

2.5 Charl算法与定量平衡聚类算法对比

Quant算法是一个典型的定量平衡聚类算法^[10],现在将Charl算法进行对比。Quant算法可以指定最后聚类结果必须是完全平衡的，即n_k≥$\left\lceil \frac{n}{K} \right\rceil ,1\le k\le K$。

为了满足定量平衡的性质，Quant算法不得不牺牲一些其他的性能，比如效率和准确率。而与之相对应的是，Charl算法使用了启发式算法来鼓励平衡聚类。在下一节，将通过试验对比Charl算法与Quant算法。

3 试验

为了验证Charl算法的有效性，首先在13个常用的聚类数据集上对比Charl算法、Lloyd算法和Quant算法的各项指标：平衡因子b、代价函数值J、聚类准确率(accuracy,ACC)、归一化互信息(normalized mutual information,NMI)和调整随机指标(adjusted rand index,ARI)、运行时间。

3.1 数据集、评估参数和设定

试验使用了具有平衡先验的13组数据集，数据集的具体属性在表 1中给出。对于RCV数据集，使用了4个种类：“C15”“ECAT”“GCAT”和“MCAT”中的样本点。对Isolet数据集，使用了Isolet1子数据集。

这些数据集的属性差异较大。这些数据集中，样本点个数的范围是165~70 000个，样本点的维度范围是617~29 992。数据集涉及到不同的领域，包括图像、音频和文本。数据集中类的个数也不相同，范围是4~68。因为Charl算法适用的问题域是具有平衡先验的，所以选择这些数据集的一个重要特征是，它们的类别全部是平衡的(b≤0.2)。

所有的样本点都进行行归一化，即，x_i→$\frac{{{x}_{i}}}{\left\| {{x}_{i}} \right\|}$。考虑到主成分分析(principal component analysis,PCA)已经被广泛应用于聚类数据的预处理中，所以在运行聚类算法之前，均对数据进行PCA预处理，保留累计贡献率大于90%时对应的主成分个数^[2]。每种算法在每个数据集上均运行100次，并报告100次运行后算法平均表现。在运行Lloyd算法和Charl算法时，设置最大的迭代次数是100(T=100)，并且为了保证Charl算法和Lloyd算法的公平对比，两种算法使用相同的初始化值。在Quant算法中，根据文献^[10]的建议，首先采样大小为$\left\lfloor \frac{n}{K} \right\rfloor $的样本点子集，并且使用Lloyd算法对样本点子集进行聚类。 在得到聚类结果之后，首先通过Hungarian算法^[22]找到聚类结果的类别和真正类别之间的最佳二分图匹配。 在获得最佳匹配之后，使用ACC、 NMI、 ARI等指标来评价聚类结果^[22-24]。

3.2 Charl算法聚类结果

首先，对比Charl算法和Lloyd算法的聚类结果的平衡程度。表 2中列出了进行100次聚类试验得到的平均平衡因子b和平均代价函数。在本研究中，使用↓表示所比较的算法的表现劣于Charl的算法(例如更不平衡或者是代价函数值更高)，使用↑表示所比较的算法优于Charl算法。

试验结果显示，在全部数据集中，Charl算法都比Lloyd算法得到了更低的平衡因子b。因此，Charl算法可以得到更加平衡的聚类结果。本研究计算了Charl算法和Lloyd算法的平衡因子的相对减少程度

$\frac{{{b}_{\text{Lloyd}}}-{{b}_{\text{Charl}}}}{{{b}_{\text{Lloyd}}}}.$

相对减少程度从3.6%到35.9%不等，平均为20%。这说明了Charl算法达到比Lloyd算法更平衡的聚类结果。

表 2可以看出，与Lloyd算法和Charl算法相比，Quant算法能得到明显更平衡的聚类结果。这是由于Quant算法的设计目标是实现一个定量平衡的聚类算法。但是这并不意味着Quant算法是一个较好的聚类算法。

3.3 Charl算法代价函数值

出乎本研究意料的是，在8个数据集上，Charl算法获得了比Lloyd算法更低的代价函数值 J(如表 2所示)。因为设计Lloyd算法的初衷是得到尽可能低的代价函数值J，Charl算法则不是。但是试验结果表明，Charl算法却能得到更低的代价函数。造成这一现象的原因为： Lloyd算法只能收敛到一个局部最优解，而在p较大时，Charl算法并不能保证最后一定收敛。在有平衡先验的问题域上，Charl算法更倾向于收敛到有较低代价函数值J的局部最优解，因为Charl算法更倾向于寻找平衡的聚类结果。而Lloyd算法忽略了数据集本身的平衡先验知识，可能收敛到一个不平衡的解，也就意味着得到了较高的代价函数值J。更详细的原因将会在未来的工作继续深入研究。

在所有的数据集上，Charl算法都比Quant算法获得了更低的代价函数值，这是因为Quant算法必须牺牲聚类结果的质量来保证聚类结果的严格平衡。

3.4 Charl算法的聚类性能

Charl算法不仅能获得更平衡的聚类结果，还可以得到较好的聚类性能，例如ACC、NMI和ARI。表 3统计了运行100次Charl、Lloyd和Quant算法得到的平均ACC、NMI和ARI值。

对于ACC指标，在所有的13个数据集上，Charl算法都优于Lloyd算法。Charl算法的平均ACC比Lloyd算法高4.4%。总体来说，在平衡聚类问题中，Charl算法比Lloyd算法有更高的ACC，而且Charl算法在和Quant算法相比较时也有足够的竞争力。Charl算法在10个数据集上胜过Quant算法，并且ACC相对提高了2%以上。在另外3个数据集上，Quant算法虽然优于Charl算法，但是没有明显的提升。

对于NMI和ARI指标，Charl算法分别在11个数据集和10个数据集上优于Lloyd算法，分别在12个数据集和13个数据集上优于Quant算法。

简而言之，对比多次重复试验时聚类算法所能达到的最优聚类性能，Charl算法的性能仍然明显优于另外两种算法。

最后值得注意的是，虽然平衡先验存在于很多问题域上，但是在聚类过程中要求严格的平衡，得到的结果往往要劣于本研究提出的Charl算法得到的结果。

3.5 Charl算法的实际应用

Charl算法作为一种平衡聚类算法，可以应用到许多场景中。最近邻问题是计算机视觉中一个基础且重要的问题。在数据维数很大的时候，由于得到精确的最近邻点的代价太大，通常使用近似算法得到近似最近邻点。而Product Quantization算法^[25]是解决近似最近邻的一种常用算法。

本研究使用的3个数据集分别为SIFT10K、SIFT1M、GIST1M。采取Product Quantization算法^[25]建议的参数，将高维向量划分为8个子向量。对每个子向量组聚类时，类的个数设为256。每个算法在所有数据集上都运行100次，其中Charl算法使用的初始聚类中心和Lloyd算法相同。评价指标为recall@R。其中R值表示返回的近似最近邻的个数，recall指标表示返回的近似最近邻是真正的最近邻的比例。

通过表 4可知，在全部数据集上，Charl算法都能比Lloyd算法得到更高的召回率，从而进一步验证了Charl算法的有效性。

4 结论与展望

本研究通过结合Lloyd算法和Hartigan-Wong算法，提出一种新的定性平衡聚类算法Charl，来获得平衡的聚类结果。在Charl算法中，本研究提出了一种新的动态调整机制，使得Charl算法的聚类结果中不同的类有着基本相同的大小。在13个聚类数据集上进行的试验表明:Charl算法能达到更好的平衡性和聚类性能。同时，本研究把Charl算法与一种严格平衡的聚类算法进行对比，发现Charl算法能得到更低的代价函数值和更好的聚类性能。

目前，在某些问题域上，可能存在输出类别大小分布的先验知识，本研究对于此类问题的聚类比较感兴趣，将会继续研究这方面的工作。例如医疗诊断中，待聚类的图像正负样本个数差异很大，研究中期待得到一个较为不平衡的聚类结果。另外，也考虑继续拓展Charl算法来满足某些定量平衡的要求。