(二)三维点云课程---KMeans介绍

KMeans作为一种三维点云中的聚类方法，由于原理简单，使用起来方便。现在就KMeans介绍它的原理

1.KMeans原理推导

推导KMeans原理之前，先说明一下下面用到的符号。三维数据集$x_1,x_2,...,x_N,x_n\in R^D$；聚类中心${\mu }_k,k=1,...,k$,k代表聚类的个数，KMeans中需要认为设定；二进制变量 γ n k ∈ { 0 , 1 } \gamma_{nk} \in \{0,1\} γnk​∈{0,1},其中1表示三维点中属于k类，0不属于k类。

其实KMeans的目的就是最小化损失函数$J$
$$J=\sum _{n=1}^N\sum _{k=1}^K{r}_{nk}{\Vert {\mathbf{x}}_n-{\mathbf{\mu }}_k\Vert }^2$$
那么怎么最小化损失函数呢，就要用到EM算法。

1.1E(expection)步

在E步就是固定${\mu }_k$,通过关心$r_{nk}$来减小$J$。因为三维点云数据是互相独立的，因此我们可以对与每一个点$n$单独进行优化，并且$r_{nk}$取值只有0，1。因此只需分配$n^{th}$数据指向最近的群集中心，这样的话$||x_n-{\mu }_k|{|}_2$才会最小。数学表达式如下
$$r_{nk}=\{\begin{array}{ll}1& \text{ if }k=\arg {\min }_j{\Vert {\mathbf{x}}_n-{\mathbf{\mu }}_j\Vert }^2\\ 0& \text{ otherwise }\end{array}$$
其中$argmi{n}_j||x_n-{\mu }_j|{|}_2$表示$j$取什么值时，会最小化$||x_n-{\mu }_j|{|}_2$。其实公式优点晦涩难懂，通过代码深入理解一下

#E-step过程
for idx,point in enumerate(data):  #data表示数据点 维数N*2
    #x_norm=np.linalg.norm(x, ord=None, axis=None, keepdims=False)
    #x: 表示矩阵（也可以是一维）,ord：范数类型默认为二范数，axis：处理类型，axis=1表示按行向量处理，求多个行向量的范数，axis=0表示按列向量处理，求多个列向量的范数，axis=None表示矩阵范数。
    diff=np.linalg.norm(old_center-point,axis=1)
    #numpy.argmin表示最小值在数组中所在的位置
    labels[np.argmin(diff)].append(idx)

具体分析一下

 diff=np.linalg.norm(old_center-point,axis=1)

其中old_center表示更新前聚类的中心点，就是上述的${\mu }_k$。通过np.linalg.norm()函数可以求解数据点距离中心点的欧式距离，维数为N*1。这也就是$||x_n-{\mu }_j|{|}_2$

 labels[np.argmin(diff)].append(idx)

其中numpy.argmin() 表示最小值在数组中所在的位置(labels在代码之前定义为labels=[[] for i in range(self.k_)] )， labels是一个二维数组，行数=聚类数，每一行里面存储的属于这个类的原始数据点的索引，那么此时$r_{nk}=1$,同理不属于这个类，$r_{nk}=0$。

1.2M(maximization)步

在M步,固定$r_{nk}$,$J$是一个关于${\mu }_k$的二次函数,那么通过求导来求解极小值。对$J$关于${\mu }_k$求解一阶导
$$\begin{gathered} J=\sum _{n=1}^N\sum _{k=1}^K{r}_{nk}{\Vert {\mathbf{x}}_n-{\mathbf{\mu }}_k\Vert }^2=\sum _{n=1}^N(r_{n1}||x_n-{\mu }_1|{|}^2+...+r_{nK}||x_n-{\mu }_K|{|}^2) \\ \frac{dJ}{d{\mu }_k}=2\sum _{n=1}^N{r}_{nk}\left(x_n-{\mu }_k\right)=0 \end{gathered}$$
推导出
$$u_k=\frac{{\sum }_n{r}_{nk}{x}_n}{{\sum }_n{r}_{nk}}$$
此时通过${\mu }_k$的公式，可以得出一个类的中心点等于属于这个类的数据点的平均值，编程中也是这么实现。

#M-step
for i in range(self.k_):
    points=data[labels[i],:]
    centers[i]=points.mean(axis=0)  #axis=0是对列进行计算

centers就是新得到的中心点。

---

2.KMeans的原理总结

• 人工给定一个K，这里有个关于起始中心点的技巧，可以选择原始数据的任意K个点作为起始中心点，加快算法效率
• E步: 每个点属于哪个类，n个点k个中心，就是求n个点到k个中心的距离
• M步: 知道每一个点被分配到了一个类之后，更新这个类中心点的位置
• 迭代E步和M步，直到算法收敛。

注意KMeans算法收敛的条件是①算法的中心点不在移动或者变化特别小，或者②$r_{nk} $不在变化

效果如下：(b)为E步，找到$r_{nk}$，©为M步，找到${\mu }_k$,其他的类似。

---

3.KMeans完整代码

import numpy as np
import random
import matplotlib.pyplot as plt

class K_Means(object):
    #KMeans初始化
    def __init__(self, n_clusters=2, tolerance=0.0001, max_iter=600):
        self.k_ = n_clusters                     #分组数
        self.tolerance_ = tolerance              #中心点误差,用于提前终止KMeans
        self.max_iter_ = max_iter                #最大迭代次数
        self.fitted=False                        #KMeans完成的标志
        self.centers=np.array([])                #中心点的存储结构

    #KMeans的核心算法
    def fit(self, data):
        # 随机的在0-data.shape[0]范围内选取self.k_个值作为中心点
        centers=data[random.sample(range(data.shape[0]),self.k_)]
        old_center=np.copy(centers)
        print("old_center:",old_center)
        #二维数组，行数=聚类数，每一行里面存储的属于这个类的原始数据点的索引
        labels=[[] for i in range(self.k_)]

        for iter_ in range(self.max_iter_):
            #E-step过程
            for idx,point in enumerate(data):
                #x_norm=np.linalg.norm(x, ord=None, axis=None, keepdims=False)
                #x: 表示矩阵（也可以是一维）,ord：范数类型默认为二范数，axis：处理类型，axis=1表示按行向量处理，求多个行向量的范数，axis=0表示按列向量处理，求多个列向量的范数，axis=None表示矩阵范数。
                diff=np.linalg.norm(old_center-point,axis=1)
                #numpy.argmin表示最小值在数组中所在的位置
                labels[np.argmin(diff)].append(idx)
            #M-step
            for i in range(self.k_):
                points=data[labels[i],:]
                centers[i]=points.mean(axis=0)  #axis=0是对列进行计算
            #提前终止KMeans的步骤
            if np.sum(np.abs(centers-old_center))<self.tolerance_*self.k_:
                break
            #更新KMeans的中心点
            old_center=np.copy(centers)  #copy为深拷贝，即拷贝前后完全不相干
        # True表示KMeans已结束
        self.fitted=True
        # 存储KMeans最终的中心点
        self.centers=centers

    #通过最终的中心点，根据欧式距离预测每个点的分类
    def predict(self, p_datas):
        result = []
        if not self.fitted:
            print ("Unfitted")
            return result
        for point in p_datas:
            diff=np.linalg.norm(self.centers-point,axis=1)
            result.append(np.argmin(diff))
        return result

## 读取数据的函数
def read_txt(path):
    filename = path  # txt文件和当前脚本在同一目录下，所以不用写具体路径
    pos = []
    Efield = []
    with open(filename, 'r') as file_to_read:
        while True:
            lines = file_to_read.readline()  # 整行读取数据
            if not lines:
                break
                pass
            p_tmp, E_tmp = [float(i) for i in lines.split(",")]  # 将整行数据分割处理，如果分割符是空格，括号里就不用传入参数，如果是逗号， 则传入‘，'字符。
            pos.append(p_tmp)  # 添加新读取的数据
            Efield.append(E_tmp)
            pass
    pos = np.array(pos)  # 将数据从list类型转换为array类型。
    Efield = np.array(Efield)
    x=np.vstack(pos)    #将pos按照pos.shape[1]*pos.shape[0]进行排列  N*1
    y=np.vstack(Efield) #N*1
    """
    >>> a=np.array([[1,2,3],[4,5,6]])
    >>> b=np.array([[11,21,31],[7,8,9]])
    >>> np.concatenate((a,b),axis=0)
    array([[ 1,  2,  3],
           [ 4,  5,  6],
           [11, 21, 31],
           [ 7,  8,  9]])
    >>> np.concatenate((a,b),axis=1)  #axis=1表示对应行的数组进行拼接
    array([[ 1,  2,  3, 11, 21, 31],
           [ 4,  5,  6,  7,  8,  9]])
    """
    X=np.concatenate((x,y),axis=1) #axis=1表示对应行的数组进行拼接，即行数不变，增加列数,N*2
    return X


##初始化画布为4行2列
fig, ax = plt.subplots(4, 2)

## 可视化函数
## 输入：x 原始数的路径
#        i 画布的行
#        j 画布的列
#        n_clusters 聚类的个数，认为指定
def show_result(x,i,j,n_clusters):
    ax[i][j].scatter(x[:, 0], x[:, 1], s=20, c="b", marker='o')      #蓝色的"o"表示原始数据
    ax[i][j].set_xlabel('x')
    ax[i][j].set_ylabel('y')

    k_means = K_Means(n_clusters)
    k_means.fit(x)
    cat = k_means.predict(x)
    list_max = max(cat)
    if list_max == 1:   #两类
        for idx, value in enumerate(cat):
            if value == 0:
                ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="r", marker='*')  # 红色的"*"表示一类
            elif value == 1:
                ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="g", marker='+')  # 绿色的"+"表示一类
    elif list_max == 2:  #三类
        for idx, value in enumerate(cat):
            if value == 0:
                ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="r", marker='*')  # 红色的"*"表示一类
            elif value == 1:
                ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="g", marker='+')  # 绿色的"+"表示一类
            elif value == 2:
                ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="y", marker='^')  # 黄色的"^"表示一类
    ax[i][j+1].set_xlabel('x')
    ax[i][j+1].set_ylabel('y')


