告别K-Means！用Python手撸Science上的DPC算法，搞定任意形状数据聚类

密度峰值聚类DPC用Python突破传统K-Means的局限当面对螺旋形、环形或交叉分布的数据集时许多数据科学从业者都有过这样的经历反复调整K-Means参数却始终无法获得理想的聚类效果。这正是2014年发表在《Science》上的密度峰值聚类算法(DPC)要解决的核心问题。与传统方法不同DPC通过创新的密度峰值概念能够自动识别任意形状数据中的自然聚类中心无需预先指定簇数量。1. 为什么需要超越K-MeansK-Means算法自1967年提出以来因其简单高效成为最广泛使用的聚类工具。但它存在两个根本性局限球形分布假设K-Means基于样本与簇中心的欧式距离进行划分隐含假设数据呈球形分布。这在处理螺旋形、月牙形等复杂结构时效果显著下降。K值依赖需要人工指定聚类数量K而实际应用中K往往未知。错误的K值会导致完全偏离真实结构的聚类结果。# 典型K-Means在复杂数据上的失败案例 from sklearn.cluster import KMeans import matplotlib.pyplot as plt # 生成螺旋数据集 theta np.linspace(0, 4*np.pi, 200) r np.linspace(0.5, 1, 200) x r * np.cos(theta) y r * np.sin(theta) spiral np.column_stack([x, y]) # K-Means聚类 kmeans KMeans(n_clusters3).fit(spiral) plt.scatter(spiral[:,0], spiral[:,1], ckmeans.labels_) plt.title(K-Means在螺旋数据上的表现) plt.show()相比之下DPC算法具有三大突破性优势形状无关性基于局部密度而非距离可识别任意形状的簇自动中心发现通过密度峰值自动确定聚类中心无需预设K值离群点鲁棒性对噪声数据不敏感能自然区分核心点与边缘点2. DPC算法核心原理拆解DPC建立在两个直观而强大的假设基础上假设一聚类中心点的局部密度应显著高于周围邻居假设二聚类中心应与更高密度的点保持较远距离这两个假设转化为两个核心计算指标2.1 局部密度(ρ)计算局部密度衡量数据点周围的样本密集程度DPC提供两种计算方式计算方式公式适用场景截断核ρᵢ ∑ⱼχ(dᵢⱼ - d_c)离散型数据高斯核ρᵢ ∑ⱼexp(-(dᵢⱼ/d_c)²)连续型数据其中d_c是关键参数——截断距离通常选择使每个点平均有1%-2%的邻居落在d_c范围内。def compute_density(dists, dc, methodgaussian): 计算局部密度 if method cutoff: rho np.sum(dists dc, axis1) - 1 # 排除自身 else: # gaussian rho np.sum(np.exp(-(dists/dc)**2), axis1) - 1 return rho2.2 相对距离(δ)计算对于每个点iδᵢ定义为若是最高密度点δᵢ max(dᵢⱼ)其他点δᵢ min(dᵢⱼ) where ρⱼ ρᵢdef compute_delta(dists, rho): 计算相对距离 N len(rho) delta np.zeros(N) nearest_higher np.zeros(N, dtypeint) # 按密度降序排序索引 rho_order np.argsort(-rho) for i, idx in enumerate(rho_order): if i 0: # 最高密度点 delta[idx] np.max(dists[idx]) nearest_higher[idx] -1 else: # 找出密度更高的点 higher_mask rho rho[idx] if np.any(higher_mask): delta[idx] np.min(dists[idx, higher_mask]) nearest_higher[idx] np.argmin(dists[idx, higher_mask]) else: delta[idx] np.max(dists[idx]) nearest_higher[idx] -1 return delta, nearest_higher2.3 决策图与中心选择将ρ和δ绘制为决策图(Decision Graph)聚类中心表现为右上角的离群点——同时具有高ρ和高δ。实际应用中可以设定ρ和δ的双阈值或选择ρ×δ乘积最大的前K个点作为中心。3. 完整DPC实现与可视化下面我们实现完整的DPC流程并在螺旋数据集上进行测试class DensityPeakClustering: def __init__(self, dc_percent2.0): self.dc_percent dc_percent def fit(self, X): # 计算距离矩阵 self.dists np.sqrt(((X[:, None] - X) ** 2).sum(axis2)) # 确定截断距离dc N X.shape[0] position int(N * (N - 1) * self.dc_percent / 100) self.dc np.sort(self.dists.ravel())[position N] # 计算局部密度 self.rho compute_density(self.dists, self.dc) # 计算相对距离 self.delta, self.nearest_higher compute_delta(self.dists, self.rho) return self def find_centers(self, autoTrue, KNone): if auto: # 自动阈值法 rho_thresh (np.min(self.rho) np.max(self.rho)) / 2 delta_thresh (np.min(self.delta) np.max(self.delta)) / 2 self.centers np.where((self.rho rho_thresh) (self.delta delta_thresh))[0] else: # 指定K值法 product self.rho * self.delta self.centers np.argsort(-product)[:K] return self.centers def assign_clusters(self): K len(self.centers) self.labels -np.ones(len(self.rho), dtypeint) # 标记中心点 for i, center in enumerate(self.centers): self.labels[center] i # 按密度降序分配标签 for idx in np.argsort(-self.rho): if self.labels[idx] -1: self.labels[idx] self.labels[self.nearest_higher[idx]] return self.labels可视化展示# 生成测试数据 from sklearn.datasets import make_moons X, _ make_moons(n_samples300, noise0.05) # DPC聚类 dpc DensityPeakClustering(dc_percent2).fit(X) centers dpc.find_centers(autoTrue) labels dpc.assign_clusters() # 绘制结果 plt.figure(figsize(12,5)) plt.subplot(121) plt.scatter(dpc.rho, dpc.delta, ck, s10) plt.scatter(dpc.rho[centers], dpc.delta[centers], cr, s50) plt.xlabel(Density (ρ)) plt.ylabel(Delta (δ)) plt.title(Decision Graph) plt.subplot(122) plt.scatter(X[:,0], X[:,1], clabels, cmapviridis, s20) plt.scatter(X[centers,0], X[centers,1], cred, markerx, s100) plt.title(Clustering Result) plt.show()4. 实战优化与参数调优虽然DPC理论优雅但实际应用中需要注意几个关键点4.1 截断距离dc的选择dc直接影响密度计算结果建议采用以下策略百分比法保持1%-2%的邻居比例默认方法k近邻法取每个点到第k近邻的距离中位数网格搜索在候选dc值中选择使聚类结果最稳定的值def find_optimal_dc(dists, percent_range(1,3)): 寻找最佳dc值 N dists.shape[0] results [] for percent in np.linspace(*percent_range, 10): position int(N * (N - 1) * percent / 100) dc np.sort(dists.ravel())[position N] rho compute_density(dists, dc) delta, _ compute_delta(dists, rho) product rho * delta # 使用轮廓系数评估稳定性 centers np.argsort(-product)[:3] labels assign_clusters(rho, centers, nearest_higher) if len(np.unique(labels)) 1: score silhouette_score(dists, labels, metricprecomputed) results.append((percent, dc, score)) return max(results, keylambda x: x[2])4.2 处理大规模数据原始DPC的O(N²)复杂度限制了其在大型数据集的应用可通过以下方法优化子采样先在小样本上确定dc和中心再扩展到全数据集KD树/球树加速最近邻搜索分布式计算将距离矩阵计算分布到多台机器from sklearn.neighbors import KDTree def fast_compute_density(X, dc, k30): 使用KDTree加速密度计算 tree KDTree(X) counts tree.query_radius(X, rdc, count_onlyTrue) return counts - 1 # 排除自身4.3 多维度数据预处理高维数据中距离度量可能失效建议使用PCA等降维技术保留主要特征采用马氏距离考虑特征相关性对连续和离散特征分别处理from sklearn.decomposition import PCA # 高维数据预处理 pca PCA(n_components0.95) X_reduced pca.fit_transform(X_highdim) dpc DensityPeakClustering().fit(X_reduced)5. DPC进阶应用场景DPC的特殊优势使其在多个领域大放异彩5.1 异常检测高δ低ρ的点往往是异常值这种方法比传统阈值法更可靠def find_anomalies(rho, delta, alpha0.1): 基于DPC的异常检测 rho_norm (rho - np.min(rho)) / (np.max(rho) - np.min(rho)) delta_norm (delta - np.min(delta)) / (np.max(delta) - np.min(delta)) anomaly_scores delta_norm - rho_norm threshold np.percentile(anomaly_scores, 100*(1-alpha)) return anomaly_scores threshold5.2 层次聚类扩展通过调整dc值可以自然获得层次聚类结构从小dc开始识别细粒度簇逐渐增大dc观察簇的合并过程构建聚类树状图5.3 图像分割应用将像素的坐标和颜色作为特征DPC可实现自动图像分割def image_segmentation(image, dc_percent0.5): 基于DPC的图像分割 h, w image.shape[:2] coords np.indices((h, w)).transpose(1,2,0).reshape(-1,2) colors image.reshape(-1, 3) features np.hstack([coords/10, colors]) # 坐标权重降低 dpc DensityPeakClustering(dc_percent).fit(features) centers dpc.find_centers(autoTrue) labels dpc.assign_clusters() return labels.reshape(h, w)在实际项目中DPC算法特别适合处理那些传统方法难以应对的复杂结构数据。我曾在一个客户行为分析项目中使用DPC成功识别出5种独特的用户群体这些群体在传统的RFM模型中被混为一谈。算法的自动中心发现功能帮助我们发现了两个从未考虑过的重要客户细分。