Science上发表的简单快速的聚类方法
2014-08-12
工作以后发现自己学习和研究的时间变得少得可怜。
前两周因为一个同事的交流，关注了一下canopy辅助Kmeans聚类确定簇数目。然后想起最近很火的一篇Science文章：Clustering by fast search and find of density peak,据说非常简单而优美。然后上网上搜了一下，评论的文章也就那样转来转去，其实就是把人家论文拿来翻译一下，有些关键点根本没讲清楚，真不知道翻译者是不是自己实现过那个算法。
我之所以对这个算法感兴趣，主要是因为看到论文中可以识别那么奇形怪状的点簇，然后又只有两个指标，好像很有道理又很好算的样子。没想到被坑惨了，我用了差不多两个星期，偶尔下班后有时间看论文、写代码，才把这个简单的算法实现下来。其中依然还有一个参数需要手工调整，就是delta的阈值（下面有讲）。
这个算法，其实是对所有坐标点，基于相互距离，提出了两个新的属性，一是局部密度rho，即与该点距离在一定范围内的点的总数，二是到更高密度点的最短距离delta。作者提出，类簇的中心是这样的一类点：它们被很多点围绕（导致局部密度大），且与局部密度比自己大的点之间的距离也很远。
聚类过程，这篇文章讲的还不错，首先我们要把所有点的两个属性算出来，可以画出一个平面决策图；根据决策图，选出rho和delta都大的点作为聚类中心；选定中心后，让周围的点采取“跟随”策略，归类到密度比自己大的最近邻居所在的簇。 我更想说的是两点，第一点关于聚类中心的选择，我翻看了几篇文章（包括作者原文），都没有明确指出一个很好地自动确定中心数目的方法，较多的做法是画出决策图后进行人工选定范围。我的做法是，按rho排序，然后根据决策图人工对delta取一个合适的值，大于这个delta值的，才能被选为中心。这种做法很人肉，让我一直有没完结的感觉，继续探索一下有没有自动确定这个阈值的方法。
第二点是关于rho的计算，其实论文中只提到一个计算公式，是通过截断距离做线性判断，即rho=sigma(sign(dij-dc)),这个计算方法对一般的球状簇，如图1，图2，有不错的效果，而且计算快速，但是对图3的异形图（类簇形状并不呈球状分布），效果就不好了。这时候翻看作者给出的matlab源码，了解到还可以使用高斯核(Gaussian Kernel)函数来定义局部密度，引入以后就完全handle了异形的问题，见图4，我觉得高斯核函数确实是强大（但是为什么呢？）。具体实现可以参看python代码（两种rho的计算都有）。
DATA来源:http://cs.joensuu.fi/sipu/datasets/
[CODE:flame_GaussianKernel.py]
import math
def getDistance(pt1, pt2):
tmp = pow(pt1[0]-pt2[0],2) + pow(pt1[1]-pt2[1],2)
return pow(tmp,0.5)
def ChooseDc(dc_percent,points,dis,distance):
avgNeighbourNum = dc_percent*len(points)
maxd = 0
for i in range(0,len(points)):
for j in range(i+1,len(points)):
pt1 = points[i]
pt2 = points[j]
d = getDistance(pt1,pt2)
dis.append(d)
distance[i,j] = d
dis.append(d)
distance[j,i] = d
if d>maxd:
maxd = d
dis.sort()
return dis[int(avgNeighbourNum*len(points)*2)]
def drawOriginGraph(pl,points,cl,colorNum):
x = [xx for (xx,yy) in points]
y = [yy for (xx,yy) in points]
cm = pl.get_cmap("RdYlGn")
for i in range(len(points)):
pl.plot(x[i],y[i],'o',color=cm(cl[i]*1.0/colorNum))
def drawDecisionGraph(pl,rho, delta,cl,colorNum):
cm = pl.get_cmap("RdYlGn")
for i in range(len(rho)):
pl.plot(rho[i], delta[i],'o',color=cm(cl[i]*1.0/colorNum))
pl.xlabel(r'$\rho$')
pl.ylabel(r'$\delta$')
def Cluster():
#=========Load Data=========
InputFileName = "flame"
OutputFileName = InputFileName + "_out"
suffix = ".txt"
Fin = open(InputFileName+suffix,"r")
Fout = open(OutputFileName+suffix,"w")
points = []
for line in Fin.readlines():
data = line.split()
if len(data)==3:
a = float(data[0])
b = float(data[1])
points.append((a,b))
#=========Calculating=========
#-----choose dc-----
dc_percent = 0.015
dis = []
distance = {}
dc = ChooseDc(dc_percent,points,dis,distance)
print("dc:"+str(dc))
#-----cal rho:"Cut off" kernel
'''
rho = [0 for i in range(len(points))]
for i in range(0,len(points)):
for j in range(i+1,len(points)):
dij = getDistance(points[i],points[j])
if dij<dc:
rho[i] += 1
rho[j] += 1
'''
#-----cal rho:"Gaussian Kernel"
rho = [0 for i in range(len(points))]
for i in range(0,len(points)):
for j in range(i+1,len(points)):
dij = getDistance(points[i],points[j])
rho[i] += math.exp(-(dij/dc)*(dij/dc))
rho[j] += math.exp(-(dij/dc)*(dij/dc))
rho_list =[(rho[i],i) for i in range(len(rho))]
rho_sorted = sorted(rho_list, reverse=1)
print("Highest rho:",rho_sorted[0])
maxd = dis[-1]
delta = [maxd for i in range(len(points))]
nneigh = [-1 for i in range(len(points))]
for ii in range(1,len(rho_sorted)):
for jj in range(0,ii):
id_p1 = rho_sorted[ii][1] #get point1's id
id_p2 = rho_sorted[jj][1] #get point2's id
if (distance[id_p1,id_p2]<delta[id_p1]):
delta[id_p1] = distance[id_p1,id_p2]
nneigh[id_p1] = id_p2
#assignment
cl = [-1 for i in range(len(points))]
colorNum = 0
for ii in range(len(rho_sorted)):
id_p = rho_sorted[ii][1]
if (cl[id_p] == -1 and delta[id_p]>7.0):
cl[id_p] = colorNum
colorNum += 1
else:
if (cl[id_p] == -1 and cl[nneigh[id_p]!=-1]):
cl[id_p] = cl[nneigh[id_p]]
print(colorNum)
import pylab as pl
fig1 = pl.figure(1)
pl.subplot(121)
drawOriginGraph(pl,points,cl,colorNum)
pl.subplot(122)
drawDecisionGraph(pl,rho,delta,cl,colorNum)
pl.show()
for i in range(len(points)):
Fout.write(str(i)+","+str(rho[i])+","+str(delta[i])+"\n")
#Assign Cluster
Fin.close()
Fout.close()
if __name__=="__main__":
Cluster()