BPMF(Bayesian Probabilistic Matrix Factorization)によるレコメンド

こんにちは、リブセンスで機械学習関係の仕事をしている北原です。
弊社の転職ナビアプリには求人をレコメンドする機能が実装されていて、求人の好みを回答すると各ユーザーに合った求人がレコメンドされるようになっています。このサービスではいくつかのレコメンドアルゴリズムが使われているのですが、その中にBPMF(Bayesian Probabilistic　Matrix Factorization)というアルゴリズムがあります。基本的な問題をフルベイズで扱っている典型的なベイズ手法なのですが、使いどころが難しいのか、使われているのをあまり見たことがありません。そこで、今回はこのBPMFを紹介しようと思います。

アプリの求人レコメンド

レコメンドに限らず機械学習では、やりたいことや使えるデータの種類、特徴に応じて適切なアルゴリズムを使うことが大事です。BPMFを使った背景として、まず簡単に求人レコメンドの目的とアプリの特徴について説明します。

転職ナビは成功報酬型の転職サイトであるため、掲載求人数が多いという特徴があります。働きたい仕事や会社の雰囲気などは人によってさまざまなので、多くの求人があるほうが自分に合った仕事探しがしやすいです。一方で、スマホでは画面が小さいため検索で多くの求人をチェックするのは大変です。そこで、求人の好みを回答するだけで自分に合った求人がレコメンドされると便利です。

レコメンドの精度はユーザーの回答が多いほど向上するので、使い続けるほど適切な求人がレコメンドされやすくなります。そのため、長く使い続けてもらいたいのですが、一般にアプリユーザーは利用回数が少ない段階で離脱することも多く、弊社アプリも例外ではありません。使い始めやたまにしか使わないユーザーにも高精度なレコメンドを行えるようにする必要があります。

データの特徴として、アクティブユーザー数は比較的少なく、求人数もECサイトの商品数などと比べれば少ないと言えるでしょう。そうなると、使うレコメンドアルゴリズムも大規模データを扱う場合とは異なるものが適している場合があります。今回は計算量が多少多くても構わないので、データが少なくても精度が高くなるものが必要です。それに、できれば調整パラメータが少なくシステムへの導入がしやすいと楽です。

BPMFは今回の目的に合っていて、計算量は多いものの評価数が少ないユーザーの精度が高いという特徴があります。正則化項の係数を調整する必要もありませんし、ユーザーの求人評価データだけを使うので既存のレコメンドアルゴリズムをそのまま置き換えやすく、システムへの導入もしやすいという事情がありました。また、階層ベイズが実システムでどれくらい有効なのかについての技術的興味もありました。そんなわけで、BPMFを使ってみることにしました。

PMF

BPMFはベースとなったPMF(Probabilistic Matrix Factorization)というモデルと比較するとわかりやすいので、まずPMFについて説明しようと思います。PMFはMF(Matrix Factorization)を確率モデルで表現したもので、MFと同じく協調フィルタリングアルゴリズムの一つです。協調フィルタリングやMFについては、検索すると多くの記事が見つかるのでここでは説明を省きます。

PMFでは、ユーザー数をN$N$、アイテム数をM$M$、因子数をD$D$とし、評価行列をR∈RN×M$R\in {\mathbb{R}}^{N\times M}$、分解後のユーザー因子行列をU∈RD×N$U\in {\mathbb{R}}^{D\times N}$、アイテム因子行列をV∈RD×M$V\in {\mathbb{R}}^{D\times M}$としたとき

p(R|U,V,α)p(U|αU)p(V|αV)===∏i=1N∏j=1M[N(Rij|UTiVj,α−1)]Iij　　∏i=1NN(Ui|0,α−1UI)∏j=MNN(Vj|0,α−1VI)(1)(2)(3)$$\begin{array}{}\text{(1)}& p(R|U,V,\alpha )& =& \prod _{i=1}^N\prod _{j=1}^M[\mathcal{N}(R_{ij}|U_i^T{V}_j,{\alpha }^{-1}){]}^{I_{ij}} \text{(2)}& p(U|{\alpha }_U)& =& \prod _{i=1}^N\mathcal{N}(U_i|0,{\alpha }_U^{-1}I)\text{(3)}& p(V|{\alpha }_V)& =& \prod _{j=M}^N\mathcal{N}(V_j|0,{\alpha }_V^{-1}I)\end{array}$$

