SIAM-2010《Making $k$-means even faster》

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

核心思想

论文的核心思想是通过改进Elkan的$k$-means加速算法，提出一种更快、更高效的精确$k$-means聚类算法，特别适用于低维和中维数据（维度最高约50）。该算法利用三角不等式和距离界（distance bounds）避免不必要的点-中心距离计算，尤其通过引入一个新颖的单下界（lower bound）来显著减少最内层循环的执行频率（实验中超过80%），从而加速Lloyd’s $k$-means算法，同时保持结果的精确性。算法在时间复杂度和内存开销上均优于其他加速方法（如$k$-d树和Elkan算法），且实现简单，适合并行化和大规模数据集处理。

目标函数

$k$-means算法的目标函数是最小化点到其分配中心的平方距离和，即$k$-means准则。数学表达为：

$$J=\sum _{i=1}^n\Vert x(i)-c(a(i)){\Vert }^2$$

其中：

• $x(i)$：数据集中的第$i$个数据点， i ∈ { 1 , … , n } i \in \{1, \dots, n\} i∈{1,…,n}。
• $c(j)$：第$j$个聚类中心， j ∈ { 1 , … , k } j \in \{1, \dots, k\} j∈{1,…,k}。
• $a(i)$：数据点$x(i)$分配到的聚类中心的索引。
• $n$：数据点总数。
• $k$：聚类中心数。
• $\Vert \cdot \Vert$：欧几里得距离。

目标是通过迭代优化，使得每个数据点分配到距离最近的中心，同时更新中心位置以最小化$J$。

目标函数的优化过程

论文提出的算法基于Lloyd’s算法，通过引入距离界和三角不等式优化目标函数的计算过程。优化过程分为以下步骤：

1. 初始化：

   • 使用$k$-means++算法选择初始中心$c(j)$，以提高收敛质量。
   • 为每个数据点$x(i)$初始化：
     • 上界$u(i)$：$x(i)$到其分配中心$c(a(i))$的距离上界。
     • 下界$l(i)$：$x(i)$到第二近中心（非$c(a(i))$）的距离下界。
     • 分配索引$a(i)$：记录$x(i)$所属的中心。
   • 初始化中心相关结构，如中心移动距离$p(j)$和中心间最小距离$s(j)$。

2. 迭代过程：

   • 更新中心间距离：计算每个中心$c(j)$到其他最近中心的距离$s(j)$，时间复杂度为$O(d{k}^2)$，其中$d$是数据维度。
   • 数据点分配：
     • 对于每个数据点$x(i)$，检查是否可以跳过最内层循环（遍历所有中心的距离计算）：
       • 如果$u(i)\le s(a(i))/2$（基于三角不等式）或$u(i)\le l(i)$（新颖下界），则无需重新计算距离，保持当前分配$a(i)$。
       • 否则，重新计算$x(i)$到其分配中心$c(a(i))$的精确距离，更新$u(i)$，再次检查上述条件。
       • 如果仍需计算，则调用POINT-ALL-CTRS子程序，遍历所有中心，更新$a(i)$、$u(i)$和$l(i)$，时间复杂度为$O(kd)$。
   • 更新中心：
     • 使用增量更新（delta updates）方法，根据分配变化调整中心向量和$c^{\prime }(j)$（分配点的向量和）和计数$q(j)$，避免重新计算所有点的均值。
     • 更新中心$c(j)=c^{\prime }(j)/q(j)$。
   • 更新距离界：
     • 上界$u(i)$增加$c(a(i))$的移动距离$p(a(i))$。
     • 下界$l(i)$减少所有中心的最大移动距离${\max }_{j}p(j)$。
     • 保证$u(i)\ge \Vert x(i)-c(a(i))\Vert$，$l(i)\le {\min }_{j\ne a(i)}\Vert x(i)-c(j)\Vert$。

3. 收敛判断：

   • 当中心不再移动或分配不再变化时，算法收敛，输出最终聚类结果。

通过上述优化，算法显著减少了最内层循环的执行频率（实验中高达80%以上），从而降低总计算时间。

主要贡献点

