Bo

density_peaks_clustering_algorithm

发表于 更新于 分类于 阅读次数： Valine： 本文字数： 9.9k 阅读时长 ≈ 9 分钟

处理数据时学习了一种新的聚类方法–密度峰值聚类。记录一下学习过程。

密度峰值聚类说白了就是4步：

第1步，求每个点的密度rho。点的密度就是，以点为中心，以dc为半径，画一个小圆圈，数数里面几个点，圆圈中点的个数就是点的密度。（还可以用高斯核密度求点的密度，求出来的密度是连续型的）

第2步，计算每个点的delta。假设有一个点x，求比点x的密度大的且距离点x最近的那个点y，那么点x与点y之间的距离，就是点x的delta，就这样遍历所有点，把每个点的delta都求出来（注意，delta是距离，谁和谁的距离？x和y的距离，y是谁？y就是比x密度大，且距离x最近的那个点，要满足两个条件，密度比x大且距离最近）

第3步，每个点的密度rho和delta都求出来了，以rho为横坐标，delta为纵坐标，画个二维图，图中右上角的那几个点就是聚类中心，也就是rho和delta都很大的那几个点。（为什么？因为聚类中心有个特点，密度很大，且与密度比它大的点的距离也很大）

第4步，找到聚类中心了，就可以扩展聚类簇了，按照rho从大到小的顺序进行聚类扩展。

就是这么简单！这个算法有很多缺点，比如不适合高维，第1步中的半径dc难以选择，效率太低等等，因此成百上千篇的论文都在优化这个算法。

1、背景介绍

　　密度峰值算法（Clustering by fast search and find of density peaks）由Alex Rodriguez和Alessandro Laio于2014年提出，并将论文发表在Science上。Science上的这篇文章《Clustering by fast search and find of density peaks》主要讲的是一种基于密度的聚类方法，基于密度的聚类方法的主要思想是寻找被低密度区域分离的高密度区域。 密度峰值算法(DPCA)基于这样的假设：（1）类簇中心点的密度大于周围邻居点的密度；（2）类簇中心点与更高密度点之间的距离相对较大。因此，DPCA主要有两个需要计算的量：第一，局部密度；第二，与高密度点之间的距离。

2、局部密度

　　数据对象x{_{i}}的局部密度img定义为：img

其中，dist_{cutoff}表示截断距离img，这个公式的含义是说找到与第img个数据点之间的距离小于截断距离img的数据点的个数，并将其作为第i个数据点真的密度。

img

3、定义聚类中心距离

　　密度峰聚类算法的巧妙之处：就是在于聚类中心距离 δi的选定。根据局部密度的定义，我们可以计算出上图中每个点的密度，依照密度确定聚类中心距离 δi。

  1. 首先将每个点的密度从大到小排列： ρi > ρj > ρk > ….；密度最大的点的聚类中心距离与其他点的聚类中心距离的确定方法是不一样的;

  2. 先确定密度最大的点的聚类中心距离–i点是密度最大的点，它的聚类中心距离δiδi等于与i点最远的那个点n到点i的直线距离 d(i,n)；

  3. 再确定其他点的聚类中心距离——其他点的聚类中心距离是等于在密度大于该点的点集合中，与该点距离最小的的那个距离。例如i、j、k的密度都比n点的密度大，且j点离n点最近，则n点的聚类中心距离等于d(j,n);

  4. 依次确定所有的聚类中心距离δ
    img

4、聚类效果

　　将所有点的聚类中心密度都统计出来后，将其值按 δi和pi作为坐标轴作图可以得到如图所示结果。可以看到图中1,10两个聚类中心同时远离坐标轴。普通点则是靠近p轴，异常点靠近 δ轴。

img

5、代码

​ 根据以上方法试着写出代码

1、计算距离矩阵

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 计算距离矩阵
def cal_dis(data):
    distance_all = np.random.rand(len(data), len(data))

    for i in range(len(data)):
        for n in range(len(data)):
            distance_all[i][n] = np.sqrt(np.square(abs(data[i][0]) - abs(data[n][0])) +
                                         np.square(abs(data[i][1]) - abs(data[n][1])) +
                                         np.square(data[i][2] - data[n][2]))
            # np.square(data[i][3]-data[n][3])+
            # np.square(data[i][4]-data[n][4])+
            # np.square(data[i][5]-data[n][5])+
            # np.square(data[i][6]-data[n][6])+
            # np.square(data[i][7]-data[n][7]))
    # print('距离矩阵：\n', distance_all, '\n')
    return distance_all

