密度峰值聚类CFDP算法实战：原理、Python实现与调参指南

CFDP（Clustering 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, metric='euclidean') 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), axis=1) - 1.0 return rho

注意减1.0是因为距离矩阵对角线上自己到自己的距离是0，exp(-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, dtype=int) 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_n=10): 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, dtype=int) 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, dtype=bool) 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, dc=None, center_indices=None, auto_center=True): dist_matrix = compute_distance_matrix(X) if dc is None: triu = dist_matrix[np.triu_indices(dist_matrix.shape[0], k=1)] 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²)存储，n=20000时就是3.2GB的float64，基本不可用。我在处理十万级样本时直接用Mini-Batch思路：先随机抽样做一个“骨架”CFDP，得到聚类中心，然后把其余样本分配给最近的中心或最近的高密度骨架点。这个近似方法效果不差，在图像特征聚类上实测稳定。

更优雅的做法是引入KNN优化。密度计算只考虑每个点最近的K个邻居，不需要全距离矩阵；δ计算也只需要和密度高于当前点的邻居比较，维护一个近邻堆就行。算法复杂度能降到O(nK log n)，K取20~50时内存占用直线下降。楼主如果后续要上大规模数据，建议直接KNN版起步。

5. 实战调参经验与后续扩展思路

5.1 不同场景下的参数参考表

我把在不同数据上试过的参数组合整理成一个速查表，可以作为起点再微调。

表里的分位数经验值都是基于“距离矩阵上三角”的。有几次我直接在原始特征上设d_c，效果一塌糊涂，后来统一改用距离矩阵分位数才稳定。

5.2 算法变体和我的使用体会

CFDP后续有很多变体值得关注。基于KNN局部密度的DPC变体解决了一部分参数敏感问题；把δ替换为“到密度更高点的最短路径长度”的版本能处理流形数据；还有与密度聚类结合的层次化方法，能输出聚类层次而不是单一划分。

我在实际项目中用得最多的是给CFDP接一个“粗聚类+细分类”的流程：先用CFDP找出几个大簇，然后在每个簇内部再跑一次CFDP或者KMeans，这样能在无监督情况下识别出子结构，而且比直接指定簇数更自然。

5.3 最后聊一点个人心得

用CFDP这几年，最深的体会是“聚类结果要先看决策图，再谈指标”。很多同学一上来就调参、跑指标，忽略了决策图本身就是最直观的质量诊断工具。如果决策图右上角是模糊一片，说明数据没有清晰的聚类结构，再怎么调参也白搭。另外，CFDP虽然参数少，但每个参数的背后都有明确物理含义，理解和调整起来比KMeans那种纯数学优化直观得多。

如果只是需要一个快速聚类工具，现在很多库也有成熟实现，但我建议至少手动实现一次，真正理解密度和距离的关系以后，你用任何聚类库都会更有底气。