背景與原理：
聚類問題與分類問題有一定的區別，分類問題是對每個訓練數據，我給定了類別的標簽，現在想要訓練一個模型使得對於測試數據能輸出正確的類別標簽，更多見於監督學習；而聚類問題則是我們給出了一組數據，我們並沒有預先的標簽，而是由機器考察這些數據之間的相似性，將相似的數據聚為一類，是無監督學習的一個典型應用。
而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方法，當屬於某個類的數據點過少就刪掉這個類，而當屬於某個類的數據過多、分散程度較大時將這個類分成兩個子類