層次聚類演算法原理及實現
聚類
聚類是對點集進行考察並按照某種距離測度將他們聚成多個“簇”的過程。聚類的目標是使得同一簇內的點之間的距離較短,而不同簇中點之間的距離較大。
一、聚類演算法介紹
層次法聚類和點分配法聚類。
1.1 點、空間和距離
點集是一種適合於聚類的資料集,每個點都是某空間下的物件。一般意義上,空間只是點的全集,也就是說資料集中的點從該集合中抽樣而成。特別地,歐式空間下的點就是實數向量。向量的長度就是空間的維度數,而向量的分量通常稱為所表示點的座標(coordinate)。
能夠進行聚類的所有空間下都有一個距離測度,即給出空間下任意兩點的距離。一般的歐氏距離(點的座標在各個維度上差值的平方和的算術平方根)可以作為所有歐式空間下的距離測度。
現代聚類問題可能並不這麼簡單。他們可能牽涉非常高維的歐式空間或者根本不在歐式空間下聚類。比如,基於文件間高頻區分詞的共現情況來依據主題對文件聚類。而按照電影愛好者所喜歡的電影類別對他們進行聚類。
1.2 聚類策略
按照聚類演算法使用的兩種不同的基本策略,可以將聚類演算法分成兩類。
1) 層次(hierarchical)或凝聚式(agglomerative)演算法。
這類演算法一開始將每個點都看成簇。簇與簇之間按照接近度(closeness)來組合,接近度可以按照“接近”的不同含義採用不同的定義。當進一步的組合導致多個原因之下的非期望結果時,上述組合過程結束。比如停止條件為:達到預先給定的簇數目
2) 點分配過程演算法。按照某個順序依次考慮每個點,並將它分配到最適合的簇中。該過程通常有一個短暫的初始簇估計階段。一些變形演算法允許臨時的簇合並或分裂的過程,或者當點為離群點時允許不將該點分配到任何簇中。
聚類演算法也可以使用如下方式分類:
1) 歐式空間下,我們可以將點集合概括為其質心,即點的平均。而在非歐式空間下根本沒有質心的概念,因此需要其他的簇概括方法。
2) 演算法是否假設資料足夠小的能夠放入記憶體?或者說資料是否必須主要存放在二次儲存器?由於不能將所有簇的所有點都同時放入記憶體,所以我們將簇的概括表示存放在記憶體中也是必要的。
1.3 維數災難
“災難”的一個體現是,在高維空間下,幾乎所有點對之間的距離都差不多相等。另一個表現是,幾乎任意的兩個向量之間都近似正交。
1. 高維空間下的距離分佈
一個d維歐式空間,假設在一個單位立方體內隨機選擇n個點,每個點可以表示成[x1,x2,…,xd],其每個xi都是0到1之間。假定d非常大,兩個隨機點[x1,x2,…,xd]和[y1,y2,…,yd]之間的歐式距離為
上述基於隨機資料的論證結果表明,在這麼多所有距離近似相等的點對之中發現聚類是很難的一件事。
2. 向量之間的夾角
在d維空間的隨機點A,B,C,其中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
- #include <iostream>
- #include <vector>
- #include <algorithm>
- #include <fstream>
- usingnamespace std;
- constint iniClusNum = 12;
- constint 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++;
- }
- elseif(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;
- }
- elseif(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;
- }
- elseif(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");
- }
#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;輸出如下,採用的是輸入資料的索引