The 19th Annual Conference of the Japanese Society for Artificial Intelligence, 2005

1F3-01

局所尤度推定に基くノイジーデータからの大規模分散クラスタリング
Large-scale Distributed Clustering from Noisy Data using Localized Likelihood Estimation
佐久間

淳∗1

小林 重信∗1Jun Sakuma∗1

Shigenobu Kobayashi

東京工業大学 大学院 総合理工学研究科

Tokyo Institute of Technology, Interidisciplinary Graduate School of Science and Engineering

This paper presents a new clustering algorithm that can be easily extended for distributed computation. Our proposal algorithm, Distributed Clustering based on Mixture Modeling (DCMM), assumes that the noise is distributed
following a uniform distribution and the each cluster is estimated as a component of a mixture distribution. Basedon this assumption, the estimation is conducted by each component in a serial manner, that is, the computation canbe conducted on multiple CPUs separately. The eﬀectiveness of the proposal is veriﬁed by the numerical examples
using artiﬁcially generated data sets.

1.

はじめに

に基づき発生したと考える. ここで, k は混合数 (クラスタ数), αj は αj > 0, j αj = 1 なる混合係数である. X のオンラインで入手可能なデータの増加に伴い, データが潜在的に有する知識を獲得するためのツールへの要求が高まっている. クラスタリングは教示情報を用いずに, 類似度の高い複数のデータを一つのクラスタに割り当てる方法であり, データ
セットの大域的な構造を俯瞰することで, データセットが非明
示的に有する知識の獲得を容易にする. 大規模情報源からのス
ケーラブルなクラスタリングアルゴリズムは活発に研究されいているが [1], たとえデータ数 N について O(N ) の計算量が実現されているとしても N が極めて多い場合にはアルゴリズムの並列分散化が必須である.
また実問題で利用されるデータは通常どのクラスタにも属さないであろうノイジーデータを含むが、このようなデータセッ
トを扱うためには, 確率分布に基づくクラスタリングアルゴリ
ズムが有効であると考えられる. 確率分布に基づくアルゴリズ
ムは理論的背景が明確であるが, 必ずしも現実のデータに則してアルゴリズムが設計されているとは言いがたく, 大規模デー
タにおいて時間計算量および空間計算量が問題となる. たとえばクラスタ数の決定には正則化項の導入によるモデル選択が用いられるが, クラスタ数が数百や数千に達するような状況ではこれは計算量の面からも精度の面からも適切ではない.
本稿では EM 法を利用した混合分布推定に基いて, あるクラスタに属さないデータは一様分布に従うとの仮定の下で 1 つのクラスタを推定するアルゴリズムを提案する. このアルゴリズ
ムは, (1) 並列実行が可能なため分散化が容易である, (2) デー
タに多量のノイズが混入している場合でもその影響を受けにくい, (3) 終了条件が必ずしも明確ではないが, クラスタ数を前もって与える必要がない, 等の利点を持つ. 提案手法の有効性は人工的に生成したデータセットを用いた計算機実験によって検証される.

2.

p(x; θ) に対する尤度 (= Ni=1 p(xi ; θ)) を最大にするようなパラメータを θM L とする (最尤推定). このとき, 確率モデルに基づくクラスタリングでは, データ x について最大の確率密度
pj (x; θM L) を与える j をデータ x が属するクラスタと決定する. 混合分布の最尤推定には一般に EM 法 [3] を用いる.
[EM]
1. 初期パラメータ θ(0) を設定, t ← 0.
2. Q(θ|θ (t) ) = E{log p(X, Z|θ)|X, θ (t) } を計算 (E-step).
3. θ(t+1) = arg maxθ Q(θ|θ (t) ) とする (M-step).
4. 収束条件を満たせば終了. そうでなければ t ← t + 1 として 2 へ.
上記の EM 法は対数尤度の局所収束性を保証する.

3.1

