K-Means 概要及其实现代码

本文先整体介绍 K-Means 算法的原理及特点，然后从模型、策略、算法三个要素切入介绍这个聚类算法，紧接着讲了几点 K-Means 算法缺点的解决方案；之后从 UDF 和框架两方面给出 K-Means 在 Python 中的实现方法；最后引入一个实际案例进行实操。

一、K-Means 概要

K-Means 是基于样本集合划分的聚类算法。它将样本集合划分为 k 个子集，构成了 k 个类，并将 n 个样本分到 k 个类中去，每个样本到其所属类的中心的距离最小。

1.1 特点

K-Means ：

1. 是基于划分的聚类，属于硬聚类方法；
2. 类别数k事先指定（如何指定？见下文）；
3. 距离或相似度的度量是 欧式距离平方 ，以中心或样本的均值表示类别；
4. 用于优化的目标函数：样本和其所属类的中心之间的距离的总和（SSE）；
5. 算法是迭代算法，不能保证得到全局最优。

K-Means 的优缺点：

• 优点：
  • 原理简单，实现容易；
  • 计算速度相对较快；
• 缺点：
  • 类别数 k 需要事先指定；
  • 初始中心的选择会直接影响聚类结果；
  • 不能保证收敛到全局最优（局部最优）；
  • 对噪声和异常点敏感；
  • 对于非凸数据集或类别规模差异太大的数据效果较差（适用于球状类别）；
  • 大规模的数据集上效率较低；
  • 只能使用连续变量进行聚类。

1.2 模型、策略、算法

1.2.1 模型

K-Means 模型的本质是一个划分函数，输入样本，输出样本对应的类别。
假设数据集有 n 行观测，拟聚为 k 类。
$l = C(i) \tag{1} \quad i\in(1,2,…,n);l\in(1,2,…,k)$
函数 $C$ 是一个多对一的划分函数，输入样本 $i$，输出类别 $l$。

1.2.2 策略

相似的样本被聚到同类，即让相同类中样本相似程度最小，故 K-Means 就是求解最优化的问题。
目标函数（csdn 里对 Latex语法的支持度不是很好，就用了截图）：

其中 $W(C)$ 是损失函数，表示样本与其所属类的中心之间的距离总和（用这来表示相同类中的样本相似程度）。
那怎么让这个损失函数最小呢？
这是个组合优化问题，n 个样本分到 k 类的分法是指数级别的。现实中解决这种 NP 困难问题，一般采用迭代的方式求解。

1.2.3 算法

上面说了 K-Means 是求解最优化的问题，现实中会采用迭代的方式求解。

K-Means 算法：
拟定分为 k 类，确定距离度量为 欧式距离平方，通常还需要设置迭代停止条件（如最大迭代次数、误差）。

1. 初始化质心，随机选用数据集中的 k 个样本作为初始类的中心（这些中心点代表了其所属的类）；
2. 计算各样本点到各中心的距离，将各样本点划分到离它们最近的类；
3. 重新计算各类的中心，再重复步骤2。
   ① 如果迭代收敛（重新计算中心后的划分结果与之前一致）或符合迭代停止条件，则输出 $C^*$ （最小化的误差平方和）和$l$ （划分的结果）；
   ② 否则重复步骤2、3。

1.3 算法缺点的解决方案

实际运用中建议跑多次，最好再试试其他聚类算法，结合业务意义与误差平方法两个维度选取最佳的结果。

1.3.1 类别数 k 需要事先指定

K-Means 聚类前需要事先指定类别数 $k$，实际运用中最优的 $k$ 值是未知的，所以通常比较难确定。
解决方案：

尝试使用不同的 $k$ 值聚类，根据聚类结果（通常用类的平均直径衡量），绘制出 类别数 $k$- 平均直径 的折线图，找到拐点，再结合实际业务在拐点或其附近确定相应的 $k$ 值。

1.3.2 聚类结果依赖于初始中心的选择

