CFDPClustering by fast search and find of density peaks这个算法中文圈子里一般叫密度峰值聚类2014年发表在Science上。我在实际项目里用过不少聚类方法KMeans、DBSCAN、层次聚类都折腾过但遇到簇形状不规则、密度差异大的数据时CFDP往往是最先能给出合理结果的那一个。这篇不打算复述论文而是从一个实现者的角度把算法拆开讲清楚它到底在做什么、代码怎么写、参数怎么调以及我踩过的那些坑。适合正在做无监督学习、数据挖掘或者想找一个能处理任意形状簇聚类的同学参考。1. 先想清楚CFDP在干什么密度、距离和决策图1.1 聚类中心的两条判定规则CFDP的核心思路很简单就两条规则。第一聚类中心本身一定是局部密度极高的点也就是它周围邻居足够密。第二这个中心和比它密度更高的其他点要离得够远换句话说真正好的聚类中心是被一群低密度点包围的身边没有更高密度的点紧贴着。这两条规则听起来朴素但实现出来效果意外地好。KMeans要预先指定K对簇形状还敏感DBSCAN靠epsilon和minPts两个参数密度不均时经常一个参数管不了全局。CFDP不需要预知簇数量聚类中心可以通过决策图人工选或者用阈值自动选。它把“聚类中心”这个概念落到了一个可量化的二维坐标上——每个样本都有两个属性局部密度ρ和到更高密度点的最小距离δ。这里我打个比方。你把每个样本想象成一座山上的登山者ρ就是这个人位置的人气周围围了多少人。δ就是他抬眼望过去离自己最近的、人气比自己高的人隔了多远。真正的“山大王”两边都突出人气爆棚而且周围没有人气比他更高的邻居只有一群比他矮的仰望着他。而那些人气一般但δ很大的点通常是离群点孤零零地待着旁边没有人气更高的人所以δ也大但ρ不行。1.2 局部密度ρ和距离δ的数学定义密度ρ的定义有两种方式对应两种核函数。一种是截断核ρ_i Σ_j χ(d_ij - d_c)其中χ(x) 1当x 0否则为0。大白话就是数一数有多少个样本到点i的距离小于截止距离d_c这个数量就是密度。另一种是高斯核ρ_i Σ_j exp(-(d_ij / d_c)^2)。截断核是硬计数高斯核是软加权。我一开始图省事直接用截断核后来发现它会导致大量样本的ρ值相等——因为距离小于d_c的样本数是离散整数尤其在小数据集上一大半点的ρ都是同一个数。而ρ一相等后面算δ就要额外处理并列情况非常麻烦。高斯核平滑得多每个点即使距离很近也会有不同的ρ值代码逻辑可以少一个分支。δ的定义稍微绕一点。对每个点i找到所有密度比它高的点取这些点中到i最近的那个距离就是δ_i。对于密度最高的那个点它没有“更高密度的点”可以找惯例上直接取它到所有其他点的最大距离保证它在决策图里稳坐右上角。这套逻辑把聚类中心的选择从“看数据分布”变成了“看决策图”本质上是一个降维和可视化问题。每次我跑完CFDP第一步都是画ρ-δ散点图看右上角有几个点孤立地飘着那些就是候选中心。1.3 决策图怎么读中心怎么定决策图的横轴是ρ纵轴是δ。两类点会出现在右上角一类是真中心ρ高、δ也高另一类是离群点ρ低但δ极高。区别就在ρ上所以选中心时不要只看δ大的要盯着右上角——“又高又右边”的区域。实际选中心时我一般先画出决策图人工数一下右上角有几个明显的点然后在代码里把它们的索引挑出来。如果想自动选可以用γ ρ × δ这个标量γ越大越像中心把γ排序后从大到小看通常曲线会有一个明显的“肘部”肘部之前的就是聚类中心。这个思路和找KMeans的K值时看肘部图一样不完美但够用。2. 动手前必须搞清的三个关键细节2.1 截断距离d_c到底怎么选d_c是整个算法里最敏感的参数没有之一。论文给的经验是取所有样本两两距离排序后的第1%到2%分位数我实测下来这个范围对大多数二维数据集都还行但三维以上的数据就不一定了尤其是尺度差异大或者存在高维稀疏特征时需要重新调。我自己的习惯是先算所有样本间的距离矩阵用np.percentile(distances, 1.5)得到一个初始值然后跑一遍聚类看决策图和分配结果如果簇数量明显不合预期就在1%到5%之间二分搜索。注意d_c选的过大密度全都差不多高中心选不出来选得过小每个点都是孤立峰直接变成离群点大杂烩。还有一个小技巧如果数据各维度的物理意义和量纲差很多一定要先做标准化。不然欧式距离被数值大的维度主导d_c等于白设。2.2 密度并列时δ的计算陷阱截断核会带来一个经典问题密度并列。多个点ρ值相同排序时它们的相对顺序会影响后续δ的分配。我最初实现时直接按ρ降序排密度相同的点之间会互相“找”导致δ计算出现随机性聚类结果每次都略微不同。解决办法有两个。最简单的就是改用高斯核ρ值几乎不会并列。如果坚持用截断核就需要在排序时增加一个次要关键字比如按原始索引升序保证结果可复现。我建议新手直接上高斯核省掉这个麻烦。另一个细节是找“密度更高的点”时严格比较是ρ_j ρ_i不是≥。如果用了≥就变成密度最高的点也把自己算进“更高密度集合”δ会变成0整个决策图就废了。用截断核时尤其容易踩这个坑。2.3 噪声点和边界点的密度阈值策略CFDP本身不是为噪声设计的好手但它论文里给了一个精巧的边界策略。核心思想每个簇的边界点就是那些到其他簇样本距离小于d_c的点这些边界点的密度最大值作为该簇的密度阈值。分配完成后密度低于这个阈值的样本被标记为噪声。我在实现时是按簇分别算阈值的因为各簇密度差异很大全局一个阈值不公平。具体做法是遍历所有样本对如果距离小于d_c且标签不同那这两个点都算边界点然后每个簇取这些边界点中密度最大的那个值作为阈值。最后凡是某簇内部密度低于阈值的点统统置为噪声标签-1。这个方法有一个弱点是簇重叠区域宽时阈值会被拉得偏高把正常边缘样本误杀。所以我一般把噪声标记当成参考真正做性能评估时对比有噪声和没噪声两种版本的聚类指标。3. Python实现从距离矩阵到聚类结果3.1 距离矩阵计算与内存预判CFDP第一步就是算全样本的距离矩阵n个样本就是n×n的矩阵。n在1万以内都能忍超过1万就得考虑分块或者KNN近似了否则光矩阵就占几百MB甚至几个GB。我用SciPy的pdist算成对距离再转成方阵。注意pdist返回的是扁平数组顺序是按行展开的squareform能直接还原成方阵。import numpy as np from scipy.spatial.distance import pdist, squareform def compute_distance_matrix(X): # X: shape (n_samples, n_features) dist_array pdist(X, metriceuclidean) dist_matrix squareform(dist_array) np.fill_diagonal(dist_matrix, 0.0) return dist_matrix这里的fill_diagonal不是必须的因为pdist默认对角线是0方阵化之后也能对上但为了防止浮点误差把对角变成长尾小值我习惯手动归零。3.2 密度ρ和最小距离δ的实现密度计算很简单高斯核版本一行搞定。d_c我一般取距离矩阵上三角部分排除0的排序后一个分位数。def compute_rho(dist_matrix, dc): n dist_matrix.shape[0] # 高斯核密度 rho np.sum(np.exp(-(dist_matrix / dc) ** 2), axis1) - 1.0 return rho注意减1.0是因为距离矩阵对角线上自己到自己的距离是0exp(-0)是1要抠掉自身。δ的计算稍微讲究。我先按ρ降序得到索引然后从第二个开始遍历当前点的“高密度邻居集合”就是排在它前面所有点取它们与当前点的距离最小值。def compute_delta(dist_matrix, rho): n dist_matrix.shape[0] rho_order np.argsort(-rho) # 密度从高到低 delta np.zeros(n) nearest_higher np.full(n, -1, dtypeint) max_rho_idx rho_order[0] delta[max_rho_idx] np.max(dist_matrix[max_rho_idx]) nearest_higher[max_rho_idx] -1 for i in range(1, n): idx rho_order[i] higher_indices rho_order[:i] dist_to_higher dist_matrix[idx, higher_indices] min_pos np.argmin(dist_to_higher) delta[idx] dist_to_higher[min_pos] nearest_higher[idx] higher_indices[min_pos] return delta, nearest_higher这个循环天然是O(n²)每层循环最多n次。n大于5000时这个纯Python循环会有点慢可以向量化预先按ρ排序重排距离矩阵然后用np.min在切片上批量算但代码可读性会变差我一般在n小于1万时就用这个写法再大就考虑优化。3.3 绘制决策图并选定聚类中心决策图就是ρ和δ的散点。我通常把ρ放横轴、δ放纵轴然后手动观察右上角。为了方便我会算γ ρ * δ把γ排序打印出来人工看拐点。def plot_decision_graph(rho, delta, top_n10): gamma rho * delta sorted_idx np.argsort(-gamma) print(Top candidates (index, rho, delta, gamma):) for i in range(top_n): idx sorted_idx[i] print(idx, round(rho[idx], 4), round(delta[idx], 4), round(gamma[idx], 4)) # 画散点 import matplotlib.pyplot as plt plt.scatter(rho, delta) plt.xlabel(rho) plt.ylabel(delta) plt.show()选定中心的方法从sorted_idx里取前K个。K就是你在决策图上数出来的中心个数。我见过有人用二阶差分自动找拐点效果取决于数据的清晰程度噪声大时容易多选或少选还是人工看决策图最可靠。3.4 样本分配与噪声标记的完整流程中心定了之后剩下所有样本按密度从高到低依次归入“最近更高密度点”所在的簇。这个策略是CFDP的精华它等于把每个非中心的样本沿着“密度上升路径”一路挂到最近的峰值上。def assign_labels(dist_matrix, rho, nearest_higher, centers): n dist_matrix.shape[0] labels np.full(n, -1, dtypeint) for c_id, c_idx in enumerate(centers): labels[c_idx] c_id rho_order np.argsort(-rho) for idx in rho_order: if labels[idx] ! -1: continue labels[idx] labels[nearest_higher[idx]] return labels理论上nearest_higher指向的节点一定比当前点密度高而且由于密度关系有向无环最终一定指向某个中心节点。所以按密度降序分配能保证nearest_higher的标签已经在之前分配好。这里有个隐含假设密度相同的点之间不能互为nearest_higher否则可能死循环。高斯核的密度基本是浮点数不会完全相等所以安全。然后做噪声标记用边界点策略。def label_noise(dist_matrix, rho, labels, dc): n dist_matrix.shape[0] boundary_rho_max {} for i in range(n): for j in range(i 1, n): if dist_matrix[i, j] dc and labels[i] ! labels[j] and labels[i] ! -1 and labels[j] ! -1: rho_bound min(rho[i], rho[j]) if labels[i] not in boundary_rho_max: boundary_rho_max[labels[i]] rho_bound else: boundary_rho_max[labels[i]] max(boundary_rho_max[labels[i]], rho_bound) if labels[j] not in boundary_rho_max: boundary_rho_max[labels[j]] rho_bound else: boundary_rho_max[labels[j]] max(boundary_rho_max[labels[j]], rho_bound) noise_mask np.zeros(n, dtypebool) for idx in range(n): lab labels[idx] if lab -1: continue if lab in boundary_rho_max and rho[idx] boundary_rho_max[lab]: noise_mask[idx] True result labels.copy() result[noise_mask] -1 return result噪声标记这段代码是O(n²)的双层循环n一大地狱就来了。我实际用的优化版是先构建KDTree只查d_c半径内的邻居能快一个数量级。数据集不大时用这个直白版没问题正好方便理解算法流程。完整调用入口def cfdp_clustering(X, dcNone, center_indicesNone, auto_centerTrue): dist_matrix compute_distance_matrix(X) if dc is None: triu dist_matrix[np.triu_indices(dist_matrix.shape[0], k1)] dc np.percentile(triu, 1.5) rho compute_rho(dist_matrix, dc) delta, nearest_higher compute_delta(dist_matrix, rho) if center_indices is None: # 人工/自动选择中心返回索引列表 gamma rho * delta # 简化自动选gamma前5或者按拐点这里示意取排序后前3 center_indices list(np.argsort(-gamma)[:3]) labels assign_labels(dist_matrix, rho, nearest_higher, center_indices) labels label_noise(dist_matrix, rho, labels, dc) return labels, rho, delta, dc如果你有聚类中心的人工标注比如先画决策图再回来填索引把center_indices传进去就行。4. 常见问题与排查方法实录4.1 聚类中心“多出来”几个怎么办这是我最常遇到的问题。决策图右上角不止有真中心还会有几个离群点它们δ很大但ρ很小。把它们选进去之后簇的数量就虚高了每个离群点成了独立簇。排查方法是看γ排名后面几个候选的ρ值。真中心的ρ一定显著大于周围邻居离群点的ρ通常就在整个数据集的低分位段。我通常会打印候选点的ρ值凡是ρ小于整体密度中位数的直接排除掉。如果自动选中心就用这个规则做硬过滤。但也要注意一种特殊情况小簇的真中心ρ可能确实不高。如果数据集里同时存在极密的大簇和很稀疏的小簇全局中位数阈值会误删小簇中心。这种时候建议画一下密度直方图再决定阈值放哪里。4.2 d_c的微小变化导致聚类结果剧烈变化d_c对截断核尤其致命。我试过同一个数据集d_c从1.2改到1.5聚类簇数从4跳到7。根本原因是截断核的密度计算是离散跳变的d_c一过某个阈值好多点的ρ同时变排序瞬间重排δ也跟着全变。应对策略就是改用高斯核它连续可导d_c的微小变化只会带来ρ的平滑变化聚类结果稳定很多。另外我强烈建议对d_c做一次小网格搜索比如1%、1.5%、2%、3%分位数分别跑一遍看决策图和簇数是否稳定。如果结果在某个范围内变化很大说明数据本身对密度尺度敏感需要先检查是否做了标准化或者是否有离群特征干扰。4.3 高维数据效果变差怎么办维度升高后所有样本两两距离都趋向于差不多大这是高维欧式距离的普遍问题。CFDP依赖距离排序和密度峰值遇到这个情况会直接退化ρ差距小δ随机决策图上没有明显峰值。我的经验是先用t-SNE或UMAP把维度压到2~5维再跑CFDP。这样既能保留局部结构又让距离度量恢复分辨力。但注意降维本身也有信息损失如果业务上不能接受就试试换距离度量比如余弦距离或者基于相关性的距离。另外也可以用KNN距离作为核心距离只考虑每个点前K个近邻的密度贡献能缓解高维稀疏问题。4.4 样本量大时内存和速度双双告急距离矩阵是O(n²)存储n20000时就是3.2GB的float64基本不可用。我在处理十万级样本时直接用Mini-Batch思路先随机抽样做一个“骨架”CFDP得到聚类中心然后把其余样本分配给最近的中心或最近的高密度骨架点。这个近似方法效果不差在图像特征聚类上实测稳定。更优雅的做法是引入KNN优化。密度计算只考虑每个点最近的K个邻居不需要全距离矩阵δ计算也只需要和密度高于当前点的邻居比较维护一个近邻堆就行。算法复杂度能降到O(nK log n)K取20~50时内存占用直线下降。楼主如果后续要上大规模数据建议直接KNN版起步。5. 实战调参经验与后续扩展思路5.1 不同场景下的参数参考表我把在不同数据上试过的参数组合整理成一个速查表可以作为起点再微调。数据特点距离度量d_c分位数密度核建议二维稠密点云欧式1%~2%高斯核标准化后直接跑二维含噪声欧式2%~3%高斯核配合边界噪声过滤高维文本/图像特征余弦1%按余弦距离高斯核先降维或只用KNN近邻密度密度差异极大的簇欧式3%~5%高斯核考虑分簇分别设置d_c稀疏高维曼哈顿2%截断核谨慎建议降维表里的分位数经验值都是基于“距离矩阵上三角”的。有几次我直接在原始特征上设d_c效果一塌糊涂后来统一改用距离矩阵分位数才稳定。5.2 算法变体和我的使用体会CFDP后续有很多变体值得关注。基于KNN局部密度的DPC变体解决了一部分参数敏感问题把δ替换为“到密度更高点的最短路径长度”的版本能处理流形数据还有与密度聚类结合的层次化方法能输出聚类层次而不是单一划分。我在实际项目中用得最多的是给CFDP接一个“粗聚类细分类”的流程先用CFDP找出几个大簇然后在每个簇内部再跑一次CFDP或者KMeans这样能在无监督情况下识别出子结构而且比直接指定簇数更自然。5.3 最后聊一点个人心得用CFDP这几年最深的体会是“聚类结果要先看决策图再谈指标”。很多同学一上来就调参、跑指标忽略了决策图本身就是最直观的质量诊断工具。如果决策图右上角是模糊一片说明数据没有清晰的聚类结构再怎么调参也白搭。另外CFDP虽然参数少但每个参数的背后都有明确物理含义理解和调整起来比KMeans那种纯数学优化直观得多。如果只是需要一个快速聚类工具现在很多库也有成熟实现但我建议至少手动实现一次真正理解密度和距离的关系以后你用任何聚类库都会更有底气。