Counterfactualを知りたい

about Counterfactual ML, Recommender Systems, Learning-to-Rank, Causal Inference, and Bandit Algorithms

因果推論で検索システムを問い直す(2)

はじめに

ランキング学習のシリーズ記事の第二弾です*1.

前回の記事ではUnbiased Learning-to-Rankと呼ばれる, clickというimplicit feedbackを用いて relevanceに対して最適なスコアリング関数を学習するための損失関数を設計する方法について議論しました. その中で紹介したのがexamination parameter  \theta_kの逆数によって損失を重み付けするInverse Propensity Weighting (IPW)と呼ばれる方法でした.

IPWはclickデータのみから真の損失関数を不偏推定することができるという嬉しさがあった反面, 肝心のexamination parameterの推定方法に関しては Result Randomizationの紹介のみに留まっていました.

本記事では, ユーザー体験を著しく害したり, KPIに大きな打撃をもたらすため実務的にはなかなか行うことができないRandomizaitonではなくて offlineのclickデータからexamination parameterを推定する方法を2つ紹介します.

目次

Unbiased Learning-to-Rankのおさらい

Notationの導入の意味も込めて, 簡単にUnbiased LTRの定式化とIPW損失のおさらいをします. 前回記事を読んでいただいている方は流していただいて大丈夫かと思います.

 \mathcal{L}_{U}=\{(q, d, r) \}を各queryとdocumentにrelevance labelが付与されたデータ集合とします. また, あるqueryとデータ集合において共起しているdocumentの集合を \Omega_{q}=\left\{d |(q, d) \in \mathcal{L}_{U} \right\}で表しておきます.

ここで, ランキングアルゴリズムが最小化したい理想的な損失は次の通りでした.


\begin{aligned}
M=\sum_{(q, d, r) \in \mathcal{L}_{U}} r \cdot f\left(q, d, \Omega_{q}\right)
\end{aligned}


ここで,  fはランキングに関連する関数で, 例えばARPやDCGが用いられます. しかし, 通常はrelevance labelが得られないので, 代わりにより安価に手に入るclickを用いて損失を設計することを考えます.

Position-Based Modelに基づくとrelevanceとclickの関係はexaminationという変数を媒介する形で次のように表せます.


\begin{aligned} 
(1) \quad & C =E \cdot R \\ \\
(2) \quad & \underbrace{ \mathbb{P}(C=1 | q, d, k)}_{Click} =\underbrace{ \mathbb{P}(E=1 | k)}_{Examination} \cdot \underbrace{ \mathbb{P}(R=1 | q, d)}_{Relevance} = \theta_{k} \cdot \gamma_{q, d}
\end{aligned}


ここで, 両辺をexamination parameter  \theta_kで割ることによって右辺にrelevance parameter  \gamma_{q, d}が分離されることを利用して IPW損失は次のように定義されます.


\begin{aligned}
\hat{M}_{I P W} &=\sum_{(q, d, k, c) \in \mathcal{L}} \frac{c}{\theta_{k}} \cdot f\left(q, d, \Omega_{q}\right)
=\sum_{(q, d, k, c=1) \in \mathcal{L}} \frac{1}{\theta_{k}} \cdot f\left(q, d, \Omega_{q}\right)
\end{aligned}


IPWに基づいた損失関数 ( \hat{M}_{IPW})は真の損失に対してunbiasedであることが示されるため, 真の損失の代替として正当性を持つという話をしました.

EM Algorithm

 \theta_{k}を推定するための一つ目のアプローチは,  \theta_k \gamma_{q, d}を潜在変数と見てEMアルゴリズムでパラメータを推定するというアプローチです.

Naive EM

まず, パラメータ推定に際して最大化したい対数尤度は次のように表されます.


\begin{aligned}
\log P(\mathcal{L})=\sum_{(q, d, k, c) \in \mathcal{L}} c \log \theta_{k} \gamma_{q, d}+(1-c) \log \left(1-\theta_{k} \gamma_{q, d}\right)
\end{aligned}


 \theta_k \cdot \gamma_{q, d}がclick確率を表すというPBMを思い出せば直感的なはずです.

このまま愚直にE-step, M-stepを書き下すと次のようになります.

E-step: その時点までに得られているパラメータを用いてclickが発生する・しない場合のrelevanceとexaminationの条件付き確率を推定します.


\begin{aligned}
& \mathbb{P}(E=1, R=1 | C=1, q, d, k)=1 \\ 
& \mathbb{P}(E=1, R=0 | C=0, q, d, k)=\frac{\theta_{k}^{(t)}\left(1-\gamma_{q, d}^{(t)}\right)}{1-\theta_{k}^{(t)} \gamma_{q, d}^{(t)}} \\ 
& \mathbb{P}(E=0, R=1 | C=0, q, d, k)=\frac{\left(1-\theta_{k}^{(t)}\right) \gamma_{q, d}^{(t)}}{1-\theta_{k}^{(t)} \gamma_{q, d}^{(t)}}  \\ 
& \mathbb{P}(E=0, R=0 | C=0, q, d, k) = \frac{\left(1-\theta_{k}^{(t)} \right) \left(1-\gamma_{q, d}^{(t)}\right)}{1-\theta_{k}^{(t)} \gamma_{q, d}^{(t)}}
\end{aligned}


M-step: E-stepで推定された条件付き確率を用いて, パラメータを更新します.


\begin{aligned}
& \theta_{k}^{(t+1)}=\frac{\sum_{c, q, d, k^{\prime}} \mathbb{I}_{k^{\prime}=k} \cdot(c+(1-c) P(E=1 | c, q, d, k))}{\sum_{c, q, d, k^{\prime}} \mathbb{I}_{k^{\prime}=k}} \\ \\
& \gamma_{q, d}^{(t+1)}=\frac{\sum_{c, q^{\prime}, d^{\prime}, k} \mathbb{I}_{q^{\prime}=q, d^{\prime}=d} \cdot(c+(1-c) P(R=1 | c, q, d, k))}{\sum_{c, q^{\prime}, d^{\prime}, k} \mathbb{I}_{q^{\prime}=q, d^{\prime}=d}}
\end{aligned}


とEMの更新式を書き下しましたが, 最後の \gamma_{q, d}の更新式で問題が発生します. この更新は, 同じqueryとdocumentのペアについてのデータが一定数存在しないとうまくいきませんが, データ数の少ないペアが大半で更新はあまりうまくいかないと考えられます. また, プライバシーの問題であるqueryについてどのdocumentをclickしたかの情報が手に入らないこともあるらしいです (論文にそう書いてありました.).

Regression-based EM

この \gamma_{q, d}はqueryとdocumentの組み合わせ分だけ存在し, M-stepでの更新がうまくいかない問題に対し, Wang et al. (2018) はRegression-EMと呼ばれるテクニックでこの問題を解消しました.

Regression-EMは, M-stepでのrelevance parameterの推定を一つの学習器で行ってしまおうというアイデアです. 手順の理解には, 擬似コードを見るのが良いと思います.

f:id:usaito:20190525070334p:plain
[Wang et al. (2018)]のAlgorithm 1

擬似コード中の重要な手順の説明をします.

  • (1). 学習器  F(\cdot)を初期化
  • (3). E-step: それまでに得られている \{ \theta_k \}  \{ \gamma_{q, d} \}を用いて条件付きclick確率を推定
  • (4)-(8). E-stepで得られた条件付き確率を用いてrelevance labelをサンプリングし, relevance推定器の学習用のデータを生成
  • (9). 一回前の推定結果とcontextをinputとして, サンプリングしたrelevance labelに対してfit
  • (10)-(11).  \theta_kはNaive EMと同じ式で更新,  \gamma_{q, d}は予測器の出力によって更新

このようにRegression-based EMでは, M-stepにおけるqueryとdocumentの膨大な組み合わせのパラメータ更新式を一つの学習器(論文では, GBDT)で代替することで解決を図っています. 事実, このRegression-based EMを用いることにより前回紹介したResult Randomizationを用いた時よりも良い検索ランキングを提示できるようになったという実験結果が出ました ([Wang et al. (2018)]).

Dual Learning Algorithm

Examination parameter  \theta_kをofflineで推定するもう一つの方法に, Dual Learning Algorithm (DLA)と呼ばれる手法があります. これは, IPWと同じようなアイデアによってExamination parameterの推定のためのunbiasedな損失をclickデータから作れるというアイデアに基づいています.

PBMの仮定のもとで, IPWは次のように両辺を \theta_kで割った形を活用していました.


\begin{aligned}
\frac{\mathbb{P}(C=1 | q, d, k)}{\mathbb{P}(E=1 | k)} = \underbrace{\mathbb{P}(R=1 | q, d)}_{Relevance}
\end{aligned}


[Ai et al. (2018) ]では, これと同じことが両辺を \gamma_{q, d}で割った形にも活用できるよね, というアイデアIRW (Inverse Relevance Weighting)と名付けました.


\begin{aligned}
\frac{\mathbb{P}(C=1 | q, d, k)}{\mathbb{P}(R=1 | q, d)} = \underbrace{\mathbb{P}(E=1 | k)}_{Examination}
\end{aligned}

この操作により, 右辺にexamination parameterを分離できています. 今, examination推定のための真の損失を次のように設定するとします.


\begin{aligned}
M=\sum_{(q, d, e) \in \mathcal{L}_{U}} e \cdot f\left(q, d, \Omega_{q}\right)
\end{aligned}


しかし, ユーザーがdocumentをexamineしたかどうかなんてデータに残っておらず e (Examinationの実現値) がわからないため, この損失は計算不可能です. ここで, 先のIRWを用いると, clickしたかどうかの情報から, Examination parameter推定のための損失に対する不偏推定量を作ることができます.


\begin{aligned}
\hat{M}_{I R W}=\sum_{(q, d, k, c) \in \mathcal{L}_{U}} \frac{c}{\gamma_{q, d}} \cdot f\left(q, d, \Omega_{q}\right)
\end{aligned}


実際, IPWの時と同様にして,


\begin{aligned} 
\mathbb{E}\left[ \hat{M}_{I R W}\right] 
&=\sum_{(q, d, k, c) \in \mathcal{L}} \frac{\mathbb{E} [C ]}{\gamma_{q, d}} \cdot f\left(q, d, \Omega_{q}\right) \\ 
&=\sum_{(q, d, k, c) \in \mathcal{L}} \frac{\theta_{k} \cdot \gamma_{q, d}}{\gamma_{q, d}} \cdot f\left(q, d, \Omega_{q}\right) \\ 
&=\sum_{(q, d, k, c) \in \mathcal{L}} \theta_{k} \cdot f\left(q, d, \Omega_{q}\right) \\ 
&=\sum_{(q, d, k, c) \in \mathcal{L}} \mathbb{P}(E=1 | k) \cdot f\left(q, d, \Omega_{q}\right)
= \mathbb{E}_{E} \left[ M \right] 
\end{aligned}


であることから不偏性が確かめられます.

さて, DLAはこの事実を活用して次のlist-wiseな損失関数をもつ2つの関数 f_{\phi}, \, g_{\psi}を学習します.


\begin{aligned} 
\hat{M}_{I P W}
 & =-\sum_{(q, d, k, c=1) \in \mathcal{L}_{U}} \frac{1}{\theta_{k}} \cdot \log \frac{e^{f_{\phi}\left(x_{q, d}\right)}}{\sum_{d' \in \Omega_{q}} e^{f_{\phi}\left(x_{q, d'}\right)}} \\
\hat{M}_{I R W}
 & =-\sum_{(q, d, k, c=1) \in \mathcal{L}_{U}} \frac{1}{\gamma_{q, d}} \cdot \log \frac{e^{g_{\psi}(k)}}{\sum_{k'} e^{g_{\psi}(k')}}
\end{aligned}


これらの損失を最小化することで,  f_{\phi}, \, g_{\psi}をソフトマックス関数に通すことで, 各パラメータを推定するように学習してほしいという気持ちが込められています.


\begin{aligned} 
\gamma_{q, d}
  \approx \frac{e^{f_{\phi}\left(x_{q, d}\right)}}{\sum_{d' \in \Omega_{q}} e^{f_{\phi}\left(x_{q, d'}\right)}}, \;
\theta_{k}
  \approx \frac{e^{g_{\psi}(k)}}{\sum_{k'} e^{g_{\psi}(k')}}
\end{aligned}


