AAAI-2016《Approximate K-Means++ in Sublinear Time》

推荐深蓝学院的《深度神经网络加速：cuDNN 与 TensorRT》，课程面向就业，细致讲解CUDA运算的理论支撑与实践，学完可以系统化掌握CUDA基础编程知识以及TensorRT实战，并且能够利用GPU开发高性能、高并发的软件系统，感兴趣可以直接看看链接：深蓝学院《深度神经网络加速：cuDNN 与 TensorRT》

---

核心思想

论文的核心思想是针对传统$k$-means++算法在大规模数据集上的初始化效率瓶颈，提出了一种基于马尔可夫链蒙特卡洛（MCMC）采样的快速初始化算法，称为K-MC²。传统$k$-means++通过$D^2$采样选择初始中心以提高聚类质量，但需要$k$次全数据集扫描，时间复杂度为$\Theta (nkd)$，对大规模数据（如百万或亿级点）计算成本过高。K-MC²通过用MCMC方法近似$D^2$采样，显著降低计算复杂度至亚线性（$\mathcal{O}(k^{3}d{\log }^{2}n\log k)$），且在非病态数据集（满足温和分布假设）下保留$k$-means++的$\mathcal{O}(\log k)$近似保证。该算法特别适用于需要快速初始化的场景（如在线聚类、mini-batch $k$-means或coreset构造），并通过理论分析和实验验证其高效性和实用性。

---

目标函数

$k$-means聚类的目标是通过最小化量化误差（即平方误差和，Sum of Squared Error, SSE）将数据集 X = { x 1 , x 2 , … , x n } \mathcal{X} = \{x_1, x_2, \ldots, x_n\} X={x1​,x2​,…,xn​}（$x_i\in {\mathbb{R}}^d$）划分为$k$个簇。给定中心集合 C = { c 1 , c 2 , … , c k } C = \{c_1, c_2, \ldots, c_k\} C={c1​,c2​,…,ck​}，量化误差定义为：

$${\varphi }_C(\mathcal{X})=\sum _{x\in \mathcal{X}}\mathrm{d}(x,C{)}^2=\sum _{x\in \mathcal{X}}\underset{c\in C}{\min }\Vert x-c{\Vert }_2^2$$

其中，$\mathrm{d}(x,C)={\min }_{c\in C}\Vert x-c{\Vert }_2$表示点$x$到最近中心的欧几里得距离。令${\varphi }_{OPT}^k(\mathcal{X})$表示最优$k$中心解的量化误差，一个解$C$被称为$\alpha$近似解，若满足：

$${\varphi }_C(\mathcal{X})\le \alpha \cdot {\varphi }_{OPT}^k(\mathcal{X})$$

$k$-means++通过$D^2$采样初始化，期望上达到$\mathcal{O}(\log k)$近似，即：

$$\mathbb{E}[{\varphi }_{C_k}(\mathcal{X})]\le 8({\log }_{2}k+2){\varphi }_{OPT}^k(\mathcal{X})$$

K-MC²的目标是近似$k$-means++的初始化过程，保留相同的目标函数和$\mathcal{O}(\log k)$近似保证，同时降低初始化时间复杂度。

---

目标函数的优化过程

K-MC²算法通过MCMC近似$D^2$采样来优化初始中心选择，随后结合标准的$k$-means迭代（Lloyd算法）最小化${\varphi }_C(\mathcal{X})$。优化过程分为以下步骤：

1. 初始中心选择（均匀采样）：

   • 从数据集$\mathcal{X}$中均匀随机选择第一个中心$c_1$，加入中心集合 C 1 = { c 1 } C_1 = \{c_1\} C1​={c1​}。
   • 时间复杂度：$\mathcal{O}(1)$。

