基于GSA算法改進的K均值聚類

婁 奧，姚敏立，袁 丁

(火箭軍工程大學 作戰保障學院，陜西 西安 710025)

0 引 言

傳統K均值聚類存在對初始聚心敏感、全局搜索能力弱和人工設定聚類數等缺點。為改善K均值聚類的性能，歐陽浩等[1]引入小生境和禁忌算法的思想，田詩宵等[2]則提出基于密度峰值優化的算法改進方案，這些方法顯著提升了正確率，但未考慮算法穩定性的問題。Goel等[3]和Pizzuti等[4]引入智能算法思想，把進化算子應用到迭代尋優的過程中，提高了算法魯棒性，但仍存在參數較多、運算繁瑣等問題。

萬有引力搜索算法(gravitational search algorithms，GSA)是伊朗科學家Esmat Rashed等提出的一種啟發式群智能算法。有研究結果表明，當優化基準測試函數時，GSA算法的尋優精度和速度等都明顯優于遺傳算法、模擬退火算法等其它智能算法[5]。

本文針對傳統K均值聚類的缺點，提出一種基于GSA算法的改進K均值聚類。該算法以均方誤差作為GSA的適應度函數，全局搜索聚類質量最優的初始聚類中心；引入種群成熟度因子避免GSA陷入局部最優；設置戴維森堡丁指數為聚類質量評價指標，確定最佳的聚類數。本文將改進后K均值聚類與工程實踐中常用的最大最小距離聚類、傳統K均值聚類、K-means++聚類等算法作實驗對比。結果驗證，該算法具有更好的聚類效果和穩定性。

1 GSA算法

牛頓第二定律指出，自然界任何兩個物體之間都是相互吸引的，它們之間存在的力稱作引力。引力的大小與兩物體質量的乘積成正比，與距離的平方成反比。

GSA算法作為一種通用型優化算法，吸收牛頓第二定律的特點，在求解優化問題時，不僅搜索粒子位置，還考慮粒子質量。粒子質量用來評價粒子位置的優劣，粒子位置越好，質量越大。粒子之間相互吸引并向質量較大的粒子位置方向移動，通過迭代運算，質量最大的粒子位置成為最優解。

(1)

(2)

其中，fitnessi(t) 表示算法第t次迭代時粒子i的適應度值。對于最小值優化問題，最優適應度best_fit(t) 和最差適應度worst_fit(t) 分別為

(3)

(4)

用于最大值優化問題時則與之相反。

第t次迭代時，定義粒子i和粒子j在第d維的相互吸引力大小為

(5)

其中，Mi(t) 和Mj(t) 分別表示粒子i和粒子j的慣性質量；ε為常量；G(t) 表示第t次迭代時的引力系數；Rij(t) 表示粒子i和粒子j之間的歐式距離；計算公式分別為

(6)

(7)

其中，G0是引力系數初值；a是引力系數的衰減因子，一般為定值；T表示最大迭代次數。

第t次迭代時，粒子i在第d維所受的總作用力為

(8)

其中，rand1表示一個[0,1]區間內的隨機數；Kbest初始值為N， 隨迭代次數增加線性減小至1，定義為

(9)

其中，final_per表示對其它粒子產生作用力的粒子百分比。

第t次迭代時，粒子i在第d維上的加速度為

(10)

每次迭代進化時，粒子i的速度和位置由以下公式進行更新

(11)

(12)

其中，rand2表示一個[0,1]區間內的隨機數。

2 K均值聚類算法

K均值聚類是無監督的硬聚類算法，把樣本點到中心點的某種距離作為優化目標，目的是使簇內各個對象的相似度盡可能高，且各簇之間相似度盡可能小。

K均值聚類通常把歐式距離作為相似度測度，取誤差的平方和作為聚類準則函數。在算法運行時，K均值聚類通常有兩種迭代終止條件供選擇：①聚類準則函數值盡可能?。虎诰垲愔行狞c位置未更新。

K均值聚類的具體流程如下：

input: 聚類中心個數k；N個樣本點組成的數據集

步驟1 隨機選擇k個樣本點作為初始聚類中心。

步驟2 根據最小距離準則，把剩余樣本點賦給k個簇。

步驟3 計算每個簇內所有點的平均值，得到新的聚類中心集合。

步驟4 循環步驟2和步驟3，直至聚類中心不變或誤差平方和減小至閾值。

output:k個兩兩之間相似度很低的簇

3 基于GSA算法改進的K均值聚類

經典K均值聚類隨機選取初始聚心和人工選擇聚類數，這嚴重影響算法的搜索精度和穩定性?？紤]GSA算法的全局尋優能力強，而K均值聚類又有局部精確解的搜索能力，將兩者有機結合，優勢互補，是本文算法改進的研究方向。