混合分布推定に基づく分散クラスタリング

Hard assignment: データ xi が j 番目の分布から発生した確率を p(zij |xi) とする. この条件付確率について, 常にp(zij |xi) ∈ {0, 1}, j p(zij |xi ) = 1 が満たされる. これはあるデータ xi が分布 pj から発生する確率がつねに {0, 1} であると考えることに相当し, k-mean 法などと同等の制約である.
本論文ではこの性質を Hard assignment と呼ぶ.

クラスタリングの対象となるデータ集合を X = {x1 , ..., xN },N をデータ数とする. モデルパラメータ θ で特徴づけられた確率分布を p(x; θ) とする. このときデータ集合 X が, 混合分布

この仮定に基づき 1 式の混合分布推定に基づく分散クラス
タリング (Distributed Clustering based on Mixture Model,
DCMM) のアルゴリズムの概要を Algorithm 1 に示す.
DCMM は主に”LLEM”, ”混合比推定”および”クラスタの検定の 3 つのパートで構成されており, これらの詳細は次章で説明する. DCMM の特徴を以下に説明する.

kαj pj (x; θj )分散クラスタリング法の提案

EM 法では E-step において Q(θ|θ (t) ) の計算に隠れ変数に関する条件付確率 p(zij |xi ; θj ) が必要であるが, その計算に iと j に関する和計算が必要であるため, 少なくとも O(kN ) の計算量を要する. 多くのクラスタがオーバーラップするようにデータが分布する場合, これらの計算は避けがたいが, 各クラスタがある程度独立して分布している場合, p(zij |xi ; θj ) はある j について 1 に近い値をとり, その他については 0 に近い値をとることから, N 個の全データに対する計算は実際には不要であることが多い. そこで本研究では確率モデルにおいて以下の性質を仮定したうえで条件付尤度の計算を分散化する.

混合分布と EM 法

p(x|θ) =

3.

(1)

j=1

連絡先: 佐久間 淳，東京工業大学 大学院 総合理工学研究科，神奈川県横浜市緑区長津田町 4259 G5-801，045-924-5544，
045-924-5442，jun@fe.dis.titech.ac.jp1

The 19th Annual Conference of the Japanese Society for Artificial Intelligence, 2005
Algorithm 1 DCM M

Algorithm 2 LLEM

1: X, α0 を入力
2: j←0
3: while 収束条件が満たされない do
4:αj ← α05:while αj , θj が収束しない do
6:
αj を固定し, pj (x; θj ) における θj を推定 (LLEM)
7:
θj を固定し, αj を推定 (混合比推定)
8:end while9:
pj (x; θj ) をクラスタとして受理するかどうかを検定
10:
if 検定が真 then11:
Cj ← X において密度 pj (x; θj ) が高い αj N 個のデータ
12:
X ← X \ Cj
13:j ← j+1
14:
end if
15: end while
16: Cj , αj , θj (j = 0, ...) を出力

1: X を入力
2: t → 0(0)
3: θj および P̃ (Z|Φ)(0) を初期化
(t)4: while θj
5:
6:

i

上 記 の E-step，M-step の 両 方 に お い て ，具 体 的 な(t+1)の求め方は示されていない．以下P̃ (Z|Φ)(t+1) および θjに，両 step におけるその実現方法を考察する．
E-step では K の最大化は困難であるが，増大化であれば簡単な局所探索等のヒューリスティックによって以下のように実現可能である．
(t)

1. 固定された θj において，下式を満たすような部分集合
(t+1)
Dj
を生成

(t)

=

(t)

=

1

if

P̃ (zij |Φ)(t+1)

=

0

else

xi ∈ D j

(t)

3.3

混合比推定

上述した LLEM は局所尤度に関する局所収束性を保証する
が, αj が固定値として構成されている. この値を推定するために, あるモデル pj に着目したときに, pj 以外のモデルについては PU ( ; u) と同一視して，(4)xi ∈DjΣj Qj (θj |θj ).

