欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

K-Means 概要及其实现代码

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

本文先整体介绍 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)i(1,2,,n);l(1,2,,k)(1)l = C(i) \tag{1} \quad i\in(1,2,…,n);l\in(1,2,…,k)
函数 CC 是一个多对一的划分函数,输入样本 ii,输出类别 ll

1.2.2 策略

相似的样本被聚到同类,即让相同类中样本相似程度最小,故 K-Means 就是求解最优化的问题。
目标函数(csdn 里对 Latex语法的支持度不是很好,就用了截图):
K-Means 概要及其实现代码
其中 W(C)W(C) 是损失函数,表示样本与其所属类的中心之间的距离总和(用这来表示相同类中的样本相似程度)。
那怎么让这个损失函数最小呢?
这是个组合优化问题,n 个样本分到 k 类的分法是指数级别的。现实中解决这种 NP 困难问题,一般采用迭代的方式求解。

1.2.3 算法

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

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

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

1.3 算法缺点的解决方案

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

1.3.1 类别数 k 需要事先指定

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

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

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

在指定了类别数 kk 后,需要指定 kk 个初始中心。
K-Means 初始中心是从样本点随机选取的。但是在数据量增大、类别数增多时,随机选定初始中心常会出现:一个较大子集有多个中心点,而其他多个较小的子集共用一个中心点的情况(如下图所示)。
K-Means 概要及其实现代码
原因是初始中心距离太近。那选取距离最远的点作为初始中心岂不是更简单直接,但是这样又可能会因为数据集中的离群点干扰影响聚类效果,所以这里引入 K-Means++ 算法,该算法对初始化中心环节做了改进,从源头改善局部最优解的现象。
解决方案:

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

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

还有一种解决方案来确定初始化聚类中心。
先用层次聚类对样本进行聚类,得到想要的 kk 个类时停止,然后从每个类中选取一个与中心距离最近的点,用这些点作为 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 相似度量函数

给定样本集合XXXX是m维实数向量空间RmR^m中点的集合,其中xi,xjXx_i,x_j \in X

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

注:当SS位单位矩阵时,马氏距离就是欧式距离.

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 值越小时,越满足「不同划分中的样本差异尽可能的大」

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

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)  # 将数据可视化展示

K-Means 概要及其实现代码
这个数据很干净,直观的就可以看出分为5个类别。
2.选择类别数 k

plot_k_vs_inertia_udf(df)

plot_k_vs_inertia(df)

K-Means 概要及其实现代码
左图是 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")

K-Means 概要及其实现代码
可以看到不论是自己写的函数(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