从零实现K-means聚类：Python3实战指南与可视化解析

初识K-means：聚类算法的核心思想

当我们需要对无标签数据进行分组时，K-means算法就像一位经验丰富的分类专家。想象你面前摆着数百种未知的植物标本，K-means能自动将它们分成K个具有相似特征的群体，而无需预先知道任何分类标准。这种无监督学习算法通过迭代寻找数据中的自然分组，广泛应用于客户细分、图像压缩、异常检测等领域。

K-means的魅力在于其数学简洁性与实践有效性的完美结合。算法核心仅需三步：

1. 随机初始化K个中心点（质心）
2. 将每个数据点分配到最近的中心点形成簇
3. 重新计算每个簇的新中心点

# 基础K-means伪代码展示
def k_means(data, k):
    # 1. 随机选择k个初始中心点
    centroids = initialize_centroids(data, k)
    
    while not converged:
        # 2. 分配每个点到最近的中心点
        clusters = assign_points_to_clusters(data, centroids)
        
        # 3. 重新计算中心点
        new_centroids = calculate_new_centroids(clusters)
        
        # 判断是否收敛
        converged = check_convergence(centroids, new_centroids)
        centroids = new_centroids
    
    return clusters, centroids

环境准备与数据加载

1.1 工具库配置

现代Python生态为机器学习提供了强大支持。我们将使用以下核心库：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score

关键工具对比：

工具库	用途	优势
NumPy	数值计算	高效的数组操作
Matplotlib	数据可视化	丰富的图表类型
scikit-learn	机器学习	完善的算法实现

1.2 Iris数据集探索

Iris数据集是机器学习界的"Hello World"，包含三种鸢尾花的四个特征：

# 加载并查看数据集
iris = load_iris()
X = iris.data
y = iris.target
feature_names = iris.feature_names
target_names = iris.target_names

print(f"特征矩阵形状: {X.shape}")
print(f"特征名称: {feature_names}")
print(f"目标类别: {target_names}")

注意：虽然我们使用有标签的Iris数据集进行演示，但K-means作为无监督算法在实际应用中并不需要这些标签信息

数据预处理是确保算法效果的关键步骤：

# 数据标准化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 可视化特征分布
plt.figure(figsize=(12, 5))
for i in range(X.shape[1]):
    plt.subplot(1, 4, i+1)
    plt.hist(X[:, i], bins=20)
    plt.title(feature_names[i])
plt.tight_layout()
plt.show()

K-means算法实现详解

2.1 核心算法实现

让我们从零构建K-means算法的完整实现：

class KMeans:
    def __init__(self, n_clusters=3, max_iter=300, tol=1e-4):
        self.n_clusters = n_clusters
        self.max_iter = max_iter
        self.tol = tol
        
    def fit(self, X):
        # 初始化质心
        n_samples = X.shape[0]
        random_indices = np.random.choice(n_samples, self.n_clusters, replace=False)
        self.centroids = X[random_indices]
        
        for _ in range(self.max_iter):
            # 分配样本到最近的质心
            distances = self._compute_distances(X)
            self.labels_ = np.argmin(distances, axis=1)
            
            # 计算新质心
            new_centroids = np.array([
                X[self.labels_ == k].mean(axis=0) 
                for k in range(self.n_clusters)
            ])
            
            # 检查收敛
            centroid_shift = np.linalg.norm(new_centroids - self.centroids)
            if centroid_shift < self.tol:
                break
                
            self.centroids = new_centroids
            
        return self
    
    def _compute_distances(self, X):
        return np.array([
            np.linalg.norm(X - centroid, axis=1) 
            for centroid in self.centroids
        ]).T
    
    def predict(self, X):
        distances = self._compute_distances(X)
        return np.argmin(distances, axis=1)

2.2 关键数学原理

K-means本质上是优化问题，目标是最小化簇内平方和（WCSS）：

$$ J = \sum_{i=1}^{k} \sum_{x \in C_i} ||x - \mu_i||^2 $$

其中：

• $C_i$ 是第i个簇
• $\mu_i$ 是第i个簇的质心
• $||x - \mu_i||$ 表示数据点与质心的欧氏距离

距离计算优化技巧：

# 欧氏距离的向量化计算
def euclidean_distance(x, y):
    return np.sqrt(np.sum((x - y)**2, axis=1))

# 更高效的距离矩阵计算
pairwise_diff = X[:, np.newaxis, :] - self.centroids
distances = np.sqrt(np.sum(pairwise_diff**2, axis=2))

确定最佳K值：肘部法则与轮廓系数

3.1 肘部法则实践

寻找最佳簇数是K-means应用中的关键挑战。肘部法则通过观察WCSS随K值变化的拐点来确定：

wcss = []
k_range = range(1, 11)

for k in k_range:
    kmeans = KMeans(n_clusters=k)
    kmeans.fit(X_scaled)
    wcss.append(kmeans.inertia_)  # inertia_属性即WCSS

plt.figure(figsize=(10, 6))
plt.plot(k_range, wcss, 'bo-')
plt.xlabel('Number of clusters (K)')
plt.ylabel('Within-Cluster Sum of Squares (WCSS)')
plt.title('Elbow Method for Optimal K')
plt.xticks(k_range)
plt.grid()
plt.show()

3.2 轮廓系数分析

轮廓系数结合了簇内的凝聚度和簇间的分离度：

silhouette_scores = []
for k in range(2, 11):
    kmeans = KMeans(n_clusters=k)
    preds = kmeans.fit_predict(X_scaled)
    score = silhouette_score(X_scaled, preds)
    silhouette_scores.append(score)

