机器学习聚类算法之HDBSCAN

先前的文章中介绍了基于密度的聚类方法DBSCAN，今天要学习的是 HDBSCAN。单从名字上看，两者必然存在一定的关系。我们先来看看官方的介绍：

HDBSCAN – Hierarchical Density-Based Spatial Clustering of Applications with Noise. Performs DBSCAN over varying epsilon values and integrates the result to find a clustering that gives the best stability over epsilon. This allows HDBSCAN to find clusters of varying densities (unlike DBSCAN), and be more robust to parameter selection.

从介绍中我们可以知道是 DBSCAN 算法与基于层次聚类算法结合而来的。DBSCAN 算法的原理是：对于聚类中的每个对象，在给定的半径邻域内的数据对象必须超过某个阀值。其算法简洁，对噪声点不敏感，而且可以发现任意形状的簇，但还是存在不足之处：

• 由于需要在整个数据空间构建树，算法需要很大的 IO 开销
• 算法输入参数没有一个很完美的科学标准来作为参考，这就使得人为干扰的因素变得很大，参数选取略有偏差对于聚类的效果有时会呈现出完全不同的效果

HDBSCAN 算法是对OPTICS 算法的一种改进，但并不是没有缺点。比如其对于边界点的处理方面效果却不是很理想。

HDBSCAN 算法原理

关于算法原理，这里直接引用官方文档的内容进行说明。首先创建一些测试点：

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import sklearn.datasets as data
%matplotlib inline
sns.set_context('poster')
sns.set_style('white')
sns.set_color_codes()
plot_kwds = {'alpha': 0.5, 's': 80, 'linewidths': 0}

moons, _ = data.make_moons(n_samples=50, noise=0.05)
blobs, _ = data.make_blobs(n_samples=50, centers=[(-0.75, 2.25), (1.0, 2.0)], cluster_std=0.25)
test_data = np.vstack([moons, blobs])
plt.scatter(test_data.T[0], test_data.T[1], color='b', **plot_kwds)

要清楚的解释 HDBSCAN 的原理，最好的办法是在实际使用过程中去了解中间发生了什么。

HDBSCAN 的使用方式

import hdbscan

clusterer = hdbscan.HDBSCAN(min_cluster_size=5, gen_min_span_tree=True)
clusterer.fit(test_data)

上述代码非常的简单，但中间可以把它拆成如下几个步骤：

• 根据密度/稀疏度对空间进行变换
• 建立距离加权图的最小生成树
• 构造连接组件的簇层次结构
• 根据最小的簇大小压缩簇层次结构
• 从压缩树中提取稳定的簇

变换空间

为了找到簇，我们希望在一片稀疏的噪音海洋中找到密度更高的孤岛。聚类算法的核心是单链接聚类，它对噪声非常敏感：一个位于错误位置的单个噪声数据点可以充当岛屿之间的桥梁，将它们粘合在一起。显然，我们希望我们的算法对噪声是鲁棒的，所以我们需要找到一种方法，以帮助”降低海平面”之前运行一个单一的连接算法。

我们如何在不进行聚类的情况下描述“海洋”和“陆地”?我们只要能够得到一个密度的估计，我们就可以把密度较低的点看作是“海洋”。这里的目标不是完全区分”海洋”和”陆地”，只是为了使我们的簇核心对噪音更加健壮。因此，鉴于”海洋”的定义，我们希望降低海平面。就实际目的而言，这意味着使”海洋”中的点彼此之间和”陆地”之间的距离更远。

然而，这只是设想。它在实践中是如何工作的？我们需要一个非常低成本的密度估计，最简单的是到 kth 最近邻距离。将其称为为针对点 x 的参数 k 定义的核心距离（定义为当前点到其第 k 近的点的距离），并表示为 $core_k(x)$：

$$core_k(x)=d(x,N^{\mathcal{k}}(x))$$

现在我们需要一种方法，以低密度(相应的高核心距离)分散点。要做到这一点，简单的方法是定义一个新的点之间的距离度量，我们将调用相互可达距离。我们将相互可达距离定义如下：

$$d_{\mathrm{mreach-}k}(a,b)=\max\{\mathrm{core}_k(a),\mathrm{core}_k(b),d(a,b)\}$$

