Python實現(xiàn)的Kmeans++算法實例

更新時間：2014年04月26日 10:38:34 作者：

這篇文章主要介紹了Kmeans和kmeans++算法,講解了Kmeans算法的缺點和kmeans++算法的實現(xiàn)思路,以及Python和matlab中實現(xiàn)的Kmeans++算法,需要的朋友可以參考下

1、從Kmeans說起

Kmeans是一個非?；A(chǔ)的聚類算法，使用了迭代的思想，關(guān)于其原理這里不說了。下面說一下如何在matlab中使用kmeans算法。

創(chuàng)建7個二維的數(shù)據(jù)點：

復(fù)制代碼代碼如下:

x=[randn(3,2)*.4;randn(4,2)*.5+ones(4,1)*[4 4]];

使用kmeans函數(shù)：

復(fù)制代碼代碼如下:

class = kmeans(x, 2);

x是數(shù)據(jù)點，x的每一行代表一個數(shù)據(jù)；2指定要有2個中心點，也就是聚類結(jié)果要有2個簇。 class將是一個具有70個元素的列向量，這些元素依次對應(yīng)70個數(shù)據(jù)點，元素值代表著其對應(yīng)的數(shù)據(jù)點所處的分類號。某次運行后，class的值是：

復(fù)制代碼代碼如下:

這說明x的前三個數(shù)據(jù)點屬于簇2，而后四個數(shù)據(jù)點屬于簇1。 kmeans函數(shù)也可以像下面這樣使用：

復(fù)制代碼代碼如下:

>> [class, C, sumd, D] = kmeans(x, 2)

class =
     2
     2
     2
     1
     1
     1
     1

C =
4.0629 4.0845
-0.1341 0.1201

sumd =
1.2017
0.2939

D =
   34.3727    0.0184
   29.5644    0.1858
   36.3511    0.0898
    0.1247   37.4801
    0.7537   24.0659
    0.1979   36.7666
    0.1256   36.2149

class依舊代表著每個數(shù)據(jù)點的分類;C包含最終的中心點，一行代表一個中心點；sumd代表著每個中心點與所屬簇內(nèi)各個數(shù)據(jù)點的距離之和；D的每一行也對應(yīng)一個數(shù)據(jù)點，行中的數(shù)值依次是該數(shù)據(jù)點與各個中心點之間的距離，Kmeans默認使用的距離是歐幾里得距離（參考資料[3]）的平方值。kmeans函數(shù)使用的距離，也可以是曼哈頓距離（L1-距離），以及其他類型的距離，可以通過添加參數(shù)指定。

kmeans有幾個缺點（這在很多資料上都有說明）：

1、最終簇的類別數(shù)目（即中心點或者說種子點的數(shù)目）k并不一定能事先知道，所以如何選一個合適的k的值是一個問題。
2、最開始的種子點的選擇的好壞會影響到聚類結(jié)果。
3、對噪聲和離群點敏感。
4、等等。

2、kmeans++算法的基本思路

kmeans++算法的主要工作體現(xiàn)在種子點的選擇上，基本原則是使得各個種子點之間的距離盡可能的大，但是又得排除噪聲的影響。以下為基本思路：

1、從輸入的數(shù)據(jù)點集合（要求有k個聚類）中隨機選擇一個點作為第一個聚類中心
2、對于數(shù)據(jù)集中的每一個點x，計算它與最近聚類中心(指已選擇的聚類中心)的距離D(x)
3、選擇一個新的數(shù)據(jù)點作為新的聚類中心，選擇的原則是：D(x)較大的點，被選取作為聚類中心的概率較大
4、重復(fù)2和3直到k個聚類中心被選出來
5、利用這k個初始的聚類中心來運行標準的k-means算法

假定數(shù)據(jù)點集合X有n個數(shù)據(jù)點，依次用X(1)、X(2)、……、X(n)表示，那么，在第2步中依次計算每個數(shù)據(jù)點與最近的種子點（聚類中心）的距離，依次得到D(1)、D(2)、……、D(n)構(gòu)成的集合D。在D中，為了避免噪聲，不能直接選取值最大的元素，應(yīng)該選擇值較大的元素，然后將其對應(yīng)的數(shù)據(jù)點作為種子點。

