医療職からデータサイエンティストへ

統計学、機械学習に関する記事をまとめています。

k-meansとk-means++を視覚的に理解する~Pythonにてスクラッチから~

k-means(k平均法)は教師なし学習の中でもとても有名なアルゴリズムの一つです。例えば、顧客のデータから顧客を購買傾向によってグループ分けしたり、商品の特性からいくつかのグループに分けたりと使用法は様々です。

そんなk-measですが、実は中学生でも知っている点と点の間の距離を使うだけのアルゴリズムで成り立っているので、簡単に実装することができます。

今回は、k-means とk-meansの弱点を克服したk-means++をPythonを使って実装していきます。

k-meansの仕組み

今回は2つの変数からサンプルをクラスタリングすることを想定してみましょう。グラフはこんな感じになります。

import numpy as np
import matplotlib.pyplot as plt
% matplotlib inline

#データの生成
np.random.seed(123)
x1 = np.r_[np.random.normal(size=20,loc=1,scale=2),np.random.normal(size=20,loc=8,scale=2),np.random.normal(size=20,loc=15,scale=2),np.random.normal(size=20,loc=25,scale=2)]
x2 = np.r_[np.random.normal(size=20,loc=15,scale=2),np.random.normal(size=20,loc=1,scale=2),np.random.normal(size=20,loc=20,scale=2),np.random.normal(size=20,loc=0,scale=2)]
X = np.c_[x1,x2]

#可視化
plt.scatter(X[:,0],X[:,1],c="black",s=10,alpha=0.5)
plt.show()

k-means

みたところ4つのグループに分けられそうですね。

さて、k-meansで最初に行うことは分割するグループの数を決めることです。今回は見た目からグループ数を4としましょう。多次元で可視化できない場合や、グループ境界が曖昧な場合のグループの決め方はまた後程書いていきますのでご安心ください。

K=4

グループ数を決めたら、お次はそのグループの中心となる座標を決めます。本来はランダムに4つの点を選んで、それを中心座標とするのですが、今回はここにしましょう。

centers = np.array([[0,5],[5,0],[10,15],[20,10]])
>>array([[ 0,  5],
       [ 5,  0],
       [10, 15],
       [20, 10]])

k-means

そして、すべての点に一番近い中心点のグループのラベルを付けます。

例えば、矢印の点であれば、赤色の中心点が一番近いので、赤ラベルがつきます。この時に中学校でならった2点間の距離の公式を使うのです。

k-means

#Xのサンプル数だけ空のラベルを作る
idx = np.zeros(X.shape[0])

#距離の二乗が一番近い中心点のインデックスを返す。
for i in range(X.shape[0]):
    idx[i] = np.argmin(np.sum((X[i,:] - centers)**2,axis=1))

Pythonのnp.sumにaxis=1をつけることで、行方向に足し算を行って、np.argminで一番小さい値となるインデックスを返します。すべてのラベルを付け終えるとこんな感じになります。

k-means

ここから今度は中心点を移動させます。どこに移動させるかというと、それぞれのグループの平均値をとったところを次の中心点の移動先とするのです。例えば、黄色グループの中心座標はここになります。

X[idx==3,:].mean(axis=0)
>>array([23.87552823,  2.6349639 ])

k-means

このようにすべての中心点を移動させるとこんな感じになります。

#中心の移動
for k in range(k):
    centers[k,:] = X[idx==k,:].mean(axis=0)

k-means

あとは先ほど同様に、中心点を使ってラベルを付け直し、中心点を移動させ、ラベルを付け直し.....というのを中心が移動しなくなるまで繰り返すというのが、k-meansになります。簡単ですね。実際にどのように変化していくかみてみましょう。

idx = np.zeros(X.shape[0])
centers = np.array([[0,5],[5,0],[10,15],[20,10]])

plt.subplot(2, 2, 1)
plt.scatter(X[:,0],X[:,1],c="black",s=10,alpha=0.5)
plt.text(x=22,y=20,s="1",size=15)
plt.scatter(centers[:,0],centers[:,1],color=["r","b","g","orange"])

