Scipy 的 cut_tree() 不返回请求的簇数，并且使用 scipy 和 fastcluster 获得的链接矩阵不匹配
PDR*_*DRX 3 python debugging hierarchical-clustering data-analysis scipy
我正在使用fastcluster包与scipy.cluster.hierarchy模块函数进行凝聚层次聚类 ( AHC ) 实验，在 中，我发现cut_tree()函数的令人费解的行为。Python 3
我毫无问题地对数据进行聚类，并Z使用linkage_vector()with获得链接矩阵method=ward。然后，我想切割树状图树以获得固定数量的簇（例如 33），并且我使用 正确地执行此操作cut_tree(Z, n_clusters=33)。（请记住，AHC 是一种确定性方法，生成连接所有数据点的二叉树，这些数据点位于树的叶子；您可以在任何级别查看这棵树，以“查看”您最终想要的集群数量；所有cut_tree() 的作用是返回一组从 0 到 n_clusters - 1 的“n_cluster”整数标签，归因于数据集的每个点。）
我在其他实验中已经做过很多次，并且总是得到我请求的集群数量。问题是，对于这个数据集，当我要求cut_tree()33 个簇时，它只给我 32 个。我不明白为什么会出现这种情况。这可能是一个错误吗？您知道 的任何错误吗cut_tree()？我尝试调试这种行为，并使用 scipy 的links()函数执行相同的聚类实验。将生成的链接矩阵作为输入，cut_tree()我没有得到意外数量的簇作为输出。我还验证了两种方法输出的链接矩阵不相等。
我使用的 [数据集]由 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底部）。括号内的数字表示该分支中包含的人口；不带括号的数字是树的叶子（数字是 X 矩阵中向量的索引）。人们可以看到唯一的区别（在绘制的水平上）是用红色方块标记的分支，它们在 0 距离处合并，因为它们包含重叠的向量。
因此，我再次运行了前面代码中所示的聚类过程，但这次仅包含包含至少有一个副本的 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()。如何解决这个问题？
使用的库版本：scipy'0.19.1'，numpy'1.13.3'，fastcluster'1.1.24'
编辑：它也发布在这里： https: //github.com/scipy/scipy/issues/7977。
scipy.cluster.hierarchy.linkage()首先，我不关心 、 、fastcluster.linkage()和三种方法的输出差异fastcluster.linkage_vector()。差异可能因以下两个原因而产生：
- \n
由于浮点运算的限制，簇距离存在微小差异（代数等效公式产生不同的结果）。
\n除了算术差异之外，算法还可以以不同的方式解决平局问题。在这种情况下， A \xe2\x80\x9ctie\xe2\x80\x9c 是两对簇(A,B)、(C,D) ， A和B以及C和D之间的距离完全相同。正如OP所写，例如在过程开始时有许多距离为0的点对，算法可以按任何顺序对它们进行聚类。
\n
其次，为什么scipy.cluster.hierarchy.cut_tree()没有产生所需数量的簇是这里真正有趣的问题。从代数上来说，\xe2\x80\x9cWard\xe2\x80\x9d 聚类方法无法在树状图中产生所谓的\xe2\x80\x9cinversions\xe2\x80\x9d。（当一步中的聚类距离大于下一步中的聚类距离时，就会发生反转。）也就是说，逐步树状图的第三列中的距离序列（ 的输出linkage()）理想地是单调递增序列\xe2\x80\x9cWard\xe2\x80\x9d 方法。然而，由于浮点不准确，在 OP 提供的数据集中，步骤 35 中的聚类距离报告为 2.2e\xe2\x88\x9216 fastcluster.linkage_vector()，但在下一步 36 中为 0.0。
不幸的是， SciPycut_tree()不能很好地处理这些反转。即使它们发生在树状图中的\xe2\x80\x9c深处\xe2\x80\x9d（如本例中），这也会使整个合并过程的其余部分的算法感到困惑，并产生错误形成的簇的效果。（乍一看，我认为树状图不仅 \xe2\x80\x9ccut\xe2\x80\x9d 处于错误的级别，而且簇实际上是错误的。不过，我还没有完全分析情况。）
这甚至更令人不安，因为反转不仅是由于当前情况下的数值不准确而发生的，而且 \xe2\x80\x9ccentroid\xe2\x80\x9d 和 \xe2\x80\x9cmedian\xe2\x80\x9d 聚类方法也会产生它们即使有完美的算术，这种情况也经常发生。
\n\n最后，一个不完美的解决方案：为了解决这个问题，替换SciPy中循环中的两行标记行cut_tree()：
for i, node in enumerate(nodes):\n idx = node.pre_order()\n this_group = last_group.copy()\n # ------------- Replace this:\n this_group[idx] = last_group[idx].min()\n this_group[this_group > last_group[idx].max()] -= 1\n # -------------\n if i + 1 in cols_idx:\n groups[np.where(i + 1 == cols_idx)[0]] = this_group\n last_group = this_group\n通过以下几行：
\n\n unique_idx = np.unique(last_group[idx])\n this_group[idx] = unique_idx[0]\n for ui in unique_idx[:0:-1]:\n this_group[this_group > ui] -= 1\n（查看SciPy源代码了解上下文。）
\n\n但是，我宁愿建议重新实施cut_tree()从头开始重新实现该方法，因为当前的实现过于复杂，并且在运行时复杂性方面可以更有效地执行任务。