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

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

Joint Topic Modelを実装した

LDAの簡単な拡張になっている Joint Topic Model を実装した。青いトピックモデル本で紹介されてた。この本はいろんなモデルが載ってるのでいいね。

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

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

実装したあとで気づいたけど既にnzw君が実装して実験してたのでこちらも参考に。

nzw0301.github.io

Joint Topic Model

Joint Topic Model (JTM) はLDAとほとんど同じなんだけど、文書に付加情報(カテゴリとか)がついてる場合、それも使うことができる。

どんな付加情報を扱えるかというと、基本的にはカテゴリ変数だけ。生成過程を見ると分かるように、付加情報の生成にはカテゴリカル分布が使われているので、何かしらの数値(レビュー文書のレーティングとか)が付加情報としてついていてもこのモデルではうまく扱えない(理論上は)。

このモデルが応用されている例として、対訳文書のモデリングをしてる以下の論文がある。

http://dirichlet.net/pdf/mimno09polylingual.pdf

この論文では文書として言語Aの文書、付加情報として言語Bの文書をモデルに与えて、両方の情報を用いてトピック推定をしている。JTMは複数の付加情報を同時にモデリングできるので、言語C, D, E, …の文書を同時にモデルに与えることも出来る。

実装

周辺化ギブスサンプリングで実装。ほとんどLDAと同じ。Z(トピック割り当て)のサンプリングの式がちょっと違うだけ。

Joint Topic Models

実験

JTMとLDAを比較してみた。PCもしくはFOODというカテゴリが与えられている文書がうまくモデリング出来るかを比較した。結果を見ると、JTMはカテゴリ情報をうまく使って、全く同じ単語が出現している文書にもうまくトピックを割り当てることが出来ている。LDAではカテゴリ情報が与えられていないので当然そういうことは出来ない。

JTM experiments

