でかいチーズをベーグルする

でかいチーズはベーグルすべきです。

Stochastic Block Model を Edward で実装する

前回の記事で Edward を使ってみたらすごく良かったので、もう一回遊んでみる。今回はグラフクラスタリングによく使われる Stochastic Block Model (SBM) を実装する。

前回の記事はこれ。

yamaguchiyuto.hatenablog.com

ちなみにプルリク送ったらマージされたので、コードはリポジトリedward/examples/stochastic_block_model.py にある。

github.com

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}} を以下のように生成する。

$$ z_{i} \sim Categorical(\gamma) $$

ここで、{\gamma} は K 次元のベクトルで、各クラスタにどれくらいの割合でノードが所属するかを表すパラメータ。Categorical 分布のパラメータなので、当然 {\sum_{k}^{K} \gamma_{k} = 1} となっている必要がある。{z_{i}} は 1 から K の整数で、例えば {z_{i} = 2} ならノード i がクラスタ2に属すことを表す。

つぎに、1 から N の i と j のすべての組み合わせについて、隣接行列 {X} の要素を以下のように生成する。

$$ X_{ij} \sim Bernoulli(\Pi_{z_{i}z_{j}}) $$

{\Pi} は上に書いたような K x K 行列で、各クラスタ間にどれくらいの確率でエッジが存在するかを表すパラメータ。Bernoulli 分布のパラメータなので、{\Pi} の要素はすべて0以上1以下である必要がある。

生成過程は以上。とてもシンプル。2つのパラメータ {\gamma}{\Pi} を決めてやればなんらかの隣接行列 {X} が確率的に生成される。また、それと同時に各ノードがどのクラスタに属すかを表す {z} も生成される。 つまり、ノードのクラスタリングが得られる。

ちなみに、簡単のため上では {z}{X} の生成過程しか書かなかったけど、多くの場合、{\gamma}{\Pi} に次のように事前分布を入れる。

$$ \gamma \sim Dirichlet(\alpha) $$

1 から K のすべての k, l の組み合わせについて、

$$ \Pi_{kl} \sim Beta(a, b) $$

ここで、{\alpha} は Dirichlet 分布のパラメータで、{a}{b} は Beta 分布のパラメータ。

モデルが定義できたらパラメータを推定したくなるのが世の常なわけで、それを Edward つかってやっちゃいましょうってのが今回の記事の趣旨。普通は VB やら MCMC やらを頑張って導出して頑張って実装するんだけど、Edward を使うとその辺をすっ飛ばして楽に実装できる。ちなみに、MCMCの導出は毎度おなじみ『関係データ学習』に詳しく載ってるので参考にしてもらえるといいとおもう。

関係データ学習 (機械学習プロフェッショナルシリーズ)

関係データ学習 (機械学習プロフェッショナルシリーズ)

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

モデルを定義する。上に書いた生成過程をそのままコードにしたような感じで、非常に書きやすい。{\gamma}{\Pi} にも事前分布を入れている。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 いいわ。

Probabilistic Matrix Factorization を導出して Edward で実装する

Edward っていう確率モデリングのためのライブラリがよさげって話を聞いたので入門してみたら良かったという話。せっかくなので、行列分解を確率モデルとして定義した Probabilistic Matrix Factorization を実装してみた。

Edward – Home

行列分解 (Matrix Factorization)

前にも書いた気がするけど、行列分解ってのは N x M 行列 X を、適当な K に対して N x K 行列 U と M x K 行列 V(の転置)との積に分解する手法のこと。つまり、

 {\mathbf{X} \simeq \mathbf{U}\mathbf{V}^T}

となるような U と V 見つければOK。ここで、{\mathbf{X}}{\mathbf{U}\mathbf{V}^T} が近くなる({\simeq} になる)というのは例のごとく二乗誤差で評価する。つまり、

{
L(\mathbf{U}, \mathbf{V}) = \sum_{i = 1}^{N} \sum_{j = 1}^{M} \left( x_{ij} - u_{i}^T v_{j} \right)^2
}

が最小となるような U と V を求める。{u_{i}} は U の i 番目の(K次元)行ベクトルで、{v_{j}} は V の j 番目の(K次元)行ベクトルを表す。要素ごとに書かずに行列ノルムで書くと、

