ExpoMFの紹介と実装

  • タイトル:Modeling User Exposure in Recommendation
  • 著者: Dawen Liang, Laurent Charlin, James McInerney, David M. Blei
  • 所属:Columbia University, McGill University
  • 発表会議:WWW 2016

[link]

はじめに

こんにちは、さとうみなとです。2019年最後は前々から読もうと思っていたExpoMFで締めたいと思います。最近userの閲覧確率を考慮したMatrix Factorizationでも解説されたりしていましたが、自分なりに読んで実装まで落とし込みたいと思いました。

論文概要

Implicit-feedback(e.g. クリックした/していない)を用いた推薦システムについて通常ユーザーがconsumeしていないアイテムを含む全てのアイテムが考慮されています。 しかし、ユーザーはすべてのアイテムを知覚することができないので、直感的にはおかしな仮定に思えます。 本論文では、アイテムへのユーザーのexposureを組み込む新しい確率的アプローチを提案しました。

平たく言えば下記のようなモデルを仮定しています。[参考]

yuiクリックされたかどうか=auiアイテムへユーザーがexposeしたかどうか×ruirelevance\underbrace{y_{ui}}_{\text{クリックされたかどうか}} = \underbrace{a_{ui}}_{\text{アイテムへユーザーがexposeしたかどうか}} \times \underbrace{r_{ui}}_{\text{relevance}}

つまり、アイテムへユーザーがexposeする、且つrelevanceがあった場合にはじめて yui=1y_{ui} = 1 となります。逆にrelevanceがあってもexposeされな(i.e. ユーザーがアイテムを認知しな)ければ yui=0y_{ui} = 0 となります。

Exposure Matrix Factorization (提案手法)

Notation

  • uu:ユーザー。1,,U1, \cdots , U
  • ii:アイテム。1,,I1, \cdots , I
  • A=auiA = \\{ a_{ui} \\}:ユーザー uu がアイテム ii にexposeされたかどうかを表す行列。
  • Y=yuiY = \\{ y_{ui} \\}:ユーザー uu がアイテム ii をクリックしたかどうかを表す行列。

Model Description

次のような確率モデルを仮定する。

θuN(0,λθ1IK)βiN(0,λβ1IK)yuiaui=1N(θuTβi,λy1)yuiaui=0δ0\begin{aligned} \boldsymbol{\theta}_{u} & \sim \mathcal{N}\left(\mathbf{0}, \lambda_{\theta}^{-1} I_{K}\right) \\ \boldsymbol{\beta}_{i} & \sim \mathcal{N}\left(\mathbf{0}, \lambda_{\beta}^{-1} I_{K}\right) \\ y_{u i} | a_{u i}=1 & \sim \mathcal{N}\left(\boldsymbol{\theta}_{u}^{\mathsf{T}} \boldsymbol{\beta}_{i}, \lambda_{y}^{-1}\right) \\ y_{u i} | a_{u i}=0 & \sim \delta_{0} \end{aligned}

ここで、δ0\delta_{0}p(yui=0aui=0)=1p\left(y_{u i}=0 | a_{u i}=0\right)=1を表す。 つまり、そもそもユーザーにアイテムがexposeされなければ必ず yui=0y_{ui} = 0 となる。 また、μui\mu_{ui} はexposureの事前確率。
ユーザー uu がアイテム ii にexposeされたかどうかについてベルヌーイ分布を仮定して

p(aui,yuiμui,θu,βi,λy1)=p(auiμui,θu,βi,λy1)×p(yuiaui,μui,θu,βi,λy1)=Bernoulli(auiμui)×(N(yuiθuTβi,λy1)auiI[yui=0](1aui))\begin{aligned} p\left(a_{u i}, y_{u i} | \mu_{u i}, \boldsymbol{\theta}_{u}, \boldsymbol{\beta}_{i}, \lambda_{y}^{-1}\right) & = p\left(a_{u i} | \mu_{u i}, \boldsymbol{\theta}_{u}, \boldsymbol{\beta}_{i}, \lambda_{y}^{-1}\right) \times p\left(y_{u i} | a_{u i}, \mu_{u i}, \boldsymbol{\theta}_{u}, \boldsymbol{\beta}_{i}, \lambda_{y}^{-1}\right) \\ & = \operatorname{Bernoulli}\left(a_{u i} | \mu_{u i}\right) \times \left( \mathcal{N}\left(y_{u i} | \boldsymbol{\theta}_{u}^{\mathsf{T}} \boldsymbol{\beta}_{i}, \lambda_{y}^{-1}\right)^{a_{u i}} * \mathbb{I}\left[y_{u i}=0\right]^{\left(1-a_{u i}\right)} \right) \\ \end{aligned}