式中，$d(a,b)$ 是 a 与 b 的原始距离。在该式中密集点（核心距离较低）彼此保持相同的距离，但较稀疏的点被推开，以使其核心距离至少远离任何其他点。这实际上”降低了海平面”，稀疏的”海洋”指向外界，而”陆地”则没有受到影响。这里需要注意的是，这显然取决于 k 的选择，较大的 k 值将更多的点解释为处于“海洋”中。所有这些用一张图片来说都比较容易理解，我们使用 k 值为 5，然后对于给定的一个点，我们可以画一个核心距离的圆，作为与第六个最近邻接触的圆(包括点本身)，如下所示:

再选择另外一个点，我们可以做同样的事情，这一次用一组不同的邻居（其中一个甚至包含我们选择的第一个点）：

我们可以再用另一组六个最近邻，和另一个半径略有不同的圆：

现在，如果我们想知道蓝点和绿点之间的相互可达距离，我们可以先画一个箭头，给出绿点和蓝点之间的距离：

它穿过蓝色的圆圈，但不是绿色的圆圈——绿色的核心距离大于蓝色和绿色之间的距离。因此，我们需要将蓝色和绿色之间的相互可达距离标记为大于等于绿色圆的半径。另外，从红色到绿色的相互反应距离就是从红色到绿色的距离，因为这个距离大于两个核心距离：

一般来说，有潜在的理论来证明，相互可达距离作为一种变换，可以很好地允许单链接聚类更接近水平集的层次结构，无论我们采样的点的实际密度分布是什么。

建立最小生成树

现在我们在数据上有了一个新的相互可达性度量，我们希望开始在稠密数据上寻找孤岛。当然，密集区域是相对的，不同的岛屿可能有不同的密度。从概念上讲，我们将要做的是: 将数据看作一个加权图，其中数据点为顶点，任意两点之间的边的权重等于这些点之间的相互可达距离。

现在考虑一个阈值，从高开始，逐步降低。删除任何重量超过该阈值的边。当我们删除边时，我们将开始断开图形的连接组件。最终，我们将在不同的阈值水平上得到一个连接组件的层次结构(从完全连接到完全不连接)。在实践中，这是非常低效的：我们有$n^2$个边，并且不期望连接的组件算法运算那么多次。正确的做法是找到一个最小的边集合，这样从集合中删除任何边都会导致组件断开。幸运的是，图论为我们提供了这样一个东西: 图的最小生成树。

我们可以通过Prim算法 非常有效地构建最小生成树树 – 我们一次构建一条边，总是添加最小的权重边，将当前的树连接到树中还没有的顶点。您可以看到下面构造的HDBSCAN树。注意这是相互可达距离的最小生成树，它不同于图中的纯距离。在这个例子中，k值为5。

clusterer.minimum_spanning_tree_.plot(edge_cmap='viridis', edge_alpha=0.6, node_size=80, edge_linewidth=2)

构建簇层次结构

给定最小生成树，下一步是将其转换为连接组件的层次结构。这很容易以相反的顺序完成：根据距离对树的边进行排序(按增加的顺序)，然后遍历，为每条边创建一个新的合并的簇。这里唯一困难的部分是确定每个将2个簇接在一起的边，但可以通过联合查找 数据结构很容易实现。我们可以将结果视为树状图，如下图所示：

clusterer.single_linkage_tree_.plot(cmap='viridis', colorbar=True)

这就把我们带到了单链接停止的点。但是我们需要更多，簇层次结构虽然是好的，但是我们真正需要的是一组扁平的簇。我们可以通过在上面的图中画一条水平线，然后选择它穿过的簇来做到这一点。这就是DBSCAN 的实现逻辑（将水平剪枝上的任何单例簇声明为噪声）。问题是，我们如何知道界限在哪里？DBCSCAN只是将其作为一个(非常不直观的)参数。另外，我们真正想要处理的是可变密度的簇，任何剪枝的选择都是相互可达距离的选择，因此是单一的固定密度水平。理想情况下，我们希望能够在不同的地方剪枝来选择我们的簇。这就是HDBSCAN下一步要做的地方，它将创建与健壮的单链接不同的内容。