2. MCMC近似$D^2$采样（K-MC²初始化）：

   • 对于剩余的$k-1$个中心（$i=2,\dots ,k$）：
     • 初始化一个马尔可夫链，均匀随机选择起点$x$，计算其到当前中心集$C_{i-1}$的平方距离$d_x=\mathrm{d}(x,C_{i-1}{)}^2$。
     • 运行$m$步Metropolis-Hastings算法：
       • 均匀随机选择候选点$y$，计算$d_y=\mathrm{d}(y,C_{i-1}{)}^2$。
       • 计算接受概率：
         $$\pi =\min \left(1,\frac{p(y)}{p(x)}\right)=\min \left(1,\frac{d_y}{d_x}\right)$$
       • 以概率$\pi$接受$y$（设置$x\leftarrow y$，$d_x\leftarrow d_y$），否则保留$x$。
     • 取马尔可夫链第$m$步的状态$x_m$作为新中心$c_i$，更新 C i = C i − 1 ∪ { c i } C_i = C_{i-1} \cup \{c_i\} Ci​=Ci−1​∪{ci​}。
   • 时间复杂度：每步计算距离需要$\mathcal{O}(kd)$，$m$步共$\mathcal{O}(mkd)$，$k-1$次迭代共$\mathcal{O}(m{k}^{2}d)$。
   • 理论保证：马尔可夫链的平稳分布为$D^2$采样分布$p(x)=\frac{\mathrm{d}(x,C_{i-1}{)}^2}{{\sum }_{x^{\prime }\in \mathcal{X}}\mathrm{d}(x^{\prime },C_{i-1}{)}^2}$，总变差距离随$m$几何收敛：
     $$\Vert {\tilde{p}}_m-p{\Vert }_{TV}=\mathcal{O}\left({\left(1-\frac{1}{\gamma }\right)}^m\right),\gamma =\underset{x\in \mathcal{X}}{\max }p(x)$$
     选择$m=\mathcal{O}({\gamma }^{\prime }\log \frac{k}{ϵ})$，可使总变差距离$\Vert {\tilde{p}}_m-p{\Vert }_{TV}\le \frac{ϵ}{k-1}$，其中${\gamma }^{\prime }={\max }_{C\subset \mathcal{X},|C|\le k}{\max }_{x\in \mathcal{X}}n\frac{\mathrm{d}(x,C{)}^2}{{\sum }_{x^{\prime }\in \mathcal{X}}\mathrm{d}(x^{\prime },C{)}^2}$。

3. 数据集假设与亚线性复杂度：

   • 假设数据集$\mathcal{X}$从分布$F$独立同分布采样，满足：
     • (A1)：$F$具有有限方差和指数尾（如高斯、指数、拉普拉斯分布），则：
       $$\alpha =\frac{{\max }_{x\in \mathcal{X}}\mathrm{d}(x,\mu (\mathcal{X}){)}^2}{\frac{1}{n}{\sum }_{x^{\prime }\in \mathcal{X}}\mathrm{d}(x^{\prime },\mu (\mathcal{X}){)}^2}=\mathcal{O}({\log }^{2}n)$$
     • (A2)：$F$为非退化分布（如超球面上的近似均匀分布），则：
       $$\beta =\frac{{\varphi }_{OPT}^1(\mathcal{X})}{{\varphi }_{OPT}^k(\mathcal{X})}=\mathcal{O}(k)$$
   • 因此，${\gamma }^{\prime }\le 4\alpha \beta =\mathcal{O}(k{\log }^{2}n)$，链长$m=\mathcal{O}(k{\log }^{2}n\log k)$，总复杂度为：
     $$\mathcal{O}(k^{3}d{\log }^{2}n\log k)$$
     这与$n$亚线性相关。

4. $k$-means迭代：

   • 使用K-MC²选择的初始中心$C_k$，运行Lloyd算法：
     • 将每个点分配到最近中心，更新簇。
     • 计算每个簇的质心作为新中心。
     • 重复直到收敛或达到最大迭代次数。
   • 时间复杂度：每迭代$\mathcal{O}(nkd)$，通常迭代次数较少。

5. 理论保证：

   • K-MC²的量化误差满足：
     $$\mathbb{E}[{\varphi }_{{\text{K-MC}}^2}]\le \mathbb{E}[{\varphi }_{\text{k-means++}}]+2ϵ\beta {\varphi }_{OPT}^k(\mathcal{X})$$
     设置$ϵ=\mathcal{O}(1/\beta )$，则：
     $$\mathbb{E}[{\varphi }_{{\text{K-MC}}^2}]\le \mathcal{O}(\log k){\varphi }_{OPT}^k(\mathcal{X})$$
     保留$k$-means++的$\mathcal{O}(\log k)$近似保证。

---

主要的贡献点

1. 亚线性复杂度的初始化算法：
   • 提出K-MC²算法，通过MCMC近似$D^2$采样，将初始化复杂度从$\Theta (nkd)$降至$\mathcal{O}(k^{3}d{\log }^{2}n\log k)$，实现亚线性时间初始化。
2. 理论保证：
   • 证明K-MC²在总变差距离下收敛到$k$-means++，并在非病态数据集（满足A1和A2假设）下保留$\mathcal{O}(\log k)$近似保证。
