備忘録

主に機械学習、統計、競技プログラミングに関する備忘録

Regression-based EMの導出

はじめに

Learning to rankにおけるポジションバイアスの求め方の一種であるRegression-based EM[Wang18]の導出についてまとめる。 より具体的には論文中のEMアルゴリズムの更新式(2)のあたりを詳しく追う。

定義

基本的には論文の表記に従う。

  • $q$ : クエリ
  • $d$ : 文書
  • $C$ : ユーザが文書をクリックしたか
  • $R$ : ユーザが文書に興味があるか
  • $E$ : ユーザが文書を確認したか
  • $k$ : 文書の位置
  • $\mathcal{L} = \left\{ \left(c, q, d, k\right)\right\}$ : ユーザのログの集合

Position Bias Model

Posigion Bias Model (PBM) ではクリックしたかを表す確率変数 $C$ が潜在変数 $E, R$ を通じて以下のように生成されると仮定する: $$ P(C = 1 \mid q, d, k) = P(E = 1 \mid k) \cdot P(R = 1 \mid q, d) $$

すなわち、$C=1$ となるのは $E=R=1$ の場合に限る。

第一項及び第二項をそれぞれ以下のように定義しておく:

$$ \begin{align} \theta_k &= P(E = 1 \mid k), \\ \gamma_{q, d} &= P(R = 1 \mid q, d). \end{align} $$

$\gamma_{q, d}$ は通常適当な特徴量 $x_{q,d}$ と予測器 $f$ を用いて $\gamma_{q,d} = f(x_{q, d})$ と計算される。しかし、ここでは機械学習モデルによる予測には立ち入らず、PBMの尤度最大化による $\gamma_{q, d}$ の最適化を考えるものとする。

PBMにおけるパラメータは $\left\{ \theta_k \right\}, \left\{ \gamma_{q, d}\right\}$ であり、最適なパラメータは対数尤度

$$ \begin{align} \log P(\mathcal{L}) &= \sum_{(c, q, d, k) \in \mathcal{L}} c \log P(C = 1 \mid q, d, k) + (1 - c) \log P(C = 0 \mid q, d, k) \\ &= \sum_{(c, q, d, k) \in \mathcal{L}} c \log \theta_k \gamma_{q, d} + (1 - c) \log (1 - \theta_k \gamma_{q, d} ) \end{align} $$

を最大化するものである。この目的関数を直接最適化するのは困難なため、EMアルゴリズムを用いて最適なパラメータを反復法で求めることになる。

EMアルゴリズムによるパラメータ推定

EMアルゴリズムの概要

観測可能な確率変数 $X$ と観測不可能な潜在変数 $Z$ が確率密度関数 $p(X, Y \mid \xi)$ に従ってるとする。 $Z$ を観測出来ないという設定の下で素朴に $\xi$ を推定しようと思えば対数尤度

$$ \log p(X^n \mid \theta) = \log \sum_{Z^n} p(X^n, Z^n \mid \xi) $$

を最大化すればよいが、この操作は一般には優しくない。EMアルゴリズムでは以下の操作によって $\hat{\xi}^{(0)}, \hat{\xi}^{(1)}, \dots$ を生成する。このように生成された列 $\hat{\xi}^{(t)}$ は局所最適解に収束することが示される。

  • Eステップ : $p(Z \mid X, \hat{\xi}^{(t)})$ を計算する
  • Mステップ : $\hat{\xi}^{(t+1)} = \mathrm{arg~max}~Q(\xi \mid \hat{\xi}^{(t)})$ を計算する ただし、 $Q$ は以下のような関数である: $$ \begin{align} Q(\xi \mid \hat{\xi}^{(t)}) &= \mathbb{E}_{Z^n \mid X^n, \hat{\xi}^{(t)}} \left[ \log p(X^n, Z^n \mid \xi) \right] \\ &= \sum_{Z^n} p(Z^n \mid X^n, \hat{\xi}^{(t)}) \log p(X^n, Z^n \mid \xi) \end{align} $$

$\theta_{k}, \gamma_{q, d}$ の推定

