从零实现KMeans聚类：用Python代码揭开算法的神秘面纱

当我在第一次接触机器学习时，那些复杂的数学公式和抽象概念总让我望而却步。直到有一天，我决定抛开厚重的教科书，直接在Jupyter Notebook中一行行敲出KMeans算法的实现代码，才真正理解了聚类的精髓。本文将带你体验这种"做中学"的乐趣，无需死记硬背公式，通过可视化代码理解KMeans的每个关键步骤。

1. 准备工作：搭建实验环境

在开始之前，我们需要准备一个交互式的Python环境。推荐使用Jupyter Notebook，它能让我们实时看到代码执行结果和可视化效果。以下是需要安装的核心库：

pip install numpy matplotlib ipython

这些库将分别提供：

• NumPy ：高效的数值计算基础
• Matplotlib ：数据可视化工具
• IPython ：交互式编程环境

提示：如果你使用Anaconda发行版，这些库通常已经预装好了

让我们先导入必要的模块并设置随机种子，确保每次运行结果一致：

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)  # 宇宙的终极答案

2. 创建模拟数据集：理解聚类的输入

好的数据集是理解算法的第一步。我们将生成三个明显分组的二维数据点，这样在可视化时可以直观地看到聚类效果。

# 定义每个簇的参数
cluster_1 = np.random.normal(loc=[0,0], scale=0.5, size=(50,2))
cluster_2 = np.random.normal(loc=[5,5], scale=0.8, size=(70,2))
cluster_3 = np.random.normal(loc=[-4,5], scale=0.3, size=(30,2))

# 合并所有数据点
X = np.vstack([cluster_1, cluster_2, cluster_3])

# 可视化原始数据
plt.scatter(X[:,0], X[:,1], alpha=0.6)
plt.title("原始未标记数据分布")
plt.xlabel("特征1")
plt.ylabel("特征2")
plt.show()

这段代码会产生约150个数据点，分布在三个不同的区域。在真实场景中，你的数据可能来自：

• 客户购买行为数据
• 传感器采集的测量值
• 图像的特征向量

3. 手动实现KMeans核心逻辑

现在来到最激动人心的部分——从零实现KMeans算法。我们将分步骤构建这个经典的无监督学习算法。

3.1 初始化质心：算法的起点

KMeans的第一步是随机选择K个点作为初始质心。质心的选择会影响算法的收敛速度，但不会改变最终结果（在全局最优情况下）。

def initialize_centroids(X, k):
    """随机选择k个数据点作为初始质心"""
    indices = np.random.choice(X.shape[0], k, replace=False)
    return X[indices]

3.2 分配数据点到最近质心

对于每个数据点，我们计算它与所有质心的距离，并将其分配到最近的质心所在的簇。

def assign_clusters(X, centroids):
    """将每个点分配到最近的质心"""
    distances = np.linalg.norm(X[:, np.newaxis] - centroids, axis=2)
    return np.argmin(distances, axis=1)

这里使用了欧式距离，但你可以尝试其他距离度量如曼哈顿距离或余弦相似度。

3.3 更新质心位置

根据当前簇分配，重新计算每个簇的质心（即该簇所有点的平均值）。

def update_centroids(X, labels, k):
    """计算每个簇的新质心"""
    new_centroids = np.array([X[labels == i].mean(axis=0) for i in range(k)])
    return new_centroids

3.4 组合成完整算法

现在我们将这些步骤组合起来，形成完整的KMeans算法：

def kmeans(X, k, max_iter=100):
    """完整的KMeans实现"""
    # 1. 初始化质心
    centroids = initialize_centroids(X, k)
    
    for _ in range(max_iter):
        # 2. 分配簇
        labels = assign_clusters(X, centroids)
        
        # 3. 更新质心
        new_centroids = update_centroids(X, labels, k)
        
        # 检查收敛
        if np.all(centroids == new_centroids):
            break
            
        centroids = new_centroids
    
    return labels, centroids

4. 可视化聚类过程：眼见为实

为了真正理解KMeans的工作原理，让我们创建一个动态可视化，展示算法每一步的变化。

def plot_kmeans_steps(X, k=3, max_iter=10):
    """可视化KMeans的迭代过程"""
    centroids = initialize_centroids(X, k)
    
    plt.figure(figsize=(15, 10))
    
    for i in range(max_iter):
        # 分配簇
        labels = assign_clusters(X, centroids)
        
        # 绘制当前状态
        plt.subplot(2, 5, i+1)
        for cluster in range(k):
            cluster_points = X[labels == cluster]
            plt.scatter(cluster_points[:,0], cluster_points[:,1], alpha=0.6)
            plt.scatter(centroids[cluster,0], centroids[cluster,1], 
                        marker='*', s=200, c='black')
        plt.title(f'迭代 {i+1}')
        
        # 更新质心
        new_centroids = update_centroids(X, labels, k)
        
        # 检查收敛
        if np.all(centroids == new_centroids):
            break
            
        centroids = new_centroids
    
    plt.tight_layout()
    plt.show()