3. 数据集分布假设：
   • 分析${\gamma }^{\prime }$的上界，证明在常见分布（如高斯、指数、拉普拉斯）下，$\alpha =\mathcal{O}({\log }^{2}n)$，$\beta =\mathcal{O}(k)$，确保算法高效性。
4. 广泛的实验验证：
   • 在六个真实大规模数据集（USGS、CSN、KDD、BIGX、WEB、SONG）上验证K-MC²的性能，显示其在量化误差和运行时间上的竞争力。
5. 适用性扩展：
   • K-MC²适用于在线聚类、mini-batch $k$-means和coreset构造等场景，提升了$k$-means++在大规模数据处理中的实用性。

---

实验结果

实验在六个数据集上进行，数据集信息如表1所示：

数据集	点数 ($n$)	维度 ($d$)	$\alpha$	$\beta$ ($k=200$)
CSN	80,000	17	546.27	3.04
KDD	145,751	74	1267.65	1.81
USGS	59,209	3	2.68	51.67
WEB	45,811,883	5	2.33	57.09
BIGX	1,162,000	57	5.22	1.17
SONG	515,345	90	525.67	1.23

实验设置：

• 对于USGS、CSN、KDD，设置$k=200$，在全数据集上训练，评估量化误差和距离计算次数。
• 对于BIGX、WEB、SONG，设置$k=2000$，保留250,000点作为测试集，评估训练误差和泛化误差。
• 比较方法：K-MC²（链长 m ∈ { 1 , 2 , 5 , 10 , 20 , 50 , 100 , 150 , 200 } m \in \{1, 2, 5, 10, 20, 50, 100, 150, 200\} m∈{1,2,5,10,20,50,100,150,200}）、$k$-means++、RANDOM、HEURISTIC（在子集上运行$k$-means++，子集大小 s ∈ { 100 , 200 , … , 20 , 000 } s \in \{100, 200, \ldots, 20,000\} s∈{100,200,…,20,000}）、$k$-means||（$r=5$轮，过采样因子 l ∈ { 0.02 k , 0.05 k , … , 2 k } l \in \{0.02k, 0.05k, \ldots, 2k\} l∈{0.02k,0.05k,…,2k}）。
• 所有方法重复运行多次，取平均量化误差，计算95%置信区间。

主要结果（表2）：

1. K-MC² vs $k$-means++：
   • K-MC²随链长$m$增加快速接近$k$-means++的量化误差：
     • $m=20$时，相对误差（相对于$k$-means++）在USGS为2.63%，KDD为32.62%，SONG为0.75%。
     • $m=200$时，相对误差降至USGS 0.33%，KDD 1.00%，SONG 0.02%，接近$k$-means++。
   • 运行时间显著降低：
     • 在USGS上，K-MC² ($m=200$)比$k$-means++快5倍。
     • 在WEB上（$n=45.8M$），K-MC²快275.1倍。
2. K-MC² vs RANDOM：
   • RANDOM的量化误差远高于$k$-means++，如CSN上高334.50%，WEB上高105.54%。
   • K-MC²在所有数据集上显著优于RANDOM，即使$m$较小。
3. K-MC² vs HEURISTIC：
   • HEURISTIC的性能随子集大小$s$增加而改善，但总体不如K-MC²：
     • 在CSN上，HEURISTIC ($s=20,000$)相对误差为54.72%，而K-MC² ($m=200$)为6.53%。
     • 图1显示K-MC²的量化误差收敛更快，优于HEURISTIC。
4. K-MC² vs $k$-means||：
   • $k$-means||的性能随过采样因子$l$增加而接近$k$-means++，但仍逊于K-MC²：
     • 在BIGX上，K-MC² ($m=200$)相对误差0.03%，而$k$-means|| ($l=2k$)为0.33%。
     • 图2显示K-MC²在更少的距离计算下达到竞争性解。
5. 泛化误差：
   • 在BIGX、WEB、SONG上，K-MC²的测试集误差与训练集误差一致，表明其初始化具有良好的泛化能力。

总结：

• K-MC²在所有数据集上以较少的计算成本（亚线性复杂度）获得接近$k$-means++的量化误差，尤其在超大规模数据集（如WEB）上效率优势显著。
• 其性能优于RANDOM、HEURISTIC和$k$-means||，验证了理论保证的实用性。

---

算法的实现过程

K-MC²算法的核心是通过MCMC近似$k$-means++的$D^2$采样，以下是详细实现过程，结合伪代码和说明：

