kmeans算法性能改进_kmeans++算法+kmeans++优化算法+距离计算优化

***题记&＃xff1a;*我一直在路上&＃xff0c;害怕停下

在我的另一篇博客里《读sklearn源码学机器学习——kmeans聚类算法》我详细的阐述了kmeans算法的工作过程。截至目前为止&＃xff0c;还没有深入的刨析kmeans算法的工作原理(会用和知道怎么用&＃xff0c;跟理解背后深刻的数学原理是有本质区别的&＃xff0c;我对此深感敬畏)。其实kmeas算法和高斯混合算法都是em算法的具体应用。今天站在工程应用的角度&＃xff0c;刨析kmeans在工程中的应用和优化。主要是弥补我上篇博客中没有说明白的两个函数&＃xff08;elativeDist,squaredNorna&＃xff09;和一个初始化。

1、kmeans&＃43;&＃43;算法

kmeans&＃43;&＃43;是David Arthur and Sergei Vassilvitskii在一篇名为“k-means&＃43;&＃43;: The Advantages of Careful Seeding&＃xff08;http://ilpubs.stanford.edu:8090/778/1/2006-13.pdf&＃xff09;文中中首次提及的&＃xff0c;可以用来解决kmeans在初始化的时候&＃xff0c;随机选取初始点带来的负面影响。算法的工作过程如下&＃xff1a;

翻译成中文&＃xff1a;
1a、随机的选取一个样本点作为一个一个聚类中心。
2b、在选择下一个聚类中心的时候&＃xff0c;以到达现有聚类中心的距离为依据进行聚类中心得选择。离现在的聚类中心越远被选择的概率越大。
&＃xff08;一般通过轮盘赌选择&＃xff0c;概率一般用欧式距离平方占比&＃xff09;
1c、重复上述的过程直到选出所有的聚类中心。
基本依据&＃xff1a;在选择初始距离的时候&＃xff0c;样本之间离得越远越好。
kmeans&＃43;&＃43;实现过程&＃xff1a;
数据通过矩阵的形式存在matrix_data中&＃xff08;个人觉得自己的实现过程还是优雅且美妙的&＃xff09;。
&＃xff08;可以仔细体会下对轮盘赌的实现过程以及numpy相应函数的应用&＃xff0c;可以猜测np.searchsorted()的算法复杂度最高时log(n),比我们遍历选取要划算多了&＃xff09;

center&＃61;np.zeros((n_cluster,n_feature))
center[0]&＃61;matrix_data[np.random.randint(n_samples)]
for i in range(1,n_cluster):#计算每个样本点到已有中心点的距离distance_to_centers&＃61;euclideanDistance(matrix_data,center[[i for i in range(i)]],square&＃61;True)#选取据距离最小值closed_distance&＃61;np.min(distance_to_centers,axis&＃61;1)#轮盘赌denominator&＃61;closed_distance.sum()point&＃61;np.random.rand()*denominator#轮盘赌的指针be_choosed&＃61;np.searchsorted(np.cumsum(closed_distance),point)be_choosed&＃61;min(be_choosed,n_samples-1)#避免选的是最后一个而造成下标越界center[i]&＃61;matrix_data[be_choosed]


其中欧式距离的计算函数euclideanDistance实现原理在上一篇博客中证明过&＃xff0c;过程如下&＃xff1a;

import pandas as pd
import numpy as np
def rowNorms(X):#对行每个元素取平方加和return np.einsum("ij,ij->i",X,X)
def euclideanDistance(x,y,square&＃61;False):#x的每个样本与y之间的距离""":param x: 矩阵x:param y: 矩阵y:param squared: 表示是否返回二者欧式距离的平方值,很明显如果返回欧式距离的平方&＃xff0c;计算量又小一些:return: 矩阵x中的每一个样本与y中样本之间的距离""""""对于这个操作的理解一定要理解下面这个操作np.array([[1],[2],[1]])&＃43;np.array([[1,2,3],[4,5,6],[7,8,9]])Out[28]:array([[ 2, 3, 4],[ 6, 7, 8],[ 8, 9, 10]])np.array([[ 2, 3, 4],[ 6, 7, 8],[ 8, 9, 10]])&＃43;np.array([1,2,3])Out[28]:array([[ 3, 5, 7],[ 7, 9, 11],[ 9, 11, 13]])P矩阵与C矩阵每一行之间的距离可以用公式[[p1^2],[p2^2],[...]]-2PC^T&＃43;[c1^2,c2^2,c3^2...](p1^2是p1行向量平方和&＃xff0c;c1^2是c1行向量平方和)"""xx&＃61;rowNorms(x)[:,np.newaxis]#转化为列向量&＃xff0c;便于让dot的每一行都相加同一个数yy&＃61;rowNorms(y)[np.newaxis,:]#与xx同理dot&＃61;np.dot(x,y.T)res &＃61; xx &＃43; yy - 2 * dotreturn res if square else np.sqrt(res)


