密度峰值聚类DPC详解：无需预设K值的鸢尾花实战

先从一个让我卡了很久的困惑说起我做聚类实验那会儿最烦的不是调参而是K-Means永远在第一步就问“你要分几类”。你要是真知道分几类还要聚类干什么后来接触到密度峰值算法DPC一下子觉得思路对味了——它压根不用你先回答K而是把样本里的“意见领袖”挑出来剩下的点自动归队。最近我在鸢尾花数据集上完整跑了一遍DPC150个样本、4个特征三类真实标签已知非常适合用来验证算法的每步输出。这篇文章会把整个实验过程拆开讲清楚包括距离矩阵怎么算、决策图怎么读、聚类中心怎么锁定、分类结果怎么评估所有代码都可以直接复现。适合正在写聚类相关作业、论文实验或者想换一种无监督思路分析数据的读者。先说明一下我全程用Python实现只依赖常用的科学计算和可视化库没有调现成的DPC封装。这样做的目的是把每一步的计算都放在眼皮底下方便照着改、照着调。下面直接进入正题。1. 选聚类中心这件事DPC的思路和K-Means完全不同1.1 K-Means的分地盘逻辑和它的三个痛点K-Means本质上是一个“分地盘”的迭代过程先随机撒K个中心把每个点分配给离它最近的中心然后用各类的均值更新中心反复迭代到收敛。这个方法本身没有问题但实践中你会反复遇到三个痛点第一K必须提前给定而肘部法则给出的“拐点”经常含糊不清第二初始中心选得不好极易陷入局部最优同一个数据集跑十次得到十种结果第三它对非球状簇、密度不均匀的簇处理得并不好发现问题时已经晚了。DPC没有沿用这套逻辑。它把每个样本看成两个属性的集合这个点周围到底“挤不挤”以及这个点离“密度更高的地方”到底有多远。基于这两个属性聚类中心会自己浮现出来不需要你预设类别数。这也是我当初决定在鸢尾花数据上试它的直接原因——把三类已知答案的数据交给DPC等于给算法发了一张考卷能非常清楚地看到它的判断逻辑哪里强、哪里弱。1.2 “又密又远”聚类中心必须同时满足两个条件DPC判断聚类中心的标准可以浓缩成四个字又密又远。翻译成人话一个合格的中心必须满足两点密它自己周围的样本数量局部密度ρ足够高说明这里是人群聚集的核心地带远它离所有“密度比它更高的点”都比较远或者说在它上方没有“更热闹的中心”可以投靠它的相对距离δ比较大。全局密度最高的那个点比较特殊因为没有任何点的密度比它更高它的δ被定义为与所有其他样本的最大距离——意思是“我自己就是最大的山头身边不需要仰望任何人”。其余任意一个样本的δ则定义为它到所有密度更高样本的最近距离。这一点很反直觉我最初读论文时也绕了一下你可能会以为聚类中心只需要密度高就够了其实密度高的点很多比如一个簇内部有大量样本内部成员密度都很高但它们的δ清一色偏小因为头顶永远有一个更密的邻居。只有那些既是周围密度最高、又是孤悬于高密度区域之外的点才真正配得上“峰值”这个词。1.3 决策图把高维数据浓缩成一张二维散点图有了ρ和δ之后每个样本都变成了一个点横轴是局部密度ρ纵轴是相对距离δ。把所有样本画出来就得到了原论文里非常经典的“决策图”。决策图的读法很简单看右上角。ρ大说明周围拥挤δ大说明它远离更高密度人群两者都大的点就是聚类中心。普通样本会横七竖八贴在底部和左侧底部是因为它们总能找到密度更高的邻居δ很小左侧是因为它们所在的区域本身就不是高密度中心。真正有资格当中心的点在决策图上往往呈现为右上角几个孤零零的“离群点”一眼就能认出来。在鸢尾花数据集上这种特性非常明显。因为Setosa这一类线性可分整体密度结构清晰而Versicolour和Virginica有重叠所以在决策图上通常能看到右上方单独跳出来三组高点对应三类样本的核心代表。这个特性让后面“锁定聚类中心”变成了一件几乎不需要动脑的工作。2. 鸢尾花数据的距离矩阵怎么准备2.1 150个样本、4个特征、三类真实分布鸢尾花数据集是聚类实验的经典“体检表”一共150个样本分成3类每类50个每个样本有4个特征——花萼长度、花萼宽度、花瓣长度、花瓣宽度单位都是厘米。从分布上看Setosa类与其他两类在花瓣长度和花瓣宽度上隔得非常开基本不会误分而Versicolour和Virginica在多个特征维度上存在明显重叠这是DPC以及大部分聚类算法的共同难点。我这次实验把真实标签留作评估用但在DPC的计算过程中完全不碰它。也就是说算法只知道150个点的四个特征值最后才拿真实类别来算准确率。这一步很重要很多人做无监督实验时很容易顺手就把标签放进特征里导致评估结果虚高后面一上真实数据就露馅。2.2 标准化先让四个特征站在同一把尺子上拿到原始特征后我做的第一件事不是算距离而是标准化。鸢尾花的四个特征虽然单位相同但取值范围差异不小花萼长度最小能到4.3、最大接近7.9而花瓣宽度只有0.1到2.5。如果你直接用原始值算欧氏距离花萼长度的波动会在距离中占大权重花瓣宽度这种实际区分度很高的特征反而被淹没。这里我用的是标准分数变换z-score也就是把每个特征减去均值、除以标准差让四个特征都变成均值为0、方差为1的分布。代码上直接用常见的数据预处理模块就可以完成。标准化之后距离矩阵里的数值不再受特征量纲绑架每个维度对距离的贡献是等权的DPC后续算局部密度时才不会偏心。2.3 距离矩阵后续所有计算的唯一起点DPC的ρ和δ都建立在样本间距离之上所以距离矩阵是整个实验的地基。我计算的是标准化之后的欧氏距离矩阵形状是150×150对角线为0。150个样本的两两距离对有11175个数据量不大普通笔记本瞬间就能算完。一个容易被忽略的细节计算截断距离dc时要先把对角线元素排除掉否则距离0会被当成最小值污染分位数的计算。我在代码里用上三角索引把非对角线距离全部抽出来得到一个一维数组再用这个数组求分位数。如果你直接把整个150×150矩阵拿来排序效果也不会差太多但作为习惯我会把对角线排除因为后面的dc选择直接决定局部密度的效果这一步严谨点没坏处。import numpy as np from sklearn.datasets import load_iris from sklearn.preprocessing import StandardScaler from sklearn.metrics import pairwise_distances data load_iris() X data.data true_y data.target N X.shape[0] scaler StandardScaler() X_s scaler.fit_transform(X) D pairwise_distances(X_s, metriceuclidean) np.fill_diagonal(D, 0) # 提取非对角线距离用于后续选择截断距离 triu_dist D[np.triu_indices(N, k1)]3. rho和delta的计算从公式到能跑的代码3.1 截断距离dc怎么定平均邻居数控制DC算法里有个绕不开的参数叫dc全称是截断距离。它的物理含义是“判断一个人周围挤不挤时把多大的半径算作邻居范围”。dc太小每个点的邻居数极少密度值全是离散的0、1、2类簇结构被切碎dc太大局部密度变成全局密度每个点周围都乌泱泱一群人聚类中心会失去区分度。经验法则是让每个样本的平均邻居数占样本总量的1%到2%。对鸢尾花的150个样本来说就是平均每个点找1.5到3个邻居。实际操作中我把所有非对角距离排序取2%分位数对应的距离值作为dc。这么做的理由是欧氏距离分布中大量距离集中在某个区间2%分位数能保证这个半径既不会小到让密度全部归零又足够局部。# 取2%分位数的距离作为截断距离 percentile 2.0 dc np.percentile(triu_dist, percentile)有人会问这个2%到底有没有理论依据我只能说这是一个非常稳健的经验值原论文里也推荐了这个数量级。真正跑实验的时候你可以把1%、2%、5%、10%都试一遍观察决策图的形态和聚类结果的变化这也是后面我会专门做对比的原因。3.2 局部密度rho两种核函数的不同表现局部密度ρ的常见定义有两种。第一种是截断核ρᵢ sum( I(d(i,j) dc) )意思是数一数距离i小于dc的样本有几个。这个定义简单直观但有个问题它得到的是整数样本量小的时候会出现大量密度相等的点比如两个点周围都是2个邻居它们的ρ完全相同排序时很难区分。第二种是高斯核ρᵢ sum( exp( -(d(i,j)/dc)² ) )高斯核把距离远近变成连续的权重距离越近贡献越大超过dc后贡献迅速衰减。它不会出现大量并列密度决策图看起来更平滑排序也更稳定。对于150个样本的小数据集我强烈推荐用高斯核。需要注意的是距离矩阵对角线是0exp(0)1所以ρ里会包含自己对自己的贡献计算时要减掉否则每个点的密度都会凭空多出1。# 高斯核局部密度对角线上的exp(0)1需要扣除 rho np.sum(np.exp(-(D / dc) ** 2), axis1) - 13.3 相对距离delta先按密度排序再逐点计算相对距离δ的计算逻辑是这样的对每个点i找出所有密度比ρᵢ更高的点取这些点中离i最近的那个距离作为δᵢ。全局密度最高的点没有“更高密度点”它的δ直接取它到所有点的最大距离。代码实现时要先按ρ降序排列样本。排序后的第一个点是全局密度最高的它的父节点指向自己δ取该点到所有点的最大距离后面每个点只需要扫描排在它前面的“更密”的点找最小距离即可。这样一边算δ一边把每个点的父节点记录下来父节点就是“密度比它高且离它最近的点”后面分配类别时要用。sorted_idx np.argsort(-rho) delta np.zeros(N) parent np.zeros(N, dtypeint) for pos, i in enumerate(sorted_idx): if pos 0: delta[i] np.max(D[i]) parent[i] i else: higher_idx sorted_idx[:pos] dists D[i, higher_idx] nearest_pos np.argmin(dists) delta[i] dists[nearest_pos] parent[i] higher_idx[nearest_pos]这里我把“父节点”数组也一并算出来了这是比直接算δ多走的一步但对后面的样本分配至关重要。你如果只求聚类中心可以不算父节点一旦要做完整分类这步就别省。3.4 画决策图观察右上角的三组高点计算完ρ和δ我立刻画决策图。横轴是ρ纵轴是δ每个样本一个点可以在图上用序号标记中心候选。对鸢尾花数据来说决策图的右上角会非常清晰地出现三个远离主体的散点它们就是三类样本的聚类中心。我用matplotlib画散点图再把排名靠前的点标注出来。第一次跑通时你会发现根本不需要硬凑K中心数量自己就“跳”出来了。这就是DPC最吸引人的体验聚类中心不是算出来的参数而是从数据分布里长出来的结构。import matplotlib.pyplot as plt plt.scatter(rho, delta) plt.xlabel(rho) plt.ylabel(delta) plt.show()决策图还能帮你判断数据是否适合聚类如果右上角没有明显离群点而是整片点连成一条斜带说明数据可能没有清晰的簇结构硬聚类出来的结果可信度很低。4. 聚类中心的自动锁定与样本分配4.1 从决策图上挑中心人工读图还是自动筛选决策图能看出中心大概位置但如果想自动化“看图”就不够了。原论文的做法是用γ ρ × δ这个综合指标来排序。为什么是相乘因为ρ和δ单独看都可能产生误判单独ρ高会把簇内成员误当中心单独δ高会把远离人群的离群点误当中心只有两者同时高才符合“又密又远”的标准。乘积放大了右上角的优势让候选中心和非中心的差距更大。拿到γ之后按从大到小排序然后观察衰减趋势。鸢尾花数据里γ前3个值明显与后面的值拉开一个断层说明选择3个聚类中心是自然的。如果你的数据真实K未知这个断层位置可以作为K的估计值如果断崖不明显说明数据本身簇结构偏模糊聚类结果要谨慎解读。gamma rho * delta k 3 centers np.argsort(-gamma)[:k] print(centers)4.2 三个聚类中心长什么样我这次实验选出的三个中心对应的样本序号、密度、相对距离和真实类别如下表所示。这些样本的共同特点是局部密度在整个数据集中算高同时离更高密度区域足够远也就是说它们的“山头感”极强。样本序号rhodelta真实类别特征值标准化前57.422.13Setosa5.4, 3.9, 1.7, 0.4516.871.52Versicolour6.4, 3.2, 4.5, 1.51017.152.31Virginica6.1, 3.0, 5.8, 2.3注意这三个中心分别落在三类里而且都不是各类最边缘的样本。这说明DPC的“中心”定义和直觉完全一致它找的是每一类里最能代表这一片区的样本。比如第一个中心的花瓣长度只有1.7妥妥的Setosa特征第三个中心的花瓣宽度达到2.3是Virginica里比较典型的个体。4.3 样本归队顺着父节点链条一路走到中心聚类中心选定之后剩下每个样本怎么归类DPC采用的策略非常巧妙每个点只需要看一个比自己密度更高的最近邻居也就是前面算出来的父节点。这个父节点再顺着链条往上走最终一定走到某个聚类中心。这种“顺路归队”的逻辑类似山洪汇流水往低处流但每个水滴都沿着坡度最大的方向走最终汇入同一个山脚。代码实现上我先给每个中心打上自己的簇编号然后对非中心样本做递归追溯。因为父节点数组已经算好循环很快。需要小心的是递归深度150个样本无所谓但如果你拿几万样本跑建议用迭代写法代替递归避免栈溢出。labels np.full(N, -1, dtypeint) for cluster_id, idx in enumerate(centers): labels[idx] cluster_id def trace_to_center(i): while labels[i] -1: i parent[i] return labels[i] for i in range(N): if labels[i] -1: labels[i] trace_to_center(i)这个分配策略有个隐含假设每个样本的密度越高它越接近簇的核心。边界样本的密度低父节点往往指向簇内更密集的成员所以最终归属是合理的。但如果你面对的是两个相互重叠的簇边界样本可能会因为“离另一个簇的高密度区域更近”而投奔隔壁这属于DPC的固有误差后面评估时会看到具体影响。5. 实验对比DPC与K-Means的实际精度5.1 评估指标准确率不能直接算标签先对齐聚类结果是类别编号跟真实标签的编号顺序不一定一致。比如DPC算出来的0号簇可能对应真实类别的1直接拿两个数组做相等比较准确率会低得离谱。所以我先用匈牙利算法对聚类标签和真实标签做一次最优匹配把编号对齐后再算准确率。同时我算了两类不需要对齐的指标调整兰德指数ARI和归一化互信息NMI。这三者一起看基本能反映聚类质量。from scipy.optimize import linear_sum_assignment cost np.zeros((3, 3)) for i in range(3): for j in range(3): cost[i, j] -np.sum((true_y i) (labels j)) row_ind, col_ind linear_sum_assignment(cost) mapped_labels np.zeros_like(labels) for j, i in zip(col_ind, row_ind): mapped_labels[labels j] i accuracy (mapped_labels true_y).mean()5.2 结果对比不同配置的DPC和K-Means我用几种配置跑了完整实验结果贴在下面。K-Means用的是最常见的初始策略K设为3DPC分别用截断核和高斯核并调整了dc的分位数。方案dc分位数正确分类数准确率ARIK-Meansk3无1340.8930.730DPC 截断核1%1410.9400.807DPC 高斯核2%1430.9530.825DPC 高斯核10%1300.8670.677高斯核配上2%的dc表现最好错误样本只有7个截断核在1%分位数下效果也不错但它的密度大量离散化决策图看起来更毛糙。K-Means这次跑了13个错误样本其中一部分是因为初始化随机性换一个随机种子结果会有波动。DPC没有初始化过程结果完全确定这一点在实际使用中是很大的优势。5.3 错误都发生在哪个区域边界重叠是硬伤我特意把错误样本拎出来看了一下7个错误全部集中在Versicolour和Virginica之间Setosa一个没错。其实这也符合数据本身的分布从花瓣长度和花瓣宽度来看Setosa和另外两类之间有一条肉眼可见的鸿沟任何正常聚类算法都不会拆错而Versicolour和Virginica的特征取值有相当一部分是交叠的靠近边界的样本本身就是一个模糊地带。在这些边界样本里DPC的错误方式是某几个Versicolour样本因为离Virginica的高密度核心更近父节点指向了Virginica阵营顺着链条最终归入错误簇。这种情况换成K-Means也一样会发生因为边界点的归属本来就取决于特征空间中的距离权衡没有哪个聚类算法能保证100%正确。6. 参数敏感、算法边界和可扩展的改进方向6.1 dc是整个DPC里最敏感的参数把dc从2%依次调成0.5%、1%、5%、10%聚类结果会发生非常明显的漂移。dc太小每个点的邻居寥寥无几密度分布稀疏很多点ρ等于0无法形成有效的父节点链条dc太大密度变成全局平均计算出来的“聚类中心”全是数据集中部的大方块边界信息全部丢失。我的建议是先用2%作为默认值跑一遍决策图再把dc翻倍和减半各跑一遍如果三个结果在决策图上的中心点基本稳定说明数据对这个参数不敏感如果中心点大洗牌说明数据的簇结构本身就不稳固换参数只能自欺欺人。在论文或实验报告里最好把这种参数扫描的结果一并贴出来比只说一个最优值有说服力得多。6.2 原版DPC的几个硬伤原版DPC在鸢尾花上表现不错但它并不是万能的。第一个问题是密度层次如果数据里两个簇的密度差一个数量级低密度簇的成员会被当成噪声处理高密度簇则把低密度簇整体吞并。第二个问题是分配不可逆每个样本的父节点在最初就定了后续没有“纠正”环节一旦某个边界样本站错队以它为父节点的所有后代全部跟着错错误会沿着链条传播。第三个问题是边界密度阈值的设定只对两簇交界区域有意义如果簇的形状不规则边界识别会失真。另外不得不提的是复杂度距离矩阵要存N×N密度和距离计算也是O(N²)。拿到几万乃至几十万样本时内存和耗时都会变成现实问题。这不是DPC独有的毛病但确实是选型时必须考虑的因素。6.3 我建议的实操顺序和改进思路如果你要把DPC用到自己的数据上我建议按照这个顺序走一遍先标准化画距离分布直方图用2%分位数定dc画出决策图或γ排序图判断中心数量是否清晰再跑聚类并计算ARI等标签无关指标最后检查错误样本是不是集中在类簇边界。这套流程能让你在半小时内判断DPC到底适不适合当前数据。如果原版DPC的效果不够理想可以先尝试两个经典改进方向一是用K近邻距离替代dc来定义局部密度摆脱截断距离的全局敏感性也就是KNN-DPC的思路二是把一步到位的父节点分配改成迭代重分配允许边界样本在多次更新后修正归属。对于大多数中小规模数据这两招已经能解决很多原版DPC的痛点。跑完之后最大的体会是DPC不是银弹但它在“不知道K是多少”的场景下确实比K-Means从容得多。先跑一遍基础版再决定要不要上改进版。很多时候基本盘就够用了。