for j in np.arange(1,4):
    for i in range(X.shape[0]):
        idx[i] = np.argmin(np.sum((X[i,:] - centers)**2,axis=1))
        
    for k in range(K):
        centers[k,:] = X[idx==k,:].mean(axis=0)
        
    plt.subplot(2, 2, j+1)
    plt.scatter(X[idx==0,0],X[idx==0,1],color="r",s=10,alpha=0.5)
    plt.scatter(X[idx==1,0],X[idx==1,1],color="b",s=10,alpha=0.5)
    plt.scatter(X[idx==2,0],X[idx==2,1],color="g",s=10,alpha=0.5)
    plt.scatter(X[idx==3,0],X[idx==3,1],color="orange",s=10,alpha=0.5)
    plt.text(x=22,y=20,s=str(j+1),size=15)
    plt.scatter(centers[:,0],centers[:,1],color=["r","b","g","orange"])

plt.show()

k-means

4回目の更新でほとんどが理想的なクラスに属するようになっているみたいです。ここまでの処理を関数にまとめておきます。

def kmeans(X,K,centers,iter):
    #Xのサンプル数だけ空のラベルを作る
    idx = np.zeros(X.shape[0])
    for _ in range(iter):
        #距離の二乗が一番近い中心点のインデックスを返す。
        for i in range(X.shape[0]):
            idx[i] = np.argmin(np.sum((X[i,:] - centers)**2,axis=1))
        #中心の移動
        for k in range(K):
            centers[k,:] = X[idx==k,:].mean(axis=0)
        
    return idx

引数に渡すのはデータ(X)、グループ数(K)、中心の初期値(centers)、更新回数(iter)です。結果はラベルで返ってきます。

k-meansの問題点とk-means++

一見、とても簡単かつ有用なアルゴリズムにみえるk-meansにも弱点があります。それは、中心の初期値がランダムに選ばれてしまうということです。例えば、こんな少し偏ったところが中心に選ばれたとしましょう。

k-means

すると、何度更新を繰り返してもうまくクラスタリングすることができません。こんな感じになります。

k-means

このk-meansの初期値問題を解消するのがk-means++です。(なぜ+が2つなのかは分かりません...)k-means++では中心の初期値を決める際に、すべての中心点どうしが遠いところに配置されるような確率が高くなるようにします。

例えば、最初の点(青点)が以下の図のようになっていたとすると、2つ目の点(赤点)は青点からの距離が遠いところに配置してほしいですよね。

k-means

このとき、ある点と青点の距離の2乗を青点とすべての点の距離の2乗を足したもので割った値を次の赤点が選ばれる場所の確率とします。(ややこしいので下の図を参照ください)

ここでは一番遠い点が選ばれる確率は36%、一番近い点が選ばれる確率は1.4%となります。同様に3つ目の点は、ある点と赤点・青点の距離の2乗を赤点・青点とすべての点の距離の2乗を足したもので割った値が選ばれる場所の確率となります。

k-means

こうすることで、遠くにある点ほど選ばれる確率が高くなるため、すべての中心点が遠くに配置されやすくなります。あとはk-meansと同様な手順で、グルーピングしていきます。つまり、k-meansとk-means++の違いは、最初の中心点の座標をランダムに決めるか、そうでないかという違いだけです。

実際にk-means++を実装してみましょう。

#距離の2乗を保存するサンプル数 × グループ数の行列をつくる
n = X.shape[0]
distance = np.zeros(n*K).reshape(n,K)

#中心点の座標を保存するための入れ物
centers = np.zeros(2*K).reshape(K,-1)

#最初の確率は均等
pr = np.repeat(1/n,n)
#1つ目の中心点はランダムに選ばれる
centers[0,:] = X[np.random.choice(np.arange(n),1,p=pr),]
distance[:,0] = np.sum((X-centers[0,:])**2,axis=1)

#1つ目の中心点からの距離によって確率を変更
pr = np.sum(distance,axis=1)/np.sum(distance)
#確率に従って2つ目の点を選ぶ
centers[1,:] = X[np.random.choice(np.arange(n),1,p=pr),]
distance[:,1] = np.sum((X-centers[1,:])**2,axis=1)

