基于最小聚類求解k-means問題算法

王守強，朱大銘

（1. 山東交通學院 信息工程系，山東 濟南 250023；2. 山東大學 計算機科學與技術學院，山東 濟南 250100）

1 引言

給定 d維空間中的點集P，k-means聚類問題要求在空間中選取k個中心點，使P中點與其距離最近的中心點的距離平方和最小。形式化描述為

實例：點集 P ∈Rd，正整數k∈Z+。

k-means問題相當于在d維空間中計算k個中心點，以中心點為核心將給定點集P劃分為k個子集，優化目標為給定點到其所屬子集中心點的距離平方和最小。該問題是NP-Hard問題[1]。其教科書算法為Lloyd給出的啟發式算法[2,3]，Lloyd算法簡單而容易實現，但運行結果依賴于初始值，算法無法保證一個確切的求解近似度。Kanungo[4]等采用局部搜索技術給出k-means問題（9+ε）近似度算法。在執行算法前需要對點集空間結構進行劃分求得一個候選中心點集[5]，候選中心點集的求解十分復雜，算法顯得不夠實用。Song等[6]進一步證明，如果將給定點集P作為中心點的候選點集，對候選點集執行局部搜索，可使算法的近似度達到O(1)。2004年Kumar[7]等人給出求解k-means問題的(1+ε)隨機近似算法，算法的時間復雜度為

2) 對于m最小聚類問題，改進了Kumar提出的 k-means問題的(1+ε)隨機近似算法，將 Kumar給出算法的時間復雜度由 O ( 2(k/ε)O(1)dn)改進為證明改進算法求到(1+ε)近似解的成功概率至少為如果該算法運行次，則以接近于 1的概率求到該問題的(1+ε)近似解。

3) 設計出求解k-means問題的局部搜索隨機算法。從k個初始中心點的選取以及生成候選中心點集2個方面分別對Song的局部搜索算法提出新的隨機策略。Song局部搜索算法的時間復雜度為本文算法的時間復雜度為O(nk3dln(k)2)/α)。實驗結果表明：新算法所求解的精度和算法的運算時間均優于 Song給出的局部搜索算法。

2 符號標記與基本結論

定義 2對于點集 P的 m最小聚類問題，稱α=km/|P|為該問題的最小聚類參數。

定義3設是d維空間中的k個點，若存在 l個實數滿足：

d對于 d維空間中某點 x ∈ R ，如果成立，則稱 x為的凸組合點。

引理1[8]點集P的質心點為P的1-means問題最優解的中心點。

Inaba等[7]指出在歐氏空間中，只需從P中隨機選取部分點，計算取樣點的質心點，則該質心點以較大的概率成為P的 ε 近似質心點。該定理描述如下。

引理2[7]給定點集P，從P中均勻地隨機選取m個點，設S為取樣點的集合，m=|S|， c(S)為S的質心點。存在δ (0＜δ＜1)，使下述不等式成立的概率至少為1-δ。

引理3[8]給定點集P，對于任一個點x ∈ Rd，

3 基于取樣技術k-means近似算法

3.1 近似度期望值為2的算法

記 ?k2(P )為給定實例點集 P的最優解的解值。在給出該算法之前，首先探討當中心點集包含每個最優子集 Pi(1≤i≤k)中的一個點時，算法近似度的期望值。

定理 1設 P所對應的最優聚類劃分子集為如果從每個 Pi中均勻地隨機選取一個點ci，以這k個點作為中心點求到的解值記為?，則

證明設 P的最優解所對應的中心點為從每個 Pi中均勻地隨機選取一個點 ci，由所有 ci構成中心點集合記為不失一般性，假設 ci∈Pi，根據已知條件，Pi中任一點被選作中心點 ci的概率為針對點 x ∈Pi,定義{distance(x,ci)}，D(x)實際表示為點x與集合 C - { ci}中最近點的距離，則由此可得：

(根據引理3)

給定正整數m，如果P為m最小聚類問題。下面證明只需從P中隨機選取一個樣本點集S，則S以較高的概率滿足它包含每個最優子集 Pi中至少一個點，即