EMアルゴリズムの設定において $Z = (E, R), \xi = \left( \left\{ \theta_{k} \right\}, \left\{ \gamma_{q, d} \right\} \right)$ として更新式の導出を行う。

Eステップ

観測データで条件付けた下での潜在変数の条件付き確率 $P(E, R \mid C, q, d, k)$ を求める。 これは変数の定義と設定した確率モデルから直ちに求まる。

Mステップ

論文中ではMステップにおいて $E$ と $R$ それぞれの周辺確率について $Q$ 関数を計算し、$\theta_{k}$ と $\gamma_{q, d}$ の更新式を求めているが、その部分が若干天下り的だったので詳細を確認する。 まず今回の設定で $Q$ 関数を書き下すと以下のようになる:

$$ \begin{align} Q(\xi \mid \hat{\xi}^{(t)}) = \sum_{(c, q, d, k) \in \mathcal{L}} \sum_{(E, R)}P(E, R \mid c, q, d, k) \left( c \log P(C = 1, E, R \mid q, d, k) + (1 - c) \log P(C=0, E, R \mid q, d, k) \right) \end{align} $$

ここで $\log P(C,E,R \mid q, d, k)$ は

$$ \log P(C, E, R \mid q, d, k) = \log P(C \mid E, R) + \log P(E \mid k) + \log P(R \mid q, d) $$

と分解できる。パラメータ $\xi$ のうち、$\theta_{k}$ は $\log P(E \mid k)$に、$\gamma_{q, d}$ は $\log P(R \mid q, d)$ にそれぞれ依存することに着目して、$Q$ 関数も以下のように分解する:

$$ \begin{align} Q(\xi \mid \hat{\xi}^{(t)}) &= Q_{\theta}(\theta \mid \hat{\xi}^{(t)}) + Q_{\gamma}(\gamma \mid \hat{\xi}^{(t)}) + \mathrm{const.}, \\ Q_{\theta}(\theta \mid \hat{\xi}^{(t)}) &= \sum_{(c, q, d, k) \in \mathcal{L}} \sum_{E}P(E \mid c, q, d, k) \left( c \log P(E \mid k) + (1 - c) \log P (E \mid k) \right), \\ Q_{\gamma}(\gamma \mid \hat{\xi}^{(t)}) &= \sum_{(c, q, d, k) \in \mathcal{L}} \sum_{R}P(R \mid c, q, d, k) \left( c \log P(R \mid q, d) + (1 - c) \log P(R \mid q, d) \right). \end{align} $$

次に $\theta$ の第 $k$ 要素 $\theta_{k}$ の更新式を求める。$Q_{\theta}$ を $\theta_{k}$ で偏微分した式

$$ \frac{\partial}{\partial \theta_{k}} Q_{\theta}(\theta \mid \hat{\xi}^{(t)}) = 0 $$

を満たす $\theta_{k}$ を求めれば良い。

$$ \begin{align} \frac{\partial}{\partial \theta_{k}} Q_{\theta}(\theta \mid \hat{\xi}^{(t)}) &= \frac{\partial}{\partial \theta_{k}} \sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \sum_{E}P(E \mid c, q, d, k) \left( c \log P(E \mid k) + (1 - c) \log P (E \mid k) \right) \\ &= \frac{\partial}{\partial \theta_{k}} \sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \left\{P(E=1 \mid c, q, d, k) \left( c \log P(E=1 \mid k) + (1 - c) \log P (E=1 \mid k) \right) + P(E=0 \mid c, q, d, k) \left( c \log P(E=0 \mid k) + (1 - c) \log P (E = 0 \mid k) \right) \right\} \\ &= \frac{\partial}{\partial \theta_{k}} \sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \left\{ P(E=1 \mid c, q, d, k) \left(c \log \theta_{k} + (1-c)\log \theta_{k}\right) + P(E=0 \mid c, q, d, k) \left( c \log (1 - \theta_{k} )+ (1-c) \log (1 - \theta_{k}) \right)\right\} \\ &= \frac{\partial}{\partial \theta_{k}} \sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \left\{ P(E=1 \mid c, q, d, k) \left(c \log \theta_{k} + (1-c)\log \theta_{k}\right) + P(E=0 \mid c, q, d, k) (1-c) \log (1 - \theta_{k})\right\} (\because~ c=0~\mathrm{when}~E = 0)\\ &= 0 \cdots(1) \end{align} $$