したがって,  f_{\phi}, \, g_{\psi}はお互いの損失の重み付け部分がお互いに依存する形の損失を持っています ( fの損失は \theta_kに依存しており,  gの損失は,  \gamma_{q, d}に依存していますね). これを同時に最適化することにより, end-to-endにExamination parameterの推定と検索ランキング提示のためのスコアリング関数を学習することが可能になります. (ただし, かなり初期値依存が激しそうな印象を持っています)

Dual Learning Algorithmについては以前勉強会(MLPRP)で発表させていただきました. その時のスライドを貼っておきます. 擬似コードや実験結果についてはこちらの資料に掲載してあるので, 参考にしていただけたらと思います.

speakerdeck.com

まとめ

今回は, IPWによって損失を補正するために必要なexamination parameterをofflineで推定する方法を紹介しました. Regression EMなどは推定したいパラメータ数が多い別の問題にも応用できそうな気がしています. 次回は, Position-Based Modelをより現実的なモデル化に拡張した研究について紹介しようと思っています.

参考

[Joachims et al. (2017)]: Thorsten Joachims, Adith Swaminathan, and Tobias Schnabel. 2017. Unbiased learning-to-rank with biased feedback. In Proceedings of the 10th ACM International Conference on Web Search and Data Mining (WSDM ’17).
[Wang et al. (2018)]: Xuanhui Wang, Nadav Golbandi, Michael Bendersky, Donald Metzler, and Marc Najork. 2018. Position Bias Estimation for Unbiased Learning to Rank in Personal Search. In Proceedings of the 11th ACM International Conference on Web Search and Data Mining (WSDM ’18).
[Ai et al. (2018)]: Qingyao Ai, Keping Bi, Cheng Luo, Jiafeng Guo, and W. Bruce Croft. Unbiased learning to rank with unbiased propensity estimation. In The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval (SIGIR’18).

*1:たぶん第三弾まである.

因果推論で検索システムを問い直す(1)

はじめに

ランキング学習のシリーズ記事の第一弾です.

検索システムで達成したいのは, あるqueryに対してよりrelevantなdocumentを上位に提示することです. queryとdocumentのペアに対して, relevanceが大きければ大きいほど, 大きいスコアをつけるような関数を学習することができれば, 効果的にdocumentを並べることができるでしょう. このような関数を学習するためのアルゴリズムがRanking SVMやLambdaMARTなどのランキングアルゴリズムです.

これらのランキングアルゴリズムは, 当然ながら教師ラベルとしてrelevanceを必要とします. 推薦システム的に言えばexplicit feedbackを必要とします. しかし, 学習のたびにannotationするのは高コストなので, clickというimplict feedbackからどうにかスコアリング関数を学習したいというお決まりの動機が湧き上がります. 本記事では, 因果推論で使われる手法を用いてclick情報からrelevance情報を抜き出しつつスコアリング関数を学習する研究を紹介します.

目次

ランキング学習の定式化

さて, ランキング学習を定式化します.  qをquery,  dをdocumentとした時に,  \mathcal{L}_{U}=\{(q, d, r) \}を各queryとdocumentにrelevance labelとして rが付与されたデータ集合とします. このデータ集合は, explicit feedbackとしてrelevanceラベルが得られているのでランキング学習を行うための理想的なデータセットです. また, 便宜上あるqueryとデータ集合において共起しているdocumentの集合を \Omega_{q}=\left\{d |(q, d) \in \mathcal{L}_{U} \right\}で表しておきます.

ここで, Ranking SVMやLambdaMARTなどのランキングアルゴリズムは次の損失関数を最小化することを目指して設計されています.


\begin{aligned}
M=\sum_{(q, d, r) \in \mathcal{L}_{U}} r \cdot f\left(q, d, \Omega_{q}\right)
\end{aligned}


ここで,  fはランキングに関連する関数で, 具体的には次のような損失が用いられることが多いと思います.


\begin{aligned}
A R P &=\sum_{(q, d, r) \in \mathcal{L}_{U}} r \cdot \operatorname{rank}\left(q, d, \Omega_{q}\right) \\ 
D C G &=\sum_{(q, d, r) \in \mathcal{L}_{U}} r \cdot \frac{1}{\operatorname{rank}\left(q, d, \Omega_{q}\right)} 
\end{aligned}


ここで,  \operatorname{rank}\left(q, d, \Omega_{q}\right)は,  d \Omega_qにおける順位を表しています. ARP (Average Relevance Position)は, relevantなdocumentのランキングの総和が小さくなってほしい, DCG (Discounted Cumulative Gain)は, relevantなdocumentのランキングの逆数の総和が小さくなってほしい, という気持ちを反映させています.

これらの損失を具体的にデータから計算できるならば, 最先端のアルゴリズムを用いて効率的にスコアリング関数を学習することができるでしょう. しかし普通そんな理想的な状況*1 はあり得ないので, 次章以降ではより現実的な損失関数の設計方法に関して議論します.

Position-Based Model (PBM)

さて, 現実的には先に定義した \mathcal{L}_{U}のようなexplicit feedbackに基づくデータは得られないはずです. よって, ここからはより現実的なimplicit feedback (click情報) に基づいたデータ集合 \mathcal{L}=\{(q, d, k, c) \}を用いた損失関数の設計方法を考えます. ここで,  k qに対して dを提示した時のpositionを表し,  cはその時にclickが発生したか否かを表します.

さて, ナイーブに考えると既存のアルゴリズムを用いるために, 次のようなrelevanceをclickに入れ替えただけの損失関数を用いることが選択肢の一つとして思いつきます.


\begin{aligned}
\hat{M}_{naive}=\sum_{(q, d, k, c) \in \mathcal{L}} c \cdot f\left(q, d, \Omega_{q}\right)
\end{aligned}


前章で定義した損失関数 Mにおけるrelevanceをclickに入れ替えただけです. しかし, このようなNaive損失をランキングアルゴリズムで最適化してもrelevanceについて最適なランキングを出力するスコアリング関数を得ることはできません. この事実を, 最も単純なclick生成モデルの一つであるPosition-Based Model (PBM)を引き合いに出して説明します.

PBMでは, clickとrelevanceの間に次のような関係性を仮定します.


\begin{aligned} 
(1) \quad & C =E \cdot R \\ \\
(2) \quad & \underbrace{ \mathbb{P}(C=1 | q, d, k)}_{Click} =\underbrace{ \mathbb{P}(E=1 | k)}_{Examination} \cdot \underbrace{ \mathbb{P}(R=1 | q, d)}_{Relevance} = \theta_{k} \cdot \gamma_{q, d}
\end{aligned}


ここで,  Eはexaminationを表す2値確率変数です. このモデル化により, 私たちは次の3つの仮定を置いたことになります.

  1. Queryとdocumentがrelevantかつexaminedのときにのみ, またその時は必ずclickが発生
  2. Relevanceはqueryとdocumentのみに依存し, positionには依存しない
  3. Examinationはpositionのみに依存し, queryとdocumentには依存しない

非常にきつい仮定のように思う一方で, ある程度現実を反映しているかつ単純であるがゆえに思考の整理を助けるモデル化だと思います. このモデル化によって起こる事象を視覚的に説明したのが次の図になります.

f:id:usaito:20190521051747p:plain

この例によると, 1と2番目に提示されたdocumentはexamineされているので, clickはすなわちrelevanceを表しています. 一方で, 9と10番目に提示されたdocumentはもはやexamineされないので, clickはrelevanceを意味しません. つまり, clickされていなくてもrelevantなqueryとdocumentのペアは存在するはずで, 先ほどのナイーブな損失の設計ではそのようなペアの存在を無視してしまうのです.

(この状況ではclickしなかったとしてもそれはirrelevanceを意味しないので, まさにimplicit feedbackといっていい状況でしょう.)

Inverse Propensity Weighting (IPW)

さて, PBMに基づいた議論によりNaive損失を最適化していても, 私たちが本当に得たいスコアリング関数を得ることはできないことがわかりました. しかし, 因果推論の文脈で多用されるInverse Propensity Weighting (IPW)と呼ばれる方法で損失関数を補正することで, 使用可能なデータから真の損失関数をよりよく近似することが可能になります.

IPWに基づいた損失関数は次のように定義されます.


\begin{aligned}
\hat{M}_{I P W} &=\sum_{(q, d, k, c) \in \mathcal{L}} \frac{c}{\theta_{k}} \cdot f\left(q, d, \Omega_{q}\right)
=\sum_{(q, d, k, c=1) \in \mathcal{L}} \frac{1}{\theta_{k}} \cdot f\left(q, d, \Omega_{q}\right)
\end{aligned}


つまり, 各データに対する局所損失をそのデータのexamination parameter ( \theta_k)の逆数で重み付けすることで得られる損失関数のことです. また, 損失の最後の表現を見ると,  c=1についてsummationを取っているためimplict feedbackに基づくデータ集合から計算可能であることがわかります. さて, IPWに基づいた損失関数 ( \hat{M}_{IPW})は次のように真の損失に対してunbiasedであることが示されます.


\begin{aligned} 
\mathbb{E}\left[ \hat{M}_{I P W}\right] 
&=\sum_{(q, d, k, c) \in \mathcal{L}} \frac{\mathbb{E} [C ]}{\theta_{k}} \cdot f\left(q, d, \Omega_{q}\right) \\ 
&=\sum_{(q, d, k, c) \in \mathcal{L}} \frac{\theta_{k} \cdot \gamma_{q, d}}{\theta_{k}} \cdot f\left(q, d, \Omega_{q}\right) \\ 
&=\sum_{(q, d, k, c) \in \mathcal{L}} \gamma_{q, d} \cdot f\left(q, d, \Omega_{q}\right) \\ 
&=\sum_{(q, d, k, c) \in \mathcal{L}} \mathbb{P}(R=1 | q, d) \cdot f\left(q, d, \Omega_{q}\right)
= \mathbb{E}_{R} \left[ M \right] 
\end{aligned}


よって, IPW損失は少なくともbiasの意味でナイーブな損失を用いるよりも真の損失をよく近似するということがわかります.

補足
試しに, Naive損失の期待値を取ってみると,


\begin{aligned} 
\mathbb{E}\left[ \hat{M}_{naive}\right] 
&=\sum_{(q, d, k, c) \in \mathcal{L}} \mathbb{E} [C ] \cdot f\left(q, d, \Omega_{q}\right) \\ 
&=\sum_{(q, d, k, c) \in \mathcal{L}} \theta_{k} \cdot \gamma_{q, d} \cdot f\left(q, d, \Omega_{q}\right) \\ 
&=\sum_{(q, d, k, c) \in \mathcal{L}} \; \underbrace{ \mathbb{P}(E=1 | k)}_{(a)} \cdot \mathbb{P}(R=1 | q, d) \cdot f\left(q, d, \Omega_{q}\right)
\end{aligned}


ここで, (a)の部分は過去のデータのexamination parameterでrelevance最適化を考える上では不要な重みと言えます. そして, この重みは過去のランキング方策に依存します. なぜなら, examinaiton parameterはdocumentの提示positionのみに依存し, positionはデータが集められたときに走っていたランキング方策によって決められているからです. したがって, Naive損失を用いることはすなわち暗黙のうちに過去のランキング方策に似せるような重みを含んだ損失を最適化していることになります.

Examination Parameterの推定方法

前章では, IPW損失を用いることでNaive損失よりも良い損失関数を設計できることを示しました. しかし, IPW損失を計算するためには,  \theta_kが必要です. そして, このパラメータは残念ながら未知です*2. よって, データから推定する必要があります. 本記事では一番率直な方法として, Result Randomizationと呼ばれる方法を紹介します.

ここで,  \mathcal{R}_kを上位N番目までの提示結果をランダムにしたときに得られた k (\leqq N)番目のpositionのqueryとdocumentの集合とします. このとき,


\begin{aligned} 
\mathbb{E} \left[ C \, | \, k \right] 
&=\int_{q, d \in \mathcal{R}_{k}} \mathbb{E} \left[ C | q, d, k \right]  \, P(q, d) \\ 
&=\int_{q, d \in \mathcal{R}_{k}} \theta_{k} \cdot \gamma_{q, d} \, P(q, d) \\ 
&=\theta_{k} \cdot \int_{q, d \in \mathcal{R}_{k}} \gamma_{q, d} \, P(q, d) \\ 
& \propto \theta_{k} \end{aligned}


したがって, randomizationによって得られたデータにおいてk番目のpositionのclick確率はexamination確率に比例します. よって, IPW損失を計算する前段でResult Randomizationによって得られたデータのclick確率からexamination parameterの比がわかるのでそれを重みづけに使用すれば良いということになります.