如何選擇值較大的元素呢，下面是一種思路（暫未找到最初的來源，在資料[2]等地方均有提及，筆者換了一種讓自己更好理解的說法）：把集合D中的每個元素D(x)想象為一根線L(x)，線的長度就是元素的值。將這些線依次按照L(1)、L(2)、……、L(n)的順序連接起來，組成長線L。L(1)、L(2)、……、L(n)稱為L的子線。根據(jù)概率的相關(guān)知識，如果我們在L上隨機選擇一個點，那么這個點所在的子線很有可能是比較長的子線，而這個子線對應(yīng)的數(shù)據(jù)點就可以作為種子點。下文中kmeans++的兩種實現(xiàn)均是這個原理。

3、python版本的kmeans++

在http://rosettacode.org/wiki/K-means%2B%2B_clustering 中能找到多種編程語言版本的Kmeans++實現(xiàn)。下面的內(nèi)容是基于python的實現(xiàn)（中文注釋是筆者添加的）：

復(fù)制代碼代碼如下:

from math import pi, sin, cos
from collections import namedtuple
from random import random, choice
from copy import copy

try:
    import psyco
    psyco.full()
except ImportError:
    pass

FLOAT_MAX = 1e100

class Point:
    __slots__ = ["x", "y", "group"]
    def __init__(self, x=0.0, y=0.0, group=0):
        self.x, self.y, self.group = x, y, group

def generate_points(npoints, radius):
points = [Point() for _ in xrange(npoints)]

    # note: this is not a uniform 2-d distribution
    for p in points:
        r = random() * radius
        ang = random() * 2 * pi
        p.x = r * cos(ang)
        p.y = r * sin(ang)

return points

def nearest_cluster_center(point, cluster_centers):
    """Distance and index of the closest cluster center"""
    def sqr_distance_2D(a, b):
        return (a.x - b.x) ** 2 + (a.y - b.y) ** 2

min_index = point.group
min_dist = FLOAT_MAX

    for i, cc in enumerate(cluster_centers):
        d = sqr_distance_2D(cc, point)
        if min_dist > d:
            min_dist = d
            min_index = i

return (min_index, min_dist)

'''
points是數(shù)據(jù)點，nclusters是給定的簇類數(shù)目
cluster_centers包含初始化的nclusters個中心點，開始都是對象->(0,0,0)
'''

def kpp(points, cluster_centers):
cluster_centers[0] = copy(choice(points)) #隨機選取第一個中心點
d = [0.0 for _ in xrange(len(points))] #列表，長度為len(points)，保存每個點離最近的中心點的距離

    for i in xrange(1, len(cluster_centers)): # i=1...len(c_c)-1
        sum = 0
        for j, p in enumerate(points):
            d[j] = nearest_cluster_center(p, cluster_centers[:i])[1] #第j個數(shù)據(jù)點p與各個中心點距離的最小值
            sum += d[j]

sum *= random()

        for j, di in enumerate(d):
            sum -= di
            if sum > 0:
                continue
            cluster_centers[i] = copy(points[j])
            break

for p in points:
p.group = nearest_cluster_center(p, cluster_centers)[0]

'''
points是數(shù)據(jù)點，nclusters是給定的簇類數(shù)目
'''
def lloyd(points, nclusters):
cluster_centers = [Point() for _ in xrange(nclusters)] #根據(jù)指定的中心點個數(shù)，初始化中心點，均為(0,0,0)

# call k++ init
kpp(points, cluster_centers) #選擇初始種子點

# 下面是kmeans
lenpts10 = len(points) >> 10

    changed = 0
    while True:
        # group element for centroids are used as counters
        for cc in cluster_centers:
            cc.x = 0
            cc.y = 0
            cc.group = 0

        for p in points:
            cluster_centers[p.group].group += 1 #與該種子點在同一簇的數(shù)據(jù)點的個數(shù)
            cluster_centers[p.group].x += p.x
            cluster_centers[p.group].y += p.y

        for cc in cluster_centers:    #生成新的中心點
            cc.x /= cc.group
            cc.y /= cc.group

        # find closest centroid of each PointPtr
        changed = 0 #記錄所屬簇發(fā)生變化的數(shù)據(jù)點的個數(shù)
        for p in points:
            min_i = nearest_cluster_center(p, cluster_centers)[0]
            if min_i != p.group:
                changed += 1
                p.group = min_i

        # stop when 99.9% of points are good
        if changed <= lenpts10:
            break

for i, cc in enumerate(cluster_centers):
cc.group = i

return cluster_centers

def print_eps(points, cluster_centers, W=400, H=400):
Color = namedtuple("Color", "r g b");

    colors = []
    for i in xrange(len(cluster_centers)):
        colors.append(Color((3 * (i + 1) % 11) / 11.0,
                            (7 * i % 11) / 11.0,
                            (9 * i % 11) / 11.0))

