Stochastic Block Model を Edward で実装する
前回の記事で Edward を使ってみたらすごく良かったので、もう一回遊んでみる。今回はグラフクラスタリングによく使われる Stochastic Block Model (SBM) を実装する。
前回の記事はこれ。
ちなみにプルリク送ったらマージされたので、コードはリポジトリの edward/examples/stochastic_block_model.py にある。
Stochastic Block Model
Stochastic Block Model (SBM) はグラフの確率的生成モデルの一つ。グラフの確率的生成モデルと言うと、有名どころでは Erdos-Renyi model とか Barabasi-Albert model とかあるけど、そういうやつ。
SBM がどういうモデルなのか、すごく簡単に説明する。ものすごくざっくり言うと、SBM は以下の2つの仮定を置いたモデルだといえる。
- すべてのノードは、K 個あるクラスタのうち(必ず)どれか一つに属している。
- あるノード i とノード j の間にエッジが存在する確率は、i と j がどのクラスタに属しているか のみ によって定まる。
めちゃくちゃシンプルなモデルなんだけど、これが結構注目されていて、いろいろ理論的な解析がされていたり、グラフクラスタリングに応用されていたり、拡張モデルがいろいろ提案されていたりする。
例を使って説明してみよう。例えば、グラフ上に3つのクラスタがあり、それぞれのクラスタ間には以下の確率でエッジが張られるということがわかっているとしよう。
$$ \left( \begin{array}{ccc} 0.9 & 0.2 & 0.2 \\ 0.1 & 0.3 & 0.7 \\ 0.1 & 0.8 & 0.1 \end{array} \right) $$
これは例えば、(2,3)要素を見ると、2番目のクラスタに属すノードから3番目のクラスタに属すノードへエッジが張られる確率は 0.7 であるということを表している。上にも書いたように、SBMではクラスタ2に属すどんなノードからも、クラスタ3に属すノードへエッジが張られる確率は等しく 0.7 であるという仮定を置く。 ここが SBM で最も重要。グラフ上のノードはいろいろな性質(例えばTwitterユーザなら性別、年齢、趣味、居住地などなど)を持つわけだけど、それらはすべて無視して、そのノードがどのクラスタに属すかのみによって張られるエッジが決まる。ここが SBM がシンプルたる所以。
さて、実際の生成過程を見ていこう。ノード数は N 、クラスタ数は K とする。
まず、1 から N の i について、クラスタ割りあて を以下のように生成する。
$$ z_{i} \sim Categorical(\gamma) $$
ここで、 は K 次元のベクトルで、各クラスタにどれくらいの割合でノードが所属するかを表すパラメータ。Categorical 分布のパラメータなので、当然 となっている必要がある。 は 1 から K の整数で、例えば ならノード i がクラスタ2に属すことを表す。
つぎに、1 から N の i と j のすべての組み合わせについて、隣接行列 の要素を以下のように生成する。
$$ X_{ij} \sim Bernoulli(\Pi_{z_{i}z_{j}}) $$
は上に書いたような K x K 行列で、各クラスタ間にどれくらいの確率でエッジが存在するかを表すパラメータ。Bernoulli 分布のパラメータなので、 の要素はすべて0以上1以下である必要がある。
生成過程は以上。とてもシンプル。2つのパラメータ と を決めてやればなんらかの隣接行列 が確率的に生成される。また、それと同時に各ノードがどのクラスタに属すかを表す も生成される。 つまり、ノードのクラスタリングが得られる。
ちなみに、簡単のため上では と の生成過程しか書かなかったけど、多くの場合、 と に次のように事前分布を入れる。
$$ \gamma \sim Dirichlet(\alpha) $$
1 から K のすべての k, l の組み合わせについて、
$$ \Pi_{kl} \sim Beta(a, b) $$
ここで、 は Dirichlet 分布のパラメータで、 と は Beta 分布のパラメータ。
モデルが定義できたらパラメータを推定したくなるのが世の常なわけで、それを Edward つかってやっちゃいましょうってのが今回の記事の趣旨。普通は VB やら MCMC やらを頑張って導出して頑張って実装するんだけど、Edward を使うとその辺をすっ飛ばして楽に実装できる。ちなみに、MCMCの導出は毎度おなじみ『関係データ学習』に詳しく載ってるので参考にしてもらえるといいとおもう。
- 作者: 石黒勝彦,林浩平
- 出版社/メーカー: 講談社
- 発売日: 2016/12/07
- メディア: 単行本(ソフトカバー)
- この商品を含むブログ (1件) を見る
Edward で実装
さて、いよいよ Edward で実装する。Edward そのものの使い方については本家チュートリアルをみてもらうなり、前回の記事をみてもらうなりしてほしい。
まず必要なものをインポートする。適当に乱数シードも決めておく。
import edward as ed import numpy as np import tensorflow as tf from sklearn.metrics.cluster import adjusted_rand_score from edward.models import Bernoulli, Multinomial, Beta, Dirichlet, PointMass ed.set_seed(42)
そんで次にデータを用意する。edward/examples/data/
以下に必要なデータが置いてあるので、それを読んでいる。この example コードではかの有名な Zachary’s karate club のデータを使う。
def build_dataset(label_filepath, graph_filepath): Z = np.loadtxt(label_filepath, dtype=np.int) N = Z.shape[0] X = np.zeros((N, N)) for line in open(graph_filepath, 'r'): src, dst = map(int, line.strip().split(' ')) X[src, dst] = 1 return X, Z # DATA label_filepath = 'data/karate_labels.txt' graph_filepath = 'data/karate_edgelist.txt' X_data, Z_true = build_dataset(label_filepath, graph_filepath) N = X_data.shape[0] # number of vertices K = 2 # number of clusters
モデルを定義する。上に書いた生成過程をそのままコードにしたような感じで、非常に書きやすい。 と にも事前分布を入れている。tf.matmul
は行列積で、tf.transpose
は転置をしている。
# MODEL gamma = Dirichlet(concentration=tf.ones([K])) Pi = Beta(concentration0=tf.ones([K, K]), concentration1=tf.ones([K, K])) Z = Multinomial(total_count=1., probs=gamma, sample_shape=N) X = Bernoulli(probs=tf.matmul(Z, tf.matmul(Pi, tf.transpose(Z))))
モデルが定義できたらいよいよパラメータを推定する。今回は簡単に MAP を EM algorithm で推定する。PointMass
はそのパラメータを点推定するよということを宣言していると思っていい。
# INFERENCE (EM algorithm) qgamma = PointMass(params=tf.nn.softmax(tf.Variable(tf.random_normal([K])))) qPi = PointMass(params=tf.nn.sigmoid(tf.Variable(tf.random_normal([K, K])))) qZ = PointMass(params=tf.nn.softmax(tf.Variable(tf.random_normal([N, K])))) inference = ed.MAP({gamma: qgamma, Pi: qPi, Z: qZ}, data={X: X_data}) n_iter = 100 inference.initialize(n_iter=n_iter) tf.global_variables_initializer().run() for _ in range(inference.n_iter): info_dict = inference.update() inference.print_progress(info_dict) inference.finalize()
パラメータが推定できたので、最後にモデルを評価する。評価指標には Adjusted Rand Index (ARI) を使う。
# CRITICISM Z_pred = qZ.mean().eval().argmax(axis=1) print("Result (label filp can happen):") print("Predicted") print(Z_pred) print("True") print(Z_true) print("Adjusted Rand Index =", adjusted_rand_score(Z_pred, Z_true))
結果は以下のようになる。クラスタリング結果と正解が完全に一致していれば ARI は 1 になるんだけど、一つ間違っていたので ARI = 0.88 になった。まぁそれでも SBM はこのデータをかなり正しくクラスタリングできたと言えると思う。
Result (label filp can happen): Predicted [0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 0 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1] True [0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1] Adjusted Rand Index = 0.882257541356
まとめ
やっぱ Edward いいわ。