1.基本Kmeans演算法[1]

選擇K個點作為初始質心
repeat
	將每個點指派到最近的質心，形成K個簇
	重新計算每個簇的質心
until 簇不發生變化或達到最大迭代次數

時間複雜度：O(tKmn)，其中，t為迭代次數，K為簇的數目，m為記錄數，n為維數

空間複雜度：O((m+K)n)，其中，K為簇的數目，m為記錄數，n為維數

2.注意問題

（1）K如何確定

        kmenas演算法首先選擇K個初始質心，其中K是使用者指定的引數，即所期望的簇的個數。這樣做的前提是我們已經知道資料集中包含多少個簇，但很多情況下，我們並不知道資料的分佈情況，實際上聚類就是我們發現數據分佈的一種手段，這就陷入了雞和蛋的矛盾。如何有效的確定K值，這裡大致提供幾種方法，若各位有更好的方法，歡迎探討。

1.與層次聚類結合[2]

         經常會產生較好的聚類結果的一個有趣策略是，首先採用層次凝聚演算法決定結果粗的數目，並找到一個初始聚類，然後用迭代重定位來改進該聚類。

2.穩定性方法[3]

        穩定性方法對一個數據集進行2次重取樣產生2個數據子集，再用相同的聚類演算法對2個數據子集進行聚類，產生2個具有k個聚類的聚類結果，計算2個聚類結果的相似度的分佈情況。2個聚類結果具有高的相似度說明k個聚類反映了穩定的聚類結構，其相似度可以用來估計聚類個數。採用次方法試探多個k，找到合適的k值。

3.系統演化方法[3]

         系統演化方法將一個數據集視為偽熱力學系統，當資料集被劃分為K個聚類時稱系統處於狀態K。系統由初始狀態K=1出發，經過分裂過程和合並過程，系統將演化到它的穩定平衡狀態Ki，其所對應的聚類結構決定了最優類數Ki。系統演化方法能提供關於所有聚類之間的相對邊界距離或可分程度，它適用於明顯分離的聚類結構和輕微重疊的聚類結構。

4.使用canopy演算法進行初始劃分[4]

          基於Canopy Method的聚類演算法將聚類過程分為兩個階段
         Stage1、聚類最耗費計算的地方是計算物件相似性的時候，Canopy Method在第一階段選擇簡單、計算代價較低的方法計算物件相似性，將相似的物件放在一個子集中，這個子集被叫做Canopy ，通過一系列計算得到若干Canopy，Canopy之間可以是重疊的，但不會存在某個物件不屬於任何Canopy的情況，可以把這一階段看做資料預處理；
          Stage2、在各個Canopy 內使用傳統的聚類方法(如K-means)，不屬於同一Canopy 的物件之間不進行相似性計算。
從這個方法起碼可以看出兩點好處：首先，Canopy 不要太大且Canopy 之間重疊的不要太多的話會大大減少後續需要計算相似性的物件的個數；其次，類似於K-means這樣的聚類方法是需要人為指出K的值的，通過Stage1得到的Canopy 個數完全可以作為這個K值，一定程度上減少了選擇K的盲目性。

         其他方法如貝葉斯資訊準則方法（BIC）可參看文獻[5]。

（2）初始質心的選取

          選擇適當的初始質心是基本kmeans演算法的關鍵步驟。常見的方法是隨機的選取初始質心，但是這樣簇的質量常常很差。處理選取初始質心問題的一種常用技術是：多次執行，每次使用一組不同的隨機初始質心，然後選取具有最小SSE（誤差的平方和）的簇集。這種策略簡單，但是效果可能不好，這取決於資料集和尋找的簇的個數。           第二種有效的方法是，取一個樣本，並使用層次聚類技術對它聚類。從層次聚類中提取K個簇，並用這些簇的質心作為初始質心。該方法通常很有效，但僅對下列情況有效：（1）樣本相對較小，例如數百到數千（層次聚類開銷較大）；（2）K相對於樣本大小較小            第三種選擇初始質心的方法，隨機地選擇第一個點，或取所有點的質心作為第一個點。然後，對於每個後繼初始質心，選擇離已經選取過的初始質心最遠的點。使用這種方法，確保了選擇的初始質心不僅是隨機的，而且是散開的。但是，這種方法可能選中離群點。此外，求離當前初始質心集最遠的點開銷也非常大。為了克服這個問題，通常該方法用於點樣本。由於離群點很少（多了就不是離群點了），它們多半不會在隨機樣本中出現。計算量也大幅減少。           第四種方法就是上面提到的canopy演算法。