(t+1)

P̃ (zij |Φ)(t+1)

また M-step における K(P̃ (Z|Φ)(t) |θj ) の増大化法は，確率モデルおよび事前分布の具体的な形状に依存しているが, デー(t)タセット Dj についての最尤推定で最大化することができる.
上記, E-step, M-step に基づいて, LLEM 法は混合モデルのなかの一つのモデルに関する局所尤度について，局所収束 (局所最大化) が保証される．証明は [5] を参照されたい．

pj (xj , zij ; θj ) (3)log pj (xi , zij ; θj )jlog Pj (X, Z; θj )
(t)
x∈Dj

2. P̃ (Z|Φ)(t+1) を以下のように設定:

j

=

(t)

log Pj (X, Z; θj ) ≥
(t+1)
x∈Dj

この計算に i と j の和計算が必要であるため, EM アルゴリ
ズムは少なくとも O(kN ) の計算量を要することに注意されたい. ここで隠れ変数の条件付確率について, 常に p(zij |xi ) ={0, 1}, j p(zij |xi ) = 1 が満たされる (Hard assignment 仮定) ならば, Q 関数はQ(θ|θ (t) )(t)) ≥ K(P̃ (Z|Φ)(t+1) , θj ) となるよのように, 各分布毎に和の形式で書くことができる. ただし，
I(i) はデータ xi に対して zij = 1 なるインデックス j を返す
関数，Dj = {xi |I(i) = j} である．これを条件付対数局所尤
度 (Conditional Log Localized Likelihood) と呼ぶ．
条 件付 対 数局 所 尤度 が コ ンポ ー ネン ト 毎に 個 別に 計 算可能である性質を用いて，混合分布推定のための局所尤度
EM(Localized Liklihood EM, LLEM) を構築する [5]．分布は p0 , ..., pj−1 の順番に推定されるものとし, コンポーネント
pj を推定する LLEM のアルゴリズムを説明する. ここでは分
布 pj に hard assignment されるデータ数を αj N = cj 個と固定されるものとする．また，P̃ (Z|Φ) を hard assignment 仮定を満たす関数とする．
K(P̃ (Z|Φ), θj ) = Z P̃ (Z|Φ) log pj (X; |θj ) とすると，コンポーネント pj (x; |θj ) における LLEM のアルゴリズムは
Algorithm 2 のように示される．

Z
(t)

(t+1)

