1. 程式人生 > >層次聚類演算法原理及實現

層次聚類演算法原理及實現

聚類

聚類是對點集進行考察並按照某種距離測度將他們聚成多個“簇”的過程。聚類的目標是使得同一簇內的點之間的距離較短,而不同簇中點之間的距離較大。

一、聚類演算法介紹

     層次法聚類和點分配法聚類。

1.1     點、空間和距離

點集是一種適合於聚類的資料集,每個點都是某空間下的物件。一般意義上,空間只是點的全集,也就是說資料集中的點從該集合中抽樣而成。特別地,歐式空間下的點就是實數向量。向量的長度就是空間的維度數,而向量的分量通常稱為所表示點的座標(coordinate)

能夠進行聚類的所有空間下都有一個距離測度,即給出空間下任意兩點的距離。一般的歐氏距離(點的座標在各個維度上差值的平方和的算術平方根)可以作為所有歐式空間下的距離測度。

現代聚類問題可能並不這麼簡單。他們可能牽涉非常高維的歐式空間或者根本不在歐式空間下聚類。比如,基於文件間高頻區分詞的共現情況來依據主題對文件聚類。而按照電影愛好者所喜歡的電影類別對他們進行聚類。

1.2     聚類策略

按照聚類演算法使用的兩種不同的基本策略,可以將聚類演算法分成兩類。

1)   層次(hierarchical)或凝聚式(agglomerative)演算法。

這類演算法一開始將每個點都看成簇。簇與簇之間按照接近度(closeness)來組合,接近度可以按照“接近”的不同含義採用不同的定義。當進一步的組合導致多個原因之下的非期望結果時,上述組合過程結束。比如停止條件為:達到預先給定的簇數目

,或者使用簇的緊密度測度方法,一旦兩個小簇組合之後得到簇內的點分散的區域較大就停止簇的構建。

2)   點分配過程演算法。按照某個順序依次考慮每個點,並將它分配到最適合的簇中。該過程通常有一個短暫的初始簇估計階段。一些變形演算法允許臨時的簇合並或分裂的過程,或者當點為離群點時允許不將該點分配到任何簇中。

聚類演算法也可以使用如下方式分類:

1)   歐式空間下,我們可以將點集合概括為其質心,即點的平均。而在非歐式空間下根本沒有質心的概念,因此需要其他的簇概括方法

2)   演算法是否假設資料足夠小的能夠放入記憶體?或者說資料是否必須主要存放在二次儲存器?由於不能將所有簇的所有點都同時放入記憶體,所以我們將簇的概括表示存放在記憶體中也是必要的。

1.3     維數災難

“災難”的一個體現是,在高維空間下,幾乎所有點對之間的距離都差不多相等。另一個表現是,幾乎任意的兩個向量之間都近似正交。

1.        高維空間下的距離分佈

一個d維歐式空間,假設在一個單位立方體內隨機選擇n個點,每個點可以表示成[x1,x2,…,xd],其每個xi都是01之間。假定d非常大,兩個隨機點[x1,x2,…,xd][y1,y2,…,yd]之間的歐式距離為

                                                                                                                    

 上述基於隨機資料的論證結果表明,在這麼多所有距離近似相等的點對之中發現聚類是很難的一件事。

2.        向量之間的夾角

d維空間的隨機點ABC,其中d很大。AC可以在任意位置,而B處於座標原點。那麼夾角ABC的餘弦為:

                                                                                                 

d不斷增長時,分母會隨d線性增長,但是分子卻是隨機值之和,可能為正也可能為負。分子期望為0,分子最大值為 。所以對於很大的d而言,任意兩個向量的夾角餘弦值幾乎肯定接近為0,即夾角近似度等於90度。

推論為:如果dAB = d1, dBC=d2,dAC≈ 。

二、層次聚類

首先考慮歐式空間下的層次聚類。該演算法僅可用於規模相對較小的資料集。層次聚類用於非歐式空間時,還有一些與層次聚類相關的額外問題需要考慮。因此,當不存在簇質心或者說簇平均點時,可以考慮採用簇中心點(clustroid)來表示一個簇。

2.1     歐式空間下的層次聚類

首先,每個點看作一個簇,通過不斷的合併小簇而形成大簇。我們需要提前確定

(1)         簇如何表示?

(2)         如何選擇哪兩個簇進行合併?

(3)         簇合並何時結束?

對於歐式空間,(1)通過簇質心或者簇內平均點來表示簇。對於單點的簇,該點就是簇質心。可以初始化簇數目為歐式空間點的數目Cnumber=n。簇之間的距離為質心之間的歐式距離,

2)選擇具有最短距離(或者其他方式)的兩個簇進行合併。