#以下同様
pr = np.sum(distance,axis=1)/np.sum(distance)
centers[2,:] = X[np.random.choice(np.arange(n),1,p=pr),]
distance[:,2] = np.sum((X-centers[2,:])**2,axis=1)

pr = np.sum(distance,axis=1)/np.sum(distance)
centers[3,:] = X[np.random.choice(np.arange(n),1,p=pr),]
distance[:,3] = np.sum((X-centers[3,:])**2,axis=1)


#可視化
plt.scatter(X[:,0],X[:,1],c="black",s=10,alpha=0.5)
plt.scatter(centers[:,0],centers[:,1],color=["r","b","g","orange"])

k-means

4つの点がうまく離れたところに配置されたようです(点は確率的に選ばれるので、全てが離れた場所に選ばれるとは限りません)。 分かりやすいように分けて書いてありますが、forループで書いても問題ありません。ポイントはnp.sum(distance,axis=1)/np.sum(distance)で、点と点の距離の2乗が遠いほど大きな値になるように計算して、np.random.choiceのpに渡すことで、選ばれやすさに重みづけをしているところにあります。k-means++も関数にしておきましょう。

def kmeansplus(X,K,inter):
    n = X.shape[0]
    distance = np.zeros(n*K).reshape(n,K)
    centers = np.zeros(X.shape[1]*K).reshape(K,-1)

    #最初の確率は均等
    pr = np.repeat(1/n,n)
    #1つ目の中心点はランダムに選ぶ
    centers[0,:] = X[np.random.choice(np.arange(n),1,p=pr),]
    distance[:,0] = np.sum((X-centers[0,:])**2,axis=1)
    
    for k in np.arange(1,K):
        pr = np.sum(distance,axis=1)/np.sum(distance)
        centers[k,:] = X[np.random.choice(np.arange(n),1,p=pr),]
        distance[:,k] = np.sum((X-centers[k,:])**2,axis=1)
    
    idx = kmeans(X,K,centers,inter)
    
    return idx

そして実際にクラスタリングしてみます。

idx = kmeansplus(X,K,10)

plt.scatter(x1[idx==0],x2[idx==0],color="r",s=10,alpha=0.5)
plt.scatter(x1[idx==1],x2[idx==1],color="b",s=10,alpha=0.5)
plt.scatter(x1[idx==2],x2[idx==2],color="g",s=10,alpha=0.5)
plt.scatter(x1[idx==3],x2[idx==3],color="orange",s=10,alpha=0.5)
plt.show()

k-means

確率的に選ばれるので、たまに失敗しますが、全くランダムに中心点が選ばれるよりは、うまくクラスタリングすることができます。


補足:なぜ最大値ではなく確率で選ぶのか?

k-means++は最初の中心点を距離が遠いほど確率的に選ばれやすくするアルゴリズムです。では、なぜ最も遠い点ではなく、確率的に選ぶのでしょうか?それは、外れ値に対応するためです。もし、仮に外れ値があった場合には、どの点からも距離が遠くなるので、必然的に最初の中心点として選ばれてしまいます。これを防ぐためにk-means++では、確率によって遠い中心点を選ぶというアルゴリズムを採用しています。


クラスの数を決めるエルボー法

以上でk-meansとk-measn++の実装は終わりなのですが、ここで高次元のためデータを可視化できなかったり、クラスの境界が曖昧だった場合に適切だと思われるクラスの数を決める方法があります。それがエルボー法です。エルボー法とはその名の通りです。

k-meansはグループの中心点からの距離が最小になるようにクラスタリングされていきます。ということは、そのグループ中心からの距離の合計を、全てのグループで足し合わせた値が小さければ、よいクラスタリングができたと言えそうです。これをクラスタ内誤差平方和(SSE)と呼びます。

