DBSCAN & OPTICS

我如期来更新啦!!!聚类算法是很常用的一种算法,不过最常见的就是KMeans了,虽然很多人都会用,不过讲道理,如果是调用现成机器学习库里面的KMeans的话,我敢保证90%的人答不上来具体的是什么算法。相信我,库里的KMeans跟教科书讲的那个随机取初始点的KMeans不是一个算法哟~

因为KMeans依赖K,但是我怎么知道K要用多少呢?另外,KMeans受限于算法本身,对于球状的数据效果较好,但是不规则形状的就不行了。这种情况下,相对而言,基于密度的聚类算法就比较好用了。sklearn里面现在是放了一个DBSCAN,下一版会更新OPTICS。刚好最近都用了,这里把DBSCAN跟OPTICS算法复现一遍。

DBSCAN

DBSCAN算法的原始论文是96年的这篇《A Density-Based Algorithm for Discovering Clusters in Large Spatial Database with Noise》

DBSCAN是一种基于密度的聚类算法,也就是说,密度高的区域自动聚成一类。这样一来,我们就避免了人为去设定群组数量的问题,算法可以自动发现群组数量。另外用这种方法,如果一个sample不在高密度区域,就有可能被判定为异常值,那么也可以拿来作为异常值检验的方法。

DBSCAN的思路非常简单,有两个参数,一个是\(\varepsilon\),另一个是minimum points。这里首先定义DBSCAN的几个核心概念,一个是\(\varepsilon\)-neighborhood,另一个是core object,还有就是reachable。

首先是\(\varepsilon\)-neighborhood。DBSCAN开始的时候会随机选取一个初始点,然后按照\(\varepsilon\)的距离寻找临近点,这些点的集合叫做\(\varepsilon\)-neighborhood。参数里面的minimum points就是限定这个点的集合的。minimum points限定了这个集合最小需要包含多少样本点,\(\varepsilon\)则是限定了要用多大的范围去框定这些样本点。这里有个小细节要注意,那就是,算neighborhood的时候,中心点自己是算进来的。

那么这样圈完一波neighborhood后,我们会将符合样本数量大于等于minpts的中心点叫做core object。而core object跟这些neighbor就是reachable的。

然后为了让算法运作起来,只要neighborhood这个集合里面有点,我们就不断重复这样圈地的动作,然后把中心点从集合中拿掉,直到neighborhood为空。

那我们很自然就会想到,一定会圈到一些点,它们不是core object,但是也在集合里面。这些点我们叫做border,也就是说,这些点是这个类的边界了。那么我们就很自然会想到,也有一些点压根不在neighborhood里面,也不是core的点,这些点就是noise。那既然样本点多了两个状态,reachable的情况也就变得多了,如果是直接可以在neighborhood里面找到的,我们叫做directly-reachable;如果通过neighborhood一层层找,最后找到的,我们叫density-reachable。

可以看图说话:

这个图里面,A是core,B、C是border,N是noise。A跟B、C都是density-reachable。但是有没有发现,B是没法返回去找到A的。所以这种reachable是有方向的。

如果对SNA有点了解的朋友就知道,就是一度人脉和N度人脉,但是是一个有向图。

那么算法思路理清了,代码就好写了,这里我就用最常用的欧氏距离了。有人可能会想,如果还是欧氏距离,那跟KMeans还有什么分别,都是画圈圈嘛!请回想一个微积分,只要圈画的够小,就能做出各种形状来。下面是核心部分的代码,详细的可以去看我的notebook

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
def _euclidean_dist(p, q):
return np.sqrt(np.power(p - q, 2).sum())


def _eps_reachable(p, q, eps):
return _euclidean_dist(p, q) < eps


def _eps_neighborhood(data, point_id, eps):
neighborhood = []
point = data[point_id]
for q in range(len(data)):
if _eps_reachable(point, data[q], eps):
neighborhood.append(q)
return neighborhood


def _is_core_object(neighborhood, minpts):
if len(neighborhood) >= minpts:
return True
else:
return False


def dbscan(data, eps, minpts):
class_id = 0
class_label = np.full(len(data), -1)
for p_id in range(len(data)):
if class_label[p_id] == -1:
neighborhood = _eps_neighborhood(data, p_id, eps)
if _is_core_object(neighborhood, minpts):
class_label[neighborhood] = class_id
neighborhood.remove(p_id)
while len(neighborhood) > 0:
current_p = neighborhood[0]
current_neighborhood = _eps_neighborhood(data, current_p, eps)
for n in current_neighborhood:
if class_label[n] == -1:
class_label[n] = class_id
neighborhood.append(n)
neighborhood = neighborhood[1:]
class_id += 1
return class_label