根据数据的个数，建立n*n阶矩阵，在这里距离值是数据内各维度数据之差的平方，维度越高可能精度越不准确，还有很多论文提出有关高维的密度峰值聚类方法，可以加以参考，这里就用了三个维度主要数据。

2、计算密度数组

1
2
3
4
5
6
7
8
9
10
11
def cal_den(distance_all, min_distance):
    point_density = np.random.rand(distance_all.shape[0])
    for i in range(distance_all.shape[0]):
        tmp = 0
        for n in range(distance_all.shape[0]):
            if 0 < distance_all[i][n] < min_distance:
                tmp = tmp + 1
        point_density[i] = tmp
    # print('点密度数组：\n', point_density, '\n')

    return point_density

计算密度数组需要输入两个参数，距离矩阵和最小距离，这里的最小距离是以这个点为中心，以多长距离为半径画圆，在这个范围内的数据的个数就是这个中心点的密度。我参考的一个最小距离的值是在距离矩阵内最大值和最小值的平均值。

3、计算密度最大的点的聚类中心距离

1
2
3
4
5
6
7
8
# 计算密度最大的点的聚类中心距离
def get_max_dis(distance_all, point_density):
    point_density = point_density.tolist()
    a = max(point_density)  # 寻找密度最大的点
    b = point_density.index(a)  # 寻找密度最大的索引
    c = max(distance_all[b])  # 密度最大点的最大中心距离

    return c

之前说先确定密度最大的点的最大中心距离，因为这个点没有比它密度还要大的点，所以将这个点另行计算。

4、到点密度大于自身最近点的距离

1
2
3
4
5
6
7
8
9
10
# 到点密度大于自身最近点的距离
def get_min_dis(point, i, distance_all, point_density):
    min_distance = []
    if point:
        for j in point:
            min_distance.append(distance_all[i][j])
        return min(min_distance)
    else:
        max_distance = get_max_dis(distance_all, point_density)
        return max_distance

point参数来自下方的传参，就是一个列表，列表中存着密度比第i个点大的点的数据。如果point是个空列表，说明这个点是密度最大点，所以用get_max_dis方法返回这个点的聚类中心距离。

5、计算各点到聚类中心的距离

1
2
3
4
5
6
7
8
9
10
11
# 计算各点到聚类中心的距离(dist)
def get_each_dis(distance_all, point_density):
    nn = []
    for i in range(len(point_density)):
        aa = []
        for j in range(len(point_density)):
            if point_density[i] < point_density[j]:
                aa.append(j)
        ll = get_min_dis(aa, i, distance_all, point_density)
        nn.append(ll)
    return nn

这里就是根据密度数组，找出每一个数据密度比自身大的点，把这些点加入到一个列表中，这里aa就是存放比自己密度大点的列表，如果自身是密度最大点，那么传到get_min_dis中的aa就是一个空的列表，最后将每个点到聚类中心的距离返回。

6、根据特征分类

1
2
3
4
5
6
7
8
9
10
11
# 分簇，根据需要划分的特征结果数找到数据密度中最大的n个点
def clustering(point_density, delta):
    mul1 = (point_density * delta).tolist()
    mul2 = mul1.copy()
    mul2.sort(reverse=True)
    v = np.random.rand(feature_num)
    for i in range(feature_num):
        for j in range(len(mul1)):
            if mul1[j] == mul2[i]:
                v[i] = mul1.index(mul1[j])
    return v

寻找最大的feature_num个密度中心点，寻找的方法是$\rho$和$\delta$的乘积，$\rho$越大，说明这个点周围的点越多，$\delta$越大，说明这个点离比它密度大的点距离越远，所以根据两者乘积，就可以基本上得到几个分类中心，我这里将数据分为了四类，也就是feature_num=4

7、输出每个点的分类

1
2
3
4
5
6
7
8
9
# 归类，判断所有点与V中点的距离，距离最近的划为一类，输出每个数据的类别
def classify(distance_all, v):
    c = np.random.rand(distance_all.shape[0])
    for i in range(distance_all.shape[0]):
        dist = []
        for j in range(len(v)):
            dist.append(distance_all[i][int(v[j])])
        c[i] = dist.index(min(dist)) + 1
    return c