在指定了类别数 $k$ 后，需要指定 $k$ 个初始中心。
K-Means 初始中心是从样本点随机选取的。但是在数据量增大、类别数增多时，随机选定初始中心常会出现：一个较大子集有多个中心点，而其他多个较小的子集共用一个中心点的情况（如下图所示）。

原因是初始中心距离太近。那选取距离最远的点作为初始中心岂不是更简单直接，但是这样又可能会因为数据集中的离群点干扰影响聚类效果，所以这里引入 K-Means++ 算法，该算法对初始化中心环节做了改进，从源头改善局部最优解的现象。
解决方案：

K-Means++ 算法
在确认了类别数 k 的前提下：

1. 从数据集中随机选择一个观测作为某个类的中心。计算每个观测点到该点的欧氏距离平方 $D_x$（并放入一个数组中），并将这些距离求和得到 $D_{total}$；
2. 从 $[0, D_{total})$ 区间随机选择一个值 $R$，依次重复计算 $R = R - D_x$ ，直到 $R≤0$ ，选择此时 $D_x$ 对应的点作为下一个中心点；
3. 重复第 2 步 ，直至 $k$ 个初始聚类中心都确定。

还有一种解决方案来确定初始化聚类中心。
先用层次聚类对样本进行聚类，得到想要的 $k$ 个类时停止，然后从每个类中选取一个与中心距离最近的点，用这些点作为 K-Means 聚类的初始中心点。
但是这种方法在数据集大时效率很低。

1.3.3 大规模的数据集上效率较低

毫无疑问的是 K-Means 算法的计算量是随着数据量递增的，在面对大规模的数据集时，K-Means的效率较低。
解决方案：

Mini-Batch K-Means：
算法原理：在每次迭代过程中，从数据集中随机抽取一部分以形成小批量的数据集，用该子集计算距离、更新类中心点。

1.3.4 只能使用连续变量进行聚类

分类变量可以通过独热编码（最常用）、虚拟编码等其他形式，作为连续变量输入。

二、K-Means 实现

使用工具：Python.

2.1 UDF 实现

2.1.1 相似度量函数

给定样本集合$X$，$X$是m维实数向量空间$R^m$中点的集合，其中$x_i,x_j \in X$。

• 闵氏距离：
  $d_{ij} = (\sum_{k=1}^m |x_{ik} - x_{jk}|^p)^\frac{1}{p} \tag{1}$
• 马氏距离：
  $d_{ij} = [(x_i - x_j)^T S^{-1} (x_i - x_j)] ^ \frac{1}{2} \tag{2}$

注：当$S$位单位矩阵时，马氏距离就是欧式距离.

import numpy as np
import pandas as pd

def dist_measure(x, y, typ='MIN', p=2, X=''):
    """
    距离度量函数（衡量观测与观测之间的相似程度）
	
    输入：
    x, y -- 两个观测点（numpy 矩阵对象）
    typ  -- 'MIN'表示闵氏距离；'MA'表示马氏距离。默认闵氏距离（也可用np.linalg.norm(x-y, ord=2)实现）
    p    -- 闵氏距离的参数：1表示曼哈顿距离；2表示欧氏距离；'inf'表示切比雪夫距离
    X    -- 马氏距离的参数：数据集的矩阵，要求样本量n（行数）大于变量数m（列数），必须指定
    返回：
    dist -- 指定的距离（float对象）
    """
    if typ == 'MIN':
        if type(p) == int:
            dist = np.power(np.sum(np.power(np.abs(x - y), p)), 1/p)
            return dist
        elif p == 'inf':
            dist = np.max(np.abs(x - y))    # 当 p 取无穷大时，是切比雪夫距离
            return dist
        else:
            print('错误：闵氏距离的p要求大于1的整数，输入不符合要求!')
            return None
    elif typ == 'MA':
        x = x.T
        y = y.T    # 行向量（样品）的转置
        Sigma = np.cov(X.T)    # np.cov(M) M的行是变量，列是样品，故需要转置
        try:
            Sigma_inv = np.linalg.inv(Sigma)
            dist = np.power((x - y).T @ Sigma_inv @ (x - y), 1/2)
            # 矩阵乘法：
            # matrix中@ * dot等价（np.multiply表示数量积）；
            # array中只能使用@和dot（普通的*和np.multiply表示数量积）
            return dist
        except Exception as e:
            print('错误：协方差矩阵不是半正定矩阵，即样本矩阵的行数小于等于列数.')
            return None
    else:
        print('错误：输入方法错误！只能是MIN或MA.')
        return None