以上から対数尤度は

logp(aui,yuiμui,θu,βi,λy1)=logBernoulli(auiμui)+auilogN(yuiθuTβi,λy1)+(1aui)logI[yui=0]\begin{aligned} \log p\left(a_{u i}, y_{u i} | \mu_{u i}, \boldsymbol{\theta}_{u}, \boldsymbol{\beta}_{i}, \lambda_{y}^{-1}\right) = & \log \operatorname{Bernoulli}\left(a_{u i} | \mu_{u i}\right) \\ & + a_{u i} \log \mathcal{N}\left(y_{u i} | \boldsymbol{\theta}_{u}^{\mathsf{T}} \boldsymbol{\beta}_{i}, \lambda_{y}^{-1}\right) \\ & + \left(1-a_{u i}\right) \log \mathbb{I}\left[y_{u i}=0\right] \end{aligned}

となる。 また、今回はexposureを表す方法としてアイテムのpopularityを用いる方法を採用。共役事前分布として下記のベータ分布を仮定する。

μiBeta(α1,α2)\mu_{i} \sim \operatorname{Beta}\left(\alpha_{1}, \alpha_{2}\right)

Inference

論文中ではEMアルゴリズムでの学習方法を紹介。

E-step

Eステップでは、positive feedbackの観測されなかった全てのユーザー、アイテムの組み合わせ に対するexposureの期待値を計算します。 (positive feedbackの観測された u,iu, i のペアは必ず aui=1a_{u i} = 1 なので。)

E[auiθu,βi,μui,yui=0]=aui0,1aui×p(auiyui=0)=p(aui=1yui=0)=p(yui=0aui=1)×p(aui=1)p(yui=0)=p(yui=0aui=1)×p(aui=1)p(yui=0,aui=0)+p(yui=0,aui=1)=p(aui=1)μui×p(yui=0aui=1)p(aui=1)×p(yui=0aui=1)+p(aui=0)=1p(aui=1)=1μui×p(yui=0aui=0)=1=μuiN(0θuTβi,λy1)μuiN(0θuTβi,λy1)+(1μui)\begin{aligned} \mathbb{E}\left[a_{u i} | \boldsymbol{\theta}_{u}, \boldsymbol{\beta}_{i}, \mu_{u i}, y_{u i}=0\right] & = \sum_{a_{ui} \in 0,1} a_{ui} \times p(a_{ui} | y_{ui} = 0)\\ & = p(a_{ui} = 1 | y_{ui} = 0)\\ & = \frac{p(y_{ui} = 0 | a_{ui} = 1) \times p(a_{ui} = 1)}{p(y_{ui} = 0)} \\ & = \frac{p(y_{ui} = 0 | a_{ui} = 1) \times p(a_{ui} = 1)}{p(y_{ui} = 0, a_{ui} = 0) + p(y_{ui} = 0, a_{ui} = 1)} \\ & = \frac{\overbrace{p(a_{ui} = 1)}^{\mu_{ui}} \times p(y_{ui} = 0 | a_{ui} = 1) }{p(a_{ui} = 1) \times p(y_{ui} = 0 | a_{ui} = 1) + \underbrace{p(a_{ui} = 0)}_{=1-p(a_{ui}=1)=1-\mu_{ui}} \times \underbrace{p(y_{ui} = 0 | a_{ui} = 0)}_{=1}}\\ & = \frac{\mu_{u i} \cdot \mathcal{N}\left(0 | \boldsymbol{\theta}_{u}^{\mathsf{T}} \boldsymbol{\beta}_{i}, \lambda_{y}^{-1}\right)}{\mu_{u i} \cdot \mathcal{N}\left(0 | \boldsymbol{\theta}_{u}^{\mathsf{T}} \boldsymbol{\beta}_{i}, \lambda_{y}^{-1}\right)+\left(1-\mu_{u i}\right)} \end{aligned}

