K-Means 聚类算法的 CUDA 实现

K-Means 聚类算法的 CUDA 实现

本文将介绍 K - means 聚类算法的 CUDA C 实现。K - means 聚类是一种硬聚类算法，这意味着每个数据点都被分配到一个聚类中（而不是以不同概率分配到多个聚类）。该算法从随机的聚类分配开始，并在以下两个步骤之间迭代：

1. 根据最近的质心（通过某种距离度量）将数据点分配到聚类中。
2. 根据上一步的新聚类分配更新质心。

最终聚类分配会收敛，得到最终结果。在这个 CUDA 实现中，这两个步骤中的每一个都将并行执行。此实现是针对一维情况，其中每个数据点是一个标量，但它可以很容易地扩展到多维情况。

将数据点分配到聚类

第一步是将数据点分配到它们最近的质心。这一步不难并行化，因为每个数据点的距离计算可以独立进行。

__global__ void kMeansClusterAssignment(float *d_datapoints, int *d_clus)
{
    // 获取此数据点的索引
    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
    // 边界检查
    if (idx >= N) return;
    // 找到此数据点最近的质心
    float min_dist = INFINITY;
    int closest_centroid = 0;
    for (int c = 0; c < K; ++c)
    {
        float dist = distance(d_datapoints[idx], d_centroids[c]);//计算距离
        if (dist < min_dist)
        {
            min_dist = dist;
            closest_centroid = c;
        }
    }
    // 为这个数据点/线程分配最近的聚类 ID
    d_clust_assn[idx] = closest_centroid;
}

上述代码实际上与串行情况非常相似。这里线程索引 idx 对应于数据点的索引。其余代码相当简单明了。计算数据点 d_datapoints[idx] 与每个质心 d_centroids[c] 之间的距离，然后将最近的质心分配给该数据点。这段代码的一个缺点是质心和数据点是从全局内存中读取的，速度有点慢。

这里使用的距离函数就是一维情况下的欧几里得距离。

__device__ float distance(float x1, float x2)
{
    return sqrt((x2 - x1) * (x2 - x1));
}

更新质心

下一步是根据上一步计算出的聚类分配重新计算质心。这一步并行化要棘手得多，因为质心计算依赖于其聚类中的所有其他数据点。然而，依赖分布式数据集计算单个输出值的操作仍然可以并行化，这被称为规约（reductions）。下面的图展示了求和规约。

一般的过程是对输入数组进行分区，并行地对每个分区求和，然后合并分区，并重复这个过程，直到所有分区都被合并，最终得到求和结果。这种并行化使得复杂度为对数级，而不是像串行情况那样为线性。

质心重新计算也可以看作是一种规约操作，代码如下所示。

__global__ void kMeansCentroidUpdate(float *d_datapoints, int *d_clust_a)
{
    // 获取网格级线程的索引
    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
    // 边界检查
    if (idx >= N) return;
    // 获取块级线程的索引
    const int s_idx = threadIdx.x;
    // 将数据点和相应的聚类分配放入共享内存
    __shared__ float s_datapoints[TPB];//TPB为每个线程块中的线程数量，即数据点的数量等于线程块中的线程数量
    s_datapoints[s_idx] = d_datapoints[idx];//把全局内存中数据点复制到共享内存中，共享内存速度更快
    __shared__ int s_clust_assn[TPB];//聚类分配
    s_clust_assn[s_idx] = d_clust_assn[idx];//将全局内存中的聚类复制到共享内存中
    __syncthreads();
    // 每个块中索引为 0 的线程对所有数据进行求和
    if (s_idx == 0)
    {	
        float b_clust_datapoint_sums[K] = { 0 };//临时存储每个聚类中数据点总和的数组
        int b_clust_sizes[K] = { 0 };//每个聚类中数据点的数量
        for (int j = 0; j < blockDim.x; ++j)
        {
            int clust_id = s_clust_assn[j];//获取s_datapoints[j]所属的聚类编号
            b_clust_datapoint_sums[clust_id] += s_datapoints[j];//将数据点s_datapoints[j]的值累加
            b_clust_sizes[clust_id] += 1;//统计聚类的数据点数量
        }
        // 现在将求和结果添加到全局质心和全局聚类大小数组
        for (int z = 0; z < K; ++z)
        {
        	//d_centroids 用于存储所有聚类的质心信息（在当前步骤中，存储的是所有线程块处理后每个聚类的数据点总和）。
        	//累加数据点总和到全局质心数组
            atomicAdd(&d_centroids[z], b_clust_datapoint_sums[z]);
            //累加数据点数量到全局聚类大小数组
            atomicAdd(&d_clust_sizes[z], b_clust_sizes[z]);
        }
    }
    __syncthreads();
    // 目前质心只是求和结果，所以除以大小以得到实际质心
    if (idx < K)
    {
        d_centroids[idx] = d_centroids[idx] / d_clust_sizes[idx];
    }
}

在上述代码中，s_idx 指的是块级线程（即，如果使用多个块，多个线程将具有相同的 s_idx）。接下来的代码：

__shared__ float s_datapoints[TPB];
s_datapoints[s_idx] = d_datapoints[idx];
__shared__ int s_clust_assn[TPB];
s_clust_assn[s_idx] = d_clust_assn[idx];

是规约的分区步骤，其中全局数据点（d_datapoints）和聚类分配（d_clust_assn）变量被读入线程 idx 所在块的共享内存中。

接下来的代码仅在每个块的第一个线程中执行（即 if(s_idx==0)...）。这个线程创建了两个本地变量——b_clust_datapoint_sums，它是一个长度为 K 的数组，用于保存分配到每个聚类的数据点的总和；b_clust_sizes，它也是一个长度为 K 的数组，用于保存每个聚类的大小。然后线程遍历其块的共享内存中的数据点和聚类分配，并将结果累加到这些本地变量中。

下一步是将这些中间结果合并到全局变量 d_centroids 和 d_clust_sizes 中以获得最终更新结果。这是通过 atomicAdd 函数完成的，该函数确保当多个线程写入同一位置时不会发生竞态条件。但是，我们不希望质心是求和结果，所以最后一步是将总和除以聚类大小，这可以在全局线程级别并行完成。

心得

K-means算法的cuda实现分为两步

1. 根据最近的质心（通过某种距离度量）将数据点分配到聚类中。
2. 根据上一步的新聚类分配更新质心。

质心：质心的位置，使每个数据点到它所属聚类的质心的距离之和尽可能小。

第一步并行计算各个数据点与各个质心得距离，将数据点分配给距离最近得质心。缺点是数据点与质心从全局内存读取，速度较慢。
在 CUDA 编程中，多个线程块可能会同时对全局数组进行写操作。如果不使用原子操作，多个线程同时修改同一个全局数组元素时，可能会出现数据竞争问题，导致最终结果错误。例如，两个线程同时读取 d_centroids[z] 的值，然后各自进行累加操作，再写回该位置，这样就会丢失其中一个线程的累加结果。
atomicAdd 函数确保了在同一时刻只有一个线程能够对指定的内存位置进行修改，避免了数据竞争，保证了全局数组更新的线程安全性

求实际质心过程中应该确保聚类的大小不为0