定理2如果從P中均勻地隨機選取k ln(2k)個點，記取樣點集為S，那么對于所有P，αi均成立的概率至少為1，即?iPr[|S∩P|2i

證明不失一般性，設|記定義則 nis的期望值根據Chernoff Bounds不等式得：

其中，0＜λ＜1。定義事件Ai為滿足足 nis小于λ|S |nni的事件，則：

對每個最優子集Pi，如果均成立，枚舉S中所有k個點的子集，則至少存在一個子集滿足：k個點分別來自于最優子集 P1,P2,… ,Pk中的一個點。以這k個點作為中心點，根據定理1知，所求解的近似性能比的期望值至多為 2。據此，給出下述算法。

算法1k-means問題近似算法

輸入：點集P，參數α。

2) 從S中任取k個點，以這k個點作為中心點，對P進行一次劃分，計算劃分后的值?。

3) 重復2)的操作，枚舉所有可能k個點，選擇? 最小的值所對應的k個中心點作為算法最終解。

4) 返回所求的中心點。

定理3算法1至少以1/2的概率求到近似度期望值為2的解。

證明根據定理1及定理2，該定理結論顯然成立。

算法在求得k個中心點時，通過枚舉取樣點集中所有k個點的子集找出最小解。枚舉子集的個數為 C|kS|，由于將S值代入該式可得枚舉子集個數為由于每次將實例點集P中的點分配給k個中心點時，時間復雜度為 O(nkd)。因此，算法 1的時間復雜度為

3.2 k-means聚類(1+ε)近似算法

定理4對滿足m最小聚類問題實例點集P，如果從P中均勻地隨機選取8 k ln(4k)個點，記取樣εα點集為S，則S包含每個最優子集Pi中至少2/ε個點的概率大于等于3/4。

證明證明過程類似于定理 2。記 ni=|Pi|，定義則 nis的期望值根據Chernoff Bounds不等式得：

定義事件Ai為 nis小于λ|S |nni的事件，則：

由于|Pk|≥m，根據最小聚類參數定義α=km/|P|可得成立。

引理4給定點集P，設均屬于P中同一個最優子集Pi中的點，如果x∈P并且 x是由所組成凸組合的點，則x∈Pi

證明由于 x是屬于中凸組合內的點，根據凸組合定義，存在l個大于等于0的數值滿 足 ：其 中設 ci為 Pi的質心點，則：

算法2k-means問題(1+ε)近似算法

輸入：點集P，參數ε，每個子集最少聚類個數m。

1) 計算參數α＝(mk)/|P|，置C=φ。

3) 調用函數 Cost=Online-k-means(S, C, ε)。

4) 返回代價Cost及集合C。

函數 Online-k-means(·)的實現采用遞歸方法，從i＝1開始，首先從S中枚舉所有2/ε個點的子集，則至少存在一個子集 Si′滿足 Si′?Si，由 Si′求出 Si′′并進一步計算 Si′∪ Si′的質心點 ci′，然后從S中刪除屬于集合 Si′∪ Si′中的所有點。刪除 Si′∪ Si′中的所有點的目的在于當求解下一個子集時，可以減少枚舉子集的個數，從而提高算法的運行效率。函數的實現過程描述如下。

函數 Online-k-means(S, C, ε)

輸入：樣本點集S，已求中心點集C，參數ε 值。

1) If |C|=k then

2) 以C作為中心點，計算解值Sum。

3) if Sum＜MinCost then MinCost=Sum，保存C和MinCost。

4) If |C|＜k and |S|＜2/ε 返回。

5) Repeat。

6) 從S中取2/ε個點,設點的集合為S′。

7) 從S-S′中找出滿足S′凸組合中的點，設點集為 S ''。

9) Online-k-means(S,C, ε)。

10) Until 窮舉完畢S中所有2/ε個點的子集。

11) 返回MinCost的解值。

引理 5如果對每個最優子集 Pi，均成立，則函數 Online-k-means(·)能夠從 S中求出一個子集 Si′，滿足成立。

證明記由上述第 6)步知，算法是從 S中枚舉所有 2/ε個點，則枚舉子集中必存在一個子集S′屬于Si。上述第7）步求出S中滿足S′的凸組合的點集 S′，根據引理 4，，因此并且成立。

引理 6對每個最優子集 Pi，如果成立，則函數 Online-k-means(·)求出 k-means問題(1+ε)近似解的成功概率不小于(1/2)k。

證明根據引理5，對每個最優子集Pi，函數Online-k-means(·)從 S中求出一個子集 Si′，Si′滿足成立。以 Si′的質心點作為 Pi的中心點，根據引理2，該質心點至少以1/2的概率滿足Pi的 1-means問題的(1+ε)近似解。因此，算法求解 k-means問題(1+ε)近似解的成功概率至少為(1/2)k。

定理5對于滿足m最小聚類的點集P，算法2至少以2的概率求出該問題的(1+ε)近似算法，算