max_x = max_y = -FLOAT_MAX
min_x = min_y = FLOAT_MAX

    for p in points:
        if max_x < p.x: max_x = p.x
        if min_x > p.x: min_x = p.x
        if max_y < p.y: max_y = p.y
        if min_y > p.y: min_y = p.y

    scale = min(W / (max_x - min_x),
                H / (max_y - min_y))
    cx = (max_x + min_x) / 2
    cy = (max_y + min_y) / 2

print "%%!PS-Adobe-3.0\n%%%%BoundingBox: -5 -5 %d %d" % (W + 10, H + 10)

    print ("/l {rlineto} def /m {rmoveto} def\n" +
           "/c { .25 sub exch .25 sub exch .5 0 360 arc fill } def\n" +
           "/s { moveto -2 0 m 2 2 l 2 -2 l -2 -2 l closepath " +
           "   gsave 1 setgray fill grestore gsave 3 setlinewidth" +
           " 1 setgray stroke grestore 0 setgray stroke }def")

    for i, cc in enumerate(cluster_centers):
        print ("%g %g %g setrgbcolor" %
               (colors[i].r, colors[i].g, colors[i].b))

        for p in points:
            if p.group != i:
                continue
            print ("%.3f %.3f c" % ((p.x - cx) * scale + W / 2,
                                    (p.y - cy) * scale + H / 2))

print ("\n0 setgray %g %g s" % ((cc.x - cx) * scale + W / 2,
(cc.y - cy) * scale + H / 2))

print "\n%%%%EOF"

def main():
npoints = 30000
k = 7 # # clusters

    points = generate_points(npoints, 10)
    cluster_centers = lloyd(points, k)
    print_eps(points, cluster_centers)

main()

上述代碼實現(xiàn)的算法是針對二維數(shù)據(jù)的，所以Point對象有三個屬性，分別是在x軸上的值、在y軸上的值、以及所屬的簇的標識。函數(shù)lloyd是kmeans++算法的整體實現(xiàn)，其先是通過kpp函數(shù)選取合適的種子點，然后對數(shù)據(jù)集實行kmeans算法進行聚類。kpp函數(shù)的實現(xiàn)完全符合上述kmeans++的基本思路的2、3、4步。

4、matlab版本的kmeans++

復(fù)制代碼代碼如下:

function [L,C] = kmeanspp(X,k)
%KMEANS Cluster multivariate data using the k-means++ algorithm.
%   [L,C] = kmeans_pp(X,k) produces a 1-by-size(X,2) vector L with one class
%   label per column in X and a size(X,1)-by-k matrix C containing the
%   centers corresponding to each class.

% Version: 2013-02-08
% Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)

L = [];
L1 = 0;