# 测试：求观测x与观测y的欧式距离平方
x = np.mat([1, 2])
y = np.mat([2, 1])
print(dist_measure(x, y))   # 1.414 ✔

2.1.2 初始化中心函数

下方 UDF 支持 random 和 kmeans++ 两种初始化中心的方法：

# K-Means++ 支持函数（得到每个样本距离最近中心点的距离及距离之和）
def get_sum_distance(centers, dataset):
    """
    参数:
    centers -- 中心点集合
    dataset    -- 数据集

    返回:
    np.sum(dis_lst) -- 样本距离最近中心点的距离之和
    dis_lst         -- 每个样本距离最近中心点的距离（列表形式存放）
    """
    dis_lst = np.array([])
    for each_data in dataset:
        distances = np.array([])
        for each_center in centers:
            temp_distance = dist_measure(each_data, each_center)  # 计算样本和中心点的欧式距离
            distances = np.append(distances, temp_distance)
        lab = np.min(distances)
        dis_lst = np.append(dis_lst, lab)
    return np.sum(dis_lst), dis_lst


def init_center(dataset, k, typ=2, random_state=42):
    """
    初始化质心（初始质心不能相同）

    输入：
    dataset   -- 需要聚类的数据集（DataFrame或Matrix）
    k         -- 事先指定的聚类的数（int对象）
    typ       -- 类型1时是"random"，即从数据集中从随机抽取 k 行作为初始质心；
                 类型2时是"kmeans++"

    返回：
    centroids -- 初始质心（k行m列的np.array对象，m由数据集决定）

    """
    n = dataset.shape[0]    # 行（观测）
    m = dataset.shape[1]    # 列（变量）
    if typ == 1:
        rng = np.random.RandomState(random_state)
        indices = rng.randint(n, size=k)
        centroids = np.matrix(dataset[indices])
    elif typ == 2:
        dataset = np.array(dataset)
        # 先从数据集随机选一个观测作为某一类的中心
        rng = np.random.RandomState(random_state)
        p = rng.randint(0, len(dataset))
        first_center = dataset[p]
        center_lst = []
        center_lst.append(first_center)       

        for i in range(k-1):
            sum_dis, dis_lst = get_sum_distance(center_lst, dataset)
            r = np.random.randint(0, sum_dis)
            for j in range(len(dis_lst)):
                r = r - dis_lst[j]
                if r <= 0:
                    center_lst.append(dataset[j])
                    break
                else:
                    pass
        centroids = np.mat(center_lst)
    
    return centroids

2.1.3 K-Means 聚类函数

