聚类方法小全（二）

最新推荐文章于 2024-08-09 06:54:59 发布

花咲溶溶月

最新推荐文章于 2024-08-09 06:54:59 发布

阅读量96

点赞数 1

分类专栏： Python 点云文章标签：聚类机器学习支持向量机

本文链接：https://blog.csdn.net/apple_51522252/article/details/132840184

版权

Python 同时被 2 个专栏收录

6 篇文章 0 订阅

订阅专栏

点云

3 篇文章 0 订阅

订阅专栏

4 谱聚类

谱聚类是一种利用无向图进行聚类的方法。实际上，就是让相似的点之间进行连线，从而形成一个无向图进行分析的方法，具体如下图所示。

建一个图 $G = (V, E)$ ，存放方式用的是邻接矩阵 $w$ 。如果两点不相连，则 $w_{ij} = 0$ ，否则 $w_{ij} = d(i, j)$ ， $d$ 函数的基本特性就是离得越近，函数值越大，可以是反比例，也可以是用很大的数去减掉距离，总之只要满足要求就可以。特别地， $w_{ii} = 0$ 。

4.1 建图方法

我们首先要知道，每一个点我们都设为是一个节点。那么建图方法主要分为以下三种：

** $\epsilon$ -邻居建图。**每一个节点而言，在一定半径范围内的点都连边，令 $w_{ij} = d(v_i, v_j)$
**k-邻居建图。**实际上就是KNN，K邻近的点就连边，令 $w_{ij} = d(v_i, v_j)$
- 需要满足 $v_i$ 是 $v_j$ 的K邻近点，同时 $v_j$ 也是 $v_i$ 的K邻近点，才可以连边。
- 满足 $v_i$ 是 $v_j$ 的K邻近点或者 $v_j$ 是 $v_i$ 的K邻近点其中之一就连边。
**全连接。**顾名思义，两点之间全部连边。

4.2 一些数学概念

**度矩阵 $D$ 。**这是一个对角矩阵，其中 $d_{ii}$ 就等于 $w$ 矩阵第 $i$ 行的和，也就是 $d_{ii} = \sum_{j = 1}^nw_{ij}$ 。它实际上就是说这个点连向其他点的权重和，也就是总共出去（进来）了多少。
拉普拉斯矩阵 $L$ （未标准化的）。 $L = D - W$
标准化拉普拉斯矩阵。
- $L_{sym} = D^{-1/2}WD^{-1/2}$
- $L_{rw} = D^{-1}L = I - D^{-1}W$

4.3 无标准化的谱聚类步骤

首先建图，生成邻接矩阵 $W\in\R^{n*n}$
计算未标准化的拉普拉斯矩阵 $L = D - W$
计算拉普拉斯矩阵 $L$ 的 $k$ （聚类数量）个最小的特征向量 $v_1, v_2, ..., v_k$
这些特征向量作为列向量拼接成一个新的矩阵 $V\in\R^{n * k}$
之后，将矩阵 $V$ 中的每一行作为向量进行分析，即 $y_i$ 为矩阵 $V$ 的第 $i$ 行
使用K-Means算法将 $\{y_i \in \R^k\}$ 进行聚类，分为 $C_1, C_2, ..., C_k$
最终的聚类结果就是 $A_1, A_2, ..., A_k$ ，其中 $A_i = \{j|y_j \in C_i\}$ 。实际上就是说一个类里面有 $y$ 几就说明点几属于这个类。

4.4 标准化的谱聚类步骤

和4.3一样，只不过第二步和第三步有一点区别：