（3）距離的度量

          常用的距離度量方法包括：歐幾里得距離和餘弦相似度。兩者都是評定個體間差異的大小的。歐幾里得距離度量會受指標不同單位刻度的影響，所以一般需要先進行標準化，同時距離越大，個體間差異越大；空間向量餘弦夾角的相似度度量不會受指標刻度的影響，餘弦值落於區間[-1,1]，值越大，差異越小。但是針對具體應用，什麼情況下使用歐氏距離，什麼情況下使用餘弦相似度？           從幾何意義上來說，n維向量空間的一條線段作為底邊和原點組成的三角形，其頂角大小是不確定的。也就是說對於兩條空間向量，即使兩點距離一定，他們的夾角餘弦值也可以隨意變化。感性的認識，當兩使用者評分趨勢一致時，但是評分值差距很大，餘弦相似度傾向給出更優解。舉個極端的例子，兩使用者只對兩件商品評分，向量分別為(3,3)和(5,5)，這兩位使用者的認知其實是一樣的，但是歐式距離給出的解顯然沒有餘弦值合理。[6]

（4）質心的計算

         對於距離度量不管是採用歐式距離還是採用餘弦相似度，簇的質心都是其均值，即向量各維取平均即可。對於距離對量採用其它方式時，這個還沒研究過。

（5）演算法停止條件

         一般是目標函式達到最優或者達到最大的迭代次數即可終止。對於不同的距離度量，目標函式往往不同。當採用歐式距離時，目標函式一般為最小化物件到其簇質心的距離的平方和，如下：
         當採用餘弦相似度時，目標函式一般為最大化物件到其簇質心的餘弦相似度和，如下：

（6）空聚類的處理

           如果所有的點在指派步驟都未分配到某個簇，就會得到空簇。如果這種情況發生，則需要某種策略來選擇一個替補質心，否則的話，平方誤差將會偏大。一種方法是選擇一個距離當前任何質心最遠的點。這將消除當前對總平方誤差影響最大的點。另一種方法是從具有最大SSE的簇中選擇一個替補的質心。這將分裂簇並降低聚類的總SSE。如果有多個空簇，則該過程重複多次。另外，程式設計實現時，要注意空簇可能導致的程式bug。

3.適用範圍及缺陷

           Kmenas演算法試圖找到使平凡誤差準則函式最小的簇。當潛在的簇形狀是凸面的，簇與簇之間區別較明顯，且簇大小相近時，其聚類結果較理想。前面提到，該演算法時間複雜度為O(tKmn)，與樣本數量線性相關，所以，對於處理大資料集合，該演算法非常高效，且伸縮性較好。但該演算法除了要事先確定簇數K和對初始聚類中心敏感外，經常以區域性最優結束，同時對“噪聲”和孤立點敏感，並且該方法不適於發現非凸面形狀的簇或大小差別很大的簇。

4.實現

#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include <math.h>
#include <stdlib.h>
#define k 3//簇的數目
using namespace std;
//存放元組的屬性資訊
typedef vector<double> Tuple;//儲存每條資料記錄

int dataNum;//資料集中資料記錄數目
int dimNum;//每條記錄的維數

//計算兩個元組間的歐幾裡距離
double getDistXY(const Tuple& t1, const Tuple& t2) 
{
	double sum = 0;
	for(int i=1; i<=dimNum; ++i)
	{
		sum += (t1[i]-t2[i]) * (t1[i]-t2[i]);
	}
	return sqrt(sum);
}

//根據質心，決定當前元組屬於哪個簇
int clusterOfTuple(Tuple means[],const Tuple& tuple){
	double dist=getDistXY(means[0],tuple);
	double tmp;
	int label=0;//標示屬於哪一個簇
	for(int i=1;i<k;i++){
		tmp=getDistXY(means[i],tuple);
		if(tmp<dist) {dist=tmp;label=i;}
	}
	return label;	
}
//獲得給定簇集的平方誤差
double getVar(vector<Tuple> clusters[],Tuple means[]){
	double var = 0;
	for (int i = 0; i < k; i++)
	{
		vector<Tuple> t = clusters[i];
		for (int j = 0; j< t.size(); j++)
		{
			var += getDistXY(t[j],means[i]);
		}
	}
	//cout<<"sum:"<<sum<<endl;
	return var;

}
//獲得當前簇的均值（質心）
Tuple getMeans(const vector<Tuple>& cluster){
	
	int num = cluster.size();
	Tuple t(dimNum+1, 0);
	for (int i = 0; i < num; i++)
	{
		for(int j=1; j<=dimNum; ++j)
		{
			t[j] += cluster[i][j];
		}
	}
	for(int j=1; j<=dimNum; ++j)
		t[j] /= num;
	return t;
	//cout<<"sum:"<<sum<<endl;
}

