K-Means算法

一、前言

        聚类是观察式学习，在聚类前可以不知道类别甚至不给定类别数量，是无监督学习的一种。相比于分类，聚类不依赖预定义的类和类标号的训练实例, 分类是示例式学习，要求分类前明确各个类别，并断言每个元素映射到一个类别。很多时候上述条件得不到满足，尤其是在处理海量数据的时候，如果通过预处理使得数据满足分类算法的要求，则代价非常大，这时候可以考虑使用聚类算法。

       所谓聚类问题，就是给定一个元素集合D，其中每个元素具有n个可观察属性，使用某种算法将D划分成k个子集，要求每个子集内部的元素之间相异度尽可能低，而不同子集的元素相异度尽可能高。其中每个子集叫做一个簇。

二、相异度的计算

        在介绍聚类问题前先说说相异度的计算，主要包含标量、二元变量、分类变量、序数变量、向量共六种相异度的计算。

        对于第一个也是最常用的一个标量相异度的计算，也就是求点群中心的计算，主要有三种：

1、闵可夫斯基距离公式（MinkowskiDistance）

其中λ可以随意取值，可以是负数，也可以是正数，或是无穷大。

下面两种是在第一种的基础上λ分别取2和1的情况。

2、欧几里德距离公式（Euclidean Distance）

3、曼哈顿距离公式（CityBlock Distance）

下图是三个距离公式逼近中心的方式：

（1）Minkowski Distance  （2）Euclidean Distance （3） CityBlock Distance

第一个图以星形的方式，第二个图以同心圆的方式，第三个图以菱形的方式。其中第一个λ在0-1之间。

 

三、K-Means均值算法

    先看一副K=2的均值算法的过程图，其中灰色是种子点：

K-Means的算法步骤如下：

1.  随机在图中取K（这里K=2）个种子点。

2.  然后对图中的所有点求到这K个种子点的距离，假如点Pi离种子点Si最近，那么Pi属于Si点群。（上图中，我们可以看到A，B属于上面的种子点，C，D，E属于下面中部的种子点）

3.  接下来，我们要移动种子点到属于他的“点群”的中心。（见图上的第三步）

4.  然后重复第2）和第3）步，直到，种子点没有移动（我们可以看到图中的第四步上面的种子点聚合了A，B，C，下面的种子点聚合了D，E）。

算法演示：

http://home.dei.polimi.it/matteucc/Clustering/tutorial_html/AppletKM.html

鼠标左键是初始化点，右键初始化“种子点”，然后勾选“ShowHistory”可以看到一步一步的迭代。

k-means算法一个有趣的应用示例：中国男足近几年到底在亚洲处于几流水平？

参考：http://www.cnblogs.com/leoo2sk/archive/2010/09/20/k-means.html

亚洲15只球队在2005年-2010年间大型杯赛的战绩（由于澳大利亚是后来加入亚足联的，所以这里没有收录）。

