6.2. 吴恩达机器学习-Kmeans

1. 理论基础

1. 算法

K-means 是我们最常用的基于欧氏距离的聚类算法，其认为两个目标的距离越近，相似度就越大。

1.1 牧师-村民模型

K-means有一个著名的解释：牧师—村民模型：
  有四个牧师去郊区布道，一开始牧师们随意选了几个布道点，并且把这几个布道点的情况公告给了郊区所有的村民，于是每个村民到离自己家最近的布道点去听课。
  听课之后，大家觉得距离太远了，于是每个牧师统计了一下自己的课上所有的村民的地址，搬到了所有地址的中心地带，并且在海报上更新了自己的布道点的位置。
  牧师每一次移动不可能离所有人都更近，有的人发现A牧师移动以后自己还不如去B牧师处听课更近，于是每个村民又去了离自己最近的布道点……
  就这样，牧师每个礼拜更新自己的位置，村民根据自己的情况选择布道点，最终稳定了下来。
我们可以看到该牧师的目的是为了让每个村民到其最近中心点的距离和最小。

1.2 算法步骤

所以K-means的算法步骤为：
  1. 随机初始化 k 个聚类中心 $a=a_1,a_2,...,a_k$
  2. 针对数据集中每个样本 $x_i$ 计算它到 k 个聚类中心的距离并将其分到距离最小的聚类中心所对应的类中
  3. 针对每个类别 $a_j$，重新计算它的聚类中心 $a_j=\frac{1}{|c_i|}{\sum }_{x\in c_i}x$
  4. 重复上面 2 3 两步操作，直到达到某个中止条件（迭代次数、最小误差变化等）

1.3 复杂度

我们先看下伪代码：

获取数据（n 个 m 维的数据）
随机生成 k 个 m 维的数据
while(t)
	for(int i = 0; i < n; i++)
		for(int j = 0; j < k; j++)
			计算点 i 到类 j 的距离

	for(int i = 0; i < k; i++)
		1. 找出所有属于自己这一类的所有数据点
		2. 把自己的坐标修改为这些数据点的中心点坐标
end
1
2
3
4
5
6
7
8
9
10
11

时间复杂度：$O(tknm)$，其中，t为迭代次数，k为簇的数目，n为样本点数，m为样本点维度
空间复杂度：$O(m(n+k))$，其中，k为簇的数目，m为样本点维度，n为样本点数

2. 优缺点

2.1 优点

• 容易理解，聚类效果不错，虽然是局部最优，但往往局部最优就够了
• 处理大数据集的时候，该算法可以保证较好的伸缩性
• 当簇近似高斯分布的时候，效果非常不错
• 算法复杂度低

2.2 缺点

• K 值需要人为设定，不同 K 值得到的结果不一样
• 对初始的簇中心敏感，不同选取方式会得到不同结果
• 对异常值敏感
• 样本只能归为一类，不适合多分类任务
• 不适合太离散的、样本类别不平衡的、非凸形状的分类

3. 算法调优与改进

  针对 K-means 算法的缺点，我们可以有很多种调优方式：如数据预处理（去除异常点），合理选择 K 值，高位映射等。以下将简单介绍：

3.1 数据预处理

  K-means 的本质是基于欧式距离的数据划分算法，均值和方差大的维度将对数据的聚类产生决定性影响。所以未做归一化处理和统一单位的数据是无法直接参与运算和比较的。常见的数据预处理方式有：数据归一化，数据标准化。
  此外，离群点或者噪声数据会对均值产生较大的影响，导致中心偏移，因此我们还需要对数据进行异常点检测。

3.2 合理选择 K 值

  K 值的选取对 K-means 影响很大，这也是 K-means 最大的缺点，常见的选取 K 值的方法有：手肘法、Gap statistic方法。