void print(const vector<Tuple> clusters[])
{
	for(int lable=0; lable<k; lable++)
	{
		cout<<"第"<<lable+1<<"個簇："<<endl;
		vector<Tuple> t = clusters[lable];
		for(int i=0; i<t.size(); i++)
		{
			cout<<i+1<<".(";
			for(int j=0; j<=dimNum; ++j)
			{
				cout<<t[i][j]<<", ";
			}
			cout<<")\n";
		}
	}
}

void KMeans(vector<Tuple>& tuples){
	vector<Tuple> clusters[k];//k個簇
	Tuple means[k];//k箇中心點
	int i=0;
	//一開始隨機選取k條記錄的值作為k個簇的質心（均值）
	srand((unsigned int)time(NULL));
	for(i=0;i<k;){
		int iToSelect = rand()%tuples.size();
		if(means[iToSelect].size() == 0)
		{
			for(int j=0; j<=dimNum; ++j)
			{
				means[i].push_back(tuples[iToSelect][j]);
			}
			++i;
		}
	}
	int lable=0;
	//根據預設的質心給簇賦值
	for(i=0;i!=tuples.size();++i){
		lable=clusterOfTuple(means,tuples[i]);
		clusters[lable].push_back(tuples[i]);
	}
	double oldVar=-1;
	double newVar=getVar(clusters,means);
	cout<<"初始的的整體誤差平方和為："<<newVar<<endl; 
	int t = 0;
	while(abs(newVar - oldVar) >= 1) //當新舊函式值相差不到1即準則函式值不發生明顯變化時，演算法終止
	{
		cout<<"第 "<<++t<<" 次迭代開始："<<endl;
		for (i = 0; i < k; i++) //更新每個簇的中心點
		{
			means[i] = getMeans(clusters[i]);
		}
		oldVar = newVar;
		newVar = getVar(clusters,means); //計算新的準則函式值
		for (i = 0; i < k; i++) //清空每個簇
		{
			clusters[i].clear();
		}
		//根據新的質心獲得新的簇
		for(i=0; i!=tuples.size(); ++i){
			lable=clusterOfTuple(means,tuples[i]);
			clusters[lable].push_back(tuples[i]);
		}
		cout<<"此次迭代之後的整體誤差平方和為："<<newVar<<endl; 
	}

	cout<<"The result is:\n";
	print(clusters);
}
int main(){

	char fname[256];
	cout<<"請輸入存放資料的檔名： ";
	cin>>fname;
	cout<<endl<<" 請依次輸入: 維數 樣本數目"<<endl;
	cout<<endl<<" 維數dimNum: ";
	cin>>dimNum;
	cout<<endl<<" 樣本數目dataNum: ";
	cin>>dataNum;
	ifstream infile(fname);
	if(!infile){
		cout<<"不能開啟輸入的檔案"<<fname<<endl;
		return 0;
	}
	vector<Tuple> tuples;
	//從檔案流中讀入資料
	for(int i=0; i<dataNum && !infile.eof(); ++i)
	{
		string str;
		getline(infile, str);
		istringstream istr(str);
		Tuple tuple(dimNum+1, 0);//第一個位置存放記錄編號，第2到dimNum+1個位置存放實際元素
		tuple[0] = i+1;
		for(int j=1; j<=dimNum; ++j)
		{
			istr>>tuple[j];
		}
		tuples.push_back(tuple);
	}

	cout<<endl<<"開始聚類"<<endl;
	KMeans(tuples);
	return 0;
}

這裡是隨機選取的初始質心，以鳶尾花的資料集為例，原資料集中1-50為一個簇，51-100為第二個簇，101到150為第三個簇：
第一次執行結果 SSE=97.5905

第二次執行結果 SSE=98.1404

。。。

第五次執行結果 SSE=123.397

         由於初始質心是隨機選取的，前兩次還算正常，執行到第五次時，第一個簇基本包括了後51-150個記錄，第二個簇和第三個簇包含了第1-50個記錄，可能的原因就是隨機選擇初始點時，有兩個初始點都選在了1-50個記錄中。

參考：

[1]Pang-Ning Tan等著，《資料探勘導論》，2011

[2]Jiawei Han等著，《資料探勘概念與技術》，2008

[3]聚類分析中類數估計方法的實驗比較

[4]http://www.cnblogs.com/vivounicorn/archive/2011/09/23/2186483.html

[5]一種基於貝葉斯資訊準則的k均值聚類方法

[6]http://www.zhihu.com/question/19640394?nr=1&noti_id=8736954