请考虑这样一个问题：如何在超大型数据集中识别异常数据项，例如，如何识别可能具有欺骗性的信用卡交易、有风险的贷款应用程序等等。 检测异常数据的一种方法是将数据项分组为类似的聚类，然后在每个聚类中寻找在某种意义上与该聚类中的其他数据项不同的数据项。

有许多不同的聚类分析算法。 k 平均值算法是其中一种使用最久且最广泛的算法。 在本文中，我将介绍 k 平均值算法的工作原理并提供一个完整的 C# 演示程序。 现在有很多独立的数据聚类分析工具，为何要从头开始创建 k 平均值聚类分析代码呢？ 现有的聚类分析工具可能很难或无法集成到软件系统之中，它们可能无法通过定制处理异常情况，而且这些工具可能存在版权或其他知识产权问题。 阅读完本文后，您将能够尝试进行 k 平均值聚类分析并掌握向 .NET 应用程序添加聚类分析功能的基本知识。

为了体会 k 平均值聚类分析并了解我要在本文中讨论的问题，最好先看一下 图 1 。 演示程序首先创建了一个虚拟数据集，其中包含 20 个数据项。 在聚类分析术语中，数据项有时也称为元组。 这里，每个元组表示一个人，并有两个数值属性值：身高（以英寸为单位）和体重（以磅为单位）。 k 平均值算法的一个局限性在于它仅适用于数据元组完全是数字的情况。

现在，这些虚拟数据已载入内存中的数组。 接着，将聚类的数目设置为三个。 虽然有一些高级聚类分析方法可建议要使用的最佳聚类数目，但一般来说，数据聚类分析是一个探索性的过程，往往需要通过反复试验才能找到要使用的最佳聚类数目。 您稍后将看到，k 平均值聚类分析是一个迭代过程。 演示程序有一个变量 maxCount，该变量用于限制主聚类分析循环的执行次数。 这里，该值随意设为 30。

接着，演示程序在后台使用 k 平均值算法将每个数据元组放入三个聚类中的一个。 可以使用许多方法对聚类分析进行编码。 在本例中，聚类分析由一个 int 数组定义，其中数组索引表示元组，而关联的数组值表示从 0 开始的聚类 ID。 因此，在 图 1 中，元组 0 (65.0, 220.0) 分配给聚类 0，元组 1 (73.0, 160.0) 分配给聚类 1，元组 2 (59.0, 110.0) 分配给聚类 2，元组 3 (61.0, 120.0) 也分配给聚类 2，等等。 可以看到，有八个元组分配给聚类 0，五个元组分配给聚类 1，七个元组分配给聚类 2。

随后，演示程序显示按聚类分组的数据。 如果仔细查看聚类数据，可以看出聚类 0 可称为大体重人员聚类，聚类 1 可称为高个人员聚类，而聚类 2 可称为矮个人员聚类。 演示程序通过分析分配给聚类 0 的元组，根据某种标准得出结论：元组 5 (67.0, 240.0) 是最异常的元组。

在随后几部分内容中，我将逐步为您介绍生成 图 1 所示屏幕快照的代码，这样，您将能够根据自己的需要修改该代码。 本文假定您至少具有 C 系列语言的中级编程技能，但对数据聚类分析知识不作要求。 演示程序虽以 C# 编写，但并未采用 OOP 样式，因此如果您愿意，将该演示程序重构为其他语言应当不会太困难。 在本文中，我将提供演示程序的所有源代码。 您也可以从 archive.msdn.microsoft.com/mag201302kmeans 获取源代码。

k 平均值算法

至少在理论上来说，k 平均值算法是相当简单的。 但您将会发现，有些实现细节会有些棘手。 k 平均值算法的核心概念是质心。 在数据聚类分析中，一组数据元组的质心就是该组中最具代表性的那个元组。 通过举例可以很好地解释这个原理。 假设您有三个身高-体重元组，类似 图 1 所示：

[a] (61.0, 100.0)
[b] (64.0, 150.0)
[c] (70.0, 140.0)
其中哪个元组最具代表性？ 一种方法是计算数学平均值元组，然后将最接近该平均值元组的元组选为质心。 因此，在本例中，平均值元组为：
[m] = ((61.0 + 64.0 + 70.0) / 3, (100.0 + 150.0 + 140.0) / 3)
    = (195.0 / 3, 390.0 / 3)
    = (65.0, 130.0)