1. 新颖的下界设计：

   • 提出每个数据点仅维护一个下界$l(i)$，表示到第二近中心的距离下界，而非Elkan算法中的$k$个下界（每个点到每个中心的下界）。这减少了内存开销（从$O(nk)$降至$O(n)$）和更新开销，同时仍能有效跳过最内层循环（实验中跳跃率达80%以上）。

2. 高效的循环跳跃：

   • 通过结合$u(i)\le s(a(i))/2$和$u(i)\le l(i)$两个条件，算法在低维和中维数据上跳跃最内层循环的频率远高于Elkan算法，尤其在均匀分布数据上表现优异。

3. 低内存开销：

   • 内存开销为$O(3n+2k)$，远低于Elkan算法的$O(nk+k^2)$和$k$-d树的$O(nd)$，适合大规模数据集和内存受限环境。

4. 实现简单且易并行化：

   • 算法无需复杂的数据结构（如$k$-d树），实现简单，迭代式设计避免递归调用，易于并行化处理大规模数据集。

5. 实验验证：

   • 在低维（2-32维）和中维（50-56维）数据集上，算法显著优于Lloyd’s算法（平均快6.1倍）、$k$-d树（快2.8倍）和Elkan算法（快2.4倍），尤其在均匀分布数据上表现突出。

6. 互补性与扩展性：

   • 算法与Elkan算法互补，前者适合低维和中维数据，后者适合高维数据。论文提出未来可结合两者，动态选择下界数量以适应不同维度。

实验结果

1. 数据集：

   • 合成数据集：均匀分布的超立方体数据，维度为2、8、32、128，数据点数$n=1250000$。
   • 应用数据集：birch（$n=100000$，$d=2$）、covtype（$n=150000$，$d=54$）、kddcup（$n=95412$，$d=56$）、mnist50（$n=60000$，$d=50$）。

2. 时间性能（表2）：

   • 在2-32维数据上，Hamerly算法总时间和每迭代时间均优于其他算法。例如，在8维均匀随机数据集（$k=500$），Hamerly算法总时间为199.6秒，远低于Lloyd’s（13854.4秒）、$k$-d树（373.5秒）和Elkan（187.9秒）。
   • 在50-56维数据上，Hamerly算法在小$k$时表现最佳，大$k$时与Elkan算法相当。
   • 平均每迭代时间，Hamerly算法比Lloyd’s快6.1倍，比$k$-d树快2.8倍，比Elkan快2.4倍。

3. 最内层循环跳跃率（表3）：

   • Hamerly算法跳跃率远高于Elkan算法。例如，在8维均匀随机数据上，Hamerly跳跃率为0.88，Elkan仅为0.01；在应用数据集上，Hamerly跳跃率在0.82-0.94之间，Elkan为0.18-0.52。
   • 即使在128维数据上，Hamerly跳跃率仍达0.83，显示其鲁棒性。

4. 内存使用（表4）：

   • Hamerly算法内存使用量最低，接近Lloyd’s算法。例如，在$n=1250000$，$d=2$，$k=500$时，Hamerly使用14.7 MB，远低于Elkan的1205.2 MB和$k$-d树的32.1 MB。

5. 消融研究（Lesion Study）：

   • 移除下界$l(i)$（No-Lower）导致性能显著下降，证明下界是加速的关键。
   • 移除$s$（No-S）在高维数据上表现较好，因避免了$O(k^2)$中心间距离计算。
   • 移除上界$u(i)$（No-Upper）在低维数据上效果较好，但高维下因精确距离计算成本高而变慢。
   • 移除增量更新（No-Delta）在大$k$时表现较好，因频繁分配变化使增量更新开销增加。

算法实现过程详细解释

以下详细描述Hamerly算法的实现过程，基于论文中的Algorithm 1-5。

1. 数据结构

• 中心相关：
  • $c(j)$：第$j$个中心，维度$d$的向量。
  • $c^{\prime }(j)$：分配到第$j$个中心的点向量和。
  • $q(j)$：分配到第$j$个中心的点数。
  • $p(j)$：中心$c(j)$最近一次移动的距离。
  • $s(j)$：$c(j)$到其他最近中心的距离。
