#はじめに
非階層クラスタリングの一つにxmeansという手法があります。xmeansとはkmeansを拡張した手法であり、クラスタ数をBIC(ベイズ情報量規準)を用いて自動的に決定してくれます。今回は会社での業務を通じて、xmeansの弱点を発見したので記事にしようと思います。
なお、xmeansについては以下の記事を参考にしました。
#Xmeansの弱点
##xmeans実行
xmeansを用いたときに生じた問題を再現するため、以下のようにシミュレーションデータを作成します。ポイントは 60個の要素をもつベクトル $z$ の最初の20個が $0$ であることです。
#基礎計算系のライブラリ
import numpy as np
from scipy import stats
#クラスタリングに用いるライブラリ
from sklearn.cluster import KMeans
%matplotlib inline
#シミュレーションデータ作成
np.random.seed(0)
x = np.array([np.random.normal(loc, 0.3, 20) for loc in [1, 1, -1]]).flatten() #ランダムな60個の数を生成
y = np.array([np.random.normal(loc, 0.3, 20) for loc in [1, 1, -1]]).flatten() #ランダムな60個の数を生成
z = np.r_[np.zeros(20), np.array([np.random.normal(loc, 0.3, 20) for loc in [1, -1]]).flatten()] #0を20個とランダムな40個の数を生成
data = np.c_[x, y, z]
作成したデータの散布図は下記です。
このdataに対して、@deaikeiさんが作成したスクリプト(クラスタ数を自動推定するX-means法を調べてみた)を用いてxmeansを実行します。
(この記事わかりやすくて本当にすごいですね。Respectです！！)
#xmeansのクラス
class XMeans:
"""
x-means法を行うクラス
"""
def __init__(self, k_init = 2, **k_means_args):
"""
k_init : The initial number of clusters applied to KMeans()
"""
self.k_init = k_init
self.k_means_args = k_means_args
def fit(self, X):
"""
x-means法を使ってデータXをクラスタリングする
X : array-like or sparse matrix, shape=(n_samples, n_features)
"""
self.__clusters = [] #clustersにはもうこれ以上分割しないデータが格納される
clusters = self.Cluster.build(X, KMeans(self.k_init, **self.k_means_args).fit(X))
self.__recursively_split(clusters)
self.labels_ = np.empty(X.shape[0], dtype = np.intp)
for i, c in enumerate(self.__clusters):
self.labels_[c.index] = i
self.cluster_centers_ = np.array([c.center for c in self.__clusters])
self.cluster_log_likelihoods_ = np.array([c.log_likelihood() for c in self.__clusters])
self.cluster_sizes_ = np.array([c.size for c in self.__clusters])
return self
def __recursively_split(self, clusters):
"""
引数のclustersを再帰的に分割する
clusters : list-like object, which contains instances of 'XMeans.Cluster'
'XMeans.Cluster'のインスタンスを含むリスト型オブジェクト
"""
for cluster in clusters:
if cluster.size <= 3:
self.__clusters.append(cluster) #仮説：クラスタに含まれるデータの個数が3個以内の場合はこれ以上分割しない
continue #cluster.size<=3だったら次のループへいく。そうではなかったら次の処理へ進む
k_means = KMeans(2, **self.k_means_args).fit(cluster.data)
c1, c2 = self.Cluster.build(cluster.data, k_means, cluster.index)
beta = np.linalg.norm(c1.center - c2.center) / np.sqrt(np.linalg.det(c1.cov) + np.linalg.det(c2.cov))
alpha = 0.5 / stats.norm.cdf(beta)
bic = -2 * (cluster.size * np.log(alpha) + c1.log_likelihood() + c2.log_likelihood()) + 2 * cluster.df * np.log(cluster.size)
if bic < cluster.bic():
self.__recursively_split([c1, c2]) #更なる分割を行う
else:
self.__clusters.append(cluster) #これ以上分割しない
class Cluster:
"""
k-means法によって生成されたクラスタに関する情報を持ち、尤度やBICの計算を行うクラス
"""
@classmethod
def build(cls, X, k_means, index = None):
if index is None:
index = np.array(range(0, X.shape[0]))
labels = range(0, k_means.get_params()["n_clusters"])
return tuple(cls(X, index, k_means, label) for label in labels) #ここでクラスター結果を入れたインスタンスを生成し、クラスタ数分の「インスタンスを返す
# index: Xの各行におけるサンプルが元データの何行目のものかを示すベクトル
def __init__(self, X, index, k_means, label):
self.data = X[k_means.labels_ == label]
self.index = index[k_means.labels_ == label]
self.size = self.data.shape[0]
self.df = self.data.shape[1] * (self.data.shape[1] + 3) / 2
self.center = k_means.cluster_centers_[label]
self.cov = np.cov(self.data.T)
def log_likelihood(self):
return sum(stats.multivariate_normal.logpdf(x, self.center, self.cov) for x in self.data)
def bic(self):
return -2 * self.log_likelihood() + self.df * np.log(self.size)
#x-means実行!
x_means = XMeans(random_state = 1).fit(data)
実行結果は...
LinAlgError: singular matrix
というエラーを吐いてしまいました。
##原因
このスクリプトがエラーを吐いた理由は、xmeansがBICを用いて計算しているところにあります。BICを各クラスに対して計算するときに、各クラスに含まれるデータの共分散を計算し、その逆行列を求める必要があります。今回は、 $z$ に大量の $0$ が含まれていたために分割後のクラスタについて、共分散の逆行列を計算できなくなり、エラーを吐いてしまいました。なお、@deaikeiさんが作成したスクリプトのBICの計算方法は、
クラスター数を自動決定するk-meansアルゴリズムの拡張について
を参考にされたのではないかと思いました。(違っていたらごめんなさい！)
つまり、xmeansはある次元が説明力を持たない(=ある次元について一定の値が並ぶ)ときにBICを計算できなくなり、実行不可能になるということがわかりました。
#pyclusteringを用いたXmeans
##xmeans実行
さて、xmeansの限界がわかったところでpythonの一般公開されているライブラリ(pyclustering)を用いて、同様のシチュエーションでxmeansを実行するとどうなるのか試してみました。
ちなみにpyclusteringの公式ドキュメントはこちらです⇒pyclustering.cluster.xmeans.xmeans Class Reference
#クラスタリングに用いるライブラリ
from sklearn.cluster import KMeans
from pyclustering.cluster.xmeans import xmeans
#pyclusteringを用いたx-meansの実行
x_means_pyclustering = xmeans(
data = data,
)
x_means_pyclustering.process()
実行結果はエラー！ではなく、ちゃんと実行されてしまいました。実行した結果は下図です。
実行結果は4クラスになっており、意図した結果とは若干異なっております。しかし、一応分類されてしまいました。逆行列の問題はどこに行ってしまったのでしょうか...?
##pyclusteringは実行できた理由
答えは、クラスタ内で等分散を仮定しているからです。
具体的に式を比較してみます。
まず、あるクラスタに含まれる標本を $x_i (i=1, \cdots, N_k)$とします。ただし、$x_i$はベクトルとし、$N_k$ は $k$ 番目のクラスタに含まれるサンプルの個数です。このとき、クラスター数を自動決定するk-meansアルゴリズムの拡張についてで用いている共分散行列は、$R$を標本の次元として
V = \frac{1}{N_k-R}\sum_{i=1}^{N_k}(\boldsymbol{x}_i-\bar{\boldsymbol{x}})(\boldsymbol{x}_i-\bar{\boldsymbol{x}})^\top \\
となります。ただし、
\bar{\boldsymbol{x}} = \frac{1}{N_k}\sum_{i=1}^{N_k}\boldsymbol{x}_i
です。この場合は、ベクトル内の複数ある要素のうち値の変化がない要素が $x_i$ の中にあると $V$がランク落ちしてしまい、逆行列が計算できなくなってしまいます。一方で、pyclusteringのxmeansでは、共分散を
\hat{\sigma}^2 = \frac{1}{N_k}\sum_{i=1}^{N_k}(\boldsymbol{x}_i-\bar{\boldsymbol{x}})^\top(\boldsymbol{x}_i-\bar{\boldsymbol{x}}) \\
V = \hat{\sigma}^2 \boldsymbol{I}_R
としています。ただし、$\boldsymbol{I}_R$は$R$次元の単位行列です。このときは $\hat{\sigma}^2 = 0$でない限りは逆行列が計算できます。これは推測ですが、原著論文においては後者の等分散の仮定をおいていることから、pyclusteringは原著論文にしたがったのではないかと思われます。
※ちなみに、$\hat{\sigma}^2 = 0$ となるのは、クラスタ内の全ての点が同一点のときです。pyclusteringのxmeansはスクリプトをみたところ、クラスタ内に同一点が存在した場合は分割をやめるようになっていました。
#まとめ
今回記事に載せたスクリプトには以下の特徴がありました。
- @deaikeiさん作成のスクリプトは、等分散の仮定をおかずBICを計算するため、ある次元について一定の値が並ぶデータについては計算できない。データに求める仮定は緩い。
- pyclusteringのxmeansは等分散の仮定をおいてBICを計算するため、ある次元について一定の値が並ぶデータについても計算できるということが分かりました。データに求める仮定が厳しい。
今後xmeansを用いてクラスタリングをする場合は、データの特徴を掴んで適切なスクリプトを選択しようと思います！
#参考文献
- Qiitaの記事
- 論文
- pyclustering公式ドキュメント