Main Content

本页的翻译已过时。点击此处可查看最新英文版本。

kmeans

k 均值聚类

说明

示例

idx = kmeans(X,k) 执行 k 均值聚类,以将 n×p 数据矩阵 X 的观测值划分为 k 个聚类,并返回包含每个观测值的簇索引的 n×1 向量 (idx)。X 的行对应于点,列对应于变量。

默认情况下,kmeans 使用欧几里德距离平方度量,并用 k-means++ 算法进行簇中心初始化。

示例

idx = kmeans(X,k,Name,Value) 进一步按一个或多个 Name,Value 对组参数所指定的附加选项返回簇索引。

例如,指定余弦距离、使用新初始值重复聚类的次数或使用并行计算的次数。

示例

[idx,C] = kmeans(___)k×p 矩阵 C 中返回 k 个簇质心的位置。

示例

[idx,C,sumd] = kmeans(___)k×1 向量 sumd 中返回簇内的点到质心距离的总和。

示例

[idx,C,sumd,D] = kmeans(___) 在 n×k 矩阵 D 中返回每个点到每个质心的距离。

示例

全部折叠

对数据进行 k 均值聚类,然后绘制簇区域。

加载 Fisher 鸢尾花数据集。使用花瓣长度和宽度作为预测变量。

load fisheriris
X = meas(:,3:4);

figure;
plot(X(:,1),X(:,2),'k*','MarkerSize',5);
title 'Fisher''s Iris Data';
xlabel 'Petal Lengths (cm)'; 
ylabel 'Petal Widths (cm)';

Figure contains an axes. The axes with title Fisher's Iris Data contains an object of type line.

较大的簇似乎被分成两个区域:一个较低方差区域,一个较高方差区域。这可能表明较大的簇是两个重叠的簇。

对数据进行聚类。指定 k = 3 个簇。

rng(1); % For reproducibility
[idx,C] = kmeans(X,3);

idx 是与 X 中的观测值对应的预测簇索引的向量。C 是包含最终质心位置的 3×2 矩阵。

使用 kmeans 计算从每个质心到网格上各点的距离。为此,将质心 (C) 和网格上的点传递给 kmeans,并实现算法的一次迭代。

x1 = min(X(:,1)):0.01:max(X(:,1));
x2 = min(X(:,2)):0.01:max(X(:,2));
[x1G,x2G] = meshgrid(x1,x2);
XGrid = [x1G(:),x2G(:)]; % Defines a fine grid on the plot

idx2Region = kmeans(XGrid,3,'MaxIter',1,'Start',C);
Warning: Failed to converge in 1 iterations.
    % Assigns each node in the grid to the closest centroid

kmeans 显示一条警告,指出算法未收敛,这是您应预料到的,因为软件只实现了一次迭代。

绘制簇区域。

figure;
gscatter(XGrid(:,1),XGrid(:,2),idx2Region,...
    [0,0.75,0.75;0.75,0,0.75;0.75,0.75,0],'..');
hold on;
plot(X(:,1),X(:,2),'k*','MarkerSize',5);
title 'Fisher''s Iris Data';
xlabel 'Petal Lengths (cm)';
ylabel 'Petal Widths (cm)'; 
legend('Region 1','Region 2','Region 3','Data','Location','SouthEast');
hold off;

Figure contains an axes. The axes with title Fisher's Iris Data contains 4 objects of type line. These objects represent Region 1, Region 2, Region 3, Data.

随机生成样本数据。

rng default; % For reproducibility
X = [randn(100,2)*0.75+ones(100,2);
    randn(100,2)*0.5-ones(100,2)];

figure;
plot(X(:,1),X(:,2),'.');
title 'Randomly Generated Data';

Figure contains an axes. The axes with title Randomly Generated Data contains an object of type line.

数据中似乎有两个簇。

将数据分成两个簇,并从五个初始化中选择最佳排列。显示最终输出。

opts = statset('Display','final');
[idx,C] = kmeans(X,2,'Distance','cityblock',...
    'Replicates',5,'Options',opts);
Replicate 1, 3 iterations, total sum of distances = 201.533.
Replicate 2, 5 iterations, total sum of distances = 201.533.
Replicate 3, 3 iterations, total sum of distances = 201.533.
Replicate 4, 3 iterations, total sum of distances = 201.533.
Replicate 5, 2 iterations, total sum of distances = 201.533.
Best total sum of distances = 201.533

