[論文メモ]Factorization Machines

## 読んだ論文
Steffen Rendle, <i>Factorization Machineas</i>, in ICDM, 2010
Link: https://www.csie.ntu.edu.tw/~b97053/paper/Rendle2010FM.pdf

## TL;DR
交差項をfactorized parametersを用いてモデル化した手法

## この手法でやりたきこと
- SVMが得意としないsparseなデータセットを用いた予測タスク
- 他のfactorized model（ex. Matrix Factorization, SVD++）が得意としない通常の予測タスク

## 仕組み

$y \in \mathbb{R} \rightarrow T, T = \(0,1\)$の分類タスクを考える
Facotorization Machines（以下FM）は以下のように定式化される。

```math
\hat{y}(\boldsymbol{x}) := w_0 + \sum_{i=1}^n w_i x_i + \sum_{i=1}^n \sum_{j=i+1}^n <\boldsymbol{v}_i, \boldsymbol{v}_j> x_i x_j

```

それぞれの記号の定義は以下の通り

- $w_0$ ... バイアス
- $\boldsymbol{w} \in \mathbb{R}^n$ ... パラメータ
- $\boldsymbol{v}_i$ ... i番目の特徴量に対するパラメータ（k次元）   
　 　\* \<v, v\>はドット積を表す  

FMの中でも一番キモなのはこのドット積を用いている二項目。
交差項一つ一つに対して重みを設定しようとするとそもそも大変な上、レコメンドのようにsparseなデータセットの場合データセット中頑張っても1回ぐらいしかでてこずSVMのような手法での学習は難しい。
そこで交差項に対する重みをfactorizeすることで交差項間の類似性を取り入れ、未知の入力パターンに対しても対応できるようにした（と私は理解した）。

論文中でも触れられているがpositive definiteな行列$W$について、$k$が十分に大きい場合$W = V \cdot V^T$となるような$k$次元のベクトル$V$が存在するとのこと。
ただしレコメンドの場合はデータセットがsparseなためあまり大きくならないらしい（あとで試す）。

## 計算上の工夫

パラメータをfactorizeしている二項目をそのまま素直に計算しようとすると$O(kn^2)$の計算コストがかかってしまう。
しかしながらきちんと式を変形することで$O(kn)$までコストを減らすことができる。ということで式変形する。

```math
\begin{align}

\sum_{i=1}^n \sum_{j=i+1}^n <\boldsymbol{v}_i, \boldsymbol{v}_j> x_i x_j &= \frac{1}{2}\sum_{i=1}^n \sum_{j=1}^n <\boldsymbol{v}_i, \boldsymbol{v}_j> x_i x_j - \frac{1}{2}\sum_{i=1}^n \sum_{j=1}^n <\boldsymbol{v}_i, \boldsymbol{v}_i> x_i x_i　\\
&= \frac{1}{2}\sum_{i=1}^n \sum_{j=1}^n \sum_{f=1}^k v_{i,f} v_{j,f} x_i x_j - \frac{1}{2}\sum_{i=1}^n \sum_{j=1}^n \sum_{f=1}^k v_{i,f} v_{i,f} x_i x_i　\\
&= \frac{1}{2} \sum_{f=1}^k \Bigl\{ \bigl(\sum_{i=1}^n v_{i,f}x_i \bigl) \bigl(\sum_{j=1}^n v_{j,f}x_j \bigl) - \sum_{i=1} v_{i,f}^2 x_i^2 \Bigr\} \\
&= \frac{1}{2} \sum_{f=1}^k \Bigl\{ \bigl(\sum_{i=1}^n v_{i,f}x_i \bigl)^2  - \sum_{i=1} v_{i,f}^2 x_i^2 \Bigr\}

\end{align}
```
一つ目の変形はもともと$n \times n$の行列のうち半分（からさらに対角成分除いた部分）だけ計算していたものをいったん$n \times n$全てを計算し、対角成分を除いた後半分にするという操作。
二つ目の変形は$\boldsymbol{v}$を$k$について展開しただけ。ここから一項目を$i, j$それぞれについてまとめ、二項目を二乗の形に直す。
三つ目の式の一項目をよく見ると$i$にしても$j$にしても$n$個という同じ数の要素があるのでこの部分もまとめられる。
ということで最後の式が導かれる。

## SVMとの比較
ここでは結論だけをまとめる

- linear kernelの場合FMにおける$d=1$のモデルと同じ
- polinomial kernelの場合交差項のサンプルが極端に少なく十分な学習が望めない（論文中ではlinear kernelと大差ない言ってる）

## 実装

今回は[fastFM](https://github.com/ibayer/fastFM)というライブラリを使用し実装。
コードの全体は[こちらのGitHubリポジトリ](https://github.com/wmichi/recommend_notebooks)にあげている[^1]。

#### データセットの準備
データセットには[MovieLens20M](https://grouplens.org/datasets/movielens/20m/)を用いている。
二値分類の問題にするためネガティブサンプリングした上でデータセットを以下のような形に変更[^2]。

```python
import numpy as np
import pandas as pd
from fastFM import sgd
from sklearn.feature_extraction import DictVectorizer
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

df = pd.read_csv('preprocessed_dataset.csv')
df.head() 
```

| user | item | rating |
|--:|--:|--:|
| 52224 | 74738 | -1 |
| 3179 | 5596 | 1 |
| 121962 | 3383 | -1 |
| 121626 | 111628 | -1 |
| 36358 | 4752 | -1 |

{0,1}ではなく{+1,-1}で作成している点に注意。
このデータセットを{user: xxx, item: zzz}という辞書のリストに直す。
yは一次元にreshapeする（.reshape(-1,1)とかすると実行中にコケる）。

```python
X_list = []

for row in df.itertuples(index=False, name=None):
    X_list.append({"user": str(int(row[0])), "item": str(int(row[1]))})
y = np.array(df.iloc[:,-1]).reshape(-1,)
```

#### 学習

データセットを投入用にDictVectorizerで変換。
このとき全てのtrain/testに全てのインデックスが入っていないと予測のときにコケる。
ということでここでX全体を作成してからtrain/testを分割しなくてはいけない。

```python
v = DictVectorizer()
X = v.fit_transform(X_list)

X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.2, random_state=0)
```

学習はいつものsklearn的なやり方でOK。
先の式でいう$k$はここでいうrank。
先程までの式では出てこなかったが正則化項も設けられる。

```python
fm = sgd.FMClassification(n_iter=param_dict['n_iter'], 
                          init_stdev=param_dict['init_stdev'], 
                          l2_reg_w=param_dict['l2_reg_w'],
                          l2_reg_V=param_dict['l2_reg_v'], 
                          rank=param_dict['rank'], 
                          step_size=param_dict['step_size'])
fm.fit(X_train, y_train)

fm.fit(X_train, y_train)
```

param_dictの部分、それぞれのパラメータを以下のように設定している。
rankだけ{10,20,40,80}の4通り試行。

```
n_iter: 1000000
init_stdev: 0.1
l2_reg_w: 0.1
l2_reg_V: 0.01
step_size: 0.1

--
rank: {10, 20, 40, 80}
```

そして予測する

```python
y_pred = fm.predict(X_test)
accuracy = accuracy_score(y_test, np.round(y_pred))

rankname = 'rank-' + str(param_dict['rank'])
print('Accuracy:')
print(rankname, '{:.4f}'.format(accuracy))
```

```
>>> Accuracy:
>>> rank-10: 0.7538
>>> rank-20: 0.7538
>>> rank-40: 0.7540
>>> rank-80: 0.7538
```

rank=40のときが結果として一番よかったけどそこまで大きな差ではないな、というのは正直な印象。
他のパラメータ（特に正則化項の強さ）ももっといじると面白そう。

ということで今回はここまで。
お粗末！

### 参考資料
Steffen Rendle, <i>Factorization machines with libFM</i>, in ACM, Trans.  Intell.  Syst.  Technol., 2012
Link: https://www.csie.ntu.edu.tw/~b97053/paper/Rendle2010FM.pdf

##### （自分用メモ。ただのメモ）
factorizeするパラメータの次元が$k=3$の場合二項目は以下のようなイメージになる。

```math
\begin{align}
\sum_{i=1}^n \sum_{j=u+1}^n <\boldsymbol{v}_i, \boldsymbol{v}_j> x_i x_j &= \sum_{i=1}^n \sum_{j=u+1}^n \sum_{f=1}^k v_{i,f} v_{j,f} x_i x_j \\
&= \sum_{i=1}^n \sum_{j=u+1}^n \bigl\{v_{i,f}^1 v_{j,f}^1 + v_{i,f}^2 v_{j,f}^2 + v_{i,f}^3 v_{j,f}^3 \bigr\} x_i x_j
\end{align}
```
いってしまえばEmbeddingと発想が同じなのだろうか。

ちなみにlibfm使った方が早いと参考資料の論文で言ってたが、個人的にMovieLens20MぐらいならLibfm使わなくても十分早く感じた。

[^1]: Train/Eval/Test分割の都合上ここで記載しているコードとは一部異なっている（本筋には影響ない）。
[^2]: 少し具体的に説明すると適当な$(user, item)$の組み合わせを乱数で作成し、そのうち元のデータセットに存在する$(user, item)$の組み合わせのものを除外したデータセットを負例として使用。それに対してもとのデータセットを正例として使用したというもの。