というモデルを考えます。Iij$I_{ij}$はi$i$番目のユーザーがj$j$番目のアイテムを評価している時のみ1、それ以外は0となる変数で、評価済みのデータのみを扱うようにするために導入されています。N(x|μ,α−1)$\mathcal{N}(x|\mu ,{\alpha }^{-1})$は、平均をμ$\mu$、分散をα−1${\alpha }^{-1}$とする正規分布です。R$R$の各成分はU$U$とV$V$の積を平均とした正規分布で表した素直なモデルになっています。

これをMAP推定すると、正則化項のついた通常のMFと同じ更新式が得られます。MFの確率モデル解釈が得られるので、MFの正則化項の係数がPMFの分散に対応することなどがわかって理解は進みますが、基本的にはMFと同じです。

MFを使うときには因子数や正則化係数を適切な値にしないと精度が上がりません。パラメータ調整は面倒なので、これらのパラメータも推定できるとうれしいです。

また、U$U$やV$V$の分散に注目すると、PMFでは各因子の分散は全て同一で因子間は無相関であることを仮定していることがわかります。因子の分散が全て同じというのは仮定として強すぎる気がします。例えば、因子が年収と福利厚生に関するものだった場合、年収を気にする人のばらつきと福利厚生を気にする人のばらつきは違いそうです。因子同士の関係も、理想的には無相関であってほしいところですが、因子数の選び方によっては相関があってもおかしくありません。

BPMFのモデル

BPMFは、PMFにおけるU$U$とV$V$のパラメータにも事前分布を設定したモデルです。分布の分散を一般的な分散共分散行列Λ$\mathrm{\Lambda }$とし、事前分布に共役事前分布を使うことでギブスサンプリングを使えるようにしたところが特徴です。

p(R|U,V,α)$p(R|U,V,\alpha )$はPMFと同じで、U$U$、V$V$についてはより一般的な表現

p(U|μU,ΛU)p(V|μV,ΛV)==∏i=1NN(Ui|μU,Λ−1U)∏j=MNN(Vj|μV,Λ−1V)(4)(5)$$\begin{array}{}\text{(4)}& p(U|{\mu }_U,{\mathrm{\Lambda }}_U)& =& \prod _{i=1}^N\mathcal{N}(U_i|{\mu }_U,{\mathrm{\Lambda }}_U^{-1})\text{(5)}& p(V|{\mu }_V,{\mathrm{\Lambda }}_V)& =& \prod _{j=M}^N\mathcal{N}(V_j|{\mu }_V,{\mathrm{\Lambda }}_V^{-1})\end{array}$$

に拡張されています。Λ$\mathrm{\Lambda }$がポイントです。これをみるとBPMFはPMFの自然な一般化であることがわかると思います。

μU${\mu }_U$, ΛU${\mathrm{\Lambda }}_U$, μV${\mu }_V$, ΛV${\mathrm{\Lambda }}_V$の事前分布は、自由度ν0${\nu }_0$、D×D$D\times D$のスケール行列W0$W_0$のWishart分布W(Λ|W0,ν0)$\mathcal{W}(\mathrm{\Lambda }|W_0,{\nu }_0)$を使って

p(μU,ΛU|μ0,ν0,W0)p(μV,ΛV|μ0,ν0,W0)==N(μU|μ0,(β0ΛU)−1)W(ΛU|W0,ν0)　　N(μV|μ0,(β0ΛV)−1)W(ΛV|W0,ν0)(6)(7)$$\begin{array}{}\text{(6)}& p({\mu }_U,{\mathrm{\Lambda }}_U|{\mu }_0,{\nu }_0,W_0)& =& \mathcal{N}({\mu }_U|{\mu }_0,({\beta }_0{\mathrm{\Lambda }}_U{)}^{-1})\mathcal{W}({\mathrm{\Lambda }}_U|W_0,{\nu }_0) \text{(7)}& p({\mu }_V,{\mathrm{\Lambda }}_V|{\mu }_0,{\nu }_0,W_0)& =& \mathcal{N}({\mu }_V|{\mu }_0,({\beta }_0{\mathrm{\Lambda }}_V{)}^{-1})\mathcal{W}({\mathrm{\Lambda }}_V|W_0,{\nu }_0)\end{array}$$

