从零到一用Matlab亲手实现Kmeans聚类不只是调包如果你刚开始接触数据分析面对一堆散乱的数据点感到无从下手或者你听说过“聚类”这个概念但总觉得它被包裹在一层神秘的技术面纱之后那么这篇文章就是为你准备的。我们不会停留在“调用一个函数”的层面而是要亲手从零开始用Matlab构建一个完整的Kmeans聚类算法。这个过程远比直接使用kmeans内置函数更有价值——它能让你真正理解数据是如何被“看见”和“组织”的而不仅仅是得到一个结果。无论你是理工科的学生还是希望将数据分析技能应用于自己领域的初学者跟随本文的步骤你将在理解原理的同时获得一份可以随时修改、调试和应用于自己项目的实战代码。1. 聚类与Kmeans用算法为数据绘制地图当我们面对一个陌生的数据集时第一步往往是“观察”。但如果数据维度很高或者样本量巨大人眼就无能为力了。聚类算法的核心思想就是让机器自动发现数据中内在的分组结构将相似的数据点归为同一类不相似的点区分开。想象一下你有一堆未分类的新闻文章聚类可以帮助你自动分出体育、科技、财经等板块或者在电商用户行为数据中区分出“价格敏感型”、“品质追求型”等不同客群。这是一种无监督学习我们不需要事先告诉机器哪些数据属于哪一类完全由数据自己“说话”。在众多聚类算法中Kmeans以其思想直观、实现相对简单、效率较高的特点成为了最经典和广泛使用的入门算法。它的目标很明确将n个数据点划分到k个簇中使得每个数据点到其所属簇的中心点质心的距离平方和最小。这个“距离”通常采用欧氏距离也就是我们最直观的“直线距离”。Kmeans这个名字本身就揭示了其核心流程K代表最终要形成的簇的个数means代表均值即用簇内所有点的均值来代表这个簇。整个算法就像一个不断迭代优化的过程初始化随机选择K个点作为初始的簇中心质心。分配计算每个数据点到各个质心的距离并将其分配给距离最近的质心所在的簇。更新重新计算每个簇的质心即该簇所有点的均值。迭代重复步骤2和3直到质心的位置不再发生显著变化或者达到预设的迭代次数。这个过程听起来简单但其中蕴含着几个关键问题也是我们动手实现时需要仔细考虑的K值怎么选初始质心怎么选更好算法一定会收敛吗我们将在后续章节逐一拆解。2. 搭建你的Matlab实验环境与理解数据在开始写代码之前确保你的Matlab环境已经就绪。本文的代码适用于Matlab R2016a及之后的版本大部分基础函数都是通用的。我们首先来准备和观察我们的“实验材料”——数据。2.1 数据准备构造一个可视化的二维数据集为了最直观地展示聚类过程我们使用一个二维数据集。这样每一个数据点都可以在平面上用一个(x, y)坐标表示聚类结果可以用不同颜色的散点图清晰呈现。下面我们手动构造一组数据它混合了四个相对集中的点群。% 清空环境关闭所有图形窗口确保一个干净的起点 clc; clear; close all; % 构造二维数据集这里我们手动输入了两组数据模拟四个簇 % 第一列是x坐标第二列是y坐标 data [ 90, 33; 35, 71; 52, 62; 83, 34; 64, 49; 24, 48; 49, 46; 92, 69; 99, 56; 45, 59; 19, 28; 38, 14; 1, 55; 71, 41; 56, 39; 97, 78; 63, 23; 32, 99; 3, 68; 34, 30; 33, 87; 55, 85; 75, 43; 84, 88; 53, 2; 15, 47; 88, 50; 66, 77; 41, 22; 51, 76; 39, 94; 78, 11; 67, 80; 65, 51; 25, 6; 40, 7; 77, 72; 13, 36; 69, 90; 29, 96; 14, 44; 54, 61; 87, 70; 47, 60; 44, 75; 58, 74; 8, 63; 68, 40; 81, 81; 31, 4 ]; % 可视化原始数据 figure(Position, [100, 100, 800, 400]) % 设置图形窗口位置和大小 subplot(1, 2, 1) % 创建1行2列的子图当前操作第1个 scatter(data(:, 1), data(:, 2), 40, b, filled, MarkerEdgeColor, k, LineWidth, 0.5); grid on; % 显示网格便于观察 xlabel(特征一); ylabel(特征二); title(原始数据散点图, FontSize, 12, FontWeight, bold); axis equal; % 设置坐标轴比例相等防止图形变形运行这段代码你会在左侧子图中看到50个蓝色的散点。虽然数据是人工构造的但肉眼已经能大致分辨出有几团点聚集在一起。这就是我们接下来要让算法去自动识别的结构。提示在实际项目中你的数据可能来自CSV文件、数据库或传感器。可以使用readtable,xlsread或load函数加载。数据预处理如标准化、归一化对于Kmeans至关重要因为该算法对特征的量纲敏感。我们这里为了简化使用了构造的、尺度相近的数据。2.2 核心挑战如何确定K值在我们动手实现算法主循环之前必须先面对Kmeans的第一个也是最重要的超参数K值即簇的个数。在上面的数据中我们“知道”构造了4个簇所以设K4。但在真实未知数据中K是未知的。有几种常用的方法可以帮助我们选择K肘部法则计算不同K值下所有数据点到其所属簇质心的距离平方和称为畸变程度。随着K增大这个值会下降。当K增加到真实簇数时再增加K带来的收益畸变程度下降幅度会骤减图像上会出现一个“肘部”拐点。轮廓系数它结合了内聚度一个点与同簇其他点的平均距离和分离度一个点与最近其他簇中所有点的平均距离。轮廓系数越接近1说明聚类效果越好。我们可以计算不同K值下的平均轮廓系数选择使其最大化的K。为了不给初学者增加过多负担本文后续将预设K4进行算法实现。但请你务必记住在实际应用中盲目指定K值是非常危险的很可能得到没有意义的结果。我们会在第4章探讨如何将K值选择自动化。3. 亲手编写Kmeans核心函数从原理到代码现在进入最核心的部分我们不使用Matlab自带的kmeans函数而是自己从头实现一个my_kmeans。这个过程会让你对算法的每一个细节都了然于胸。3.1 函数骨架与初始化策略首先我们定义函数的输入输出。输入是数据矩阵data每行一个样本每列一个特征和簇的数目K。输出是每个样本所属簇的索引idx以及最终的K个质心centroids。function [idx, centroids] my_kmeans(data, K, max_iters, tolerance) % MY_KMEANS 自定义K均值聚类实现 % 输入: % data - M x N 矩阵M个样本N个特征 % K - 期望的簇数量 % max_iters - 最大迭代次数可选默认100 % tolerance - 质心变化容差可选默认1e-4 % 输出: % idx - M x 1 向量每个样本的簇标签1到K % centroids - K x N 矩阵最终各个簇的质心坐标 % 设置默认参数 if nargin 3 max_iters 100; end if nargin 4 tolerance 1e-4; end [M, N] size(data); % M样本数N特征维数 % --- 初始化质心采用“K-means”改进策略 --- centroids zeros(K, N); % 1. 随机选择第一个质心 rng(default); % 固定随机种子使结果可复现 first_idx randi(M); centroids(1, :) data(first_idx, :); % 2. 选择后续质心距离已有质心越远的点被选中的概率越大 for i 2:K distances zeros(M, i-1); for j 1:(i-1) % 计算每个点到已有质心的欧氏距离 diff data - centroids(j, :); distances(:, j) sqrt(sum(diff.^2, 2)); end % 取每个点到其最近质心的距离 min_distances min(distances, [], 2); % 将距离平方作为概率分布 prob_distribution min_distances.^2 / sum(min_distances.^2); % 根据概率分布随机选择一个点作为新质心 cum_prob cumsum(prob_distribution); new_idx find(cum_prob rand(), 1, first); centroids(i, :) data(new_idx, :); end % --- 初始化完成 ---这里没有采用最简单的完全随机选择而是实现了**K-means**初始化。这种方法能显著降低算法对初始质心的敏感性提高收敛速度和最终结果的质量。其核心思想是让初始质心彼此尽可能远离。3.2 主循环分配与更新的艺术初始化完成后算法进入迭代循环将每个点分配给最近的质心然后根据新的分配重新计算质心。% 主迭代循环 for iter 1:max_iters % --- 分配步骤计算每个点到所有质心的距离 --- distances_to_centroids zeros(M, K); for k 1:K diff data - centroids(k, :); distances_to_centroids(:, k) sqrt(sum(diff.^2, 2)); % 欧氏距离 end % 找到每个样本距离最近的质心索引 [~, idx] min(distances_to_centroids, [], 2); % idx 是 Mx1 向量 % --- 更新步骤计算新的质心 --- new_centroids zeros(K, N); for k 1:K % 找出所有属于第k簇的样本点 members (idx k); if sum(members) 0 % 防止空簇 new_centroids(k, :) mean(data(members, :), 1); else % 如果出现空簇重新随机初始化该质心一种处理策略 new_centroids(k, :) data(randi(M), :); warning(迭代 %d: 簇 %d 为空已重新初始化其质心。, iter, k); end end % --- 收敛判断质心移动是否足够小 --- centroid_shift sqrt(sum((new_centroids - centroids).^2, 2)); % 每个质心的移动距离 max_shift max(centroid_shift); if max_shift tolerance fprintf(算法在 %d 次迭代后收敛。\n, iter); break; end % 更新质心准备下一次迭代 centroids new_centroids; end if iter max_iters fprintf(达到最大迭代次数 %d算法可能未完全收敛。\n, max_iters); end end这个循环体现了Kmeans的本质分配通过计算欧氏距离并取最小值完成了“物以类聚”。更新计算簇内均值得到了新的“聚类的中心”。收敛判断当质心几乎不动时认为算法已经找到了一个稳定的划分。注意代码中包含了空簇处理。在随机初始化或迭代过程中有可能某个簇会失去所有样本点。这里我们采用了一种简单的策略重新随机选择一个数据点作为该簇的新质心。更复杂的策略可以是将距离其质心最远的点分配过来。3.3 可视化与结果评估让聚类结果“看得见”算法跑完了我们怎么知道它好不好对于二维数据可视化是最直接的检验方式。同时我们引入一个定量的评估指标轮廓系数。% 调用自定义函数进行聚类 K 4; [idx, final_centroids] my_kmeans(data, K); % 1. 可视化聚类结果 subplot(1, 2, 2) % 操作第2个子图 colors lines(K); % 获取K种区分度高的颜色 for k 1:K cluster_points data(idx k, :); scatter(cluster_points(:, 1), cluster_points(:, 2), 60, ... filled, MarkerFaceColor, colors(k, :), ... MarkerEdgeColor, k, LineWidth, 0.8); hold on; end % 绘制最终质心 plot(final_centroids(:, 1), final_centroids(:, 2), kh, ... MarkerSize, 12, MarkerFaceColor, y, LineWidth, 2); grid on; xlabel(特征一); ylabel(特征二); title(sprintf(自定义Kmeans聚类结果 (K%d), K), FontSize, 12, FontWeight, bold); legend(簇 1, 簇 2, 簇 3, 簇 4, 质心, Location, best); axis equal; hold off; % 2. 计算并显示轮廓系数 silhouette_values silhouette(data, idx); avg_silhouette mean(silhouette_values); fprintf(本次聚类的平均轮廓系数为%.4f\n, avg_silhouette); % 3. 绘制轮廓系数图 figure; silhouette(data, idx); title(sprintf(各样本轮廓系数分布 (平均: %.3f), avg_silhouette), FontSize, 11);运行后你会在第二个子图中看到数据点被清晰地分成了四种颜色黄色的菱形代表最终收敛的质心。轮廓系数图则提供了更细致的评估每个样本的轮廓系数被绘制成水平条形理想情况下大部分条形都较长且为正值越接近1越好。一个较高的平均轮廓系数例如大于0.5通常意味着聚类结构较为清晰。4. 超越基础解决实战中的关键问题与优化一个能运行的Kmeans只是起点。要让它在真实数据上发挥作用我们必须考虑更多。4.1 自动化选择最佳K值肘部法则实现如前所述K值不能靠猜。下面我们实现肘部法则让数据自己“告诉”我们可能的K值范围。% 尝试不同的K值计算其畸变程度Within-Cluster-Sum of Squared Errors, WCSS max_K_to_try 8; wcss zeros(max_K_to_try, 1); % 存储每个K对应的WCSS for k 1:max_K_to_try [idx_temp, centroids_temp] my_kmeans(data, k); % 计算WCSS每个点到其所属质心距离的平方和 sum_sq_dist 0; for i 1:k cluster_points data(idx_temp i, :); if ~isempty(cluster_points) dist sqrt(sum((cluster_points - centroids_temp(i, :)).^2, 2)); sum_sq_dist sum_sq_dist sum(dist.^2); end end wcss(k) sum_sq_dist; end % 绘制肘部法则图 figure; plot(1:max_K_to_try, wcss, bo-, LineWidth, 1.5, MarkerSize, 8, MarkerFaceColor, b); xlabel(簇的数量 K); ylabel(WCSS (畸变程度)); title(肘部法则寻找最佳K值); grid on; hold on; % 可以尝试标注可能的“肘部”例如K4 plot(4, wcss(4), r*, MarkerSize, 15, LineWidth, 2); text(4, wcss(4), 可能的肘部 (K4), VerticalAlignment, bottom); hold off;观察生成的折线图WCSS会随着K增大而单调递减。当K小于真实簇数时每增加一个簇WCSS会大幅下降当K超过真实簇数后下降幅度会变得平缓。那个拐点像胳膊肘一样对应的K值就是肘部法则建议的值。从我们的构造数据图中你应该能清晰地看到在K4处有一个明显的转折。4.2 处理高维与量纲问题数据标准化我们的示例数据是二维且量纲相近的。现实中的数据比如一个人的“年龄20-60”和“收入5000-100000”量级差异巨大。Kmeans基于欧氏距离量级大的特征会主导距离计算导致聚类结果失真。标准化是必须的预处理步骤。常用的标准化方法有Z-score标准化将数据转换为均值为0标准差为1的分布。(x - mean(x)) / std(x)Min-Max归一化将数据缩放到[0, 1]区间。(x - min(x)) / (max(x) - min(x))下面演示如何在聚类流程中加入Z-score标准化% 假设data_raw是我们的原始数据 data_raw your_high_dim_data; % 替换为你的实际数据 % Z-score标准化 data_mean mean(data_raw, 1); data_std std(data_raw, 0, 1); % 0表示使用N-1计算标准差样本标准差 % 防止标准差为0导致除零错误 data_std(data_std 0) 1; data_normalized (data_raw - data_mean) ./ data_std; % 对标准化后的数据进行聚类 K 4; [idx_norm, centroids_norm] my_kmeans(data_normalized, K); % 注意得到的质心centroids_norm是在标准化空间中的。 % 如果需要原始空间的质心可以进行逆变换 centroids_original_space centroids_norm .* data_std data_mean;4.3 与内置函数对比验证与理解差异我们实现了自己的版本那么和Matlab内置的kmeans函数相比结果和性能如何呢进行对比是一个很好的学习方式。% 使用内置kmeans函数需要Statistics and Machine Learning Toolbox if license(test, statistics_toolbox) rng default; % 固定随机种子确保比较公平 [idx_builtin, centroids_builtin] kmeans(data, K, Display, iter, MaxIter, 100); % 比较聚类结果的一致性使用调整兰德指数ARI % 注意ARI是衡量两个聚类结果一致性的指标值域[-1,1]1表示完全一致 % 这里需要安装或实现ARI函数为简化我们直接比较质心位置和轮廓系数 fprintf(\n--- 对比结果 ---\n); fprintf(自定义函数平均轮廓系数: %.4f\n, avg_silhouette); silhouette_builtin silhouette(data, idx_builtin); fprintf(内置函数平均轮廓系数: %.4f\n, mean(silhouette_builtin)); fprintf(\n自定义函数最终质心:\n); disp(final_centroids); fprintf(内置函数最终质心:\n); disp(centroids_builtin); % 可视化内置函数的结果 figure; for k 1:K cluster_pts data(idx_builtin k, :); scatter(cluster_pts(:, 1), cluster_pts(:, 2), 60, ... filled, MarkerFaceColor, colors(k, :), ... MarkerEdgeColor, k, LineWidth, 0.8); hold on; end plot(centroids_builtin(:, 1), centroids_builtin(:, 2), kh, ... MarkerSize, 12, MarkerFaceColor, y, LineWidth, 2); grid on; xlabel(特征一); ylabel(特征二); title(Matlab内置kmeans函数聚类结果, FontSize, 12, FontWeight, bold); legend(簇 1, 簇 2, 簇 3, 簇 4, 质心, Location, best); axis equal; else fprintf(未检测到Statistics and Machine Learning Toolbox无法运行内置kmeans函数进行对比。\n); end通过对比你可能会发现两种实现的结果略有不同这很正常因为Kmeans对初始质心敏感且内置函数可能采用了更复杂的默认设置如更多的重复启动以寻找全局最优。关键在于我们的自定义函数揭示了算法的核心逻辑而内置函数提供了经过高度优化的工业级实现。5. 从示例到实战将你的代码应用于新数据集掌握了核心代码和优化技巧后最后一步就是将其应用于你自己的数据。这里提供一个清晰的流程模板。第一步准备数据将你的数据加载到Matlab工作区并整理成一个数值矩阵my_data其中每一行是一个样本每一列是一个特征。检查是否有缺失值NaN需要进行处理如删除或填充。第二步探索与预处理绘制散点图对于二维或三维数据或查看基本统计信息。决定是否需要进行标准化/归一化强烈建议进行。使用肘部法则或轮廓系数法在你的数据上尝试不同的K值选择一个合理的K。第三步运行聚类% 假设你的数据已经准备好 K 5; % 根据第二步确定 max_iters 150; tolerance 1e-5; [idx_my_project, centroids_my_project] my_kmeans(my_data, K, max_iters, tolerance);第四步分析与解释可视化如果维度允许。计算轮廓系数等指标评估聚类质量。分析每个簇的质心坐标理解每个簇代表的特征。例如在客户分群中一个质心可能代表“高价值、低活跃度”客户另一个代表“低价值、高活跃度”客户。将聚类标签idx_my_project作为新的一列合并回原始数据进行进一步的群体差异分析。一个常见问题算法不收敛或结果不稳定这通常是由于初始质心选择不佳或数据本身结构模糊导致的。除了使用K-means初始化一个实用的技巧是多次运行算法选择WCSS最小或轮廓系数最大的一次结果作为最终输出。num_runs 10; best_wcss inf; best_idx []; best_centroids []; for run 1:num_runs [idx_temp, centroids_temp] my_kmeans(my_data, K); % 计算本次运行的WCSS wcss_temp 0; for i 1:K members my_data(idx_temp i, :); if ~isempty(members) wcss_temp wcss_temp sum(sum((members - centroids_temp(i, :)).^2, 2)); end end % 保留最优结果 if wcss_temp best_wcss best_wcss wcss_temp; best_idx idx_temp; best_centroids centroids_temp; end end fprintf(经过%d次运行最优WCSS为%.2f\n, num_runs, best_wcss);至此你已经拥有了一个从原理理解、代码实现、问题解决到实战应用的全套Kmeans聚类工具。真正的掌握源于动手尝试和解决实际数据中的问题。不妨现在就找一个你感兴趣的数据集用这份代码去探索其中隐藏的群体结构吧。