代码之家 › 专栏 › 技术社区 › PDRX

Scipy的cut_tree()不返回请求的簇数,并且使用Scipy和fastcluster获得的链接矩阵不匹配

hierarchical-clustering data-analysis scipy debugging python

PDRX · 技术社区 · 7 年前

AHC )使用 fastcluster scipy.cluster.hierarchy 模块功能,in Python 3 ,我发现 cut_tree()

我毫无问题地对数据进行聚类,得到一个链接矩阵, Z linkage_vector() 具有 method=ward . 然后,我想切割dendogram树以获得固定数量的聚类(例如33),我使用 cut_tree(Z, n_clusters=33) . (请记住,AHC是一种确定性方法,生成一个连接所有数据点的二叉树,这些数据点位于树的叶子上;您可以在任何级别查看此树,以“查看”最终想要的簇数;cut\u tree()所做的一切就是返回一组从0到n\u clusters-1的“n\u cluster”整数标签,归属于数据集的每个点。)

我在其他实验中已经做过很多次了,我总是能得到我要求的集群数量。问题是,当我问这个数据集时 cut_tree() 对于33个集群,它只给了我32个。我不明白为什么会这样。可能是虫子吗?你知道有什么错误吗 ? 我试图调试这种行为,并使用scipy的 linkage() 作用将得到的连杆矩阵作为输入 cut_tree() 我没有得到意外数量的集群作为输出。我还验证了两种方法输出的链接矩阵不相等。

[dataset] 我使用的是10680个向量,每个向量有20个维度。检查以下实验:

import numpy as np
import fastcluster as fc
import scipy.cluster.hierarchy as hac
from scipy.spatial.distance import pdist

### *Load dataset (10680 vectors, each with 20 dimensions)*
X = np.load('dataset.npy')

### *Hierarchical clustering using traditional scipy method*
dists = pdist(X)
Z_1 = hac.linkage(dists, method='ward')

### *Hierarchical clustering using optimized fastcluster method*
Z_2 = fc.linkage_vector(X, method='ward')

### *Comparissons*

## Are the linkage matrices equal?
print("Z_1 == Z_2 ? ", np.allclose(Z_1, Z_2))

## Is scipy's cut_tree() returning the requested number of clusters when using Z_2?
print("Req.\tGot\tequal?")
for i in range(1,50):
    cut = hac.cut_tree(Z_2, i)
    uniq = len(np.unique(cut))
    print(i,"\t",uniq,"\t",i==uniq)

## The same as before, but in condensed form. When requesting cut_tree() for clusters
#  in the range [1,50] does it return wrong results at some point?
print("Any problem cutting Z_1 for n_clusters in [1,50]? ", not np.all([len(np.unique(
                                      hac.cut_tree(Z_1, i)))==i for i in range(1,50)]))
print("Any problem cutting Z_2 for n_clusters in [1,50]? ", not np.all([len(np.unique(
                                      hac.cut_tree(Z_2, i)))==i for i in range(1,50)]))

#Output:
#
#Z_1 == Z_2 ?  False
#
#Req.    Got     equal?
#1        1       True
#2        2       True
#3        3       True
#4        4       True
#5        5       True
#6        6       True
#7        7       True
#8        8       True
#9        9       True
#10       10      True
#11       11      True
#12       12      True
#13       13      True
#14       14      True
#15       15      True
#16       16      True
#17       17      True
#18       18      True
#19       19      True
#20       20      True
#21       21      True
#22       22      True
#23       23      True
#24       24      True
#25       25      True
#26       26      True
#27       27      True
#28       28      True
#29       29      True
#30       30      True
#31       31      True
#32       32      True
#33       32      False
#34       33      False
#35       34      False
#36       35      False
#37       36      False
#38       37      False
#39       38      False
#40       39      False
#41       40      False
#42       41      False
#43       42      False
#44       43      False
#45       44      False
#46       45      False
#47       46      False
#48       47      False
#49       48      False
#
#Any problem cutting Z_1 for n_clusters in [1,50]?  False
#Any problem cutting Z_2 for n_clusters in [1,50]?  True

您可能已经注意到,数据集包含37个向量,其中至少有一个精确副本,并且计算所有副本时,数据集中总共有55个向量,其中至少有一个副本。

Z_1 在顶部和 Z_2

因此,我再次运行了如前一段代码所示的聚类过程,但这次只使用了包含至少有一个副本的55个向量的数据子集。我获得 X_subset 使用:

