1、從Kmeans說起
Kmeans是一個非常基礎的聚類算法，使用了迭代的思想，關于其原理這里不說了。下面說一下如何在matlab中使用kmeans算法。
創建7個二維的數據點：
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依舊代表著每個數據點的分類;C包含最終的中心點，一行代表一個中心點；sumd代表著每個中心點與所屬簇內各個數據點的距離之和；D的每一行也對應一個數據點，行中的數值依次是該數據點與各個中心點之間的距離，Kmeans默認使用的距離是歐幾里得距離（參考資料[3]）的平方值。kmeans函數使用的距離，也可以是曼哈頓距離（L1-距離），以及其他類型的距離，可以通過添加參數指定。
kmeans有幾個缺點（這在很多資料上都有說明）：
1、最終簇的類別數目（即中心點或者說種子點的數目）k并不一定能事先知道，所以如何選一個合適的k的值是一個問題。
2、最開始的種子點的選擇的好壞會影響到聚類結果。
3、對噪聲和離群點敏感。
4、等等。
2、kmeans++算法的基本思路
kmeans++算法的主要工作體現在種子點的選擇上，基本原則是使得各個種子點之間的距離盡可能的大，但是又得排除噪聲的影響。 以下為基本思路：
1、從輸入的數據點集合（要求有k個聚類）中隨機選擇一個點作為第一個聚類中心
2、對于數據集中的每一個點x，計算它與最近聚類中心(指已選擇的聚類中心)的距離D(x)
3、選擇一個新的數據點作為新的聚類中心，選擇的原則是：D(x)較大的點，被選取作為聚類中心的概率較大
4、重復2和3直到k個聚類中心被選出來
5、利用這k個初始的聚類中心來運行標準的k-means算法
假定數據點集合X有n個數據點，依次用X(1)、X(2)、……、X(n)表示，那么，在第2步中依次計算每個數據點與最近的種子點（聚類中心）的距離，依次得到D(1)、D(2)、……、D(n)構成的集合D。在D中，為了避免噪聲，不能直接選取值最大的元素，應該選擇值較大的元素，然后將其對應的數據點作為種子點。
如何選擇值較大的元素呢，下面是一種思路（暫未找到最初的來源，在資料[2]等地方均有提及，筆者換了一種讓自己更好理解的說法）： 把集合D中的每個元素D(x)想象為一根線L(x)，線的長度就是元素的值。將這些線依次按照L(1)、L(2)、……、L(n)的順序連接起來，組成長線L。L(1)、L(2)、……、L(n)稱為L的子線。根據概率的相關知識，如果我們在L上隨機選擇一個點，那么這個點所在的子線很有可能是比較長的子線，而這個子線對應的數據點就可以作為種子點。下文中kmeans++的兩種實現均是這個原理。
3、python版本的kmeans++
在http://rosettacode.org/wiki/K-means%2B%2B_clustering 中能找到多種編程語言版本的Kmeans++實現。下面的內容是基于python的實現（中文注釋是筆者添加的）：
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是數據點，nclusters是給定的簇類數目
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個數據點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是數據點，nclusters是給定的簇類數目
'''
def lloyd(points, nclusters):
cluster_centers = [Point() for _ in xrange(nclusters)] #根據指定的中心點個數，初始化中心點，均為(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 #與該種子點在同一簇的數據點的個數
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 #記錄所屬簇發生變化的數據點的個數
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()
上述代碼實現的算法是針對二維數據的，所以Point對象有三個屬性，分別是在x軸上的值、在y軸上的值、以及所屬的簇的標識。函數lloyd是kmeans++算法的整體實現，其先是通過kpp函數選取合適的種子點，然后對數據集實行kmeans算法進行聚類。kpp函數的實現完全符合上述kmeans++的基本思路的2、3、4步。
4、matlab版本的kmeans++
% 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)是數據集合X的數據點的數目，C是中心點的集合
L = ones(1,size(X,2));
for i = 2:k
D = X-C(:,L); %-1
D = cumsum(sqrt(dot(D,D,1))); %將每個數據點與中心點的距離，依次累加
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的第二個參數表示返回的索引的數目
[~,L] = max(bsxfun(@minus,2*real(C'*X),dot(C,C,1).')); %碉堡了，這句，將每個數據點進行分類。
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
這個函數的實現有些特殊，參數X是數據集，但是是將每一列看做一個數據點，參數k是指定的聚類數。返回值L標記了每個數據點的所屬分類，返回值C保存了最終形成的中心點（一列代表一個中心點）。測試一下：
>> [L, C] = kmeanspp(x',2)
L =
2 2 2 1 1 1 1
C =
4.4839 0.3699
3.9913 0.1175
好了，現在開始一點點理解這個實現，順便鞏固一下matlab知識。
unique函數用來獲取一個矩陣中的不同的值，示例：
rand是返回在(0,1)這個區間的一個隨機數。在注釋%-1所在行，C被擴充了，被擴充的方法類似于下面：
C中第二個參數的元素1，其實是代表C的第一列數據，之所以在值2時候出現Index exceeds matrix dimensions.的錯誤，是因為C本身沒有第二列。如果C有第二列了：
關于bsxfun函數，官方文檔指出：
假定數據點是2維的，某個數據點為(x1,y1)，某個中心點為(c1,d1)，那么通過bsxfun(@minus,2real(C'X),dot(C,C,1).')的計算，數據點與中心點的距離為2c1x1 + 2d1y1 -c1.^2 - c2.^2，可以變換為x1.^2 + y1.^2 - (c1-x1).^2 - (d1-y1).^2。對于每一列而言，由于是數據點與各個中心點之間的計算，所以可以忽略x1.^2 + y1.^2，最終計算結果是歐幾里得距離的平方的相反數。這也說明了使用max的合理性，因為一個數據點的所屬簇取決于與其距離最近的中心點，若將距離取相反數，則應該是值最大的那個點。
新聞熱點
疑難解答
圖片精選
