1. 项目概述：从数据分堆到模式发现
如果你手头有一堆数据，比如几百个客户的消费记录，或者一堆图片的像素值，想看看它们能不能自然地分成几个“小团体”，这时候Kmeans算法就是你工具箱里最趁手的那把“瑞士军刀”。它不是什么高深莫测的黑科技，核心思想简单到可以用一个生活场景来比喻：假设你有一大袋混合的彩色玻璃珠（红、绿、蓝），你的任务是不用看颜色标签，仅凭珠子本身的颜色，把它们分成三堆。你会怎么做？一个很自然的做法是，先随便指定三个珠子作为“中心点”（比如一个偏红的、一个偏绿的、一个偏蓝的），然后把袋子里剩下的每个珠子，都跟这三个中心点比一比颜色，看它跟谁最像（距离最近），就把它归到那一堆里去。等所有珠子都分完堆，你再重新计算每一堆珠子的平均颜色，把这个新的平均颜色作为这一堆新的“中心点”。接着，再用新的中心点去重新给所有珠子分堆……如此反复，直到中心点不再明显移动，分堆结果稳定下来。Kmeans干的就是这个“自动分堆”的活儿，它在数据挖掘、图像分割、客户细分、异常检测等场景中应用极广。今天，我就结合自己多次在数学建模竞赛和实际项目中应用Kmeans的经验，不仅带你用MATLAB手把手实现这个算法，更要把算法里那些容易踩坑的细节和调参心得掰开揉碎了讲清楚。
2. Kmeans算法核心原理拆解：不只是“算距离”
很多人对Kmeans的理解停留在“计算距离、找中心、再分配”的循环上，这没错，但要想用好它，必须深入理解其数学本质和几个关键假设。
2.1 算法目标：最小化“类内”差异
Kmeans算法的终极目标，是让同一个簇（类）内的数据点尽可能相似，不同簇之间的数据点尽可能不同。用数学语言说，就是要最小化所有数据点到其所属簇中心的距离平方和。这个指标叫做误差平方和（Sum of Squared Errors, SSE），也叫“畸变（Distortion）”。其公式为：
SSE = Σ(i=1 to k) Σ(x in Ci) ||x - μi||²
其中，k是簇的个数，Ci是第i个簇的集合，μi是第i个簇的中心点（均值向量），||x - μi||是数据点x到中心μi的距离（通常是欧氏距离）。
为什么是距离的平方？使用平方项有两个主要原因：一是数学上便于求导和优化（在推导中心点更新公式时，对SSE求关于中心点的导数，平方项会让导数形式更简洁，直接得到均值）；二是它对远离中心点的“异常值”给予更大的惩罚，使得算法对簇的形状更敏感，倾向于形成紧凑的球形簇。这也是Kmeans一个重要的局限性：它隐含地假设每个簇都是凸形的、各向同性的（在各个方向上方差相近）。如果你的数据是拉长的带状或嵌套的环形，标准Kmeans可能就力不从心了。
2.2 标准流程与两大核心步骤
Kmeans是一个迭代优化算法，其标准流程可以精炼为两个交替进行的步骤：
分配步骤（Assignment Step）：固定簇中心
{μ1, μ2, ..., μk}，将每个数据点x分配到离它最近的簇中心所属的簇中。即，为每个数据点i寻找簇标签c(i)：c(i) = argmin(j) ||x(i) - μj||²这一步决定了当前的簇划分。更新步骤（Update Step）：固定簇划分
{C1, C2, ..., Ck}，重新计算每个簇的中心点，即该簇内所有数据点的均值向量：μj = (1/|Cj|) * Σ(i in Cj) x(i)这一步根据新的成员关系，优化簇中心的位置。
这两个步骤循环进行，直到满足停止条件。停止条件通常有两种：一是簇中心的位置在连续两次迭代中变化小于一个很小的阈值（如1e-5）；二是分配的簇成员不再发生变化；三是达到预设的最大迭代次数（防止不收敛或震荡）。
注意：这个算法保证每次迭代后SSE都不会增加（通常会减少），最终会收敛到一个局部最优解。但它不能保证找到全局最优解，最终结果严重依赖于初始簇中心的选取。这是实操中需要应对的第一个关键点。
2.3 距离度量的选择：欧氏距离是唯一答案吗？
在MATLAB的kmeans函数中，默认使用平方欧氏距离。但在不同场景下，距离度量需要慎重选择。
- 欧氏距离（‘sqeuclidean’）：最常用，适用于连续数值型特征，且各特征量纲和重要性相当时。它对坐标系旋转是不变的，但受量纲影响大，因此数据标准化（如Z-score标准化）通常是预处理必备步骤。
- 城市街区距离（‘cityblock’）：即曼哈顿距离，对异常值不如欧氏距离敏感。在某些维度特征有明确物理意义，且差异是线性叠加的场景下可能更合适。
- 余弦距离（‘cosine’）：衡量的是向量方向的差异，而非绝对距离。在文本挖掘（如文档词向量）、高维稀疏数据（如用户-物品评分矩阵）中非常有效，因为它只关注模式（哪些维度有值），而不关心数值的绝对大小。
- 相关距离（‘correlation’）：基于皮尔逊相关系数，衡量两个向量变化趋势的相似性。在时间序列分析、基因表达数据分析中常用。
在MATLAB中，可以通过kmeans(data, k, ‘Distance‘, ‘cosine‘)这样的参数来指定。选择距离度量前，一定要思考你的数据特征和业务目标。例如，做客户画像，如果特征包括“年收入（万元）”和“每周购物频率（次）”，直接算欧氏距离，年收入的微小波动（如1万元）会完全压倒购物频率的显著变化（如5次）。必须先标准化。
3. MATLAB实战编程：从函数调用到自己动手
MATLAB提供了高度优化的kmeans函数，但为了真正理解算法，我们分两步走：先学习如何正确使用官方函数并解读结果，再自己动手实现一个简化版，加深对每个环节的理解。
3.1 使用内置kmeans函数：参数详解与结果分析
我们先在一个经典的鸢尾花（Iris）数据集上演示。这个数据集包含150个样本，每个样本有4个特征（花萼长宽、花瓣长宽），真实类别有3种。我们用Kmeans尝试对其进行无监督聚类。
% 1. 加载并准备数据（这里用内置数据集为例，实际中可能是你自己的数据矩阵） load fisheriris; % 加载鸢尾花数据集，变量名为meas（150x4数据）和species（150x1标签） data = meas; % 我们的特征数据 % 2. 数据标准化（强烈建议！） data_zscore = zscore(data); % Z-score标准化，使每个特征均值为0，标准差为1 % 3. 确定簇数量k（这里我们先假设知道是3，实际中需要方法确定，后文会讲） k = 3; % 4. 调用kmeans函数，进行多次重复以降低初始值影响 opts = statset(‘Display‘, ‘final‘); % 设置显示最终结果 [idx, C, sumd, D] = kmeans(data_zscore, k, ‘Options‘, opts, ‘Replicates‘, 10, ‘MaxIter‘, 1000); % 参数解释： % - data_zscore: 输入数据矩阵，每行一个样本，每列一个特征。 % - k: 簇的个数。 % - ‘Options‘, opts: 算法选项，这里设置显示最终迭代信息。 % - ‘Replicates‘, 10: 重复运行算法10次，每次使用不同的随机初始中心，返回SSE最小的那次结果。这是**避免糟糕局部最优的必备操作**。 % - ‘MaxIter‘, 1000: 最大迭代次数。 % - 输出: % idx: 一个150x1的向量，每个元素是对应样本的簇索引（1,2,3）。 % C: 一个k x 4的矩阵，每一行是一个簇的中心点坐标。 % sumd: 一个1 x k的向量，每个元素是该簇内所有点到中心点的距离之和（不是平方和）。 % D: 一个150 x k的矩阵，D(i,j)是第i个样本到第j个簇中心的距离。 % 5. 可视化结果（以降维后的二维散点图为例） % 使用主成分分析(PCA)将四维数据降至二维以便可视化 [coeff, score] = pca(data_zscore); figure; gscatter(score(:,1), score(:,2), idx); % 根据聚类结果idx着色 hold on; % 将簇中心也投影到PCA空间 C_pca = C * coeff(:,1:2); % 注意：中心点C是在标准化后的原始空间，需要乘以PCA系数得到投影坐标 plot(C_pca(:,1), C_pca(:,2), ‘kx‘, ‘MarkerSize‘, 15, ‘LineWidth‘, 3); % 用黑色‘x‘标出中心 title(‘K-means聚类结果（PCA降维可视化）‘); xlabel(‘第一主成分‘); ylabel(‘第二主成分‘); legend(‘Cluster 1‘, ‘Cluster 2‘, ‘Cluster 3‘, ‘Cluster Centers‘); % 6. 评估聚类效果（与真实标签对比，仅用于有标签数据的验证） % 计算调整兰德指数(Adjusted Rand Index, ARI)或归一化互信息(NMI) % 这里使用简单的混淆矩阵和调整兰德指数（需要下载外部函数或使用Statistics and Machine Learning Toolbox） if exist(‘species‘, ‘var‘) % 将文本标签转为数字 [~, ~, true_labels] = unique(species); % 计算调整兰德指数 (需要 stats toolbox) % ari = adjustedrand(idx, true_labels); % fprintf(‘调整兰德指数(ARI)为：%.4f\n‘, ari); % 简易计算混淆矩阵 confusion_mat = confusionmat(true_labels, idx); disp(‘混淆矩阵（行：真实类别， 列：预测簇）:‘); disp(confusion_mat); end运行这段代码，你会在命令窗口看到最终的迭代信息，并得到一张聚类结果图。从混淆矩阵可以直观看出聚类结果与真实类别的匹配程度。请注意，Kmeans是无监督学习，它发现的“簇”不一定对应真实的“类”，顺序也是任意的。图中簇中心（黑色X）大致位于各簇的“重心”位置。
3.2 手撕Kmeans：实现一个简易版理解迭代过程
自己实现一遍是理解算法细节的最佳途径。我们来实现一个基础版本，包含随机初始化和基本的迭代循环。
function [idx, centers, SSE_history] = my_kmeans(data, k, max_iters, tol) % 简易版Kmeans实现 % 输入： % data - m x n 矩阵，m个样本，n个特征 % k - 簇数量 % max_iters - 最大迭代次数（可选，默认100） % tol - 中心点变化容忍度（可选，默认1e-5） % 输出： % idx - m x 1 向量，每个样本的簇标签（1到k） % centers - k x n 矩阵，最终簇中心 % SSE_history - 记录每次迭代的SSE，用于观察收敛 if nargin < 3 max_iters = 100; end if nargin < 4 tol = 1e-5; end [m, n] = size(data); SSE_history = zeros(max_iters, 1); % 1. 初始化：随机选择k个样本作为初始中心点（这是一种简单策略） rng(‘default‘); % 设置随机种子，使结果可复现 random_indices = randperm(m, k); centers = data(random_indices, :); % k x n old_centers = centers; for iter = 1:max_iters % 2. 分配步骤：计算每个样本到所有中心的距离，并分配标签 distances = zeros(m, k); for i = 1:k % 计算data中所有行与centers(i,:)的欧氏距离（向量化操作，高效） % 利用 (a-b)^2 = a^2 + b^2 - 2ab data_sq = sum(data.^2, 2); % m x 1 center_sq = sum(centers(i,:).^2); cross_term = data * centers(i,:)‘; distances(:, i) = data_sq + center_sq - 2*cross_term; % 平方欧氏距离 end [~, idx] = min(distances, [], 2); % 找到每个样本距离最小的中心索引 % 3. 更新步骤：重新计算每个簇的中心（均值） for i = 1:k members = (idx == i); % 逻辑索引，找出属于簇i的样本 if sum(members) > 0 % 防止空簇 centers(i, :) = mean(data(members, :), 1); else % 如果出现空簇，处理策略：随机选择一个样本作为新中心 warning(‘簇 %d 为空，将随机重新初始化该中心.‘, i); centers(i, :) = data(randi(m), :); end end % 4. 计算当前SSE并记录 current_SSE = 0; for i = 1:k members = data(idx == i, :); if ~isempty(members) diff = members - centers(i, :); % 差值矩阵 current_SSE = current_SSE + sum(sum(diff.^2, 2)); % 累加平方和 end end SSE_history(iter) = current_SSE; % 5. 检查收敛条件：中心点变化是否小于容忍度 center_shift = sqrt(sum((centers - old_centers).^2, 2)); % 每个中心移动的距离 if max(center_shift) < tol fprintf(‘迭代在 %d 步后收敛。\n‘, iter); SSE_history = SSE_history(1:iter); % 截断记录 break; end old_centers = centers; % 更新旧中心 end if iter == max_iters warning(‘达到最大迭代次数 %d，可能未完全收敛。\n‘, max_iters); end end你可以用同样的鸢尾花数据测试这个函数，并与内置函数的结果对比。这个简易实现揭示了几个关键点：
- 初始化的随机性：我们用了最简单的“随机选点”法，效果不稳定。
- 空簇处理：在更新步骤中，如果某个簇失去了所有成员（空簇），我们必须有处理策略（如代码中的随机重初始化），否则算法会出错。
- 距离计算的向量化：我们使用了向量化运算
data * centers(i,:)‘来计算点积，这比用循环逐样本计算快得多。这是MATLAB编程的性能关键。 - 收敛判断：我们监控中心点的最大移动距离和SSE的变化。
实操心得：自己实现的版本有助于理解，但在生产环境或严肃分析中，务必使用MATLAB内置的
kmeans函数。它不仅经过高度优化，速度极快，还内置了更鲁棒的初始化方法（如‘kmeans++‘）和并行计算支持，并且妥善处理了各种边界情况（如空簇）。自己写的版本更适合教学和调试思想。
4. 关键问题与调参实战：让Kmeans真正为你所用
会用函数只是开始，解决实际问题时，以下几个问题才是真正的挑战。
4.1 如何确定最佳的簇数量k？
这是Kmeans应用中最经典、最没有银弹的问题。我们不知道数据应该分成几堆。以下是几种常用方法：
- 肘部法则（Elbow Method）：原理是随着k增大，SSE会下降（因为每个簇更精细）。我们希望找到一个k，使得再增加k，SSE的下降幅度突然变缓，这个拐点像“肘部”。在MATLAB中实现：
% 肘部法则示例 max_k = 10; sse_values = zeros(max_k, 1); for k = 1:max_k [~, ~, sumd] = kmeans(data_zscore, k, ‘Replicates‘, 5, ‘Display‘, ‘off‘); sse_values(k) = sum(sumd); % 注意：内置kmeans返回的sumd是距离和，不是平方和。对于肘部法则，趋势一致即可。 % 更精确的SSE需要自己计算：sse = sum(sumd.^2)? 不，sumd已经是距离和。这里我们用 sum(sumd) 作为SSE的代理指标。 end figure; plot(1:max_k, sse_values, ‘bo-‘); xlabel(‘簇数量 k‘); ylabel(‘SSE (或距离和)‘); title(‘肘部法则‘); grid on;观察曲线，寻找那个“肘点”。但很多时候曲线是平滑的，没有明显的肘部，这就需要结合其他方法。
- 轮廓系数（Silhouette Coefficient）：它结合了簇内的凝聚度和簇间的分离度。对于每个样本i，计算：
a(i)= i到同簇其他样本的平均距离（凝聚度）。b(i)= i到其他所有簇中，样本i到该簇所有样本平均距离的最小值（分离度）。- 样本i的轮廓系数
s(i) = (b(i) - a(i)) / max(a(i), b(i))。s(i)范围在[-1, 1]，越接近1说明聚类越好。所有样本s(i)的均值即为整体轮廓系数。我们可以计算不同k下的平均轮廓系数，取最大值对应的k。
% 轮廓系数示例 (需要 Statistics and Machine Learning Toolbox) silhouette_values = zeros(max_k-1, 1); % k从2开始 for k = 2:max_k idx = kmeans(data_zscore, k, ‘Replicates‘, 5, ‘Display‘, ‘off‘); silhouette_vals = silhouette(data_zscore, idx); silhouette_values(k-1) = mean(silhouette_vals); end figure; plot(2:max_k, silhouette_values, ‘rs-‘); xlabel(‘簇数量 k‘); ylabel(‘平均轮廓系数‘); title(‘轮廓系数法‘); grid on; [best_sil, best_k_idx] = max(silhouette_values); best_k = best_k_idx + 1; fprintf(‘轮廓系数建议的最佳k值为：%d (系数=%.4f)\n‘, best_k, best_sil);- 间隙统计量（Gap Statistic）：比较实际数据的SSE与随机参考数据（如均匀分布）的SSE的差异。当实际数据的对数SSE与随机期望的对数SSE之差（间隙）最大时，对应的k较优。MATLAB没有内置函数，需要自己实现或找第三方代码。
- 业务理解与可视化：有时，数学指标不如业务直觉。将数据用PCA或t-SNE降维到2D/3D可视化，观察不同k下的聚类结果，结合你对数据背景的理解（比如，你知道客户大概分高、中、低三档价值），做出最终决策。
注意事项：没有一种方法绝对可靠。我通常的做法是“三管齐下”：先看肘部图和轮廓系数图，得到一个候选k的范围（比如2-5），然后把这些k值的结果都做出来，通过降维可视化观察簇的分离情况和合理性，最后结合业务目标拍板。在数学建模中，需要清晰阐述你选择k的方法和理由。
4.2 初始化策略：Kmeans++ 为什么更好？
随机初始化可能导致算法收敛到差的局部最优，或者迭代次数增多。Kmeans++是一种智能初始化方法，基本思想是让初始中心点彼此尽可能远离。步骤：
- 随机选择第一个中心点。
- 对于每个数据点，计算它到已选中心点的最短距离
D(x)。 - 按照
D(x)²的概率分布，随机选择下一个中心点（距离越远的点被选中的概率越大）。 - 重复2-3步，直到选出k个中心点。
MATLAB的kmeans函数默认使用的就是‘kmeans++‘初始化（从R2019a开始？早期版本需指定‘Start‘, ‘plus‘）。它能显著提高聚类质量和解的稳定性。在你自己的实现中，强烈建议加入Kmeans++初始化，代码稍复杂，但效果提升明显。
4.3 处理不同尺度与异常值
数据标准化是Kmeans的标配，而非可选。如果特征量纲不同（如年龄和收入），量级大的特征会主导距离计算。常用标准化方法：
- Z-score标准化：
(x - mean)/std。最常用，适用于特征大致符合正态分布。 - Min-Max归一化：
(x - min)/(max - min)。将数据缩放到[0,1]区间，但对异常值敏感。
异常值会严重拉偏簇中心的位置。应对方法：
- 预处理时剔除或缩尾：使用箱线图或3σ原则识别并处理异常值。
- 使用更鲁棒的距离度量：如曼哈顿距离。
- 考虑其他聚类算法：如DBSCAN，它能自动识别噪声点。
4.4 评估聚类结果：没有真实标签怎么办？
在没有真实标签的无监督学习中，评估本身就是个难题。除了前述用于确定k的内部指标（如轮廓系数、Davies-Bouldin Index），还可以：
- 簇内相似性/簇间分离性：计算簇内平均距离（要小），簇间中心点距离（要大）。
- 稳定性分析：对数据做子采样，多次运行聚类，看样本的簇分配是否稳定。一致性高的聚类结果更可靠。
- 业务逻辑校验：将聚类结果（每个簇的样本）拿给业务专家看，根据他们的经验判断分群是否合理、是否有解释性。在数学建模论文中，对每个簇进行特征画像（描述该簇客户的典型特征）是必不可少的环节。
5. 数学建模与进阶应用场景
在数学建模竞赛中，Kmeans很少单独作为最终模型，它通常是特征工程、数据预处理或更大模型中的一个环节。
5.1 特征工程与降维后聚类
高维数据直接聚类会遭遇“维数灾难”，且计算距离的意义会减弱。常见做法是：
- 先用PCA（主成分分析）或t-SNE进行降维，保留主要信息，在低维空间（如2-3维）进行聚类和可视化。
- 将聚类得到的簇标签作为一个新的类别特征，加入到后续的分类或回归模型中。这相当于让模型学习到数据内部的一种分组结构信息。
% 示例：PCA降维后聚类 [coeff, score, latent] = pca(data_normalized); explained_variance = cumsum(latent)./sum(latent); % 选择累积贡献率>95%的主成分数量 n_components = find(explained_variance >= 0.95, 1); data_pca = score(:, 1:n_components); % 在降维后的数据上做Kmeans k = 3; [idx, C] = kmeans(data_pca, k, ‘Replicates‘, 10); % 此时可视化非常方便 figure; gscatter(data_pca(:,1), data_pca(:,2), idx); hold on; plot(C(:,1), C(:,2), ‘kx‘, ‘MarkerSize‘, 15, ‘LineWidth‘, 3); title(‘PCA降维后Kmeans聚类‘);5.2 图像分割与颜色量化
Kmeans在图像处理中一个经典应用是颜色量化（压缩颜色数量）和图像分割。将每个像素的RGB（或Lab）颜色值作为三维特征数据，进行聚类。聚类后，每个像素被赋予其簇中心的颜色值，从而实现用少数几种颜色代表整个图像。
% 简易图像颜色量化示例 img = imread(‘peppers.png‘); img_double = im2double(img); % 转为双精度 [m, n, d] = size(img_double); % 将图像重塑为 m*n 行，3列（RGB）的矩阵 pixel_list = reshape(img_double, m*n, d); k = 16; % 希望压缩到16种颜色 [idx, color_table] = kmeans(pixel_list, k, ‘Replicates‘, 3, ‘MaxIter‘, 200); % 用簇中心颜色替换每个像素的颜色 quantized_pixels = color_table(idx, :); quantized_img = reshape(quantized_pixels, m, n, d); figure; subplot(1,2,1); imshow(img); title(‘原始图像‘); subplot(1,2,2); imshow(quantized_img); title(sprintf(‘颜色量化 (k=%d)‘, k));5.3 时序数据与轨迹聚类
对于时间序列数据，不能直接对原始序列聚类，因为长度可能不同。需要先进行特征提取，如计算序列的统计特征（均值、方差、趋势）、或使用动态时间规整（DTW）作为距离度量，再应用Kmeans。对于轨迹数据（如车辆GPS点），常需要先进行轨迹预处理（分段、采样、对齐），提取方向、速度、曲率等特征，再进行聚类以发现常见的移动模式。
5.4 与层次聚类的结合：二分Kmeans
标准Kmeans需要指定k，且对初始值敏感。一种改进是二分Kmeans：
- 开始时将所有数据视为一个簇。
- 选择当前所有簇中SSE最大的那个簇，对其用Kmeans（k=2）进行二分裂。
- 重复步骤2，直到达到预设的簇数量k。 这种方法有时能产生更稳定的层次化聚类结果，并且可以在分裂过程中通过SSE下降幅度来决定何时停止，从而自动确定k。
6. 常见陷阱、调试技巧与性能优化
即使理解了原理，实际编码和调试中还是会遇到各种问题。这里记录几个我踩过的坑和解决方法。
6.1 空簇（Empty Cluster）问题
在迭代过程中，可能某个簇会失去所有成员。内置的kmeans函数有稳健的处理机制。如果你自己实现，必须添加检查：
- 策略一（简单但可能不稳定）：随机选择一个数据点作为该簇的新中心。
- 策略二（更优）：选择距离当前所有中心最远的那个数据点作为新中心，或者选择SSE贡献最大的那个数据点所在的簇，将其分裂。
6.2 迭代不收敛或震荡
如果设置了MaxIter但算法在达到最大次数前仍未收敛（中心点还在较大幅度变化），可能原因：
- 数据有异常值：异常值会“吸引”中心点，导致中心点不断追逐异常值。检查数据，进行清洗。
- k值选择过大：过于细分的簇可能导致边界点在不同簇间反复横跳。尝试减小k，或使用轮廓系数检查。
- 距离度量不合适：尝试更换距离度量（如从‘sqeuclidean‘换为‘cityblock‘）。
- 增加
MaxIter：对于复杂数据集，可能需要更多迭代。但通常100-1000次足够，再多可能意味着问题本身。
6.3 性能优化：处理大数据集
当数据量很大（数十万以上）时，Kmeans的计算（尤其是距离计算）会成为瓶颈。优化策略：
- 使用内置函数并开启并行：MATLAB的
kmeans支持并行计算。确保Parallel Computing Toolbox已安装，并在调用时使用‘Options‘, statset(‘UseParallel‘, true)。对于多核CPU，加速效果显著。 - 数据采样：如果数据量极大，可以先进行随机采样，在样本上确定k和初始中心，然后再用全部数据进行一次精细聚类（初始中心设为样本聚类得到的中心）。
- 使用更快的算法变种：如Mini-Batch K-Means，它每次迭代只使用一个数据子集来更新中心，速度更快，适用于海量数据。MATLAB未直接提供，但可以自己实现或寻找工具箱。
- 向量化与矩阵运算：如我们手写代码示例所示，避免在循环内对单个样本计算距离，利用MATLAB的矩阵运算能力。
6.4 结果的可复现性
Kmeans的结果依赖于随机初始化。为了确保结果可复现：
- 设置随机数种子：在运行
kmeans前，使用rng(seed)（如rng(42)）固定随机数生成器状态。 - 使用
‘Replicates‘参数：即使设置了种子，单次运行仍可能陷入局部最优。设置‘Replicates‘, 10（或更多）让算法运行多次并返回最佳结果，这是获得稳定、高质量结果的黄金准则。虽然计算时间增加，但结果可靠性大幅提升。
最后，记住Kmeans是一个工具，它有明确的适用场景（球形簇、簇大小密度相近）和局限。当你的数据不符合这些假设时，不要强行使用。不妨试试层次聚类、DBSCAN（基于密度）或高斯混合模型（GMM）。理解算法的“为什么”和“怎么样”，远比记住函数调用更重要。在数学建模中，清晰阐述你选择Kmeans的理由、参数确定的依据、以及对结果局限性的讨论，往往比单纯跑出一个结果更能赢得评委的青睐。