plt.figure(figsize=(10, 6))
plt.plot(range(2, 11), silhouette_scores, 'go-')
plt.xlabel('Number of clusters (K)')
plt.ylabel('Silhouette Score')
plt.title('Silhouette Analysis for Optimal K')
plt.xticks(range(2, 11))
plt.grid()
plt.show()

评估指标对比：

方法	优点	局限性
肘部法则	直观易懂	拐点有时不明显
轮廓系数	综合考虑内外距离	计算复杂度较高
间隔统计	适合不同形状的簇	实现较复杂

聚类结果可视化与解读

4.1 二维特征投影

虽然Iris有四个特征，但我们可以通过降维可视化：

from sklearn.decomposition import PCA

# 降维到2D
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

# 聚类并可视化
kmeans = KMeans(n_clusters=3)
clusters = kmeans.fit_predict(X_scaled)

plt.figure(figsize=(12, 6))
plt.scatter(X_pca[:, 0], X_pca[:, 1], c=clusters, cmap='viridis', s=50)
plt.scatter(kmeans.centroids[:, 0], kmeans.centroids[:, 1], c='red', s=200, marker='X')
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.title('K-means Clustering Results (PCA Projection)')
plt.colorbar(label='Cluster')
plt.grid()
plt.show()

4.2 平行坐标图

多维度可视化技术能更好展示高维数据：

from pandas.plotting import parallel_coordinates

# 创建包含聚类结果的数据框
df = pd.DataFrame(X_scaled, columns=feature_names)
df['cluster'] = clusters

plt.figure(figsize=(12, 6))
parallel_coordinates(df, 'cluster', colormap='viridis')
plt.title('Parallel Coordinates Plot of Iris Features by Cluster')
plt.grid()
plt.show()

高级优化与实战技巧

5.1 K-means++初始化

原始K-means对初始质心敏感，K-means++提供了更好的初始化策略：

def initialize_centroids_plus(X, k):
    centroids = [X[np.random.randint(X.shape[0])]]
    
    for _ in range(1, k):
        distances = np.array([min([np.linalg.norm(x - c)**2 for c in centroids]) for x in X])
        probabilities = distances / distances.sum()
        next_centroid = X[np.random.choice(X.shape[0], p=probabilities)]
        centroids.append(next_centroid)
    
    return np.array(centroids)

5.2 处理不同尺度特征

当特征尺度差异大时，标准化至关重要：

# 对比标准化前后的聚类效果
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
kmeans_raw = KMeans(n_clusters=3).fit(X)
plt.scatter(X[:, 0], X[:, 1], c=kmeans_raw.labels_)
plt.title('Clustering on Raw Data')

plt.subplot(1, 2, 2)
kmeans_scaled = KMeans(n_clusters=3).fit(X_scaled)
plt.scatter(X_scaled[:, 0], X_scaled[:, 1], c=kmeans_scaled.labels_)
plt.title('Clustering on Scaled Data')

plt.tight_layout()
plt.show()

5.3 评估聚类质量

除了轮廓系数，还有多种评估方法：

from sklearn.metrics import calinski_harabasz_score, davies_bouldin_score

# 计算多种指标
def evaluate_clustering(X, labels):
    print(f"Silhouette Score: {silhouette_score(X, labels):.3f}")
    print(f"Calinski-Harabasz Index: {calinski_harabasz_score(X, labels):.3f}")
    print(f"Davies-Bouldin Index: {davies_bouldin_score(X, labels):.3f}")

evaluate_clustering(X_scaled, clusters)

常见问题与解决方案

6.1 空簇问题处理

当某个簇失去所有成员时，可采取以下策略：

def safe_mean(cluster_points):
    if len(cluster_points) == 0:
        return np.random.rand(X.shape[1])  # 随机生成新质心
    return cluster_points.mean(axis=0)

6.2 局部最优解

多次运行取最优结果是常用策略：

best_score = -1
best_model = None

for _ in range(10):
    model = KMeans(n_clusters=3)
    model.fit(X_scaled)
    score = silhouette_score(X_scaled, model.labels_)
    
    if score > best_score:
        best_score = score
        best_model = model

print(f"Best Silhouette Score: {best_score:.3f}")

6.3 大数据集优化

对于大规模数据，可考虑Mini-Batch K-means：

from sklearn.cluster import MiniBatchKMeans

mbkmeans = MiniBatchKMeans(n_clusters=3, batch_size=100)
mbkmeans.fit(X_scaled)

扩展应用与进阶方向

7.1 图像压缩

K-means可用于颜色量化，减少图像中的颜色数量：

from sklearn.datasets import load_sample_image

china = load_sample_image("china.jpg")
X_img = china.reshape(-1, 3) / 255.0

kmeans_img = KMeans(n_clusters=64).fit(X_img)
compressed = kmeans_img.centroids[kmeans_img.labels_].reshape(china.shape)

plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.imshow(china)
plt.title('Original Image')

plt.subplot(1, 2, 2)
plt.imshow(compressed)
plt.title('Compressed Image (64 colors)')
plt.show()

7.2 时间序列聚类

通过适当距离度量，K-means可应用于时间序列：

from scipy.spatial.distance import cdist

def dtw_distance(x, y):
    """动态时间规整距离"""
    # 实现略
    pass

class TimeSeriesKMeans(KMeans):
    def _compute_distances(self, X):
        return np.array([
            [dtw_distance(x, centroid) for centroid in self.centroids]
            for x in X
        ])

7.3 与其他算法结合

K-means常作为特征工程步骤：

from sklearn.pipeline import Pipeline
from sklearn.ensemble import RandomForestClassifier

pipeline = Pipeline([
    ('kmeans', KMeans(n_clusters=10)),
    ('classifier', RandomForestClassifier())
])