以下簡単のため $ p = P(E=1 \mid c, q, d, k)$ とおいて変形を進める。

$$ \begin{align} (1) &\iff \sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \frac{\partial}{\partial \theta_{k}} \left\{ p\left(c \log \theta_{k} + (1-c)\log \theta_{k}\right) + (1-p)(1-c)\log (1 - \theta_{k}) \right\} = 0 \\ &\iff \sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \left\{ p\left(\frac{c}{\theta_{k}} + \frac{1-c}{\theta_{k}}\right) - \frac{(1-p)(1-c)}{1 - \theta_{k}} \right\} = 0 \\ &\iff \sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \left\{ c\frac{p}{\theta_{k}} + (1-c) \left( \frac{p}{\theta_{k}} - \frac{1-p}{1-\theta_{k}} \right) \right\} = 0\\ &\iff \sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \left\{ cp(1 - \theta_{k}) + (1-c) \left( p(1 - \theta_{k}) - (1-p) \theta_{k} \right)\right\} = 0 \\ &\iff \theta_{k} \sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \left(cp + (1-c)\right) = \sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \left(cp + (1-c)p\right) \end{align} $$

ここで、$c=1$ のとき $p = P(E=1 \mid c, q, d, k) = 1$より $cp + 1 - c = 1\cdot 1 + 1 - 1 = 1$, $c=0$ のとき $cp + 1 - c = 1$ より、左辺は $\theta_{k} \sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k')$ となる。右辺に関しても $cp = c$ となることから $\sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \left(c + (1-c)p\right)$ と整理出来る。以上より、$\theta_{k}$ に関する更新式

