任意の距離関数で k-means 法を実行する
任意の距離関数を用いて k-means 法によるクラスタリングを行うプログラムを作成しました。コードは以下のとおりです。MATLAB で実装しました。fminunc 関数を利用しているため、実行には MATLAB 本体のほかに Optimization Toolbox が必要です。
function [idx, C] = mykmeans(X, k, distance) % データ数 (X は n 行 p 列の行列。n 個の p 次元データを表す) n = size(X, 1); % データからランダムに k 個を選択する方法でセントロイドを初期化する C = X(randperm(n), :); C = C(1:k, :); % fminunc のオプション設定 options = optimoptions('fminunc', 'Algorithm', 'quasi-newton', 'Display', 'off'); % 前回の割り当てを保持する変数 (収束判定用) previdx = zeros(n, 1); % 最大反復回数を 100 回として k-means を実行する for iteration = 1:100 % 各データを最も近いクラスタに割り当てる (E-step) d = zeros(n, k); for ci = 1:k d(:, ci) = arrayfun(@(xi) distance(X(xi, :), C(ci, :)), 1:n); end [o, idx] = min(d, [], 2); % 割り当てが前回と変わらなければ終了する if all(idx == previdx) return; end previdx = idx; % セントロイドをクラスタ内距離の総和が最小になる位置に設定する (M-step) for ci = 1:k C(ci, :) = fminunc(@(c) sumofdistance(X(find(idx == ci), :), c, distance), C(ci, :), options); end end end % X の各点から C への距離の総和を計算する関数 function d = sumofdistance(X, C, distance) d = sum(arrayfun(@(xi) distance(X(xi, :), C), 1:size(X, 1))); end
K-means 法は、クラスタリングのアルゴリズムの一つです。上記のコードにコメントで示したように、次の処理を交互に繰り返すことでデータをクラスタリングします。
後半の処理で、関数の最小化問題を解く必要があります。通常使われるユークリッド距離では、クラスタに割り当てられたデータの中心をセントロイドとしたときに距離の総和が最小になるので、最小化問題の解が簡単に求まります。この性質はユークリッド距離以外では一般に成り立たないので、他の距離関数を使う場合には、距離関数に応じて適切な計算を行わなければいけません。
作成したプログラムは、非線形の最小化問題を解くソルバー (fminunc 関数) を利用してセントロイドを求めます。したがって、距離関数に合わせた最小化問題を自分で解かなくても k-means クラスタリングを実行できます。
実行例として、緯度経度で表されたデータのクラスタリングを試してみます*1。以下のような 16 個のデータを作成します。
latlng = [ ... -85, -135; -85, -45; -85, 45; -85, 135; ... # 南極付近 85, -135; 85, -45; 85, 45; 85, 135; ... # 北極付近 -10, -10; -10, 10; 10, -10; 10, 10; ... # 赤道, 経度 0 度付近 -10, -170; -10, 170; 10, -170; 10, 170 ... # 赤道, 経度 180 度付近 ];
MATLAB でプロットしてみます。分かりやすいように陸地を重ねてみました。赤色の点が 16 個のデータの位置を表します。
figure; geoshow('landareas.shp', 'FaceColor', [0.5 1.0 0.5]); geoshow(latlng(:, 1), latlng(:, 2), ... 'DisplayType', 'point', 'Marker', 'o', 'MarkerEdgeColor', 'red', 'MarkerFaceColor', 'red');
まず、MATLAB の kmeans 関数でユークリッド距離を用いたクラスタリングを実行します。クラスタ数を 4 として kmeans を実行した結果をプロットします。十字の印は各クラスタのセントロイドです。緯度経度をユークリッド平面上の座標と見なしているので、図のような結果になってしまいました*2。
[idx, C] = kmeans(latlng, 4); figure; geoshow('landareas.shp', 'FaceColor', [0.5 1.0 0.5]); colors = { 'red', 'blue', 'magenta', 'black' }; for i = 1:4 X = latlng(find(idx == i), :); geoshow(X(:, 1), X(:, 2), ... 'DisplayType', 'point', 'Marker', 'o', 'MarkerEdgeColor', colors{i}, 'MarkerFaceColor', colors{i}); geoshow(C(i, 1), C(i, 2), ... 'DisplayType', 'point', 'Marker', '+', 'MarkerEdgeColor', colors{i}, 'MarkerFaceColor', colors{i}); end
次に、冒頭のプログラムでクラスタリングを実行します。緯度経度で表された二点間の距離は、Mapping Toolbox に含まれる distance 関数で求められます。先ほどの実行例で kmeans 関数を呼び出していたところを、次のように変更します。図のように、期待した結果が得られました。
% distance(lat1, lng1, lat2, lng2) で二点間の距離を計算する [idx, C] = mykmeans(latlng, 4, @(a, b) distance(a(1), a(2), b(1), b(2)));
余談ですが、セントロイドの初期配置によっては、図の結果が得られないこともあります。今回は簡単のためデータからランダムに選ぶ方法で実装しましたが、100 回の試行で、図のような分け方になった回数は 36 回でした。次のようなコードで k-means++*3 を実装したところ、100 回の試行中 66 回で図の結果が得られました。
function C = kmeanspp(X, k, distance) n =size(X, 1); % 最初のクラスタはデータからランダムに選択する C(1, :) = X(randi(k), :); for i = 2:k % 各データから (これまでに決まった) 最も近いクラスタへの距離を計算する d = zeros(n, i-1); for ci = 1:i-1 d(:, ci) = arrayfun(@(xi) distance(X(xi, :), C(ci, :)), 1:n); end d = min(d, [], 2); % 距離に比例する確率でデータから次のセントロイドを選択する % (既に選ばれたデータは距離 0 になるので、わざわざ除外する必要はない) d = cumsum(d); r = rand * d(end); C(i, :) = X(sum(r > d) + 1, :); end end
*1:この実行例を試すには Mapping Toolbox が必要です。
*2:セントロイドの初期配置によって結果は変わります。南極と北極のクラスタが、経度 0 度付近の点を二つずつ取り合う結果になることもありました。
*3:David Arthur and Sergei Vassilvitskii. 2007. k-means++: the advantages of careful seeding. In Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms (SODA ‘07). Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1027-1035.