3.1 聚類質量評價指標

本文引入戴維森堡丁指數(Davies-Bouldin index，DBI)作為聚類質量評價指標獲取最佳聚類數。有研究結果表明，DBI對聚類結果中數據成員變動的敏感性很高，可作為評價聚類數優劣的依據[6]。

DBI指數的實質是度量每個簇最大相似度的均值，其計算公式如下所示

(13)

式中：Si指第i個簇內數據到中心點的平均距離，代表各粒子的位置分散程度，如式(14)所示。N是簇的總數

(14)

其中，Ai代表第i個簇的中心點，Xij代表第i個簇內第j個樣本點，Ti表示第i個簇內樣本點總數。

式(13)中Rij定義為第i個和第j個簇的中心點之間的歐式距離

(15)

DBI指數越小，說明不同簇的相似度越小，即聚類質量越好。

3.2 種群成熟度因子

GSA算法雖然在迭代初期有較高的搜索速度，但臨近收斂時極易陷入局部最優，效率也會大大降低。由文獻[7]可知，可以依據種群多樣性的變化尋找算法收斂的臨界點。常見的判斷種群多樣性大小的方式有平均粒距和粒子之間適應度的差異[8]。但這兩種方法都未考慮粒子未來的運動趨勢，所以并不完善。

根據模糊理論，粒子之間越相似則種群多樣性越低[9]。粒子之間的相似程度可以用它目前位置和歷史最優位置的相似程度來表示。當種群內粒子的相似程度大于閾值時，說明該區域粒子分布較密，種群可能陷入局部最優或早熟收斂。基于此，本文設置種群成熟度因子判斷GSA算法是否早熟收斂：

Y(t)=(Y1(t)Y2(t)Y3(t)…Yi(t)…YN(t))T

(16)

(17)

對矩陣Y(t)作歸一化處理，得到矩陣Y′(t)

(18)

Y′(t) 矩陣中每個行向量可看作一個模糊集合，表示粒子當前位置和歷史最優位置在不同維度上的隸屬度。Y′(t) 中任意兩個模糊集合Y′a(t) 和Y′c(t) 的相似度可以用貼進度σac(t) 表示，記為

(19)

設δ(t) 為種群成熟度因子，它反映的是種群中粒子的平均相似程度，如式(20)所示

(20)

δ(t) 的取值區間為[0,1]。δ(t) 越接近1且大小穩定時，說明在第t次迭代時種群的多樣性越差，算法可能已收斂陷入局部最優。

3.3 適應度函數

在智能優化算法中，適應度函數體現算法尋優的目標和方向[10]。本文使用均方誤差(mean square error，MSE)作為GSA算法的適應度函數，即改進后聚類的目標函數。均方誤差值越小，說明簇內相似度越高，聚類效果越好。本文中均方誤差的定義如下所示

(21)

其中，Ai表示第i個聚類中心；Xj是簇Ci內的任意樣本點；k是簇的個數。N是數據集樣本點的總數。

3.4 改進后K均值聚類的算法流程

為改善傳統K均值聚類對初值敏感而易落入局部最優的問題，本文采用具有較強全局搜索能力的GSA算法優化聚類中心。

GSA是單目標優化算法，求解的目標為全局極值，因此本文視一個聚類中心的集合為GSA的優化對象單元。為方便運算，本文建立粒子編碼模型表示一個聚類中心集合。具體方法如下：設種群中粒子X由聚類中心ai(i=1,2,3…,n) 組成，用一維數組表示，長度為n×m，m是樣本向量的維數。則種群中粒子X的編碼結構為

(22)

改進后K均值聚類主要包括3部分。第一部分，設置初始聚類數n0為2，應用GSA全局搜索適應度函數值最小的初始聚心，再引入標準K均值聚類迭代生成簇和更新聚心；第二部分，遞增聚類數直至其上限nmax， 對不同的n重復GSA和K均值聚類的工作，最終計算nmax-1個聚類結果的DBI值；第三部分，取最小的DBI值對應的n為最佳聚類數，對應的聚類結果為最佳聚類。

改進后K均值聚類的流程如下所示：

input:N個樣本點組成的數據集

步驟2 從數據集中隨機選擇n個樣本作為初始聚類中心，構成一個粒子，粒子形式見式(22)。重復L次，得到L個粒子。對每個粒子使用K均值求聚類，得到L個簇；

步驟3 初始化每個粒子的參數初值；

步驟4 計算所有粒子的慣性質量、所受的作用力和加速度；

步驟5 更新所有粒子的速度和編碼值；