と設定されます。これはNormal-Wishart分布で、多次元正規分布の共役事前分布であることが知られています。このように設定することで、事後分布p(μU,ΛU|U)$p({\mu }_U,{\mathrm{\Lambda }}_U|U)$、p(μV,ΛV|V)$p({\mu }_V,{\mathrm{\Lambda }}_V|V)$が容易に計算できるようになっています。ベイズに馴染みのない方には複雑に見えるかもしれませんが、事前分布に共役事前分布を使って計算を容易にさせるのはよく用いられるテクニックです。

BPMFの利点としては

• 分散Λ$\mathrm{\Lambda }$も推定されるので正則化の調整が不要であること
• 評価数の少ないユーザーの推定精度が高いこと

が挙げられます。一方で、欠点としては、計算時間が比較的長いことが挙げられます。これらの特徴は階層ベイズを使った機械学習でよく見られるものと同じです。

MFを一般化したFM(Factorization Machine)でもMCMCを使う方法がありますが、こちらでは分散行列を対角化し共役事前分布にNormal-Gamma分布を使っています。分散行列の逆行列演算はBPMFの計算量が多い一因となっているので、因子間の相関が小さい場合はFMと同じやり方でもよいかもしれません。

BPMFの実装

BPMFはギブスサンプリングを使って実装できます。条件付き分布p(U|R,V,μU,ΛU)$p(U|R,V,{\mu }_U,{\mathrm{\Lambda }}_U)$、p(V|R,U,μV,ΛV)$p(V|R,U,{\mu }_V,{\mathrm{\Lambda }}_V)$、p(μU,ΛU|U)$p({\mu }_U,{\mathrm{\Lambda }}_U|U)$、p(μV,ΛV|V)$p({\mu }_V,{\mathrm{\Lambda }}_V|V)$が計算できればギブスサンプリングを使うことができます。ここで共役事前分布がフル活用されます。

それでは条件付き分布を導出してみましょう。ユーザー同士、アイテム同士は条件付き独立と考えると

p(U|R,V,μU,ΛU)p(V|R,U,μV,ΛV)==∏i=1Np(Ui|R,V,μU,ΛU)∏j=1Mp(Vj|R,U,μV,ΛV)(8)(9)$$\begin{array}{}\text{(8)}& p(U|R,V,{\mu }_U,{\mathrm{\Lambda }}_U)& =& \prod _{i=1}^{N}p(U_i|R,V,{\mu }_U,{\mathrm{\Lambda }}_U)\text{(9)}& p(V|R,U,{\mu }_V,{\mathrm{\Lambda }}_V)& =& \prod _{j=1}^{M}p(V_j|R,U,{\mu }_V,{\mathrm{\Lambda }}_V)\end{array}$$

と分解できるので、Ui$U_i$について考えてみます。ベイズの定理と式1$\text{1}$、4$\text{4}$より

p(Ui|R,V,μU,ΛU,α)∼=p(R|Ui,V,α)p(Ui|μU,ΛU)∏j=1M[N(Rij|UTiVj,α−1)]IijN(Ui|μU,Λ−1U)(10)(11)$$\begin{array}{}\text{(10)}& p(U_i|R,V,{\mu }_U,{\mathrm{\Lambda }}_U,\alpha )& \sim & p(R|U_i,V,\alpha )p(U_i|{\mu }_U,{\mathrm{\Lambda }}_U)\text{(11)}& & =& \prod _{j=1}^M[\mathcal{N}(R_{ij}|U_i^T{V}_j,{\alpha }^{-1}){]}^{I_{ij}}\mathcal{N}(U_i|{\mu }_U,{\mathrm{\Lambda }}_U^{-1})\end{array}$$

となります。平均が未知、分散が既知で共役事前分布に正規分布を使っているので、事後分布も正規分布になります。指数部分に着目して整理すると事後分布は

p(Ui|R,V,μU,ΛU,α)Λ∗iμ∗i===N(Ui|μ∗i,Λ∗i−1)ΛU+α∑j=1M[VjVTj]IijΛ∗i−1(α∑j=1M[VjRij]Iij+ΛUμU)(12)(13)(14)$$\begin{array}{}\text{(12)}& p(U_i|R,V,{\mu }_U,{\mathrm{\Lambda }}_U,\alpha )& =& \mathcal{N}(U_i|{\mu }_i^{\ast },{{\mathrm{\Lambda }}_i^{\ast }}^{-1})\text{(13)}& {\mathrm{\Lambda }}_i^{\ast }& =& {\mathrm{\Lambda }}_U+\alpha \sum _{j=1}^M[V_j{V}_j^T{]}^{I_{ij}}\text{(14)}& {\mu }_i^{\ast }& =& {{\mathrm{\Lambda }}_i^{\ast }}^{-1}(\alpha \sum _{j=1}^M[V_j{R}_{ij}{]}^{I_{ij}}+{\mathrm{\Lambda }}_U{\mu }_U)\end{array}$$

