★ FEATURED ARTICLE 简介这份MATLAB聚类算法实践包面向机器学习初学者与需要快速上手无监督聚类的开发者围绕k-means、DBSCAN与基于密度峰值的ADPC三种算法展开帮助读者在MATLAB环境中完成数据分组与簇结构探索。压缩包共23个文件约81KB以14个m脚本为核心辅以5个mat数据文件、3个txt样本与1个md说明分别承担算法实现、数据集加载与运行指引等用途。资源按算法分目录组织k-means、dbscan与adpc各自独立便于对照阅读与逐模块调试。目前已有183人学习。读者可借此理解k-means的迭代中心更新、DBSCAN的核心对象与噪声判定以及ADPC通过密度峰值确定簇中心的思路并直接替换自带数据集验证效果是入门聚类算法、对比不同方法特性的轻量实践素材。1. 三种聚类算法放在同一份 MATLAB 工程里到底该怎么跑很多人第一次拿到「k-means、DBSCAN、密度峰值聚类」打包在一起的 MATLAB 实践文件第一反应是双击脚本看结果然后被一堆散落的.m文件和.mat数据搞懵。这个标题真正要解决的问题不是「聚类是什么」而是在同一套数据上把三种思路完全不同的聚类算法跑通、对比、并知道各自什么时候会翻车。k-means 靠距离均值切分DBSCAN 靠密度连通性找任意形状簇密度峰值聚类DPC靠局部密度和最近高密度距离画决策图选中心点。三者适合谁做图像分割、点云分组、用户分群的工程师只要手头有 MATLAB 和一份二维或低维特征数据就能直接复现。下面按「先立住原理、再动手跑、最后避坑」的顺序拆开讲所有代码都是最小可运行版本参数含义和调法一并给出。2. 三种聚类算法的核心差异与选型依据2.1 k-means 的适用边界球形簇、预设 K、对初值敏感k-means 的目标函数是最小化簇内平方误差和数学上等价于在每轮迭代中把样本分配给最近质心再更新质心为簇内均值。它的隐含假设是簇呈凸形且大小相近一旦数据是长条形、环形或密度差异大结果就会强行把相邻簇切开。MATLAB 里最直接的调用是kmeans函数但很多人不知道它默认用 k-means 初始化重复次数Replicates默认是 1数据稍复杂就容易陷局部最优。选型判断很简单如果你能接受「簇是球状、K 值靠肘部法或轮廓系数确定」k-means 是速度最快的基线。它的计算复杂度是 O(nKd)每轮迭代n 是样本数K 是簇数d 是维度。对于几万条二维数据MATLAB 单次运行通常在秒级。2.2 DBSCAN 的密度定义eps 与 MinPts 如何决定簇的生死DBSCAN 不预设簇数而是通过两个参数——邻域半径eps和最小点数MinPts——把样本分为核心点、边界点和噪声点。核心点的 eps 邻域内至少包含 MinPts 个样本所有密度可达的核心点连成同一个簇。它的优势是能发现任意形状簇并识别噪声代价是对参数极其敏感eps 稍大两个邻近簇会合并eps 稍小大量样本被判为噪声。MATLAB 没有内置 DBSCAN 函数截至常见版本通常做法是自己写距离矩阵加区域查询或者用knnsearch加速。判断 eps 的一个实用方法是画 k-距离曲线对每个点计算到第 k 近邻的距离升序排列后找拐点拐点对应的距离就是 eps 的合理起点。MinPts 一般取维度加一或维度乘二二维数据常用 4 到 10。2.3 密度峰值聚类的决策图局部密度 ρ 与距离 δ 的联合判据密度峰值聚类DPC的核心思想是簇中心同时具有较高的局部密度和较大的到更高密度点的最近距离。对每个样本 i计算局部密度 ρ_i常用截断核或高斯核和 δ_i到任何密度更高点的最小距离。把所有点画在 ρ-δ 决策图上右上角的离群点就是簇中心。确定中心后其余点沿密度梯度分配给最近的更高密度点所属簇。DPC 的优点是能识别非球形簇且不需要迭代优化缺点是截断距离dc的选取直接影响 ρ 的分布。经验做法是让平均邻居数占样本总数的 1% 到 2%即dc取所有成对距离升序排列的第 1% 到 2% 分位数。MATLAB 实现时距离矩阵是 O(n²) 内存样本超过一万就要考虑稀疏化或分块计算。2.4 三种算法在同一份数据上的选型对照维度k-meansDBSCAN密度峰值聚类是否需要预设簇数是否否看决策图能发现的簇形状凸形/球形任意形状任意形状对噪声的处理全部分配无噪声概念显式标记噪声低密度点归入远端簇主要敏感参数K、初始化eps、MinPtsdc、中心选取计算复杂度O(nKd) 每轮O(n log n) 带索引O(n²) 距离矩阵适合场景快速基线、大样本密度不均、含噪簇中心明显、中小样本这张表不是让你背而是遇到新数据时先问三个问题簇数知不知道形状是不是球形噪声多不多答案决定你从哪个算法开始试。3. 在 MATLAB 里把三种算法跑通的最小工程3.1 生成可复现的二维测试数据不要用随机数直接跑否则每次结果不同没法对比。下面这段代码生成三簇高斯数据加一圈环形数据再混入均匀噪声覆盖球形和非球形两种情形。% 生成测试数据三簇高斯 一圈环形 噪声 rng(42); % 固定随机种子保证可复现 n1 200; n2 200; n3 200; % 三簇高斯均值分开 X1 randn(n1,2) * 0.6 [0, 0]; X2 randn(n2,2) * 0.6 [5, 5]; X3 randn(n3,2) * 0.6 [0, 6]; % 环形数据半径 4角度均匀 theta linspace(0, 2*pi, 300); X_ring [4*cos(theta), 4*sin(theta)] randn(300,2)*0.15; % 均匀噪声 X_noise rand(80,2) * 10 - 2; X [X1; X2; X3; X_ring; X_noise]; % 标准化到零均值单位方差避免量纲影响距离 X (X - mean(X)) ./ std(X); figure; scatter(X(:,1), X(:,2), 10, filled); title(原始数据分布);逻辑说明rng(42)锁定随机种子保证你和我跑出的图一致。标准化这一步很多人跳过但只要两列特征量纲不同距离计算就会被大量纲主导。参数上高斯簇的标准差 0.6 控制簇的松散程度环形半径 4 决定它离高斯簇多远噪声比例约 10% 用来检验 DBSCAN 的噪声识别能力。3.2 k-means 的肘部法与轮廓系数确定 K直接跑kmeans(X, 3)只是猜 K。正确做法是遍历 K 从 2 到 8记录簇内平方和和轮廓系数找拐点。K_range 2:8; sse zeros(size(K_range)); sil zeros(size(K_range)); for i 1:length(K_range) K K_range(i); [idx, C, sumd] kmeans(X, K, Replicates, 5, Distance, sqeuclidean); sse(i) sum(sumd); sil(i) mean(silhouette(X, idx)); end figure; subplot(1,2,1); plot(K_range, sse, -o); xlabel(K); ylabel(SSE); title(肘部法); subplot(1,2,2); plot(K_range, sil, -o); xlabel(K); ylabel(轮廓系数); title(轮廓系数);逻辑说明Replicates设为 5 表示用 5 次不同初始化取最优避免单次陷入局部最优。Distance选sqeuclidean与目标函数一致。看结果时SSE 曲线拐点对应的 K 和轮廓系数峰值对应的 K 如果一致就可以确定不一致时优先看轮廓系数因为它衡量的是簇间分离度。参数上silhouette对大规模数据计算慢样本超过两万建议抽样计算。3.3 手写 DBSCAN距离矩阵加区域查询MATLAB 没有内置 DBSCAN下面是一个紧凑实现核心是regionQuery找邻域再用队列扩展簇。function [idx, isNoise] dbscan_simple(X, eps, MinPts) n size(X,1); D pdist2(X, X); % 成对距离矩阵 visited false(n,1); idx zeros(n,1); % 0 表示未分类 clusterId 0; for i 1:n if visited(i), continue; end visited(i) true; neighbors find(D(i,:) eps); if numel(neighbors) MinPts idx(i) -1; % 暂时标为噪声 else clusterId clusterId 1; idx expandCluster(i, neighbors, clusterId, D, eps, MinPts, idx, visited); end end isNoise idx -1; end function idx expandCluster(~, neighbors, cid, D, eps, MinPts, idx, visited) idx(neighbors(1)) cid; k 1; while k numel(neighbors) j neighbors(k); if ~visited(j) visited(j) true; nb2 find(D(j,:) eps); if numel(nb2) MinPts neighbors [neighbors, nb2]; %#okAGROW end end if idx(j) 0 || idx(j) -1 idx(j) cid; end k k 1; end end逻辑说明pdist2生成 n×n 距离矩阵样本超过一万时内存会吃紧可以改用knnsearch只查 eps 邻域。expandCluster用队列式扩展把密度可达的点不断加入当前簇。参数上eps决定邻域半径MinPts决定核心点门槛。调用时[idx, noise] dbscan_simple(X, 0.5, 5);其中 eps 需要根据数据尺度调标准化后通常从 0.3 到 0.8 试。3.4 密度峰值聚类的决策图与中心选取DPC 的关键是算 ρ 和 δ然后画决策图。下面用高斯核计算局部密度避免截断核在 dc 附近的不连续。function [rho, delta, ord] dpc_core(X, dc) n size(X,1); D pdist2(X, X); % 高斯核局部密度 rho sum(exp(-(D./dc).^2), 2) - 1; % 减去自身 [~, ord] sort(rho, descend); delta zeros(n,1); nneigh zeros(n,1); delta(ord(1)) max(D(ord(1),:)); % 密度最高点的 delta 取最大距离 for i 2:n idx ord(i); higher ord(1:i-1); [delta(idx), pos] min(D(idx, higher)); nneigh(idx) higher(pos); end end逻辑说明dc是截断距离控制密度计算的邻域范围。rho用高斯核累加离得越近贡献越大。delta对每个点找密度比它高的最近点距离密度最高点的 delta 取全局最大距离。nneigh记录每个点的父节点后续分配簇时沿父节点回溯。参数上dc建议取所有成对距离的第 1% 到 2% 分位数可以用dc prctile(pdist(X), 2);快速估计。决策图上ρ 和 δ 都大的点就是中心通常选 δ 明显高于其余点的几个。3.5 三种结果的可视化对比与指标计算跑完三种算法后把标签画在同一张图里对比再算调整兰德指数ARI看一致性。% 假设 idx_kmeans, idx_dbscan, idx_dpc 已得到 figure; subplot(1,3,1); gscatter(X(:,1), X(:,2), idx_kmeans); title(k-means); subplot(1,3,2); gscatter(X(:,1), X(:,2), idx_dbscan); title(DBSCAN); subplot(1,3,3); gscatter(X(:,1), X(:,2), idx_dpc); title(DPC); % 调整兰德指数需要真实标签或两两对比 ari_kd rand_index(idx_kmeans, idx_dbscan); fprintf(k-means 与 DBSCAN 的 ARI: %.3f\n, ari_kd);逻辑说明gscatter按标签着色能直观看出簇形状差异。ARI 衡量两种划分的一致性取值 0 到 1越接近 1 越一致。注意 DBSCAN 的噪声标签是 -1算 ARI 前要统一处理可以只对非噪声点计算。参数上如果三种算法结果差异很大不要急着怀疑代码先回到数据分布看哪种假设更匹配。4. 参数调优与结果验证的实操细节4.1 k-means 的 K 值选择肘部法之外还要看业务约束肘部法在 SSE 曲线上找拐点但拐点有时不明显。更稳的做法是结合轮廓系数和业务可解释性。比如用户分群K3 可能对应高、中、低价值K5 可能细分出流失预警群这时候业务含义比数学指标更重要。MATLAB 的evalclusters函数可以一次性比较多种准则。eva evalclusters(X, kmeans, CalinskiHarabasz, KList, 2:8); fprintf(Calinski-Harabasz 推荐 K: %d\n, eva.OptimalK);逻辑说明CalinskiHarabasz准则计算簇间离散度与簇内离散度的比值越大越好。KList指定搜索范围。这个函数内部会多次运行 k-means样本大时耗时明显建议先抽样。4.2 DBSCAN 的 eps 选择k-距离曲线找拐点k-距离曲线是 DBSCAN 调参最实用的工具。对每个点计算到第 k 近邻的距离升序排列拐点处的距离就是 eps 的合理值。k 5; % 与 MinPts 一致 [knn_dist, ~] knnsearch(X, X, K, k1); k_dist sort(knn_dist(:, end), ascend); figure; plot(k_dist); xlabel(样本按距离排序); ylabel(第 k 近邻距离); title(k-距离曲线拐点处即 eps 候选);逻辑说明knnsearch的K设为 k1 是因为第一个最近邻是自身。取最后一列即第 k 近邻距离。曲线急剧上升的拐点对应 eps拐点之前是密集区之后是稀疏区。参数上k 通常取 MinPts二维数据从 4 开始试。如果曲线没有明显拐点说明数据密度连续变化DBSCAN 可能不是最佳选择。4.3 DPC 的 dc 选取分位数法与平均邻居数dc 太大所有点密度趋同决策图区分度下降dc 太小密度估计噪声大。常用做法是让平均邻居数占总数的 1% 到 2%。D pdist(X); dc prctile(D, 2); % 取距离分布的 2% 分位数 fprintf(dc %.4f平均邻居数约 %.1f\n, dc, sum(D dc) / size(X,1));逻辑说明prctile取分位数2% 是经验起点。sum(D dc) / n估算平均每个点的邻居比例理想在 1% 到 2% 之间。如果比例过高减小分位数过低则增大。DPC 对 dc 的敏感度低于 DBSCAN 对 eps 的敏感度但仍建议在决策图上验证中心点是否稳定。4.4 用调整兰德指数和轮廓系数做交叉验证没有真实标签时用内部指标有真实标签时用 ARI 或归一化互信息。MATLAB 的evalclusters支持多种内部准则但只针对部分算法。自己实现 ARI 也不复杂。function ari rand_index(labels1, labels2) % 只保留非噪声点 valid labels1 0 labels2 0; l1 labels1(valid); l2 labels2(valid); n numel(l1); % 列联表 C accumarray([l1, l2], 1); sum_comb (x) sum(x .* (x - 1) / 2); a sum_comb(sum(C, 2)); b sum_comb(sum(C, 1)); c sum_comb(C(:)); expected a * b / (n * (n - 1) / 2); max_idx (a b) / 2; ari (c - expected) / (max_idx - expected); end逻辑说明accumarray构建列联表统计两种划分中同簇和异簇的配对情况。ARI 取值 1 表示完全一致0 表示随机水平。参数上valid过滤掉 DBSCAN 的噪声点避免 -1 标签干扰。如果三种算法 ARI 都在 0.5 以上说明数据本身簇结构清晰如果低于 0.2说明算法假设与数据不匹配需要换方法或重新审视特征。5. 避坑与排查三种聚类算法最容易翻车的五个地方5.1 k-means 结果每次不一样现象同一份数据跑两次kmeans标签顺序和簇边界都变了。原因默认Replicates1初始化随机且簇标签编号本身不固定。解决设Replicates为 5 到 10并用idx的簇内均值排序后重新编号保证可复现。另外rng固定种子只影响数据生成不影响kmeans内部初始化要复现需在调用前再设一次rng。5.2 DBSCAN 把大部分点判成噪声现象idx里 -1 占比超过一半。原因eps 太小或 MinPts 太大核心点数量不足。解决先画 k-距离曲线确认 eps 拐点再把 MinPts 降到维度加一。如果数据本身密度差异极大考虑对数据做标准化或对数变换让密度更均匀。标准化后 eps 通常在 0.3 到 0.8 之间超过 1 说明数据尺度没处理好。5.3 DPC 决策图右上角点太多选不出中心现象ρ-δ 图上多个点 δ 值接近无法判断谁是中心。原因dc 过大导致密度趋同或者数据本身没有明显密度峰值。解决减小 dc 分位数到 1%重新计算如果仍然分不开说明数据不适合 DPC改用 DBSCAN 或谱聚类。另一个检查点是看 γ ρ × δ 的排序取前几个作为中心候选比肉眼判断更稳。5.4 距离矩阵内存溢出现象样本数超过一万时pdist2报内存不足。原因n×n 距离矩阵是 O(n²) 内存两万样本就是 3.2 GBdouble 类型。解决DBSCAN 改用knnsearch只查 eps 邻域不存全矩阵DPC 用分块计算或随机采样估计 dc再对全量数据用knnsearch找近邻。MATLAB 的knnsearch支持 KD-tree低维数据效率很高。5.5 标准化之后结果反而变差现象对数据做零均值单位方差标准化后聚类结果不如原始数据。原因标准化会放大低方差特征的权重如果某个特征本身噪声大标准化后噪声被放大。解决先做特征筛选去掉方差极低或与业务无关的列再标准化。或者改用鲁棒标准化中位数和四分位距MATLAB 里可以手动实现X_robust (X - median(X)) ./ iqr(X);。6. 把三种算法封装成可复用函数与批量对比脚本走到这一步你已经能单独跑通每种算法。但实际项目里数据会换、参数会调、对比要反复做。我一般会把三种算法封装成统一接口输入数据矩阵和参数结构体输出标签和指标然后用一个批量脚本遍历参数组合。这样下次换数据只需改一行路径。function results compare_clustering(X, params) % 统一接口输入数据 X 和参数结构体 params % 输出 results 包含三种算法的标签和指标 results struct(); % k-means [idx_km, ~, sumd] kmeans(X, params.K, Replicates, 5); results.kmeans.idx idx_km; results.kmeans.sse sum(sumd); results.kmeans.sil mean(silhouette(X, idx_km)); % DBSCAN [idx_db, ~] dbscan_simple(X, params.eps, params.MinPts); results.dbscan.idx idx_db; results.dbscan.noise_ratio mean(idx_db -1); % DPC [rho, delta, ord] dpc_core(X, params.dc); gamma rho .* delta; [~, center_idx] sort(gamma, descend); centers center_idx(1:params.n_clusters); idx_dpc assign_dpc(X, rho, delta, ord, centers); results.dpc.idx idx_dpc; results.dpc.centers centers; % 交叉指标 results.ari_km_db rand_index(idx_km, idx_db); results.ari_km_dpc rand_index(idx_km, idx_dpc); end逻辑说明params结构体把 K、eps、MinPts、dc、n_clusters 集中管理换数据集时只改这个结构体。assign_dpc是分配函数沿ord顺序把每个点归到最近的更高密度点所属簇中心点自己成簇。gamma排序取前 n_clusters 个作为中心比手动看决策图更可复现。参数上n_clusters对 DPC 是必须的因为算法本身不自动确定簇数需要从决策图或 γ 排序中选。批量对比时我会写一个循环遍历 eps 和 dc 的候选值记录每种组合的 ARI 和噪声比例最后画热力图找稳定区域。eps_list 0.3:0.1:0.8; dc_list prctile(pdist(X), [1, 2, 3]); ari_matrix zeros(length(eps_list), length(dc_list)); for i 1:length(eps_list) for j 1:length(dc_list) params.eps eps_list(i); params.dc dc_list(j); params.K 3; params.MinPts 5; params.n_clusters 3; res compare_clustering(X, params); ari_matrix(i,j) res.ari_km_db; end end figure; imagesc(ari_matrix); colorbar; xlabel(dc 分位数); ylabel(eps); title(k-means 与 DBSCAN 的 ARI 随参数变化);逻辑说明imagesc画热力图颜色越亮表示两种算法一致性越高。稳定区域通常出现在 eps 和 dc 的中间值附近边缘区域一致性骤降说明参数到了临界点。这个脚本跑一遍就能看出哪些参数组合值得进一步看避免盲目试。参数上eps_list和dc_list的范围根据数据尺度调整标准化后从 0.3 到 0.8 是安全起点。最后说一个我踩过的坑不要迷信单一指标。ARI 高不代表聚类结果对业务有用轮廓系数好也不代表簇中心有物理意义。我现在的习惯是跑完指标后一定把标签映射回原始特征空间看每个簇的均值向量和样本量再决定要不要调参数。这个习惯帮我省了很多次「指标漂亮但结论荒谬」的后悔药。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