例如,有如下12個點,首先我們將每一個點看做一個簇。

                                                           

point.txt檔案

4 10
4 8
7 10
6 8
3 4
2 2
5 2
9 3
10 5
11 4
12 3
12 6

[cpp] view plain copy print?
  1. #include <iostream>
  2. #include <vector>
  3. #include <algorithm>
  4. #include <fstream>
  5. usingnamespace std;  
  6. constint iniClusNum = 12;  
  7. constint stopNum = 3;  
  8. class Point  
  9. {  
  10. public:  
  11.     double x;  
  12.     double y;  
  13.     int NumPBelong;  
  14.     Point()  
  15.     {  
  16.         x=0;  
  17.         y=0;  
  18.         NumPBelong = -1;  
  19.     }  
  20.     Point(double x1, double y1, int f=-1):x(x1),y(y1),NumPBelong(f){}  
  21.     const Point& operator=(const Point& p)  
  22.     {  
  23.         x = p.x;  
  24.         y=p.y;  
  25.         NumPBelong = p.NumPBelong;  
  26.         return *this;  
  27.     }  
  28. };  
  29. class ManagerP  
  30. {  
  31. public:  
  32.     double getDistance(const Point& p1, const Point& p2)  
  33.     {  
  34.         return sqrt(pow((p1.x-p2.x),2)+pow((p1.y-p2.y),2));  
  35.     }  
  36.     Point getMean(const Point& p1, const Point& p2)  
  37.     {  
  38.         Point p;  
  39.         p.x = (p1.x+p2.x)/2;  
  40.         p.y=(p1.y+p2.y)/2;  
  41.         return p;  
  42.     }  
  43. };  
  44. class ManagerC  
  45. {  
  46. public:  
  47.     Point Cluster[iniClusNum];  
  48.     vector<int> ClusterLast[iniClusNum];  
  49.     bool isIndexClose[iniClusNum];  
  50.     bool isIndexClose2[iniClusNum];  
  51.     void initCluster()//use txt to init, import txt
  52.     {  
  53.         ifstream  myfile ( "point.txt" ) ;  
  54.         if  ( !myfile )   
  55.         {   
  56.             cout << "cannot open file." ;   return  ;   
  57.         }  
  58.         Point p;  
  59.         int x,y;    
  60.         int i=0;  
  61.         while(!myfile.eof())  
  62.         {  
  63.             myfile>>x>>y;  
  64.             p.x=x;  
  65.             p.y=y;  
  66.             Cluster[i]=p;  
  67.             i++;  
  68.         }  
  69.         myfile.close();  
  70.     }  
  71.     void initIndexClose()  
  72.     {  
  73.             for(int i=0;i<iniClusNum;i++)  
  74.             {  
  75.                 isIndexClose[i]=false;  
  76.                 isIndexClose2[i]=false;  
  77.             }  
  78.     }  
  79.     void print()  
  80.     {  
  81.         for(int i =0;i<iniClusNum;i++)  
  82.         {  
  83.             if(ClusterLast[i].empty())  
  84.             {  
  85.                 continue;  
  86.             }  
  87.             cout<<"cluster "<<i+1<<endl;  
  88.             vector<int>::iterator ite=ClusterLast[i].begin();  
  89.             for(;ite!= ClusterLast[i].end();ite++)  
  90.             {  
  91.                 cout<<*ite<<"\t";  
  92.             }  
  93.             cout<<endl;  
  94.         }  
  95.         cout<<endl;  
  96.     }  
  97.         void ClusterAlgo()//use minheap to realize, to optimize
  98.     {  
  99.         int ClustNum = iniClusNum;  
  100.         int clus_index =0;  
  101.         while(ClustNum>stopNum)  
  102.         {  
  103.             double min=INT_MAX;  
  104.             int x=-1,y=-1;  
  105.             ManagerP mp;  
  106.             for(int i=0;i<iniClusNum;i++)  
  107.             {  
  108.                 if(isIndexClose[i])  
  109.                 {  
  110.                     continue;  
  111.                 }  
  112.                 for(int j=i+1;j<iniClusNum;j++)  
  113.                 {  
  114.                     if(isIndexClose[j])  
  115.                     {  
  116.                         continue;  
  117.                     }  
  118.                     double new_d = mp.getDistance(Cluster[i],Cluster[j]);  
  119.                     if(new_d < min)  
  120.                     {  
  121.                         min = new_d;  
  122.                         x=i;y=j;  
  123.                     }  
  124.                 }  
  125.             }  
  126.             if(x==-1 || y==-1)  
  127.             {  
  128.                 break;  
  129.             }  
  130.             Point p = mp.getMean(Cluster[x],Cluster[y]);  
  131.             //x<y    store the result
  132.             if(Cluster[x].NumPBelong==-1 && Cluster[y].NumPBelong==-1)  
  133.             {  
  134.                 cout<<"a0"<<endl;  
  135.                 ClusterLast[clus_index].push_back(x);//xchange to p, y close
  136.                 ClusterLast[clus_index].push_back(y);  
  137.                 p.NumPBelong = clus_index;  
  138.                 isIndexClose[y]=true;//y is closed
  139.                 Cluster[x]=p;//new p is open
  140.                 isIndexClose[x]=false;  
  141.                 isIndexClose2[x]=true;  
  142.                 isIndexClose2[y]=true;  
  143.                 clus_index++;  
  144.             }  
  145.             elseif(Cluster[x].NumPBelong==-1 && Cluster[y].NumPBelong!=-1)//already exists one cluster
  146.             {  
  147.                 cout<<"a1"<<endl;  
  148.                 ClusterLast[Cluster[y].NumPBelong].push_back(x);  
  149.                 isIndexClose[x]=true;//x is closed
  150.                 p.NumPBelong = Cluster[y].NumPBelong;  
  151.                 Cluster[y]=p;//new p is open
  152.                 isIndexClose2[x]=true;  
  153.             }  
  154.             elseif(Cluster[x].NumPBelong!=-1 && Cluster[y].NumPBelong==-1)  
  155.             {  
  156.                 cout<<"a2"<<endl;  
  157.                 ClusterLast[Cluster[x].NumPBelong].push_back(y);  
  158.                 isIndexClose[y]=true;//y is closed
  159.                 p.NumPBelong = Cluster[x].NumPBelong;  
  160.                 Cluster[x]=p;//new p is open
  161.                 isIndexClose2[y]=true;  
  162.             }  
  163.             elseif(Cluster[x].NumPBelong!=-1 && Cluster[y].NumPBelong!=-1)//both are clusteroid
  164.             {  
  165.                 cout<<"a3"<<endl;  
  166.                 vector<int>::iterator ite = ClusterLast[Cluster[y].NumPBelong].begin();//put y's node in x
  167.                 for(;ite!=ClusterLast[Cluster[y].NumPBelong].end();ite++)  
  168.                 {  
  169.                     ClusterLast[Cluster[x].NumPBelong].push_back(*ite);  
  170.                 }  
  171.                 ClusterLast[Cluster[y].NumPBelong].clear();  
  172.                 isIndexClose[y]=true;//y is closed
  173.                 p.NumPBelong = Cluster[x].NumPBelong;  
  174.                 Cluster[x]=p;//new p is open
  175.             }  
  176.             ClustNum--;  
  177.         }  
  178.         int total_size =0;  
  179.         for(int i=0;i<stopNum;i++)  
  180.         {  
  181.             total_size+=ClusterLast[i].size();  
  182.         }  
  183.         if(total_size<iniClusNum)  
  184.         {  
  185.             int j=0;  
  186.             for(int i=0;i<iniClusNum;i++)  
  187.             {  
  188.                 if(isIndexClose2[i]==false)  
  189.                 {  
  190.                     ClusterLast[stopNum-1-j].push_back(i);  
  191.                     j++;  
  192.                 }  
  193.             }  
  194.         }  
  195.     }  
  196. };  
  197. int main()  
  198. {  
  199.     ManagerC M;  
  200.     M.initCluster();  
  201.     M.initIndexClose();  
  202.     M.ClusterAlgo();  
  203.     M.print();  
  204.     system("pause");  
  205. }  