那么现在，三个元组中的哪一个最接近 (65.0, 130.0)？ 有多种方法可以定义最近元组。 最常用的方法是使用欧氏距离，这也是本演示程序中所用的方法。 从字面上定义，两个元组之间的欧氏距离是指元组每个组分之差的平方和的平方根。 还是举例说明最清楚。 元组 (61.0, 100.0) 与平均值元组 (65.0, 130.0) 之间的欧氏距离为：
dist(m,a) = sqrt((65.0 - 61.0)^2 + (130.0 - 100.0)^2)
          = sqrt(4.0^2 + 30.0^2)
          = sqrt(16.0 + 900.0)
          = sqrt(916.0)
          = 30.27
dist(m,b) = sqrt((65.0 - 64.0)^2 + (130.0 - 150.0)^2)
          = 20.02
dist(m,c) = sqrt((65.0 - 70.0)^2 + (130.0 - 140.0)^2)
          = 11.18
由于三个距离中的最小距离是数学平均值与元组 [c] 之间的距离，因此元组 [c] 就是这三个元组的质心。 您可能希望使用两个元组之间距离的不同定义来试用演示程序，以便了解这些定义对最终生成的聚类分析有何影响。
在确定聚类质心的概念后，k 平均值算法就比较简单了。 伪代码如下：
assign each tuple to a randomly selected cluster
compute the centroid for each cluster
loop until no improvement or until maxCount
  assign each tuple to best cluster
   (the cluster with closest centroid to tuple)
  update each cluster centroid
   (based on new cluster assignments)
end loop
return clustering
如果您在网上搜索，可以找到一些很好的 k 平均值算法在线动画来播放。 图 2 中的图像显示了演示程序所生成的聚类分析。 每个聚类中带圆圈的数据项就是聚类质心。
图 2 聚类数据和质心

程序的整体结构
图 3 列出了图 1 所示演示程序的整体结构（有少许改动）。 我使用 Visual Studio 2010 新建了一个名为 ClusteringKMeans 的 C# 控制台应用程序；任何最近版本的 Visual Studio 也都应适用。 在解决方案资源管理器窗口中，我将文件 Program.cs 重命名为 ClusteringKMeans­Program.cs，这会自动重命名模板生成的类。 我删除了文件顶部不需要的 using 语句。
图 3 程序的整体结构
using System;
namespace ClusteringKMeans
  class ClusteringKMeansProgram
    static void Main(string[] args)
        Console.WriteLine("\nBegin outlier data detection demo\n");
        Console.WriteLine("Loading all (height-weight) data into memory");
        string[] attributes = new string[] { "Height", "Weight" };
        double[][] rawData = new double[20][];
        rawData[0] = new double[] { 65.0, 220.0 };
        rawData[1] = new double[] { 73.0, 160.0 };
        rawData[2] = new double[] { 59.0, 110.0 };
        rawData[3] = new double[] { 61.0, 120.0 };
        rawData[4] = new double[] { 75.0, 150.0 };
        rawData[5] = new double[] { 67.0, 240.0 };
        rawData[6] = new double[] { 68.0, 230.0 };
        rawData[7] = new double[] { 70.0, 220.0 };
        rawData[8] = new double[] { 62.0, 130.0 };
        rawData[9] = new double[] { 66.0, 210.0 };
        rawData[10] = new double[] { 77.0, 190.0 };
        rawData[11] = new double[] { 75.0, 180.0 };
        rawData[12] = new double[] { 74.0, 170.0 };
        rawData[13] = new double[] { 70.0, 210.0 };
        rawData[14] = new double[] { 61.0, 110.0 };
        rawData[15] = new double[] { 58.0, 100.0 };
        rawData[16] = new double[] { 66.0, 230.0 };
        rawData[17] = new double[] { 59.0, 120.0 };
        rawData[18] = new double[] { 68.0, 210.0 };
        rawData[19] = new double[] { 61.0, 130.0 };
        Console.WriteLine("\nRaw data:\n");
        ShowMatrix(rawData, rawData.Length, true);
        int numAttributes = attributes.Length;
        int numClusters = 3;
        int maxCount = 30;
        Console.WriteLine("\nk = " + numClusters + " and maxCount = " + maxCount);
        int[] clustering = Cluster(rawData, numClusters, numAttributes, maxCount);
        Console.WriteLine("\nClustering complete");
        Console.WriteLine("\nClustering in internal format: \n");
        ShowVector(clustering, true);
        Console.WriteLine("\nClustered data:");
        ShowClustering(rawData, numClusters, clustering, true);
        double[] outlier = Outlier(rawData, clustering, numClusters, 0);
        Console.WriteLine("Outlier for cluster 0 is:");
        ShowVector(outlier, true);
        Console.WriteLine("\nEnd demo\n");
      catch (Exception ex)
        Console.WriteLine(ex.Message);
    } // Main
    // 14 short static method definitions here