しかし、小さければ良いというものでもなく、例えばクラスの数がサンプル数分だけあればSSEは0になります。そこで、クラスの数を増やしていった時にこのSSEがガクっと下がったところ(肘)を最適なクラス数としようというのがエルボー法になります。

SSEを計算するために先ほどのkmeans関数に少し手を加えます。

def kmeans(X,K,centers,iter):
    idx = np.zeros(X.shape[0])
    #クラスごとのsseを格納するための入れ物
    sse = np.zeros(K)
    for _ in range(iter):
        
        for i in range(X.shape[0]):
            idx[i] = np.argmin(np.sum((X[i,:] - centers)**2,axis=1))
       
        for k in range(K):
            centers[k,:] = X[idx==k,:].mean(axis=0)
            sse[k] = np.sum((X[idx==k,:]-centers[k,:])**2)#sseの計算
    sse_sum = np.sum(sse)#全てのクラスで合計する
    return idx,sse_sum

新たにsseを計算する処理を加えました。これでクラスの数を変化させた時にsseがどのようになるかをみてみましょう。クラス数は1から10までを試してみます。

K = 10
sse_vec = np.zeros(K)
for k in range(K):
    idx , sse = kmeansplus(X,k+1,10)
    sse_vec[k] = sse

plt.plot(np.arange(1,11),sse_vec)
plt.plot(np.arange(1,11),sse_vec,"bo")
plt.show()

k-means

最後のクラス10の場合はおそらく初期値選びに失敗していますが、これを見る限りクラス数が4つの時にsseがガクッと肘のようになっています。ここからクラス数4が適切だと判断できるのです。

scikit-learnを使ったk-meansとk-means++

さて、ここまでは学習のために1から実装していきましたが、おそらく実務ではscikit-learnを使うことが多いでしょう。ついでにscikit-learnを使った方法も紹介しておきます。

from sklearn.cluster import KMeans

kmeans_mod = KMeans(n_clusters=4, # クラスター数
            init='k-means++',       # 中心の設定
            n_init=10,               # 異なる初期値を用いたk-meansの実行回数 
            max_iter=10,            # 最大イテレーション回数  default: '300'
            tol=1e-04,              # 収束と判定するための相対的な許容誤差 default: '1e-04'
            random_state=0)   

sklernのkmeansでは色々なパラメーターを設定できます。最低限n_clusterさえ指定すれば、デフォルト設定で動いてくれます。ちなみにデフォルトはkmeans++になります。もしどうしても通常のk-meansがいい!という場合はinitに"random"と指定します。n_initは中心の初期値選びに失敗することを考慮して、指定した回数だけk-meansを実行し、sseが一番小さくなる初期値を選びます(ありがたや)。また、今回は実装しなかったtolはk-meansの中心がどれぐらいの範囲で動かなくなったら終了させるかも指定できます。

モデルを指定したらあとはfit_predictとすれば

idx = kmeans_mod.fit_predict(X)

#可視化
plt.scatter(X[idx==0,0],X[idx==0,1],color="r",s=10,alpha=0.5)
plt.scatter(X[idx==1,0],X[idx==1,1],color="b",s=10,alpha=0.5)
plt.scatter(X[idx==2,0],X[idx==2,1],color="g",s=10,alpha=0.5)
plt.scatter(X[idx==3,0],X[idx==3,1],color="orange",s=10,alpha=0.5)
plt.scatter(kmeans_mod.cluster_centers_[:,0],kmeans_mod.cluster_centers_[:,1],color=["r","b","g","orange"])
plt.show()

k-means

とても簡単ですね!また、sseにアクセスするには

kmeans_mod.inertia_
>>722.0352358740058

とします。

まとめ

k-meansはとてもシンプルなアルゴリズムですが、応用範囲が広く直感的にも分かりやすいのが特徴です。sklernを使えば数行で実行できてしまうのも魅力的ですね!

※本記事は筆者が個人的に学んだことをまとめた記事なります。所属する組織の意見・見解とは無関係です。また、数学の記法や詳細な理論、用語等で誤りがあった際はご指摘頂けると幸いです。

参考

k-meansの最適なクラスター数を調べる方法 - Qiita