kmeans&＃43;&＃43;的优化

该方法David Arthur and Sergei Vassilvitskii没有发表出来&＃xff0c;但是在机器学习开源社区上被用到&＃xff0c;即在用kmeas&＃43;&＃43;计算距离的过程中&＃xff0c;每一轮再增加一个最大迭代次数&＃xff0c;比如通过轮盘赌一次选取n_local_trials个&＃xff0c;然后再依次选取能做到使距离平方和最小的聚类中心点&＃xff08;总距离平方和最小中的离现有中心点最大的点&＃xff09;。n_local_trals的选取原则为&＃xff08;2&＃43;log(聚类中心点数)&＃xff09;。
实现过程如下&＃xff08;引用的sklearn开源社区源码&＃xff08;https://github.com/scikit-learn/scikit-learn&＃xff09;&＃xff09;&＃xff1a;

n_samples, n_features &＃61; X.shape
centers &＃61; np.empty((n_clusters, n_features), dtype&＃61;X.dtype)
n_local_trials &＃61; 2 &＃43; int(np.log(n_clusters))
center_id &＃61; random_state.randint(n_samples)
centers[0] &＃61; X[center_id]
closest_dist_sq &＃61; euclidean_distances(centers[0, np.newaxis], X, Y_norm_squared&＃61;x_squared_norms,squared&＃61;True)
current_pot &＃61; closest_dist_sq.sum()
for c in range(1, n_clusters):rand_vals &＃61; random_state.random_sample(n_local_trials) * current_potcandidate_ids &＃61; np.searchsorted(stable_cumsum(closest_dist_sq),rand_vals)np.clip(candidate_ids, None, closest_dist_sq.size - 1,out&＃61;candidate_ids)distance_to_candidates &＃61; euclidean_distances(X[candidate_ids], X, Y_norm_squared&＃61;x_squared_norms, squared&＃61;True)np.minimum(closest_dist_sq, distance_to_candidates,out&＃61;distance_to_candidates)candidates_pot &＃61; distance_to_candidates.sum(axis&＃61;1)best_candidate &＃61; np.argmin(candidates_pot)current_pot &＃61; candidates_pot[best_candidate]closest_dist_sq &＃61; distance_to_candidates[best_candidate]best_candidate &＃61; candidate_ids[best_candidate]centers[c] &＃61; X[best_candidate]return centers

距离计算的优化

假设有4个样本3个特征如下&＃xff1a;

再比如想要把样本聚类为2类&＃xff0c;这样类的聚类中线为1和2。将上面数据存为矩阵&＃xff1a;
x&＃61;np.array([[1,2,3],[4,5,6],[7,8,9],[10,11,12]])
通过euclideanDistance计算所有样本点到两个聚类中心的距离&＃xff0c;如下&＃xff1a;

euclideanDistance(x,x[[0,1]])
Out[4]:
array([[ 0. , 5.19615242],[ 5.19615242, 0. ],[10.39230485, 5.19615242],[15.58845727, 10.39230485]])


结果显示1到样本1的距离0&＃xff0c;到样本2的距离为5.19615242&＃xff0c;2到样本1的距离为5.19615242到样本2的距离为0&＃xff0c;依次类推我们就得到了各个样本到各个聚类中心的距离&＃xff0c;并以此将各个样本贴上不同的类别标签上面的类别标签为[1,2,2,2]。
所以我们再选择的过程中直接比较的是行向量的相对大小&＃xff0c;那如果我们在求||x-x[[0,1]]||的时候行向量同时减去一个||x||^2哪&＃xff1f;分类的结果其实是不变的&＃xff0c;但是对于这里的距离求取过程却可以简化为如下&＃xff1a;

def relativeDist(x,y):#x的每个样本与y之间的相对距离"""我们知道&＃xff0c;如果单纯选出x行向量与y之间最小的距离&＃xff0c;完全可以同时减去xx也就是(x-y)^2-x^2"""yy&＃61;rowNorms(y)[np.newaxis,:]dot&＃61;np.dot(x,y.T)res&＃61;-2*dot&＃43;yyreturn res


该方法在计算距离的时候每次减少的计算量与x的大小成正比&＃xff0c;减少了计算量。
上面例子用relativeDist求取的结果为&＃xff1a;

relativeDist(x,x[[0,1]])
Out[5]:
array([[ -14, 13],[ -50, -77],[ -86, -167],[-122, -257]])


分类结果依然是[1,2,2,2]。