由上面得到了分出来的四类中心点，这块代码就是将除了这四个中心点外的其他点，判断它们离哪一个中心点更近，将它与离其最近的中心点划为一类。

8、计算最大最小距离

1
2
3
4
5
6
7
8
9
10
11
12
13
# 计算最大最小距离
def max_min_dis(distance_all):
    max_dis = 0
    min_dis = 999
    for i in range(distance_all.shape[0]):
        for j in range(distance_all.shape[0]):
            if distance_all[i][j] == 0:
                continue
            if distance_all[i][j] > max_dis:
                max_dis = distance_all[i][j]
            if distance_all[i][j] < min_dis:
                min_dis = distance_all[i][j]
    return max_dis, min_dis

这里的代码可有可无，就是计算距离矩阵中距离的最大差，方便计算第二步中的密度数组。

源代码

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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
import matplotlib.pyplot as plt
import numpy as np

feature_num = 4     # 根据特征分成四组


# 读取文件将数据存到一个列表
def read_file(file_name):
    file = open(file_name, 'r')
    array = []
    for line in file:
        array.append([int(x) for x in line.split()])  # int将字符串转换成整型，方便后续计算
    file.close()
    return array


# 将距离矩阵归一化  可以直接遍历除以行最大值，我这里的方法有点画蛇添足
def norm(data):
    data = np.array(data).astype(float)
    max_array = np.zeros(data.shape[1])
    for i in range(data.shape[0]):
        for j in range(data.shape[1]):
            if abs(data[i][j]) > max_array[j]:
                max_array[j] = abs(data[i][j])
    for i in range(data.shape[0]):
        for j in range(data.shape[1]):
            data[i][j] = data[i][j] / max_array[j]
    return data


# 计算距离矩阵和密度数组
def cal_dis(data):
    distance_all = np.random.rand(len(data), len(data))

    for i in range(len(data)):
        for n in range(len(data)):
            distance_all[i][n] = np.sqrt(np.square(abs(data[i][0]) - abs(data[n][0])) +
                                         np.square(abs(data[i][1]) - abs(data[n][1])) +
                                         np.square(data[i][2] - data[n][2]))
            # np.square(data[i][3]-data[n][3])+
            # np.square(data[i][4]-data[n][4])+
            # np.square(data[i][5]-data[n][5])+
            # np.square(data[i][6]-data[n][6])+
            # np.square(data[i][7]-data[n][7]))
    # print('距离矩阵：\n', distance_all, '\n')
    return distance_all


def cal_den(distance_all, min_distance):
    point_density = np.random.rand(distance_all.shape[0])
    for i in range(distance_all.shape[0]):
        tmp = 0
        for n in range(distance_all.shape[0]):
            if 0 < distance_all[i][n] < min_distance:
                tmp = tmp + 1
        point_density[i] = tmp
    # print('点密度数组：\n', point_density, '\n')

    return point_density


# 计算密度最大的点的聚类中心距离
def get_max_dis(distance_all, point_density):
    point_density = point_density.tolist()
    a = max(point_density)  # 寻找密度最大的点
    b = point_density.index(a)  # 寻找密度最大的索引
    c = max(distance_all[b])  # 密度最大点的最大中心距离

    return c


# 到点密度大于自身最近点的距离
def get_min_dis(point, i, distance_all, point_density):
    min_distance = []
    if point:
        for j in point:
            min_distance.append(distance_all[i][j])
        return min(min_distance)
    else:
        max_distance = get_max_dis(distance_all, point_density)
        return max_distance


# 计算各点到聚类中心的距离(dist)
def get_each_dis(distance_all, point_density):
    nn = []
    for i in range(len(point_density)):
        aa = []
        for j in range(len(point_density)):
            if point_density[i] < point_density[j]:
                aa.append(j)
        ll = get_min_dis(aa, i, distance_all, point_density)
        nn.append(ll)
    return nn


# 分簇，根据需要划分的特征结果数找到数据密度中最大的n个点
def clustering(point_density, delta):
    mul1 = (point_density * delta).tolist()
    mul2 = mul1.copy()
    mul2.sort(reverse=True)
    v = np.random.rand(feature_num)
    for i in range(feature_num):
        for j in range(len(mul1)):
            if mul1[j] == mul2[i]:
                v[i] = mul1.index(mul1[j])
    return v