其中包括两次世界杯和一次亚洲杯。对数据做了如下预处理：对于世界杯，进入决赛圈则取其最终排名，没有进入决赛圈的，打入预选赛十强赛赋予40，预选赛小组未出线的赋予50。对于亚洲杯，前四名取其排名，八强赋予5，十六强赋予9，预选赛没出现的赋予17。这样做是为了使得所有数据变为标量，便于后续聚类。

      下面先对数据进行[0,1]规格化，下面是规格化后的数据：

      接着用k-means算法进行聚类。设k=3，即将这15支球队分成三个集团。

      现抽取日本、巴林和泰国的值作为三个簇的种子，即初始化三个簇的中心为A：{0.3, 0, 0.19}，B：{0.7,0.76, 0.5}和C：{1, 1, 0.5}。下面，计算所有球队分别对三个中心点的相异度，这里以欧氏距离度量。下面是我用程序求取的结果：

      从做到右依次表示各支球队到当前中心点的欧氏距离，将每支球队分到最近的簇，可对各支球队做如下聚类：

      中国C，日本A，韩国A，伊朗A，沙特A，伊拉克C，卡塔尔C，阿联酋C，乌兹别克斯坦B，泰国C，越南C，阿曼C，巴林B，朝鲜B，印尼C。

      第一次聚类结果：

      A：日本，韩国，伊朗，沙特；

      B：乌兹别克斯坦，巴林，朝鲜；

      C：中国，伊拉克，卡塔尔，阿联酋，泰国，越南，阿曼，印尼。

      下面根据第一次聚类结果，调整各个簇的中心点。

      A簇的新中心点为：{(0.3+0+0.24+0.3)/4=0.21,(0+0.15+0.76+0.76)/4=0.4175, (0.19+0.13+0.25+0.06)/4=0.1575} = {0.21, 0.4175,0.1575}

      用同样的方法计算得到B和C簇的新中心点分别为{0.7, 0.7333, 0.4167}，{1, 0.94, 0.40625}。

      用调整后的中心点再次进行聚类，得到：

      第二次迭代后的结果为：

      中国C，日本A，韩国A，伊朗A，沙特A，伊拉克C，卡塔尔C，阿联酋C，乌兹别克斯坦B，泰国C，越南C，阿曼C，巴林B，朝鲜B，印尼C。

      结果无变化，说明结果已收敛，于是给出最终聚类结果：

      亚洲一流：日本，韩国，伊朗，沙特

      亚洲二流：乌兹别克斯坦，巴林，朝鲜

      亚洲三流：中国，伊拉克，卡塔尔，阿联酋，泰国，越南，阿曼，印尼

      看来数据告诉我们，说国足近几年处在亚洲三流水平真的是没有冤枉他们，至少从国际杯赛战绩是这样的。

      其实上面的分析数据不仅告诉了我们聚类信息，还提供了一些其它有趣的信息，例如从中可以定量分析出各个球队之间的差距，例如，在亚洲一流队伍中，日本与沙特水平最接近，而伊朗则相距他们较远，这也和近几年伊朗没落的实际相符。另外，乌兹别克斯坦和巴林虽然没有打进近两届世界杯，不过凭借预算赛和亚洲杯上的出色表现占据B组一席之地，而朝鲜由于打入了2010世界杯决赛圈而有幸进入B组，可是同样奇迹般夺得2007年亚洲杯的伊拉克却被分在三流，看来亚洲杯冠军的分量还不如打进世界杯决赛圈重啊。其它有趣的信息，有兴趣的朋友可以进一步挖掘。

 

四、K均值算法的三个主要问题及解决方案

存在的问题：

1、    不能保证找到定位聚类中心的最佳方案，但是可以保证收敛到某个解决方案；

2、    无法指出应该使用多少个类别，K的值是初始时给定的，但是，对于同一个数据集，不同的初始种子点个数K可能会产生完全不同的效果，甚至是不合理的；

3、    K均值假设控件的协方差矩阵不会影响结果，或者已经归一化。

解决方法：

1、    多运行几次K均值，每次初始的聚类中心点都不一样，最后选择方差最小的那个结果；

2、    首先将类别数设为1，然后提高类别数，即K值，每次聚类的时候都采用1的解决方法多运行几次。一般情况下，方差会逐渐下降，直到到达一个拐点，保存K值，即为最佳K值；

3、    将数据乘以逆协方差。

 

五、K-Means算法的改进——K-Means++算法

K-Means++算法步骤：

1.  先从我们的数据库随机挑个随机点当“种子点”。

2.  对于每个点，我们都计算其和最近的一个“种子点”的距离D(x)并保存在一个数组里，然后把这些距离加起来得到Sum(D(x))。

3.  然后，再取一个随机值，用权重的方式来取计算下一个“种子点”。这个算法的实现是，先取一个能落在Sum(D(x))中的随机值Random，然后用Random -= D(x)，直到其<=0，此时的点就是下一个“种子点”。

4.  重复第（2）和第（3）步直到所有的K个种子点都被选出来。

5.  进行K-Means算法。

 

代码的实现：

 

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
typedef struct { double x, y; int group; } point_t, *point;
 
double randf(double m)
{
	return m * rand() / (RAND_MAX - 1.);
}
 
point gen_xy(int count, double radius)
{
	double ang, r;
	point p, pt = malloc(sizeof(point_t) * count);
 
	/* note: this is not a uniform 2-d distribution */
	for (p = pt + count; p-- > pt;) {
		ang = randf(2 * M_PI);
		r = randf(radius);
		p->x = r * cos(ang);
		p->y = r * sin(ang);
	}
 
	return pt;
}
 
inline double dist2(point a, point b)
{
	double x = a->x - b->x, y = a->y - b->y;
	return x*x + y*y;
}
 
inline int
nearest(point pt, point cent, int n_cluster, double *d2)
{
	int i, min_i;
	point c;
	double d, min_d;
 
#	define for_n for (c = cent, i = 0; i < n_cluster; i++, c++)
	for_n {
		min_d = HUGE_VAL;
		min_i = pt->group;
		for_n {
			if (min_d > (d = dist2(c, pt))) {
				min_d = d; min_i = i;
			}
		}
	}
	if (d2) *d2 = min_d;
	return min_i;
}
 
