K-means聚类算法原理和C++实现

给定训练集$\{x^{(1)},...,x^{(m)}\}$，想把这些样本分成不同的子集，即聚类，$x^{(i)}\in\mathbb{R^{n}}$，但是这是个无标签数据集，也就是说我们再聚类的时候不能利用标签信息，所以这是一个无监督学习问题。

k-means聚类算法的流程如下：

1. 随机初始化聚类中心$\mu_{1},\mu_{2},...,\mu_{k}\in\mathbb{R}^{n}$

2. a. 对与每一个聚类中心，计算所有样本到该聚类中心的距离，然后选出距离该聚类中心最近的几个样本作为一类；

　　$c^{(i)}:=\arg\min_{j}||x^{(i)}-\mu_{j}||^{2}$

 

　　这个公式的意思是，某个样本 i 属于哪一类，取决于该样本距离哪一个聚类中心最近，步骤a就是利用这个规则实现。 

　 b. 对上面分成的k类，根据类里面的样本，重新估计该类的中心：

　　　$\mu_{j}:=\frac{\sum_{i=1}^{m}1\{c^{(i)}=j\}x^{(j)}}{\sum_{i=1}^{m}1\{c^{(i)}=j\}}$

　　对于新的聚类中心，重复a，这里1{...}是一个真值判断，例如1{3=2}=0,1{3=3}=1.

　　

　 c. 重复a和b直至收敛

但是k-means真的能保证收敛吗？k-means的目的是选出聚类中心和每一类的样本，定义失真函数：

$J(c,\mu)=\sum_{i=1}^{m}||x^{(i)}-\mu_{c^{(i)}}||^{2}$

 

这个函数衡量的是某个聚类的中心与该类中所有样本距离的平方和，根据上面k-means的算法，可以看出，a 是固定聚类中心，选择该类的样本，b 是样本固定，调整聚类中心，即每次都是固定一个变量，调整另一个变量，所以k-means完全是在针对失真函数 J 坐标上升，这样，J 必然是单调递减，所以J的值必然收敛。在理论上，这种方法可能会使得k-means在一些聚类结果之间产生震荡，即几组不同的 c 和 μ 有着相同的失真函数值，但是这种情况在实际情况中很少出现。

由于失真函数是一个非凸函数，所以坐标上升不能保证该函数全局收敛，即失真函数容易陷入局部收敛。但是大多数情况下，k-means都可以产生不错的结果，如果担心陷入局部收敛，可以多运行几次k-means(采用不同的随机初始聚类中心)，然后从多次结果中选出失真函数最小的聚类结果。

 

下面是一个简单k-means的C++代码，对{1, 2, 3, 11, 12, 13, 21, 22, 23}这9个样本值聚类：

#include<iostream>
#include<cmath>
#include<vector>
#include<ctime>
using namespace std;
typedef unsigned int uint;

struct Cluster
{
    vector<double> centroid;
    vector<uint> samples;
};
double cal_distance(vector<double> a, vector<double> b)
{
    uint da = a.size();
    uint db = b.size();
    if (da != db) cerr << "Dimensions of two vectors must be same!!\n";
    double val = 0.0;
    for (uint i = 0; i < da; i++)
    {
        val += pow((a[i] - b[i]), 2);
    }
    return pow(val, 0.5);
}
vector<Cluster> k_means(vector<vector<double> > trainX, uint k, uint maxepoches)
{
    const uint row_num = trainX.size();
    const uint col_num = trainX[0].size();

    /*初始化聚类中心*/
    vector<Cluster> clusters(k);
    uint seed = (uint)time(NULL); 
 33     for (uint i = 0; i < k; i++)
    {
        srand(seed);
        int c = rand() % row_num;
        clusters[i].centroid = trainX[c];
        seed = rand();
    }

    /*多次迭代直至收敛，本次试验迭代100次*/
    for (uint it = 0; it < maxepoches; it++)
    {
        /*每一次重新计算样本点所属类别之前，清空原来样本点信息*/
        for (uint i = 0; i < k; i++)
        {
            clusters[i].samples.clear();
        }
        /*求出每个样本点距应该属于哪一个聚类*/
        for (uint j = 0; j < row_num; j++)
        {
            /*都初始化属于第0个聚类*/    
            uint c = 0;
            double min_distance = cal_distance(trainX[j],clusters[c].centroid);
            for (uint i = 1; i < k; i++)
            {
                double distance = cal_distance(trainX[j], clusters[i].centroid);
                if (distance < min_distance)
                {
                    min_distance = distance;
                    c = i;
                }
            }
            clusters[c].samples.push_back(j);
        }

        /*更新聚类中心*/
        for (uint i = 0; i < k; i++)
        {
            vector<double> val(col_num, 0.0); 
            for (uint j = 0; j < clusters[i].samples.size(); j++)
            {
                uint sample = clusters[i].samples[j];
                for (uint d = 0; d < col_num; d++)
                {
                    val[d] += trainX[sample][d];
                    if (j == clusters[i].samples.size() - 1)
                        clusters[i].centroid[d] = val[d] / clusters[i].samples.size();
                }
            }
        }
    }
    return clusters;
}

int main()
{
    vector<vector<double> > trainX(9,vector<double>(1,0));
    //对9个数据{1 2 3 11 12 13 21 22 23}聚类
    double data = 1.0;
    for (uint i = 0; i < 9; i++)
    {
        trainX[i][0] = data;
        if ((i+1) % 3 == 0) data += 8;
        else data++;
    }

    /*k-means聚类*/
    vector<Cluster> clusters_out = k_means(trainX, 3, 100);

    /*输出分类结果*/
    for (uint i = 0; i < clusters_out.size(); i++)
    {
        cout << "Cluster " << i << " :" << endl; 

        /*子类中心*/
        cout << "\t" << "Centroid: " << "\n\t\t[ ";
        for (uint j = 0; j < clusters_out[i].centroid.size(); j++)
        {
            cout << clusters_out[i].centroid[j] << " ";
        }
        cout << "]" << endl;

        /*子类样本点*/
        cout << "\t" << "Samples:\n";
        for (uint k = 0; k < clusters_out[i].samples.size(); k++)
        {
            uint c = clusters_out[i].samples[k];
            cout << "\t\t[ ";
            for (uint m = 0; m < trainX[0].size(); m++)
            {
                cout << trainX[c][m] << " ";
            }
            cout << "]\n";
        }
    }
    return 0;
}

 

下面是4次运行结果：

由于数据简单，容易看出第一次和第是三次结果是理想的，而第二次和第四次都是较差出的聚类结果，即上面说的失真函数陷入了局部最优，所以在实践中多次运行，取出较好的聚类结果。