となることがわかります。Vj$V_j$についても同様なので、他のパラメータを固定すれば、U$U$、V$V$は正規分布を使ってサンプリングできます。

ハイパーパラメータについては共役事前分布が設定されているので、Normal-Wishart分布を共役事前分布としたときの事後分布の公式が使えて

p(μU,ΛU|U,μ0,ν0,W0)μ∗0β∗0[W∗0]−1U¯¯¯¯S¯¯¯======N(μU|μ∗0,β∗0ΛU−1)W(ΛU|W∗0,ν∗0)β0μ0+NU¯¯¯¯β0+Nβ0+NW−10+NS¯¯¯+β0Nβ0+N(μ0−U¯¯¯¯)(μ0−U¯¯¯¯)T1N∑i=1NUi1N∑i=1NUIUTi(15)(16)(17)(18)(19)(20)$$\begin{array}{}\text{(15)}& p({\mu }_U,{\mathrm{\Lambda }}_U|U,{\mu }_0,{\nu }_0,W_0)& =& \mathcal{N}({\mu }_U|{\mu }_0^{\ast },{{\beta }_0^{\ast }{\mathrm{\Lambda }}_U}^{-1})\mathcal{W}({\mathrm{\Lambda }}_U|W_0^{\ast },{\nu }_0^{\ast })\text{(16)}& {\mu }_0^{\ast }& =& \frac{{\beta }_0{\mu }_0+N\overline{U}}{{\beta }_0+N}\text{(17)}& {\beta }_0^{\ast }& =& {\beta }_0+N\text{(18)}& {[W_0^{\ast }]}^{-1}& =& W_0^{-1}+N\overline{S}+\frac{{\beta }_{0}N}{{\beta }_0+N}({\mu }_0-\overline{U})({\mu }_0-\overline{U}{)}^T\text{(19)}& \overline{U}& =& \frac{1}{N}\sum _{i=1}N{U}_i\text{(20)}& \overline{S}& =& \frac{1}{N}\sum _{i=1}^N{U}_I{U}_i^T\end{array}$$

となります。μV${\mu }_V$、ΛV${\mathrm{\Lambda }}_V$も同様です。これでμU${\mu }_U$、ΛU${\mathrm{\Lambda }}_U$、μV${\mu }_V$、ΛV${\mathrm{\Lambda }}_V$も他のパラメータを固定すればサンプリングできることがわかりました。

必要な条件付き分布がわかったので、これらを使ってギブスサンプリングを実装します。実装は以下のようになります。

def bpmf_gibbs_sampling(R, U, V, N, M, D, n_sample):
    # 初期値(BPMFの論文と同じ)
    beta0 = 2
    mu0 = 0
    nu0 = D
    W0 = np.identity(D)
    alpha = 2

    for t_ in range(n_sample - 1):
        # sample lam_u
        S_bar = np.sum([np.outer(U[t_, i, :], U[t_, i, :].T) for i in range(N)], axis=0) / N
        U_bar = np.sum(U[t_], axis=0) / N
        W0_ast = inv(inv(W0) + N * S_bar +
                     (beta0 * N / (beta0 + N)) * np.outer(mu0 - U_bar, (mu0 - U_bar).T))
        lam_u = wishart.rvs(df=nu0 + N, scale=W0_ast)

        # sample mu_u
        mu0_ast =(beta0 * mu0 + N * U_bar) / (beta0 + N)
        mu_u = multivariate_normal(mu0_ast, inv((beta0 + N) * lam_u))

        # sample lam_v
        S_bar = np.sum([np.outer(V[t_, :, j], V[t_, :, j].T) for j in range(M)], axis=0) / M
        V_bar = np.sum(V[t_], axis=1) / M
        W0_ast = inv(inv(W0) + M * S_bar +
                     (beta0 * M / (beta0 + M)) * np.outer(mu0 - V_bar, (mu0 - V_bar).T))
        lam_v = wishart.rvs(df=nu0 + M, scale=W0_ast)

        # sample mu_v
        mu0_ast = (beta0 * mu0 + M * V_bar) / (beta0 + M)
        mu_v = multivariate_normal(mu0_ast, inv((beta0 + M) * lam_v))

        # sample U
        for i in range(N):
            V_VT_I = np.sum([np.outer(V[t_, :, j], V[t_, :, j].T)
                             for j in R.getrow(i).nonzero()[1]], axis=0)
            lam_ast_inv = inv(lam_u + alpha * V_VT_I)

            V_R_I = np.sum([V[t_, :, j] * r for j, r in zip(R.getrow(i).nonzero()[1],
                                                            R.getrow(i).data[0])], axis=0)
            mu_ast = lam_ast_inv.dot((alpha * V_R_I + lam_u.dot(mu_u.T)).T)

            U[t_ + 1, i, :] = multivariate_normal(mu_ast, lam_ast_inv)

        # sample V
        for j in range(M):
            U_UT_I = np.sum([np.outer(U[t_ + 1, i, :], U[t_ + 1, i, :].T)
                             for i in R.getcol(j).nonzero()[0]], axis=0)
            lam_ast_inv = inv(lam_v + alpha * U_UT_I)

            U_R_I = np.sum([U[t_ + 1, i, :] * r for i, r in zip(R.getcol(j).nonzero()[0],
                                                                R.getcol(j).data)], axis=0)
            mu_ast = lam_ast_inv.dot((alpha * U_R_I + lam_v.dot(mu_v.T)).T)

            V[t_ + 1, :, j] = multivariate_normal(mu_ast, lam_ast_inv)

    return U, V