K(P̃ (Z|Φ)(t+1) , θj

7:
t→t+1
8: end while(t)(t)
9: θj および K(P̃ (Z|Φ)(t) , θj ) = 1 なるデータを出力

E-step で出現する Q 関数は混合分布の場合以下のように書き下せる．
Q(θ|θ (t) ) =p(Z|X; θ(t) ) log p(X, Z; θ)(2)p(zij |xi; θj ) log

(t)

(t+1)
を決定 (M-step)
うに，θjLocalized Liklelihood EM=(t)

K(P̃ (Z|Φ)(t+1) , θj ) > K(P̃ (Z|Φ)(t) , θj ) となるようにP̃ (Z|Φ)(t+1) を決定 (E-step)

推定の分散化: Hard Assinment 仮定によって Q(θ|θ (t) ) の計算を各分布毎に局所化し, j αj pj (x; θj ) の推定を全分布の同時推定ではなく, 各分布毎に局所的に推定する (step. 5-8).
あるクラスタの推定には必ずしも全データセットは不要であり, 興味のある領域に存在するデータセットについてアルゴリ
ズムを実行すれば十分である. データセットを適宜分割することで, O(N k) の N を減少させることができる. また推定の分散化によって, O(N k) の k について, 分布毎に異なる CPU を割り当てることで O(N ) に減少させることができる.
クラスタ数の事前の決定が不要: DCMM ではクラスタ数のかわりに混合係数 α0 を初期値として与える. 混合係数はそのクラスタが保持するデータ数の全データ数に対する割合である. 混合係数を固定した上でモデルで pj (x; θj ) を推定し (step6), 推定パラメータ θj の下で αj を推定する (step 7). Step 6
および 7 の詳細については 3.2 および 3.3 をそれぞれ参照されたい.
ノイズ対応のための密度分布に関する検定: EM 法は局所最適性しか保証しないため, Step 6 および Step 7 によって求められたクラスタが望ましい形状をしているとは限らない. そこで Step 9 ではクラスタの形状について, そのクラスタがノイズかどうか, および, 密度の形状が単峰であるかどうかをコル
モゴロフ検定および二項検定により検定する.
本章では任意の分布形の混合分布 j αj pj に基づく分散ク
ラスタリング法を説明し, 次章では分布 pj が正規分布を取る場
合の具体的なアルゴリズムについて説明する. 以下, Step 6,7を中心に分散クラスタリング法のアルゴリズムの理論的バックグラウンドとその手順を説明する. 検定部分 (Step 9) の詳細については [4] を参照されたい.

3.2

が収束 dop(x; θ) = αj pj (x; θj ) + (1 − αj )pU (x; u),

(5)

2

(6)

The 19th Annual Conference of the Japanese Society for Artificial Intelligence, 2005
なる確率モデルを考えることにする。ここで，pU はデータの領域全体にわたる一様分布，u は一様分布におけるモデルパ
ラメータである. pj についての LLEM が実行されていると考えれば, θj はすでに推定済みである. u を固定した場合, 尤度
は極値を持つならば, σ に関する一次元最適化によって望ましい σ を得ることができる. たとえば u を正規分布の 2σ 線∗1 を，クラスタの境界と考え，その外は一様分布であるとみなせば発散せず安定した推定が可能でる. 逆に，u を決めたうえでσ を探索し，なお σ が発散するならば，それは正規分布をなさない，と結論することができる. 適切な σ が決定されれば
pj (x; θj ) > u となるようなデータを収集し, これに基づいて
αj を決定することができる.

4.

図 1: クラスタ数 k = 32, 混入ノイズ数 r = 800, データ次
元数 d = 2 において生成したデータセットに DCMM を適用した際の推定の様子
50

30
20

x2

10

-20
-30

前章では任意の分布形 pj における DCMM の構成法を説明した. 本節では最も基本的な多次元正規混合分布におけるDCMM の具体的なアルゴリズムについて説明する.
多次元正規分布はモデルパラメータ平均値ベクトル µ および分散共分散行列 Σ によって特徴つけられる. d 次元正規分布は下式のように記述される.

1(2π)d/2 |Σ|

exp
1/2

-40
-50
-50

-20

-10

0

10

20

30

40

50

1: X, α0 , s を入力
2: µ, Σ を初期化
3: t ← 0, s ← α0 N
4: while 局所尤度が未収束 do
5:xi ∈ X を P (x; µ, Σ) における密度に基づきソートし，上位 s

個を選択 (Lists )Lists から s 個の要素を密度に基づきルーレット選択し，これを Lists とする
7:
L ← xi ∈Lists P (xi ; µ, Σ)(局所尤度)
8:
if L > Lbest then9:
Lists から µ(t) , Σ(t) を最尤推定
10:µ ← µ(t) , Σ ← Σ(t)11:
L ← Lbest
12:
else
13:
Step 6 へ戻る
14:
end if
15:t←t+1
16: end while
17: µ, Σ, Lists を出力
6:

k

(7)

j=1

となる. ただし j αj + αU = 1 である. 本章ではデータ
X が与えられたときに αj , µj , Σj を推定するアルゴリズムを考える.
混合正規分布における LLEM: LLEM では，混合比を固定し局所尤度を極大化するアルゴリズムである. 入力の α0は対象とするモデルに対する混合比の初期値である. s は
s > αj N なるパラメータであり, E-step における Dj の探索において, hard assignment するデータの候補数である. 実験
的には s = 4α0 N 程度で十分である.

評価基準として以下に定義する accuracy と error を用いる．
データにはその発生元となった分布があらわすラベルがついている (アルゴリズムにはこの情報は不可視である). クラスタリ
ングの結果, そのクラスタにおいて最も多くのデータが持つラベルを majority と呼ぶ.

混合正規分布における混合比推定: クラスタの平均ベクト
ル µ および固有ベクトル e1 , ..., ed を固定し, 固有値 a1 , ..., akを変化させ混合比 αj を推定する．固有値を LLEM における推定値の rmin (< 1) 倍から rmax (> 1) について局所尤度 L を
1 次元探索し, 最大となる r∗ を探索する．

• True Positive(TP) : あるクラスタに属するデータの中の, majority のデータ数• False Positive(FP) : あるクラスタに属するデータの中の, 非majority のデータ数• True Negative(TN) : あるラベルを有するデータの中の, 非majority のデータ数

上記のアルゴリズムを組み合わせることで, 正規混合分布における DCMM が構成される. クラスタ数 k = 16, 混入ノ
イズ数 r = 800, データ次元数 d = 2 において生成したデータセットに DCMM を適用した際の推定の様子を図 1 に示す.

5.

-30

Algorithm 3 混合正規分布における LLEM

1
− (x − µ)T Σ−1 (x − µ)2

αj P (x; µj , Σj ) + αU PU (x; u)

-40

x1

固有ベクトル e1 , ..., ed について行列 E = { 1 , ..., k } とし,Σ の固有値 a1 , ..., ad について対角行列 A = diag(a1 , ..., ad )とすると, Σ = E T AE となる.
前章同様, 一様分布を PU = (x; u) とすると, 複数の正規分
布と一様分布の混合分布はPi(x; θ) =0
-10

混合正規分布における DCMM

Pi(x; θ) =’data_d2_k32_r800.dat’’contour_d2_k32_r800.dat’

40

このとき以下に定義する accuracy および error を評価基準とする.
P• accuracy(acc) = T PT+T
: あるラベルが，単一のクラスタにN

実験属している率P• error(err) = T PF+F
: あるクラスタについて，majority ではP

7 式においてクラスタ数 k, ノイズ混入数 r およびデータ次
元数 d を変化させて生成したデータについて実験を行う．1 ク
ラスタあたりのデータ数は [100, 200] における一様乱数である.

ないデータの含有率

accuracy は高いほど精度の高いクラスタを生成していることを表し, error は低いほど誤ったデータをクラスタのメンバとしてはいないことをあらわす. これに加えて, 発見したクラ∗1 正規分布から生成されたデータの 95%が 2σ 線内に存在する

3

The 19th Annual Conference of the Japanese Society for Artificial Intelligence, 2005
1.1
r=200
1

r=3200
30

0.4

0.8

err

0.3

r=32000.2
r=3200

0.7

0.1

0.6

0

0.5

-0.1

0

5

10

15

20

25

30

35

25
20

10

0
0

5

10

15

20

25

30

0

35

0.4

10

0.3

0.7

# of found cluster

1
0.9
err

12

0.8

0.2
0.1

0.6

0

0.5

-0.1
10 15 20 25# of dimension

30

35

10

15

20

25

30

35

10 15 20 25# of dimensiona

30

35# of cluster

0.5

5

5

# of cluster

1.10

r=200,400,800

155r=200 r=400 r=800# of cluster

accuracy# of found cluster

r=800

0.9
accuracy

35

0.5

r=400

8
6
4
2
0

0

5

10 15 20 25# of dimension

30

35

0

5

図 2: 上段は左から，クラスタ数 vs accuracy, error, 見つかったクラス数. 下段は左から，次元数 vs accuracy, error,
スタ数. データはすべて 10 試行平均値．

6.

Algorithm 4 M ixtureRatioEstimation
1: X, Lists , , i (i = 1, ..., d), rmin , rmax , rstep を入力
2: r0 ← rmin , t ← 0
3: while rt < rmax do
4:
Σt ← E T (rt I)E5:xi ∈ X を Pj ( ; , Σ) における密度に基づきソートし,

発見クラ

おわりに

本稿ではノイズモデルとして一様分布を想定した確率モデ
ルにおいて, 混合分布に基づくクラスタリングアルゴリズム
DCMM を提案した. DCMM はノイズに対するロバスト性が強く, また次元数に対する性能のスケーラビリティにも優れていることを実験により明らかにした. 今後はオンメモリでの処理が困難なギガデータに対してグリッド計算機を利用したク
ラスタリングを試みる予定である. また XML などの半構造化
データや文書データ, 木構造で記述されたルール, 関係グラフ,などのさまざまな非構造化データにおけるクラスタリングに拡張する予定である.

Pj ( ; , Σ) > d なる xi を Listt に入れる6:
Lt ← xi ∈Listt P ( i ; , Σ)(局所尤度)
7:
rt+1 ← rt + rstep , t ← t + 1
8: end while
9: t∗ ←
Lt を最大にする t
10: r ∗ ←Lt を最大にする rt∗
T ∗11: List∗t と Σ (= E (r I)E) を出力

参考文献

スタ数を評価基準に用いる. 提案手法はクラスタ数を前もって与えないため, 発見するクラスタ数は試行によって異なる. これは真のモデルのクラスタ数に近いほど望ましい.
ロバスト性: 推定のロバスト性を見るために，クラスタ数 k
およびノイズ混入数 r を変化させて実験を行う．この実験では k = {4, 8, 16, 32}, r = {0, 200, 800, 3200} としてデータ次
元数を d = 2 とした．評価基準として accuracy, error および発見クラスタ数を用いる．図 2(上) は異なるノイズ混入数について, 真のクラスタ数に対する accuracy, error, 発見クラスタ
数をプロットしたものである. ノイズは一様にデータに混入するため, r が大きいほど accuracy は悪化する. ただし提案手法
では確率モデルが正確に推定されているため, error および発見クラスタ数はノイズの影響を受けていないことがわかる. またデータの真のクラスタ数が増えても, これらの評価基準においては性能が悪化しないことも確認できる.
スケーラビリティ: アルゴリズムのスケーラビリティを見るために，ノイズ混入数 r = 800, クラスタ数 k = 8 として, 次
元数 d を d = {2, 4, 8, 16, 32} と変化させて実験を行った (図2(下))．評価基準は前節と同様である．次元数についても, 精
度・発見クラスタ数ともに正確な推定を行っており, 高次元データにおいても性能の悪化がないことを確認できた.

[1] Domingos, P. and Hulten, G.: A General Method for Scaling Up Machne Learning Algorithms and its Application to
Clustering, Proc. of the 18th Int’l Conf. on Machine Learning, pp. 106-113 (2001)[2] Kollios, G., Gunopulos, D, Koudas, N., Berchtold, S.: Efﬁcient Biased Sampling for Approximate Clustering andOutlier Detection in Large Datasets, IEEE Transactions on
Knowledge and Data Engineering, Vol. 15, No. 5, (2003).
[3] Depmster, A. P., Laird, N.M. and Rubin, D. B. (1977),
Maximum likelihood from incomplete data via the EM algorithm, J. Roy. Statist. Soc. B, 39:1-38 (1977).
[4] 佐久間, 小林, 局所尤度推定に基づくノイジーデータからの分散
クラスタリング, 第 35 回システム工学部会研究会, (2005).
[5] 佐久間, 小林 , 依存性を有する事前分布における Localized MAPEM, 情報論的学習理論ワークショップ, pp. 289-294 (2003).

4