sklearn聚类之OPTICS算法

简介

OPTICS算法，全称是Ordering points to identify the clustering structure，是一种基于密度的聚类算法，是DBSCAN算法的一种改进。

众所周知，DBSCAN算法将数据点分为三类：

1. 核心点：若样本$x_i$的$\epsilon$邻域内至少包含了$M$个点，则为核心点
2. 边界点：若样本$x_i$的$\epsilon$邻域内包含的点数小于$M$，但在其他核心点的$\epsilon$邻域内，则为边界点
3. 噪声：既非核心点也非边界点则为噪声

这里面有两个关键参数，即$\epsilon$邻域内点的个数$M$，二者作为判定条件，用以区分核心点、边界点以及噪声。这两个参数至关重要，甚至可以直接影响聚类结果。换言之，DBSCAN受经验影响，不同的参数会导致不同的聚类结果。

OPTICS的改进思路是，降低算法对$\epsilon$值的依赖，从而引入了核心距离和可达距离的概念的概念，即对于样本$x$，如果给定M，则使得$x$成为核心点的最小邻域半径为$x$的核心距离。如果$x$是核心点，若距离$x$最近的核心点为$y$，则可达距离为$y$的核心距离与$x,y$实际距离的最大值。

OPTICS的优越之处在于，可以为每个聚类簇设置不同的核心距离与可达距离，更能在点密度的意义上，为样本提供更加个性化的聚类结果。

sklearn实现

若将OPTICS算法的聚类结果进行绘制，能更加直观地理解可达距离的作用。在sklearn中提供了OPTICS类，测试如下

from sklearn.cluster import OPTICS
import matplotlib.pyplot as plt
import numpy as np

np.random.seed(0)   # 设置随机数种子
cens = [[-5, -2], [4, -1], [1, -2], [-2, 3], [3, -2],[5, 6]]
scales =[0.8, 0.1, 0.2, 0.3, 1.6, 2]
X = np.vstack([c+s*np.random.randn(250,2) for c,s in zip(cens, scales)])

clust = OPTICS(min_samples=50, xi=0.05, min_cluster_size=0.05)

# 开始聚类
clust.fit(X)

space = np.arange(len(X))
reachability = clust.reachability_[clust.ordering_]
labels = clust.labels_[clust.ordering_]

下面对聚类结果进行可视化演示

colors = ["g.", "r.", "b.", "y.", "c."]
# 绘制可达距离
for ind, color in enumerate(colors):
    Xk = space[labels == ind]
    Rk = reachability[labels == ind]
    plt.plot(Xk, Rk, color, alpha=0.3)

# 此为噪声
plt.plot(space[labels == -1], reachability[labels == -1], "k.", alpha=0.3)
plt.plot(np.full_like(space, 2.0, dtype=float), "k-", alpha=0.5)
plt.plot(np.full_like(space, 0.5, dtype=float), "k-.", alpha=0.5)

plt.tight_layout()
plt.show()

在上图中，横坐标为点的序号，纵坐标为可达距离，不同颜色代表OPTICS聚类得到的不同的距离。可以看出不同颜色截止时对应的$y$值是不同的，说明在OPTICS聚类的过程中，对不同的聚类簇生成了不同的参数。

上图在0.5和2.0处画了两条线，如果以这两条线所在位置为$\epsilon$进行DBSCN聚类，则意味着产生不同的聚类结果。

cluster_optics_dbscan

sklearn中提供了cluster_optics_dbscan函数，可以指定统一的可达距离，并进行聚类，调用如下

from sklearn.cluster import cluster_optics_dbscan

# 可达距离为0.5或者2时的DBSCN聚类
labels = [cluster_optics_dbscan(
    reachability=clust.reachability_,
    core_distances=clust.core_distances_,
    ordering=clust.ordering_,
    eps=0.5,
) for eps in [0.5, 2]]

然后可以对比一下这三种不同聚类方案的结果

import matplotlib.gridspec as gridspec

fig  = plt.figure()
ax = fig.subplots(1, 3)

labels = [clust.labels_] + labels
infos = [
    "Automatic Clustering\nOPTICS",
    "Clustering at 0.5 epsilon cut\nDBSCAN",
    "Clustering at 2 epsilon cut\nDBSCAN"
]

for i in range(3):
    for k, c in enumerate(colors):
        Xk = X[labels[i] == k]
        ax[i].plot(Xk[:, 0], Xk[:, 1], c, alpha=0.3)
    noise = X[labels[i] == -1].T
    indNoise = labels[i] == -1
    ax[i].plot(noise[0], noise[1], "k+", alpha=0.1)
    ax[i].set_title(infos[i])

绘图结果如下