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

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

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