def K_Means(df, k=5, typ=2, random_state=-1, print_typ=False):
    """
    K-均值聚类函数
    
    不足：
    	 - 不能设置迭代停止条件
    	 - 在 n 行的数据集要分为 n 类时可能存在问题
    输入：
    df             --  数据集（DataFrame对象，各变量只能是数值型）
    k              --  类别数（默认为 5）
    typ            --  初始化中心的方法（默认是2。1表示随机抽 k 个样本作为中心，2是 kmeans++ 算法的初始化方法）
    random_state   --  随机种子（用来固定聚类结果）
    print_typ      --  是否打印出每次迭代后中心更新的情况（默认关闭）
    
    返回：
    centroids      --  最后的类中心（numpy.matrix 对象，shape是 n * m）
    cluster_result --  第一列是各观测对应的类别（np.squeeze(np.array(cluster_result[:, 0]))）
                       第二列是各观测到各自类中心的距离（cluster_result[:, 1].sum()）
    
    """
    if random_state == -1:
        random_state = np.random.randint(0, 1000)  # 若不指定随机数种子，则每次划分的结果都可能不一样
    
    dataset = np.mat(df.values)
    n = dataset.shape[0]    # n条观测
    centroids = init_center(dataset, k, typ, random_state)    # 初始聚类中心（k行m列）
    cluster_result = np.mat(np.zeros((n, 2)))  # 初始化聚类结果(n行2列，第一列存放n行观测的类别，第二列存放最小距离（欧式距离平方）)
    
    cluster_changed = True 
    t = 0    # 迭代次数（质心更新次数）

    while cluster_changed:
        cluster_changed = False
        # 1.寻找最近的质心
        for i in range(n):    # n是总观测数
            min_index = -1; min_dist = np.inf    # 初始化：min_dist是每一条观测的最近距离，min_index是每一条观测所属类别（聚类中心的索引）
            for j in range(k):  # k是总类别数，n>=k
                dist_ij = dist_measure(dataset[i, :], centroids[j, :])
                if dist_ij < min_dist:
                    min_index = j; min_dist = dist_ij
            if cluster_result[i, 0] != min_index: cluster_changed = True    # 只要有一条观测的类别发生变化，循环就不停止
            cluster_result[i, :] = min_index, min_dist

        # 2.更新各类质心的位置
        if print_typ == True:
            print(f'-------------\n第{t}次更新后质心的位置：\n{centroids}')
        for clust in range(k):
            points_in_clust = dataset[np.nonzero(cluster_result[:, 0] == clust)[0]]  # np.nonzero返回数组 a 中非零元素的索引值数组(tuple)
            centroids[clust] = np.mean(points_in_clust, axis=0)
        t += 1
#    print(f'--------共迭代了{t-1}次')  
    return centroids, cluster_result

2.1.4 k 值的选择

• 当 K 值越大时，越满足「同一划分的样本之间的差异尽可能的小」;
• 当 K 值越小时，越满足「不同划分中的样本差异尽可能的大」

选择合适的 $k$ 是很重要的，参考 1.3.1 ，绘制类别数 $k$ 与误差（这里定义各观测和其最近的类中心距离之和）的折线图，找到拐点。具体结合业务意义，选择拐点或拐点附近的点作为选定的 $k$ 值。

def plot_k_vs_inertia_udf(df, k_min = 1, k_max = 10):
    """
    K-Means UDF k 值选取函数
    
    具体是绘制出类别数k与误差（这里定义各观测和其最近的类中心距离之和）的折线图
    
    """
    index = []  # 横坐标数组（k的取值）
    inertia = []  # 纵坐标数组（误差，即各观测和其最近的类中心距离之和）

    for k in range(k_min, k_max+1):
        centroids, cluster_result = K_Means(df, k=k, typ=2)
        index.append(k)
        inertia.append(cluster_result[:, 1].sum())

    # 绘制折线图
    plt.plot(index, inertia, "-o")
    plt.grid()

2.2 Sklearn 实现

class sklearn.cluster.KMeans 重要参数与类属性：

参数（一般只需要调整类别数 n_clusters 即可）：

• n_clusters = 8; # 默认聚成 8 类
• init = ‘k-means++’ : ‘k-means++’/‘random’/ndarray，初始类中心位置
  • ‘k-means++’ : 采用优化后的算法确定类中心
  • ‘random’ : 随机选取k个案例作为初始类中心
  • ndarray : (n_clusters, n_features)格式提供的初始类中心位
• precompute_distances = ‘auto’ : {‘auto’, True, False} 是否预先计算距离，分析速度更快，但需要更多内存
  ‘auto’ : 如果n_samples*n_clusters > 12 million，则不事先计算距离
• algorithm = ‘auto’ : ‘auto’, ‘full’ or ‘elkan’，具体使用的算法
  ‘full’ : 经典的EM风格算法
  ‘elkan’ : 使用三角不等式，速度更快，但不支持稀疏数据
  ‘auto’ : 基于数据类型自动选择

常用用法：