{
L(\mathbf{U}, \mathbf{V}) = || \mathbf{X} - \mathbf{U}\mathbf{V}^T ||^2_2
}

こうなる。これを最小化すればいいんだけど、これもまた例のごとく正則化したほうがいいので、L2-正則化項を足してやると、

{
L_{\lambda}(\mathbf{U}, \mathbf{V}) = || \mathbf{X} - \mathbf{U}\mathbf{V}^T ||^2_2 + \lambda \left( || \mathbf{U} ||_2^2 + || \mathbf{V} ||_2^2  \right)
}

こうなる。{\lambda}正則化パラメータ。これをどう最小化するかはたぶん Web 上にたくさん説明があると思うのでここでは説明しない。

Probabilistic Matrix Factorization(論文

天下り的に説明した行列分解だけど、実は確率的生成モデルとして解釈することができる。以下のモデルを考えよう。

(1) 1 から N のすべての i について、平均0、分散 {1/\lambda}正規分布から {u_{i}} を生成する。

{u_{i} \sim N(\cdot | \mathbf{0}, 1/\lambda)}

(2) 1 から M のすべての j について、平均0、分散 {1/\lambda}正規分布から {v_{j}} を生成する。

{v_{j} \sim N(\cdot | \mathbf{0}, 1/\lambda)}

(3) 1 から N の i 、1 から M の j のすべての組み合わせについて、平均 {u_{i}^T v_{j}} 、分散 1 の正規分布から {x_{ij}} を生成する。

{x_{ij} \sim N(\cdot | u_{i}^T v_{j}, 1)}

与えられたデータ {\mathbf{X}} がこのモデルから生成されたと考えて、{\mathbf{U}}{\mathbf{V}} をMAP推定してみよう。{\mathbf{X}}{\lambda} が与えられたもとでの {\mathbf{U}}{\mathbf{V}} の事後分布は

{
p(\mathbf{U}, \mathbf{V} | \mathbf{X}, \lambda) \propto p(\mathbf{U}, \mathbf{V}, \mathbf{X} | \lambda) = p(\mathbf{X} | \mathbf{U}, \mathbf{V}) p(\mathbf{U} | \lambda) p(\mathbf{V} | \lambda)
}

と書ける。これの対数をとったものを {\mathbf{U}, \mathbf{V}} の関数として {J_{\lambda}(\mathbf{U}, \mathbf{V})} と置くと、

{
J_{\lambda}(\mathbf{U}, \mathbf{V}) = \log p(\mathbf{X} | \mathbf{U}, \mathbf{V}) + \log p(\mathbf{U} | \lambda) + \log p(\mathbf{V} | \lambda)
}

となる。各項について見ていこう。まず第一項は以下のように書ける。

{
\log p(\mathbf{X} | \mathbf{U}, \mathbf{V}) = \sum_{i = 1}^{N} \sum_{j = 1}^{M} \log N(x_{ij} | u_{i}^T v_{j}, 1)
}

{
= - \frac{1}{2} \sum_{i = 1}^{N} \sum_{j = 1}^{M} \left( x_{ij} - u_{i}^T v_{j} \right)^2 + (const)
}

{
= - \frac{1}{2} || \mathbf{X} - \mathbf{U}\mathbf{V}^T ||^2_2 + (const)
}

第二項は以下のように書ける。

{
\log p(\mathbf{U} | \lambda) = \sum_{i = 1}^{N} \sum_{k = 1}^{K} \log N(u_{ik} | 0, 1/\lambda)
}

{
= - \frac{\lambda}{2} \sum_{i = 1}^{N} \sum_{k = 1}^{K} u_{ik}^2 + (const)
}

{
= - \frac{\lambda}{2} ||\mathbf{U}||_2^2 + (const)
}

第三項も同様に、以下のように書ける。

{
\log p(\mathbf{V} | \lambda) = - \frac{\lambda}{2} ||\mathbf{V}||_2^2 + (const)
}

以上より、

{
J_{\lambda}(\mathbf{U}, \mathbf{V}) = - \frac{1}{2} || \mathbf{X} - \mathbf{U} \mathbf{V}^T ||^2_2 - \frac{\lambda}{2} ||\mathbf{U}||_2^2 - \frac{\lambda}{2} ||\mathbf{V}||_2^2 + (const)
}

MAP推定を求めるためには {J_{\lambda}} を最大化する {\mathbf{U}, \mathbf{V}} を求めれば良いので、定数項と定数倍を除去して改めて {J_{\lambda}} と置くと、

{
J_{\lambda}(\mathbf{U}, \mathbf{V}) = - || \mathbf{X} - \mathbf{U} \mathbf{V}^T ||^2_2 - \lambda \left( ||\mathbf{U}||_2^2 +
 ||\mathbf{V}||_2^2 \right)
}

となる。ここまでくればわかると思うけど、{J_{\lambda}} を最大化する {\mathbf{U}, \mathbf{V}} を求めることは、{L_{\lambda}} を最小化する {\mathbf{U}, \mathbf{V}} を求めることと等しい。つまり、

上記の確率的生成モデルにおける {\mathbf{U}, \mathbf{V}} の MAP推定と、行列分解の解 {\mathbf{U}, \mathbf{V}} は一致する。

ということが言える。

ちなみに、行列とかグラフとかその辺のデータに対する手法はここでも示したように互いに密接に関わり合っていて、調べていくとすごく面白い。今までまとまった資料ってあんまりなかった(と思う)んだけど、MLPシリーズでいい本が出てるので興味ある人はそれを読んでみると楽しいと思う。

関係データ学習 (機械学習プロフェッショナルシリーズ)

関係データ学習 (機械学習プロフェッショナルシリーズ)

Edward で実装

長くなったけど、行列分解を確率モデルとして解釈できたので、Edward で実装できる。Edward での実装の流れは、

確率モデルの定義 → データを流し込んで推論 → (評価)

という感じになる。

コード(jupyter notebook)はここに置いた。

順を追って説明していく。まず必要なものをインポートする。あと適当に乱数シードを決めておく。

import tensorflow as tf
import edward as ed
import numpy as np

from edward.models import Normal

ed.set_seed(42)

つぎに実験につかうトイデータセットを作る。ここまでは Edward は全く関係ない。U と V を適当に作って、それをもとにデータ X_data と y_data を生成する。実装の都合上、X_data には観測された {x_{ij}} のインデックスのペア (i, j) を格納し、y_data にはそれに対応する値を格納している。また、{\lambda = 1} とした。

なんで普通に行列を作らないでインデックス X_data とそれに対応する値 y_data を別々に作るかと言うと、値が観測された要素のみを扱うことができるため。値が観測されていない要素のインデックスは X_data に含めなければ良い。scipy.sparseで観測されたところにだけ値を入れればいいじゃんと思う人もいるかもしれないけど、それだと値が入っていないところは0が観測されたということを意味する。「0が観測される」ということと、「値が観測されていない」ということは全く違う。

def build_dataset(U, V, n_samples):
    N = U.shape[0]
    M = V.shape[0]
    K = U.shape[1]
    X = []
    sampled = set()
    
    for _ in range(n_samples):
        i = np.random.randint(N)
        j = np.random.randint(N)
        while (i, j) in sampled:
            i = np.random.randint(N)
            j = np.random.randint(N)
        sampled.add((i, j))
        X.append([i, j])
        
    X = np.array(X)
    y = np.random.normal(loc = np.sum(U[X[:, 0]] * V[X[:, 1]], axis=1), scale = 1.0)
    
    return X, y

def train_test_split(X, y, test_ratio):
    n_samples = X.shape[0]
    test_indices = np.random.binomial(1, test_ratio, size = n_samples)
    return X[test_indices==0], y[test_indices==0], X[test_indices==1], y[test_indices==1]
    
N = 30
M = 40
K = 5  # number of latent dimensions
n_samples = 300
U_true = np.random.normal(loc = 0.0, scale = 1.0, size = (N, K))
V_true = np.random.normal(loc = 0.0, scale = 1.0, size = (M, K))
X_data, y_data = build_dataset(U_true,  V_true, n_samples)

test_ratio = 0.1
X_train, y_train, X_test, y_test = train_test_split(X_data, y_data, test_ratio)

ここからようやく Edward を使っていく。まずデータを流し込む “プレースホルダ” を定義する。ここでは X_ph がそれに当たる。X_ph には tf.int32 型でサイズ None x 2 のデータが流し込まれてくるよということを宣言している。None はサイズが事前にはわからないことを意味する。具体的には、推論時には X_data に格納されているインデックスの array を流し込むことになる。

# Prepare a placeholder

X_ph = tf.placeholder(tf.int32, [None, 2])

モデルを定義する。U と V についてはたぶん見たとおりで、平均0、分散1の正規分布から生成されるよということを定義している。y についてはすこしややこしいんだけど、やってることは同じで、平均が tf.reduce_sum(tf.gather(U, X_ph[:, 0]) * tf.gather(V, X_ph[:, 1]), axis=1) で分散が1の正規分布から生成されると定義している。tf.ones_likeで loc と scale が同じサイズになるようにしている。詳細はtensofflowのドキュメントを読んでもらうとして、ざっくり言うと、tf.reduse_sumnp.sumみたいなもんで、tf.gatherはarrayのsliceみたいなもの。まぁつまり上の方で説明したとおりの分布から y (行列の要素の値)が生成されるということを定義している。

# Building the model

U = Normal(loc = tf.zeros([N, K]), scale = tf.ones([N, K]))
V = Normal(loc = tf.zeros([M, K]), scale = tf.ones([M, K]))
y = Normal(loc = tf.reduce_sum(tf.gather(U, X_ph[:, 0]) * tf.gather(V, X_ph[:, 1]), axis=1), scale = tf.ones_like(tf.reduce_sum(tf.gather(U, X_ph[:, 0]) * tf.gather(V, X_ph[:, 1]), axis=1)))

Variational inference のためのモデルを定義する。ここでは真のモデル(U と V)と同様に qU も qV も正規分布とする。平均と分散は適当に正規乱数で初期化しておく。tf.nn.softplusを通すことで、分散が負の値にならないようにしている。

# Building the variational model

qU = Normal(loc = tf.Variable(tf.random_normal([N, K])), scale = tf.nn.softplus(tf.Variable(tf.random_normal([N, K]))))
qV = Normal(loc = tf.Variable(tf.random_normal([M, K])), scale = tf.nn.softplus(tf.Variable(tf.random_normal([M, K]))))

ようやく準備が整ったのでパラメータを推論する。いま求めたいのは U と V なんだけど、それらを直接求めるんじゃなくて、適当な分布(ここでは qU と qV)を用意して、それらが求めたい U と V にできるだけ近づくように、用意した適当な分布のパラメータを求めていく。

詳細はドキュメントなりなんなりを読んでもらうとして、KLqpは求めたい分布(ここでは U と V)と、近似した分布(ここでは qU と qV)の間のKLダイバージェンスが小さくなるようにがんばる。この辺を説明するとまた長くなるのでなんか頑張ってるということにしておこう。

見ればわかると思うけど、KLqpには求めたい分布と近似した分布を対応させた辞書と、プレースホルダと実際のデータを対応させた辞書を渡す。そんであとはrun()するだけ。n_iterは iteration の回数。

# Inference

inference = ed.KLqp({U: qU, V: qV}, data = {X_ph: X_train, y: y_train})
inference.run(n_iter = 1000)

推論も終わったので、最後に評価する。ここではとりあえず適当に推論前と推論後の分布それぞれを使って Mean squared error を計算してみる。y_priorには推論前のyをそのまま入れて、y_postには推論後の qU と qV を使って新たに正規分布を作る。あとは Edward が用意してくれてる評価用の関数を使うだけ。

結果を見ると、正しく推論されてるっぽいことが分かる。

# Evaluation

y_prior = y
y_post = Normal(loc = tf.reduce_sum(tf.gather(qU, X_ph[:, 0]) * tf.gather(qV, X_ph[:, 1]), axis=1), scale = tf.ones_like(tf.reduce_sum(tf.gather(qU, X_ph[:, 0]) * tf.gather(qV, X_ph[:, 1]), axis=1)))

print("Mean squared error (prior): ", ed.evaluate("mean_squared_error", data = {X_ph: X_test, y_prior: y_test}))
print("Mean squared error (posterior): ", ed.evaluate("mean_squared_error", data = {X_ph: X_test, y_post: y_test}))

# 結果
Mean squared error (prior):  6.11875
Mean squared error (posterior):  3.83919

まとめ

試行錯誤しながらいろんなモデルを試したいときにいい感じに使えそう。つまり良い。

研究と開発のはざま

博士を取ってからの3年間はアカデミックで仕事をしていたけど、4月から民間企業に移ることにした。転職するかどうかそうとう悩んだわけだけど、その時に研究についていろいろと考えたのでちょっと書いてみたい。

工学では研究と開発の違いなんて無い

誰もが納得するような明確な違いは無いと思う。あったら教えて欲しい。

実際、ほとんどの研究(論文)が何かを開発しているし、それが世の中の役に立たないと評価されない。逆に、ほとんどのプロダクトには新規性(もしくは他のプロダクトとの差異)がある。確かに、工学において研究と呼ばれるものは開発と呼ばれるものより平均的には基礎的なことをやっているとは思うけど、とはいえ明確な線引は出来ない。

研究を神格化しすぎる風潮があるんじゃないかな。

「研究者です」というと「すごい!」と言われることがよくある。そう言ってもらえるのは嬉しいけど、べつにすごくないよ。99%の研究者は世の中に何のインパクトも与えない研究をやって、何のインパクトも無い論文を書いてる。それってすごいのか? 一方で世の中に多大なインパクトを与えるプロダクトを開発している人だってたくさんいる。要するに何が言いたいかというと、「研究は開発よりもすごい」なんてことは一切ないし、まぁそもそも違いもないよねということ。

研究とはなにか

さて、それなら研究とはなんなのか。それは自分の捉え方次第。自分が研究だと思ったらそれは研究だ。客観的な線引は出来ないんだから、当然そうなる。

じゃあ自分自身にとって研究とはなにかと考えると、「自分がワクワクするような成果をだすこと」だ。ワクワクしなければ意味は無いし、成果が出なければ意味は無い。少なくとも自分にとっては。

わくわくするには、自分の解きたい問題、答えたい質問を見つけなくてはならない。わくわくするにはその問題を解決しなくてはならない。わくわくするにはその解決法について自分が納得しなくてはならない。自分が納得するには、他人も納得させないといけない。というか、他人が納得しない結果に自分が納得するわけがない。他人が納得しないということは、どこかに論理の飛躍がある。

要するに、「自分がわくわくする問題を見つけて、それを他人も納得するような形で解く」というのが自分の中での研究だと思う。

ビッグクエスチョンがないと研究は厳しい

自分がワクワクして、解くのに5年くらいかかる問題がないと研究を続けるのは難しい。

ビッグクエスチョンがないと、適当に問題を見つけて、あるいは誰かからもらって解決するただの問題解決屋になってしまう。問題を解決するだけならどこにいても出来る。すぐに解決できる問題に取り組むのも研究職の利点を活かしていない。すぐに解決できるなら土日に趣味でやれば良い。

研究職の良いところは、自由が多いこと*1

せっかく自由に自分がやることを決められるんだから、面白いことをやりたい。せっかく自由に使える時間があるんだから、自分の好きなことに使いたい。自分がワクワクするような問題にじっくりと取り組む。こんなに楽しい仕事はないと思う。実際、面白い問題を見つけてそれに取り組んでいる時は本当に楽しかった。

ビッグクエスチョンを見つけるには世の中のことを深く知る必要がある。

何が自分にとって面白いのかなんて、結局のところこれまで自分が作ってきた価値観の中でしか考えられない。考えても考えても思いつかないなら、ずっと思いつかないはず。視野が狭い。視野を広げていけばこれまでは思いもつかなかったような面白いことが見つかるかもしれない。視野を広げるって相当難しいと思うけど、少なくとも新しいことをやれば視野が広がって面白いことが見つかる可能性は高くなると思う。

まとめ

自分が心の底から解きたいと思える問題があるなら研究職につくのは素晴らしいと思う。自分にこんなに裁量がある仕事はそうそうない。でもそういう問題が見つかっていないなら一歩踏みとどまったほうが良い。それがないなら研究と開発の間に大きな違いはない。

*1:自由なんてないという研究者の皆さん、大変お疲れ様です。

Author Topic Model の導出と実装

またまた引き続き青いトピックモデル本から。今回は Author Topic Model を導出して実装してみる。とりあえずこのシリーズは一旦今回で最後。

トピックモデル (機械学習プロフェッショナルシリーズ)

トピックモデル (機械学習プロフェッショナルシリーズ)

出典は以下の論文。これまで実装してきたモデルと比べるとずば抜けて有名っぽい。

https://arxiv.org/ftp/arxiv/papers/1207/1207.4169.pdf

Author Topic Model

Author Topic Model (ATM) は文書に付加情報として著者情報が付いているデータのモデリングをするのに使われる*1。一つの文書に複数(一人以上)の著者がいるときに、文書中のそれぞれの単語についてどの著者が書いたのかを推定して割り当てていく。

簡単な例で説明する。スポーツと音楽と政治という3つのトピックについて書かれた文書 d があるとする(あるかな?)。また、文書 d は Alice と Bob によって書かれたとする。Alice はスポーツと音楽について執筆することが多く、Bob は政治について執筆することが多いことは分かっているとしよう。ATM は、文書 d の一つ一つの単語について、Alice と Bob のどちらが担当して書いたのか(生成したのか)を割り当てる。スポーツもしくは音楽っぽい単語なら、その単語は Alice が書いたと判定されるし、政治っぽい単語なら Bob が書いたと判定されることになる。ここまでの説明では Alice と Bob の専門(何についてよく書くか)は分かっているとしたが、実はATMはそれもデータから推定する。

こうすることで何が嬉しいかというと、Bob が政治についてよく書くということがデータから推定されれば、Bob が書いた他の文書も政治について書かれていそうだという推定が出来るようになる。まさに著者情報という付加情報をうまく使ってトピック割り当てを行っている。

ATM のグラフィカルモデルは以下の通り。

f:id:yamaguchiyuto:20170323162459p:plain

LDA と同様に、w は単語、z はトピックID、θはトピック分布、φはトピックごとの単語分布。また、a は著者集合で、y は単語ごとの著者割り当てを表す。S は著者の異なり数(全著者数)で、K はトピック数。LDA では文書ごとにトピック分布があったが、ATM では著者ごとにトピック分布がある。

生成過程は基本的に LDA と同じだけど、y と z の生成だけが異なる。y はその文書の著者集合の中から一様ランダムに選ばれる。つまり、文書 d の著者が {M_d} 人いるとすると、n 番目の単語 {w_{dn}} の著者割り当ては、

{
y_{dn} \sim \frac{1}{M_d}
}

となる。

{y_{dn}=j} が生成されたとして、文書 d の j 番目の著者を {a_{dj}} と表すことにすると、単語 {w_{dn}} のトピック割り当ては、

{
z_{dn} \sim Mul(\cdot | \theta_{a_{dj}})
}

となる。つまり、著者 {a_{dj}} のトピック分布からトピックが生成される。

導出

周辺化ギブスサンプリングのサンプリング式を導出する。ATM は各単語について隠れ変数が2つ({z_{dn}}{y_{dn}})があるので、その2つを同時にサンプリングする。いつも通り右肩の添字はその集合から添字に対応する要素を抜くことを表すとすると、

{
p(z_{dn}=k,y_{dn}=j | w_{dn}=v, A, Z^{dn}, Y^{dn}, W^{dn}, \alpha, \beta)
}

{
\propto \int p(z_{dn}=k,y_{dn}=j, w_{dn}=v, A, Z^{dn}, Y^{dn}, W^{dn}, \Theta, \Phi | \alpha, \beta) d\Theta d\Phi
}

{
\propto \int p(w_{dn}=v | \phi_k) p(\phi_k | W^{dn}, Z^{dn}, \beta) d\phi_k
}

{
\times \int p(z_{dn}=k | \theta_{a_{dj}}) p(\theta_{a_{dj}} | Z^{dn}, Y^{dn}, \alpha) d\theta_{a_{dj}}
}

{
= \frac{n_{kv}^{dn} + \beta}{\sum_v (n_{kv}^{dn} + \beta)} \frac{n_{a_{dj}k}^{dn} + \alpha}{\sum_k (n_{a_{dj}k}^{dn} + \alpha)}
}

となる。添字がごちゃごちゃしててわかりづらいんだけど、 {n_{kv}^{dn}} は、{z_{dn}} を除いて、トピック k が単語 v に割り当てられた数を表し、{n_{a_{dj}k}^{dn}}{z_{dn}} を除いて、著者 {a_{dj}} が割り当てられた単語に対してトピック k が割り当てられた数を表す。

実装

z と y を同時にサンプリングするので、そこの実装だけ注意。各イテレーションごとにサンプルされた Z と Y を用いて尤度を計算して、それが最も大きかった Z と Y を保存しておく。

Author Topic Model

実験

架空の著者「山田」と「鈴木」が書いた文書という想定の人工データを使って実験してみた。山田はコンピュータに関するトピックについて多く書いていて、鈴木は食べ物に関するトピックについて多く書いている。

結果を見ると、 “apple” という単語はコンピュータトピックにも食べ物トピックにも両方登場するんだけど、著者情報を使うことによって正しくトピックが割り当てられている。また、最後の文書は山田と鈴木の共著ということになっているが、文書中の単語はコンピュータに関するものだけなので、割り当てられた著者も山田だけということになった。鈴木のギフトオーサーシップが発覚した!

ATM experiment

*1:著者情報じゃなくても適用は可能。

Noisy Correspondence Topic Model の導出と実装

さらに引き続き青いトピックモデル本から。今回はノイズ有り対応トピックモデル (Noisy Correspondence Topic Model; NCTM) を導出して実装する。

トピックモデル (機械学習プロフェッショナルシリーズ)

トピックモデル (機械学習プロフェッショナルシリーズ)

出典は以下の論文。

http://www.kecl.ntt.co.jp/as/members/iwata/nips2009.pdf

Noisy Correspondence Topic Model

このモデルは Correspondence Topic Model (CTM) の拡張になっていて、CTM と同様に付加情報を考慮しながら文書のモデリングが出来る。どう拡張されているかというと、付加情報の中からノイズとノイズでないものを切り分けることが出来る。

付加情報のノイズについて、論文中でも取り上げられているソーシャルブックマークの例を使って説明する。ソーシャルブックマークでは、Webページ(文書)にタグ(付加情報)を付けることが出来る。タグの中には、文書のトピックを表すようなもの(例えば技術系のブログ記事に “TECH” とか)と、文書のトピックを表さないようなもの(例えば、"あとで読む" とか)の2種類がある。前者は文書のモデリングに役立ちそうだけど、後者は役立ちそうにないので、このモデルではノイズと呼ぶ。

ノイズでない付加情報は特定のトピックに偏って出てくるけど、ノイズである付加情報は偏り無く、全てのトピックにまんべんなく出てくる。つまり、「ノイズとは特定のトピックに偏り無く出現する付加情報」であると、(ある意味)定義していることになる。これ重要。

で、どう切り分けるかだけど、各付加情報について、それがノイズかノイズでないかを表す確率変数 {r \in {0,1}} を導入する。r が 0 のときは対応する付加情報はノイズで、1 のときはノイズでないとする。そんでその r が 0 なのか 1 なのかをデータから推定する。この r をスイッチ変数と呼ぶ。

グラフィカルモデルは以下。

f:id:yamaguchiyuto:20170322064202p:plain

CTM から変わったのはスイッチ変数の部分だけで、r は各付加情報ごとにある。生成過程も CTM とほぼおなじだけど、r は λ をパラメータとするベルヌーイ分布から生成され、λ は η をパラメータとするベータ分布から生成される。一件ごちゃごちゃしてるけど、拡張の系譜をたどると理解しやすい。

なぜスイッチ変数をデータから推定できるか?

トピックモデルを1ミリもわかってなかったときはなんでこの変数をデータから推定できるのか全くわからなかったんだけど、今は1ミリくらいは分かっているので直感的に説明してみる。トピックモデル(LDAから派生したモデル群)は基本的には単語のハードクラスタリングを行っている。言い換えると、「それぞれの単語は特定のトピック(クラスタ)から出現した」という仮定を入れている。 つまり、全てのトピックからまんべんなく出現する単語というのはモデル上有りえない。

じゃあどうするか。上記のクラスタリングをする前に、その単語が「特定のトピックから出現した」のか、「トピックに関係なく出現した」のか判別する。前者と判別されたらその単語がどのトピックから出現したのかをクラスタリングする。後者と判別されたらしない。当然、ノイズは「トピックに関係なく出現した」としたほうが、つまり r=0 としたほうが、モデルのもっともらしさは高くなるので、そのように推定されることになる。

導出

周辺化ギブスサンプリングのサンプリング式を導出する。Z のサンプリングは CTM と同じなので省略。前回と同様に右肩の添字はその集合から添字に対応する要素を抜くことを表す。

Y のサンプリング

Y のサンプリングも CTM とほぼ同じなので導出は省略するけど、 文書 d の m 番目の付加情報のトピック {y_{dm}} をサンプリングするときに、{r_{dm}} が 0 か 1 かでサンプリングの式が変わる。

{r_{dm}=1} のとき、つまりノイズでないときは CTM の時と同じ。

{
p(y_{dm} = k | x_{dm} = u, r_{dm} = 1, Z, W, X^{dm}, Y^{dm}, R^{dm}, \alpha, \beta, \gamma, \eta)
}

{
\propto n_{dk} \frac{m^{dm}_{ku}+\gamma}{\sum_u (m^{dm}_{ku}+\gamma)}
}

{n_{dk}} は文書 d の単語のうちトピック k が割り当てられたものの数を表し、{m^{dm}_{ku}}{x_{dm}} を除いて、トピック k が割り当てられた付加情報 u の数を表す。

{r_{dm}=0} のとき、つまりノイズであるときは CTM の時と少し違う。

{
p(y_{dm} = k | x_{dm} = u, r_{dm} = 0, Z, W, X^{dm}, Y^{dm}, R^{dm}, \alpha, \beta, \gamma, \eta)
}

{
\propto n_{dk}
}

つまり、付加情報 {x_{dm}} がノイズであるときは、その文書においてトピック k が割り当てられた単語数のみが考慮される。ちなみに、{x_{dm}} がノイズであるときは {y_{dm}} が他の変数のサンプリング式に影響をあたえることはないので、サンプリングする必要はない。

R のサンプリング

{r_{dm}} のサンプリングをするときは、{r_{dm}=0} の事後分布と {r_{dm}=1} の事後分布を両方計算して2つの合計が 1 になるように正規化すれば良い。

{
p(r_{dm} = 1 | x_{dm} = u, y_{dm} = k, Z, W, X^{dm}, Y^{dm}, R^{dm}, \alpha, \beta, \gamma, \eta)
}

{
\propto \int p(r_{dm} = 1, x_{dm} = u, y_{dm} = k, Z, W, X^{dm}, Y^{dm}, R^{dm}, \Psi, \lambda | \alpha, \beta, \gamma, \eta) d\Psi d\lambda
}

{
\propto \int p(x_{dm} = u | r_{dm} = 1, y_{dm} = k, \psi_k)p(\psi_k|X^{dm}, Y^{dm}, R^{dm}, \gamma) d\psi_k
}

{
\times \int p(r_{dm} = 1 | \lambda )p(\lambda|R^{dm}, \eta) d\lambda
}

{
= \int \psi_{ku} p(\psi_k | X^{dm}, Y^{dm}, R^{dm}, \gamma)d\psi_k \times \int \lambda p(\lambda|R^{dm},\eta) d\lambda
}

{
\propto \frac{M^{dm}_{ku} + \gamma}{\sum_u (M^{dm}_{ku} + \gamma)} (M^{dm}_1 + \eta)
}

{M^{dm}_{ku}}{x_{dm}} を除いて、トピック k が付加情報 u に割り当てられた回数を表し、{M^{dm}_1}{x_{dm}} を除いて、ノイズでない付加情報の数を表す。

また同様に、

{
p(r_{dm} = 0 | x_{dm} = u, y_{dm} = k, Z, W, X^{dm}, Y^{dm}, R^{dm}, \alpha, \beta, \gamma, \eta)
}

{
\propto \frac{M^{dm}_{0u} + \gamma}{\sum_u (M^{dm}_{0u} + \gamma)} (M^{dm}_0 + \eta)
}

{M^{dm}_{0u}}{x_{dm}} を除いて、付加情報 u がノイズであると判別された回数を表し、{M^{dm}_0}{x_{dm}} を除いて、ノイズである付加情報の数を表す。

実装

{r_{dm}}が 0 か 1 かでサンプリングとかの処理を変える必要があるので、その辺を注意。

Noisy Correspondence Topic Model

実験

CTM と NCTM を比較した。ノイズと判別されるべきの付加情報 “TOREAD” を加えて、それがちゃんと NCTM によってノイズと判別されるかを見てみる。結果を見ると、NCTM は “TOREAD” をしっかりとノイズと判別するため、他の単語及び付加情報にうまくトピックを割り当てることが出来ている。一方で CTM は “TOREAD” にトピック 0 (ipadとかと同じトピック) を割り当ててしまっているせいで、4番目の文書の単語 “apple” にトピック 0 を割り当ててしまっている。

NCTM experiment