首先建图，生成邻接矩阵 $W\in\R^{n*n}$
计算标准化的拉普拉斯矩阵 $L' = L_{rw}$
计算拉普拉斯矩阵 $L^{'}$ 的 $k$ （聚类数量）个最小的特征向量 $v_1, v_2, ..., v_k$
这些特征向量作为列向量拼接成一个新的矩阵 $V\in\R^{n * k}$
之后，将矩阵 $V$ 中的每一行作为向量进行分析，即 $y_i$ 为矩阵 $V$ 的第 $i$ 行
使用K-Means算法将 $\{y_i \in \R^k\}$ 进行聚类，分为 $C_1, C_2, ..., C_k$
最终的聚类结果就是 $A_1, A_2, ..., A_k$ ，其中 $A_i = \{j|y_j \in C_i\}$ 。实际上就是说一个类里面有 $y$ 几就说明点几属于这个类。

4.5 归一化和不归一化（标准化）的区别

不标准化的结果更倾向于类之间点的数量接近

标准化的结果更倾向于类里面点的密度更加接近

具体如下图：

4.6 k的选定

谱聚类的好处就是可以确定 $k$ 的值。

拉普拉斯矩阵在计算特征向量的时候，顺便计算一下特征值。如果特征值有突变，那就可以确定 $k$ 为突变前的值了。

比如说上面这张图所示的情况，4和5之间有突变，所以最终可以确定 $k = 4$ 。

其实就是找最大的 $\Delta k = |\lambda_k - \lambda_{k - 1}|$

4.7 原理

简单来说，将各个点之间连线之后，要对其进行分类，实际上就是进行图的最小割操作。在这其中，最重要的就是拉普拉斯矩阵 $L$ 。我们就以未经过标准化的拉普拉斯矩阵为例进行讲解。

4.7.1 拉普拉斯矩阵特点

拉普拉斯矩阵 $L$ 具有以下特点：

对于每一个向量 $\in \R^n$ ，我们有 $f^TLf = \frac{1}{2}\sum_{i = 1}^n\sum_{i = 1}^nw_{ij}(f_i - f_j)^2$
拉普拉斯矩阵 $L$ 是对称矩阵且是半正定的
拉普拉斯矩阵 $L$ 的最小特征值为 $0$ ，对应的特征向量为矩阵 $1$
拉普拉斯矩阵 $L$ 有 $n$ 个非负特征值，并且 $\lambda_1 \leq \lambda_2 \leq ... \leq \lambda_n$

下面来对这些特点进行证明。

对于每一个向量 $\in \R^n$ ，我们有 $f^TLf = \frac{1}{2}\sum_{i = 1}^n\sum_{i = 1}^nw_{ij}(f_i - f_j)^2$ ：
$\begin{align*} f^TLf &= f^TDf - f^TWf = \sum_{i = 1}^nf_i^2d_i -\sum_{i = 1}^n\sum_{j = 1}^nf_if_jw_{ij} \\ &= \frac{1}{2}(\sum_{i = 1}^nd_if_i^2 - 2\sum_{i = 1}^n\sum_{j = 1}^nf_if_jw_{ij} + \sum_{j = 1}^nd_jf_j^2) \\ &= \frac{1}{2}(\sum_{i = 1}^n\sum_{j = 1}^nw_{ij}f_i^2 - \sum_{i = 1}^n\sum_{j = 1}^nf_if_jw_{ij} + \sum_{i = 1}^n\sum_{j = 1}^nw_{ji}f_j^2)\\ &= \frac{1}{2}\sum_{i = 1}^n\sum_{i = 1}^nw_{ij}(f_i - f_j)^2 \end{align*}$
拉普拉斯矩阵 $L$ 是对称矩阵且是半正定的：

$L = D - W$ ，且 $D$ 和 $W$ 矩阵为对称的，因此矩阵 $L$ 也是对称的
$f^TLf \geq 0$ ，因为平方肯定大于等于 $0$ ，且 $w_{ij} \geq 0$

拉普拉斯矩阵 $L$ 的最小特征值为 $0$ ，对应的特征向量为矩阵 $1$ ：

$[...,d_if_i - \sum_{j = 1}^nw_{ij}f_j,...]^T = 0·f, f = 1$

为什么是 $L f$ 呢，因为我们知道有 $\lambda f$ ，求解之后 $\lambda$ 就是特征值， $f$ 就是特征向量。而我们注意到， $d_i = \sum_{j = 1}^nw_{ij}$ ，那么实际上 $L f$ 中的每一个值 $d_if_i - \sum_{j = 1}^nw_{ij}f_j$ ，如果 $f$ 是常量，那么这个值无论 $i$ 和 $j$ 取几，它都是 $0$ ，所以最小特征值肯定是 $0$ ， $f = 1$

拉普拉斯矩阵 $L$ 有 $n$ 个非负特征值，并且 $\lambda_1 \leq \lambda_2 \leq ... \leq \lambda_n$

这个直接就是上面两条的直接结果。

4.7.2 另一个关于拉普拉斯矩阵的小特点

对于一个无向图 $G$ ，如果它的每一条边都是非负数，那么它有多少个特征值为 $0$ ，取决于它有多少个连通域。特征值为 $0$ 对应的特征向量为常数向量 $1$ 。如下图所示：

可问题是，你这也没有常熟向量 $1$ 啊？别急，看下面：

$\lambda x$

$Lx_1 + Lx_2 = \lambda_1x_1 + \lambda_2x_2$

$\lambda_1 = \lambda_2 = 0$

$L(x_1 + x_2) = \lambda(x_1 + x_2)$

因此， $x_1 + x_2$ 同样也是一个特征向量。

证明：

首先注意只有一个连通域的情况，假设 $f$ 是特征值为 $0$ 所对应的特征向量，那么就有 $f^TLf = \frac{1}{2}\sum_{i = 1}^n\sum_{i = 1}^nw_{ij}(f_i - f_j)^2 = 0 · f^T = 0$

唯一满足的情况就是 $f_i = f_j$ ，所以对于任意 $i, j$ 来说，都有 $f_i = f_j$ .

举例来说，有

而对于有多个连通域的情况来说，其实也是一样的。我们通过交换这个大矩阵的行和列，让每一个连通域的拉普拉斯矩阵凑在一起，也就是变成：
$\begin{bmatrix} L_1 & & & \\ & L_2 & & \\ & & \ddots & \\ & & & L_k \end{bmatrix}$
这样，它的特征值为 $0$ 所对应的特征向量就变成了
$\begin{bmatrix} \bold{1}_1 & 0 & \cdots & 0\\ 0 & \bold{1}_2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \bold{1}_k \end{bmatrix} \in\R^{n * k}$

4.7.3 图的割

我们定义将图分成两个连通域，它的割就是$cut(A, B) = \sum_{i \in A, j \in B}w_{ij} \$

那么让如果是分成多个连通域，那么它的割就是 $cut(A_1,..., A_k) = \sum_{i = 1}^k cut(A_i, \bar A_i)$

对于非标准化的谱聚类，我们有 $RatioCut(A_1,..., A_k) = \sum_{i = 1}^k\frac{cut(A_i, \bar A_i)}{|A_i|}$

对于标准化的谱聚类，我们有 $Ncut(A_1,..., A_k) = \sum_{i = 1}^k\frac{cut(A_i, \bar A_i)}{vol(A_i)}$

首先考虑只有两类的情况：

我们要求最小割，也就是如下的优化问题：

$\min_{A \subset V} RatioCut(A, \bar{A}) = \min_{A \subset V}(\frac{cut(A, \bar{A})}{|A|} + \frac{cut(\bar{A}, A)}{|\bar A|})\$

那么我们构造一个向量 $[f_1,...,f_n]^T \in \R^n$

$f_i=\left\{ \begin{array}{rcl} \sqrt{|\bar A|/|A|} & & {v_i \in A}\\ \sqrt{-|A|/|\bar A|} & & {v_i \in \bar{A}} \end{array} \right.$

通过我们构造的向量的这个特性，我们就可以有：

$f_i=\left\{ \begin{array}{rcl} \ v_i \in A & & {f_i \geq 0}\\ \ v_i \in \bar A & & {f_i < 0} \end{array} \right.$

然后我们利用之前所得到的结论进行推论：
$\begin{align*} f^TLf &= \frac{1}{2}\sum_{i,j = 1}^nw_{ij}(f_i - f_j)^2 \\ &= \frac{1}{2}\sum_{i \in A, j \in \bar A} w_{ij}(\frac{|\bar A|}{|A|} + \frac{|A|}{|\bar A|})^2 + \frac{1}{2}\sum_{i \in \bar A, j \in A} w_{ij}(-\frac{|\bar A|}{|A|} - \frac{|A|}{|\bar A|})^2 \\ &= cut(A, \bar A)(\frac{|\bar A|}{|A|} + \frac{|A|}{|\bar A|} + 2) \\ &= cut(A, \bar A)(\frac{|A| + |\bar A|}{|A|} + \frac{|A| + |\bar A|}{|\bar A|})\\ &= |V|·RatioCut(A, \bar A) \end{align*}$
另外，我们所设计的 $f$ 还有两个特点：

$f$ 与常数向量 $1$ 垂直
- 证明：$\sum_{i = 1} ^ n f_i = \sum_{i \in A}\frac{|\bar A|}{|A|} - \sum_{i \in \bar A}\frac{|A|}{|\bar A|} = |A|\frac{|\bar A|}{|A|} - |\bar A|\frac{|A|}{|\bar A|} = 0\$
$\sqrt{n}$
- 证明： $||f||^2 = \sum_{i = 1} ^n f_i^2 = |A|\frac{|\bar A|}{|A|} - |\bar A|\frac{|A|}{|\bar A|} = |A| + |\bar A| = n$

那么实际上，我们就是要求解一个优化问题：

$\min _{A \sub V} f^TLf, s.t.,f \perp 1, ||f|| = \sqrt n$

$f_i=\left\{ \begin{array}{rcl} \sqrt{|\bar A|/|A|} & & {v_i \in A}\\ \sqrt{-|A|/|\bar A|} & & {v_i \in \bar{A}} \end{array} \right.$

但是，第二个条件也就是我们构造的向量是很难满足的，因此我们采取近似值，抛弃第二个条件，只要求

$\min _{A \sub V} f^TLf, s.t.,f \perp 1, ||f|| = \sqrt n$

那这个很明显就可以用到瑞利商定理（具体可以看PCA章节），但别忘了我们有两个约束条件在，其中第一个 $\perp 1$ 这个，由于最小的特征值对应的特征向量一定是常数向量，所以就找特征值第二小的。

如果不近似，那么结果就是上图左边的情况，严谨地就只有两个值。但由于我们近似了，就会出现右边的情况。那应该怎么分开呢？跑K-Means就好了！

**对于有两类以上的情况：**也是同样的道理，只是 $f$ 向量变为了 $h$ 矩阵，其他的都一样。

4.8 代码

import numpy as np
import open3d as o3d
from sklearn.neighbors import kneighbors_graph
from sklearn import datasets

n_sample = 500
noisy_moon = datasets.make_moons(n_samples=n_sample, noise=.05)

def dis(a, b) :
    return (a[0] - b[0]) * (a[0] - b[0]) + (a[1] - b[1]) * (a[1] - b[1])

def kmeans(data, k, tolerance, maxIter) :
    idx = np.random.choice(len(data), k, replace=False)
    centroids = data[idx]
    cluster = np.zeros(len(data))
    dist = np.zeros(len(data))
    for i in range(maxIter) :
        for j in range(len(data)) :
            cluster[j] = -1
            dist[j] = 1e10
        for j in range(len(data)) :
            for k in range(len(centroids)) :
                tmp = dis(data[j], centroids[k])
                if tmp < dist[j] :
                    dist[j] = tmp
                    cluster[j] = k

        old_centroids = np.copy(centroids)
        for j in range(len(centroids)) :
            dataTmp = data[cluster == j]
            if len(dataTmp) > 0 :
                centroids[j][0] = np.mean(dataTmp[:, 0])
                centroids[j][1] = np.mean(dataTmp[:, 1])
            else :
                centroids[j] = data[np.random.choice(len(data), 1, replace=False)]

        distortion = np.linalg.norm(centroids - old_centroids)
        if distortion < tolerance :
            break

    return cluster, centroids

def spectral(data, k, maxIter = 100, n_neighbour = 20) :
    weight = kneighbors_graph(data, n_neighbour, mode='connectivity', include_self=False)
    weight = 0.5 * (weight + weight.T)
    W = weight.toarray()
    D = np.diag(np.sum(W, axis=1))
    L = D - W

    w, v = np.linalg.eig(L)
    idx = np.argsort(w)[:k]
    vec = v[:, idx]
    cluster, centroids = kmeans(vec, k, 0.0001, maxIter)
    return cluster

datas = np.zeros([len(noisy_moon[0]), 3])
for i in range(len(noisy_moon[0])) :
    datas[i][0] = noisy_moon[0][i][0]
    datas[i][1] = noisy_moon[0][i][1]
    datas[i][2] = 0
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(datas)

cluster = spectral(datas, 2)
color = np.zeros([len(noisy_moon[0]), 3])
for i in range(len(noisy_moon[0])) :
    if cluster[i] == 0 :
        color[i][0] = 1
    else :
        color[i][1] = 1

pcd.colors = o3d.utility.Vector3dVector(color)

o3d.visualization.draw_geometries([pcd])

5 Mean Shift

5.1 使用方法

随机选择一个半径为 $r$ 的圆
将圆的圆心移动到里面点的平均值上
重复步骤二，直到圆心不再移动
重复上面三个步骤，并去掉重叠的圆（若重叠，留着包含点最多的圆）
每个点选择最近的圆心作为它的分类

这个方法主打一个简单，所需要设计的参数只有半径 $r$ ，也不用提前知道需要几类

6 DBSCAN

6.1 使用方法

随机选择一个没有被访问过的点 $p$ ，然后找到它半径为 $r$ 范围内的邻居点
判断邻居点的数量是否大于 $min\_samples$
- 如果大于，那么我们就认为 $p$ 是一个核心点。此时，我们创建一个聚类 $C$ ，标记点 $p$ 为被访问过的，并继续步骤3
- 如果小于，那我们就标记点 $p$ 为访问过的，并标记它为噪声点
将这些邻居都标记为 $C$ 类的。如果这个点是核心点，那我们就把它作为新的点 $p$ ，重复该步骤
完成后从数据库移除聚类 $C$ ，然后返回步骤一
直到所有的点都被标记

所以这里面需要设计的参数为半径 $r$ 和 $min\_sample$ ，而且好处就是类别数量不需要提前定好。

6.2 各种点

这就是一个 $r = Eps, min\_sample = 4$ 的例子。可以看到，蓝色的点没有邻居点，因此就是噪声。

红色的点则会不停的去访问其他的点，由于这个图实在是太过清晰了，不再赘述。

6.3 代码

import numpy as np
import open3d as o3d
from sklearn.neighbors import kneighbors_graph
from sklearn import datasets
from scipy.spatial import cKDTree

n_sample = 500
noisy_moon = datasets.make_moons(n_samples=n_sample, noise=.05)

class DBSCAN() :
    def __init__(self, data, r, min_sample):
        self.r = r
        self.min_sample = min_sample
        self.data = data
        self.n, self.k = self.data.shape
        self.vis = np.zeros(self.n)
        self.kdtree = cKDTree(self.data)
        self.now = 0
        self.ans = np.zeros(self.n)

    def check(self, p):
        self.vis[p] = True
        index = self.kdtree.query_ball_point(self.data[p], self.r)
        if len(index) < self.min_sample :
            self.ans[p] = 0
            self.vis[p] = True
            return
        else :
            self.ans[p] = self.now
            for point in index :
                if not self.vis[point] :
                    self.check(point)

    def fit(self):
        while(1):
            if np.sum(self.vis) == self.n :
                break
            self.now += 1
            tmp = np.where(self.vis == False)[0]
            self.check(np.random.choice(tmp, 1, replace=False)[0])
        return

    def predict(self):
        return self.ans

dbscan = DBSCAN(noisy_moon[0], 0.1, 5)
dbscan.fit()
cluster = dbscan.predict()
print(cluster)

datas = np.zeros([len(noisy_moon[0]), 3])
for i in range(len(noisy_moon[0])) :
    datas[i][0] = noisy_moon[0][i][0]
    datas[i][1] = noisy_moon[0][i][1]
    datas[i][2] = 0
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(datas)

color = np.zeros([len(noisy_moon[0]), 3])
for i in range(len(noisy_moon[0])) :
    if cluster[i] == 0 :
        color[i][0] = 1
    elif cluster[i] == 1:
        color[i][1] = 1
    else :
        color[i][2] = 1

pcd.colors = o3d.utility.Vector3dVector(color)

o3d.visualization.draw_geometries([pcd])