$$ \theta_{k}^{(t+1)} = \frac{\sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k') \left(c + (1-c)P(E=1 \mid c, q, d, k)\right)}{\sum_{(c, q, d, k') \in \mathcal{L}} \mathbb{1}(k = k')} $$

が得られた。$\gamma_{q, d}$ の更新式も同様に求めることが出来る。

参考

[Wang18] Wang, X., Golbandi, N., Bendersky, M., Metzler, D., & Najork, M. (2018, February). Position bias estimation for unbiased learning to rank in personal search. In Proceedings of the Eleventh ACM International Conference on Web Search and Data Mining (pp. 610-618).

Robust Bayesian inference via coarsening

Robust Bayesian inference via coarsening についてまとめる。

背景

通常のベイズ推定ではあるデータの生成過程を仮定することで推論を行うが、これらの仮定が実際の生成過程と差異がある場合、僅かな差異であっても推定結果に大きな違いをもたらすことが知られている。
この論文ではこういった差異に対してロバストな結果を与えることが出来る coarsened posterior (または c-posterior)というものを提案している。

Toy Example

$X_1, \cdots, X_n$ がそれぞれ $\mathrm{Bernoulli}(\theta)$ からの i.i.d. サンプルであるとする。次に $H_{0} : \theta=1/2$ versus $H_{1} : \theta \ne 1/2$ の検定を考える。
各仮設に対する事前分布は一様に $\Pi(H_0) = \Pi(H_1) = 1/2$ であり、$H_1$ の下での $\theta$ の事前分布は $\theta \mid H_1 \sim \mathrm{Uniform}(0, 1)$ に従うものとする。
このとき、観測データ $x_{1:n}$ が実際に $\mathrm{Bernoulli}(\theta)$ からのサンプルであれば、事後分布は正しい "答え" を導くことが保証されてる。すなわち、

$$ \Pi(H_0 \mid x_{1:n}) \overset{\mathrm{a.s.}}{\rightarrow} \mathbb{1}(\theta = 1/2) $$

が成り立つ。ここで$\mathbb{1}(\cdot)$ は $E$ が真であれば $\mathbb{1}(E) = 1$、そうでないとき $\mathbb{1}(E) = 0$ となるような指示関数である。
しかし、現実的な問題として、仮定した生成過程が厳密に正しいケースは稀である。例えば $\theta = 0.51, 0.56$ と設定して事後分布を見てみると、従来の事後分布では $n$ が大きくなるにつれ所望の値から大きく乖離することが確認できる。

f:id:invariant:20200815211407p:plain

上の例の場合 $\theta = 0.5$ と $\theta=0.51$ の差も検出したいのであれば問題はないのだが、あるモデルを仮定したもとで予測を行いたい場合、真のモデルと僅かに乖離しただけで予測精度が大幅に悪化するのは避けたいため、このように使用するモデルに対してロバストな事後分布を構成することには大きなご利益がある。

coarsened posterior

本論文の筆者はモデルの設計者が決めたパラメトリックモデルにおける確率分布 $P_{\theta_I}$ と実際に観測されるデータの経験分布 $\hat{P}_{x_{1:n}}$ の間の距離(ダイバージェンス)に関する条件付け行う。具体的な定義を行うため、記号の準備を行う。
理想的なモデルを生成する分布(idealized distribution)を含むパラメトリックモデルを $\left\{P_\theta : \theta \in \Theta\right\}$ とする。特に、idealized distributionを実現するパラメータを $\theta_I$ とおく。$P_{\theta_I}$ からの現実には観測されないサンプルを $X_1, \cdots, X_n \in \mathcal{X}$ とする。次に、実際に観測されるサンプルを $x_1, \cdots, x_n \in \mathcal{X}$、経験分布を $\hat{P}_{x_{1:n}} = \frac{1}{n} \sum_{i=1}^{n} \delta_{x_{i}}$ とおく。最後に、2つの確率分布の距離を測る関数として $d(\cdot, \cdot)$ を用意する。通常のベイズ推定では $X_{1:n} = x_{1:n}$ であると仮定することを確認しておく。
上でも述べたように、一般に観測値が設計者の設定した分布に正確に従うことはありえないので、$X_{1:n}$ の仮想的な経験分布 $\hat{P}_{X_{1:n}}$ と実際の観測値の経験分布 $\hat{P}_{x_{1:n}}$ 間の距離に一定の"遊び"が入るのが自然と考えられる。
例えば両確率分布間に距離 $r$ までを許容すると決めたならば、条件 $d_n(X_{1:n}, x_{1:n}) := d(\hat{P}_{X_{1:n}}, \hat{P}_{x_{1:n}}) < r$ の下での $\theta$ の事後分布を構成するのが良さそうである。また、$r$の具体的な値を決めることは困難なので、確率変数化し、 $R \sim H$ に従うものとする。
以上をまとめるとc-posteriorは以下のように定義される。

Definition( c-posterior ) $$ \begin{align} \Pi(d\theta \mid d_n(X_{1:n}, x_{1:n})) &\propto \Pi(d\theta) \mathbb{P}(d_n(X_{1:n}, x_{1:n}) < R \mid \theta) \\ &=: \Pi(d\theta) \int_{\mathcal{X}^n} G(d_{n}(x'_{1:n}, x_{1:n})) P_\theta^{n}(dx'_{1:n}) \end{align} $$

ここで、$G(\cdot)$ は $G(r) = \mathbb{P}(R > r)$ を満たす関数である。 この積分はまず $x'_{1:n}$ が $P^{n}_\theta$ からサンプルされ、次に確率変数 $R$ が $d_n$ を上回る確率を畳み込んだものと読むとよい。これはある $\theta$ の下で観測値 $x_{1:n}$ とそれなりに近いサンプルがどれだけ生成され得るかを表すような値と解釈出来る。
このままでは具体性にかけるので、 $d, H$ を具体的に定めて計算を進める。$d$ をKLダイバージェンス、$R \sim \mathrm{Exp}(\alpha) \propto \exp(-\alpha R)$ とすると、c-posteriorは $G$ に指数分布の密度関数を代入することで、

$$ \begin{align} \Pi(d\theta \mid d_{n}(X_{1:n}, x_{1:n}) < R) &\propto \mathbb{E}\left[ \exp(-\alpha D(\hat{p}_{x_{1:n}} \mid\mid \hat{p}_{X_{1:n}} ) ) \mid \theta\right] \Pi(d\theta) \\ &\propto \exp(\zeta_n \sum_{i=1}^{n} \log p_{\theta} (x_i) ) \Pi(d\theta) \\ &= \Pi(d\theta) \prod_{i=1}^{n} p_{\theta}(x_i)^{\zeta_n}, \\ \zeta_n &= \frac{1}{1 + n / \alpha} \end{align} $$

が得られる。ただし、第1行から第2行へはそれなりの飛躍あるため、詳細はAppendixを参照。ここまで来れば通常の尤度に逆温度 $1/\xi_i$ を導入した事後分布であると解釈でき、MCMCなどの各種計算も可能であるが分かる。既存の数値計算アルゴリズムをそのまま流用してロバストな事後分布が構成出来るというのがこの論文の貢献である。

理論的結果

事後分布が真の生成分布 $P_{o}$ の変化に対してロバストであることは、事後分布が(何らかの位相において)$P_o$ に対して連続的であると言える。これはラフに言えば $P_o$ の微小な変化に対して事後分布も微小な変化をすることを要請する。
論文中では以下の定理で c-posterior が $P_o$ に関して連続であることを示している。

回帰への応用

論文中では時系列モデルと回帰モデルへの応用を紹介している。ここでは後者の結果について紹介しよう。
モデルとしてはSpike-slab modelによるスパースな回帰モデルを考える。

$$ \begin{align} W &\sim \mathrm{Beta}(r, s) \\ \beta_{j} &\sim \mathcal{N}(0, 1/L_0) \mathrm{with probability} W, \mathrm{otherwise} \beta_j = 0, j=1,\cdots,p \\ \lambda &\sim \mathrm{Gamma}(a, b)\\ Y_i \mid \beta, \lambda &\sim \mathcal{N}(\beta^\top x, 1/\lambda), i=0,\cdots,n \end{align} $$

これに対して、以下の分布に従うデータ $(x_i, y_i)$ を生成させる。

$$ y_i = \beta_{01} + \beta_{02}(x_{i2} + \frac{1}{16} x_{i2}^2) + \epsilon_{i} $$

ただし、 $\beta_{01} = -1, \beta_{02} = 4, \epsilon \sim \mathcal{N}(0, 1)$ とする。

f:id:invariant:20200816190015p:plain

このデータにして事後分布 $\pi(\beta, \lambda \mid y) \propto \pi(\beta, \lambda) p(y \mid \beta, \lambda)^{\zeta_{n}}$ を計算し、非ゼロ成分の個数と $\beta$ の各成分の累積分布関数を求めると以下のようになる。

f:id:invariant:20200816190500p:plain

まず非ゼロ成分の個数について、通常のベイズ事後分布ではサンプルサイズが大きくなるにつれ3~5と予測するようになるが、c-posteriorではサンプルサイズが増えても正しい非ゼロ成分の個数を当てていることが分かる。
さらに、事後分布の最頻値を見ると、$\beta_{1}$ が$-0.7$ 付近、$\beta_2$ が $3.5 \sim 4.0$ の間と、実際の値に近い値を取ってることが分かる。

まとめと感想

本論文で提案された事後分布は非常にシンプルな形で、既存のプログラムにも組み込みやすそうである。 数値実験の結果についても、真の分布がわからない場合により妥当な推論が可能になることが確認出来た。時間があれば実際に手元でデータを動かしてみて、その有用性を確認してみたい。
一点気になることとしては、回帰などの予測が伴うタスクの際に、c-posteriorを使うことでどの程度予測精度が改善されるかが分からないので、調べてみたい。

AtCoder Beginner Contest 151 F - Enclose All

問題へのリンク

問題文

平面上の $N$ 個の点 $(x_i, y_i)$ が与えられます。

これら全てを内部または周上に含む円の半径の最小値を求めてください。

自分の解法

条件を満たす円の中心の座標を $(s, t)$ とする。中心 $(s, t)$ の円のうち、与えられた点を全て満たすものの半径の2乗は$\max_{i=1,\cdots,N} (x_i-s)^{2} + (y_i - t)^{2} $ となる。
この半径の2乗を最小にする $(s, t)$ を求めれば良いので、以下の最適化問題を解けば良い。

$$ \min_{s, t} \max_{i=1, \cdots, N} (x_i - s)^{2} + (y_i - t)^{2} $$

ここで、各 $i$ について $(x_i - s)^{2} + (y_i - t)^{2}$ は $(s, t)$ に関する凸関数であり、2つの凸関数の $\max$ もまた凸関数となるので、$ \max_{i=1, \cdots, N} (x_i - s)^{2} + (y_i - t)^{2}$ も$(s, t)$ に関する凸関数となる。
よって $(s, t)$ それぞれについて勾配法でパラメータを更新していけば最適解が求まる。ただし、ステップサイズを減少させる必要があることに留意する。

提出コード

#include <iostream>
#include <fstream>
#include <cstdio>
#include <cmath>
#include <vector>
#include <string>
#include <set>
#include <map>
#include <stack>
#include <queue>
#include <bitset>
#include <algorithm>
#include <complex>
#include <array>
#include <iomanip>
using namespace std;

#define REP(i,n) for(int i=0; i<n; ++i)
#define FOR(i,a,b) for(int i=a; i<=b; ++i)
#define FORR(i,a,b) for (int i=a; i>=b; --i)
#define ALL(c) (c).begin(), (c).end()

typedef long long ll;
typedef vector<int> VI;
typedef vector<ll> VL;
typedef vector<long double> VD;
typedef vector<VI> VVI;
typedef vector<VL> VVL;
typedef vector<VD> VVD;
typedef pair<int,int> P;
typedef pair<ll,ll> PL;

template<typename T> void chmin(T &a, T b) { if (a > b) a = b; }
template<typename T> void chmax(T &a, T b) { if (a < b) a = b; }

int in() { int x; scanf("%d", &x); return x; }
ll lin() { ll x; scanf("%lld", &x); return x; }
#define INF 1LL<<60

int N;
vector< pair<double, double> > points;

pair<int, double> calc_max(double s, double t) {
    double r_max = 0;
    int idx = -1;
    REP(i, N) {
        if(r_max < sqrt( pow(points[i].first - s, 2) + pow(points[i].second - t, 2) )) {
            r_max = sqrt( pow(points[i].first - s, 2) + pow(points[i].second - t, 2) );
            idx = i;
        }
    }
    return make_pair(idx, r_max);
}

int main() {
    cin >> N;
    points.resize(N);
    REP(i, N) {
        cin >> points[i].first >> points[i].second;
    }

    double s = 50.0, t = 50.0, lr;
    REP(itr_outer, 50000) {
        if(itr_outer < 1000)
            lr = 1.0 / sqrt(1.0 + itr_outer);
        else if(itr_outer < 5000)
            lr = 1.0 / pow(1 + itr_outer, 1);
        else
            lr = 1.0 / pow(1.1, sqrt(1 + itr_outer));
        pair<int, double> ret = calc_max(s, t);
        double gs = 2 * (s - points[ret.first].first);
        double gt = 2 * (t - points[ret.first].second);
        s -= lr * gs;
        t -= lr * gt;
    }
    double ans = calc_max(s, t).second;
    printf("%.12f\n", ans);

    return 0;
}

感想

解説の解き方がどうすれば思い付くのか謎。
ステップ幅の調整で魔チューニングを行ってるが、もう少しシンプルに出来るかもしれない。

Item2Vec

embeddingを使って特徴量を生成したいなと思ってたところ、丁度いい論文を見つけたのでメモしておく。
Oren Barkan, Noam Koenigstein (2016). Item2Vec: Neural Item Embedding for Collaborative Filtering

概要

大量の文章から各単語のベクトル表現を獲得するWord2Vecを参考に、推薦システムなどにおけるアイテムのベクトル表現を獲得する。

Skip-gram

まずはWord2VecにおけるSkip-gramについて説明する。 $n$ 個の文が与えられるとして、$i$番目の文の $j$ 番目の単語を $w_{ij}$、文の長さを $L_i$ と書くことにする。また、ある単語の周辺の単語集合をコンテキストと呼ぶことにする。 単語は予め整数値のidに振り直されてるものとする。
Skip-gramとはある単語に注目した時に、その単語のコンテキストを推測するモデルのことを指す。 例えば She sells seashells by the seashore. という文において sell に注目した時、前後の単語である shesell を当てるようなタスクを行うことになる。
$w_{ij}$ の前後 $c$ 個の単語からなるコンテキストを $U_{ij}$ と書くことにする。
いま、$i$ 番目の文に注目したときに、コンテキストを当てることは以下の対数尤度を最大化することと等価であることが分かる。

$$ \frac{1}{L_i} \sum_{j=1}^{L_i} \sum_{w' \in U_{ij}} \log P(w' \mid w_{ij}) $$

元々の目標は各単語のベクトル表現を得ることだったので、注目する単語 $w$ に対応するベクトルを $u_{w} \in U \subset \mathbb{R}^{d}$、コンテキストの単語 $w'$ に対応するベクトルを $v_{w'} \in V \subset \mathbb{R}^{d}$ とし、 以下のモデルを採用する。

$$ P(w' \mid w) = \frac{\exp(u_{w}^{\top}v_{w'})}{\sum_{w{''} \in I_{W}}u_{w}^{\top}v_{w{''}}} $$

$I_{W}$ は単語の添字全体からなる集合である。
ここで、ある単語 $w$ のベクトルは、注目する単語として扱われる場合 $u_{w}$、コンテキストとして扱われる場合 $v_{w}$ となることに注意したい。 例えば $u_i$ と $v_i$ を区別しない場合、上の文章では she を表すベクトルと sell のベクトルの内積が $1$ に近づくことを要請することになる。しかし、shesell は意味的には近いとは言えないので、これでは目的が果たせなそうである。
そもそもこの学習タスクで目指すことは

  • she sells seashells.
  • she purchases seashells.

という2つの文章があった場合、sellpurchase のコンテキストが同じであることからこの2つは近い意味の単語である、と推論することである。
ある単語 $w$ に注目したとき、 別の単語 $w'$ がコンテキストに現れるかどうかは $u_{w}^{\top} v_{w'}$ の大小で測れば良さそうであるというのがこのモデルの気持ちである。

negative sample

目的関数に現れるソフトマックス関数の分母では単語数 $\left| I_W\right| $ に比例する計算量がかかる。$I_W$ は 105 程度のオーダになるので、この計算を毎イテレーションを行うのは現実的でない。 そこでSkip-gramではこの分母の計算量を削減するため、上記の条件付き確率を以下の式で代用している。