• 数据点相关：
  • $x(i)$：第$i$个数据点，维度$d$。
  • $a(i)$：$x(i)$分配的中心索引。
  • $u(i)$：$x(i)$到$c(a(i))$的距离上界。
  • $l(i)$：$x(i)$到第二近中心的距离下界。

2. 算法伪代码（Algorithm 1: K-means）

Algorithm 1 K-means(dataset x, initial centers c)
  INITIALIZE(c, x, q, c', u, l, a)
  while not converged do
    for j=1 to |c| do
      s(j) = min_{j'≠j} ||c(j) - c(j')||  // 更新中心间最小距离
    for i=1 to |x| do
      if u(i) > s(a(i))/2 and u(i) > l(i) then  // 检查是否可跳跃
        u(i) = ||x(i) - c(a(i))||  // 更新上界
        if u(i) > s(a(i))/2 and u(i) > l(i) then
          POINT-ALL-CTRS(x(i), c, u(i), l(i), a(i))  // 遍历所有中心
    MOVE-CENTERS(c, c', q, p)  // 更新中心
    UPDATE-BOUNDS(u, l, p)  // 更新界

3. 子程序

• INITIALIZE（Algorithm 2）：

  • 使用$k$-means++初始化中心$c$。
  • 计算初始$u(i)$、$l(i)$、$a(i)$，通过遍历所有中心确定每个点的最近和次近中心。
  • 初始化$c^{\prime }(j)$、$q(j)$、$p(j)$、$s(j)$。

• POINT-ALL-CTRS（Algorithm 3）：

  • 输入：点$x(i)$，中心$c$，当前$u(i)$、$l(i)$、$a(i)$。
  • 计算$x(i)$到每个中心的距离，更新：
    • $a(i)$：最近中心的索引。
    • $u(i)$：到最近中心的距离。
    • $l(i)$：到次近中心的距离。
  • 时间复杂度：$O(kd)$。

• MOVE-CENTERS（Algorithm 4）：

  • 如果分配变化，更新$c^{\prime }(j)$和$q(j)$：
    • 从旧中心减去$x(i)$的影响，添加到新中心。
  • 计算新中心：$c(j)=c^{\prime }(j)/q(j)$。
  • 更新$p(j)$为$c(j)$的移动距离。

• UPDATE-BOUNDS（Algorithm 5）：

  • 更新上界：$u(i)=u(i)+p(a(i))$。
  • 更新下界：$l(i)=l(i)-{\max }_{j}p(j)$。

4. 关键优化点

• 跳跃最内层循环：
  • 使用$u(i)\le s(a(i))/2$（三角不等式）和$u(i)\le l(i)$（新下界）判断是否需要计算所有点-中心距离。
  • 如果条件满足，直接保留当前分配，跳过$O(kd)$计算。
• 增量更新：
  • 仅更新分配变化的点的$c^{\prime }(j)$和$q(j)$，避免重新计算所有点的均值。
• 低内存开销：
  • 仅存储一个下界$l(i)$，内存开销为$O(n)$，而Elkan算法需要$O(nk)$。
• 避免递归：
  • 迭代式设计，相比$k$-d树的递归调用更高效。

5. 实现细节

• 避免平方根运算：仅比较平方距离以确定最近中心。
• 缓存距离：避免重复计算点-中心距离。
• 并行化：数据点相关结构（$x$、$u$、$l$、$a$）可分区，中心相关结构（$c$、$c^{\prime }$、$q$、$p$、$s$）复制并同步，适合多核处理。
• 大规模数据集：支持流式处理，数据点和界可从磁盘顺序读取。

总结

Hamerly的算法通过引入单一新颖下界和优化距离界更新，显著加速了$k$-means算法，特别在低维和中维数据上表现优异。其简单性、低内存需求和并行化潜力使其适用于大规模数据集和内存受限场景。实验结果验证了其高效性，跳跃最内层循环的频率高达80%以上，内存使用量远低于其他加速算法。未来可通过结合Elkan算法或引入树结构进一步优化，适应更广泛的应用场景。