因果推論で推薦システムを問い直す(学習アルゴリズム編)

はじめに

以前, こちらに本記事の評価指標編を書きました. 今回は, 同様の問題が推薦アルゴリズムの学習時にも発生し得ることを指摘し, その解決方法について議論します. 評価指標編を読んでいただいている方は, 重複する内容も多いのですんなり読んでいただけると思います.

目次

Toy Example

問題のイメージを持っていただくために, 評価指標編と同様の例を載せておきます.

今, Horror Lovers・Romance Loversというユーザー属性とHorror・Rommance・Dramaという3つのジャンルのみが存在するシンプルな映画推薦の系を考えます. さらに各ユーザーがそれぞれの映画に対して付すRatingがユーザー属性と映画ジャンルの組み合わせのみによって決まる(つまりRatingの付き方は高々6通り)とします. この時, 予測1 ( \hat{Y}_1) と予測2 ( \hat{Y}_2)(例えば予測1は, Horror LoversがRomance映画に対してRating 1という予測を出力するという見方です)の性能を自分たちが所持している推薦ログデータからMAEを用いて評価し, その評価(損失)に基づいて予測を最適化することを考えます.

f:id:usaito:20190416063747p:plain
[Schnabel+ 2016]のFigure 1より

ここで重要なのが評価に用いるログデータの生成過程です. この例では, ログデータが右上の過去の推薦方策と名付けられている行列のように各ユーザー属性と映画ジャンルの組み合わせごとに異なる推薦確率を持った方策によって集められたログデータであるとします. これにより, 例えば3200個のRatingが含まれるログデータについては各ユーザーと映画ジャンルの観測回数は図のようになり, 要素ごとの観測回数がかなり異なることがわかります. この観測回数の不均衡はひとえにこのデータが各ユーザー属性と映画ジャンルの組み合わせごとに異なる推薦確率を持った方策によって集められたことに起因します.

f:id:usaito:20190416063816p:plain
[Schnabel+ 2016]のFigure 1より

このような各要素ごとに観測回数(生成確率)が一様ではないデータを推薦システムの評価に使ってしまうと図のように, 真の性能は予測1の方が良い(MAE = 0.67)にも関わらず, ログデータ上で計算されたMAEは予測2の方が良くなってしまうという非常に厄介な問題が生じます. つまり, naiveに損失関数を設計してしまうと, 誤った方向に学習が進んでしまう可能性があるのです.

推薦アルゴリズムの定式化

前章の例で見た問題の原因を探るために, 推薦アルゴリズムの学習を定式化します. 今各ユーザーを u \in \{1, ..., U \} = \mathcal{U}, 各アイテムを i \in \{ 1, ...., I \} = \mathcal{I}とします. また,  Y_{u,i}をユーザー uがアイテム iに対して有する真のRatingとします. 本記事では,  Y_{u,i}がユーザーとアイテムの組み合わせによって決定的な値を取るとします(仮にRatingを確率変数として扱った場合, 以降の議論が同様に成り立つためには, 追加的な仮定が必要になります. ).

推薦アルゴリズムの学習によって達成したいのは, 次のように定義される真の損失関数 \mathcal{L}を最小化するような予測値集合 \hat{Y} = \{ \hat{Y}_{u, i} \}_{(u,i) \in \mathcal{U} \times \mathcal{I}} を得ることです.


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


ここで,  \delta (\cdot, \cdot): \mathbb{R} \times \mathbb{R} \rightarrow \mathbb{R}は真のRatingとそれに対する予測値を入力とする関数で, 例えばMSEを評価指標とするなら  \delta(x, y) = (x - y)^{2} とすれば良いです.

もちろん, 真の損失関数 \mathcal{L}を知ることができるならばそれを用いてアルゴリズムのパラメータを最適化すれば良くなんら困難は生じません. しかし, 私たちは全てのユーザーとアイテムのインタラクションを観測できるわけではなく, 手元にあるログデータから損失関数を推定して, それを信じた上で, 最適化を実行するしかありません.

ここで自分たちが観測しているRatingと観測していないそれを区別するために新たな変数 O を導入します. この Oは2値確率変数であり各ユーザーとアイテムのペアに対応して独立に存在します.  O_{u,i} = 1 ならば, ユーザー uとアイテム iのペアについてのRatingを観測していることを意味します. また,  O_{u,i}の期待値は,  u iに対して有する真のRating  Y_{u,i}が観測される確率を意味します.

以降, どのようにして真の損失関数の代理として妥当な, そして手持ちのログデータから計算可能な損失関数を構築すべきかについて議論します.

Naive Lossに存在するBias

まず最初に紹介するのが, Naiveな損失関数 \hat{\mathcal{L}}_{naive} \left( \hat{Y} \right)です. これは, 観測されているデータについての損失を単純に平均するという非常に単純な推定量と言えますが, 多くの論文ではこのNaiveな損失関数の推定量をいかにして最適化するかに力点が置かれているように思います.


\begin{aligned}
\hat{\mathcal{L}}_{naive} \left( \hat{Y} \right) =  \frac{1}{ \sum_{u,i} O_{u,i}  } \sum_{(u,i): O_{u,i} = 1} \delta \left( Y_{u,i}, \hat{Y}_{u,i} \right)
\end{aligned}


このNaiveな推定量 O_{u,i} = 1つまり, ログデータとして観測されているRatingに対する損失を単純に平均しています. ごくごく自然な損失関数の推定方法だと思いますが, 実はログデータの Y_{u,i}が観測される確率が一様ではない場合, このNaiveな損失関数の推定量 \hat{\mathcal{L}}_{naive}は真の損失関数 \mathcal{L} に対してBiasを持つことが次のように示されます.

今分母の \sum_{u, i} O_{u, i} を任意の正の定数 Nで置き換えた上で,  \hat{\mathcal{L}}_{naive}の確率変数 Oについての期待値をとります.


\begin{aligned}
 \mathbb{E}_O  \left[  \hat{\mathcal{L}}_{naive} \left( \hat{Y} \right) \right]  
& = \mathbb{E}_O  \left[  \frac{1}{ N } \sum_{(u,i): O_{u,i} = 1} \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right) \right] \\
 & = \mathbb{E}_O  \left[  \frac{1}{ N } \sum_{u,i}  O_{u,i} \cdot \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right) \right] \\
 & =  \sum_{u,i} \frac{\mathbb{E}_{O_{u,i}} \left[   O_{u,i} \right]} {N} \cdot \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right)
\end{aligned}


ここで,  \hat{\mathcal{L}}_{naive}\left( \hat{Y} \right) \propto \mathcal{L} \left( \hat{Y} \right) が成り立つためには, 任意の u iついて


\begin{aligned}
\frac{\mathbb{E}_{O_{u,i}} \left[   O_{u,i} \right]} {N} = \frac{C}{U \cdot I} 
\end{aligned}


が成り立つ必要があります ( Cは何かしらの正の定数.). しかし観測確率が一様ではない場合, この条件を満たす定数 Nは存在しません. よって,  \hat{\mathcal{L}}_{naive} \left( \hat{Y} \right) \neq \mathcal{L} \left( \hat{Y} \right) であり, naive推定量は真の評価値に対してBiasを持つ(比例関係にない)ことがわかりました. また, そのBiasは観測確率が一様ではないことに起因することも示唆されました.

UnbiasedなLossの構築

前章では, naiveな推定量がBiasを持ってしまうという悲しい事実を紹介しましたが, 評価指標の構築時と同様にIPS (Inverse Propensity Score)と呼ばれる方法でBiasの除去が可能です. ここで, Propensity Score (傾向スコア) とは u iに対して有する真のRating  Y_{u,i}が観測される確率のことです. 表記の簡略化のため, Propensity Scoreを


\begin{aligned}
P_{u,i} = \mathbb{E} \left[ O_{u,i} \right]
\end{aligned}


と置くと, IPSを用いた推定量 \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right)は, 次のように定義されます.


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


summationが O_{u,i} = 1 で条件付けられたデータに対するものであることから, この \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right) は観測可能なデータから推定可能な統計量です. 少し天下り的に与えてしまった \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right) ですが, その期待値が真の評価値 \mathcal{L} \left( \hat{Y} \right) に一致することを次のように示すことができます.


\begin{aligned}
 \mathbb{E}_O  \left[  \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right) \right] 
 & = \mathbb{E}_O  \left[  \frac{1}{ U \cdot I } \sum_{(u,i): O_{u,i} = 1} \frac{ \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right)} {P_{u,i}} \right] \\
 & = \mathbb{E}_O  \left[ \frac{1}{ U \cdot I } \sum_{u,i}  O_{u,i} \cdot \frac{ \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right)} {P_{u,i}} \right] \\
 & = \frac{1}{ U \cdot I }  \sum_{u,i} \frac{\mathbb{E}_{O_{u,i}} \left[   O_{u,i} \right]} {P_{u,i}} \cdot \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right) \\
 & = \frac{1}{ U \cdot I }  \sum_{u,i}  \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right) \\
 & = \mathcal{L}  \left( \hat{Y} \right)
\end{aligned}


よって Biasの観点で議論すると, 真の損失に対して比例関係を満たさないnaive推定量 \hat{\mathcal{L}}_{naive}ではなくて, 不偏性を満たすIPS推定量を最適化することでアルゴリズムを学習するのが良さそうと言えます.

Propensity Matrix Factorization

[Schnabel+ 2016]は, 真の損失関数に対する不偏性を満たす \hat{\mathcal{L}}_{IPS} を用いてユーザーとアイテムのembeddings ( \{ \theta_u \}_{u \in \mathcal{U}}  ,   \{ \beta_i \}_{i \in \mathcal{I}})を得るというPropensity Matrix Factorizationを提案しました. これは, 次のような式でembeddingsを逐次更新することにより, Unbiasedな勾配に基づいたembeddingsを獲得することができます. (実際には正則化項も入れます.)


\begin{aligned}
\theta_u & \leftarrow \theta_u -  \nabla_{\theta} \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right) \\
\beta_i &  \leftarrow \beta_i -  \nabla_{\beta} \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right)
\end{aligned}

簡易実験

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

実験設定

まず, 各ユーザーとアイテムに対応する20次元のベクトルを独立にサンプリングします.  \lambdaは適当に0.5としました.


\begin{aligned}
 \theta_u \sim \mathcal{N} \left( {\bf 0}, \lambda {\bf I} \right), \; \beta_i \sim \mathcal{N} \left( {\bf 0}, \lambda {\bf I} \right)
\end{aligned}


次にユーザーバイアス b_uとアイテムバイアス b_iを独立にサンプリングします. こちらも \lambdaは0.5に設定しました.


\begin{aligned}
 b_u \sim \mathcal{N} \left( 0, \lambda \right), \; b_i \sim \mathcal{N} \left( 0, \lambda \right)
\end{aligned}


これらのサンプリングされたベクトルやスカラー値を用いてRatingを次のように生成したあとで, 整数値に直し, Ratingの最小値が1, 最大値が5になるような処理を加えました. 末尾のrating_biasは, ratingの全体平均を調整するために入れています. 今回の実験では3に設定しています.


\begin{aligned}
 Y_{u,i} = \theta_u^{\top} \beta_i + b_u + b_i + rating\_bias
\end{aligned}


コードで書くとこのような感じになります.

# theta, beta, user_bias, item_biasからratingを生成.
preference = theta @ beta.T + user_bias + item_bias.T + rating_bias
rate = np.clip(np.round(preference), a_min=1., a_max=5.)

最後に傾向スコアを次のように生成します. ここでは, 真のRatingが高いほどデータが観測されやすいという考えに基づき次の生成過程を用いました. 実験では,  k=0.275, \alpha=0.2に設定しています. これでsparcityがおおよそ4%となりました.


\begin{aligned}
 P_{u,i} = k \cdot \alpha^{  \left( 5 - Y_{u,i} \right) }
\end{aligned}


上記のような過程で生成されたRatingの分布は次の通りです. 学習データのRatingの平均は, である一方で, テストデータのRatingの平均は, です. 大きなズレが生じていることが一目瞭然です.

f:id:usaito:20190416061820p:plain
TrainとTestのRating分布

実験結果

Naiveな損失を使ったMatrix FactorizationとIPSによるUnbiasedな損失を用いたPropensity Matrix Factorizationを比較します.

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

validation時の評価指標は, NaiveなMatix FactorizationはNaiveな評価指標を, Propensity Matrix FactorizatioはUnbiasedなMatrix Factorizationをそれぞれ対応する指標として用いています(評価指標について詳しくはこちら).

結果は次の通りです. 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 %

validationを対応する評価指標で行なった場合, Propensity MFが大きく上回っています. 一方で, test dataに対して最適だったときの性能はさほど大きな差はありません. 損失関数のBiasを取り除くことが性能向上に寄与することに加えて, 評価指標がいかに重要であるかを物語る結果となりました. さらにNaive MFをNaiveな評価指標で評価するという組み合わせは, baselineよりも悪い結果となってしまう可能性があることも示唆されました.

それぞれの学習の様子は次の通りです. Naive MFは途中でtest dataに対する性能が上昇し始めてしまっていることが見えます. 一方で, Propensity Matrix Factorizationの方はほぼ単調にtest dataに対する損失が減少しています. ちなみに, Propensity Matrix Factorizationの損失は観測確率の逆数で重み付けすることもあって分散が大きくなる(≒各エポックで損失の推定がぶれる)ので, 今回はバッチサイズが1024のミニバッチ学習を用いています.

f:id:usaito:20190416061647p:plain
Naive Matrix Factorizationのlearning curve

f:id:usaito:20190416061708p:plain
Propensity Matrix Factorizationのlearning curve

さいごに

本記事では, 簡単な例から入り理論事実と簡易実験により推薦システムの学習をnaiveな損失関数に基づいて行ってしまうと学習過程でBiasが生じること, 傾向スコアを用いたIPS推定量によりそのBiasを除去し予測性能を向上できることを示してきました. 実は, 2018年のRecsysでPropensity Matrix Factorizationを実験的に上回った手法が提案されています. 次回は, その手法の追試を行っていくつもりです.

参考

[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.