调用这个函数，你将看到质心如何一步步移动到簇的中心，以及数据点如何被重新分配。

5. 与Scikit-learn实现对比验证

为了验证我们的实现是否正确，让我们将其与Scikit-learn中的KMeans进行比较。

from sklearn.cluster import KMeans

# 使用我们的实现
our_labels, our_centroids = kmeans(X, k=3)

# 使用Scikit-learn
sklearn_kmeans = KMeans(n_clusters=3, random_state=42)
sklearn_labels = sklearn_kmeans.fit_predict(X)
sklearn_centroids = sklearn_kmeans.cluster_centers_

# 比较结果
print("我们的质心:\n", our_centroids)
print("Scikit-learn的质心:\n", sklearn_centroids)

如果实现正确，两者的质心位置应该非常接近（可能因初始化不同有微小差异）。

6. 算法优化与实用技巧

在实际应用中，我们还需要考虑一些优化和技巧：

6.1 处理空簇问题

有时一个簇可能会失去所有成员，导致计算新质心时出错。我们可以通过以下方式处理：

def update_centroids_safe(X, labels, k):
    """处理可能出现的空簇"""
    new_centroids = []
    for i in range(k):
        cluster_points = X[labels == i]
        if len(cluster_points) > 0:
            new_centroids.append(cluster_points.mean(axis=0))
        else:
            # 如果簇为空，重新随机初始化一个质心
            new_centroids.append(X[np.random.randint(X.shape[0])])
    return np.array(new_centroids)

6.2 多次初始化避免局部最优

KMeans对初始质心敏感，常见的解决方案是运行算法多次，选择最佳结果：

def run_kmeans_multiple_times(X, k, n_init=10):
    """多次运行KMeans，选择最佳结果"""
    best_inertia = float('inf')
    best_labels = None
    best_centroids = None
    
    for _ in range(n_init):
        labels, centroids = kmeans(X, k)
        inertia = np.sum([np.sum((X[labels == i] - centroids[i])**2) 
                         for i in range(k)])
        
        if inertia < best_inertia:
            best_inertia = inertia
            best_labels = labels
            best_centroids = centroids
    
    return best_labels, best_centroids

6.3 选择最佳K值：肘部法则

确定合适的簇数量K是使用KMeans的关键挑战。肘部法则通过观察不同K值下的惯性变化来帮助选择：

def find_optimal_k(X, max_k=10):
    """使用肘部法则寻找最佳K值"""
    inertias = []
    for k in range(1, max_k+1):
        labels, centroids = run_kmeans_multiple_times(X, k)
        inertia = np.sum([np.sum((X[labels == i] - centroids[i])**2) 
                         for i in range(k)])
        inertias.append(inertia)
    
    plt.plot(range(1, max_k+1), inertias, marker='o')
    plt.xlabel('K值')
    plt.ylabel('惯性(Inertia)')
    plt.title('肘部法则')
    plt.show()

7. 真实世界应用案例

理解算法后，让我们看几个实际应用场景：

7.1 客户细分

假设我们有一组客户的购买行为数据，包含两个特征：年消费额和购买频率。我们可以使用KMeans将客户分成不同的群体：

# 模拟客户数据
np.random.seed(1)
high_value = np.random.normal(loc=[10000, 15], scale=[2000, 3], size=(30,2))
medium_value = np.random.normal(loc=[5000, 8], scale=[1000, 2], size=(50,2))
low_value = np.random.normal(loc=[2000, 3], scale=[500, 1], size=(70,2))
customers = np.vstack([high_value, medium_value, low_value])

# 应用KMeans
labels, centroids = kmeans(customers, k=3)

# 可视化结果
plt.scatter(customers[:,0], customers[:,1], c=labels, alpha=0.6)
plt.scatter(centroids[:,0], centroids[:,1], marker='*', s=200, c='black')
plt.xlabel('年消费额')
plt.ylabel('年购买次数')
plt.title('客户细分')
plt.show()

7.2 图像压缩

KMeans还可以用于图像颜色压缩，通过减少颜色数量来减小文件大小：

from sklearn.datasets import load_sample_image

# 加载示例图片
china = load_sample_image("china.jpg")
image_array = china / 255.0  # 归一化到0-1

# 将图像数据重塑为二维数组
w, h, d = image_array.shape
image_2d = image_array.reshape(w * h, d)

# 使用KMeans找到16种主要颜色
k = 16
labels, centroids = kmeans(image_2d, k)

# 用质心颜色重建图像
compressed_image = centroids[labels].reshape(w, h, d)

# 显示原始和压缩后的图像
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.imshow(image_array)
plt.title('原始图像')
plt.subplot(1, 2, 2)
plt.imshow(compressed_image)
plt.title(f'压缩图像 ({k}色)')
plt.show()

在实现过程中，我发现KMeans虽然简单，但对初始质心非常敏感。多次运行算法并选择最佳结果可以显著提高稳定性。此外，在处理大型数据集时，可以考虑使用MiniBatchKMeans来提高计算效率。