if __name__ == '__main__':
    x_circle=read_txt("circle.txt")
    print(x_circle.shape)
    show_result(x_circle,0,0,2)

    x_moons=read_txt("moons.txt")
    show_result(x_moons,1,0,2)

    x_blobs=read_txt("varied.txt")
    show_result(x_blobs,2,0,3)

    x_aniso=read_txt("aniso.txt")
    show_result(x_aniso,3,0,3)

    plt.show()

仿真结果如下,左边为原始数据，右边是经过KMeans聚类后的结果

数据集下载：
aniso.txt: https://url23.ctfile.com/f/22628623-518637820-025a10 （访问密码：7875）
blobs.txt: https://url23.ctfile.com/f/22628623-518637821-76c706 （访问密码：7875）
circle.txt: https://url23.ctfile.com/f/22628623-518637822-97aad4 （访问密码：7875）
moons.txt: https://url23.ctfile.com/f/22628623-518637823-189bf8 （访问密码：7875）
varied.txt: https://url23.ctfile.com/f/22628623-518637824-7d006a （访问密码：7875）

---

4.KMeans的限制

1.对噪声数据比较敏感，故提出了K-Medois（中心点）的思想，参见附录
2.k值需要认为指定，并不是特别的智能
3.$r_{nk}$非1即0，缺乏不确定的指标，如果数据点位于中心之间的边界线周围，如下图。对于这个问题，参见下文的GMM模型讲解，它就解决了这个问题。

---

附录

KMedios算法的基本解释

由于数据存在噪声点，可能会拉偏中心点的位置，那么提出了KMedios的方法
$$\begin{gathered} K-Medios\tilde{J}=\sum _{n=1}^N\sum _{k=1}^K{r}_{nk}\mathcal{V}\left({\mathbf{x}}_n,{\mathbf{\mu }}_k\right) \\ K-MeansJ=\sum _{n=1}^N\sum _{k=1}^K{r}_{nk}{\Vert {\mathbf{x}}_n-{\mathbf{\mu }}_k\Vert }^2 \end{gathered}$$
其中$\mathcal{V}\left({\mathbf{x}}_n,{\mathbf{\mu }}_k\right)$可能是不可导的。那么怎么求解呢

1）E步

仅仅通过将每个点离其他点的距离之和代替欧式距离

2）M步

求解一个点到其他点的距离之和最小值作为更新后的中心点，因此KMedios中心点一定是某一个数据点，时间复杂度为$O(N_{nk}^2)$