OPTICS

那么DBSCAN本身是一个非常牛逼的算法,它解决了我们找K的问题,这样在海量群组的时候,我们不用像KMeans一样去到处尝试K的大小。但是DBSCAN有个问题,那就是这个算法只能检测一个密度。换句话说,如果现在存在一个数据集有两个类,一个类是方差小的,一个类是方差大的。且这两个群组离得不算太远。如果我们为了照顾方差大的群组将eps设得很大,minpts设得很小,那么可能把两个类聚在一起。反过来,我们就可能找不到方差大的类。

那么问题来了,有没有办法量化这个距离呢?三年后,同一组作者在DBSCAN的基础上进化出了OPTICS算法

既然是DBSCAN的进化版,所以很多概念上都是互通的,只是OPTICS算法多了几个概念,一个是core-distance,一个是reachability-distance。

我们知道,DBSCAN是不断跑马圈地的一个过程,但是我们很直观想就知道,有些密度大的地方,可能不需要\(\varepsilon\)那么大的范围就可以圈到minpts个样本,所以在OPTICS算法里面,我们将满足minpts这么多样本点的\(\varepsilon'\)叫做core-distance。而reachability-distance就是中心点与临近点的距离,但是,如果临近点落在\(\varepsilon'\)内,reachability-distance就用core-distance来替代。如下图:

那么用OPTICS的时候我们就需要定义两个列表,一个是seeds,一个是ordered result。seeds就是我们每一轮迭代时候的候选列表,而ordered result就是最终的结果。

具体的过程是这样的。我们先找到一个点,然后一样跑马圈地,接着计算reachability-distance,然后放到seeds里面从小到大排序。每次取第一个seed出来继续圈地,把被取出来的点以及这个点的reachability distance存在ordered result里面。接着就跟DBSCAN一样,不断重复,直到neighborhood为空。这样做的好处就是,我们可以量化评估每个群的密度大小。效果如下图:

那么我们又会想到,有些seeds里面的点可能随着核心点的移动,reachability distance会不断变小。因为A的core distance里可能是B,而C不在A的core distance里面,但是C在B的core distance里面。如果第一个处理的点是A,第二个处理的点是B,那C其实还是很核心的一个点。那这种时候我们就要跟一它开始的reachability distance做比较,如果新的reachability distance比原来的小,就把原来的值替换掉。

那么废话不多说,上代码,很多地方跟DBSCAN是可以复用的,我就放了一些OPTICS的核心部分,老规矩,详细见notebook

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
def _eps_neighborhood(data, point_id, eps):
neighborhood = []
rdist = []
point = data[point_id]
for q in range(len(data)):
dist = _euclidean_dist(point, data[q])
if dist < eps:
neighborhood.append((q, dist))
neighborhood = np.array(neighborhood, dtype=[('id', 'i'), ('dist', 'f')])
neighborhood = np.delete(neighborhood, np.where(neighborhood['id'] == point_id))
return neighborhood


def _update_order_seeds(neighborhood, minpts, reach_dists, processed, seeds):
if not _is_core_object(neighborhood, minpts):
return seeds
core_dist = np.sort(neighborhood['dist'])[minpts - 2]
for obj in neighborhood['id']:
if not processed[obj]:
obj_dist = neighborhood[neighborhood['id'] == obj]['dist'][0]
newRdist = max(core_dist, obj_dist)
if np.isnan(reach_dists[obj]):
reach_dists[obj] = newRdist
seeds[obj] = newRdist
elif newRdist < reach_dists[obj]:
reach_dists[obj] = newRdist
seeds[obj] = newRdist
return seeds