補足
とは言え, Result Randomizationなんかできないよ, という方がほとんどだと思います. この点に関しては,

  1. 近日中に, Result Randomizationを経ていないデータからexamination parameterを推定する方法について記事を書く予定です.
  2. 何もやらないよりも, 例えば過去のデータにおいてclickが発生したポジションの逆数でNaive損失を補正するだけでも結果が変わってくると思います. 過去のランキング方策が提示したポジションを保存しましょう!

実験結果

補正がうまくいった結果の例を一つ掲載しておきます. 詳しくは[Joachims et al. (2017)]をお読みください. WSDM2017のBest Paper Awardを受賞した論文です. この論文では, PBMに従う人工データを生成し, Naive損失とIPW損失をそれぞれRanking SVMで最適化したときに得られたスコアリング関数の性能をテストデータにおけるARPで評価しました. 結果は次の通りです.

f:id:usaito:20190521050153p:plain
[Joachims et al. (2017)]のFigure 2

横軸が大きな値をとるほど, Positionごとのexamination確率の変化が激しくなっています. 赤線がIPWによる補正を用いた時のARPですので, 特にexamination確率の変化が激しいときにNaive損失を用いたbaselineを大きく上回っていることがわかります. 実験的にもIPWによる補正は効果的であることがわかる例となっています.

まとめ

本記事では, 因果推論的なアプローチでclickしか得られていないimplicit feedbackデータからrelevanceに対して最適化されたスコアリング関数を得るための損失を設計する方法について書きました. 記事中にも書きましたが, 近日中に補正に必要なexaminaiton parameterをデータから推定する方法についての記事を書こうと思っています. また, 記事中に何か変なことを書いていたらご指摘いただけると幸いです.

参考

[Joachims et al. (2017)]: Thorsten Joachims, Adith Swaminathan, and Tobias Schnabel. 2017. Unbiased learning-to-rank with biased feedback. In Proceedings of the 10th ACM International Conference on Web Search and Data Mining (WSDM ’17).
[Wang et al. (2018)]: Xuanhui Wang, Nadav Golbandi, Michael Bendersky, Donald Metzler, and Marc Najork. 2018. Position Bias Estimation for Unbiased Learning to Rank in Personal Search. In Proceedings of the 11th ACM International Conference on Web Search and Data Mining (WSDM ’18).
[Agarwal et al. (2019)]: Aman Agarwal, Xuanhui Wang, Cheng Li, Mike Bendersky, and Marc Najork. 2019. Addressing Trust Bias for Unbiased Learning-to-Rank. In Proceedings of the 2019 World Wide Web Conference (WWW ’19)

*1:全てのqueryとdocumentのペアについてrelevanceが得られている状況

*2:あるポジションに提示したdocumentは何%の確率でexamineされるかを知ることはできません, よね.

統計的学習理論(Rademacher Complexityを用いた期待損失の導出)

はじめに

前回の記事では, 仮説集合 \mathcal{H}が有限である場合の, 仮説 h \in \mathcal{H}の予測損失の上界をHoeffding's ineqを用いて導きました. しかし, 無限仮説集合に対しては同様の方法で実用的な上界を得ることは不可能でした. したがって, 今回は無限仮説集合に対応する方法の一つであるRademacher Complexityを用いて予測損失の上界を導いてみようと思います.

目次

定式化のおさらい

統計的学習理論のモチベーションをおさらいします. 前回記事をお読みいただいている方は読み飛ばしていただいても結構です.

入力空間 \mathcal{X}に値をとる確率変数を X, 入力空間から出力空間 \mathcal{Y}への写像 f: \mathcal{X} \rightarrow \mathcal{Y}としてlabeling functionと呼びます. また, 入力が従う確率分布を Dとします ( X \sim D). 損失関数として \ell: \mathcal{Y} \times \mathcal{Y} \rightarrow \mathbb{R}_+を用いるとします. このとき, ある仮説 h: \mathcal{X} \rightarrow \mathcal{Y}の予測損失は次のように定義されます.


\begin{aligned}
R (h) = \mathbb{E}_{ X \sim D } \left[  \ell \left( h(X), f(X) \right) \right]
\end{aligned}


つまり, 分布 Dにおける予測値 h(X)の損失の期待値です.  R(h)をできるだけ小さくするような仮説 hを見つけ出すことが機械学習の目標です.

ここで, データの真の分布 Dは未知であることが問題であり, 有限の観測データから, 予測損失をできるだけ小さくする仮説を導く必要がありました. 予測損失を観測データから評価する上で最も大きな手がかりは経験損失です. データ  \left\{  X_i  \right\}_{i=1}^n が観測されたときの経験損失は次のように定義されます.


\begin{aligned}
\hat{R} (h) = \mathbb{E}_{ X \sim \hat{D} } \left[  \ell \left( h(X), f(X) \right) \right] = \frac{1}{n} \sum_{i=1}^n \ell  \left( h(X_i), f(X_i) \right)
\end{aligned}


ここで,  \hat{D}はデータの経験分布です.

また, 仮説集合 \mathcal{H}に含まれる仮説の中で予測損失と経験損失を最小化するような仮説をそれぞれ  h^* \in \arg \min_{h \in \mathcal{H}} R (h),  \hat{h} \in \arg \min_{h \in \mathcal{H}} \hat{R} (h)と表していました.

ここで私たちが導きたいのは, 観測データから計算できる経験損失を最小化する基準で得られる仮説 \hat{h}の期待損失 R ( \hat{h}) と 仮説集合 \mathcal{H}を用いたときに達成され得る最小の期待損失 R(h ^*)の差を次のように評価することです.


少なくとも 1 - \deltaの確率で次の不等式が成り立つ.


\begin{aligned}
R (\hat{h} ) \leqq   R (h^* ) + complexity \left( \mathcal{H} \right) + confidence ( \delta )
\end{aligned}



以降の目標は, 仮説集合 \mathcal{H}のcomplexity項と \deltaに依存するconfidence項を具体的に求めることです.

Rademacher Complexity

Rademacher Complexityの導入

何度か名前が出てきていますが, 実数値関数の集合の複雑さの指標として解釈されるEmpirical Rademacher ComplexityRademacher Complexityを定義します.

Empirical Rademacher Complexity: 空間 \mathcal{Z}上の実数値関数からなる集合を \mathcal{G} \subset \{ g: \mathcal{Z} \rightarrow \mathbb{R} \}とする. また, 入力点の集合を \mathcal{S} = \{ z_1, ..., z_n \}としする. さらに, 同数の+1と-1を等確率でとる独立な確率変数 (Rademacher variables)を \boldsymbol {\sigma} = \{ \sigma_1, ..., \sigma_n \}とする. このとき, 集合 \mathcal{G}のEmpirical Rademacher Complexityは次のように定義される.


\begin{aligned}
\hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{G}) =  \mathbb{E}_{ \boldsymbol {\sigma} } \left[ \sup_{g \in \mathcal{G}} \frac{1}{n} \sum_{i=1}^n \sigma_i g(z_i)  \right]
\end{aligned}


さらに, それぞれのサンプル z_iはある分布 Dに独立に従う( \mathcal{S} \sim D^{n}) とき, 次のRademacher Complexityを定義します.

Rademacher Complexity: 入力点 z_iがある分布 Dに従う確率変数であるとする. また, ある入力点のサンプル集合 \mathcal{S} = \{ z_1, ..., z_n \}が与えられたとする. このとき, 集合 \mathcal{G}のRademacher Complexityは次のように定義される.


\begin{aligned}
\mathfrak{R}_{n } (\mathcal{G}) 
& =  \mathbb{E}_{ S \sim D^n } \mathbb{E}_{ \boldsymbol {\sigma} } \left[ \sup_{g \in \mathcal{G}} \frac{1}{n} \sum_{i=1}^n \sigma_i g(z_i)  \right] \\
& =  \mathbb{E}_{ S \sim D^n }  \left [ \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{G})   \right]
\end{aligned}


Rademacher variables  \boldsymbol{ \sigma } \mathcal{S}とは独立にサンプルされます. よって, 入力点とそれに対するランダムなラベル付け (z_1, \sigma_1), ..., (z_n, \sigma_n)を 最大でどれくらい予測できてしまう関数を含むか?を測っていると解釈できます. 入力集合 \mathcal{S}とはなんら無関係なランダムノイズの組みに対してよくfittingできてしまうほど,  \mathcal{G}は複雑である, ということです.

Rademacher Complexityの推定

さて, 定義よりRademacher ComplexityはEmpirical Rademacher Complaxityをデータ集合について期待値をとったものでした. しかし, ある入力点のサンプル集合 \mathcal{S}が与えられたとき, Empirical Rademacher Complexityは必ずしもRademacher Complexityに一致するとは限りません. よって, Rademacher Complexityの推定誤差の裾確率を評価することで, 誤差の大きさを見積もってみようと思います.

Rademacher Complexity vs Empirical Rademacher Complexity: 入力空間 \mathcal{Z}上の実数値関数の集合を \mathcal{G} \subset \{ g: \mathcal{Z} \rightarrow [0, M ] \} とする. このとき, 少なくとも 1 - \delta / 2 ( \delta > 0)の確率で次の不等式が成り立つ.


\begin{aligned}
\mathfrak{R}_{n } (\mathcal{G})  \leqq  \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{G}) + M \sqrt {  \frac{ \log \frac{2}{ \delta } }{ 2n }  }
\end{aligned}


導出
 g(z_1, ..., z_n) = \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{G}) =  \mathbb{E}_{ \boldsymbol {\sigma} } \left[ \sup_{g \in \mathcal{G}} \frac{1}{n} \sum_{i=1}^n \sigma_i g(z_i)  \right]とおきます. このとき,