默认情况下,软件使用 k-means++ 分别对每次重复进行初始化。

绘制簇和簇质心。

figure;
plot(X(idx==1,1),X(idx==1,2),'r.','MarkerSize',12)
hold on
plot(X(idx==2,1),X(idx==2,2),'b.','MarkerSize',12)
plot(C(:,1),C(:,2),'kx',...
     'MarkerSize',15,'LineWidth',3) 
legend('Cluster 1','Cluster 2','Centroids',...
       'Location','NW')
title 'Cluster Assignments and Centroids'
hold off

Figure contains an axes. The axes with title Cluster Assignments and Centroids contains 3 objects of type line. These objects represent Cluster 1, Cluster 2, Centroids.

通过将 idx 传递给 silhouette,您可以确定簇之间的分离程度。

对大型数据集进行聚类可能需要大量时间,尤其是在您使用在线更新(默认设置)时。如果您拥有 Parallel Computing Toolbox™ 许可证并设置了并行计算的选项,则 kmeans 将并行运行每个聚类任务(或副本)。而且,如果 Replicates > 1,则并行计算会减少收敛时间。

从高斯混合模型中随机生成一个大型数据集。

Mu = bsxfun(@times,ones(20,30),(1:20)'); % Gaussian mixture mean
rn30 = randn(30,30);
Sigma = rn30'*rn30; % Symmetric and positive-definite covariance
Mdl = gmdistribution(Mu,Sigma); % Define the Gaussian mixture distribution

rng(1); % For reproducibility
X = random(Mdl,10000);

Mdl 是一个 30 维 gmdistribution 模型,包含 20 个分量。X 是一个 10000×30 矩阵,其数据出自 Mdl 模型。

指定并行计算的选项。

stream = RandStream('mlfg6331_64');  % Random number stream
options = statset('UseParallel',1,'UseSubstreams',1,...
    'Streams',stream);

RandStream 的输入参数 'mlfg6331_64' 指定使用乘法滞后 Fibonacci 生成器算法。options 是一个结构体数组,其字段指定控制估算的选项。

对数据进行 k 均值聚类。指定数据中有 k = 20 个簇,并增加迭代次数。通常,目标函数包含局部最小值。指定 10 个副本以帮助找到更低的局部最小值。

tic; % Start stopwatch timer
[idx,C,sumd,D] = kmeans(X,20,'Options',options,'MaxIter',10000,...
    'Display','final','Replicates',10);
Starting parallel pool (parpool) using the 'local' profile ...
connected to 6 workers.
Replicate 5, 72 iterations, total sum of distances = 7.73161e+06.
Replicate 1, 64 iterations, total sum of distances = 7.72988e+06.
Replicate 3, 68 iterations, total sum of distances = 7.72576e+06.
Replicate 4, 84 iterations, total sum of distances = 7.72696e+06.
Replicate 6, 82 iterations, total sum of distances = 7.73006e+06.
Replicate 7, 40 iterations, total sum of distances = 7.73451e+06.
Replicate 2, 194 iterations, total sum of distances = 7.72953e+06.
Replicate 9, 105 iterations, total sum of distances = 7.72064e+06.
Replicate 10, 125 iterations, total sum of distances = 7.72816e+06.
Replicate 8, 70 iterations, total sum of distances = 7.73188e+06.
Best total sum of distances = 7.72064e+06
toc % Terminate stopwatch timer
Elapsed time is 61.915955 seconds.

命令行窗口指示有六个工作进程可用。您的系统中的工作进程数量可能会有所不同。命令行窗口显示迭代次数和每个副本的最终目标函数值。输出参数包含副本 9 的结果,因为它的总距离最低。

kmeans 执行 k 均值聚类以将数据划分为 k 个簇。当您有要进行聚类的新数据集时,可以使用 kmeans 创建包含现有数据和新数据的新簇。kmeans 函数支持 C/C++ 代码生成,因此您可以生成接受训练数据并返回聚类结果的代码,然后将代码部署到设备上。在此工作流中,您必须传递训练数据,训练数据有可能相当大。为了节省设备上的内存,您可以分别使用 kmeanspdist2 来分离训练和预测。

使用 kmeans 在 MATLAB® 中创建簇,并在生成的代码中使用 pdist2 将新数据分配给现有簇。对于代码生成,定义接受簇质心位置和新数据集的入口函数,并返回最近邻簇的索引。然后,为入口函数生成代码。

生成 C/C++ 代码需要 MATLAB® Coder™。

执行 k 均值聚类

使用三种分布生成训练数据集。

rng('default') % For reproducibility
X = [randn(100,2)*0.75+ones(100,2);
    randn(100,2)*0.5-ones(100,2);
    randn(100,2)*0.75];

使用 kmeans 将训练数据分成三个簇。

[idx,C] = kmeans(X,3);

绘制簇和簇质心。

figure
gscatter(X(:,1),X(:,2),idx,'bgm')
hold on
plot(C(:,1),C(:,2),'kx')
legend('Cluster 1','Cluster 2','Cluster 3','Cluster Centroid')

Figure contains an axes. The axes contains 4 objects of type line. These objects represent Cluster 1, Cluster 2, Cluster 3, Cluster Centroid.

将新数据分配给现有簇

生成测试数据集。

Xtest = [randn(10,2)*0.75+ones(10,2);
    randn(10,2)*0.5-ones(10,2);
    randn(10,2)*0.75];

使用现有簇对测试数据集进行分类。使用 pdist2 找到距离每个测试数据点最近的质心。

[~,idx_test] = pdist2(C,Xtest,'euclidean','Smallest',1);

使用 idx_testgscatter 绘制测试数据并对测试数据加标签。

gscatter(Xtest(:,1),Xtest(:,2),idx_test,'bgm','ooo')
legend('Cluster 1','Cluster 2','Cluster 3','Cluster Centroid', ...
    'Data classified to Cluster 1','Data classified to Cluster 2', ...
    'Data classified to Cluster 3')

Figure contains an axes. The axes contains 7 objects of type line. These objects represent Cluster 1, Cluster 2, Cluster 3, Cluster Centroid, Data classified to Cluster 1, Data classified to Cluster 2, Data classified to Cluster 3.

生成代码

生成将新数据分配给现有簇的 C 代码。请注意,生成 C/C++ 代码需要 MATLAB® Coder™。

定义名为 findNearestCentroid 的入口函数,该函数接受质心位置和新数据,然后使用 pdist2 找到最近的簇。

在入口函数的函数签名后面添加 %#codegen 编译器指令(即 pragma),以指示您要为此 MATLAB 算法生成代码。添加此指令指示 MATLAB 代码分析器帮助您诊断和修复在代码生成过程中可能导致错误的违规。

type findNearestCentroid % Display contents of findNearestCentroid.m
function idx = findNearestCentroid(C,X) %#codegen
[~,idx] = pdist2(C,X,'euclidean','Smallest',1); % Find the nearest centroid

注意:如果您点击位于此页右上角的按钮,并在 MATLAB® 中打开此示例,则 MATLAB® 将打开示例文件夹。该文件夹包括入口函数文件。

使用 codegen (MATLAB Coder) 生成代码。由于 C 和 C++ 是静态类型语言,因此必须在编译时确定入口函数中所有变量的属性。要指定 findNearestCentroid 的输入的数据类型和数组大小,请使用 -args 选项传递表示具有特定数据类型和数组大小的值集的 MATLAB 表达式。有关详细信息,请参阅 Specify Variable-Size Arguments for Code Generation

codegen findNearestCentroid -args {C,Xtest}
Code generation successful.

codegen 生成 MEX 函数 findNearestCentroid_mex,扩展名因平台而异。

验证生成的代码。

myIndx = findNearestCentroid(C,Xtest);
myIndex_mex = findNearestCentroid_mex(C,Xtest);
verifyMEX = isequal(idx_test,myIndx,myIndex_mex)
verifyMEX = logical
   1

isequal 返回逻辑值 1 (true),这意味着所有输入都相等。这一比较结果确认 pdist2 函数、findNearestCentroid 函数和 MEX 函数均返回相同的索引。

您还可以使用 GPU Coder™ 生成优化的 CUDA® 代码。

cfg = coder.gpuConfig('mex');
codegen -config cfg findNearestCentroid -args {C,Xtest}

有关代码生成的详细信息,请参阅General Code Generation Workflow。有关 GPU Coder 的详细信息,请参阅Get Started with GPU Coder (GPU Coder)Supported Functions (GPU Coder)

输入参数

全部折叠

数据,指定为数值矩阵。X 的行对应于观测值,而列对应于变量。

如果 X 是数值向量,则 kmeans 将其视为 n×1 数据矩阵,而不管其方向如何。

该软件将 X 中的 NaN 视为缺失数据,并删除 X 中包含至少一个 NaN 的任何行。删除 X 的行会缩小样本大小。对于输出参数 idx 中的对应值,kmeans 函数返回 NaN

数据类型: single | double

数据中的簇数,指定为一个正整数。

数据类型: single | double

名称-值对组参数

指定可选的、以逗号分隔的 Name,Value 对组参数。Name 为参数名称,Value 为对应的值。Name 必须放在引号内。您可采用任意顺序指定多个名称-值对组参数,如 Name1,Value1,...,NameN,ValueN 所示。

示例: 'Distance','cosine','Replicates',10,'Options',statset('UseParallel',1) 指定使用余弦距离,以不同起始值重复 10 次聚类,并使用并行计算。

要在命令行窗口中显示的输出的级别,指定为以逗号分隔的对组,其中包含 'Display' 和以下选项之一:

  • 'final' - 显示最终迭代的结果

  • 'iter' - 显示每次迭代的结果

  • 'off' - 不显示任何内容

示例: 'Display','final'

p 维空间中的距离度量,用于最小化距离和,指定为以逗号分隔的对组,其中包含 'Distance''sqeuclidean''cityblock''cosine''correlation''hamming'

对于支持的各种距离度量,kmeans 分别以不同方式计算质心簇。下表总结了可用的距离度量。在公式中,x 是一个观测值(即 X 的一个行),c 是一个质心(一个行向量)。

距离度量说明公式
'sqeuclidean'

欧几里德距离的平方(默认值)。每个质心是该簇中各点的均值。

d(x,c)=(xc)(xc)

'cityblock'

绝对差之和,即 L1 距离。每个质心是该簇中各点的按成分中位数。

d(x,c)=j=1p|xjcj|

'cosine'

1 减去点之间夹角的余弦值(视为向量)。每个质心是该簇中各点的均值(在将这些点归一化为单位欧几里德长度之后)。

d(x,c)=1xc(xx)(cc)

'correlation'

1 减去点之间的样本相关性(视为值序列)。每个质心是该簇中各点的按成分均值(在将这些点中心化并归一化为零均值和单位标准差后)。

其中

d(x,c)=1(xx¯)(cc¯)(xx¯)(xx¯)(cc¯)(cc¯),

  • x¯=1p(j=1pxj)1p

  • c¯=1p(j=1pcj)1p

  • 1p 是 p 个 1 组成的行向量。

'hamming'

此度量仅适用于二类数据。

它是相异位所占的比例。每个质心是该簇中各点的按成分中位数。

d(x,y)=1pj=1pI{xjyj},

其中 I 是指示函数。

示例: 'Distance','cityblock'

在簇丢失其所有成员观测值时应采取的操作,指定为以逗号分隔的对组,其中包含 'EmptyAction' 和以下选项之一。

说明
'error'

将空簇视为错误。

'drop'

删除所有变空的簇。kmeansCD 中对应的返回值设置为 NaN

'singleton'

创建一个新簇,其中包含一个距其质心最远的点(默认值)。

示例: 'EmptyAction','error'

最大迭代次数,指定为以逗号分隔的对组,其中包含 'MaxIter' 和一个正整数。

示例: 'MaxIter',1000

数据类型: double | single

在线更新标志,指定为以逗号分隔的对组,其中包含 'OnlinePhase''off''on'

如果 OnlinePhaseon,则 kmeans 除了执行批量更新阶段之外,还会执行在线更新阶段。对于大型数据集来说,在线阶段可能很耗时,但可以保证解是距离标准的局部最小值。换句话说,软件会找到一个数据分区,在此分区中,将任一点移至不同簇都会增加总距离。

示例: 'OnlinePhase','on'

控制迭代算法以最小化拟合标准的选项,指定为以逗号分隔的对组,其中包含 'Options'statset 返回的结构体数组。结构体数组中受支持的字段指定用于控制迭代算法的选项。

下表总结了支持的字段。请注意,支持的字段需要 Parallel Computing Toolbox™。

字段说明
'Streams'

RandStream 对象或此类对象的元胞数组。如果您没有指定 Streamskmeans 将使用默认流。如果您指定 Streams,请使用单个对象,除非满足以下所有条件:

  • 您有一个打开的并行池。

  • UseParalleltrue

  • UseSubstreamsfalse

在这种情况下,请使用与并行池大小相同的元胞数组。如果并行池未打开,则 Streams 必须提供一个随机数流。

'UseParallel'
  • 如果为 trueReplicates > 1,则 kmeans 以并行方式对每次重复实现 k 均值算法。

  • 如果未安装 Parallel Computing Toolbox,则计算以串行模式进行。默认值为 false,表示串行计算。

'UseSubstreams'设置为 true 将以可重现的方式执行并行计算。默认值为 false。要以可重现方式执行计算,请将 Streams 设置为允许子流的类型:'mlfg6331_64''mrg32k3a'

为确保更具预测性的结果,请在调用 kmeans 并设置 'Options',statset('UseParallel',1) 之前,使用 parpool (Parallel Computing Toolbox) 并显式创建并行池。

示例: 'Options',statset('UseParallel',1)

数据类型: struct

使用新初始簇质心位置重复聚类的次数,指定为由 'Replicates' 和整数组成的以逗号分隔的对组。kmeans 返回 sumd 最低的解。

您可以通过提供三维数组作为 'Start' 名称-值对组参数的值来隐式设置 'Replicates'

示例: 'Replicates',5

数据类型: double | single

一种选择初始簇质心位置(即种子)的方法,指定为以逗号分隔的对组,其中包含 'Start' 以及 'cluster''plus''sample''uniform'、一个数值矩阵,或一个数值数组。下表总结了选择种子的可用选项。

说明
'cluster'

X 的一个随机 10% 子样本中的观测值数目多于 k 个时,对该子样本执行初步聚类阶段。此初步阶段本身是使用 'sample' 进行初始化的。

如果随机 10% 子样本中的观测值数目少于 k 个,则软件从 X 中随机选择 k 个观测值。

'plus'(默认值)实现用于簇中心初始化的 k-means++ 算法,依算法选择 k 个种子。
'sample'X 中随机选择 k 个观测值。
'uniform'X 范围内随机均匀选择 k 个点。Hamming 距离无效。
数值矩阵包含质心起始位置的 k×p 矩阵。Start 的行对应于种子。软件根据 Start 的第一个维度推断 k,因此您可以为 k 传入 []
数值数组包含质心起始位置的 k×p×r 数组。每页的行对应于种子。第三个维度调用聚类例程的重复。第 j 页包含副本 j 的种子集。软件根据第三个维度的大小推断副本的数量(由 'Replicates' 名称-值对组参数指定)。

示例: 'Start','sample'

数据类型: char | string | double | single

输出参数

全部折叠

簇索引,以数值列向量形式返回。idx 的行数与 X 的行数一样多,每行都表示对应观测值的簇分配。

簇质心位置,以数值矩阵形式返回。Ck×p 矩阵,其中第 j 行是簇 j 的质心。

簇内的点到质心距离的总和,以数值列向量形式返回。sumdk×1 向量,其中元素 j 是簇 j 内的点到质心距离的总和。默认情况下,kmeans 使用欧几里德距离平方(请参阅 'Distance' 度量)。

从每个点到每个质心的距离,以数值矩阵形式返回。D 是 n×k 矩阵,其中元素 (j,m) 是从观测值 j 到质心 m 的距离。默认情况下,kmeans 使用欧几里德距离平方(请参阅 'Distance' 度量)。

详细信息

全部折叠

k 均值聚类

k 均值聚类,即 Lloyd 算法 [2],是一种迭代数据划分算法,它将 n 个观测值分配给由质心定义的 k 个簇之一,其中 k 是在算法开始之前选择的。

算法的执行如下:

  1. 选择 k 个初始簇中心(质心)。例如,随机选择 k 个观测值(通过使用 'Start','sample')或使用 k-means ++ 算法进行簇中心初始化(默认值)。

  2. 计算所有观测值到每个质心的点到簇质心的距离。

  3. 有两种方法可以继续(由 OnlinePhase 指定):

    • 批量更新 - 将每个观测值分配给离质心最近的簇。

    • 在线更新 - 只要将观测值重新分配给另一质心可减少簇内点到质心距离平方和的总和,就对该观测值执行此分配。

    有关详细信息,请参阅算法

  4. 计算每个簇中观测值的平均值,以获得 k 个新质心位置。

  5. 重复步骤 2 到 4,直到簇分配不变,或达到最大迭代次数。

k-means++ 算法

k-means++ 算法使用启发式方法找到 k 均值聚类的质心种子。根据 Arthur 和 Vassilvitskii 的研究·[1],k-means++ 改进了 Lloyd 算法的运行时间和最终解的质量。

k-means++ 算法按如下方式选择种子,假设簇数为 k。

  1. 从数据集 X 中随机均匀选择一个观测值。所选观测值是第一个质心,表示为 c1

  2. 计算从每个观测值到 c1 的距离。将 cj 和观测值 m 之间的距离表示为 d(xm,cj)

  3. 从 X 中随机选择下一个质心 c2,概率为

    d2(xm,c1)j=1nd2(xj,c1).

  4. 要选择中心 j,请执行以下操作:

    1. 计算从每个观测值到每个质心的距离,并将每个观测值分配给其最近的质心。

    2. 对于 m = 1、...、n 和 p = 1、...、j - 1,从 X 中随机选择质心 j,概率为

      d2(xm,cp){h;xhCp}d2(xh,cp),

      其中 Cp 是最接近质心 cp 的所有观测值的集合,而 xm 属于 Cp

      也就是说,选择每个后续中心时,其选择概率与它到已选最近中心的距离成比例。

  5. 重复步骤 4,直到选择了 k 个质心。

Arthur 和 Vassilvitskii·[1] 通过对几个簇方向的模拟研究证明,在计算簇内点到质心距离平方和时,k-means++ 相比 Lloyd 算法能更快地收敛至更低的总和。

算法

  • kmeans 使用两阶段迭代算法来最小化点到质心的距离总和,该距离总和覆盖所有 k 个簇。

    1. 第一阶段使用批量更新,其中每次迭代都包括一次性将点重新分配给其最近邻的簇质心,然后重新计算簇质心。此阶段有时候不会收敛至局部最小值解。也就是说,将任一点移至不同簇都会增加总距离的数据分区。尤其是对小数据集来说,更容易出现这种情况。批量更新阶段速度很快,但在该阶段逼近的解可能只能作为第二阶段的起点。

    2. 第二阶段使用在线更新,即只要重新分配点可减少距离的总和,则进行重新分配,并在每次重新分配后重新计算簇质心。此阶段中的每次迭代都对所有点进行一次遍历。此阶段会收敛至一个局部最小值,尽管可能存在总距离更低的其他局部最小值。一般情况下,需要穷尽选择起始点才能求解全局最小值,但是,使用具有随机起始点的几个副本通常也会得到全局最小值解。

  • 如果 Replicates = r > 1 且 Startplus(默认值),则软件根据 k-means++ 算法选择 r 个可能不同的种子集。

  • 如果您在 Options 中启用 UseParallel 选项且 Replicates > 1,则每个工作进程会以并行方式选择种子和进行聚类。

参考

[1] Arthur, David, and Sergi Vassilvitskii. “K-means++: The Advantages of Careful Seeding.” SODA ‘07: Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms. 2007, pp. 1027–1035.

[2] Lloyd, Stuart P. “Least Squares Quantization in PCM.” IEEE Transactions on Information Theory. Vol. 28, 1982, pp. 129–137.

[3] Seber, G. A. F. Multivariate Observations. Hoboken, NJ: John Wiley & Sons, Inc., 1984.

[4] Spath, H. Cluster Dissection and Analysis: Theory, FORTRAN Programs, Examples. Translated by J. Goldschmidt. New York: Halsted Press, 1985.

扩展功能

在 R2006a 之前推出