DPC算法实战：用MATLAB搞定密度峰值聚类（附完整代码）

DPC算法实战：用MATLAB搞定密度峰值聚类（附完整代码）

密度峰值聚类（DPC）算法自2014年提出以来，因其独特的聚类思路和简洁的实现方式，在数据科学领域获得了广泛关注。与传统的K-means等算法不同，DPC不需要预先指定聚类数量，而是通过寻找数据中的自然密度峰值来确定聚类中心。这种特性使其在处理复杂形状的数据分布时表现出色，尤其适合那些类簇形状不规则、密度不均匀的真实数据集。

对于MATLAB用户而言，DPC算法的实现既是一次挑战也是一次探索数据内在结构的机会。本文将带你从零开始，一步步实现完整的DPC算法流程，包括关键参数的选择技巧、核心代码的逐行解读，以及如何通过可视化手段直观地评估聚类效果。我们不仅会复现经典算法，还会分享几个提升聚类效果的实用技巧，帮助你将这个强大的工具应用到实际项目中。

1. 环境准备与数据加载

在开始DPC算法实现之前，我们需要确保MATLAB环境配置正确，并准备好合适的数据集。对于初次接触DPC的开发者，建议从二维或三维数据集开始，这样可以直观地验证算法的效果。

首先创建一个新的MATLAB脚本文件，命名为DPC_demo.m。我们将使用经典的合成数据集进行演示，这些数据集可以很好地展示DPC算法的特性：

% 生成合成测试数据 rng(42); % 设置随机种子保证可重复性 data1 = mvnrnd([1 1], eye(2)*0.2, 200); data2 = mvnrnd([4 1], eye(2)*0.3, 150); data3 = mvnrnd([2.5 4], eye(2)*0.25, 180); X = [data1; data2; data3]; % 可视化原始数据 figure; scatter(X(:,1), X(:,2), 15, 'filled'); title('原始数据分布'); xlabel('特征1'); ylabel('特征2'); box off;

这段代码生成了三个高斯分布的二维数据簇，分别位于不同的位置。运行后会显示一个散点图，帮助我们直观理解数据的原始分布情况。

提示：在实际项目中，如果使用自己的数据集，务必先进行标准化处理。不同特征间的量纲差异会影响距离计算，进而影响聚类效果。常用的标准化方法包括Z-score标准化和Min-Max归一化。

2. 核心算法实现

DPC算法的核心在于两个关键概念的计算：局部密度(ρ)和相对距离(δ)。我们将分步骤实现这些计算，并最终完成整个聚类流程。

2.1 距离矩阵计算

首先需要计算所有数据点之间的欧氏距离矩阵。这一步是算法中最耗时的部分，特别是当数据量较大时：

function D = calculateDistanceMatrix(X) % 计算欧氏距离矩阵 n = size(X, 1); D = zeros(n, n); for i = 1:n for j = i+1:n D(i,j) = norm(X(i,:) - X(j,:)); D(j,i) = D(i,j); end end end

对于大型数据集，可以考虑使用向量化计算或MATLAB的pdist2函数来提高效率。距离矩阵的计算结果将用于后续的局部密度和相对距离计算。

2.2 局部密度计算

局部密度ρ_i表示数据点i周围的数据密集程度。DPC论文中提出了两种计算方式：

1. 截断核(Truncated kernel)：统计距离小于截断距离dc的邻居数量
2. 高斯核(Gaussian kernel)：使用高斯函数加权计算所有邻居的贡献

以下是两种方法的MATLAB实现：

function rho = calculateLocalDensity(D, dc, method) n = size(D, 1); rho = zeros(n, 1); if strcmp(method, 'truncated') % 截断核方法 for i = 1:n rho(i) = sum(D(i,:) < dc) - 1; % 减去自身 end else % 高斯核方法 for i = 1:n rho(i) = sum(exp(-(D(i,:)/dc).^2)) - 1; % 减去自身贡献 end end end

截断核方法计算速度快，适合大规模数据集；高斯核方法更平滑，适合小规模或密度变化平缓的数据集。在实际应用中，可以根据数据特点选择合适的方法。

2.3 相对距离计算

相对距离δ_i表示数据点i到比它密度高的最近点的距离。对于密度最高的点，其δ值设为最大距离：

function delta = calculateRelativeDistance(D, rho) n = size(D, 1); delta = zeros(n, 1); [~, ordrho] = sort(rho, 'descend'); % 密度最高的点 delta(ordrho(1)) = max(D(ordrho(1),:)); % 其他点 for i = 2:n idx_higher_rho = ordrho(1:i-1); delta(ordrho(i)) = min(D(ordrho(i), idx_higher_rho)); end end

2.4 决策图与聚类中心选择

有了ρ和δ值后，我们可以绘制决策图来帮助选择聚类中心：

% 计算gamma值 = rho * delta gamma = rho .* delta; % 绘制决策图 figure; subplot(1,2,1); plot(rho, delta, 'o'); xlabel('\rho'); ylabel('\delta'); title('决策图'); subplot(1,2,2); plot(sort(gamma, 'descend'), 'r-o'); xlabel('样本点排序'); ylabel('\gamma = \rho * \delta'); title('Gamma值排序');

