Hartigan-Wong のアルゴリズムを確認する

R の stats パッケージで提供されている kmeans 関数は、既定では Hartigan-Wong のアルゴリズムを利用します。通常の k-means (Lloyd のアルゴリズム) では、各データ点を最も近いクラスタに割り当てる操作を繰り返しますが、Hartigan-Wong の方法はより直接的に、量子化誤差の増分を最小化するクラスタにデータ点を割り当てる方法になっています。

Hartigan-Wong の論文*1は下記のウェブサイトにあります*2。
http://www.jstor.org/stable/2346830

また、下記の 2 本の論文*3 *4などで、このアルゴリズムについて論じられています。たとえば前者によると、Hartigan-Wong の方法で見つかる局所解は Lloyd の方法で見つかる局所解の真部分集合になっているとのことです (Theorem 2.2)*5。
http://www.jmlr.org/proceedings/papers/v9/telgarsky10a/telgarsky10a.pdf
http://ijcai.org/papers13/Papers/IJCAI13-249.pdf

後者の論文の内容が私には分かりやすく思えたので、こちらで Hartigan-Wong のアルゴリズムを紹介します。論文の式 (2) が Hartigan-Wong のアルゴリズムの本質です。この式は、クラスタにデータ点を追加することによって生じる量子化誤差の増加量を意味します。距離関数 $d(\cdot, \cdot)$ が二乗ユークリッド距離の場合には、式 (7) のように変形できます。また、同論文の Figure 1 は、アルゴリズムを擬似コードで示したものです。Figure 1 の ${\Delta}D(x, c)$ を式 (7) で具体的に計算すれば、Hartigan-Wong のアルゴリズムを実装できます。

二乗ユークリッド距離の場合に式 (2) が式 (7) のように変形されるということが、それぞれの式を見比べただけでは私には理解できませんでした。そこで、実際に式変形を行って、このことを確認してみます。

点 $x \notin C$ をクラスタ $C$ に追加したときの量子化誤差の増分は、論文の式 (2) のとおり、次式で表せます。ただし、小文字の $c$ は $x$ を追加する前のセントロイド、 $c^+$ は $x$ を追加した後のセントロイドとします。括弧内の二項が追加後の量子化誤差、最後の項が追加前の量子化誤差に相当します。

${\displaystyle {\Delta}D(x, C) = (\|x - c^+\|^2 + \sum_{x' \in C}\|x' - c^+\|^2) - \sum_{x' \in C}\|x' - c\|^2 }$

この式から $x'$ と $c^+$ を消去して $x$ と $c$ のみで表し、それが論文の (7) 式になっていることを示します。

まず、 $c^+$ は $x$ と $c$ を用いて以下のように表せます。

${\displaystyle \begin{eqnarray} c^+ &=& \frac{1}{|C^+|} (\sum_{x' \in C} x' + x) \\ &=& \frac{1}{|C|+1} \sum_{x' \in C} x' + \frac{1}{|C|+1} x \\ &=& \frac{|C|}{|C|+1} c + \frac{1}{|C|+1} x \end{eqnarray} }$

これを用いると、 ${\Delta}D(x, C)$ の第一項は以下のように変形できます。

${\displaystyle \begin{eqnarray} \|x - c^+\|^2 &=& \|x - (\frac{|C|}{|C|+1} c + \frac{1}{|C|+1} x)\|^2 \\ &=& \|\frac{|C|}{|C|+1}(x - c)\|^2 \\ &=& (\frac{|C|}{|C|+1})^2 \|x - c\|^2 \end{eqnarray} }$

残りの二項は、次の形にまとめて考えます。

${\displaystyle \sum_{x' \in C}(\|x' - c^+\|^2 - \|x' - c\|^2) }$

総和記号の内側は次のように計算できます。

${\displaystyle \begin{eqnarray} \|x' - c^+\|^2 - \|x' - c\|^2 &=& (\|x'\|^2 - 2 x' \cdot c^+ + \|c^+\|^2) - (\|x'\|^2 - 2 x' \cdot c + \|c\|^2) \\ &=& \|c^+\|^2 - \|c\|^2 - 2 x' \cdot (c^+ - c) \\ &=& (c^+ + c) \cdot (c^+ - c) - 2 x' \cdot (c^+ - c) \\ &=& (c^+ + c - 2 x') \cdot (c^+ - c) \end{eqnarray} }$

したがって、

${\displaystyle \begin{eqnarray} \sum_{x' \in C}(\|x' - c^+\|^2 - \|x' - c\|^2) &=& \sum_{x' \in C}\{(c^+ + c - 2 x') \cdot (c^+ - c)\} \\ &=& (c^+ - c) \cdot \sum_{x' \in C}(c^+ + c - 2 x') \\ &=& (c^+ - c) \cdot (|C| \; c^+ + |C| \; c - 2 \sum_{x' \in C} x') \\ &=& |C| \; (c^+ - c) \cdot (c^+ - c) \\ &=& |C| \; \|c^+ - c\|^2 \end{eqnarray} }$

となります。ここで、 $c^+$ が

${\displaystyle c^+ = \frac{|C|}{|C|+1} c + \frac{1}{|C|+1} x }$

だったことから、上式はさらに

${\displaystyle \begin{eqnarray} |C| \; \|c^+ - c\|^2 &=& |C| \; \|(\frac{|C|}{|C|+1} c + \frac{1}{|C|+1} x) - c\|^2 \\ &=& |C| \; \|\frac{1}{|C|+1}(x - c)\|^2 \\ &=& \frac{|C|}{(|C|+1)^2} \|x - c\|^2 \end{eqnarray} }$

となります。これで ${\Delta}D(x, C)$ の各項を $x$ と $c$ のみで表すことができました。

第一項と合わせて全体をまとめると、最終的に次のようになります。

${\displaystyle \begin{eqnarray} {\Delta}D(x, C) &=& (\frac{|C|}{|C|+1})^2 \|x - c\|^2 + \frac{|C|}{(|C|+1)^2} \|x - c\|^2 \\ &=& \frac{|C|}{|C|+1} \|x - c\|^2 \end{eqnarray} }$

論文の式 (7) と比較してみると、 $(0.5 / n)$ の有無が異なりますが*6、これは k-means の処理では定数項になるため無視できます。以上の計算により、式 (7) を最小化する割り当てが、量子化誤差の増分を最小にする割り当てを意味していることを示せました。

*1:J. A. Hartigan and M. A. Wong. Algorithm AS 136: A K-Means Clustering Algorithm. Journal of the Royal Statistical Society. Series C (Applied Statistics). Vol. 28, No. 1 (1979), pp. 100-108.

*2:論文のダウンロードは有料 ($29.00) です。アカウントを作成すれば、ウェブブラウザ上では無料で読めます。

*3:Matus Telgarsky and Andrea Vattani. Hartigan's Method: k-means Clustering without Voronoi. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics. 2010.

*4:Noam Slonim, Ehud Aharoni, and Koby Crammer. 2013. Hartigan's K-means versus Lloyd's K-means: is it time for a change?. In Proceedings of the Twenty-Third international joint conference on Artificial Intelligence (IJCAI '13), Francesca Rossi (Ed.). AAAI Press 1677-1684.

*5:ちなみに、この Theorem 2.2 は "a (possibly strict) subset" という表現になっていますが、同論文の introduction では "a strict subset" と書かれていて、ややニュアンスが異なるように感じます。

*6:量子化誤差や二乗ユークリッド距離の定義の微妙な違いに由来するものだと思います。