uniqs, uniqs_indices, uniqs_count = np.unique(X, axis=0, return_index=True, return_counts=True)
duplicate_rows_indices = list( set(range(len(X))) - set(uniqs_indices) )
number_of_duplicate_rows = len(X)-len(uniqs) # 37

all_duplicate_rows = set()
for i in duplicate_rows_indices:
    _rows = set(np.where(X == X[i])[0])
    for j in _rows:
        all_duplicate_rows.add(j)

rows_with_at_least_a_copy = list(all_duplicate_rows)
number_of_rows_with_at_least_a_copy = len(rows_with_at_least_a_copy)  # 55

X_subset = X[rows_with_at_least_a_copy]

#Z_1 == Z_2 ?  False
#Req.    Got     equal?
#1        1       True
#2        2       True
#3        2       False
#4        3       False
#5        4       False
#6        5       False
#7        6       False
#8        7       False
#9        8       False
#10       9       False
#11       10      False
#12       11      False
#13       12      False
#14       13      False
#15       14      False
#16       15      False
#17       16      False
#18       17      False
#19       18      False
#20       20      True
#21       21      True
#22       22      True
#23       23      True
#24       24      True
#25       25      True
#26       26      True
#27       27      True
#28       28      True
#29       29      True
#30       30      True
#31       31      True
#32       32      True
#33       33      True
#34       34      True
#35       35      True
#36       36      True
#37       37      True
#38       38      True
#39       39      True
#40       40      True
#41       41      True
#42       42      True
#43       43      True
#44       44      True
#45       45      True
#46       46      True
#47       47      True
#48       48      True
#49       49      True
#Any problem cutting Z_1 for n_clusters in [1,50]?  False
#Any problem cutting Z_2 for n_clusters in [1,50]?  True

因此,fastcluster和scipy不会返回相同的结果,如果只是因为重叠点,那么由于聚类情况的模糊性,这是可以接受的。但问题是 cut_tree() 在这些情况下,当给定通过以下方式获得的链接矩阵时,它有时不会返回请求的簇数: linkage_vector()

它也张贴在这里: https://github.com/scipy/scipy/issues/7977 .

1 回复 | 直到 7 年前

Daniel 7 年前

首先,我不关心这三种方法输出的差异 scipy.cluster.hierarchy.linkage() fastcluster.linkage() 和 fastcluster.linkage_vector() . 产生差异的原因有两个:

由于浮点算法的限制,簇距离存在微小差异(代数等价公式产生不同的结果)。
除了算术上的差异外,这些算法还可以以不同的方式解决关系。这方面的一个联系是两对集群 (A、B) (C,D) 之间的距离完全相同和和和 D

第二个问题是为什么 scipy.cluster.hierarchy.cut_tree() 这里真正有趣的问题是不能产生所需数量的簇。代数上,Ward聚类方法不能在树状图中产生所谓的反演。(当一步中的聚类距离大于下一步中的聚类距离时,会发生反转。)也就是说,逐步树状图第三列中的距离序列(输出 linkage() )理想情况下,是Ward方法的单调递增序列。然而,由于浮点不准确,在OP提供的数据集中,步骤35中的聚类距离报告为2.2e16 fastcluster。linkage_vector()

cut_tree() 不幸的是,不能很好地处理这些反演。即使它们发生在树状图的深处(如本例中),也会在合并过程的整个其余部分混淆算法,并产生错误形成的簇的影响。(乍一看,我认为树状图不仅切割到了错误的层次,而且聚类实际上是错误的。不过,我还没有完全分析过这种情况。)

最后,一个不完美的解决方案:为了解决这个问题,在SciPy的循环中替换两条标记线 cut_tree()

for i, node in enumerate(nodes):
    idx = node.pre_order()
    this_group = last_group.copy()
    # ------------- Replace this:
    this_group[idx] = last_group[idx].min()
    this_group[this_group > last_group[idx].max()] -= 1
    # -------------
    if i + 1 in cols_idx:
        groups[np.where(i + 1 == cols_idx)[0]] = this_group
    last_group = this_group

通过以下行:

    unique_idx = np.unique(last_group[idx])
    this_group[idx] = unique_idx[0]
    for ui in unique_idx[:0:-1]:
        this_group[this_group > ui] -= 1

(请看 SciPy source code

然而,我宁愿建议重新实现 cut_tree() 方法从头开始,因为当前实现过于复杂,并且任务可以在运行时复杂性方面更高效地执行。