2k+2法的時間復雜度為

證明根據定理4，對任意最優子集成立的概率不小于 3/4。當條件成立時，由引理 6知，函數 Online-k-means(·)求解k-means問題(1+ε)近似解的成功概率至少為(1/2)k。因此，算法2的成功概率為

記 r=2/ε，T(|S|)代表求函數 Online-k-means(·)中以 S作為實例點集的枚舉個數，則：由此可得T(|S|)至多為

將|S|及r值代入上式可得：枚舉所有可能k個中心點的子集個數至多為因此整個算法的時間復雜度為

3.3 局部搜索隨機算法

基于文獻[6]的局部搜索算法，本文提出k-means問題的局部搜索的隨機算法，其隨機策略主要體現如下。

1) 文獻[6]候選中心點集選自給定實例點集P，新的隨機算法則以 P的一個取樣子集作為候選中心點集，S滿足以較高的概率包含每個最優子集至少一個點。

2) 在 k個初始中心點的選取方面。文獻[6]是從候選中心點集中任取k個點，新的隨機算法則采用非均勻地隨機選取策略從P中選取k個點，使得這k個點盡可能地分屬于k個不同最優子集中的點。

初始中心點選取算法實現描述如下。

算法3k個初始中心點選取算法

輸入：點集P

1) 從P中隨機選取2個點 x1、 x2，選擇概率

2) While 選取點數≤k。

3) 從P中隨機選取一點xi(i≥3)，選取概率為

4) End While。

5) 返回選取點{x1,…, xk}。

算法3首先從P中隨機選取2個點 S = { x1,x2}，遵循兩點距離越大，被選取概率越大的原則。再依次隨機選取點 x3,… ,xk，選取第i(i＞2)個點時，遵循一個點與已選擇的點距離平方和越大，則該點被選取概率越大的原則。因此算法規定第一次選擇2個點{x1,x2}的概率表達式為設規定選擇第i個點xi的概率表達式為算法4k-means局部搜索隨機算法

輸入：點集P，最小聚類參數α。

3) Repeat。

4) 以C為中心點，對S作一次劃分，劃分子集為{S1,S2,…,Sk}。

5) For i=1 to k。

6) 對于Si中每一個點x，以C′為中心點,計算給定點集P的k-means解值?′,如果?′＜?,置

7) Next i。

9) 返回C以及解值?。

與文獻[6]局部搜索算法相比，選取k個初始中心點時，算法4中的第2)步不是從給定實例點集P中任取k個點，而是調用算法3從P中非均勻地隨機選取k個初始中心點。算法4中的第3)～8)步局部搜索時以取樣子集S作為候選中心點集。根據定理2，該取樣子集能夠較好地代表P。在執行中心點交換時，取樣子集能夠減少交換次數，從而達到降底算法時間復雜度的目的。

根據算法4第6)步，算法每次迭代時，所求P的值下降(1-ε/k)倍。記?k2(P)為給定實例點集 P的最優解值，?2(P,C)為算法所求最終解值，?2(P,C0)為算法初始中心點所對應的值，算法的迭代次數記為t。則：

由此可得：算法迭代次數

4 實驗結果

本文選用 Iris、RuspIni、Spath Postal Zone Data、Cloud 和SPAM 數據集測試本文相關算法。選用數據集 Iris、RuspIni、Spath Postal Zone Data來驗證算法1的運算結果；選用UCI中2個高維數據集Cloud以及SPAM測試局部搜索隨機算法的運算性能。表1列出5個數據集的基本屬性。

表1 選用數據集說明

4.1 近似度算法實驗

根據算法 1，在執行算法前，需要給出算法的最小聚類參數α。受篇幅所限、僅給出α=0.5作為參數值，對不同的k值進行實驗結果。表2中最優值一列來自于文獻[10～12]。由于算法的隨機性，表2中實驗結果取自算法運算 20次后所求解值的最小值，實驗結果參見表2。

表2 算法1實驗結果(α＝0.5)

在表2中，近似度一列值等于實驗結果與最優值的比值。由表2近似度一列可知：本文實驗所得到的算法近似度均小于2。通過表2實驗結果可以看出：對隨機取樣點集S，枚舉S中所有可能k個點子集作為中心點，則以較高概率獲得近似度期望值為2的解。1的常數，因此，迭代此數可簡化為O(kln(k))。由算法4的第7)步知，每次迭代時，需要涉及|S|個中心點交換；每次交換后，算法重新計算k-means解值的時間復雜度為O(nkd)。所以，整個算法4的時