Robust Large-Scale Machine Learning in the Cloud [KDD'16] を読んだ

KDD16で発表されてた論文。著者はかの有名なFactorization Machinesの人。Googleに行ってたのね。いままでとはちょっと違う研究をしてるように感じる。

論文はここから読める。

www.kdd.org

勉強会で紹介したので念のため、その時のスライドはこちら。

一言まとめ

一般化線形モデルの学習を Coordinate Descent で、めっちゃスケールさせるよ。

概要

Coordinate Descent (CD) はシングルマシン上では収束が早いアルゴリズムとして知られてるんだけど、分散には向かない。なぜなら、アルゴリズムの特性上、分散させると1つのワーカーに割り当てられる work load が小さくなってしまって、オーバヘッドを考えると分散効果が小さいから。CDはパラメータの更新を次元ごとに 逐次 繰り返していくことで収束するアルゴリズムなので、分散させるには基本的にはデータパラレルにすることになる。でもパラメータの1つの次元を更新するときの処理はそんなに重くないので、分散させるメリットがあまりないということ。

そこでこの論文では1つの次元ごとじゃなくて、次元の集合(ブロック)ごとに更新するのをくりかえす Scalable Coordinate Descent (SCD) というアルゴリズムを提案する。そうすることで分散させて1つのワーカーに割り当てる work load をある程度大きくすることが出来る。で、ブロックごとに更新しても収束するのか? という話だけど、パラメータの更新時にこれまでの値と更新後の値の間で line search してやると収束するらしい。なぜかは知らん(書いてない)。さらに、特徴次元をうまいことブロックに分割しれやれば、CDと全く同じ解が得られることも示した。

で、この SCD を Google Cloud 上に実装したらしい。論文中では実装時の細かいところの説明がされていた。この辺もさらっと書いてあるけどいろいろ試行錯誤して大変だったんだろなと思った。ポイントは、VM Preemptionがある環境を想定しているところなのかな? AWS の Spot instance とかがそれなんだけど、優先度が低いとみなされた VM(=worker) が勝手に強制終了されて、空いたリソースが別の優先度の高いVMに割り当てられてしまうという環境。

実際にGoogleが持ってる広告のデータを使って実装したらしい。サンプル数は1兆(!)で、特徴数は10億(!)のスケールのデータらしい。Google恐るべし。ちなみにこれはNetflix prizeの1万倍の規模らしい。結果はここでは省略するけどスケールしたとのこと。

無限潜在特徴モデルを実装した

引き続きノンパラベイズ。今回はノンパラベイズ本の7.4節で説明されてる無限潜在特徴モデル(Infinite latent feature model; ILFM)を実装した。

無限潜在特徴モデル

一言で言えば、潜在次元に無限次元を仮定する行列分解。

与えられたデータ行列 {Y \in \mathbb{R}^{N \times D}} をバイナリ行列 {Z \in \{0,1\}^{N \times K}} と特徴行列 {X \in \mathbb{R}^{K \times D}} の積に分解する。

バイナリ行列の要素 {z_{n,k}} はサンプル n が k 番目の潜在特徴を持つか持たないかを表している。k番目の潜在特徴は {X} のk番目の行ベクトル {x_k^T \in \mathbb{R}^D} で表されている。つまり、あるサンプルのデータベクトル {y_n} は無限個ある潜在特徴のうちのいくつかの和になっている。

通常の行列分解では分解後の行列の潜在次元数 K を予め決めなきゃいけないんだけど、このモデルはノンパラベイズなので決めなくて良い。データから勝手に決まる。さすがノンパラベイズ

モデルの説明とパラメータ推定の説明は上の本に詳しく書いてあるので割愛。

実装

特に注意すべき点はないと思うけど、強いて言うならサンプリングの事後分布を計算するときにunderflowしないように対数のまま正規化してやるくらい。

Infinite Latent Feature Model

実験

alphaを小さく設定してもどうしても真のKより大きな潜在次元数が推定されてしまう。何か間違ってるのかこういうもんなのか判断ができない。相対誤差は十分(?)小さくなってる。

ILFM experiment

無限混合ガウスモデルを実装した

ノンパラベイズ面白いね。佐藤一誠先生のノンパラメトリックベイズの本を読んで自分なりに理解できたので実装してみた。本読んで理解して、自分で導出して、実装・実験するの本当に重要。定着度がぜんぜん違う。

無限混合ガウスモデルは上の本の5.2節で説明されてるモデルで、ディリクレ過程混合ガウスモデル(Dirichlet Process Gaussian Mixture Model; DPGMM)とも呼ぶらしい。面倒くさいので以下DPGMMと書く。

DPGMM

普通の混合ガウスだと分割数(クラスタの数)K を事前に決めなきゃいけないんだけど、これが面倒くさいし、よほどそのデータに詳しいか超能力でも無い限り、いろんな K の値を試してみないといけない。面倒くさい。

そこでDPGMMの出番。予め K を決めなくてもデータから推定してくれる。便利。

感覚的にちょっと説明すると、あるデータ点のクラスタ割り当てを決めるとき、普通の混合ガウスだったらそのデータ点を生成したっぽいクラスタ(ガウシアン)に割り当てるんだけど、DPGMMだとどのクラスタからも生成されたっぽくないときは新しいクラスタが作られる。このクラスタ割り当てを何度も何度もやる(ギブスサンプリング)ことで、最終的に適切な数のクラスタが作られ、全てのデータ点が適切なクラスタに割り当てられる。

詳細は上の本を読むと良い。佐藤先生の本は本当に超わかりやすい。

実装

めっちゃ遅いけど実装した。本に書いてあるように、共分散行列は対角行列と仮定した。

Dirichlet Process Gaussian Mixture Model

実験

真のクラスタ数は3で、分散は全てのクラスタで1。クラスタ間が互いにかなり離れてるので簡単な設定になってる。

結果を見ると、クラスタ数を指定してないのに、クラスタが3つ作られて、全てのデータ点が適切なクラスタに割り当てられていることが分かる。

DPGMM experiment

Tucker分解の導出と実装

CP分解の次はTucker分解を導出して実装する。丁寧にTucker分解の導出を説明してる文献(Web含め)が全然なかったので、自分で書く。CP分解についてはある程度知ってる前提とする。CP分解についてはこちらから。

yamaguchiyuto.hatenablog.com

まとめ

  • Tucker分解とは
  • ALSでTucker分解の更新式の導出
  • PythonでTucker分解を実装
  • 人工データを使って実験

Tucker分解とは

Tucker分解は、テンソルを1つのテンソル(コアテンソルと呼ぶ)と、それぞれのモードに対して一つずつの行列に分解する。

f:id:yamaguchiyuto:20161127144753p:plain

上の図の例では、もとのテンソルのサイズは IxJxK だけど、これをコアテンソルのサイズの RxSxT (R<=I, S<=J, T<=K) まで小さくしている。また、あとで説明するけど、行列 U、V、W は全て直行行列となるように分解する。このコアテンソル {\mathcal{G}} と3つの直交行列を使ってもとのテンソル {\mathcal{X}} を近似する。

これを式で書くと、

{
x_{ijk} \simeq \sum_{rst} g_{rst} u_{ir} v_{js} w_{kt}
}

という分解をしていることになる。

この分解が何かに似ていると思った人はお目が高い。実は、コアテンソルが対角要素が全部1の対角テンソルのとき、Tucker分解はCP分解と等しくなる。つまり、CP分解はモデル的にはTucker分解の特殊ケースということになる。

ALSでTucker分解の更新式の導出

上で説明したように、Tucker分解ではテンソル {\mathcal{X}} を1つのコアテンソル {G} と、(テンソルが3階の場合は)3つの直交行列 {U}{V}{W}に分解する。分解の方針はCP分解のときと同じで、ALS (Alternating Least Square) で行く。

さて、以下のコスト関数を最小にすることを考える。

{
L = \frac{1}{2} \sum_{ijk} \left(x_{ijk} - \sum_{rst} g_{rst} u_{ir} v_{js} w_{kt} \right)^{2}
}

これはただ単にテンソル {X} の全要素について近似誤差の和の二乗を取っているだけ。ただし、問題を簡単にするために {U}{V}{W} は直交行列であるとする。

L を {g_{abc}}微分すると

{
\frac{\partial L}{\partial g_{abc}} = - \sum_{ijk} \left(x_{ijk} - \sum_{rst}g_{rst}u_{ir}v_{js}w_{kt}\right)\left(u_{ia}v_{jb}w_{kc}\right)
}

{
= - \sum_{ijk}x_{ijk}u_{ia}v_{jb}w_{kc} +\sum_{ijk}\sum_{rst}g_{rst}u_{ir}v_{js}w_{kt}u_{ia}v_{jb}w_{kc}
}

{
= - \sum_{ijk}x_{ijk}u_{ia}v_{jb}w_{kc} +\sum_{rst}g_{rst}\sum_{i}u_{ir}u_{ia}\sum_{j}v_{js}v_{jb}\sum_{k}w_{kt}w_{kc}
}

ここで、{U} は直交行列だから、{r\neq a} のとき {\sum_{i} u_{ir} u_{ia} = 0}、また {r = a} のとき、{\sum_{i} u_{ir} u_{ia} = 1} となる。よって、

{
= - \sum_{ijk}x_{ijk}u_{ia}v_{jb}w_{kc} + g_{abc}
}

これを 0 と置いて {g_{abc}} について解くと、

{
g_{abc} = \sum_{ijk}x_{ijk}u_{ia}v_{jb}w_{kc}
}

となる。これを★とする。

★を L に代入するが、その前にまず L をちょっと変形すると

{
L = \frac{1}{2} \sum_{ijk} \left(x_{ijk} - \sum_{rst} g_{rst} u_{ir} v_{js} w_{kt} \right)^{2}
}

{
= \frac{1}{2} \sum_{ijk} \left(x_{ijk}^{2} - 2x_{ijk}\sum_{rst} g_{rst} u_{ir} v_{js} w_{kt} + \left(\sum_{rst}g_{rst}u_{ir}v_{js}w_{kt}\right)^{2}\right)
}

{
= \frac{1}{2} \sum_{ijk} x_{ijk}^{2} - \sum_{ijk} x_{ijk}\sum_{rst} g_{rst} u_{ir} v_{js} w_{kt} + \frac{1}{2} \sum_{ijk} \left(\sum_{rst}g_{rst}u_{ir}v_{js}w_{kt}\right)^{2}
}

ここで、第二項を変形して★を代入すると、

{
\sum_{ijk} x_{ijk}\sum_{rst} g_{rst} u_{ir} v_{js} w_{kt} = \sum_{rst} g_{rst} \sum_{ijk} x_{ijk} u_{ir} v_{js} w_{kt}
}

{
= \sum_{rst} g_{rst} g_{rst} = \sum_{rst} g_{rst}^{2}
}

また第三項に★を代入すると、

{
\sum_{ijk} \left(\sum_{rst}g_{rst}u_{ir}v_{js}w_{kt}\right)^{2} = \sum_{ijk} \left(\sum_{rst}\sum_{abc}x_{abc}u_{ar}v_{bs}w_{ct}u_{ir}v_{js}w_{kt}\right)^{2}
}

{
= \sum_{ijk} \left(\sum_{abc}x_{abc}\sum_{r}u_{ar}u_{ir}\sum_{s}v_{bs}v_{js}\sum_{t}w_{ct}w_{kt}\right)^{2}
}

ここでまた {U} が直交行列なので、{a = i} かつ {b = j} かつ {c = k} のときの項のみが残って、

{
= \sum_{ijk} x_{ijk}^{2}
}

となる。

まとめると、Lは以下のように書ける。

{
L = \sum_{ijk} x_{ijk}^{2} - \sum_{rst} g_{rst}^{2}
}

結局、L を最小化するという問題は、{U}{V}{W}が直行という制約のもとで

{\sum_{rst} g_{rst}^{2} = \sum_{rst}\left(\sum_{ijk}x_{ijk} u_{ir}v_{js}w_{kt}\right)^{2} }

を最大化するという問題に帰着した(やっとスタート地点に立った・・・)。これを L' と置こう。

さて、ALSなので、{V}{W} を固定して、L' を最大化する {U} を求める。

今後の見通しを良くするために L' をちょっと変形する、

{
L' = \sum_{rst}\left(\sum_{ijk}x_{ijk} u_{ir}v_{js}w_{kt}\right)^{2}
}

{
= \sum_{rst}\left(\sum_{i} u_{ir} \sum_{jk} x_{ijk}v_{js}w_{kt}\right)^{2}
}

ここで {\sum_{jk} x_{ijk}v_{js}w_{kt} = x^{VW}_{ist}} と置くと、

{
= \sum_{rst}\left(\sum_{i} u_{ir} x^{VW}_{ist}\right)^{2}
}

また、添字 s と t をまとめて q と書くと*1

{
= \sum_{rq}\left(\sum_{i} u_{ir} x^{VW}_{iq}\right)^{2} = \sum_{rq}\left(u_{r}^{T} x^{VW}_{q}\right)^{2}
}

ここで、{x^{VW}_{q} = (x^{VW}_{1q}, \cdots, x^{VW}_{Iq})^{T}}{u_{r} = (u_{1r}, \cdots, u_{Ir})^{T}} とした。すべての {u_{r}} は独立に考えられるから、ある {u_{r}} について {u_{r}^{T}u_{r}=1} の制約のもとで L' を最大化する。

{
L'(r) = \sum_{q}\left(u_{r}^{T} x^{VW}_{q}\right)^{2} - \lambda_{r} (u_{r}^{T}u_{r}-1)
}

とおいて、{u_{r}}微分すると、

{
\frac{\partial L'(r)}{\partial u_{r}} = 2\sum_{q}\left(u_{r}^{T} x^{VW}_{q}\right)x^{VW}_{q} - 2\lambda_{r} u_{r}
}

{
= 2 \sum_{q}\left(x^{VW}_{q} (x^{VW}_{q})^{T}\right)u_{r} - 2\lambda_{r} u_{r}
}

{X^{VW} = (x^{VW}_{1}, \cdots, x^{VW}_{Q})} と置くと、

{
= 2X^{VW}(X^{VW})^{T}u_{r} - 2\lambda_{r} u_{r}
}

これを 0 と置くと、

{
X^{VW}(X^{VW})^{T}u_{r} = \lambda_{r} u_{r}
}

つまり、{u_{r}} は 行列 {X^{VW}(X^{VW})^{T}}固有ベクトル、また行列 {X^{VW}} の左特異ベクトルということになる。

この {u_{r}} を全ての r について互いに直交するように求めれば良いから、結局行列 {X^{VW}}特異値分解して得られる左特異ベクトルが答えとなる。

長かった・・・。

{V}{W} についても同様に求まる。まとめると、

  1. {U^{(0)}}{V^{(0)}}{W^{(0)}} を適当に初期化する
  2. {k \leftarrow 0}
  3. {X^{V^{(k)}W^{(k)}}}特異値分解し、左特異ベクトルを {U^{(k+1)}} とする
  4. {X^{W^{(k)}U^{(k+1)}}}特異値分解し、左特異ベクトルを {V^{(k+1)}} とする
  5. {X^{U^{(k+1)}V^{(k+1)}}}特異値分解し、左特異ベクトルを {W^{(k+1)}} とする
  6. {k \leftarrow k+1}
  7. 収束するまで 2-6 を繰り返す。
  8. ★により {\mathcal{G}} を求める。

実装

scipy.sparse は疎テンソルを扱えないので、データテンソルX[(indices)]=value という辞書で与えた。indicesはインデックスが入ったタプル。

gist74f1026096be91ba0a147eeaa540326c

実験

平均0、標準偏差1の正規分布でコアテンソル {\mathcal{G}} と直交行列 U、V、W を初期化し、データテンソル {\mathcal{X}} を作って、平均0、標準偏差0.1の正規ノイズを加えた。

Tucker分解が真のコアテンソル {\mathcal{G}} と直交行列 U、V、W を復元できるか試した。

max_iter=5 としたけど、一回で(ほぼ)収束した。収束するまで特異値分解し続けるのがHOOI、一回だけ特異値分解するのがHOSVDと呼ばれてるけど、この例ではHOSVDで十分っぽい。

学習した後に残差(の平均)を見てみたけど十分小さいのでOK。

gist86035070f286f595f3f1078ed2002e3a

*1:q = s*T + t とすればよい