压缩簇层次结构

簇抽取的第一步是将庞大而复杂的簇层次结构压缩到一个更小的树中。正如上面的层次结构中看到的，通常情况下簇拆分是从一个簇中分离出一个或两个点，而不是将其视为一个簇拆分为两个新的簇。为了使这个具体化，我们需要一个最小簇大小的概念，我们将它作为HDBSCAN的一个参数。一旦我们有了最小簇大小的值，我们现在就可以遍历层次结构，并在每次分割时询问是否有一个由分割创建的新簇的点数少于最小簇大小。如果我们有少于最小的簇大小的点，我们声明它是’从簇中剔除的点’，并有较大的簇保留父簇的身份。另一方面，如果拆分为两个簇，每个簇至少与最小簇大小一样大，那么我们认为簇拆分就是让这个拆分保留在树中。在遍历了整个层次结构之后，我们最终得到了一个拥有少量节点的小得多的树，每个节点都有关于该节点的簇大小如何随着不同距离减小的数据.我们可以将其可视化为一个树状图，类似于上面的树状图，用线的宽度来表示簇中的点数。但是，当点被剔除时，该宽度随线的长度而变化。我们使用最小簇大小为5的数据，结果如下：

clusterer.condensed_tree_.plot()

这里更容易查看和处理，但是，我们仍然需要挑选出簇来作为扁平簇使用。看看上面的图，你就会知道如何去做这件事。

提取簇

直观地说，我们希望选择的簇能够持续存在并且有更长的生命周期; 短命的簇可能仅仅是单链接方法的产物。在前面的图中，我们可以说，我们要选择那些簇有最大面积的情节油墨。为了创建一个平面集群，我们需要添加一个进一步的要求，如果您选择了一个簇，那么您就不能选择它的后代的任何簇。事实上，关于应该做什么的直观概念正是HDBSCAN所做的。当然，我们需要把事情形式化，使之成为一个具体的算法。

首先，我们需要一种不同于距离的度量方法来考虑簇的持久性; 作为替换我们使用$\lambda=\frac{1}{distance}$。针对给定的簇，我们可以定义值$\lambda_{birth}$和$\lambda_{death}$为当簇分离并成为它自己的簇时的lambda值，以及当簇分别拆分为较小的簇时的lambda值（如果有）。反过来，对于给定的集群，对于集群中的每个点P，我们可以将值$\lambda_p$定义为“从簇中掉出来”的lambda值，该值介于$\lambda_{birth}$和$\lambda_{death}$之间。现在，对于每个簇，稳定性计算如下：

$$s_{\mathrm{cluster}} =\sum_{p\in\mathrm{cluster}}(\lambda_p-\lambda_{\mathrm{birth}})$$

其中：

• $\lambda_{\mathrm{birth}}$：团簇形成时的$\lambda$值
• $\lambda_{\mathrm{death}}$：团簇分裂为两个子团簇时的$\lambda$值
• $\lambda_p$：从团簇中分离出去时的$\lambda$值

将所有叶节点声明为选定的簇。现在通过遍历树（反向拓扑排序顺序）。如果子簇的稳定性之和大于簇的稳定性，那么我们将簇稳定性设置为子簇稳定性之和。另一方面，如果簇的稳定性大于其子簇的总和，那么我们将簇声明为选定簇，并取消选择其所有子簇。当我们到达根节点时，我们将当前选定的簇集合称为平面簇并返回它。

好的，这很复杂，但它实际上只是在执行我们的“选择图片中最大面积的簇”，这要受我们前面解释的后代约束。我们前面解释过。我们可以通过这个算法选择压缩树图中的聚类，你会得到你想要的结果：

clusterer.condensed_tree_.plot(select_clusters=True, selection_palette=sns.color_palette())

