K-Means 是一種最經(jīng)典和常用的聚類方法。它通過多輪迭代的方式不斷更新不同類樣本的中心，計算樣本到每個中心的距離，然后更新樣本所屬的類。最終能夠把樣本劃分到 K 個類中。本案例中，我們首先使用 Python 實現(xiàn) K-Means 算法，基于一份隨機數(shù)據(jù)集，使用動畫演示聚類過程和優(yōu)化目標的變化。然后將 K-Means 應(yīng)用于圖像分割問題。最后我們還將使用一份中文新聞數(shù)據(jù)集，用 K-Means 算法進行自動新聞主題聚類，并使用柱狀圖和詞云圖對聚類結(jié)果進行可視化分析。
K-Means 算法的基本運行流程為：
1 隨機選擇 個點作為初始中心
2 Repeat：
2.1 將每個樣本指派到最近的中心，形成 個類
2.2 重新計算每個類的中心為該類樣本均值
3 直到中心不發(fā)生變化
在整個算法中，2.1 步驟運算量最大，因為該步驟需要計算每一個樣本到 個中心的距離。假設(shè)我們的數(shù)據(jù)集表示為 Pandas 的 DataFrame 格式 X。要實現(xiàn)每個樣本到中心的距離計算，可以有多種實現(xiàn)方式。最簡單的是通過 for 循環(huán)遍歷 X 的每一行，計算每一行到中心的距離。為了提高運算效率，我們可以借助 DataFrame 的 apply 函數(shù)加快運算。Numpy 是 Python 中矩矩陣運算的高效工具，我們也可以使用矩陣運算的方式來實現(xiàn)樣本到中心的距離計算。下面，我們分別來看下這三種實現(xiàn)方法。
iterrows 遍歷的方式實現(xiàn)假設(shè)我們使用歐式距離計算樣本到中心的距離。對于樣本 維樣本 到中心 的歐式距離計算公式為：
使用最簡單的方式來實現(xiàn)，先用一個函數(shù) point_dist 計算一個樣本到中心的距離。這里我們使用 Numpy 的線性代數(shù)模塊 linalg 中的 norm 方法。
import numpy as np
def point_dist(x,c): #定義距離計算函數(shù)
return np.linalg.norm(x-c)
然后使用 iterrows 方法遍歷樣本計算樣本到中心的距離，定義 k_means_iterrows 方法實現(xiàn) K-Means 算法。
def k_means1(X,k):
centers = X.sample(k).values #從數(shù)據(jù)集隨機選擇 K 個樣本作為初始化的類中心，k 行 d 列
X_labels = np.zeros(len(X)) #樣本的類別
error = 10e10
while(error > 1e-6):
for i,x in X.iterrows():#指派樣本類標簽
X_labels[i] = np.argmin([point_dist(x,centers[i,:]) for i in range(k)])
centers_pre = centers
centers = X.groupby(X_labels).mean().values #更新樣本均值，即類中心
error = np.linalg.norm(centers_pre - centers)#計算error
return X_labels, centers
用一個簡單的隨機數(shù)據(jù)集來測試時間性能。Sklearn 中的 datasets 模塊的 make_blobs 函數(shù)能夠自動生成一些供測試聚類算法的隨機數(shù)據(jù)集。它能夠根據(jù)輸入的參數(shù)生成數(shù)據(jù)集和對應(yīng)的類標簽。常用的參數(shù)如下表：
| 參數(shù) | 含義說明 |
|---|---|
| n_samples | 需要生成的樣本數(shù)量,默認值100 |
| n_features | 特征數(shù)量，默認值是2 |
| centers | 數(shù)據(jù)的中心點數(shù)量，默認值3 |
| cluster_std | 數(shù)據(jù)集的標準差，浮點數(shù)或者浮點數(shù)序列，默認值1.0 |
| center_box | 中心確定之后的數(shù)據(jù)邊界，默認值(-10.0, 10.0) |
| shuffle | 打亂樣本順序，默認值是True |
| random_state | 官網(wǎng)解釋是隨機生成器的種子 |
from sklearn import datasets
import pandas as pd
X, y = datasets.make_blobs(n_samples=5000, n_features=8, cluster_std = 0.5,centers=3,random_state=99)
X_df = pd.DataFrame(X)
在該數(shù)據(jù)集上用我們實現(xiàn)的 k_means1 方法運行 K-Means 聚類。使用 iPython 提供的魔法命令 %time 記錄運行時間。
%time labels,centers = k_means1(X_df,3) # for 循環(huán)
CPU times: user 5.16 s, sys: 39.3 ms, total: 5.19 sWall time: 5.14 s
apply 遍歷的方式實現(xiàn)提高運算效率，可以使用 DataFrame 的 apply 函數(shù)，它可以對數(shù)據(jù)框中的每一行執(zhí)行一個復(fù)雜的函數(shù)。在我們的例子中，是計算每一行與每一個中心的距離。
def k_means2(X,k):
#初始化 K 個中心，從原始數(shù)據(jù)中選擇樣本
centers = X.sample(k).values
X_labels = np.zeros(len(X)) #樣本的類別
error = 10e10
while(error > 1e-6):
#********#
X_labels = X.apply(lambda r : np.argmin([point_dist(r,centers[i,:]) for i in range(k)]),axis=1)
centers_pre = centers
centers = X.groupby(X_labels).mean().values #更新樣本均值，即類中心
error = np.linalg.norm(centers_pre - centers)#計算error
return X_labels, centers
%time labels,centers = k_means2(X_df,3) # apply 運算
CPU times: user 3.24 s, sys: 19.9 ms, total: 3.26 sWall time: 3.23 s
數(shù)據(jù)集表示成 矩陣 ，其中 為樣本數(shù)量， 為樣本的維度。 個聚類中心表示成 矩陣 ， 每一行表示一個聚類中心。樣本到 個中心的距離表示成 矩陣 。
已經(jīng)聚類中心，計算樣本到中心距離，并將樣本劃分到距離最小的類的流程如下圖所示。
使用 Numpy 實現(xiàn)上述計算流程的代碼為：
for i in range(k):
D[:,i] = np.sqrt(np.sum(np.square(X - C[i,:]),axis=1))
labels = np.argmin(D,axis=1)
得到樣本的類標簽后，聚類中心的更新流程為：1）根據(jù)類標簽對樣本進行分組；2）將聚類中心更新為每一組樣本的均值。Python 實現(xiàn)的代碼為：
C = X.groupby(labels).mean().values
現(xiàn)在我們更新 K-Means 算法的實現(xiàn)，函數(shù)名為 k_means。
import pandas as pd
import numpy as np
def k_means(X,k):
C = X.sample(k).values #從數(shù)據(jù)集隨機選擇 K 個樣本作為初始化的類中心，k 行 d 列
X_labels = np.zeros(len(X)) #記錄樣本的類別
error = 10e10 #停止迭代的閾值
while(error > 1e-6):
D = np.zeros((len(X),k)) #樣本到每一個中心的距離，n 行 k 列
for i in range(k):
D[:,i] = np.sqrt(np.sum(np.square(X - C[i,:]),axis=1))
labels = np.argmin(D,axis=1)
C_pre = C
temp_C = X.groupby(labels).mean() #更新樣本均值，即類中心
C = np.zeros((k,X.shape[1]))
for i in temp_C.index:
C[i,:] = temp_C.loc[i,:].values
if C.shape == C_pre.shape:
error = np.linalg.norm(C_pre - C)#計算error
else:
print(C.shape, C_pre.shape)
return labels, C
%time labels,centers = k_means(X_df,3) # 矩陣運算
CPU times: user 150 ms, sys: 0 ns, total: 150 msWall time: 149 ms
下面我們使用一份隨機生成的二維數(shù)據(jù)集，使用我們上一小節(jié)實現(xiàn)的 k_means 完成聚類，然后使用不同顏色標注不同類的樣本以及類中心。
color_dict = {0:'#E4007F',1:'#007979',2:'blue',3:'orange'} #洋紅，深綠，藍色，橘色
from sklearn import datasets
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline
X, y = datasets.make_blobs(n_samples=1000, n_features=2, cluster_std = 1.5,centers=4,random_state=999)
X_df = pd.DataFrame(X,columns=['x1','x2'])
labels,centers= k_means(X_df,4)
fig, ax = plt.subplots(figsize=(8, 8)) #設(shè)置圖片大小
for i in range(len(centers)):
ax.scatter(X_df[labels == i]['x1'],X_df[labels == i]['x2'],color=color_dict[i],s=50,alpha=0.4)
ax.scatter(centers[int(i),0],centers[int(i),1],color='r',s=100,marker='+')
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
Text(0, 0.5, '')
要動態(tài)展示 K-Means 聚類過程，我們需要在每一步迭代中記錄每一個類的中心，以及每一個類的樣本集合。創(chuàng)建 k_means_steps，在完成聚類的同時，將每一步迭代的每類樣本和中心返回。
def k_means_steps(X,k):
#初始化 K 個中心，從原始數(shù)據(jù)中選擇樣本
#********#
samples_list = [] #記錄每一個中間迭代中每一類樣本
centers_list = [] #記錄每一個中間迭代中每一類樣本中心
#********#
C = X.sample(k).values
labels = np.zeros(len(X)) #樣本的類別
error = 10e10
while(error > 1e-6):
D = np.zeros((len(X),k)) #樣本到每一個中心的距離
for i in range(k):
D[:,i] = np.sqrt(np.sum(np.square(X - C[i,:]),axis=1))
labels = np.argmin(D,axis=1)
C_pre = C
C = X.groupby(labels).mean().values #更新樣本均值，即類中心
#********# 記錄當前迭代地每一類的樣本集合和中心
samples,centers2 = [],[]
for i in range(k):
samples.append(X[labels == i])
centers2.append(C[i,:])
samples_list.append(samples)
centers_list.append(centers2)
#********#
if C.shape == C_pre.shape:
error = np.linalg.norm(C_pre - C)#計算error
else:
print(C.shape, C_pre.shape)
return labels, C,samples_list,centers_list #********# 返回最終的聚類結(jié)果，聚類中心，每一步的聚類結(jié)果和聚類中心
我們可以借助 matplotlib.animation 動畫模塊來實現(xiàn)下面的 init_draw 函數(shù)是動畫最開始時繪制的內(nèi)容，包含數(shù)據(jù)。update_draw 則是每次更新的內(nèi)容。
labels,centers,samples_list,centers_list= k_means_steps(X_df,4)
fig, ax = plt.subplots(figsize=(8, 8))
samples_obj = []
centers_obj = []
def init_draw(): # 展現(xiàn)樣本數(shù)據(jù)
ax.set_title('K-Means聚類過程:')
for i in range(len(centers)):
samples_obj.append(ax.scatter(samples_list[0][i]['x1'],samples_list[0][i]['x2'],color=color_dict[i],s=50,alpha=0.6))
centers_obj.append(ax.scatter(centers_list[0][i][0],centers_list[0][i][1],color='r',s=100,marker='+'))
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
def update_draw(t): # 實現(xiàn)動畫中每一幀的繪制函數(shù)，i為第幾幀
ax.set_title('K-Means聚類過程:' + str(t))
samples,centers = samples_list[t],centers_list[t]
for i in range(len(centers)):
samples_obj[i].set_offsets(samples[i])
centers_obj[i].set_offsets(centers[i])
plt.close()
#演示決策面動態(tài)變化
import matplotlib.animation as animation
from IPython.display import HTML
animator = animation.FuncAnimation(fig, update_draw, frames= range(1,len(centers_list)), init_func=init_draw,interval=2000)
HTML(animator.to_jshtml())
我們來查看隨著迭代的進行，K-Means聚類模型的優(yōu)化目標，即失真度量 的變化。對算法進行修改，記錄每一步的失真度量。
import pandas as pd
import numpy as np
def k_means_inertia(X,k):
#初始化 K 個中心，從原始數(shù)據(jù)中選擇樣本
C = X.sample(k).values
labels = np.zeros(len(X)) #樣本的類別
inertia_list = [] #*****記錄優(yōu)化目標****#
error = 10e10
while(error > 1e-6):
D = np.zeros((len(X),k)) #樣本到每一個中心的距離
for i in range(k):
D[:,i] = np.sqrt(np.sum(np.square(X - C[i,:]),axis=1))
labels = np.argmin(D,axis=1)
inertia_list.append(np.square(np.min(D,axis=1)).sum()) #****記錄當前步驟的失真度量****#
C_pre = C
temp_C = X.groupby(labels).mean() #更新樣本均值，即類中心
C = np.zeros((k,X.shape[1]))
for i in range(len(temp_C)):
C[i,:] = temp_C.loc[i,:].values
if C.shape == C_pre.shape:
error = np.linalg.norm(C_pre - C)#計算error
return labels,C,inertia_list
在隨機數(shù)據(jù)上進行聚類，并將失真度量的變化以折線圖的形式繪制出來。
X, y = datasets.make_blobs(n_samples=1000, n_features=2, cluster_std=1,centers=3,random_state=99)
X_df = pd.DataFrame(X,columns=['x1','x2'])
labels,centers,inertia_list = k_means_inertia(X_df,3)
fig, ax = plt.subplots(figsize=(9, 6)) #設(shè)置圖片大小
t = range(len(inertia_list))
plt.plot(t,inertia_list,c='#E4007F',marker='o',linestyle='dashed')
plt.xlabel('t')
plt.ylabel('inertia')
plt.xticks(t)
plt.title('K-Means算法優(yōu)化目標的變化')
Text(0.5, 1.0, 'K-Means算法優(yōu)化目標的變化')
我們先加載一張測試圖片，將圖片打印出來。使用 PIL.Image.open 方法來打開圖片，然后使用 matplotlib 中的 imshow 方法將圖片可視化。
from PIL import Image
fig, ax = plt.subplots(figsize=(6, 5)) #設(shè)置圖片大小
path = './input/timg.jpg'
img = Image.open(path)
plt.imshow(img)
plt.box(False) #去掉邊框
plt.axis('off')#不顯示坐標軸
(-0.5, 599.5, 514.5, -0.5)
將一張圖片轉(zhuǎn)換成表格形式。每一行為一個像素，三列分別為像素的 R，B，G取值。獲取圖片的每一個像素 的 RBG 值可以使用 Image 類的 getpixel 方法。
import pandas as pd
def image_dataframe(image): #將圖片轉(zhuǎn)換成DataFrame,每個像素對應(yīng)每一行，每一行包括三列
rbg_values = []
for i in range(image.size[0]):
for j in range(image.size[1]):
x,y,z= image.getpixel((i,j)) # 獲取圖片的每一個像素 (i,j)(i,j) 的 RBG 值
rbg_values.append([x,y,z])
return pd.DataFrame(rbg_values,columns=['R','B','G']),img.size[0],img.size[1]
img_df,m,n = image_dataframe(img)
打印轉(zhuǎn)換的數(shù)據(jù)框如下：
img_df.head()
| R | B | G | |
|---|---|---|---|
| 0 | 247 | 248 | 243 |
| 1 | 247 | 248 | 243 |
| 2 | 247 | 248 | 243 |
| 3 | 247 | 248 | 243 |
| 4 | 247 | 248 | 243 |
print(m,n,m*n,len(img_df))
600 515 309000 309000
使用我們在上一節(jié)實現(xiàn)的 K-Means 算法對像素進行聚類。
labels, _ = k_means(img_df,2)
將生成的灰度圖可視化，對圖像可視化使用 plt.imshow 方法。
fig, ax = plt.subplots(figsize=(6, 5)) #設(shè)置圖片大小
labels = labels.reshape((m,n))
pic_new = Image.new('L',(m,n))
#根據(jù)類別向圖片中添加灰度值
for i in range(m):
for j in range(n):
pic_new.putpixel((i,j),int(256/(labels[i][j] + 1)))
plt.imshow(pic_new)
plt.box(False) #去掉邊框
plt.axis('off')#不顯示坐標軸
(-0.5, 599.5, 514.5, -0.5)
img_from_labels ，將像素聚類類別標簽，轉(zhuǎn)換成一張灰度圖。def img_from_labels(labels,m,n):
labels = labels.reshape((m,n))
pic_new = Image.new('L',(m,n))
#根據(jù)類別向圖片中添加灰度值
for i in range(m):
for j in range(n):
pic_new.putpixel((i,j),int(256/(labels[i][j] + 1)))
return pic_new
調(diào)整聚類數(shù)量 , 將聚類得到的不同的灰度圖使用 Matplotlib 將生成的灰度圖繪制出來。
fig, ax = plt.subplots(figsize=(18, 10)) #設(shè)置圖片大小
img = Image.open(path) #顯示原圖
plt.subplot(2,3,1)
plt.title('原圖')
plt.imshow(img)
plt.box(False) #去掉邊框
plt.axis('off')#不顯示坐標軸
for i in range(2,7):
plt.subplot(2,3,i)
plt.title('k=' + str(i))
labels, _ = k_means(img_df,i)
pic_new = img_from_labels(labels,m,n)
plt.imshow(pic_new)
plt.box(False) #去掉邊框
plt.axis('off')#不顯示坐標軸
首先，我們使用 Pandas 的 read_csv 方法將中文新聞?wù)Z料讀取進來， encoding 參數(shù)和 sep 參數(shù)分別設(shè)置為 'utf8' 和 '\t' 。
import pandas as pd
news = pd.read_csv('./input/chinese_news_cutted_train_utf8.csv',sep='\t',encoding='utf8')
news.head()
| 分類 | 分詞文章 | |
|---|---|---|
| 0 | 科技 | “ 一路 聽 天下 ” 開拓 “ 無聊 經(jīng)濟 ” 在 納斯達克 風光 上市 的 分眾傳媒... |
| 1 | 旅游 | 中國 赴美 旅游 首發(fā) 團 ６ 月 １ ７ 日 啟程 報價 兩萬左右 六月 初 的 ... |
| 2 | 新聞 | 科學家 教育家 蒙特 梭利 ：給 孩子 愛 與 自由 （ 圖 ） 蒙特 梭利 和 “ ... |
| 3 | 教育 | ７ 名 大學生 作出 “ 震后 恢復(fù) 建議 ” 獲 國務(wù)院 肯定 大學生 災(zāi)后 重建... |
| 4 | 教育 | 最牛 高 考生 ：重慶 學生 稱 “ 押 中 ” 今年 作文題 最牛 鼓勁 ：數(shù)學考... |
news['分類'].value_counts()
汽車 544
科技 511
旅游 510
健康 492
文化 491
房地產(chǎn) 480
財經(jīng) 469
新聞 458
體育 437
教育 437
娛樂 370
女人 322
Name: 分類, dtype: int64
使用 sklearn.feature_extraction.text 模塊的 TfidfVectorizer，將詞列表格式的新聞轉(zhuǎn)換成向量形式。向量中的每一個維度代表字典中的一個詞，維度的取值代表詞在對應(yīng)文檔中的 TF-IDF 取值。
#加載停用詞
stop_words = []
file = open('./input/stopwords.txt')
for line in file:
stop_words.append(line.strip())
file.close()
from sklearn.feature_extraction.text import TfidfVectorizer
vectorizer = TfidfVectorizer(stop_words=stop_words,min_df=0.01,max_df=0.5,max_features=500)
news_vectors = vectorizer.fit_transform(news['分詞文章'])
k_means 方法對新聞進行聚類在真實的文本應(yīng)用中，文本表示成向量后，由于詞典的詞通常會很大，我們得到的是稀疏矩陣。例如上一步 news_vectors 就是一個稀疏存儲格式。
type(news_vectors)
scipy.sparse.csr.csr_matrix
在本案例中，為了簡便我們直接使用 todense 方法將稀疏向量轉(zhuǎn)換成稠密矩陣。(!注意，實際中很少這么操作！)
news_df = pd.DataFrame(news_vectors.todense())
現(xiàn)在我們可以直接使用第一節(jié)實現(xiàn)的 k_means 方法對新聞進行聚類，這里我們將聚類個數(shù)設(shè)置成 12 。
labels, _ = k_means(news_df,12)
在 news 中新建 labels 列保存將得到的樣本聚類結(jié)果。
news['labels'] = labels
在新聞數(shù)據(jù)集中，由于 分類 一列已經(jīng)標注了新聞的主題類別，我們可以聚類得到的每一組新聞中的主題分布，從而定性地觀察聚類效果。將聚類結(jié)果按照聚類結(jié)果進行分組，可以使用 DataFrame 的 groupby 方法。然后使用 value_counts 方法統(tǒng)計每一組的不同主題的新聞數(shù)量分布。
fig, ax = plt.subplots(figsize=(18, 12)) #設(shè)置圖片大小
cluster_topics = news.groupby(['labels'])['分類'].value_counts()
for i in range(12):
plt.subplot(3,4,i+1)
plt.title('cluster ' + str(i))
topic_dist = pd.Series(data=0,index = news['分類'].unique())
topic_dist += cluster_topics[i]
topic_dist.plot(kind='bar')
plt.xticks(rotation=45)
plt.tight_layout()
從上圖可以看到，部分分組對應(yīng)于特定主題，然而在另一上部分分組中，不同主題的新聞分布較為分散。
首先，使用 groupby 方法，根據(jù)聚類結(jié)果對新聞進行分組。
cluster_contents = news[['labels','分詞文章']].groupby(['labels'])
將 wordcloud 模塊導(dǎo)入詞云類 WordCloud，新建一個詞云對象 cloud 。為了顯示效果，構(gòu)造詞云對象時需要設(shè)置停用詞詞表。由于我們顯示的是中文，為了避免亂碼，還需要制定顯示字體。其中字體文件為 ./input/simfang.ttf 。
from wordcloud import WordCloud
stop_words_set = set(stop_words)
cloud = WordCloud(font_path = './input/simfang.ttf',background_color='white', max_words=100,stopwords=stop_words_set,width=600,height=400)
def series_to_word_list(ts):
results = []
for words in ts.values:
results.extend(words.split(' '))
return ' '.join(results)
遍歷每一個聚類的新聞樣本，得到新聞內(nèi)容的內(nèi)容，然后顯示成詞云圖。
fig, ax = plt.subplots(figsize=(18, 12)) #設(shè)置圖片大小
for cluster, group in cluster_contents:
plt.subplot(3,4,cluster+1)
plt.title('cluster ' + str(cluster))
wc = cloud.generate_from_text(series_to_word_list(group['分詞文章']))
plt.imshow(wc)
plt.axis('off')
plt.tight_layout()
可以看到，我們可以直接根據(jù)詞云圖對得到的聚類進行解釋。
本案例中我們使用首先使用三種方式實現(xiàn)了 K-Means 算法并對不同實現(xiàn)的時間性能進行了對比，結(jié)果發(fā)現(xiàn)向量化實現(xiàn)能夠大大提高運行效率。然后我們使用 K-Means 算法進行圖像分割，展示了不同 K 取值下生成的灰度圖的變化。最后，我們使用 K-Means 在一份中文新聞數(shù)據(jù)集進行了主題聚類。
本案例使用的主要 Python 工具如下：
| 工具包 | 用途 |
|---|---|
| NumPy | 矩陣運算 |
| Pandas | 數(shù)據(jù)讀取與預(yù)處理 |
| Matplotlib | 數(shù)據(jù)集可視化、聚類結(jié)果動畫 |
| Sklearn | 中文新聞的向量化 |
| Wordcloud | 繪制聚類結(jié)果詞云圖 |