#include <iostream>
#include <vector>
#include <algorithm>
#include <fstream>
using namespace std;
const int iniClusNum = 12;
const int stopNum = 3;

class Point
{
public:
	double x;
	double y;
	int NumPBelong;
	Point()
	{
		x=0;
		y=0;
		NumPBelong = -1;
	}
	Point(double x1, double y1, int f=-1):x(x1),y(y1),NumPBelong(f){}
	const Point& operator=(const Point& p)
	{
		x = p.x;
		y=p.y;
		NumPBelong = p.NumPBelong;
		return *this;
	}
};

class ManagerP
{
public:
	double getDistance(const Point& p1, const Point& p2)
	{
		return sqrt(pow((p1.x-p2.x),2)+pow((p1.y-p2.y),2));
	}
	Point getMean(const Point& p1, const Point& p2)
	{
		Point p;
		p.x = (p1.x+p2.x)/2;
		p.y=(p1.y+p2.y)/2;
		return p;
	}
};
class ManagerC
{
public:
	Point Cluster[iniClusNum];
	vector<int> ClusterLast[iniClusNum];
	bool isIndexClose[iniClusNum];
	bool isIndexClose2[iniClusNum];
	void initCluster()//use txt to init, import txt
	{
		ifstream  myfile ( "point.txt" ) ;
		if  ( !myfile ) 
		{ 
			cout << "cannot open file." ;   return  ; 
		}

		Point p;
		int x,y;  
		int i=0;
		while(!myfile.eof())
		{
			myfile>>x>>y;
			p.x=x;
			p.y=y;
			Cluster[i]=p;
			i++;
		}
		myfile.close();
	}
	void initIndexClose()
	{
			for(int i=0;i<iniClusNum;i++)
			{
				isIndexClose[i]=false;
				isIndexClose2[i]=false;
			}
	}
	void print()
	{
		for(int i =0;i<iniClusNum;i++)
		{
			if(ClusterLast[i].empty())
			{
				continue;
			}
			cout<<"cluster "<<i+1<<endl;
			vector<int>::iterator ite=ClusterLast[i].begin();
			for(;ite!= ClusterLast[i].end();ite++)
			{
				cout<<*ite<<"\t";
			}
			cout<<endl;

		}
		cout<<endl;
	}
		void ClusterAlgo()//use minheap to realize, to optimize
	{

		int ClustNum = iniClusNum;
		int clus_index =0;
		while(ClustNum>stopNum)
		{

			double min=INT_MAX;
			int x=-1,y=-1;
			ManagerP mp;
			for(int i=0;i<iniClusNum;i++)
			{
				if(isIndexClose[i])
				{
					continue;
				}
				for(int j=i+1;j<iniClusNum;j++)
				{
					if(isIndexClose[j])
					{
						continue;
					}

					double new_d = mp.getDistance(Cluster[i],Cluster[j]);
					if(new_d < min)
					{
						min = new_d;
						x=i;y=j;

					}
				}
			}
			if(x==-1 || y==-1)
			{
				break;
			}

			Point p = mp.getMean(Cluster[x],Cluster[y]);
			//x<y	store the result
			if(Cluster[x].NumPBelong==-1 && Cluster[y].NumPBelong==-1)
			{
				cout<<"a0"<<endl;
				ClusterLast[clus_index].push_back(x);//xchange to p, y close
				ClusterLast[clus_index].push_back(y);
				p.NumPBelong = clus_index;
				isIndexClose[y]=true;//y is closed
				Cluster[x]=p;//new p is open
				isIndexClose[x]=false;
				isIndexClose2[x]=true;
				isIndexClose2[y]=true;
				clus_index++;

			}
			else if(Cluster[x].NumPBelong==-1 && Cluster[y].NumPBelong!=-1)//already exists one cluster
			{
				cout<<"a1"<<endl;
				ClusterLast[Cluster[y].NumPBelong].push_back(x);
				isIndexClose[x]=true;//x is closed
				p.NumPBelong = Cluster[y].NumPBelong;
				Cluster[y]=p;//new p is open
				isIndexClose2[x]=true;
			}
			else if(Cluster[x].NumPBelong!=-1 && Cluster[y].NumPBelong==-1)
			{
				cout<<"a2"<<endl;
				ClusterLast[Cluster[x].NumPBelong].push_back(y);
				isIndexClose[y]=true;//y is closed
				p.NumPBelong = Cluster[x].NumPBelong;
				Cluster[x]=p;//new p is open
				isIndexClose2[y]=true;
			}
			else if(Cluster[x].NumPBelong!=-1 && Cluster[y].NumPBelong!=-1)//both are clusteroid
			{
				cout<<"a3"<<endl;
				vector<int>::iterator ite = ClusterLast[Cluster[y].NumPBelong].begin();//put y's node in x
				for(;ite!=ClusterLast[Cluster[y].NumPBelong].end();ite++)
				{
					ClusterLast[Cluster[x].NumPBelong].push_back(*ite);
				}
				ClusterLast[Cluster[y].NumPBelong].clear();
				isIndexClose[y]=true;//y is closed
				p.NumPBelong = Cluster[x].NumPBelong;
				Cluster[x]=p;//new p is open
							
			}
			ClustNum--;
		}
		int total_size =0;
		for(int i=0;i<stopNum;i++)
		{
			total_size+=ClusterLast[i].size();
		}
		if(total_size<iniClusNum)
		{
			int j=0;
			for(int i=0;i<iniClusNum;i++)
			{
				if(isIndexClose2[i]==false)
				{
					ClusterLast[stopNum-1-j].push_back(i);
				    j++;
				}
			}

		}
	}

};
int main()
{
	ManagerC M;
	M.initCluster();
	M.initIndexClose();
	M.ClusterAlgo();
	M.print();

	system("pause");
}
  如果仔細觀察資料的資料在座標系的分佈,可以分成3簇。於是我們使StopNum =3;輸出如下,採用的是輸入資料的索引