手肘法：

  当 K < 3 时，曲线急速下降；当 K > 3 时，曲线趋于平稳，通过手肘法我们认为拐点 3 为 K 的最佳值。
  手肘法的缺点在于需要人工看不够自动化，所以我们又有了 Gap statistic 方法，这个方法出自斯坦福大学的几个学者的论文：Estimating the number of clusters in a data set via the gap statistic$$Gap(K)=E(\log D_k)-\log D_k$$  其中 $D_k$ 为损失函数，这里 $E(\log D_k)$ 指的是 $\log D_k$ 的期望。这个数值通常通过蒙特卡洛模拟产生，我们在样本里所在的区域中按照均匀分布随机产生和原始样本数一样多的随机样本，并对这个随机样本做 K-means，从而得到一个 $D_k$。如此往复多次，通常 20 次，我们可以得到 20 个 $\log D_k$。对这 20 个数值求平均值，就得到了 $E(\log D_k)$ 的近似值。最终可以计算 Gap statistic。而 Gap statistic 取得最大值所对应的 K 就是最佳的 K。

  由图可见，当 K=3 时，Gap(K)取值最大，所以最佳的簇数是 K=3。

这里解释一下为什么 Gap(K) 取最大值所对应的 K 值就是最佳的簇数？
  关于这一点，参考论文原文图1可以得到很好的理解。当 K 小于真实的类簇数目时，随着 K 的增大，类簇间样本的划分粒度逐渐增大，则 $D_k$ 下降的幅度较大，也就是迅速下降，因此 $\log (D_k)$ 的下降幅度也较大。当 K 大于真实的类簇数目时，则 $D_k$ 和 $\log (D_k)$ 的下降幅度逐渐平缓，因此在真实的 K 值附近，$D_k$ 和$\log (D_k)$ 的下降会出现一个明显的拐点。然而，$E(\log (D_k))$ 却由于平均的原因，使得其变化幅度变得平滑，随着 K 的增大，$E(\log (D_k))$ 保持着近乎均匀的变化幅度。所以 $E(\log (D_k))-\log (D_k)$ 会在真实的 K 值时达到最大。因此，当Gap(K) 最大时所对应的就是最佳的类簇数目。
  Github 上有一个项目叫 gap_statistic，可以更方便的获取建议的类簇个数。

3.3 采用核函数

  基于欧氏距离的 K-means 假设了各个数据簇的数据具有一样的先验概率并呈现球形分布，但这种分布在实际生活中并不常见。面对非凸的数据分布形状时我们可以引入核函数来优化，这时算法又称为 核K-means 算法，是核聚类方法的一种。核聚类方法的主要思想是通过一个非线性映射，将输入空间中的数据点映射到高维的特征空间中，并在新的特征空间中进行聚类。非线性映射增加了数据点线性可分的概率，从而在经典的聚类算法失效的情况下，通过引入核函数可以达到更为准确的聚类结果。

3.4 K-means++

  我们知道初始值的选取对结果的影响很大，对初始值选择的改进是很重要的一部分。在所有的改进算法中，K-means++ 最有名。
K-means++ 算法步骤如下：
  1. 随机选取一个中心点 $a_1$
  2. 计算数据到之前 n 个聚类中心最远的距离 $D(x)$，并以一定概率 $\frac{D(x{)}^2}{\sum D(x{)}^2}$ 选择新中心点 $a_i$
  3. 重复第二步
  简单来说，K-means++ 就是选择离已选中心点最远的点。这也比较符合常理，各个聚类中心当然是离得越远越好。但是这个算法的缺点在于，难以并行化。所以 K-means || 改变取样策略，并非按照 K-means++ 那样每次遍历只取样一个样本，而是每次遍历取样 k 个，重复该取样过程 $\log (n)$ 次，则得到 $k\log (n)$ 个样本点组成的集合，然后从这些点中选取 k 个。当然一般也不需要 $\log (n)$ 次取样，5 次即可。