M-step

対数尤度の auia_{ui} についてEステップで計算した期待値

pui=E[auiθu,βi,μui,yui]p_{u i} = \mathbb{E}\left[a_{u i} | \boldsymbol{\theta}_{u}, \boldsymbol{\beta}_{i}, \mu_{u i}, y_{u i}\right]

を用いて、各パラメータ θu\boldsymbol{\theta}_uβi\boldsymbol{\beta}_i で偏微分し =0=0 とおいて求めた値で更新する。

[2020/01/05更新追記]
対数尤度は

u,ilogp(aui,yuiμui,θu,βi,λy1)+ulogp(θu)+ilogp(βi)=u,ipuilogN(yuiθuTβi,λy1)+ulogp(θu)+ilogp(βi)+Const.=u,ipui×log12πλy1exp((yuiθuTβi)2λy2)+ulogp(θu)+ilogp(βi)+Const.=u,ipui((yuiθuTβi)2λy2)+ulogp(θu)+ilogp(βi)+Const.\begin{aligned} & \sum_{u,i} \log p\left(a_{u i}, y_{u i} | \mu_{u i}, \boldsymbol{\theta}_{u}, \boldsymbol{\beta}_{i}, \lambda_{y}^{-1}\right) + \sum_u \log p(\boldsymbol{\theta}_u) + \sum_i \log p(\boldsymbol{\beta}_i) \\ & = \sum_{u,i} p_{u i} \log \mathcal{N}\left(y_{u i} | \boldsymbol{\theta}_{u}^{\mathsf{T}} \boldsymbol{\beta}_{i}, \lambda_{y}^{-1}\right) + \sum_u \log p(\boldsymbol{\theta}_u) + \sum_i \log p(\boldsymbol{\beta}_i) + \text{Const.} \\ & = \sum_{u,i} p_{u i} \times \log \frac{1}{\sqrt{2 \pi \lambda_y^{-1}}} \exp \left(- \frac{(y_{ui} - \boldsymbol{\theta}_{u}^{\mathsf{T}} \boldsymbol{\beta}_{i})^2 \lambda_y}{2}\right) + \sum_u \log p(\boldsymbol{\theta}_u) + \sum_i \log p(\boldsymbol{\beta}_i) + \text{Const.}\\ & = \sum_{u,i} p_{u i} \left( - \frac{(y_{ui} - \boldsymbol{\theta}_{u}^{\mathsf{T}} \boldsymbol{\beta}_{i})^2 \lambda_y}{2} \right) + \sum_u \log p(\boldsymbol{\theta}_u) + \sum_i \log p(\boldsymbol{\beta}_i) + \text{Const.} \end{aligned}

これについて

θuu,ipui((yuiθuTβi)2λy2)+θuulogp(θu)=θuu,ipui((yui22yuiθuTβi+(θuTβi)2)λy2)λθθu=iλypui(yuiβiβiβiTθu)λθθu=0\begin{aligned} & \frac{\partial}{\partial \boldsymbol{\theta}_{u}} \sum_{u,i} p_{ui} \left(- \frac{(y_{ui} - \boldsymbol{\theta}_{u}^{\mathsf{T}} \boldsymbol{\beta}_{i})^2 \lambda_y}{2}\right) + \frac{\partial}{\partial \boldsymbol{\theta}_{u}} \sum_u \log p(\boldsymbol{\theta}_u)\\ & = \frac{\partial}{\partial \boldsymbol{\theta}_{u}} \sum_{u,i} p_{ui} \left(- \frac{(y_{ui}^2 -2 y_{ui}\boldsymbol{\theta}_{u}^{\mathsf{T}} \boldsymbol{\beta}_{i} + (\boldsymbol{\theta}_{u}^{\mathsf{T}} \boldsymbol{\beta}_{i})^2) \lambda_y}{2}\right) - \lambda_{\theta} \boldsymbol{\theta}_u \\ & = \sum_{i} \lambda_y p_{ui} (y_{ui}\boldsymbol{\beta}_i - \boldsymbol{\beta}_i\boldsymbol{\beta}_i^{\mathsf{T}} \boldsymbol{\theta}_u) - \lambda_{\theta} \boldsymbol{\theta}_u = 0 \end{aligned}

したがって

θu(λyipuiβiβiT+λθIK)1(iλypuiyuiβi)βi(λyupuiθuθuT+λβIK)1(uλypuiyuiθu)\begin{aligned} &\boldsymbol{\theta}_{u} \leftarrow\left(\lambda_{y} \sum_{i} p_{u i} \boldsymbol{\beta}_{i} \boldsymbol{\beta}_{i}^{\mathsf{T}}+\lambda_{\theta} I_{K}\right)^{-1}\left(\sum_{i} \lambda_{y} p_{u i} y_{u i} \boldsymbol{\beta}_{i}\right)\\ &\boldsymbol{\beta}_{i} \leftarrow\left(\lambda_{y} \sum_{u} p_{u i} \boldsymbol{\theta}_{u} \boldsymbol{\theta}_{u}^{\mathsf{T}}+\lambda_{\beta} I_{K}\right)^{-1}\left(\sum_{u} \lambda_{y} p_{u i} y_{u i} \boldsymbol{\theta}_{u}\right) \end{aligned}

を得る。

[2020/01/06更新追記]
また、各 μi\mu_i ついては事後分布を計算すると、

p(μia1:U,i)p(a1:U,iμi)×p(μi)μia1i(1μi)1a1iμiaUi(1μi)1aUi×Γ(α1+α2)Γ(α1)+Γ(α2)μiα11(1μi)α21μiα1+uaui1(1μi)α2+(Uuaui)1\begin{aligned} p(\mu_i | \boldsymbol{a}_{1:U,i}) & \propto p(\boldsymbol{a}_{1:U,i} | \mu_i) \times p(\mu_i)\\ & \propto \mu_i^{a_{1i}}(1-\mu_i)^{1-a_{1i}} \cdots \mu_i^{a_{Ui}}(1-\mu_i)^{1-a_{Ui}}\times \frac{\Gamma(\alpha_1+\alpha_2)}{\Gamma(\alpha_1) + \Gamma(\alpha_2)}\mu_i^{\alpha_1-1}(1-\mu_i)^{\alpha_2-1}\\ & \propto \mu_i^{\alpha_1 + \sum_u a_{ui}-1}(1-\mu_i)^{\alpha_2 + (U - \sum_u a_{ui})-1} \end{aligned}

この式はベータ分布

Beta(α1+uaui,α2+Uuaui)\text{Beta}\left(\alpha_1 + \sum_u a_{ui}, \alpha_2 + U - \sum_u a_{ui}\right)

に従っている。auia_{ui} は実際には期待値を用いるので puip_{ui} で置き換えて

Beta(α1+upui,α2+Uupui)\operatorname{Beta}\left(\alpha_{1}+\sum_{u} p_{u i}, \alpha_{2}+U-\sum_{u} p_{u i}\right)

の最頻値

μiα1+upui1α1+α2+U2\mu_{i} \leftarrow \frac{\alpha_{1}+\sum_{u} p_{u i}-1}{\alpha_{1}+\alpha_{2}+U-2}

を用いて更新する。

擬似コード

以上を踏まえた学習の擬似コードは下記の通り。

01: input: click matrix YY
02: randomly initialize user factors θ1:U\boldsymbol{\theta}_{1:U}
03: randomly initialize item factors β1:I\boldsymbol{\beta}_{1:I}, exposure priors μ1:I\mu_{1:I}
04: loop for each iteration
05:   Compute expected exposure AA (e-step)
06:   Update user factors θ1:U\boldsymbol{\theta}_{1:U} (m-step)
07:   Update item factors β1:I\boldsymbol{\beta}_{1:I} (m-step)
08:   Update priors μi\mu_{i} (m-step)

Pythonコード

ExpoMFは著者実装があって幸いPython3系でも割と簡単に動かすことができるのですが、せっかくならなるべく数式を追って理解したいので自分でも実装してみたくなりました。

movielensのデータ取得の部分でlightfmの力を拝借しておりますので

pip install lightfm

としてインストールください。
コード中に注釈をつけておりますので擬似コードと照らし合わせてください。

import numpy as np
from lightfm.datasets import fetch_movielens

dataset = fetch_movielens(min_rating=4.0)

iterations = 3
K = 20
lambda_theta = 1e-2
lambda_beta = 1e-2
lambda_y = 1e-2
alpha_1 = 1.0
alpha_2 = 1.0
init_mu = 0.01
U, I = dataset["train"].shape
A = np.zeros(dataset["train"].shape)


# 擬似コード01行目 = click matrix
Y = dataset["train"].toarray()
Y[Y > 0] = 1.0

# 擬似コード02行目 = ユーザーの潜在ベクトルの初期化
theta = np.random.uniform(low=-0.1, high=0.1, size=(U, K)) / K

# 擬似コード03行目 = アイテムの潜在ベクトル、\muの初期化
beta = np.random.uniform(low=-0.1, high=0.1, size=(I, K)) / K
mu = np.ones(I) * init_mu

# 擬似コード04行目 = loop開始
for iteration in range(iterations):
    print(f"CURRENT ITERATION = {iteration+1}")

    # 擬似コード05行目 = Exposureの期待値計算 (e-step)
    for u in range(U):
        for i in range(I):
            if Y[u, i] == 1.0:
                A[u, i] = 1.0
            else:
                n_ui = lambda_y / np.sqrt(2*np.pi) * np.exp(- np.dot(theta[u], beta[i]) ** 2 * lambda_y ** 2)
                A[u, i] = n_ui / (n_ui + (1 - mu[i]) / mu[i])

    # 擬似コード06行目 = ユーザーの潜在ベクトルの更新 (m-step)
    for u in range(U):
        B = np.zeros((K, K))
        a = np.zeros(K)
        for i in range(I):
            B += A[u, i] * beta[i, None] * beta[i, None].T
            a += lambda_y * A[u, i] * Y[u, i] * beta[i]
        B *= lambda_y
        B += lambda_theta * np.eye(K)
        new_theta = np.linalg.solve(B, a)
        theta[u] = new_theta

    # 擬似コード07行目 = アイテムの潜在ベクトルの更新 (m-step)
    for i in range(I):
        B = np.zeros((K, K))
        a = np.zeros(K)
        for u in range(U):
            B += A[u, i] * theta[u, None] * theta[u, None].T
            a += lambda_y * A[u, i] * Y[u, i] * theta[u]
        B *= lambda_y
        B += lambda_beta * np.eye(K)
        new_beta = np.linalg.solve(B, a)
        beta[i] = new_beta

    # 擬似コード08行目 = \muの更新 (m-step)
    mu[:] = (alpha_1 + A.sum(axis=0) - 1) / (alpha_1 + alpha_2 + U - 2)

Y_test = dataset["test"].toarray()
Y_test[Y_test > 0] = 1.0

from sklearn import metrics
predicted = theta.dot(beta.T)
scores = np.zeros(U)
for u in range(U):
    fpr, tpr, thresholds = metrics.roc_curve(Y_test[u], predicted[u])
    scores[u] = metrics.auc(fpr, tpr) if len(set(Y_test[u])) != 1 else 0.0

print(f"test mean auc: {scores.mean()}")

僕の環境での実行結果は

CURRENT ITERATION = 1
CURRENT ITERATION = 2
CURRENT ITERATION = 3
test mean auc: 0.8930167190209362

となりました。

さいごに

今回はExpoMFについて簡単に紹介と実装を行ってみました。数式やコード、説明に間違いございましたらTwitterかコメント欄までお願い致します。


Written by@Minato Sato
Senior Software Engineer - Embedded AI

GitHubTwitterFacebookLinkedIn

© 2023 Minato Sato. All Rights Reserved