読者です 読者をやめる 読者になる 読者になる

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

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

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 とすればよい

クロネッカー積とvec作用素とRoth's column lemma

クロネッカー積とvec作用素は見た目簡単なんだけど、各要素のインデックスを書き下すと頭の中こんがらがってわけわからなくなるから一旦整理した。インデックスに関する記事ってほとんどないね。あとそれに関連して Roth's column lemma っていうのが便利なのでちょっと紹介する。

まとめ

クロネッカー積とは

クロネッカー積 - Wikipedia

下の式を見ればひと目で分かるけど、 IxJ 行列 {A} と KxL 行列 {B} から IKxJL 行列を作る。

{
A \otimes B = \left(
\begin{array}{ccc}
a_{11} B & \cdots & a_{1J} B \\
\vdots & \ddots & \vdots \\
a_{I1} B & \cdots & a_{IJ} B
\end{array}
\right)
}

Wikipediaの例を見ると分かりやすい。

https://wikimedia.org/api/rest_v1/media/math/render/svg/130f0defabfa409a56ff89937d40f1656fcb108f

クロネッカー積のインデックス

クロネッカー積の定義は一目見ればすぐ分かるのに、各要素の定義を見ると途端にめんどくさくなる。

{
\left[A \otimes B\right]_{(i-1)K+k,(j-1)L+l} = a_{ij}b_{kl}
}

行のインデックスを見ると、K進数みたいになってるのが分かる。 {A} の行のインデックスが1つ増えるためには、{B} の行のインデックスを一回ずつ辿らなきゃいけないって考えるとそうなることが分かる。列についても同じ。

実験

上式のインデキシングが本当に合ってるのか numpy.kron を使って確かめる。

gist6cf39c77df3bf44fe6ed5aaf98c343cc

vec作用素

vec作用素 - Wikipedia

vec作用素は、以下のように IxJ 行列 {A} の列ベクトルを全て一列に並べて1つの大きな IJ 次元列ベクトルを作る。

{
vec(A) = \left(
\begin{array}{c}
a_{1} \\
a_{2} \\
\vdots \\
a_{J}
\end{array}
\right)
}

例えば、以下のようになる

{
vec\left(
\begin{array}{ccc}
1 & 3 & 2 \\
0 & 7 & 4
\end{array}
\right) = \left(
\begin{array}{c}
1 \\
0 \\
3 \\
7 \\
2 \\
4
\end{array}
\right)
}

vec作用素のインデックス

こっちはクロネッカー積よりは簡単。

{
\left[vec(A)\right]_{(i-1)J+j} = a_{ij}
}

クロネッカー積のときと同じで、J進数みたいになってる。各列ベクトルがJ次元なので、{a_{ij}}{i} が増えるためには {j} を全部辿らなきゃならない。

Roth's column lemma

これが成り立つ。

{
vec(XYZ) = (Z^{T} \otimes X)vec(Y)
}

証明

{X} をIxJ行列、{Y} をJxK行列、{Z} をKxL行列とし、

{W = XYZ}

と置く。{W} は IxL 行列。

{W} の各要素は(とうぜん)以下のように書ける。

{
w_{il} = \sum_{pq} x_{ip}y_{pq}z_{ql}
}

クロネッカー積の要素を考えると

{
\left[Z^{T} \otimes X\right]_{(l-1)I+i,(k-1)J+j} = z_{kl}x_{ij}
}

vec作用素の要素を考えると

{
\left[vec(Y)\right]_{(j-1)K+k} = y_{jk}
}

{
\left[vec(W)\right]_{(i-1)L+l} = w_{il}
}

なので、以下の行列ベクトル積の要素を考えると、

{
\left[(Z^{T} \otimes X)vec(Y)\right]_{(i-1)L+l}
}

{
= \sum_{pq} \left[Z^{T} \otimes X\right]_{(l-1)I+i,(q-1)J+p} \left[vec(Y)\right]_{(p-1)J+q}
}

{
= \sum_{pq}z_{ql} x_{ip} y_{pq}
}

{
= \sum_{pq}x_{ip} y_{pq}z_{ql}
}

{
= w_{il} = \left[vec(W)\right]_{(i-1)L+l}
}

{W=XYZ} だから、証明完了。