\begin{aligned}
& | g(z_1, ..., z_i, .., z_n) - g(z_1, ..., z'_i, .., z_n) | \\
& =\mathbb{E}_{ \boldsymbol {\sigma} } \left[ \sup_{g \in \mathcal{G}} \left( \frac{1}{n} \sum_{j=1}^n \sigma_j g(z_j) \right) 
-  \sup_{g \in \mathcal{G}} \left( \frac{1}{n} \sum_{j=1}^n \sigma_j g(z_j) - \frac{1}{n} \sigma_i g(z_i) + \frac{1}{n} \sigma_{i} g(z_i') \right) \right]  \\
& \leqq  \mathbb{E}_{ \boldsymbol {\sigma} } \left[ \sup_{g \in \mathcal{G}} \left( \frac{1}{n} \sum_{j=1}^n \sigma_j g(z_j) \right) 
-  \sup_{g \in \mathcal{G}} \left( \frac{1}{n} \sum_{j=1}^n \sigma_j g(z_j) \right)  + \sup_{g \in \mathcal{G}} \frac{1}{n} |  \sigma_i g(z_i) - \sigma_{i} g(z_i') |   \right]  \\
& =  \mathbb{E}_{ \boldsymbol {\sigma} } \left[ \sup_{g \in \mathcal{G}} \frac{1}{n} |  \sigma_i g(z_i) - \sigma_{i} g(z_i') |   \right]
\leqq \frac{M}{n}
\end{aligned}


ここでは, 関数 gの値域が [0, M]であることとと,  \sigma \in \{-1, +1\}であることを用いました. これにて, 関数 gこちらの記事で導出したMacDiamid's ineqのboundedness conditonを満たすことがわかります. したがって, MacDiamid's ineqで c_i = M / nとおくと, 任意の \epsilon > 0に対して, 次の不等式が成り立つことがわかります.


\begin{aligned}
\mathbb{P} \left(  \mathbb{E} \left [ g ( \mathcal{S} ) \right] -   g ( \mathcal{S} ) \geqq \epsilon  \right) \Leftrightarrow
\mathbb{P} \left(  \mathfrak{R}_{n } (\mathcal{G}) -   \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{G}) \geqq \epsilon  \right) \leqq \exp \left( - \frac{ 2 n \epsilon^2 }{M^2} \right)
\end{aligned}


ここで, 左辺を \delta / 2とおくと,


\begin{aligned}
\exp \left( - \frac{ 2 n \epsilon^2 }{M^2} \right) = \frac{\delta}{2} \:  \Leftrightarrow  \: M \sqrt {  \frac{ \log \frac{2}{ \delta } }{ 2n }  }
\end{aligned}


です. よって,


\begin{aligned}
\mathbb{P} \left(  \mathfrak{R}_{n } (\mathcal{G}) -   \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{G}) \leqq  M \sqrt {  \frac{ \log \frac{2}{ \delta } }{ 2n }  } \right) 
>  1 - \frac{\delta}{2}
\end{aligned}


ですので, 所望の不等式を得ます.

Uniform law of large numbers

さて, 前章でRademacher Complexityを導入しました. 本章では, これを用いて次のUniform law of large numbersを導きます.

Uniform law of large numbers: 入力空間 \mathcal{Z}上の実数値関数の集合を \mathcal{G} \subset \{ g: \mathcal{Z} \rightarrow [0, M ] \} とする. また, ある入力点のサンプル集合を \mathcal{S} = \{ z_1, ..., z_n \}とし, 一つ一つのサンプルは独立に zと同じ分布 Dに従うとする. このとき, 少なくとも 1 - \delta ( \delta > 0)の確率で次の不等式が成り立つ.


\begin{aligned}
\mathbb{E} [g (z) ]  \leqq \frac{1}{n} \sum_{i=1}^n g (z_i) + 2  \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{G}) + 3 M \sqrt {  \frac{ \log \frac{2}{ \delta } }{ 2n }  }
\end{aligned}



導出
まず,


\begin{aligned}
\varphi( z_1, ..., z_n)  = \sup_{g \in \mathcal{G}}  \left( \mathbb{E} [g (z) ]  - \frac{1}{n} \sum_{i=1}^n g (z_i)  \right)
\end{aligned}


と置きます. ここで先ほど同様に \varphi( z_1, ..., z_n)  に対してMacDiamid's ineqを適用します.


\begin{aligned}
\varphi( z_1, ..., z_i, ..., z_n)  - \varphi( z_1, ..., z'_i, ..., z_n)  \leqq   \sup_{g \in \mathcal{G}} \frac{g(z_i) - g(z_i') }{n}  \leqq \frac{M}{n}
\end{aligned}


より,  \varphi( z_1, ..., z_n)はboundedness conditionを満たすので, 少なくとも 1 - \delta / 2の確率で, 次の不等式が成り立ちます.


\begin{aligned}
\varphi( \mathcal{S} )  \leqq \mathbb{E}_{S \sim D^n} \left[  \varphi( \mathcal{S} ) \right] + M \sqrt {  \frac{ \log \frac{2}{ \delta } }{ 2n }  }
\end{aligned}


ここで,  \varphi( z_1, ..., z_n) = \varphi( \mathcal{S} ) としました.

次に, 今導いた不等式の右辺の第1項をRademacher Complexityを用いて評価します.


\begin{aligned}
 \mathbb{E}_{S} \left[  \varphi( \mathcal{S} ) \right]  
& =  \mathbb{E}_{S} \left[ \sup_{g \in \mathcal{G}}  \left( \mathbb{E} [g (z) ]  - \frac{1}{n} \sum_{i=1}^n g (z_i)  \right)  \right]   \\
& =  \mathbb{E}_{S} \left[ \sup_{g \in \mathcal{G}}  \left( \mathbb{E}_{S'} \left [ \frac{1}{n} \sum_{i=1}^n g (z'_i)  \right ]  - \frac{1}{n} \sum_{i=1}^n g (z_i)  \right)  \right] \\
& =  \mathbb{E}_{S} \left[ \sup_{g \in \mathcal{G}}  \left( \mathbb{E}_{S'} \left[ \frac{1}{n} \sum_{i=1}^n g (z'_i)   - \frac{1}{n} \sum_{i=1}^n g (z_i)  \right ] \right)   \right] \\
& \leqq \mathbb{E}_{S, S'} \left[ \sup_{g \in \mathcal{G}}  \frac{1}{n} \sum_{i=1}^n (g (z'_i)   -  g (z_i) )    \right] 
\end{aligned}


ここで,  z_i z'_iは同一分布 Dに従います. また,  \sigma_iは+1と-1を等確率でとるRademacher variableです. このとき,  (g (z'_i)   -  g (z_i) )   \sigma_i (g (z'_i)   -  g (z_i) )  は同一分布に従います. よって,


\begin{aligned}
& \mathbb{E}_{S, S'} \left[ \sup_{g \in \mathcal{G}}  \frac{1}{n} \sum_{i=1}^n (g (z'_i)   -  g (z_i) )    \right]  \\
&= \mathbb{E}_{ \boldsymbol{\sigma}, S, S'} \left[ \sup_{g \in \mathcal{G}}  \frac{1}{n} \sum_{i=1}^n \sigma_i (g (z'_i)   -  g (z_i) )    \right] \\
& \leqq \mathbb{E}_{ \boldsymbol{\sigma}, S, S'} \left[ \sup_{g \in \mathcal{G}}  \frac{1}{n} \sum_{i=1}^n \sigma_i g (z'_i)    \right] 
+ \mathbb{E}_{ \boldsymbol{\sigma}, S, S'} \left[ \sup_{g \in \mathcal{G}}  \frac{1}{n} \sum_{i=1}^n - \sigma_i g (z_i)    \right] \\
& = \mathbb{E}_{ \boldsymbol{\sigma}, S} \left[ \sup_{g \in \mathcal{G}}  \frac{1}{n} \sum_{i=1}^n \sigma_i g (z_i)    \right] 
+ \mathbb{E}_{ \boldsymbol{\sigma}, S} \left[ \sup_{g \in \mathcal{G}}  \frac{1}{n} \sum_{i=1}^n  \sigma_i g (z_i)    \right] \\
& = 2 \mathfrak{R}_{n } (\mathcal{G})
\end{aligned}


これらの結果を合わせると, 少なくとも 1 - \delta / 2の確率で次の不等式が成り立ちます.


\begin{aligned}
\sup_{g \in \mathcal{G}}  \left( \mathbb{E} [g (z) ]  - \frac{1}{n} \sum_{i=1}^n g (z_i)  \right) =  \varphi( z_1, ..., z_n)  \leqq 2 \mathfrak{R}_{n } (\mathcal{G}) + M \sqrt {  \frac{ \log \frac{2}{ \delta } }{ 2n }  }
\end{aligned}


さらに, 前章の (Rademacher Complexity vs Empirical Rademacher Complexity) の結果を用いると, 少なくとも 1 - \delta / 2の確率で次の不等式が成り立ちます.


\begin{aligned}
 2 \mathfrak{R}_{n } (\mathcal{G}) \leqq 2 \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{G})  + 2M \sqrt {  \frac{ \log \frac{2}{ \delta } }{ 2n }  }
\end{aligned}


以上の結果を統合し, union boundを用いることで, 少なくとも 1 - \deltaの確率で次の不等式が成り立つことを得ます.


\begin{aligned}
\mathbb{E} [g (z) ]  \leqq \frac{1}{n} \sum_{i=1}^n g (z_i) + 2  \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{G}) + 3 M \sqrt {  \frac{ \log \frac{2}{ \delta } }{ 2n }  }
\end{aligned}


 L_p損失を用いた時の予測損失の上界

さてようやく本題です. ここでは, 次にように表される L_p損失を用いたときの, 予測損失の確率的上界を導出します.  p \geqq 1で, また損失は定数 M^{p}でboundされるとします.


\begin{aligned}
R (h) = \mathbb{E}_{ X \sim D } \left[  \ell_{p} \left( h(X), f(X) \right) \right] = \mathbb{E}_{ X \sim D } \left[ \left| h(X) -  f(X) \right|^p \right] 
\end{aligned}


ここで,  L_p損失と仮説の合成関数の集合 \mathcal{H}_p = \left\{ x \rightarrow | h(x) - f(x) |^p : h \in \mathcal{H} \right\}のEmpirical Rademacher Complexityを仮説集合 \mathcal{H}のEmpirical Rademacher Complexityで評価します.

まず, 先ほどの集合 \mathcal{H}_pを関数 \phi_p: x \rightarrow |x|^pと集合 \mathcal{H}' = \{x \rightarrow h(x) - f(x): h \in \mathcal{H} \}を用いて,  \mathcal{H}_p = \{ \phi_p \circ  \mathcal{H}' \}と表しておきます. 関数 \phi_pは,  pM^{p-1}-Lipschitzなので, Talagrand's lemma*1を用いると,


\begin{aligned}
 \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{H}_p) \leqq  pM^{p-1}  \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{H}')
\end{aligned}


が成り立ちます. また,


\begin{aligned}
 \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{H}') 
& = \mathbb{E}_{ \boldsymbol {\sigma} } \left[ \sup_{h \in \mathcal{H}} \frac{1}{n} \sum_{i=1}^n \sigma_i (h(x_i) - f(x_i))  \right] \\
& \leqq \mathbb{E}_{ \boldsymbol {\sigma} } \left[ \sup_{h \in \mathcal{H}} \frac{1}{n} \sum_{i=1}^n \sigma_i h(x_i)   \right] 
+ \underbrace{ \mathbb{E}_{ \boldsymbol {\sigma} } \left[ \sup_{h \in \mathcal{H}} \frac{1}{n} \sum_{i=1}^n \sigma_i f(x_i)   \right]}_{= 0} \\
& =  \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{H}) 
\end{aligned}


したがって, 結局のところ


\begin{aligned}
 \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{H}_p) \leqq  pM^{p-1}  \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{H})
\end{aligned}


なので,  L_p損失と仮説集合の合成関数の集合のEmpirical Rademacher Complexityは, 仮説集合のEmpirical Rademacher Complexityで上から評価できます. この事実と, Uniform law of large numbersにおいて,  \mathcal{G} \mathcal{H}_pとすれば, 任意の仮説 h \in \mathcal{H}に対して, 少なくとも 1 - \deltaの確率で次の不等式が成り立ちます.


\begin{aligned}
R (h) -  \hat{R} (h) = & \mathbb{E}_{X \sim D} [ \ell_{p} (h(X),  f(X ))  ] - \frac{1}{n} \sum_{i=1}^n \ell_{p} (h(X_i),  f(X_i)) \\
&  \leqq   2pM^{p-1}  \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{G}) + 3 M^p \sqrt {  \frac{ \log \frac{2}{ \delta } }{ 2n }  }
\end{aligned}


この不等式を用いれば, 少なくとも 1 - \deltaの確率で次の不等式が成り立つことを導くことができます.


\begin{aligned}
R (\hat{h} ) -   R (h^* ) 
& =  R (\hat{h} ) -   \hat{R} (\hat{h} ) + \underbrace{ \hat{R} (\hat{h} )  - \hat{R} (h^* ) }_{ \leqq 0 } + \hat{R} (h^*) -  R (h^* )  \\
& \leqq R (\hat{h} ) -   \hat{R} (\hat{h} ) +  \hat{R} (h^*) -  R (h^* )  \\
& \leqq 2 \sup_{h \in \mathcal{H}} \left| \hat{R} (h) -  R (h) \right| \\
& \leqq 2 \left(  2pM^{p-1}  \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{H}) + 3 M^p \sqrt {  \frac{ \log \frac{2}{ \delta } }{ 2n }  }  \right)
\end{aligned}


したがって, 所望の確率的上界を次のように得ることができました.


\begin{aligned}
R (\hat{h} )  \leqq  R (h^* )  + \underbrace{ 4 pM^{p-1}  \hat{\mathfrak{R}}_{ \mathcal{S} } (\mathcal{H})}_{complexity}
+   \underbrace{ 6 M^p \sqrt {  \frac{ \log \frac{2}{ \delta } }{ 2n }  } }_{confidence}
\end{aligned}


さいごに

本記事では,  L_p損失を用いた場合の予測損失の確率的上界をRademacher Complexityを用いて導出してみました. 相変わらず, 私の誤解で誤った記述をしている可能性が大いにありますので, 見つけた場合はご指摘いただけたら幸いです.

参考

[金森 (2015)] 金森敬文. 2015. 統計的学習理論. 講談社 機械学習プロフェッショナルシリーズ.
[Mohri et al. (2012)] Mohri, M.; Rostamizadeh, A.; and Talwalkar, A. 2012. Foundations of Machine Learning. MIT Press.

*1:Mohri et al. (2012)のlemma 4.2

統計的学習理論(有限仮説集合の場合の予測損失の上界)

はじめに

最近自身の研究で使うため統計的学習理論の勉強をしています. 2回に渡って基本的な内容をまとめてみます.

目次

定式化