# k 值的选取，找拐点
def plot_k_vs_inertia(df, k_min = 1, k_max = 10):
    """
    K-Means k 值选取函数（基于 Sklearn 框架）
    不足：折线图有点吃藕（丑），后续优化下
    具体是绘制出类别数k与误差（这里定义各观测和其最近的类中心距离之和）的折线图
    
    """
    index = []  # 横坐标数组（k的取值）
    inertia = []  # 纵坐标数组（误差，即各观测和其最近的类中心距离之和）

    for i in range(k_min, k_max+1):
        model = KMeans(n_clusters=i).fit(df)
        index.append(i)
        inertia.append(model.inertia_)

    # 绘制折线图
    plt.plot(index, inertia, "-o")
    plt.grid()

假如发现拐点在 k=5：

# sklearn 实现
from sklearn.cluster import KMeans
# 将 df 聚成 5 类
kmeans = KMeans(n_clusters=5, random_state=0).fit(df)
kmeans.labels_    # 聚类结果（即各行被划分的类别），(n,)的array
# kmeans.cluster_centers_  # 最后各类别的中心，array (n_clusters, n_features)
# kmeans.inertia_ # 误差（各观测与其最近中心的距离之和）

对于样本量超过1万的情形，建议使用MiniBatchKMeans，算法改进为在线增量，速度会更快。

from sklearn.cluster import MiniBatchKMeans

mini_kmeans = MiniBatchKMeans(n_clusters=5, random_state=0).fit(df)
mini_kmeans.labels_

2.3 函数效果验证

1.获取数据：用 sklearn 框架里的 make_blobs 生成测试聚类算法的玩具数据（1000行2列），这样会直观点：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import make_blobs
from sklearn.cluster import KMeans, MiniBatchKMeans
%matplotlib inline

plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False   #用来正常显示负号
plt.style.use("tableau-colorblind10")

data,label = make_blobs(n_samples=100,n_features=2,centers=3,center_box=(-10.0,10.0),random_state=None)
df, _ = make_blobs(n_samples=1000, centers=5, random_state=18)  # 生成数据
df = pd.DataFrame(df)
plt.scatter(df.iloc[:, 0], df.iloc[:, 1], s=20)  # 将数据可视化展示

这个数据很干净，直观的就可以看出分为5个类别。
2.选择类别数 k ：

plot_k_vs_inertia_udf(df)

plot_k_vs_inertia(df)

左图是 K-Means UDF 的结果，右图是 Sklearn 框架中 KMeans 的结果。拐点都是 k = 5 时，直接确定 k = 5。
3.聚类划分样本：这里因为数据集是二维的，就用散点图做呈现了，也更直观点。

fig, axes = plt.subplots(1, 2, figsize=(15, 5))

# UDF
centers1, labels1 = K_Means(df, 5) 
centers1 = np.squeeze(np.array(centers1))
labels1 = np.squeeze(np.array(labels1[:, 0]))

axes[0].scatter(df.iloc[:, 0], df.iloc[:, 1], s=20, c=labels1)
axes[0].scatter(centers1[:, 0], centers1[:, 1], s=100, marker='*', c="r")

# Sklearn
kmeans = KMeans(n_clusters=5, random_state=0).fit(df)
centers2, labels2 = kmeans.cluster_centers_, kmeans.labels_

axes[1].scatter(df.iloc[:, 0], df.iloc[:, 1], s=20, c=labels2)
axes[1].scatter(centers2[:, 0], centers2[:, 1], s=100, marker='*', c="r")

可以看到不论是自己写的函数（UDF）还是套用的框架（sklearn），都精准地将这份数据集划分成了5个类别。

注： sklearn 的拟合速度要快很多。此外由于数据比较干净，考虑点较少。实际运用时，k 值最好不光看拐点，再还需结合具体业务考虑，需要多跑几次或多用几种不同的聚类算法，以期得到更合适的结果。

三、案例实操

本文的篇幅过长，因为放在另一篇博文中：K-Means 案例实操

参考：
[1] 李航. 统计学习方法（第二版）[M]
[2] Peter Harrington. 机器学习实战 [M]
[3] *. K-Means算法
[4] 张文彤. Python数据分析-玩转数据挖掘

本文地址：https://blog.csdn.net/xienan_ds_zj/article/details/107116215