为简单起见，我使用了一种静态方法，并删除了所有错误检查功能。 演示代码的第一部分设置要进行聚类分析的身高和体重数据。 因为只有 20 个元组，所以我对数据进行了硬编码并将这些数据存储到内存中一个名为 rawData 的数组中。 一般来说，数据将存储在文本文件或 SQL 表中。 这种情况下，必须编写一个 Helper 函数来将数据加载到内存中。 如果数据源太大而无法放入计算机内存，您必须修改演示代码以循环访问外部数据源而不是数据数组。
设置原始数据后，演示程序会调用 Helper 函数 ShowMatrix 来显示数据。 接着，将值 2（身高和体重）、3 和 30 分别赋予变量 num­Attributes、numClusters 和 maxCount。 前面曾提到，maxCount 用于限制主算法处理循环中的迭代次数。 k 平均值算法往往会快速聚合，但您可能有必要使用 maxCount 值试验一下。
所有聚类分析工作均由方法 Cluster 执行。 该方法返回一个 int 数组，后者定义如何将每个元组分配给某一聚类。 完成后，演示程序会显示编码的聚类分析，还会显示按照聚类进行分组的原始数据。
演示程序最后使用方法 Outliers 分析聚类数据，确定离群元组，也就是可能异常的元组。 该方法接受聚类 ID 并返回距聚类质心（最具代表性的元组）最远的数据元组的值（由欧氏距离测量）。 在本例中，对于聚类 0（即大体重人员聚类），离群值元组是 (67.0, 240.0)，也就是体重最重的人员。
计算聚类质心
前面提到，聚类质心是在分配给聚类的元组中最具代表性的元组，而确定聚类质心的一种方法是计算数学平均值元组，然后找到最接近该平均值元组的那个元组。 Helper 方法 UpdateMeans 计算每个聚类的数学平均值元组，如图 4 所示。
图 4：方法 UpdateMeans
static void UpdateMeans(double[][] rawData, int[] clustering,
  double[][] means)
  int numClusters = means.Length;
  for (int k = 0; k < means.Length; ++k)
    for (int j = 0; j < means[k].Length; ++j)
      means[k][j] = 0.0;
  int[] clusterCounts = new int[numClusters];
  for (int i = 0; i < rawData.Length; ++i)
    int cluster = clustering[i];
    ++clusterCounts[cluster];
    for (int j = 0; j < rawData[i].Length; ++j)
      means[cluster][j] += rawData[i][j];
  for (int k = 0; k < means.Length; ++k)
    for (int j = 0; j < means[k].Length; ++j)
      means[k][j] /= clusterCounts[k]; // danger
  return;
方法 UpdateMeans 假定已存在一个数组的数组，该数组名为 means，因而无需创建该数组并返回它。 由于假定数组 means 存在，您可能希望将其设置为一个 ref 参数。 数组 means 是使用 Helper 方法 Allocate 创建的：
static double[][] Allocate(int numClusters, int numAttributes)
  double[][] result = new double[numClusters][];
  for (int k = 0; k < numClusters; ++k)
    result[k] = new double[numAttributes];
  return result;