まず, 統計的学習理論のモチベーションを述べます. 入力空間 \mathcal{X}に値をとる確率変数を X, 入力空間から出力空間 \mathcal{Y}への写像 f: \mathcal{X} \rightarrow \mathcal{Y}としてlabeling functionと呼ぶことにします. また, 入力が従う確率分布を Dとします ( X \sim D). 損失関数として \ell: \mathcal{Y} \times \mathcal{Y} \rightarrow \mathbb{R}_+を用いるとします. このとき, ある仮説 h: \mathcal{X} \rightarrow \mathcal{Y}の予測損失は次のように定義されます.


\begin{aligned}
R (h) = \mathbb{E}_{ X \sim D } \left[  \ell \left( h(X), f(X) \right) \right]
\end{aligned}


つまり, 分布 Dにおける予測値 h(X)の損失の期待値です.  R(h)をできるだけ小さくするような仮説 hを見つけ出すことが目標です.

ここで, 予測損失 R(h)がわかっていれば, 問題は随分簡単になるのですが, 残念ながらデータの真の分布 Dは未知です. よって, 有限の観測データから, 予測損失をできるだけ小さくする仮説を導く必要があります. 予測損失を観測データから評価する上で最も大きな手がかりは経験損失です. データ  \left\{ X_i \right\}_{i=1}^n が観測されたときの経験損失は次のように定義されます.


\begin{aligned}
\hat{R} (h) = \mathbb{E}_{ X \sim \hat{D} } \left[  \ell \left( h(X), f(X) \right) \right] = \frac{1}{n} \sum_{i=1}^n \ell  \left( h(X_i), f(X_i) \right)
\end{aligned}


ここで,  \hat{D}をデータの経験分布としました. これは, データ数が nのときに, 確率 1 / nで各観測データ X_iの値をとるような分布のことです.

 \mathcal{H} = \{h_1, h_2, ..., h_{} \}をある有限な仮説集合とします. このとき, この仮説集合に含まれる仮説の中で予測損失と経験損失を最小化するような仮説をそれぞれ h^*, \hat{h}と表しておきます.


\begin{aligned}
h^* \in \arg \min_{h \in \mathcal{H}} R (h), \quad \hat{h} & \in \arg \min_{h \in \mathcal{H}} \hat{R} (h)
\end{aligned}


 \hat{h}は, Empirical Risk Minimizerと呼んだりもします.  h^*, \hat{h}について, 定義より次の2つの不等式が成り立ちます.


\begin{aligned}
R (h^* )  \leqq R( \hat{h} ), \quad \hat{R} (\hat{h})  \leqq  \hat{R} (h^*)
\end{aligned}


ここで私たちが導きたいのは, 観測データから計算できる経験損失を最小化する基準で得られる仮説 \hat{h}の期待損失 R ( \hat{h}) と 仮説集合 \mathcal{H}を用いたときに達成され得る最小の期待損失 R(h ^*)の差を次のように評価することです.


少なくとも 1 - \deltaの確率で次の不等式が成り立つ.


\begin{aligned}
R (\hat{h} ) \leqq   R (h^* ) + [extra \: term]
\end{aligned}



次章では, 仮説集合 \mathcal{H} \deltaに依存する [extra term]を具体的に求めます.

有限仮説集合の場合の予測損失の上界の導出

さて, 本記事で考えている仮説集合 \mathcal{H}は有限でした. またここでは , 損失関数 \ellがある定数 Mで上からboundできるとします (例えば, 01-lossの場合は,  M=1). この場合, こちらの記事で紹介したHoeffding's ineqを用いて R (\hat{h} )  R (h^* )の差を評価することができます.

まず,


\begin{aligned}
R (\hat{h} ) -   R (h^* ) 
& =  R (\hat{h} ) -   \hat{R} (\hat{h} ) + \underbrace{ \hat{R} (\hat{h} )  - \hat{R} (h^* ) }_{ \leqq 0 } + \hat{R} (h^*) -  R (h^* )  \\
& \leqq R (\hat{h} ) -   \hat{R} (\hat{h} ) +  \hat{R} (h^*) -  R (h^* )  \\
& \leqq 2 \max_{h \in \mathcal{H}} \left| \hat{R} (h) -  R (h) \right|
\end{aligned}


ここで, 右辺の裾確率をunion boundとHoeffding's ineqを用いて次のように評価します. 途中で係数2が登場しているのは, 誤差の絶対値を評価するため両側の裾確率を考慮しているからです.


\begin{aligned}
& \mathbb{P} \left( 2 \max_{h \in \mathcal{H}} \left| \hat{R} (h) -  R (h) \right|  \geqq \epsilon  \right) \\
& \leqq \sum_{h \in \mathcal{H}} \mathbb{P} \left( \left| \hat{R} (h) -  R (h) \right| \geqq \epsilon / 2 \right) \\
&  \leqq  \sum_{h \in \mathcal{H}} 2 \exp \left( - \frac{n \epsilon^2} {2M^2} \right) \quad \because Hoeffding's \,  ineq \\
& = 2 | \mathcal{H} | \exp \left( - \frac{n \epsilon^2} {2M^2} \right)
\end{aligned}


右辺を \deltaと置いて \epsilonについて解けば,


\begin{aligned}
2 | \mathcal{H} | \exp \left( - \frac{n \epsilon^2} {2M^2} \right) = \delta \Leftrightarrow \epsilon = M \sqrt {  \frac{2}{n} \log \frac{2 | \mathcal{H}| }{ \delta}  }
\end{aligned}


よって,


\begin{aligned}
\mathbb{P} \left( 2 \max_{h \in \mathcal{H}} \left| \hat{R} (h) -  R (h) \right|  \geqq M \sqrt {  \frac{2}{n} \log \frac{2 | \mathcal{H}| }{ \delta}  }  \right)  \leqq  \delta
\end{aligned}


なので, 少なくとも 1 - \deltaの確率で次の不等式が成り立ちます.


\begin{aligned}
2 \max_{h \in \mathcal{H}} \left| \hat{R} (h) -  R (h) \right|  \leqq M \sqrt {  \frac{2}{n} \log \frac{2 | \mathcal{H}| }{ \delta}  }
\end{aligned}


以上の結果を統合すると, 最終的な目標であった次の不等式が少なくとも 1 - \deltaの確率で成り立つという結果を得ます.


\begin{aligned}
R (\hat{h} ) \leqq   R (h^* )  + M \sqrt {  \frac{2}{n} \log \frac{2 | \mathcal{H}| }{ \delta}  }
\end{aligned}


したがって, 仮説集合 \mathcal{H}が有限である場合,  R ( \hat{h} )  R (h^* )に 少なくとも \mathcal{O}_p \left( \sqrt{ \log | \mathcal{H} |  / n } \right) で収束することがわかります.

さいごに

今回は, 仮説集合が有限な場合の \hat{h}の予測損失の上界をHoeffding's ineqを用いて導出しました. しかし, 今回導いた上界は仮説集合が無限である場合 (i.e.,  |\mathcal{H}| = \infty) 実用的ではありません. (上界の第2項を見れば一目瞭然.) よって, 次はRademacher Complexityという仮説集合の複雑さの指標を用いて仮説集合が無限の場合も意味を成す期待損失の上界を導いてみます.

参考

[金森 (2015)] 金森敬文. 2015. 統計的学習理論. 講談社 機械学習プロフェッショナルシリーズ.
[Mohri et al. (2012)] Mohri, M.; Rostamizadeh, A.; and Talwalkar, A. 2012. Foundations of Machine Learning. MIT Press.

多腕バンディットアルゴリズムのリグレット解析

はじめに

活用と探索をいい具合にバランスしつつ, 報酬を最大化することを目標とする多腕バンディット問題は, web広告の最適化などへの応用が期待されることから大きな注目を集めています. この分野に関してはこれまでにいくつかの素晴らしい記事が存在しますが, その理論解析に踏み込んだものはありませんでした. 本記事では, 有名なアルゴリズムを例にして, バンディットアルゴリズムの理論解析の雰囲気をさらってみます.

目次

定式化

まずは(確率的)多腕バンディット問題を定式化します. 定式化に関しては, こちらの記事で非常にわかりやすくまとめられているため, ここでは簡単に記述します.

多腕バンディット問題では, playerに対してarmと呼ばれるいくつかの行動の選択肢が与えられます. 各armは報酬と呼ばれる未知の確率分布に従う確率変数と対応付けられています. playerは, 各時刻にいずれかのarmを選択してそのarmに紐づく報酬の実現値を観測します. このような試行を繰り返すうちに, 各armの報酬が従う確率分布のパラメータを推定(探索)しながら, 同時に得られる報酬を最大化すること(活用)がplayerの目標です.

今, armの集合を \mathcal{A} = \{1, ..., a , ..., K \}とし, ある時刻 tにあるarm  aから観測される報酬は, 区間 [0, 1]に値をとる確率変数 X_{a, t}に独立に従うとします. また, arm  aを選択することで得られる期待報酬を \mathbb{E} \left[ X_{a, t} \right] = \mu_{a} と表記します.

このとき, 多腕バンディット問題では次のように定義される期待累積リグレット(以降, リグレット)を考えます.


\begin{aligned}
R (T) 
& = \mathbb{E} \left[ \max_{a \in \mathcal{A}} \sum_{t = 1}^T  X_{a, t} - X_{a(t), t}  \right] \\
& =  \sum_{a = 1}^K \Delta_{a} \mathbb{E} \left[ T_a \right]
\end{aligned}


ここで, 総試行回数を T, 時刻 tにおいて選択されるarmを a(t), 最適なarmとの期待報酬のgapを \Delta_a , 全ての試行が終了した時点においてarm  aを選択した回数を T_aとしました. リグレットは, 最適な固定戦略  a^{*} = \arg \max_{ a \in \mathcal{A}} \mu_{a} を取ることで得られる期待累積報酬と playerの取る戦略  \{ a(t) : t = 1, 2, ..., T \} によって得られる期待累積報酬の差の合計だと解釈できます. これは, 達成し得た最大の報酬からの後悔の量と言えるのでリグレットと名付けられています.

playerはリグレット R(T)を最小化することを目指します. また, armの選択を出力するアルゴリズム (player) が達成するリグレットの上界を導出することが主要な理論的興味になります. ちなみにLai et al. (1985)によると, 期待リグレットの理論限界は \Omega (\log T)であることが知られており, アルゴリズムにはこの理論限界を達成していて欲しいということになります.

Upper Confidence Boundアルゴリズム

多腕バンディット問題の中でとてもよく知られたアルゴリズムに, Upper Confidence Bound (UCB) アルゴリズムというものがあります. これは, ある時刻においてそれまでに観測された報酬によって推定される期待報酬の上側信頼区間が最大になるarmを選択するアルゴリズムになります. 上側信頼区間を推定する際の有意水準は時刻 tに対して t^{- \alpha}程度に設定されることが多いので今回はそれに従います.

期待報酬の上側信頼区間が最大になるarmを選択するという方針の解釈は決して難しくないと思いますが, 具体的なアルゴリズムはどのように構成すれば良いのでしょうか? 実は, 以前の記事で紹介したHoeffding's ineqを活用することでこの方針を実現するアルゴリズムを導くことができます.

ここで, 時刻 tの試行が終了した時点においてarm  aを選択した回数を T_{a, t}として, 時刻 tにおける選択前までの観測報酬を用いて各armの期待報酬の推定量を次のように構築します. と言っても, ただ観測報酬を平均しているだけです.


\begin{aligned}
\hat{\mu}_{ a, T_{a, t - 1} } = \frac{1}{ T_{a, t - 1} } \sum_{s=1}^{t-1} \mathbb{I} \{ a(t) = a \} X_{a, s}
\end{aligned}


この推定量が真の期待報酬 \mu_aを大きく過小評価してしまう確率はHoeffding's ineqを用いて次のように評価できます.


\begin{aligned}
\mathbb{P} \left( \mu_a - \hat{\mu}_a \geqq \epsilon \right) \leqq  \exp \left( - 2 T_{a, t - 1} \epsilon^2 \right)
\end{aligned}


右辺を t^{-\alpha}とおいて変形を加えると,


\begin{aligned}
\mathbb{P} \left( \mu_a < \hat{\mu}_a +  \sqrt{ \frac{\alpha \log (t)}{2 T_{a, t - 1}}  } \right) > 1 - t^{- \alpha}
\end{aligned}


よって, Hoeffding's ineqを用いることにより各armの有意水準 t^{- \alpha}での上側信頼区間の上界を得ることができました. これに従い, UCBアルゴリズムは次のように定義されるUCBスコアが最大になるようなarmを各時刻において選択するようなアルゴリズムになっています.


\begin{aligned}
UCB_{a, t} = \hat{\mu}_a +  \sqrt{ \frac{\alpha \log (t - 1)}{2 T_{a, t - 1}} }
\end{aligned}


擬似コードは次の通りです.

f:id:usaito:20190501053314p:plain

リグレット上界の導出

さて, 本章では前章で紹介したUCBアルゴリズムが達成するリグレットの上界を導きます. しかしここでは, UCBアルゴリズムが理論限界を達成することを示すため簡易なリグレット上界の導出を採用します.

以降一般性を失うことなく,  a^* = 1であるとします. armの期待報酬の推定が大きく外れる場合とそうでない場合に場合分けする方針をとります. まず, 次の2つの事象 A_t B_tを考えます.

  •  A_t:  \hat{\mu}_{ a, T_{a, t} } \leqq \mu_a +  \sqrt{ \frac{\alpha \log (t)}{2 T_{a, t}} }
  •  B_t:  \hat{\mu}_{ 1, T_{1, t} } \geqq \mu_1 -  \sqrt{ \frac{\alpha \log (t)}{2 T_{1, t}} }

つまり,  A_tは時刻 tにおけるあるarmの期待報酬の推定がそこまで大きくならないという事象のこと.  B_tは時刻 tにおける最適なarmの期待報酬の推定がそこまで小さくならないという事象のことです. まず,  A_tの補事象が起こる確率を評価します.


\begin{aligned}
\mathbb{P} \left( A_t^c \right) 
& = \mathbb{P} \left(  \hat{\mu}_{ a, T_{a, t} } -  \mu_a  >   \sqrt{ \frac{\alpha \log (t)}{2 T_{a, t}} } \right) \\
& \leqq t^{- \alpha} \quad \because Hoeffding's \: ineq 
\end{aligned}


同様に,  \mathbb{P} \left( B_t^{c} \right) \leqq  t^{- \alpha}  です. よってunion boundを用いれば,


\begin{aligned}
\sum_{t = 1}^T  \mathbb{E} \left[ \mathbb{I} \{ A_t^c \cup B_t^c \} \right]  
& \leqq \sum_{t = 1}^T \left(  \mathbb{E} \left[ \mathbb{I} \{ A_t^c \} \right] + \mathbb{E} \left[ \mathbb{I} \{ B_t^c \} \right]   \right) \\
& \leqq 2 \sum_{t = 1}^T t^{- \alpha} \\
& \leqq  2 \left( 1 + \int_1^{\infty} x^{-\alpha} dx \right)  \\
& = \frac{2 \alpha}{\alpha - 1}
\end{aligned}


です.

次に, 事象 A_t, B_tが共に成り立つ場合を考えます. まず, 最適ではないarmが選択されるのは, 次の不等式が成り立つ場合です.


\begin{aligned}
\hat{\mu}_a +  \sqrt{ \frac{\alpha \log (t)}{2 T_{a, t}} } > \hat{\mu}_1 +  \sqrt{ \frac{\alpha \log (t)}{2 T_{1, t}} }
\end{aligned}


今, 事象 A_tが成り立っているので, 両辺に \sqrt{ \frac{\alpha \log (t)}{2 T_{a, t}} }を足せば,


\begin{aligned}
\hat{\mu}_{ a, T_{a, t} } + \sqrt{ \frac{\alpha \log (t)}{2 T_{a, t}} } \leqq \mu_a +  \sqrt{ \frac{2 \alpha \log (t)}{ T_{a, t}} }
\end{aligned}


です. これを用いれば,


\begin{aligned}
\mu_1 
&  \leqq \hat{\mu}_{ 1, T_{a, t} } +  \sqrt{ \frac{\alpha \log (t)}{2 T_{1, t}} } \quad \because B_t \\
& \leqq  \hat{\mu}_a +  \sqrt{ \frac{\alpha \log (t)}{2 T_{a, t}} }  \\
& \leqq \mu_a +  \sqrt{ \frac{2 \alpha \log (t)}{ T_{a, t}} }
\end{aligned}


です.  T_{a, t}について整理すると,


\begin{aligned}
T_{a, t}  \leqq 2 \Delta_a^{-2} \log (t) \leqq  2 \Delta_a^{-2} \log (T) 
\end{aligned}


を得ます. 最終的には,


\begin{aligned}
\mathbb{E} \left[ T_a \right] 
& \leqq  2  \left (   \Delta_a^{-2} \log (T)  +  \frac{\alpha}{\alpha - 1}\right)
\end{aligned}


よって,


\begin{aligned}
R(T)
& = \sum_{a \neq 1} \Delta_{a} \mathbb{E} \left[ T_a \right] \\
&=  \sum_{a \neq 1} \Delta_{a} \left(  \mathbb{E} \left[ \mathbb{I} \{ A_t \cap B_t \} \right]   +  \mathbb{E} \left[ \mathbb{I} \{ A_t^c \cup B_t^c \} \right]  \right) \\
& \leqq  2 \sum_{a \neq 1}   \left (   \frac{ \log (T) }{\Delta_a } +  \frac{\alpha}{\alpha - 1} \Delta_a \right)
\end{aligned}


というリグレットの上界が導かれます. これによりUCBアルゴリズムのリグレット上界は悪くても \mathcal{O} (\log T)のオーダーです. Lai et al. (1985)における理論限界は \Omega(\log T)でしたので, UCBアルゴリズムは理論限界と同様のオーダーを達成することがわかりました.

さいごに

今回は, 単純にarmの報酬のみから次に選択するarmを出力するようなアルゴリズムを構成する問題を扱いました. しかし, 実際のweb広告などではuserやarmのcontextに報酬の確率分布のパラメータが依存すると考えるのが自然です. このようなcontextに依存するような問題設定を文脈付きバンディット問題などと呼びます. 文脈付きバンディット問題に関しては, 以前Qiitaに紹介記事を書きましたが, 理論解析には触れられていませんでした. 近いうちに本ブログで, 文脈付きバンディット問題の理論にも触れたいと思っています.

また毎度のことですが, 私の誤解で誤った記述をしてしまっている可能性が大いにあります. 誤りを見つけた場合はご指摘いただけると幸いです.

参考

[Jamison (2017)] Kevin Jamieson. Stochastic Multi-armed bandits, regret minimization. URL: https://courses.cs.washington.edu/courses/cse599i/18wi/resources/lecture3/lecture3.pdf
[Lai et al. (1985)] T.L Lai and Herbert Robbins. Asymptotically efficient adaptive allocation rules. Adv. Appl. Math., 6(1):4–22, March 1985.
[本多, 中村 (2016)] 本多淳也, 中村篤祥. バンディット問題の理論とアルゴリズム. 講談社 機械学習プロフェッショナルシリーズ, 2016.

有用な確率不等式のまとめ

はじめに

機械学習に関連する諸分野では何かしらの統計量(期待判別誤差やリグレットなど)を上から評価したい場面が多くあります. そのような場面で大活躍するのが確率不等式と呼ばれる不等式の数々です. 今後本ブログでもこれらの不等式を多用することが予想されるため, 一度まとめておきます. いくつかの不等式は証明もします. 証明は, MLPシリーズの『統計的学習理論』のAppendix Aを参考に, 自分なりに行間を埋めてみました.

目次

Jensen's inequality

まず, 凸関数の定義を確認します.

凸関数: 関数 f: \mathbb{R}^d \rightarrow (- \infty, \infty ] が, 任意の {\bf x}, {\bf y} \in \mathbb{R}^d と任意の \alpha \in [0, 1]に対して,


\begin{aligned}
f (\alpha {\bf x} + (1 - \alpha) {\bf y}) \leqq \alpha f ({\bf x}) + (1 - \alpha) f({\bf y})
\end{aligned}

を満たす時,  fを凸関数という.


ここで, 以降のHoeffding's ineqの証明などでも用いられる期待値演算にJensen's ineqを適用した形を紹介します.

Jensen's Inequality:  Xを有限な期待値を持つ確率変数とする.  f(x) Xの値域を含む区間において凸な関数とする. この時, 次の不等式が成り立つ.


\begin{aligned}
f \left( \mathbb{E} \left[ X \right] \right) \leqq \mathbb{E} \left[ f(X)  \right]
\end{aligned}

Markov's inequality / Chebyshev's inequality

ここでは, Markov's ineqとChebyshev's ineqを紹介します. これらは統計学の教科書にもだいたい掲載されていると思います.

Marcov's Inequality: 非負確率変数 Xと任意の \epsilon > 0, p > 0に対して, 次の不等式が成り立つ.

\begin{aligned}
\mathbb{P} \left( X \geqq \epsilon \right) \leqq \frac{ \mathbb{E} \left[  X^p \right] }{\epsilon^p}
\end{aligned}


導出


\begin{aligned}
\epsilon^p \mathbb{P} \left( X \geqq \epsilon \right) 
& = \epsilon^p \int_{ X \geqq \epsilon } f(x)dx\\
& =  \int_{ X \geqq \epsilon } \epsilon^p f(x)dx\\
& \leqq  \int_{ X \geqq \epsilon } X^p f(x)dx \\
& \leqq \mathbb{E} \left[  X^p \right]
\end{aligned}



特に,  X | Z - \mathbb{E} [ Z ] |,  p=2と置いた時の形をChebyshev's ineqと呼びます. この不等式は, 確率変数 Zがその期待値から大きく外れた値をとる確率を分散を用いて評価しています.

Chebyshev's inequality: 確率変数 Zと任意の \epsilon > 0に対して, 次の不等式が成り立つ.

 
\begin{aligned}
\mathbb{P} \left( |Z - \mathbb{E} [ Z ] | \geqq \epsilon \right) \leqq \frac{ \mathbb{V} \left[  Z \right] }{\epsilon^2}
\end{aligned}

Hoeffding's inequality

ここで紹介するHoeffding's ineqは機械学習の論文で頻出なので, かなり重要です. これを示す前に, 一つ補題 (Hoeffding's lemma) を示します.

Hoeffding's lemma: 確率変数 X \mathbb{E} [X ] = 0,  a \leqq X \leqq bを満たすとき, 任意の s > 0に対して, 次の不等式が成り立つ.


\begin{aligned}
\mathbb{E} \left[ \exp (sX) \right] \leqq \mathbb{E} \left[ \exp \left( \frac{s^2 (b - a)^2}{8} \right) \right]
\end{aligned}


導出
 a \leqq x \leqq bを満たすような xについて,  0 \leqq \lambda \leqq 1を次のように置きます.


\begin{aligned}
\lambda = \frac{b - x}{b - a}
\end{aligned}


これを変形すると,  sx = s \lambda a + (1 - \lambda)bです.  \exp(\cdot)は凸関数なので,


\begin{aligned}
e^{sx} = e^{s \lambda a + (1 - \lambda)b} \leqq \lambda e^{sa} + (1 - \lambda) e^{sb} = \frac{b - x}{b - a} e^{sa} + \frac{x - a}{b - a} e^{sb}
\end{aligned}


両辺の期待値をとると,


\begin{aligned}
\mathbb{E} \left[ e^{sX} \right] = \mathbb{E} \left[ \frac{b - X}{b - a} e^{sa} + \frac{X - a}{b - a} e^{sb} \right] = \frac{b }{b - a} e^{sa} - \frac{a}{b - a} e^{sb}
\end{aligned}


 Xの期待値が0であることを用いました. ここで,  p = - a / (b - a)と置くと,  1 - p = b / (b-a)なので,


\begin{aligned}
\frac{b }{b - a} e^{sa} - \frac{a}{b - a} e^{sb} = (1 - p) e^{sa} + p e^{sb}
\end{aligned}

と表されます. また,  u = s(b - a)と置き,  uについての関数 \phi (u)を次のように定義しておきます.


\begin{aligned}
 & (1 - p) e^{sa} + p e^{sb}  = pe^{(1 - p)u} + (1 - p)e^{-pu} \\
& \phi (u)  = \log \left( pe^{(1 - p)u} + (1 - p)e^{-pu} \right)
\end{aligned}


このとき,


\begin{aligned}
\phi(0)  & = \log \left(p - (1 - p) \right) =  0 \\
\phi'(0) & = 0 \quad \left( \phi'(u) = -p + \frac{p e^u}{1 - p + p e^u}  \right) \\
\phi'' (u) & = \frac{1 - p}{1 - p + pe^u} \cdot \frac{pe^u}{1 - p + pe^u}  \\
& = \frac{1 - p}{1 - p + pe^u} \cdot \left( 1 - \frac{1 - p}{1 - p + pe^u} \right) \leqq \frac{1}{4}
\end{aligned}



したがって, テイラーの定理より,  0 \leqq v \leqq uとなる vが存在して,


\begin{aligned}
\phi(u) = \phi(0) + u \phi'(0) + \frac{u^2}{2} \phi'' (0) \leqq  \frac{u^2}{8}
\end{aligned}


 u = s(b - a)だったことを思い出すと,


\begin{aligned}
\phi(u) = \log \left( \mathbb{E} \left[ e^{sX} \right] \right) \leqq \frac{  s^2(b - a)^2}{8}
\Leftrightarrow  \mathbb{E} \left[  \exp \left( sX \right)  \right] \leqq  \exp \left( \frac{s^2 (a-b)^2}{8} \right)
\end{aligned}


これにて, Hoeffding's lemmaを得ました.



このHoeffding's lemmaを用いて, Hoeffding's ineqを導出します.

Hoeffding's Inequality: 確率変数 X_1, ..., X_nは独立でかつそれぞれが有界区間 [ a_i, b_i ]に値をとるとする. この時, 任意の \epsilon > 0に対して次の不等式が成り立つ.


\begin{aligned}
\mathbb{P} \left( \frac{1}{n} \sum_{i=1}^n X_i - \frac{1}{n} \sum_{i=1}^n \mathbb{E} \left[  X_i \right] \geqq \epsilon \right) \leqq \exp \left( - \frac{2 n^2 \epsilon^2}{ \sum_{i=1}^n (b_i - a_i) } \right) 
\end{aligned}


導出
まず,  Z_i = X_i - \mathbb{E} [ X_i ]と置く. このとき,  \mathbb{E} [Z_i ] = 0であり,  a_i \leqq Z_i \leqq b_iが成り立つので, 確率変数 Z_iはHoeffding's lemmaの仮定を満たす. ここで,


\begin{aligned}
& \mathbb{P} \left( \sum_{i=1}^n Z_i \geqq \epsilon  \right)  \\
& = \mathbb{P} \left(  \exp \left(  s \sum_{i=1}^n Z_i  \right) \geqq \exp(s \epsilon ) \right) \\
& \leqq   \exp \left( -s \epsilon \right)  \mathbb{E} \left[ \prod_{i=1}^n \exp (s Z_i) \right]  \quad \because Marcov's \: ineq \\
& =  \exp \left( -s \epsilon \right) \prod_{i=1}^n  \mathbb{E} \left[ \exp (s Z_i) \right] \quad \because independency  \\
& \leqq \exp \left( -s \epsilon \right) \prod_{i=1}^n  \exp \left( \frac{s^2 (b_i - a_i)^2} {8} \right) \quad \because Hoeffding's \: lemma \\
& = \exp \underbrace{  \left( \frac{s^2}{8} \sum_{i=1}^n (b_i - a_i)^2 - s \epsilon  \right)}_{(1)}
\end{aligned}

よりtightなboundを得るため, (1)を最小化する s =   4  \epsilon / \sum_{i=1}^{n} (b_i - a_i)^{2} を代入すると,


\begin{aligned}
\mathbb{P} \left(  \sum_{i=1}^n Z_i \geqq \epsilon  \right) \leqq  \exp \left( \frac{-2 \epsilon ^2} { \sum_{i=1}^n (b_i - a_i)^2 } \right)
\end{aligned}


ここで,  Z_i = X_i - \mathbb{E} [ X_i ] を代入し, 両辺を nで割ることで, Hoeffding's ineqを得ます.



Hoeffding's ineqは, 有限仮説集合の汎化誤差解析やバンディット問題におけるリグレット解析など至る所で出てくる印象です. 定理自体,  X_1, ..., X_nが独立であれば成り立ちますが, 実際はiidが仮定されていることが多いと思います. (iidの仮定がある場合,  \frac{1}{n} \sum_{i=1}^n \mathbb{E} \left[  X_i \right] の部分が単に \mathbb{E} \left[  X_i \right]となります.)

McDiarmid's inequality

