基于GSA算法改進(jìn)的K均值聚類

婁 奧，姚敏立，袁 丁

(火箭軍工程大學(xué) 作戰(zhàn)保障學(xué)院，陜西 西安 710025)

0 引 言

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

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

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

1 GSA算法

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

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

(1)

(2)

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

(3)

(4)

用于最大值優(yōu)化問題時則與之相反。

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

(5)

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

(6)

(7)

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

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

(8)

其中，rand1表示一個[0,1]區(qū)間內(nèi)的隨機(jī)數(shù)；Kbest初始值為N， 隨迭代次數(shù)增加線性減小至1，定義為

(9)

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

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

(10)

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

(11)

(12)

其中，rand2表示一個[0,1]區(qū)間內(nèi)的隨機(jī)數(shù)。

2 K均值聚類算法

K均值聚類是無監(jiān)督的硬聚類算法，把樣本點(diǎn)到中心點(diǎn)的某種距離作為優(yōu)化目標(biāo)，目的是使簇內(nèi)各個對象的相似度盡可能高，且各簇之間相似度盡可能小。

K均值聚類通常把歐式距離作為相似度測度，取誤差的平方和作為聚類準(zhǔn)則函數(shù)。在算法運(yùn)行時，K均值聚類通常有兩種迭代終止條件供選擇：①聚類準(zhǔn)則函數(shù)值盡可能小；②聚類中心點(diǎn)位置未更新。

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

input: 聚類中心個數(shù)k；N個樣本點(diǎn)組成的數(shù)據(jù)集

步驟1 隨機(jī)選擇k個樣本點(diǎn)作為初始聚類中心。

步驟2 根據(jù)最小距離準(zhǔn)則，把剩余樣本點(diǎn)賦給k個簇。

步驟3 計算每個簇內(nèi)所有點(diǎn)的平均值，得到新的聚類中心集合。

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

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

3 基于GSA算法改進(jìn)的K均值聚類

經(jīng)典K均值聚類隨機(jī)選取初始聚心和人工選擇聚類數(shù)，這嚴(yán)重影響算法的搜索精度和穩(wěn)定性。考慮GSA算法的全局尋優(yōu)能力強(qiáng)，而K均值聚類又有局部精確解的搜索能力，將兩者有機(jī)結(jié)合，優(yōu)勢互補(bǔ)，是本文算法改進(jìn)的研究方向。

3.1 聚類質(zhì)量評價指標(biāo)

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

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

(13)

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

(14)

其中，Ai代表第i個簇的中心點(diǎn)，Xij代表第i個簇內(nèi)第j個樣本點(diǎn)，Ti表示第i個簇內(nèi)樣本點(diǎn)總數(shù)。

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

(15)

DBI指數(shù)越小，說明不同簇的相似度越小，即聚類質(zhì)量越好。

3.2 種群成熟度因子

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

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

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

(16)

(17)

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

(18)

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

(19)

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

(20)

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

3.3 適應(yīng)度函數(shù)

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

(21)

其中，Ai表示第i個聚類中心；Xj是簇Ci內(nèi)的任意樣本點(diǎn)；k是簇的個數(shù)。N是數(shù)據(jù)集樣本點(diǎn)的總數(shù)。

3.4 改進(jìn)后K均值聚類的算法流程

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

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

(22)

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

改進(jìn)后K均值聚類的流程如下所示：

input:N個樣本點(diǎn)組成的數(shù)據(jù)集

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

步驟3 初始化每個粒子的參數(shù)初值；

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

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

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

步驟7 選出適應(yīng)度值最小的粒子，對其用K均值算法求聚類。進(jìn)入步驟9；

步驟8 判斷迭代次數(shù)k是否達(dá)到最大值kmax， 若是則進(jìn)入步驟7，否則回到步驟4；

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

步驟10 選出DBI值最小時的聚類數(shù)nbest。 當(dāng)n=nbest時得到的聚類即算法的輸出結(jié)果。

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

改進(jìn)后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 實(shí)驗驗證

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

實(shí)驗數(shù)據(jù)采用UCI公開標(biāo)準(zhǔn)數(shù)據(jù)庫中的Iris、Glass、Wine和Balance-scale等4種數(shù)據(jù)集。其中，Glass屬于多類別數(shù)據(jù)集，Wine屬于高維數(shù)據(jù)集，Balance-scale屬于多樣本數(shù)據(jù)集。每種數(shù)據(jù)集的基本情況見表1。

表1 UCI數(shù)據(jù)集基本情況

設(shè)本文算法的種群粒子總數(shù)為50。統(tǒng)計聚類后各簇的個體樣本到中心點(diǎn)的距離之和的均值，用來檢驗算法的尋優(yōu)能力和穩(wěn)定性。記錄4種算法20次獨(dú)立運(yùn)算的最小值(Min)，最大值(Max)，平均值(Ave)和標(biāo)準(zhǔn)差(Std)，結(jié)果見表2～表5。

表2 Iris數(shù)據(jù)集實(shí)驗結(jié)果對比

表3 Glass數(shù)據(jù)集實(shí)驗結(jié)果對比

表4 Wine數(shù)據(jù)集實(shí)驗結(jié)果對比

表5 Balance-scale數(shù)據(jù)集實(shí)驗結(jié)果對比

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

對于Glass數(shù)據(jù)集，k-means++聚類的最小值最小，為6.4729。但本文算法的最大值、平均值等都遠(yuǎn)遠(yuǎn)小于最大最小距離、傳統(tǒng)K均值和K-means++聚類，說明本文算法對于多類別數(shù)據(jù)集亦善于求聚類，且聚類質(zhì)量較高。

對于Wine數(shù)據(jù)集，本文算法的平均值最小，標(biāo)準(zhǔn)差也優(yōu)于K均值和K-means++聚類，說明其對多特征數(shù)據(jù)集也能取得不錯的聚類效果。

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

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

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

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

5 結(jié)束語

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