means 数组中的第一个索引表示聚类 ID，第二个索引表示属性。 例如，如果 means[0][1] = 150.33，则聚类 0 中元组的体重 (1) 值的平均值为 150.33。
方法 UpdateMeans 首先将数组 means 中的现有值清空，循环访问每个数据元组，记下每个聚类中元组的计数，累计每个属性的总和，然后将每个累计总和除以对应的聚类计数。 请注意，如果有任意聚类计数为 0，则该方法将会引发异常，因此您可能需要添加错误检查功能。
方法 ComputeCentroid（如图 5 所示）确定质心值，即最接近于指定聚类的平均元组值的那个元组的值。
图 5 方法 ComputeCentroid
static double[] ComputeCentroid(double[][] rawData, int[] clustering,
  int cluster, double[][] means)
  int numAttributes = means[0].Length;
  double[] centroid = new double[numAttributes];
  double minDist = double.MaxValue;
  for (int i = 0; i < rawData.Length; ++i) // walk thru each data tuple
    int c = clustering[i]; 
    if (c != cluster) continue;
    double currDist = Distance(rawData[i], means[cluster]);
    if (currDist < minDist)
      minDist = currDist;
      for (int j = 0; j < centroid.Length; ++j)
        centroid[j] = rawData[i][j];
  return centroid;
方法 ComputeCentroid 循环访问数据集中的每个元组，并跳过未不在指定聚类中的元组。 对于指定聚类中的每个元组，使用 Helper 方法 Distance 计算该元组与聚类平均值之间的欧氏距离。 将存储并返回最接近平均值（与之距离最小）的元组值。
方法 UpdateCentroids 调用每个聚类的 ComputeCentroid 以便为所有聚类指定质心：
static void UpdateCentroids(double[][] rawData, int[] clustering,
  double[][] means, double[][] centroids)
  for (int k = 0; k < centroids.Length; ++k)
    double[] centroid = ComputeCentroid(rawData, clustering, k, means);
    centroids[k] = centroid;
方法 UpdateCentroids 假定存在一个数组的数组，该数组名为 centroids。 数组 centroids 与数组 means 十分类似： 第一个索引表示聚类 ID，第二个索引表示数据属性。
总之，每个聚类都有一个质心，即该聚类中最具代表性的元组。 计算质心值的方法是：在每个聚类中找到最接近于该聚类中的平均值元组的那个元组。 每个数据元组都将分配给其聚类质心最接近于该元组的聚类。
Distance 函数和数据规范化
方法 ComputeCentroid 调用 Distance 方法来确定哪个数据元组最接近于聚类平均值。 如前所述，欧氏距离是测量元组与平均值之间距离的最常用方法：
static double Distance(double[] tuple, double[] vector)
  double sumSquaredDiffs = 0.0;
  for (int j = 0; j < tuple.Length; ++j)
    sumSquaredDiffs += Math.Pow((tuple[j] - vector[j]), 2);
  return Math.Sqrt(sumSquaredDiffs);
您可能希望考虑使用其他方法来定义距离。 另一种十分常见的方法是使用每个组分之差的绝对值之和。 由于欧氏距离会对差值进行平方，因此较大差值的权重更加高于较小差值。
在 k 平均值聚类分析算法中，与选择距离函数相关的另一个重要因素是数据规范化。 演示程序使用的是未规范化的原始数据。 由于元组体重值通常是像 160.0 这样的值，而元组身高值通常是像 67.0 这样的值，因此体重差异要比身高差异的影响大得多。 许多情况下，除了探究对原始数据的聚类分析之外，在聚类分析之前对原始数据进行规范化也十分有用。 有许多方法可用于对数据进行规范化。 一种常用的方法是：计算每个属性的平均值 (m) 和标准偏差 (sd)，然后计算每个属性值 (v) 的规范化值 nv = (v-m)/sd。
将每个元组分配给一个聚类
有了计算每个聚类的质心的方法后，便可编写一种方法将每个元组分配给一个聚类。 方法 Assign 如图 6 所示。
图 6 方法 Assign
static bool Assign(double[][] rawData, 
  nt[] clustering, double[][] centroids)
  int numClusters = centroids.Length;
  bool changed = false;
  double[] distances = new double[numClusters];
  for (int i = 0; i < rawData.Length; ++i)
    for (int k = 0; k < numClusters; ++k)
      distances[k] = Distance(rawData[i], centroids[k]);
    int newCluster = MinIndex(distances);
    if (newCluster != clustering[i])
      changed = true;
      clustering[i] = newCluster;
  return changed;