最後に, McDiarmid ineqを紹介します. これは, 後に紹介する予定のRademacher Complexityに関連する不等式を導くときなどに役立ちます. これを示す前に, 一つ補題 (Azuma's ineq) を示します.

Azuma's inequality: 確率変数 X_i, Z_i, V_i: i=1, ..., nに対して,  V_i X_i, ..., X_iの関数でありかつ \mathbb{E} [V_i | X_i, ..., X_{i-1} ] =0 が成り立つとする. また,  Z_iは,  X_1, ..., X_{i-1}の関数として表すことができ,  i = 1, ..., nについて Z_i \leqq V_i \leqq Z_i + c_iを満たすような定数 c_iが存在するとします. このとき, 任意の \epsilon > 0に対して, 次の不等式が成り立つ.


\begin{aligned}
\mathbb{P} \left(  \sum_{i=1}^n V_i \geqq \epsilon  \right) \leqq \exp \left(  - \frac{2 \epsilon^2}{\sum_{i=1}^m c_i^2 } \right) 
\end{aligned}


導出
まず V_iについての部分和を S_k = \sum_{i=1}^k V_iとしておきます. ここで t > 0を用いて,


\begin{aligned}
& \mathbb{P}\left(  \sum_{i=1}^n V_i \geqq \epsilon  \right) \\ 
& =  \mathbb{P} \left(  S_n \geqq \epsilon  \right)  \\
& = \mathbb{P} \left( \exp \left( t S_n \right) \geqq \exp (t \epsilon ) \right) \\
& \leqq \exp (- t \epsilon ) \mathbb{E} \left[ \exp \left( t S_n \right) \right] \quad \because Marcov's \, ineq  \\
& = \exp (-t \epsilon ) \mathbb{E}_{ X_1, ..., X_{n - 1} } \left[ \exp (t S_n)  \mathbb{E}_{X_n} \left[  \exp(t V_n) | X_1, ..., X_{n-1} \right]  \right] \\
& \leqq  \exp (-t \epsilon ) \mathbb{E}_{ X_1, ..., X_{n - 1} } \left[ \exp (t S_n) \right ] \exp \left( \frac{t^2 c_n^2}{8} \right)  \quad  \because  Hoeffding's \, lemma \\
&  \leqq  \exp  \underbrace{ \left( \frac{t^2}{8} \sum_{i=1}^n c_i^2  - t \epsilon \right) }_{(2)} 
\end{aligned}


Hoeffding's lemmaの部分は, 仮定より X_1, ..., X_{n-1}で条件付けたとき V_iの期待値が0であることからlemmaの仮定を満たします. 最後の不等式は, Hoeffding's lemmaを繰り返し用いることで得ます. 最後に, (2)を最小化する t = 4 \epsilon / \sum_{i=1}^{n} c_i^{2}を代入することで,


\begin{aligned}
\mathbb{P} \left(  \sum_{i=1}^n V_i \geqq \epsilon  \right) \leqq \exp \left(  - \frac{2 \epsilon^2}{\sum_{i=1}^m c_i^2 } \right) 
\end{aligned}


を得ます.



さて, これを用いてMcDiarmid's ineqを示します.

McDiarmid's inequality: ある集合 \mathcal{X}に値をとる独立な確率変数を X_1, ..., X_nとする. また, 関数 f: \mathcal{X}^n \rightarrow \mathbb{R}と任意の x_1, ..., x_n, x_i' \in \mathcal{X}について, 次の条件 (boundedness condition) を満たすような定数 c_1, ..., c_nが存在するとする.


\begin{aligned}
| f(x_1, ..., x_i, ..., x_n) - f(x_1, ..., x'_i, ..., x_n)  | \leqq c_i, \quad \forall i \in \{1, 2, .., n \}
\end{aligned}

この時, 任意の \epsilon > 0に対して, 次の不等式が成り立つ.


\begin{aligned}
\mathbb{P} \left(  f(X_1, ..., X_n) -  \mathbb{E} \left[ f(X_1, ..., X_n) \right] \geqq \epsilon \right)   \leqq  \exp \left(  - \frac{2 \epsilon^2}{\sum_{i=1}^n c_i^2} \right)
\end{aligned}


導出
最初に f (X_1, ..., X_n) f (S)と表しておきます. また,  V_1, ..., V_n


\begin{aligned}
V_k = \mathbb{E} \left[ f(S) | X_1, ..., X_k \right] - \mathbb{E} \left[ f(S) | X_1, ..., X_{k-1} \right]
\end{aligned}


とします ( V_1 = \mathbb{E} \left[ f(S) | X_1 \right] - \mathbb{E} \left[ f(S) \right]). ここで, 今定義した V_kAzuma's ineqの仮定を満たすことを確認します. まず,  V_kの第1項は X_1, ..., X_kで条件付けた時の期待値なので, これらの確率変数の関数です. さらに,


\begin{aligned}
\mathbb{E} [ V_k  |  X_1, ..., X_{k-1} ] 
& = \mathbb{E}_{X_k, ..., X_n} \left[ f(S) | X_1, ..., X_k \right] - \mathbb{E} \left[ f(S) | X_1, ..., X_{k-1} \right] \\
& = \mathbb{E} \left[ f(S) | X_1, ..., X_{k-1} \right] - \mathbb{E} \left[ f(S) | X_1, ..., X_{k-1} \right] \\
& = 0
\end{aligned}


ここで, boundedness conditionより


\begin{aligned}
\sup_{x, x'} \{\mathbb{E} \left[ f(S) | X_1, ..., X_{k-1}, x \right] - \mathbb{E} \left[ f(S) | X_1, ..., X_{k-1}, x' \right]   \}  \leqq c_i
\end{aligned}


を満たします. ここで,  Z_k = \inf_x  \mathbb{E} \left[ f(S) | X_1, ..., X_{k-1}, x \right] とすると,  X_k, V_k, Z_k: k=1, ..., nAzuma's ineqの仮定を満たします. 最後に,


\begin{aligned}
\sum_{i=1}^n V_i  
&  =   \left( \mathbb{E} \left[ f(S) | X_1, ..., X_{k-1} \right] - \mathbb{E} \left[ f(S) | X_1, ..., X_{k-1} \right] \right) \\
& \quad + \cdot \cdot \cdot +  \left(  \mathbb{E} \left[ f(S) | X_1 \right] - \mathbb{E} \left[ f(S) \right] \right) \\
& = f(X_1, ..., X_n) -  \mathbb{E} \left[ f(X_1, ..., X_n) \right]
\end{aligned}


と見て, これにAzuma's ineqを適用すれば, MacDiamid's ineqを得ます.



さいごに

さて, 本記事では機械学習の周辺分野でよく使われる確率不等式(特に, Hoeffding's inequalityとMacDiamid's inequality)を導出してみました. 一見どのように役立つかわからないかもしれませんが, 今後本ブログで記事にするつもり内容と密接な関連があります.
なお, 導出部分で私が勘違いして誤った記述をしている可能性が大いにあります. 誤りを見つけた場合, ご指摘いただけたら幸いです.

参考

[金森 (2015)] 金森敬文. 2015. 統計的学習理論. 講談社 機械学習プロフェッショナルシリーズ.
[Duchi] John Duchi. Probability Bounds. URL: http://www.cs.berkeley.edu/~jduchi/projects/probability_bounds.pdf.

Causal Embeddingsの解説と追試

はじめに

前回は, ログデータの観測確率が一様ではない場合に傾向スコアで補正した損失関数を用いるPropensity Matrix Factorizationを紹介しました. しかし, 2018年のRecsysにてBest Paper Awardを受賞したCausal Embeddings for Recommendation [Bonner+ 2018]で, Propensity MFを実験的に上回る手法が提案されました. 本記事では人工データを用いて, その提案手法の追試を行います.

目次

定式化のおさらい

推薦アルゴリズムの学習の定式化をおさらいします. しかし, ここでは[Bonner+ 2018]のモデル化に合わせるため, 前回記事とは少し異なる部分があることに注意してください.

各ユーザーを u \in \{1, ..., U \} = \mathcal{U}, 各アイテムを i \in \{ 1, ...., I \} = \mathcal{I}とします. また,  P_{u,i}をユーザー uがアイテム iに対して有する真のPreferenceとします.  P_{u,i}はユーザーとアイテムの組み合わせにのみ依存して決定的な値を取るとします. 推薦アルゴリズムの学習によって達成したいのは, 次のように定義される真の損失関数 \mathcal{L}を最小化するような予測値集合 \hat{P} = \{ \hat{P}_{u, i} \}_{(u,i) \in \mathcal{U} \times \mathcal{I}} を得ることです.


\begin{aligned}
\mathcal{L} \left( \hat{P} \right) =  \frac{1}{U \cdot I} \sum_{u,i} \delta \left( P_{u,i}, \hat{P}_{u,i} \right)
\end{aligned}


ここで,  \delta (\cdot, \cdot): \mathbb{R} \times \mathbb{R} \rightarrow \mathbb{R}は真のRatingとそれに対する予測値を入力とする関数でした.

さらに, 自分たちが観測できるPreferenceと観測できないそれを区別するために, Recommendationの頭文字をとって新たな確率変数 Rを導入します. この Rは2値確率変数であり各ユーザーとアイテムのペアに対応して独立に存在します.  R_{u,i} = 1 ならば, 過去にユーザー uにアイテム iが推薦されたことを意味し, そのペアについてのPreferenceが観測されることを意味します. [Bonner+ 2018]では, 観測されるPreferenceを Y_{u,i}^{obs} とした時, 次のようなモデルを想定しています.


\begin{aligned}
Y_{u,i}^{obs} = P_{u,i} \cdot R_{u,i}
\end{aligned}


このモデル化では, 過去に推薦が発生したペア( R_{u,i} = 1)については,  Y^{obs}_{u,i} = P_{u,i} となり真のPreferenceが観測されますが, それ以外のペアについては Y_{u,i}^{obs} = 0となりPreferenceが観測されないということになります.(ユーザーとアイテムのMatrixを考えた時に, Preferenceが観測されず欠損となってしまっている部分を0で表しているというイメージです.)

以降, 推薦有無を表す確率変数 R_{u,i}は議論の中で重要な役割を担うので, その期待値にnotationを用意しておきましょう.


\begin{aligned}
\mathbb{E} \left [ R_{u,i} \right ] = \pi_c (u,i)
\end{aligned}


 \pi_c (u,i)はログデータが集められた際に走っていた推薦方策を表し, logging policyと呼びます. これを用いると私たちが学習の際に用いることができる(真の損失関数の推定に用いることができる)学習データセット  \mathcal{S}を次のように表すことができます.


\begin{aligned}
\mathcal{S} = \left\{ \left( u, i, Y_{u,i}^{obs} \right):  R_{u,i} = 1  \right\}
\end{aligned}

Radomized Policy下でのNaive Loss

前章で定義した真の損失関数 \mathcal{L}を知ることができるならばそれ最適化するパラメータを得れば良いのですが, それは不可能です. よって手元にあるログデータ \mathcal{S}から損失関数を推定して, それを信じた上で, 最適化を実行する必要があります.

まずNaiveな損失関数 \hat{\mathcal{L}}_{naive} \left( \hat{P} \right)は次のように定義されます. これは, 過去に推薦が発生したデータについての損失を単純に平均するという非常に単純な推定量です.


\begin{aligned}
\hat{\mathcal{L}}_{naive} \left( \hat{P} \right) =  \frac{1}{ | \mathcal{S} |  } \sum_{(u,i) \in \mathcal{S}} \delta \left( Y^{obs}_{u,i}, \hat{P}_{u,i} \right)
\end{aligned}


[Bonner+ 2018]は, このNaive推定量  \hat{\mathcal{L}}_{naive} がある条件下で望ましい性質を満たすことに注目します. まず, ある条件とは次のようなものです.


\begin{aligned}
 \pi_c (u,i) = \pi_{rand} = \frac{C \cdot | \mathcal{S} |}{U \cdot I}  \quad (constant)
\end{aligned}


つまり, logging policyがユーザーとアイテムに非依存で, 一様な定数であるという要求です. このようなlogging policyを特にrandomized policyと呼ぶことにし,  \pi_{rand}で表します. 定数部分は実はなんでも良いのですが, ここでは後の式変形が綺麗な形になる定数を置いています. さて, このlogging policyがrandomized policyであるという仮定が成り立っているときに, Naive推定量  \hat{\mathcal{L}}_{naive}の期待値をとってみましょう.


\begin{aligned}
 \mathbb{E}_R  \left[  \hat{\mathcal{L}}_{naive} \left( \hat{P} \right) \right ]  
& = \mathbb{E}_R  \left[ \frac{1}{ | \mathcal{S} |  } \sum_{(u,i) \in \mathcal{S}} \delta \left( Y^{obs}_{u,i}, \hat{P}_{u,i} \right) \right] \\
& = \mathbb{E}_R  \left[ \frac{1}{ | \mathcal{S} |  } \sum_{(u,i) \in \mathcal{U} \times \mathcal{I} } R_{u,i} \cdot  \delta \left( P_{u,i}, \hat{P}_{u,i} \right) \right] \\
& = \frac{1}{ | \mathcal{S} |  } \sum_{(u,i) \in \mathcal{U} \times \mathcal{I}} \mathbb{E}_{R_{u,i}}  \left[ R_{u,i}  \right] \cdot  \delta \left( P_{u,i}, \hat{P}_{u,i} \right)\\
& =  \frac{1}{ | \mathcal{S} |  } \sum_{(u,i) \in \mathcal{U} \times \mathcal{I}} \pi_{rand} \cdot  \delta \left( P_{u,i}, \hat{P}_{u,i} \right) \\
& =  \frac{1}{ | \mathcal{S} |  } \sum_{(u,i) \in \mathcal{U} \times \mathcal{I}} \left(  \frac{C \cdot | \mathcal{S} |}{U \cdot I}  \right)  \cdot  \delta \left( P_{u,i}, \hat{P}_{u,i} \right) \\
& = C \left(  \frac{1}{U \cdot I} \sum_{u,i} \delta \left( P_{u,i}, \hat{P}_{u,i} \right)  \right) \\
& = C \cdot \mathcal{L} \left( \hat{P} \right) \\
& \propto \mathcal{L}  \left( \hat{P} \right)
\end{aligned}


これにより, randomized policyで集められたログデータを用いて計算されるNaiveな損失関数は真の損失関数に比例します. よって,  \hat{\mathcal{L}}_{naive}を最適化するような学習は妥当であると言えそうです.

Causal Embeddings (CausE)

[Bonner+ 2018]の提案手法であるCausEは, 前章で紹介したNaive推定量の性質を活用します. CausEは学習データが次のように表せることを要求します.


\begin{aligned}
 \mathcal{S} = \mathcal{S}_t  \cup \mathcal{S}_c
\end{aligned}


ここで,  \mathcal{S}_t はrandomized policyによって集められた学習データで,  \mathcal{S}_c u iに依存するlogging policyによって集められた学習データを表します. ただし,  \mathcal{S}_t \mathcal{S}_cよりも要素数の少ない集合で良いです. また, それぞれ異なるembeddingsによって構成される2つの予測値集合 \hat{P}_t = \{ \theta_{u, t}^{\top} \beta_{i, t}  \}_{ (u,i) \in \mathcal{U} \times \mathcal{I} } ,  \hat{P}_c = \{ \theta_{u, c}^{\top} \beta_{i, c}  \}_{ (u,i) \in \mathcal{U} \times \mathcal{I} } を用いた次の損失関数をCausEは最適化します.


\begin{aligned}
 \hat{\mathcal{L}}_{CausE}  
 & = \underbrace{ \frac{1}{| \mathcal{S}_t |}  \sum_{ (u, i) \in S_t } \delta \left( P_{u,i}, \theta_{u, t}^{\top} \beta_{i, t}   \right) }_{(1)}
 + \underbrace{ \frac{1}{| \mathcal{S}_c |}  \sum_{ (u, i) \in S_c } \delta \left( P_{u,i},  \theta_{u, c}^{\top} \beta_{i, c}   \right) }_{(2)} \\
& + \underbrace{ \lambda  \left( \sum_{u \in \mathcal{U}}  \| \theta_{u, t} - \theta_{u, c}  \|_p + \sum_{i \in \mathcal{I}}  \| \beta_{i, t} - \beta_{i, c}  \|_p \right)}_{(1)}
\end{aligned}


(1)はrandomized policyによって構成された \mathcal{S}_tに対するNaiveな損失(Biasなし), (2)はlogging policyによって構成された \mathcal{S}_cに対するNaiveな損失(Biasあり), (3)は(2)によって得られるembeddings ( \theta_{u, c}, \beta_{i, c} ) が望ましい損失である(1)から得られるembeddings ( \theta_{u, t}, \beta_{i, t} ) と乖離することに対してかける正則化で,  \lambdaはその強さを調節するハイパーパラメータです. ( \theta_{u, c}, \beta_{i, c} を更新するときのみ, (3)の正則化を考慮に入れるアルゴリズムになっています.)

確かにこの損失関数からは, Biasがあるが考慮できるデータが多い損失(2)から得られたembeddingsが, Biasのない損失(1)から得られたembeddingsと離れないようにしようというお気持ちが読み取れます. しかし, この正則化について理論的な考察が論文でなされているわけではなく, どこかもどかしい気もします.

さて, 最終的な予測値集合は次のように作ります.(これが実験的に最も良かったらしいです.)


\begin{aligned}
 \hat{P}_{CausE} = \{  \theta_{u, c}^{\top} \beta_{i, c}  \}_{ (u,i) \in \mathcal{U} \times \mathcal{I} }
\end{aligned}

簡易実験

CausEとPropensity Matrix Factorizationの性能を同様の人工データを用いて比較します.

実験設定

設定に関しては前回の実験と同じです.

しかし, 前回記事ではあえて書いていなかったのですが学習データの約3%はrandomized policy ( \pi_{rand}) によって生成されています. CausEは, この少量の \mathcal{S}_tを明示的に活用することで精度向上を目指します.

また, 各手法について, Optimizerはエポック数1024のMini-batch Momentum, 学習率には初期値を0.1としたexponential decayを施しています. CausEの \mathcal{S}_t \mathcal{S}_cで得られたembeddingsの乖離に対してかける正則化はL1 normを用いました(L2 normよりも若干良かった).

実験結果

評価方法は前回と同様で, 各手法の25,000エポックの内, validation dataの評価が最も良かったエポックのtest dataに対するMSEの相対改善度を性能評価として用います. ただし, test dataのratingの平均で全埋めしたときのMSEを相対改善度のベースラインとしています. ([Bonner+ 2018]と同じような評価方法です.)

validation時の評価指標は, NaiveなMatix FactorizationはNaiveな評価指標を, Propensity Matrix FactorizatioはUnbiasedな評価指標をそれぞれ対応する指標として用いました. CausEに関しては, validationの際も少量の S_tを用いています.

結果は次の通りです. Best on Valはvalidation dataの評価指標が最も良かったときのエポックのtest dataに対する性能で, Best on Testは, test dataを対する真の性能が最も良かったときのエポックのtest dataに対する性能です. この2つの差が大きい場合, validationに使っていた評価指標が信頼できるものではなかった, ということです.

Best on Val Best on Test
Naive MF - 66.05 % 18.71 %
Propensity MF 16.22 % 19.31 %
Causal Embeddings (CausE) 23.02 % 23.11 %

Naive MFとPropensity MFは, 前回記事と同様の結果を掲載しています. CausEは, 少量の \mathcal{S}_tを明示的に活用することで精度を向上できていることがわかります. CausEの学習の様子は次のようになりました. 青い線は \mathcal{S}_tに対するNaiveな損失ですが, 学習データの3%ほどから計算しているのでエポック間でのブレが大きくなっています.

f:id:usaito:20190416073151p:plain
CausEのlearning curve

さいごに

これまではexplicit feedback recommendationのBias除去の話題について触れてきました. これからはimplicit feedback recommendationに関する同様の話題についても触れていきたいです.

参考

[Bonner+ 2019] Stephen Bonner and Flavian Vasile. Causal embeddings for recommendation: An Extended Abstract. arXiv:1904.05165. 2019.
[Bonner+ 2018] Stephen Bonner and Flavian Vasile. Causal embeddings for recommendation. In Proceedings of the 12th ACM Conference on Recommender Systems, pages 104– 112. ACM, 2018.
[Schnabel+ 2016] Tobias Schnabel, Adith Swaminathan, Ashudeep Singh, Navin Chandak, and Thorsten Joachims. Recommendations as treatments: Debiasing learning and evaluation. In Proceedings of the 33rd International Conference on International Conference on Machine Learning - Volume 48, ICML’16, pages 1670–1679, 2016.