Algorithm K-MC²
输入: 数据集 𝒳 = {x_1, x_2, ..., x_n}，中心数 k，链长 m
输出: 初始中心集合 C_k = {c_1, c_2, ..., c_k}1. 初始化:均匀随机选择 c_1 ∈ 𝒳C_1 ← {c_1}2. 循环选择 k-1 个中心:for i = 2 to k do// 初始化马尔可夫链均匀随机选择 x ∈ 𝒳计算 d_x ← d(x, C_{i-1})^2// 运行 m 步 Metropolis-Hastingsfor j = 2 to m do均匀随机选择候选点 y ∈ 𝒳计算 d_y ← d(y, C_{i-1})^2计算接受概率 π ← min(1, d_y / d_x)生成 u ← Unif(0, 1)if u < π thenx ← y, d_x ← d_yC_i ← C_{i-1} ∪ {x}end3. 返回 C_k

实现细节：

1. 均匀采样：
   • 使用随机数生成器从$\mathcal{X}$中均匀选择点，确保初始中心$c_1$和马尔可夫链起点$x$的随机性。
2. 距离计算：
   • 计算点$x$到中心集$C_{i-1}$的距离$\mathrm{d}(x,C_{i-1})={\min }_{c\in C_{i-1}}\Vert x-c{\Vert }_2$，平方后得$d_x$。
   • 每步仅需计算一个点到$C_{i-1}$的距离，复杂度为$\mathcal{O}(kd)$，避免$k$-means++中对所有点的全扫描。
3. Metropolis-Hastings采样：
   • 使用独立均匀提议分布，目标分布为$D^2$采样分布$p(x)$。
   • 接受概率$\pi =\min (1,\frac{d_y}{d_x})$无需计算归一化常数${\sum }_{x^{\prime }\in \mathcal{X}}\mathrm{d}(x^{\prime },C_{i-1}{)}^2$，简化计算。
   • 每次迭代生成均匀随机数$u\in [0,1]$，若$u<\pi$，接受新点$y$。
4. 链长选择：
   • 链长$m$根据${\gamma }^{\prime }$和期望总变差距离$ϵ$确定。在实践中，实验测试了 m ∈ { 1 , 2 , … , 200 } m \in \{1, 2, \ldots, 200\} m∈{1,2,…,200}，发现$m=100$至$200$已接近$k$-means++性能。
   • 理论上，$m=\mathcal{O}(k{\log }^{2}n\log k)$保证亚线性复杂度。
5. 复杂度分析：
   • 每步MCMC计算$\mathcal{O}(kd)$，$m$步共$\mathcal{O}(mkd)$，$k-1$次迭代共$\mathcal{O}(m{k}^{2}d)$。
   • 总复杂度为$\mathcal{O}(k^{3}d{\log }^{2}n\log k)$，不依赖$n$，适合大规模数据集。
6. 后续$k$-means：
   • K-MC²输出初始中心$C_k$后，可直接输入Lloyd算法，执行标准$k$-means迭代。
   • 实验未直接评估Lloyd迭代，但理论保证K-MC²初始化的质量接近$k$-means++，后续迭代收敛性相似。

技术考虑：

• 随机种子：实验中多次运行不同随机种子以平均量化误差，确保结果稳健。
• 并行化：K-MC²的MCMC采样本质上是顺序的，但每轮的$k-1$次采样可部分并行（如在多核CPU上并行运行不同链）。
• 数据存储：对于大规模数据集（如WEB，$n=45.8M$），需高效存储和访问数据点，可能使用分布式文件系统或内存映射。
• 数值稳定性：距离计算使用浮点运算，确保$d_y/d_x$不会因数值溢出导致错误。

与$k$-means++的对比：

• $k$-means++的$D^2$采样需要计算所有点的概率$p(x)$，每次迭代扫描全数据集，复杂度$\Theta (nkd)$。
• K-MC²通过MCMC采样，仅需计算$m$个点的距离，采样分布随$m$增加接近$p(x)$，大幅降低计算量。

---

总结

K-MC²算法通过MCMC近似$D^2$采样，解决了$k$-means++在大规模数据集上的初始化效率问题，在亚线性时间复杂度下保留了$\mathcal{O}(\log k)$近似保证。其理论分析基于温和的分布假设，实验结果在多个真实数据集上验证了其高效性和竞争力。算法实现简单，适合多种聚类场景，为大规模数据聚类提供了实用工具。