$$ P(w' \mid w) = \mathrm{sigmoid}(u_{w}^{\top}v_{w'}) \prod_{l=1}^{N} \mathbb{E}_{w{''} \sim P_{n}}\left[ \mathrm{sigmoid}(u_{w}^{\top}v_{w{''}})\right] $$ ここで、$P_{n}(\cdot)$ は単語集合 $I_W$ 上の確率分布、$N$ は定数(論文中では $N=15$)である。

Word2VecからItem2Vecへ

Skip-gramをアイテムのembeddingの学習に適用する。
各ユーザ $i$ が選択/購入したアイテムのリストを $w_{ij}$ とすれば、

文 ↔ ユーザ
単語 ↔ アイテム

といったアナロジーが使えそうなことに気づく。 ここでItem2Vecの際、$w_{ij}$ が時系列順に並んでる訳ではなければコンテキスト $U_{ij}$ は $w_{ij}$ を除くユーザ $i$ が選択したアイテム全体としてもよいことに留意する。実際元論文ではコンテキストをアイテム全体として学習を行っている。

おわり

Item2Vecの説明といいつつほとんどSkip-gramの説明になってしまった。
Word2Vecの学習にはContinuous Bag of WordsとSkip-gramの2つがあるが、性能は後者の方が良いらしい。 Item2Vecのときは $U = V$ としても良い気がするが、論文中では特に言及されていない。