Factorization MachinesでThompson Samplingしてみる



 ## Logistic Contextual Sampling
[[Chapelle,2011]](https://proceedings.neurips.cc/paper/2011/file/e53a0a2978c28872a4505bdb51db06dc-Paper.pdf)を参照
- ロジスティック回帰 + Thompson Sampling
- 回帰係数に事前分布を導入
- 尤度はベルヌーイ分布
- 係数の事後分布をラプラス近似で求める
- 係数をサンプリング(Thompson Sampling)して、それを使って予測値を計算する
- 最も最大の予測値を採用

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1793316/02bbc88b-3f79-67f3-132b-4ab4572827a4.png)


 ## contexual bandit (Factorization Machines)


FM版で考えてみる.
- データ数が$N$ ,特徴量次元が$D$, 相互作用項$v$の次元は$k$
- $ X = \lbrace \mathbf{x}_1, \mathbf{x}_2, ... , \mathbf{x}_N \rbrace$
- $ \mathbf{x}_i =  \lbrace \mathbf{x} _{i1}, \mathbf{x} _{i2}, ... , \mathbf{x} _{iD} \rbrace $
- $ \mathbf{V} =  \lbrace \mathbf{v}_1, \mathbf{v}_2, ... , \mathbf{v}_D \rbrace $
- $ \mathbf{v}_i =  \lbrace \mathbf{v} _{i1}, \mathbf{v} _{i2}, ... , \mathbf{v} _{ik} \rbrace $

一般的なFMの予測式
```math
FM(\mathbf{x}_i) = w_0 + <\mathbf{w}, \mathbf{x}_i> + \sum_{m = 1}^D\sum_{n = m+1}^D<\mathbf{v}_m, \mathbf{v}_n>x_{im}x_{in}
```
事前分布を導入
```math
\displaylines{
p(w_i) \sim \textit{N}(0, q_w^{-1}) \\
p(v_{ik}) \sim \textit{N}(0, q_v^{-1})
}
```
尤度
$\mu_i = FM(\mathbf{x}_i)$として、
```math
P(\mathbf{y} | \mathbf{X},\mathbf{w},\mathbf{V}) = \prod_i^{N}\mu_i^{y_i}(1 - \mu_i)^{1-y_i}
```
よって事後分布は
```math
\begin{align}
P(\mathbf{w},\mathbf{V}| \mathbf{X},y)
&\propto P(\mathbf{y} | \mathbf{X},\mathbf{w},\mathbf{V})P(\mathbf{w},\mathbf{V})\\
&=\prod \mu_i^{y_i}(1 - \mu_i)^{1-y_i}\prod N(0, q_w)\prod N(0, q_v)
\end{align}
```

正規加項は計算不可能なため近似手法に頼るしかない(MCMC, 変分推論, ラプラス近似)

## ラプラス近似
任意の正規化項のわからない分布を無理やり正規分布に近似する手法
確率分布$P(\boldsymbol{\theta})$の対数を$\boldsymbol{\theta}_0$の周りで二次の項までテイラー展開
```math
lnP(\boldsymbol{\theta}) \approx lnP(\boldsymbol{\theta}_0) + \nabla lnP(\boldsymbol{\theta})^T|_{\boldsymbol{\theta} = \boldsymbol{\theta}_0}(\boldsymbol{\theta} - \boldsymbol{\theta}_0) + \frac{1}{2}(\boldsymbol{\theta} - \boldsymbol{\theta}_0)^T\nabla ^2 lnP(\boldsymbol{\theta}) |_{\boldsymbol{\theta} = \boldsymbol{\theta}_0}(\boldsymbol{\theta} - \boldsymbol{\theta}_0)
```
ここで、この$\theta_0$は$\nabla lnP(\boldsymbol{\theta})$を0にする点とすると、
```math
lnP(\boldsymbol{\theta}) \approx lnP(\boldsymbol{\theta}_0) + \frac{1}{2}(\boldsymbol{\theta} - \boldsymbol{\theta}_0)^TH(\boldsymbol{\theta} - \boldsymbol{\theta}_0) \\
,H = \nabla ^2 lnP(\boldsymbol{\theta}) |_{\boldsymbol{\theta} = \boldsymbol{\theta}_0}
```
```math
P(\boldsymbol{\theta}) = exp(lnP(\boldsymbol{\theta})) \propto exp\left(\frac{1}{2}(\boldsymbol{\theta} - \boldsymbol{\theta}_0)^TH(\boldsymbol{\theta} - \boldsymbol{\theta}_0)\right)\\
N(\boldsymbol{\theta} | \boldsymbol{\mu},\Sigma) \propto exp\left(-\frac{1}{2}(\boldsymbol{\theta} - \boldsymbol{\mu})^T\Sigma^{-1}\boldsymbol{\theta} - \boldsymbol{\mu})\right) \\

```
このようにして事後分布を$N( \boldsymbol{\theta}_0,  -H^{-1})$で近似する.
今回は$ \boldsymbol{\theta} = \lbrace \mathbf{w}, \mathbf{V}\rbrace$であり、近似分布の平均$\boldsymbol{\theta}_0$は事後分布を最大にする$\boldsymbol{\theta} $, つまりMAP推定量でありこれは損失関数をクロスエントロピーとした、L2ノルム正則化付きの推定量である.


## FM + ラプラス近似
近似したい確率分布は$P(\mathbf{w},\mathbf{V}| \mathbf{X},y)$
この分布をラプラス近似していく.

分散共分散は対角行列と仮定をおいている.


対数変換すると, $\mu_i = FM(\mathbf{x}_i)$として
```math
lnP(\mathbf{w},\mathbf{V}| \mathbf{X},y) \propto  \sum_i^N y_iln\mu_i +  (1-y_i)ln(1 - \mu_i) - \sum_i^D \frac{q_w}{2}w_i^2 - \sum_i^D  \sum_j^k \frac{q_v}{2} v_{ij}^2
```



この対数事後確率の勾配を0にする点、つまり最大にする点だが、それは数式を見ると正則化項付きのクロスエントロピー損失関数の最小化と等価であることわかる. 

### 学習時
1. 勾配法などで事後確率が最大となる点を探索し、平均とする <- これは正則化項付きFMの学習と等価
1. 共分散行列を計算 <- これは解析的に可能(FMが二階微分可能なため)
     ```math
     \displaylines{
     q_i^{-1} = \frac{\partial^2 logP(\mathbf{w}, \mathbf{V} | \mathbf{X},\mathbf{y})}{\partial w_{m}^2} = q_w + \sum_i^N \mu_i(1 - \mu_i)x_{im}^2\\
     q_i^{-1} = \frac{\partial^2 logP(\mathbf{w}, \mathbf{V} | \mathbf{X},\mathbf{y})}{\partial v_{mn}^2} = q_v + \sum_i^N \mu_i(1 - \mu_i)\left(\left(\sum_j^D v_{jn}x_{ij} \right)x_{im} - v_{mn}x_{im}^2\right)^2
     }
    ```

既存のコードに対して、更新式を2つ追加すれば良いので楽です。

### 予測時
1. モデルから与えられる平均、共分散行列からパラメータをサンプリング
1. FMの式でそのまま計算

分散共分散行列の非対角成分をゼロとしているため、正規分布からのサンプリングにコレスキー分解などをする必要がなくサンプリングにかかる計算量も$O(kD)$で済みます。Factorizatoin Machinesのメリットである線形の計算量を落とすことなく推論が可能です。


## 参考

https://tech.nextroll.com/blog/data-science/2017/03/06/thompson-sampling-bayesian-factorization-machines.html