# 归类，判断所有点与V中点的距离，距离最近的划为一类，输出每个数据的类别
def classify(distance_all, v):
    c = np.random.rand(distance_all.shape[0])
    for i in range(distance_all.shape[0]):
        dist = []
        for j in range(len(v)):
            dist.append(distance_all[i][int(v[j])])
        c[i] = dist.index(min(dist)) + 1
    return c


# 计算最大最小距离
def max_min_dis(distance_all):
    max_dis = 0
    min_dis = 999
    for i in range(distance_all.shape[0]):
        for j in range(distance_all.shape[0]):
            if distance_all[i][j] == 0:
                continue
            if distance_all[i][j] > max_dis:
                max_dis = distance_all[i][j]
            if distance_all[i][j] < min_dis:
                min_dis = distance_all[i][j]
    return max_dis, min_dis


def main():
    array = read_file('data.txt')
    array = norm(array)
    distance_all = cal_dis(array)
    max_dis, min_dis = max_min_dis(distance_all)
    point_density = cal_den(distance_all, 0.5 * (max_dis + min_dis))
    nn = get_each_dis(distance_all, point_density)
    v = clustering(point_density, nn)
    c = classify(distance_all, v)
    print(c)
    x = []
    y = []
    z = []
    for i in range(len(array)):
        x.append(array[i][0])
        y.append(array[i][1])
        z.append(array[i][2])
    fig = plt.figure(0)
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(x, y)
    plt.figure(1)
    plt.scatter(point_density.tolist(), nn)
    plt.figure(2)
    plt.scatter(range(0, 500), c.tolist())
    plt.show()


if __name__ == '__main__':
    main()

数据集实例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
-65 96 67 15 11 73 1863   1105
3  87 30 15 16 73 1830   1535
-64    49 8  15 19 73 1006   1992
-34    93 12 15 12 73 1477   1614
9  139    29 15 11 73 1954   1796
178    174    40 15 18 73 1622   1007
-48    160    65 15 18 73 1508   1198
169    124    5  15 12 73 1837   1005
-59    -159   69 15 18 73 1324   1702
74 -174   17 15 12 73 1518   1178
-176   177    30 15 15 73 1647   1909
152    -137   37 15 13 73 1531   1606
72 -67    74 15 10 73 1164   1780
-10    -51    17 15 16 73 1979   1814
-21    -80    69 15 14 73 1026   1371
111    -83    7  15 10 73 1048   1093
-167   57 12 15 16 73 1117   1547
45 -89    55 15 14 73 1981   1986
83 -107   32 15 16 73 1701   1204

我这里是自己创建的数据集，就是根据数据中的某几维的不同状态判断属于那种状态，就是给大家一个参考，可以试一下。

image-20210511222630043

image-20210511222639063

其他技术点

sklearn 中 make_blobs模块的使用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
import numpy as np

if __name__ == '__main__':

    # make_blobs函数是为了聚类产生数据集
   data, label = make_blobs(n_features=2, n_samples=100, centers=3, random_state=3, cluster_std=[0.8, 2, 5])
    # n_features : 设置每个样本有几个特征值，默认值是2
    # n_samples : 取多少个样本数，默认值100
    # centers : 样本中心点有几个，默认值是3
    # random_state : 设置随机数种子，防止每次生成的数据都修改。默认是np.random。
    # cluster_std : 每个类别的方差。默认值是1.0   # shuffle : 洗乱，默认值是True
    # center_box : 中心确认之后的数据边界，默认值（-10.0， 10.0）

    # 进行样本绘制    plt.scatter(data[:, 0], data[:, 1], c=label)    plt.show()

img

绘制三维图

见参考资料绘制三维图

参考资料

Reducing the Dimensionality of Data with Neural Networks – G. E. Hinton and R. R. Salakhutdinov

基于密度峰值的聚类(DPCA)

密度峰值聚类算法MATLAB程序

python绘制三维图

总结

我这里写的只是一个很简单的一个demo，希望大佬们看个热闹，应该存在着不少理解错误的地方，代码冗余度还很高。因为还没有看关于高维数据的密度峰值聚类，也没使用标记过的数据集，以上数据集都是我自己生成的，所以分类结果也不知道误差是多大。原理应该是没错。以后有机会，会继续完善这部分代码，使用上高维数据，也会用到更加真实的数据集进行测试。