步驟6 計算種群成熟度因子δ(t)。 若超過閾值δ′或與δ(t-1) 相差小于0.001，則進入步驟7，否則進入步驟8；

步驟7 選出適應度值最小的粒子，對其用K均值算法求聚類。進入步驟9；

步驟8 判斷迭代次數k是否達到最大值kmax， 若是則進入步驟7，否則回到步驟4；

步驟9 計算DBI指數，令n=n+1， 返回步驟2。若n>nmax， 則進入步驟10；

步驟10 選出DBI值最小時的聚類數nbest。 當n=nbest時得到的聚類即算法的輸出結果。

output: 兩兩之間相似度很低的簇。

改進后K均值聚類的偽代碼如下所示：

The improved k-means clustering based on GSA

(1) Set the initial value and upper limit of cluster number.

(2) Letnequals 2,selectnsamples randomly to form particles, repeatLtimes.

(3) Identifier particles according to Eq.(22)

(4) Initialize parameters of GSA.

(5) Whilen

(6) Forkfrom 1 tokmaxstep 1

(9) Calculateδ(t) according to Eq.(16)～Eq.(20)

(10) Ifδ(t)>δ′ or |δ(t)-δ(t-1)|≤0.001 then

(11) Break

(12) End if

(13) End for

(14) Cluster the particle withbest_fit(t) according to Eq.(3) and Eq.(21)

(15) Calculate DBI according to Eq.(13)～Eq.(15)

(16)n=n+1

(17) End while

(18) Electnbestwith minimun DBI

(19) Return the cluster whenn=nbest

4 實驗驗證

本文在Intel(R)Core(TM)i3-2100CPU@3.10GHz，內存12 G，Windows7系統和Matlab2014b的環境下，對最大最小距離聚類、K均值聚類、K-means++聚類和本文改進后K均值聚類進行仿真實驗驗證。

實驗數據采用UCI公開標準數據庫中的Iris、Glass、Wine和Balance-scale等4種數據集。其中，Glass屬于多類別數據集，Wine屬于高維數據集，Balance-scale屬于多樣本數據集。每種數據集的基本情況見表1。

表1 UCI數據集基本情況

設本文算法的種群粒子總數為50。統計聚類后各簇的個體樣本到中心點的距離之和的均值，用來檢驗算法的尋優能力和穩定性。記錄4種算法20次獨立運算的最小值(Min)，最大值(Max)，平均值(Ave)和標準差(Std)，結果見表2～表5。

表2 Iris數據集實驗結果對比

表3 Glass數據集實驗結果對比

表4 Wine數據集實驗結果對比

表5 Balance-scale數據集實驗結果對比

對于Iris和Balance-scale數據集，本文算法的最小值、最大值、平均值等指標都最小，說明各簇內個體的間距很小，聚類中心的選取合適。本文算法的標準差接近0，僅次于最大最小距離聚類，說明其性能穩定，對多樣本數據集魯棒性較好。

對于Glass數據集，k-means++聚類的最小值最小，為6.4729。但本文算法的最大值、平均值等都遠遠小于最大最小距離、傳統K均值和K-means++聚類，說明本文算法對于多類別數據集亦善于求聚類，且聚類質量較高。

對于Wine數據集，本文算法的平均值最小，標準差也優于K均值和K-means++聚類，說明其對多特征數據集也能取得不錯的聚類效果。

為進一步檢驗本文改進后K均值聚類的有效性，記錄4種算法20次獨立運算的平均正確率，見表6。

表6 算法平均正確率對比/%

由表6可知，與傳統K均值聚類對比，本文算法對4個數據集測試的平均正確率分別提高0.91、1.40、2.25和2.88個百分點。這充分驗證了本文的K均值聚類改進策略的有效性。

綜上所述，本文改進后K均值聚類對于多樣本、多類別和高維數據集等都有良好的搜索性能和穩定性，與最大最小距離、傳統K均值和K-means++等聚類對比，其聚類結果更理想。

5 結束語

本文提出一種基于GSA算法的改進K均值聚類。采用粒子編碼策略，把聚類中心集合作為GSA算法的尋優對象，克服標準K均值聚類對初始聚心敏感的問題。引入種群成熟度因子，避免算法陷入局部最優。引入聚類質量評價指標，確定最佳聚類數。同時算法把均方誤差作為適應度函數，指引全局尋優方向。實驗結果表明，本文算法具有更高的正確率和更好的穩定性，對不同類型的數據集都有理想的聚類效果。進一步提高算法的聚類正確率及拓展在一些領域的應用，如改進徑向基神經網絡辨識算法，將是下一步的研究方向。