方法 Assign 接受质心值的数组并循环访问每个数据元组。 对于每个数据元组，将计算其与每个聚类质心的距离并将该距离存储在一个名为 distances 的局部数组中，其中该数组的索引表示聚类 ID。 随后，Helper 方法 MinIndex 确定数组 distances 中具有最小距离值的索引，即聚类中具有与元组最接近的质心的聚类 ID。
Helper 方法 MinIndex 如下：
static int MinIndex(double[] distances)
  int indexOfMin = 0;
  double smallDist = distances[0];
  for (int k = 0; k < distances.Length; ++k)
    if (distances[k] < smallDist)
      smallDist = distances[k]; indexOfMin = k;
  return indexOfMin;
在 Assign 中，如果计算出的聚类 ID 不同于存储在数组 clustering 中的现有聚类 ID，则将更新数组 clustering，并且切换一个布尔型标志以指示的聚类中至少发生了一项更改。 此标志将用于确定何时停止主算法循环，即超过最大迭代数或聚类中未发生任何更改时。
k 平均值算法的这一实现假定至少总是有一个数据元组分配给每个聚类。 如图 6 所示，方法 Assign 并不能避免聚类未获分配元组的情况。 实际上，这通常不是问题。 但要避免这一错误情况却有些棘手。 我通常使用的方法是创建一个名为 centroidIndexes 的数组，该数组与 centroids 数组配合使用。 前面曾提到，数组 centroids 存放质心值，例如，在图 2 中 (61.0, 120.0) 是聚类 2 的质心。 数组 centroidIndexes 存放关联的元组索引，例如 [3]。 然后，在 Assign 方法中，第一步是将存放质心值的数据元组分配给每个聚类，然后该方法循环访问其余每个元组并将其分配给一个聚类。 这种方法可保证每个聚类都至少有一个元组。
Cluster 方法
如图 7 所示，方法 Cluster 是调用所有 Helper 和子 Helper 方法来实际执行数据聚类分析的高级例程。
图 7 Cluster 方法
static int[] Cluster(double[][] rawData, int numClusters,
  int numAttributes,  int maxCount)
  bool changed = true;
  int ct = 0;
  int numTuples = rawData.Length;
  int[] clustering = InitClustering(numTuples, numClusters, 0);
  double[][] means = Allocate(numClusters, numAttributes);
  double[][] centroids = Allocate(numClusters, numAttributes);
  UpdateMeans(rawData, clustering, means);
  UpdateCentroids(rawData, clustering, means, centroids);
  while (changed == true && ct < maxCount)
    ++ct;
    changed = Assign(rawData, clustering, centroids);
    UpdateMeans(rawData, clustering, means);
    UpdateCentroids(rawData, clustering, means, centroids);
  return clustering;
主 while 循环反复将每个数据元组分配给聚类，计算每个聚类的新元组平均值，然后使用新的平均值计算每个聚类的新质心值。 当聚类分配中不再发生更改或达到某一最大计数时，该循环便退出。 由于 means 数组仅用于计算质心，因此您可能希望通过在方法 UpdateCentroids 中调用 UpdateMeans 来重构聚类。
在开始处理循环之前，将使用方法 InitClustering 初始化 clustering 数组。
static int[] InitClustering(int numTuples, 
  int numClusters, int randomSeed)
  Random random = new Random(randomSeed);
  int[] clustering = new int[numTuples];
  for (int i = 0; i < numClusters; ++i)
    clustering[i] = i;
  for (int i = numClusters; i < clustering.Length; ++i)
    clustering[i] = random.Next(0, numClusters);
  return clustering;