void kpp(point pts, int len, point cent, int n_cent)
{
#	define for_len for (j = 0, p = pts; j < len; j++, p++)
	int i, j;
	int n_cluster;
	double sum, *d = malloc(sizeof(double) * len);
 
	point p, c;
	cent[0] = pts[ rand() % len ];
	for (n_cluster = 1; n_cluster < n_cent; n_cluster++) {
		sum = 0;
		for_len {
			nearest(p, cent, n_cluster, d + j);
			sum += d[j];
		}
		sum = randf(sum);
		for_len {
			if ((sum -= d[j]) > 0) continue;
			cent[n_cluster] = pts[j];
			break;
		}
	}
	for_len p->group = nearest(p, cent, n_cluster, 0);
	free(d);
}
 
point lloyd(point pts, int len, int n_cluster)
{
	int i, j, min_i;
	int changed;
 
	point cent = malloc(sizeof(point_t) * n_cluster), p, c;
 
	/* assign init grouping randomly */
	//for_len p->group = j % n_cluster;
 
	/* or call k++ init */
	kpp(pts, len, cent, n_cluster);
 
	do {
		/* group element for centroids are used as counters */
		for_n { c->group = 0; c->x = c->y = 0; }
		for_len {
			c = cent + p->group;
			c->group++;
			c->x += p->x; c->y += p->y;
		}
		for_n { c->x /= c->group; c->y /= c->group; }
 
		changed = 0;
		/* find closest centroid of each point */
		for_len {
			min_i = nearest(p, cent, n_cluster, 0);
			if (min_i != p->group) {
				changed++;
				p->group = min_i;
			}
		}
	} while (changed > (len >> 10)); /* stop when 99.9% of points are good */
 
	for_n { c->group = i; }
 
	return cent;
}
 
void print_eps(point pts, int len, point cent, int n_cluster)
{
#	define W 400
#	define H 400
	int i, j;
	point p, c;
	double min_x, max_x, min_y, max_y, scale, cx, cy;
	double *colors = malloc(sizeof(double) * n_cluster * 3);
 
	for_n {
		colors[3*i + 0] = (3 * (i + 1) % 11)/11.;
		colors[3*i + 1] = (7 * i % 11)/11.;
		colors[3*i + 2] = (9 * i % 11)/11.;
	}
 
	max_x = max_y = -(min_x = min_y = HUGE_VAL);
	for_len {
		if (max_x < p->x) max_x = p->x;
		if (min_x > p->x) min_x = p->x;
		if (max_y < p->y) max_y = p->y;
		if (min_y > p->y) min_y = p->y;
	}
	scale = W / (max_x - min_x);
	if (scale > H / (max_y - min_y)) scale = H / (max_y - min_y);
	cx = (max_x + min_x) / 2;
	cy = (max_y + min_y) / 2;
 
	printf("%%!PS-Adobe-3.0\n%%%%BoundingBox: -5 -5 %d %d\n", W + 10, H + 10);
	printf( "/l {rlineto} def /m {rmoveto} def\n"
		"/c { .25 sub exch .25 sub exch .5 0 360 arc fill } def\n"
		"/s { moveto -2 0 m 2 2 l 2 -2 l -2 -2 l closepath "
		"	gsave 1 setgray fill grestore gsave 3 setlinewidth"
		" 1 setgray stroke grestore 0 setgray stroke }def\n"
	);
	for_n {
		printf("%g %g %g setrgbcolor\n",
			colors[3*i], colors[3*i + 1], colors[3*i + 2]);
		for_len {
			if (p->group != i) continue;
			printf("%.3f %.3f c\n",
				(p->x - cx) * scale + W / 2,
				(p->y - cy) * scale + H / 2);
		}
		printf("\n0 setgray %g %g s\n",
			(c->x - cx) * scale + W / 2,
			(c->y - cy) * scale + H / 2);
	}
	printf("\n%%%%EOF");
	free(colors);
#	undef for_n
#	undef for_len
}
 
#define PTS 100000
#define K 11
int main()
{
	int i;
	point v = gen_xy(PTS, 10);
	point c = lloyd(v, PTS, K);
	print_eps(v, PTS, c, K);
	// free(v); free(c);
	return 0;
}

聚类结果：