python机器学习——kmeans聚类算法
背景与原理:
聚类问题与分类问题有一定的区别,分类问题是对每个训练数据,我给定了类别的标签,现在想要训练一个模型使得对于测试数据能输出正确的类别标签,更多见于监督学习;而聚类问题则是我们给出了一组数据,我们并没有预先的标签,而是由机器考察这些数据之间的相似性,将相似的数据聚为一类,是无监督学习的一个典型应用。
而k-means算法则是非常常见的聚类算法,其思想是如果我们想把这些数据聚为k类,那么我们预先选择k个中心,然后计算每个数据点与这k个中心之间的“距离”(也就是这个数据点与这个中心的“相似度”),那么非常自然地,每个数据点应当被划分进离他距离最近的那个中心点对应的类。
但是这是最优的聚类方法吗?如果我们初始选的k个点很糟糕,其实有些数据点离这k个点都很远,而这种划分只是一种“矮子里面拔将军”的划分,因此可能会把两个差异巨大的点划分到一个聚类里面去,因此我们需要迭代上述过程,即当我们选中了一些数据点聚为一类之后,我们取这些数据点的“质心”作为新的中心点,这样我们会得到k个新的中心点,然后我们重复上述过程,直到中心点不再移动为止。
那么我们就要解决几个问题:
一.我们如何度量“距离”或“相似度”?
距离的度量其实有很多方式,对于两个数据点$(x_{1},x_{2},...,x_{n}),(y_{1},y_{2},...,y_{n})$,常见的距离度量有如下的方式:
欧氏距离:$d(x,y)=\sqrt{\sum_{i=1}^{n}(x_{i}-y_{i})^{2}}$
曼哈顿距离:$d(x,y)=\sum_{i=1}^{n}|x_{i}-y_{i}|$
切比雪夫距离:$d(x,y)=max_{i=1}^{n}|x_{i}-y_{i}|$
余弦距离:$d(x,y)=\cos \theta =\dfrac{\sum_{i=1}^{n}x_{i}y_{i}}{\sqrt{\sum_{i=1}^{n}x_{i}^{2}} \sqrt{\sum_{i=1}^{n}y_{i}^{2}}}$
相关系数:$\rho_{XY}=\dfrac{Cov(X,Y)}{\sqrt{D(X)} \sqrt{D(Y)}}=\dfrac{E((X-EX)(Y-EY))}{\sqrt{D(X)} \sqrt{D(Y)}}$
在这里我们选用欧氏距离来度量点之间的距离,选取SSE(误差平方和)作为损失函数,即设我们有$k$类,第$i$类的中心点为$c_{i}$,那么损失函数为:
$J(c)=\sum_{i=1}^{k}\sum_{x\in C_{i}}d(x,c_{i})^{2}$
那么此时我们每次迭代过程中选取的新中心点是什么呢?
我们对第$i$个中心点$c_{i}$求偏导:
$\dfrac{\partial J(c)}{\partial c_{i}}=\dfrac{\partial \sum_{x\in C_{i}}|x-c_{i}|^{2}}{\partial c_{i}}=-2\sum_{x\in C_{i}} (x-c_{i})$
而取得极小值时,上述偏导数为零,即:
$c_{i}=\dfrac{\sum_{x\in C_{i}}x}{|C_{i}|}$
也即新的中心点应该是聚在这一类里的所有点的算术平均值
设计、调整与评价:
如何选取初始聚类中心?
(1)凭经验直接选取
(2)将数据随机分成k类,计算每类中心作为初始聚类中心
(3)求以每个数据点为球心,某个半径内的特征点个数,选取密度最大的特征点为第一个聚类中心,然后在离这个聚类中心距离大于某个距离$d$的特征点中选取另一个密度最大的特征点,以此类推直至选出k个点
(4)用距离最远的k个点作为初始中心
(5)n较大时,先随机选出一部分聚成k类,再将这k个中心作为初始聚类中心
如何选取聚类个数?
(1)按聚类目标(比如手写数字识别)等先验知识确定k
(2)让k从小到大增加,那么损失函数显然在减少,选择损失函数下降的拐点对应的k
如何评价聚类效果?
一般我们可以用几个指标来评价,比如$P$值(纯度)和$F$值
所谓纯度,是指如果我们已知每个数据点的类别,我们不妨假设一共有$k$类,那么对于聚出的第$r$类,其纯度$P(S_{r})=\dfrac{max_{i=1}^{k}n_{ri}}{n_{r}}$,所谓$n_{ri}$,就是被聚在第$r$类的数据中心原本属于第$i$类的数据点个数,而$n_{r}$就是整个被聚出的第$r$类的元素个数,而我们整个聚类过程的总纯度为:
$P=\sum_{r=1}^{k}\dfrac{n_{r}}{n}P(S_{r})$
其中$n$为所有数据点的总数。
而所谓$F$值,是准确率(precision)和召回率(recall)的调和平均值。
这里要解决的其实是一个问题:我们已知原数据有$k$类,而我们聚出了$k$个类,那...我们怎么把聚出的$k$个类和预先想分出的$k$类对应起来?
举个例子:手写数字识别,我们要识别手写的0~9十个数字,而我们聚出了十个类,那我们要怎么知道每个类对应的是哪个数字呢?
那么我们定义$precision(i,r)=\dfrac{n_{ri}}{n_{r}}$,即如果我们认为聚出的第$r$类对应于原数据中的第$i$类,那么其准确率即为这个类里确实属于第$i$类的数据占比
而$recall(i,r)=\dfrac{n_{ri}}{n_{i}}$,其中$n_{i}$表示在原标签中属于第$i$类的数据个数,即如果我们认为聚出的第$r$类对应于原数据中的第$i$类,那么召回率即为这个类里确实属于第$i$类里的数据占所有第$i$类数据的占比
那么我们定义$f(i,r)=\dfrac{2*precision(i,r)*recall(i,r)}{precision(i,r)+recall(i,r)}$,而整个数据集的$F$定义为:
$F=\sum_{i=1}^{k}\dfrac{n_{i}}{n}f(i,r)$
当然,这个值取决于如何对聚类和原始类别对应,因此我们想最大化这个值,我们就要使用一种二分图匹配算法,常用的是KM算法。
KM算法可以查看别的介绍,这里简要介绍下其功能:我们把聚出的k类和原始分好的k类分在两侧,那么这可以看做一个二分图模型,而我们构造的聚类与分好类的对应关系就是一种二分图的匹配,而这个匹配过程的要求是我们要最大化$F$值,那么如果我们设聚出的第$r$类和原始的第$i$类之间的边权为$n_{i}f(i,r)$,我们进行的就是二分图最佳匹配,而这个匹配可以用KM算法计算出来。
代码实现:
from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt plt.figure() img = plt.imread('./cat.jpeg') plt.imshow(img) def kmeans_iteration(l): oril=[] for i in l: oril.append(i) flag=0 for i in dic: p=0 mind=10000000 for j in range(0,16): d=pow((i[0]-l[j][0]),2)+pow((i[1]-l[j][1]),2)+pow((i[2]-l[j][2]),2) if d<mind: p=j mind=d if dic[i]!=p: flag=1 dic[i]=p if flag==0: return l else:for i in range(0,16): cnt=0 r=0 b=0 g=0 for j in dic: if dic[j]==i: r+=j[0] b+=j[1] g+=j[2] cnt+=1 r/=cnt b/=cnt g/=cnt l[i]=(r,b,g) for i in range(0,16): d=pow((oril[i][0]-l[i][0]),2)+pow((oril[i][1]-l[i][1]),2)+pow((oril[i][2]-l[i][2]),2) if d>1: flag=0 if flag==1: return l else: return kmeans_iteration(l) templ=[] for i in range(1,17): templ.append((i*10,i*10,i*10)) retl=kmeans_iteration(templ) re=np.zeros((1080,1080,3)) for i in range(0,1080): for j in range(0,1080): p=0 mind=1000000 for k in range(0,16): d=pow((img[i][j][0]-retl[k][0]),2)+pow((img[i][j][1]-retl[k][1]),2)+pow((img[i][j][2]-retl[k][2]),2) if d<mind: p=k mind=d for k in range(0,3): re[i][j][k]=retl[p][k] plt.imshow(re/255)
这段代码是kmeans的一个手写实现,实现了将一张图片压缩至16色,原理是将所有颜色点聚成16类,每个颜色用其聚类中心取代
这是原始图片:
这是处理后的图片:
可以看到压缩效果还是很不错的
from sklearn import datasets, preprocessing from sklearn.decomposition import PCA from sklearn.cluster import KMeans, MeanShift import pandas as pd import numpy as np from matplotlib import pyplot as plt %matplotlib inline X = pd.read_csv('./train_X.csv') # 为了方便起见,这里只采用前6000个MNIST数据 y = pd.read_csv('./train_y.csv') X, y = np.array(X), np.array(y) print(X.shape) print(y.shape) pca2d = PCA(n_components=2) X_std = preprocessing.scale(X) # 数据标准化 X_2d = pca2d.fit_transform(X_std)# 数据降维至两维便于可视化 plt.scatter(X_2d[:, 0], X_2d[:, 1], c=y) y_pred_std = KMeans(n_clusters=10, random_state=9).fit_predict(X_std) plt.scatter(X_2d[:,0],X_2d[:,1],c=y_pred_std) l2=[] for i in range(0,10): l2.append(dict()) for j in range(0,10): l2[i][j]=0 for i in range(0,6000): l2[y_pred_std[i]][y[i][0]]+=1 P2=0 for i in range(0,10): p=0 for j in range(1,10): if l2[i][j]>l2[i][p]: p=j P2+=l2[i][p]/6000 print(P2)
这个代码展示了使用sklearn里面的KMeans包直接进行kmeans聚类,而同样还进行了一个PCA降维,这个降维的过程主要是用来可视化,即把聚类结果花在一张图上,同时这里还计算了P值
而如果实现类别的对应,我们可以这样写:
dic1=dict()#ni dic3=dict()#nr graph2=[] for i in range(0,10): dic1[i]=0 dic3[i]=0 graph2.append([]) for i in range(0,6000): dic1[y[i][0]]+=1 dic3[y_pred_std[i]]+=1 for i in range(0,10): for j in range(0,10): graph2[i].append(2*l2[j][i]*dic1[i]/(dic1[i]+dic3[j])) def find_path(graph,i, visited_left, visited_right, slack_right): visited_left[i] = True for j, match_weight in enumerate(graph[i]): if visited_right[j]: continue gap = label_left[i] + label_right[j] - match_weight if abs(gap)<1e-3 : visited_right[j] = True if j not in T or find_path(graph,T[j], visited_left, visited_right, slack_right): T[j] = i S[i] = j return True else: slack_right[j] = min(slack_right[j], gap) return False def KM(graph): m = len(graph) for i in range(m): # 重置辅助变量 slack_right = [float('inf') for _ in range(m)] while True: visited_left = [False for _ in graph] visited_right = [False for _ in graph] if find_path(graph,i,visited_left,visited_right, slack_right): break d = float('inf') for j, slack in enumerate(slack_right): if not visited_right[j] and slack < d: d = slack for k in range(m): if visited_left[k]: label_left[k] -= d if visited_right[k]: label_right[k] += d return S, T label_left, label_right = [max(g) for g in graph2], [0 for _ in graph2] S, T = {}, {} visited_left = [False for _ in graph2] visited_right = [False for _ in graph2] slack_right = [float('inf') for _ in graph2] KM(graph2) ans=0 for i in S: ans+=graph2[i][S[i]] print(ans/6000)
这段代码用KM算法计算了上述手写数字识别kmeans的F值。
小结与优化:
kmeans算法的优点是显著的:算法简单易于实现;如果类密集且类与类之间区别明显时聚类效果很好;算法复杂度为$O(Nkt)$,对大数据集而言相对高效
但是其缺点也是显著的:结果与初始质心的选取有关(如果初始质心选取的不好迭代次数会很多,同时效果可能会比较差);必须预先给出要聚类的个数k作为超参数(因此需要调参);对噪声和孤立数据点敏感,少量这样的数据就会对平均值产生较大的影响;不适合发现非凸形状的聚类;在大数据集上收敛比较慢;可能达到局部极小值。
常见的改进有kmeans++算法,即初始选取聚类中心时要求聚类中心离得越远越好;ISODATA算法(迭代自组织数据分析法):可以调整的kmeans方法,当属于某个类的数据点过少就删掉这个类,而当属于某个类的数据过多、分散程度较大时将这个类分成两个子类