聚类中心通常是那些在决策图右上角的点，即同时具有高ρ值和高δ值的点。在实际操作中，可以通过观察γ值的下降趋势来确定聚类中心的数量。

3. 完整DPC函数实现

现在我们将上述步骤整合成一个完整的DPC函数：

function [rho, delta, gamma, cluster] = DPC(X, dc, method) % 计算距离矩阵 D = calculateDistanceMatrix(X); % 计算局部密度 rho = calculateLocalDensity(D, dc, method); % 计算相对距离 delta = calculateRelativeDistance(D, rho); % 计算gamma值 gamma = rho .* delta; % 选择聚类中心 [~, ordgamma] = sort(gamma, 'descend'); nCluster = input('请输入聚类中心数量: '); centers = ordgamma(1:nCluster); % 分配非中心点到最近的更高密度点所属的类 [~, ordrho] = sort(rho, 'descend'); cluster = zeros(size(X,1), 1); cluster(centers) = 1:nCluster; for i = 1:length(ordrho) if cluster(ordrho(i)) == 0 idx_higher_rho = find(rho > rho(ordrho(i))); [~, min_idx] = min(D(ordrho(i), idx_higher_rho)); cluster(ordrho(i)) = cluster(idx_higher_rho(min_idx)); end end end

这个函数接受数据矩阵X、截断距离dc和计算方法method作为输入，返回局部密度rho、相对距离delta、gamma值和最终的聚类标签cluster。

4. 结果可视化与评估

聚类完成后，我们需要评估结果的质量。对于二维数据，可以直接绘制聚类结果：

function ShowClusterResult(X, cluster) figure; gscatter(X(:,1), X(:,2), cluster); title('DPC聚类结果'); xlabel('特征1'); ylabel('特征2'); box off; end

对于更高维的数据，可以考虑使用降维技术（如PCA或t-SNE）来可视化聚类结果。

评估聚类质量的常用指标包括：

• 轮廓系数(Silhouette Coefficient)：衡量一个样本与同簇其他样本的相似度
• Calinski-Harabasz指数：簇间离散度与簇内离散度的比值
• Davies-Bouldin指数：簇间距离与簇内直径的比值

以下是计算轮廓系数的示例代码：

% 计算轮廓系数 silhouette_values = silhouette(X, cluster); avg_silhouette = mean(silhouette_values); disp(['平均轮廓系数: ', num2str(avg_silhouette)]);

5. 参数选择与调优技巧

DPC算法的性能很大程度上依赖于截断距离dc的选择。以下是几种确定dc的方法：

1. 经验法则：dc通常选择使平均每个点的邻居数为数据总量的1%-2%
2. k-距离图：绘制每个点到其第k近邻的距离的排序图，选择拐点处的距离值
3. 自适应方法：通过优化某些指标（如轮廓系数）来自动确定dc

以下是k-距离图的实现方法：

function PlotKneeDistance(D) % 计算每个点到第k近邻的距离 k = round(size(D,1)*0.02); % 取数据量的2% knn_dist = zeros(size(D,1),1); for i = 1:size(D,1) sorted_dist = sort(D(i,:)); knn_dist(i) = sorted_dist(k+1); % +1因为包含自身距离为0 end knn_dist = sort(knn_dist, 'descend'); % 绘制k-距离图 figure; plot(knn_dist, 'b-o'); xlabel('样本点排序'); ylabel(['到第', num2str(k), '近邻的距离']); title('k-距离图（用于选择dc）'); grid on; end

在实际应用中，我发现结合k-距离图和轮廓系数来选择dc通常能得到不错的效果。首先通过k-距离图确定dc的大致范围，然后在这个范围内微调dc值，选择使轮廓系数最大的那个。

6. 高级技巧与改进思路

虽然经典DPC算法已经表现不错，但在某些场景下仍有改进空间。以下是几个实用的改进方向：

1. 局部密度计算的改进：

   • 使用k近邻(KNN)思想计算局部密度
   • 引入共享近邻(SNN)概念提高密度估计的鲁棒性

2. 距离度量的改进：

   • 对于不同量纲的特征，使用马氏距离代替欧氏距离
   • 对于高维数据，考虑使用余弦相似度等更适合的度量

3. 自动确定聚类中心：

   • 通过分析γ值的下降趋势自动确定聚类数量
   • 使用统计方法（如峰谷检测）识别显著的聚类中心

4. 样本分配策略改进：

   • 使用KNN思想优化非中心点的分配过程
   • 引入概率分配机制处理边界点

以下是使用KNN改进局部密度计算的示例代码：

function rho = calculateLocalDensityKNN(D, k) n = size(D, 1); rho = zeros(n, 1); for i = 1:n [~, idx] = sort(D(i,:)); rho(i) = 1 / mean(D(i, idx(2:k+1))); % 取k近邻距离的倒数作为密度 end end

在实际项目中，根据数据特性和应用需求选择合适的改进策略非常重要。例如，在处理高维文本数据时，余弦相似度通常比欧氏距离更合适；而在处理密度差异很大的数据集时，KNN-based的局部密度计算方法往往效果更好。