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

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

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

まとめ

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