def optics(data, eps, minpts):
processed = np.array([False] * len(data))
core_dists = np.full(len(data), np.nan)
reach_dists = np.full(len(data), np.nan)
ordered_res = []
seeds = {}
for p_id in range(len(data)):
if not processed[p_id]:
ordered_res.append(p_id)
processed[p_id] = True
neighbors = _eps_neighborhood(data, p_id, eps)
if _is_core_object(neighbors, minpts):
core_dists[p_id] = np.sort(neighbors['dist'])[minpts - 2]
seeds = _update_order_seeds(neighbors, minpts, reach_dists, processed, seeds)
while len(seeds) > 0:
nextId = sorted(seeds.items(), key=operator.itemgetter(1))[0][0]
del seeds[nextId]
processed[nextId] = True
ordered_res.append(nextId)
newNeighbors = _eps_neighborhood(data, nextId, eps)
if _is_core_object(newNeighbors, minpts):
core_dists[nextId] = np.sort(newNeighbors['dist'])[minpts - 2]
seeds = _update_order_seeds(newNeighbors, minpts, reach_dists, processed, seeds)
# 这里只是我的一个操作,强迫症看不惯有个nan存在。
reach_dists[ordered_res[0]] = core_dists[ordered_res[0]]
return ordered_res, reach_dists, core_dists


def cluster_optics_dbscan(data, eps_reachable, eps, minpts):
orders, reach_dists, core_dists = optics(data, eps, minpts)
n_samples = len(data)
labels = np.zeros(n_samples, dtype=int)
far_reach = reach_dists > eps_reachable
near_core = core_dists <= eps_reachable
labels[orders] = np.cumsum(far_reach[orders] & near_core[orders])
labels[far_reach & ~near_core] = -1
return labels, orders, reach_dists, core_dists

OPTICS的优点就是,不管是什么形状的密度,基本上都可以把这个凹槽给跑出来,但是问题就是最后的这个抽取群组小算法。目前我还没找到一个比较好的方法来自动抽取,如果是按照论文里面的分层抽取,我试过会抽的太细,如果是按照论文里面的DBSCAN来抽,就是我实现的这个,不过是一刀切的方式,太复杂的样本效果就不好了。目前还在探索用其他平滑方法来替代,有突破再来更新。

试了两种平滑方式,最后的做法是,先用一维高斯平滑,然后用max filter抹掉因为高斯平滑最大值向左漂的问题。

HDBSCAN

然后又发现了一个新的聚类算法,类似OPTICS算法,叫HDBSCAN,顾名思义,就是H就是hierarchical。事实上就是OPTICS算法的一个改进版本。

与OPTICS类似,HDBSCAN一样要算core distance和reachability distance。这个算法跟OPTICS是一样的。但是reachability distance的公式是:\(d_{\text{reach}-k}(a, b) = \max\{\text{core}_k(a), \text{core}_k(b), d(a, b)\}\),这里的\(d(a, b)\)是两个点用距离公式算的距离。如下图所示:

蓝绿两点直接的reachability distance就是绿点的core distance,红绿两点的就是红绿两点的距离。

那么这样一来,就可以得到每个点之间的reachability distance。然后就可以用这个距离作为权重来绘制一个带权重的连接图。Beyond Hartigan Consistency这篇论文认为,无论是那种概率密度分布,这种相互距离可以更好的表示单链接聚类的层次结构。

然后可以用Prim算法构建最小生成树,然后转化为层次结构,就类似层次聚类了。这个可以参考一下基于最小生成树的层次聚类。然后通过设定最小簇大小,对这个树进行剪枝。也就从下面的第一个图变成第二个图:

之后就是抽取簇。上面剪枝后的树其实已经有了结果,但是我们希望抽的簇可以自动抽取,且是稳定的。因此我们需要一种度量方式来衡量簇的稳定性。我们定义一个值\(f(x) = \lambda = \frac{1}{\text{distance}_{\text{core}}}\)。这里的distance就是一个点。所以用\(\lambda_{\text{birth}}\)\(\lambda_{\text{death}}\)分别表示簇生成和簇分裂时候的\(\lambda\)大小,也就是\(\lambda_{\text{min}}C_i\)\(\lambda_{\text{max}}C_i\),那么簇的稳定性就是\(S = \sum_{p \in \text{cluster}} (\lambda_p - \lambda_{\text{birth}})\)\(\lambda_p = \lambda_{\text{max}}(x, C) = \min(f(x), \lambda_{\text{max}}C)\)每轮优化的方向就是让这个稳定性最大,同时满足最小簇大小。

那么HDBSCAN跟DBSCAN比的话,DBSCAN实际上就是上面剪枝前的hierarchical tree切一刀,而HDBSCAN会自适应地去寻找合适的划分。

乞讨码