while length(unique(L)) ~= k

    % The k-means++ initialization.
    C = X(:,1+round(rand*(size(X,2)-1))); %size(X,2)是數(shù)據(jù)集合X的數(shù)據(jù)點的數(shù)目，C是中心點的集合
    L = ones(1,size(X,2));
    for i = 2:k
        D = X-C(:,L); %-1
        D = cumsum(sqrt(dot(D,D,1))); %將每個數(shù)據(jù)點與中心點的距離，依次累加
        if D(end) == 0, C(:,i:k) = X(:,ones(1,k-i+1)); return; end
        C(:,i) = X(:,find(rand < D/D(end),1)); %find的第二個參數(shù)表示返回的索引的數(shù)目
        [~,L] = max(bsxfun(@minus,2*real(C'*X),dot(C,C,1).')); %碉堡了，這句，將每個數(shù)據(jù)點進行分類。
    end

    % The k-means algorithm.
    while any(L ~= L1)
        L1 = L;
        for i = 1:k, l = L==i; C(:,i) = sum(X(:,l),2)/sum(l); end
        [~,L] = max(bsxfun(@minus,2*real(C'*X),dot(C,C,1).'),[],1);
    end

end

這個函數(shù)的實現(xiàn)有些特殊，參數(shù)X是數(shù)據(jù)集，但是是將每一列看做一個數(shù)據(jù)點，參數(shù)k是指定的聚類數(shù)。返回值L標記了每個數(shù)據(jù)點的所屬分類，返回值C保存了最終形成的中心點（一列代表一個中心點）。測試一下：

復(fù)制代碼代碼如下:

>> x=[randn(3,2)*.4;randn(4,2)*.5+ones(4,1)*[4 4]]
x =
   -0.0497    0.5669
    0.5959    0.2686
    0.5636   -0.4830
    4.3586    4.3634
    4.8151    3.8483
    4.2444    4.1469
    4.5173    3.6064

>> [L, C] = kmeanspp(x',2)
L =
     2     2     2     1     1     1     1
C =
    4.4839    0.3699
    3.9913    0.1175

好了，現(xiàn)在開始一點點理解這個實現(xiàn)，順便鞏固一下matlab知識。

unique函數(shù)用來獲取一個矩陣中的不同的值，示例：

復(fù)制代碼代碼如下:

>> unique([1 3 3 4 4 5])
ans =
     1     3     4     5
>> unique([1 3 3 ; 4 4 5])
ans =
     1
     3
     4
     5

所以循環(huán) while length(unique(L)) ~= k 以得到了k個聚類為結(jié)束條件，不過一般情況下，這個循環(huán)一次就結(jié)束了，因為重點在這個循環(huán)中。

rand是返回在(0,1)這個區(qū)間的一個隨機數(shù)。在注釋%-1所在行，C被擴充了，被擴充的方法類似于下面：

復(fù)制代碼代碼如下:

>> C =[];
>> C(1,1) = 1
C =
     1
>> C(2,1) = 2
C =
     1
     2
>> C(:,[1 1 1 1])
ans =
     1     1     1     1
     2     2     2     2
>> C(:,[1 1 1 1 2])
Index exceeds matrix dimensions.

C中第二個參數(shù)的元素1，其實是代表C的第一列數(shù)據(jù)，之所以在值2時候出現(xiàn)Index exceeds matrix dimensions.的錯誤，是因為C本身沒有第二列。如果C有第二列了：

復(fù)制代碼代碼如下:

>> C(2,2) = 3;
>> C(2,2) = 4;
>> C(:,[1 1 1 1 2])
ans =
     1     1     1     1     3
     2     2     2     2     4

dot函數(shù)是將兩個矩陣點乘，然后把結(jié)果在某一維度相加：

復(fù)制代碼代碼如下:

>> TT = [1 2 3 ; 4 5 6];
>> dot(TT,TT)
ans =
    17    29    45
>> dot(TT,TT,1 )
ans =
    17    29    45

<code>cumsum</code>是累加函數(shù)：

復(fù)制代碼代碼如下:

>> cumsum([1 2 3])
ans =
     1     3     6
>> cumsum([1 2 3; 4 5 6])
ans =
     1     2     3
     5     7     9

max函數(shù)可以返回兩個值，第二個代表的是max數(shù)的索引位置：

復(fù)制代碼代碼如下:

>> [~, L] = max([1 2 3])
L =
     3
>> [~,L] = max([1 2 3;2 3 4])
L =
     2     2     2

其中~是占位符。

關(guān)于bsxfun函數(shù)，官方文檔指出：

復(fù)制代碼代碼如下:

C = bsxfun(fun,A,B) applies the element-by-element binary operation specified by the function handle fun to arrays A and B, with singleton expansion enabled

其中參數(shù)fun是函數(shù)句柄，關(guān)于函數(shù)句柄見資料[9]。下面是bsxfun的一個示例：

復(fù)制代碼代碼如下:

>> A= [1 2 3;2 3 4]
A =
     1     2     3
     2     3     4
>> B=[6;7]
B =
     6
     7
>> bsxfun(@minus,A,B)
ans =
    -5    -4    -3
    -5    -4    -3

對于：

復(fù)制代碼代碼如下:

[~,L] = max(bsxfun(@minus,2*real(C'*X),dot(C,C,1).'));

max的參數(shù)是這樣一個矩陣，矩陣有n列，n也是數(shù)據(jù)點的個數(shù)，每一列代表著對應(yīng)的數(shù)據(jù)點與各個中心點之間的距離的相反數(shù)。不過這個距離有些與眾不同，算是歐幾里得距離的變形。

假定數(shù)據(jù)點是2維的，某個數(shù)據(jù)點為(x1,y1)，某個中心點為(c1,d1)，那么通過bsxfun(@minus,2real(C'X),dot(C,C,1).')的計算，數(shù)據(jù)點與中心點的距離為2c1x1 + 2d1y1 -c1.^2 - c2.^2，可以變換為x1.^2 + y1.^2 - (c1-x1).^2 - (d1-y1).^2。對于每一列而言，由于是數(shù)據(jù)點與各個中心點之間的計算，所以可以忽略x1.^2 + y1.^2，最終計算結(jié)果是歐幾里得距離的平方的相反數(shù)。這也說明了使用max的合理性，因為一個數(shù)據(jù)點的所屬簇取決于與其距離最近的中心點，若將距離取相反數(shù)，則應(yīng)該是值最大的那個點。

您可能感興趣的文章: