A Simple Linear Time (1 + ε)-Approximation Algorithm for k-Means Clustering in Any Dimensions

第三十九次写博客，本人数学基础不是太好，如果有幸能得到读者指正，感激不尽，希望能借此机会向大家学习。由于K-Means算法本身的时间复杂度很高，特别是在处理大数据集时，因此这篇文章主要介绍一种具有线性时间复杂度的近似K-Means算法，其他有关于原型聚类算法的讨论可以移步到该类算法的导航页《原型聚类算法综述（原型聚类算法开篇）》。

K-Means问题定义？

给定数据集P∈Rn×dP\in{\Bbb{R}^{n\times{d}}}P∈Rn×d、质心（centroid）集合K∈Rk×dK\in{\Bbb{R}^{k\times{d}}}K∈Rk×d，KKK中的每个样本点都是某个数据子集的质心，那么K-Means问题可以表示为，在给定数据集PPP以及聚类个数kkk时，求得可以最小化下式的质心集合，
∑p∈Pd(p,K)2\sum_{p\in{P}}{d\left({\bf{p}},K\right)^{2}}p∈P∑d(p,K)2
其中，d(p,K)d\left({\bf{p}},K\right)d(p,K)是指样本点p\bf{p}p与集合KKK中距离其最近的质心之间的距离。
给定数据集P∈Rn×dP\in{\Bbb{R}^{n\times{d}}}P∈Rn×d、中位点（median）集合K∈Rk×dK\in{\Bbb{R}^{k\times{d}}}K∈Rk×d，KKK中的每个样本点都是某个数据子集的中位点，那么K-Median问题可以表示为，在给定数据集PPP以及聚类个数kkk时，求得可以最小化下式的中位点集合，
∑p∈Pd(p,K)\sum_{p\in{P}}{d\left({\bf{p}},K\right)}p∈P∑d(p,K)
其中，d(p,K)d\left({\bf{p}},K\right)d(p,K)是指样本点p\bf{p}p与集合KKK中距离其最近的中位点之间的距离。

本文提出的算法？

本文基于“随机采样”（Random Sampling）技术，得到了一种针对K-Means问题的、时间复杂度为线性的（具体来说是与输入大小呈线性的，即O(2(k/ε)O(1)nd)O\left(2^{\left(k/\varepsilon\right)^{O\left(1\right)}}nd\right)O(2(k/ε)O(1)nd)）“(1+ε)−\left(1+\varepsilon\right)-(1+ε)−近似”（(1+ε)−\left(1+\varepsilon\right)-(1+ε)−approximate）算法。针对于2-Means聚类问题，本文还提出一个在恒定不变的概率下，以O(2(1/ε)O(1)nd)O\left(2^{\left(1/\varepsilon\right)^{O\left(1\right)}}nd\right)O(2(1/ε)O(1)nd)时间复杂度运行的(1+ε)−\left(1+\varepsilon\right)-(1+ε)−近似算法。

(k,ε)−\left(k,\varepsilon\right)-(k,ε)−不可化简与(k,ε)−\left(k,\varepsilon\right)-(k,ε)−可化简？

首先，令Δ(P,K)\Delta\left(P,K\right)Δ(P,K)表示K-Means问题中的目标函数，令Δk(P)\Delta_{k}\left(P\right)Δk(P)表示该问题的最优解，引出下述定义，
定义1 如果数据集PPP满足不等式
Δk−1(P)≥(1+32ε)Δk(P)\Delta_{k-1}\left(P\right)\geq{\left(1+32\varepsilon\right)\Delta_{k}\left(P\right)}Δk−1(P)≥(1+32ε)Δk(P)
那么就称PPP是“(k,ε)−\left(k,\varepsilon\right)-(k,ε)−不可化简”的，否则称其为“(k,ε)−\left(k,\varepsilon\right)-(k,ε)−可化简”的。
直观上讲，可化简是指能否找到一个最优解，这个解比当前最优解的簇个数小，并且更接近数据的真实聚类结构。

1-Means问题的一些性质？

首先给出下述定义，

┌定义2 给定样本点集合PPP，PPP的“质心”可以表示为c(P)c\left(P\right)c(P)，即c(P)=∑p∈Pp∣P∣c\left(P\right)=\frac{\sum_{\bf{p}\in{P}}{\bf{p}}}{|P|}c(P)=∣P∣∑p∈Pp（时间复杂度为O(nd)O\left(nd\right)O(nd)），其中∣P∣|P|∣P∣表示PPP中所含样本点个数，那么对于任意ddd维向量x\bf{x}x，下述等式成立
Δ(P,x)=Δ(P,c(P))+∣P∣⋅Δ(c(P),x)\Delta\left(P,\bf{x}\right)=\Delta\left(P,c\left(P\right)\right)+|P|\cdot{\Delta\left(c\left(P\right),\bf{x}\right)}Δ(P,x)=Δ(P,c(P))+∣P∣⋅Δ(c(P),x)┘

根据上述定义，我们可以得到下述结论，

┌结论1 对于任何数据集PPP，其1-Kmeans问题的最优解一定是其质心c(P)c\left(P\right)c(P)，因为这可以使得定义2中的等式右边的第二项为零。┘

上述结论可以扩展到K-Means问题中，对于簇质心集合K={x1,x2,…,xk}K=\{{\bf{x}}_{1},{\bf{x}}_{2},\dots,{\bf{x}}_{k}\}K={x1,x2,…,xk}，数据集P={P1,P2,…,Pk}P=\{P_{1},P_{2},\dots,P_{k}\}P={P1,P2,…,Pk}中子集PiP_{i}Pi是所有隶属于第iii个簇的样本点集合，由上述结论可知，xi{\bf{x}}_{i}xi是数据子集PiP_{i}Pi的质心。下面给出一个引理，

┌引理1 在原始数据集PPP上进行独立的均匀随机采样，得到样本大小为mmm的数据子集TTT，那么对于任意δ>0\delta\gt{0}δ>0，如下不等式以至少1−δ1-\delta1−δ的概率成立，
Δ(P,c(T))<(1+1δm)Δ1(P)\Delta\left(P,c\left(T\right)\right)\lt{\left(1+\frac{1}{\delta{m}}\right)\Delta_{1}\left(P\right)}Δ(P,c(T))<(1+δm1)Δ1(P)┘

直观来讲，原始数据集上通过随机采样得到的，某个较小的子集的质心是原始数据集质心的一个较好的近似。令δ=12\delta=\frac{1}{2}δ=21、m=2εm=\frac{2}{\varepsilon}m=ε2，根据上述引理，“通过随机采样得到的数据子集TTT的质心是Δ1(P)\Delta_{1}\left(P\right)Δ1(P)的(1+ε)−\left(1+\varepsilon\right)-(1+ε)−近似”成立的概率至少为12\frac{1}{2}21。因此，通过随机采样可以迅速得到1-Means问题最优解的较好的近似解。

假设数据集P′P'P′是PPP的某个不确定的子集，即∣P′∣≥β∣P∣|P'|\geq{\beta|P|}∣P′∣≥β∣P∣（0≤β≤10\leq{\beta\leq{1}}0≤β≤1），数据集SSS是在PPP上通过随机采样得到的、大小为4βε\frac{4}{\beta\varepsilon}βε4的数据子集，令S′S'S′表示SSS的任意一个大小为2ε\frac{2}{\varepsilon}ε2的子集，我们为每个S′S'S′计算Δ(P′,c(S′))\Delta\left(P',c\left(S'\right)\right)Δ(P′,c(S′))，下述引理表示，SSS中必定存在一个子集S′S'S′，使得c(S′)c\left(S'\right)c(S′)是P′P'P′的1-Means最优解的较好近似。

┌引理2 （超集采样引理（Superset Sampling Lemma））下述不等式以恒定不变的概率成立，
min⁡S′:S′⊂S,∣S′∣=2εΔ(P′,c(S′))≤(1+ε)Δ1(P′)\min_{S':S'\subset{S},|S'|=\frac{2}{\varepsilon}}{\Delta\left(P',c\left(S'\right)\right)}\leq{\left(1+\varepsilon\right)\Delta_{1}\left(P'\right)}S′:S′⊂S,∣S′∣=ε2minΔ(P′,c(S′))≤(1+ε)Δ1(P′)┘

根据上述采样方案得到的数据集SSS以恒定不变的概率至少包含P′P'P′中的2ε\frac{2}{\varepsilon}ε2个样本点，因此根据引理1得知上述引理成立。

2-Means问题的一个线性时间复杂度算法？

在这一部分，本文给出了一个针对于2-Means聚类问题的算法，假设对于数据集P∈Rn×dP\in{\Bbb{R}^{n\times{d}}}P∈Rn×d，c1{\bf{c}}_1c1和c2{\bf{c}}_2c2分别是两个簇质心，P1P_{1}P1和P2P_{2}P2分别是这两个质心所对应的簇集合，首先不失一般性的假设∣P1∣≥∣P∣/2|P_{1}|\geq{|P|/2}∣P1∣≥∣P∣/2，那么根据引理2，我们可以通过在集合PPP上通过采样得到大小为O(1ε)O\left(\frac{1}{\varepsilon}\right)O(ε1)的数据集SSS，以此得到c1{\bf{c}}_1c1的较好近似c1′{\bf{c}}_{1}'c1′，这一步不需要依赖于任何假设。下面假设Qi′Q_{i}'Qi′是PPP中与c1′{\bf{c}}_{1}'c1′距离最远的前n2i−1\frac{n}{2^{i-1}}2i−1n（i=1,2,…i=1,2,\dotsi=1,2,…）个样本点集合，α=ε/64\alpha=\varepsilon/64α=ε/64（在证明“当PPP满足(k,α)−\left(k,\alpha\right)-(k,α)−可化简时满足结论2”的过程中给出的假设条件），那么∣P2∣≥α∣Qi′∣/2|P_{2}|\geq{\alpha|Q_{i}'|/2}∣P2∣≥α∣Qi′∣/2，根据引理2，我们可以通过在Qi′Q_{i}'Qi′上通过采样得到大小为O(1α2)O\left(\frac{1}{\alpha^2}\right)O(α21)的数据集SSS，以此得到c2{\bf{c}}_2c2的较好近似c2′{\bf{c}}_2'c2′，由于并不知道iii的确切取值，因此需要对所有可能的iii值进行迭代，这样整个算法的时间复杂度是O(2(1/ε)O(1)nd)O\left(2^{\left(1/\varepsilon\right)^{O\left(1\right)}}nd\right)O(2(1/ε)O(1)nd)。
并通过一系列证明，得到了下述关于该算法的结论，

┌结论2 给定数据集P∈Rn×dP\in{\Bbb{R}^{n\times{d}}}P∈Rn×d，采用上述算法可以的得到2-Means问题最优解的一个(1+ε)−\left(1+\varepsilon\right)-(1+ε)−近似解（以某一恒定不变的概率），并且该算法的时间复杂度是O(2(1/ε)O(1)nd)O\left(2^{\left(1/\varepsilon\right)^{O\left(1\right)}}nd\right)O(2(1/ε)O(1)nd)。┘

下面简要对证明过程进行解释，当数据集PPP满足(k,α)−\left(k,\alpha\right)-(k,α)−可化简时（其中α=ε/64\alpha=\varepsilon/64α=ε/64），可以将该问题转化为1-Means问题并通过简单的求取PPP的质心作为最优解，这时时间复杂度是O(nd)O\left(nd\right)O(nd)，因此满足结论2。当数据集PPP不满足(k,α)−\left(k,\alpha\right)-(k,α)−可化简（即满足(k,α)−\left(k,\alpha\right)-(k,α)−不可化简）时，文章首先给出假设d(c1,c2)=td\left({\bf{c}}_1,{\bf{c}}_2\right)=td(c1,c2)=t，并给出如下引理

┌引理3 d(c1,c1′)≤t/4d\left({\bf{c}}_1,{\bf{c}}_1'\right)\leq{t/4}d(c1,c1′)≤t/4┘

这一引理是根据定义2中的等式，使用反证法证明（由于d(c1,c1′)≤t/4d\left({\bf{c}}_1,{\bf{c}}_1'\right)\leq{t/4}d(c1,c1′)≤t/4时，违背PPP的(k,α)−\left(k,\alpha\right)-(k,α)−不可化简性质）得到的，根据该引理可以认为以c1′{\bf{c}}_1'c1′为中心的球B(c1′,t/4)B\left({\bf{c}}_1',t/4\right)B(c1′,t/4)在c1{\bf{c}}_1c1代表的聚类簇中，即B(c1′,t/4)B\left({\bf{c}}_1',t/4\right)B(c1′,t/4)中的这些点均属于P1P_1P1，如下图所示（二维情况下的简单示意图），

图1 示例图

下面令P′=P1′+P2P'=P_1'+P_2P′=P1′+P2（上图阴影部分），其中P1′=P1−B(c1′,t/4)P_1'=P_1-B\left({\bf{c}}_1',t/4\right)P1′=P1−B(c1′,t/4)，并给出结论3，

┌结论3 ∣P2∣≥α∣P1′∣|P_2|\geq{\alpha|P_1'|}∣P2∣≥α∣P1′∣┘

这一结论依旧是使用反证法进行证明（由于∣P2∣≥α∣P1′∣|P_2|\geq{\alpha|P_1'|}∣P2∣≥α∣P1′∣时，违背PPP的(k,α)−\left(k,\alpha\right)-(k,α)−不可化简性质），并且用到了引理3中的不等式。最后将该结论的边界进行适当放缩，并根据引理2，可以得到“在P′P'P′上运行超集采样引理可以得到c2{\bf{c}}_2c2的一个较好近似c2′{\bf{c}}_2'c2′”，将该结论与“在集合PPP上通过采样得到大小为O(1ε)O\left(\frac{1}{\varepsilon}\right)O(ε1)的数据集SSS，以此得到c1{\bf{c}}_1c1的较好近似c1′{\bf{c}}_1'c1′”结合可以得到2-Means问题的一个(1+ε)−\left(1+\varepsilon\right)-(1+ε)−近似算法，因此，结论2得证。

一种线性时间复杂度的K-Means聚类算法？

首先给出该算法的整体框架（如下图所示），该算法的输入是数据集PPP、聚类个数kkk以及近似率ε\varepsilonε，该算法考虑所有可能的聚类个数i=1,2,…i=1,2,\dotsi=1,2,…，并将其中最小的Δi(P)\Delta_{i}\left(P\right)Δi(P)对应的聚类簇作为最终的聚类结果，下面对Irred-k-means过程进行详细的讲解。

图2 算法框架

如下图所示，是Irred-k-means过程的伪代码，该过程的输入分别是，原始数据集PPP中未被指派簇隶属的样本点集合QQQ、还未被发现的聚类个数mmm、聚类个数kkk、所有已经被发现的聚类的质心集合CCC（∣C∣=k−m|C|=k-m∣C∣=k−m）、近似因子α\alphaα以及当前阶段的整体成本（即将P−QP-QP−Q中的样本点指派给CCC所产生的成本），该过程最终产生具有最低整体成本的kkk个聚类簇。

图3 Irred-k-means伪代码

在步骤2中，主要是通过对数据集QQQ运行引理2所示的随机采样方法，得到一个新的质心，具体方法是在QQQ中随机采样得到一个大小为O(k/α2)O\left(k/{\alpha^2}\right)O(k/α2)的样本集合SSS，针对SSS中所有大小为O(1/α)O\left(1/\alpha\right)O(1/α)的子集S′S'S′求得集合的质心c\bf{c}c，然后将c\bf{c}c加入到集合CCC中，并将mmm减一，需要注意的是，这一步并不会将c\bf{c}c代表的聚类集合从数据集QQQ中剔除。然后将步骤2中的参数传递给步骤3，该步骤将QQQ中的样本点根据其与CCC中各质心的距离之和进行升序排序，并取其中的前∣Q∣/2|Q|/2∣Q∣/2个样本点组成集合UUU，然后将UUU中所有的样本点指派给CCC中距离其最近的质心，最终将UUU从数据集QQQ中剔除，并更新当前整体成本，至此，“发现一个新的质心”的过程才算结束，将更新后的参数带入步骤2中继续发现新的质心。当所有质心已经被找到，即m=0m=0m=0时，就将QQQ中的剩余样本点指派给CCC中的各个质心，并更新整体成本，将该聚类结果加入到备选中（步骤1），最终在这些备选中选择一个整体成本最小的聚类结果作为该过程的输出（步骤4）。
通过上述伪代码可以看出，步骤2、3与本文之前提出的2-Means聚类问题的算法类似，均利用了“超集采样引理”。

参考资料

【1】 Kumar, A. , Y. Sabharwal , and S. Sen . “A simple linear time (1 + ϵ)-approximation algorithm for k-means clustering in any dimensions.” IEEE Symposium on Foundations of Computer Science IEEE, 2004.