4.2 改進局部搜索算法實驗

為檢驗改進后局部搜索算法性能，本文將數據集Cloud以及SPAM分別應用在文獻[6]局部搜索算法以及本文算法4上進行測試。表3給出2種算法的實驗結果，由于初始中心點選取隨機性，將程序執行20次。表2中運算值所對應的兩列是指兩算法運算20次后所求的最小值，而表2中運算時間所對應的2列值則指程序20次運算后的平均時間。

為便于比較和描述，表中 2運算值=算法實際運行結果/點集個數。

表3 局部搜索算法實驗結果

相對于文獻[6]，定義算法 4的改進值=[1-(算法4)/(文獻[6]算法)]×100%。對于Cloud數據集，當k=10、25、50時，算法2第2）步的運算時間分別改進了96.75%、93.32%以及89.21%，因此算法時間得到顯著提高。除運算時間外，局部搜索隨機算法的運算值也均小于文獻[6]的算法。對于SPAM數據，當 k=10、25、50時，算法所求的解值分別改進了51.91%、87.38%和91.63%，而算法的運算時間則分別改進了98.98%、97.98%以及96.39%。由表3可以看出，針對給定的實例點集P，采用局部搜索隨機算法，無論是算法所求值的精度還是算法的運算時間均優于文[6]所給出的算法。

5 結束語

本文分析了基于最小聚類k-means問題的隨機近似算法，利用取樣技術，給出了該問題近似度期望值為2的隨機算法。同時探討了該子問題的(1+ε)近似算法的求解方案，將 Kumar所給出的(1+ε)近似方案的時間復雜度進行了改進，并分析了算法的成功概率。利用取樣技術，本文設計局部搜索隨機算法。最后，選取了部分實驗數據，對算法近似度以及局部搜索隨機算法進行了驗證。但本文還有幾個問題需要進行探討。1) 本文近似算法是通過枚舉樣本點集中部分點的求到的，顯然算法的時間復雜度高，能否不需要進行組合，而是通過某種策略直解在取樣子集中找出滿足給定條件的某些點？2) 該算法能否進一步提高成功概率？3) 如何找出k個初始中心點，使得這k個點以較高的概率分別來自于k個不同最優子集中的一個點。

[1] DRINEAS P, FRIEZE A, KANNAN R, et al. Clustering large graphs via the singular value decomposition[J]. Machine Learning, 2004,56(1-3)∶9-33.

[2] MACQUEEN J B. Some methods for classification and analysis of multivariate observations[A]. Proceedings of the 5th Berkeley Symposium on Mathematical Statistics and Probability[C]. California, USA,1967. 281-297.

[3] LLOYD S P. Least squares quantization in PCM [J]. IEEE Transactions on Information Theory, 1982, 28(2)∶129-137.

[4] KANUNGO T, MOUNT D M, NETANYAHU N, et al. A local search approximation algorithm for k-means clustering[J]. Computational Geometry, 2004, 28∶ 89-112.

[5] MATOUSEK J. On approximate geometric k-clustering[J]. Discrete and Computational Geometry, 2000, 24∶ 61-84.

[6] SONG M J, RAJASEKARAN S. Fast k-means algorithms with constant approximation[A]. Proceedings of the 16th Annual International Symposium on Algorithms and Computation[C]. Sanya, Hainan,China, 2005,1029-1038.

[7] INABA M, KAOTH N, IMAI H. Application of weighted voronoi diagrams and randomization to variance-based k-clustering(extended abstract)[A]. Proceedings of the tenth annual symposium on Computational Geometry[C]. Stony Brook, New York, USA, 1994. 332-339.

[8] KUMAR A, SABHARWAL Y, SEN S. A sample linear time (1+ε)algorithm for k-means clustering in any dimensions[A]. Proceedings of the 45th IEEE Symposium on the Foundations of Computer science[C].Washington, DC, USA, 2004. 454-462.

[9] MOTAWANI R, RAGHAVAN P. Randomized Algorithms[M]. Cambridge University Press, Cambridge, UK, 1995.

[10] RUSPINI E H. Numerical methods for fuzzy clustering[J]. Inform Science,1970, 2(3)∶319-350.

[11] SPAETH H. Cluster Analysis Algorithms for Data Reduction and Classification of Objects[M]. John Wiley & Sones, 1980.

[12] PENG J M, XIA Y. A New Theoretical Framework for k-Means-Type Clustering[R]. McMaster University, Advanced Optimization Laboratory, Tech Rep∶ ADVOL2004/06, 2004.