3.5 ISODATA

  ISODATA 的全称是迭代自组织数据分析法。它解决了 K 值需要预先人为设定这一缺点。而当遇到高维度、海量的数据集时，人们往往很难准确地估计出 K 的大小。ISODATA 就是针对这个问题进行了改进，它的思想也很直观：当属于某个类别的样本数过少时把这个类别去除，当属于某个类别的样本数过多、分散程度较大时把这个类别分为两个子类别。

4. 收敛证明

  我们先来看一下 K-means 算法的步骤：先随机选择初始节点，然后计算每个样本所属类别，然后通过类别再更新初始化节点。
  我们需要知道的是 K-means 聚类的迭代算法实际上是 EM 算法。EM 算法解决的是在概率模型中含有无法观测的隐含变量情况下的参数估计问题。在 K-means 中的隐变量是每个类别所属类别。K-means 算法迭代步骤中的 每次确认中心点以后重新进行标记 对应 EM 算法中的 E 步 求当前参数条件下的 Expectation。而 根据标记重新求中心点 对应 EM 算法中的 M 步 求似然函数最大化时（损失函数最小时）对应的参数。
首先我们看一下损失函数的形式：$$J=\sum _{i=1}^C\sum _{j=1}^N{r}_{ij}\cdot v(x_j,{\mu }_i)$$
其中： v ( x j , μ i ) = ∣ ∣ x j − μ i ∣ ∣ 2 , r n k = { 1 i f x n ∈ k 0 e l s e v(x_j, \mu_i)=||x_j-\mu_i||^2, r_{nk}=\left\{

\right. v(xj​,μi​)=∣∣xj​−μi​∣∣2,rnk​={10​ifxn​∈kelse​
为了求极值，我们令损失函数求偏导数且等于 0：$$\frac{\partial J}{\partial {\mu }_k}=2\sum _{i=1}^N{r}_{ik}(x_i-{\mu }_k)=0$$

k 是指第 k 个中心点，于是我们有：$${\mu }_k=\frac{{\sum }_{i=1}^N{r}_{ik}{x}_i}{{\sum }_{i=1}^N{r}_{ik}}$$
可以看出，新的中心点就是所有该类的质心。
EM 算法的缺点就是，容易陷入局部极小值，这也是 K-means 有时会得到局部最优解的原因。

2. 代码实现

1. 数据可视化

import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import scipy.io as sio

mat1 = sio.loadmat('../data/ex7data1.mat')

data1 = pd.DataFrame(mat1.get('X'), columns=['X1', 'X2'])
data1.head()
1
2
3
4
5
6
7
8
9

sns.set(context="notebook", style="white")
sns.lmplot('X1', 'X2', data=data1, fit_reg=False)
plt.show()
1
2
3

mat2 = sio.loadmat('../data/ex7data2.mat')
data2 = pd.DataFrame(mat2.get('X'), columns=['X1', 'X2'])
data2.head()
1
2
3

sns.lmplot('X1', 'X2', data=data2, fit_reg=False)
plt.show()
1
2

2. 二维K-means

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd
import scipy.io as sio

mat = sio.loadmat('../data/ex7data2.mat')
data2 = pd.DataFrame(mat.get('X'), columns=['X1', 'X2'])

sns.set(context="notebook", style="white")
sns.lmplot('X1', 'X2', data=data2, fit_reg=False)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12

1. 实现K-means算法

# 将数据集中的样本加上分类标签
def combine_data_C(data, C):
    data_with_c = data.copy()
    data_with_c['C'] = C
    return data_with_c


# k-means fn --------------------------------
# 随机初始化k个中心点
def random_init(data, k):
    """choose k sample from data set as init centroids
    Args:
        data: DataFrame
        k: int
    Returns:
        k samples: ndarray
    """
    return data.sample(k).iloc[:,:].values

# 将样本点x分类
def _find_your_cluster(x, centroids):
    """find the right cluster for x with respect to shortest distance
    Args:
        x: ndarray (n, ) -> n features
        centroids: ndarray (k, n)
    Returns:
        k: int
    """
    
    distances = np.apply_along_axis(func1d=np.linalg.norm,  # 自定义函数，这里是求范数（即距离）
                                    axis=1,             # 坐标轴
                                    arr=centroids - x)  # 输入的数组
    return np.argmin(distances)

# 将数据进行分类
def assign_cluster(data, centroids):
    """assign cluster for each node in data
       return C ndarray
    """
    return np.apply_along_axis(lambda x: _find_your_cluster(x, centroids), axis=1, arr=data.iloc[:,:].values)

# 更新中心点
def new_centroids(data, C):
    data_with_c = combine_data_C(data, C)

    return data_with_c.groupby('C', as_index=False).mean().sort_values(by='C').drop('C', axis=1).iloc[:,:].values

# 定义损失函数
def cost(data, centroids, C):
    m = data.shape[0]

    expand_C_with_centroids = centroids[C]

    distances = np.apply_along_axis(func1d=np.linalg.norm,
                                    axis=1,
                                    arr=data.iloc[:,:].values - expand_C_with_centroids)
    return distances.sum() / m

# 对各个参数进行迭代
def _k_means_iter(data, k, epoch=100, tol=0.0001):
    """one shot k-means
       with early break
    """
    centroids = random_init(data, k)
    cost_progress = []

    for i in range(epoch):
        print('running epoch {}'.format(i))

        C = assign_cluster(data, centroids)
        centroids = new_centroids(data, C)
        cost_progress.append(cost(data, centroids, C))

        if len(cost_progress) > 1:  # early break
            if (np.abs(cost_progress[-1] - cost_progress[-2])) / cost_progress[-1] < tol:
                break

    return C, centroids, cost_progress[-1]

# k-means算法
def k_means(data, k, epoch=100, n_init=10):
    """do multiple random init and pick the best one to return
    Args:
        data (pd.DataFrame)
    Returns:
        (C, centroids, least_cost)
    """

    tries = np.array([_k_means_iter(data, k, epoch) for _ in range(n_init)])

    least_cost_idx = np.argmin(tries[:, -1])

    return tries[least_cost_idx]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

random_init(data2, 3)
1

2. 聚类任务

# 初始化 3 个中心点
init_centroids = random_init(data2, 3)

# 画出那 3 个中心点所在位置
x = np.array([1,1])

fig, ax = plt.subplots(figsize=(6, 4))
ax.scatter(x=init_centroids[:, 0], y=init_centroids[:, 1])

for i, node in enumerate(init_centroids):
    ax.annotate(f'{i}: ({node[0]}, {node[1]})', node)
    
ax.scatter(x[0], x[1], marker='x', s=200)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14

_find_your_cluster(x, init_centroids)    # 测试一下样本 x 属于哪一类
1

第一轮聚类任务

# 将数据聚类
C = assign_cluster(data2, init_centroids)    
data_with_c = combine_data_C(data2, C)
data_with_c.head()
1
2
3
4

sns.lmplot('X1', 'X2', hue='C', data=data_with_c, fit_reg=False)
plt.show()
1
2

3. 更新中心点

new_centroids(data2, C)

final_C, final_centroid, _ = _k_means_iter(data2, 3)
data_with_c = combine_data_C(data2, final_C)

sns.lmplot('X1', 'X2', hue='C', data=data_with_c, fit_reg=False)
plt.show()
1
2
3
4
5
6
7

计算损失值

cost(data2, final_centroid, final_C)    # 0.794176363371587

best_C, best_centroids, least_cost = k_means(data2, 3)
least_cost            # 0.794176363371587

data_with_c = combine_data_C(data2, best_C)
sns.lmplot('X1', 'X2', hue='C', data=data_with_c, fit_reg=False)
plt.show()
1
2
3
4
5
6
7
8

4. 对比sklearn K-means

from sklearn.cluster import KMeans

sk_kmeans = KMeans(n_clusters=3)
sk_kmeans.fit(data2)
sk_C = sk_kmeans.predict(data2)

data_with_c = combine_data_C(data2, sk_C)
sns.lmplot('X1', 'X2', hue='C', data=data_with_c, fit_reg=False)
plt.show()
1
2
3
4
5
6
7
8
9

3. K-means用于图像压缩

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd
from skimage import io

pic = io.imread('../data/bird_small.png') / 255
io.imshow(pic)
1
2
3
4
5
6
7
8

pic.shape      # (128, 128, 3)

# serialize data
data = pic.reshape(128 * 128, 3)

def k_means(data, k, epoch=100, n_init=10):
    """do multiple random init and pick the best one to return
    Args:
        data (pd.DataFrame)
    Returns:
        (C, centroids, least_cost)
    """

    tries = np.array([_k_means_iter(data, k, epoch) for _ in range(n_init)])

    least_cost_idx = np.argmin(tries[:, -1])

    return tries[least_cost_idx]

# support fn --------------------------------
def combine_data_C(data, C):
    data_with_c = data.copy()
    data_with_c['C'] = C
    return data_with_c


# k-means fn --------------------------------
def random_init(data, k):
    """choose k sample from data set as init centroids
    Args:
        data: DataFrame
        k: int
    Returns:
        k samples: ndarray
    """
    return data.sample(k).iloc[:,:].values


def _find_your_cluster(x, centroids):
    """find the right cluster for x with respect to shortest distance
    Args:
        x: ndarray (n, ) -> n features
        centroids: ndarray (k, n)
    Returns:
        k: int
    """
    distances = np.apply_along_axis(func1d=np.linalg.norm,  # this give you l2 norm
                                    axis=1,
                                    arr=centroids - x)  # use ndarray's broadcast
    return np.argmin(distances)


def assign_cluster(data, centroids):
    """assign cluster for each node in data
    return C ndarray
    """
    return np.apply_along_axis(lambda x: _find_your_cluster(x, centroids),
                               axis=1,
                               arr=data.iloc[:,:].values)


def new_centroids(data, C):
    data_with_c = combine_data_C(data, C)

    return data_with_c.groupby('C', as_index=False).\
                       mean().\
                       sort_values(by='C').\
                       drop('C', axis=1).\
                       iloc[:,:].values


def cost(data, centroids, C):
    m = data.shape[0]

    expand_C_with_centroids = centroids[C]

    distances = np.apply_along_axis(func1d=np.linalg.norm,
                                    axis=1,
                                    arr=data.iloc[:,:].values - expand_C_with_centroids)
    return distances.sum() / m


def _k_means_iter(data, k, epoch=100, tol=0.0001):
    """one shot k-means
    with early break
    """
    centroids = random_init(data, k)
    cost_progress = []

    for i in range(epoch):
        print('running epoch {}'.format(i))

        C = assign_cluster(data, centroids)
        centroids = new_centroids(data, C)
        cost_progress.append(cost(data, centroids, C))

        if len(cost_progress) > 1:  # early break
            if (np.abs(cost_progress[-1] - cost_progress[-2])) / cost_progress[-1] < tol:
                break

    return C, centroids, cost_progress[-1]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

C, centroids, cost = k_means(pd.DataFrame(data), 16, epoch=10, n_init=3)
1

sklearn K-means

from sklearn.cluster import KMeans

model = KMeans(n_clusters=16, n_init=100)
model.fit(data)
centroids = model.cluster_centers_   # (16, 3)
C = model.predict(data)    # (16384,)
centroids[C].shape      # (16384, 3)

compressed_pic = centroids[C].reshape((128, 128, 3))

fig, ax = plt.subplots(1, 2)
ax[0].imshow(pic)
ax[1].imshow(compressed_pic)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14