既然我们有了簇，那么根据 sklearn api 将其转化为簇标签就足够简单了。任何不在所选簇中的点只是一个噪声点(并被分配为 -1)。不过，我们可以做得更多：对于每个簇，我们都有该簇中每个点 P 的 $\lambda_p$；如果我们简单地规范化这些值（因此它们的范围从 0 到 1），那么我们就可以度量集群中每个点的簇成员资格的强度。HDBSCAN 库将此作为 cluster 对象的概率属性返回。因此，有了标签和成员强度，我们就可以制作标准图，基于簇标签为点选择颜色，并根据成员强度对颜色进行去饱和（并使未聚集的点纯灰色）。

palette = sns.color_palette()
cluster_colors = [sns.desaturate(palette[col], sat)
if col >= 0 else (0.5, 0.5, 0.5) for col, sat in
zip(clusterer.labels_, clusterer.probabilities_)]
plt.scatter(test_data.T[0], test_data.T[1], c=cluster_colors, **plot_kwds)

这就是 HDBSCAN 的工作原理。这可能看起来有点复杂，但最终每个部分实际上是非常简单的，可以很好地优化。希望通过更好地理解 HDBSCAN 的逻辑和一些实现细节，您会感到有动力去尝试它。该库将继续开发，并将为新的思路提供基础，包括一种近乎无参数的持久密度聚类算法和一种新的半监督聚类算法。

参考链接：https://hdbscan.readthedocs.io/en/latest/how_hdbscan_works.html

HDBSCAN 使用实例

import numpy as np
import pandas as pd
import hdbscan
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from math import pi, cos, sin, atan2, sqrt


def get_centroid(cluster):
    x = y = z = 0
    coord_num = len(cluster)
    for coord in cluster:
        lat = coord[0] * pi / 180
        lon = coord[1] * pi / 180

        a = cos(lat) * cos(lon)
        b = cos(lat) * sin(lon)
        c = sin(lat)

        x += a
        y += b
        z += c
    x /= coord_num
    y /= coord_num
    z /= coord_num
    lon = atan2(y, x)
    hyp = sqrt(x * x + y * y)
    lat = atan2(z, hyp)
    return [lat * 180 / pi, lon * 180 / pi]


df = pd.read_excel("test.xlsx")

hotel_df = df[['latitude', 'longitude']]
hotel_df = hotel_df.dropna(axis=0, how='any')
hotel_coord = hotel_df.values

hotel_dbsc = hdbscan.HDBSCAN(metric="haversine", min_cluster_size=int(len(hotel_df)/50)).fit(np.radians(hotel_coord))
hotel_df['labels'] = hotel_dbsc.labels_
hotel_df['probab'] = hotel_dbsc.probabilities_
hotel_df.loc[hotel_df['probab'] < 0.5, 'labels'] = -1 #HDBSCAN 边界可能存在问题，将置信度< 0.5 的设为为噪音点

cluster_list = hotel_df['labels'].value_counts(dropna=False)
center_coords = []
for index, item_count in cluster_list.iteritems():
    if index != -1:
        df_cluster = hotel_df[hotel_df['labels'] == index]
        center_coord = get_centroid(df_cluster[["latitude", "longitude"]].values)
        center_lat = center_coord[0]
        center_lon = center_coord[1]
        center_coords.append(center_coord)
center_coords = pd.DataFrame(center_coords, columns=['latitude', 'longitude'])
print(center_coords)

#可视化
fig, ax = plt.subplots(figsize=[20, 12])
facility_scatter = ax.scatter(hotel_df['longitude'], hotel_df['latitude'], c=hotel_df['labels'], cmap=cm.Dark2,
edgecolor='None',
alpha=0.7, s=120)
centroid_scatter = ax.scatter(center_coords['longitude'], center_coords['latitude'], marker='x', linewidths=2,
c='k', s=50)
ax.set_title('Facility Clusters & Facility Centroid', fontsize=30)
ax.set_xlabel('Longitude', fontsize=24)
ax.set_ylabel('Latitude', fontsize=24)
ax.set_xlim(120, 122)
ax.set_ylim(30, 33)
ax.legend([facility_scatter, centroid_scatter], ['Facilities', 'Facility Cluster Centroid'], loc='upper right',
fontsize=20)
plt.show()

项目地址：

• https://github.com/lmcinnes/hdbscan
• https://github.com/scikit-learn-contrib/hdbscan

官方文档：

• https://hdbscan.readthedocs.io/en/latest/index.html