同じようなコードをUとVについて書いているので長く感じるかもしれませんが、条件付き分布をそのままコードに落としていることがわかると思います。実際に運用されているものは高速化や安定化のため、もう少し複雑になっていますが基本的には同じです。実際に動作させる関数呼び出しは以下。

呼び出しコード(クリックで展開)

import numpy as np
from numpy.linalg import inv
from numpy.random import multivariate_normal
from scipy.stats import wishart
from scipy.sparse import lil_matrix, coo_matrix

r = np.array([
    [0, 0, 7],
    [0, 1, 6],
    [0, 2, 7],
    [0, 3, 4],
    [0, 4, 5],
    [0, 5, 4],
    [1, 0, 6],
    [1, 1, 7],
    [1, 3, 4],
    [1, 4, 3],
    [1, 5, 4],
    [2, 1, 3],
    [2, 2, 3],
    [2, 3, 1],
    [2, 4, 1],
    [3, 0, 1],
    [3, 1, 2],
    [3, 2, 2],
    [3, 3, 3],
    [3, 4, 3],
    [3, 5, 4],
    [4, 0, 1],
    [4, 2, 1],
    [4, 3, 2],
    [4, 4, 3],
    [4, 5, 3]
])

r = coo_matrix((r[:, 2], (r[:, 0], r[:, 1])))
r = lil_matrix(r)

np.random.seed(1)
N = 5
M = 6
D = 3
n_sample = 300

U = np.zeros((n_sample, N, D))
V = np.zeros((n_sample, D, M))
U[0, :, :] = np.random.rand(N, D)
V[0, :, :] = np.random.rand(D, M)
U, V = bpmf_gibbs_sampling(r, U, V, N, M, D, n_sample)

burn_in = 100
p = np.empty((N, M))
for u in range(N):
    for i in range(M):
        p[u, i] = np.mean(np.sum(U[burn_in:, u, :] * V[burn_in:, :, i], axis=1))

print(p)

A/Bテストのようなしっかりした比較は行えていないものの、BPMFの前に使っていたMFと比較すると、ポジティブ評価の割合は4割程度向上していました。思ったより効果があったようです。

最後に

今回はBPMFというレコメンドアルゴリズムの紹介をしました。最近ではStanやEdwardなどベイズ手法を簡単に使える環境が整ってきたので、実システムでベイズを使う場面も増えるのではないでしょうか。今後はStanなどを使った事例や他のアルゴリズムの紹介もしていこうと思います。

今回はアルゴリズムの紹介をしましたが、実際の業務ではアルゴリズム以外の部分も開発する必要があるので、そちらの開発に工数を取られて、すぐに導入できる簡単なアルゴリズムしか使えないということも多いのではないでしょうか。弊社も似たような案件を個別に開発を行っていたので、ある事業でうまくいった方法を別事業で使おうとすると１から開発が必要でした。このような問題に対処すべく、機械学習を各事業で容易に導入できるようにする基盤を開発するプロジェクトを現在進めています。近いうちに紹介できると思うので乞うご期待。