InitClustering 方法首先将元组 0 到 num­Clusters-1 分别分配给聚类 0 到 numClusters-1，这样每个聚类最初至少分配有一个元组。 其余元组将分配给随机选择的聚类。
历来对 k 平均值聚类分析初始化的研究之多有些出人意料，因此除了此处提供的方法之外，您可能也希望试用一下其他方法。 许多情况下，k 平均值算法所生成的最终聚类分析取决于初始化聚类分析的方式。
查找异常数据
使用数据聚类分析的一种方法是只探究不同的聚类，然后查找异常或出人意料的结果。 另一种可能性是在聚类中查找异常数据元组。 演示程序将检查聚类 0，使用名为 Outlier 的方法在该聚类中找到离聚类质心最远的元组，该方法如图 8 所示。
图 8 Outlier 方法
static double[] Outlier(double[][] rawData, int[] clustering,
  int numClusters, int cluster)
  int numAttributes = rawData[0].Length;
  double[] outlier = new double[numAttributes];
  double maxDist = 0.0;
  double[][] means = Allocate(numClusters, numAttributes);
  double[][] centroids = Allocate(numClusters, numAttributes);
  UpdateMeans(rawData, clustering, means);
  UpdateCentroids(rawData, clustering, means, centroids);
  for (int i = 0; i < rawData.Length; ++i)
    int c = clustering[i];
    if (c != cluster) continue;
    double dist = Distance(rawData[i], centroids[cluster]);
    if (dist > maxDist)
      maxDist = dist;
      Array.Copy(rawData[i], outlier, rawData[i].Length);
  return outlier;
在初始化 means 和 centroids 数组之后，方法 Outlier 循环访问指定聚类中的每个元组，计算该元组与聚类质心的欧氏距离，然后返回与质心值距离最远的元组的值。 您也可以考虑返回最远数据元组的索引。
还有许多其他方法可用于检查聚类数据的异常情况。 例如，您可能希望确定每个元组与其已分配的聚类质心之间的平均距离，或者可能希望检查聚类质心彼此之间的距离。
为完整起见，这里提供一些简化的显示例程。 通过代码下载可获得更丰富的版本。 如果使用这些简化的例程，您必须在 Main 方法中修改其调用。 要显示原始数据、平均值和质心，您可以使用：
static void ShowMatrix(double[][] matrix)
  for (int i = 0; i < numRows; ++i)
    Console.Write("[" + i.ToString().PadLeft(2) + "]  ");
    for (int j = 0; j < matrix[i].Length; ++j)
      Console.Write(matrix[i][j].ToString("F1") + "  ");
    Console.WriteLine("");
要显示 clustering 数组，您可以使用：
static void ShowVector(int[] vector)
  for (int i = 0; i < vector.Length; ++i)
    Console.Write(vector[i] + " ");
  Console.WriteLine("");
要显示 outlier 的值，您可以使用：
static void ShowVector(double[] vector)
  for (int i = 0; i < vector.Length; ++i)
    Console.Write(vector[i].ToString("F1") + " ");
  Console.WriteLine("");
此外，要显示按聚类分组的原始数据，您可以使用：
static void ShowClustering(double[][] rawData, 
  int numClusters, int[] clustering)
  for (int k = 0; k < numClusters; ++k) // Each cluster
    for (int i = 0; i < rawData.Length; ++i) // Each tuple
      if (clustering[i] == k)
        for (int j = 0; j < rawData[i].Length; ++j)
          Console.Write(rawData[i][j].ToString("F1") + " ");
        Console.WriteLine("");
    Console.WriteLine("");
数据聚类分析与数据分类紧密相关，有时也会与之混淆。 聚类分析是一种不受监控的方法，它将数据项分组在一起而预先不知道分组结果。 聚类分析通常是一个探索性的过程。 相比之下，分类是一种受监控的方法，它需要在训练数据中指定已知组，随后将每个数据元组放入其中某个组中。 分类通常用于预测目的。
本文中提供的代码和说明应该已经给您提供了足够的信息，使您可以体验 k 平均值数据聚类分析，创建完全可自定义的独立聚类分析工具，或者向 .NET 应用程序添加聚类分析功能而无需依赖任何外部因素。 除了 k 平均值算法以外，还有许多其他聚类分析算法。我将在以后的 MSDN 杂志文章中介绍其中一些算法，包括数据熵最小化、类别实用工具和 Naive Bayes 推理等。
James McCaffrey博士 供职于 Volt Information Sciences, Inc.，在该公司他负责管理对华盛顿州雷蒙德市沃什湾 Microsoft 总部园区的软件工程师进行的技术培训。 他参与过多项 Microsoft 产品的研发工作，其中包括 Internet Explorer 和 MSN Search。 他是《.NET Test Automation Recipes》(Apress, 2006) 的作者，您可以通过以下电子邮箱地址与他联系：jammc@microsoft.com。
衷心感谢以下技术专家对本文的审阅： Darren Gehring