BPMFのモデル
BPMFは、PMFにおける\(U\)と\(V\)のパラメータにも事前分布を設定したモデルです。分布の分散を一般的な分散共分散行列\(\Lambda\)とし、事前分布に共役事前分布を使うことでギブスサンプリングを使えるようにしたところが特徴です。
\(p(R|U, V, \alpha)\)はPMFと同じで、\(U\)、\(V\)についてはより一般的な表現
\begin{eqnarray}
p(U| \mu_U, \Lambda_U) &=& \prod_{i=1}^{N} \mathcal{N}(U_i | \mu_U, \Lambda_U^{-1}) \label{eq:bpmf_p_U} \\
p(V| \mu_V, \Lambda_V) &=& \prod_{j=M}^{N} \mathcal{N}(V_j | \mu_V, \Lambda_V^{-1}) \label{eq:bpmf_p_V}
\end{eqnarray}
に拡張されています。\(\Lambda\)がポイントです。これをみるとBPMFはPMFの自然な一般化であることがわかると思います。
\(\mu_U\), \(\Lambda_U\), \(\mu_V\), \(\Lambda_V\)の事前分布は、自由度\(\nu_0\)、\(D \times D\)のスケール行列\(W_0\)のWishart分布\(\mathcal{W}(\Lambda | W_0, \nu_0)\)を使って
\begin{eqnarray}
p(\mu_U, \Lambda_U| \mu_0, \nu_0, W_0) &=&
\mathcal{N}(\mu_U | \mu_0, (\beta_0 \Lambda_U)^{-1}) \mathcal{W}(\Lambda_U | W_0, \nu_0) \label{eq:bpmf_p_muU_lamU} \\
p(\mu_V, \Lambda_V| \mu_0, \nu_0, W_0) &=&
\mathcal{N}(\mu_V | \mu_0, (\beta_0 \Lambda_V)^{-1}) \mathcal{W}(\Lambda_V | W_0, \nu_0) \label{eq:bpmf_p_muV_lamV}
\end{eqnarray}
と設定されます。これはNormal-Wishart分布で、多次元正規分布の共役事前分布であることが知られています。このように設定することで、事後分布\(p(\mu_U, \Lambda_U | U)\)、\(p(\mu_V, \Lambda_V | V)\)が容易に計算できるようになっています。ベイズに馴染みのない方には複雑に見えるかもしれませんが、事前分布に共役事前分布を使って計算を容易にさせるのはよく用いられるテクニックです。
BPMFの利点としては
- 分散\(\Lambda\)も推定されるので正則化の調整が不要であること
- 評価数の少ないユーザーの推定精度が高いこと
が挙げられます。一方で、欠点としては、計算時間が比較的長いことが挙げられます。これらの特徴は階層ベイズを使った機械学習でよく見られるものと同じです。
MFを一般化したFM(Factorization Machine)でもMCMCを使う方法がありますが、こちらでは分散行列を対角化し共役事前分布にNormal-Gamma分布を使っています。分散行列の逆行列演算はBPMFの計算量が多い一因となっているので、因子間の相関が小さい場合はFMと同じやり方でもよいかもしれません。
BPMFの実装
BPMFはギブスサンプリングを使って実装できます。条件付き分布\(p(U | R, V, \mu_U, \Lambda_U)\)、\(p(V | R, U, \mu_V, \Lambda_V)\)、\(p(\mu_U, \Lambda_U | U)\)、\(p(\mu_V, \Lambda_V | V)\)が計算できればギブスサンプリングを使うことができます。ここで共役事前分布がフル活用されます。
それでは条件付き分布を導出してみましょう。ユーザー同士、アイテム同士は条件付き独立と考えると
\begin{eqnarray}
p(U| R, V, \mu_U, \Lambda_U) &=& \prod_{i=1}^N p(U_i | R, V, \mu_U, \Lambda_U) \\
p(V| R, U, \mu_V, \Lambda_V) &=& \prod_{j=1}^M p(V_j | R, U, \mu_V, \Lambda_V)
\end{eqnarray}
と分解できるので、\(U_i\)について考えてみます。ベイズの定理と式\ref{eq:p_R}、\ref{eq:bpmf_p_U}より
\begin{eqnarray}
p(U_i| R, V, \mu_U, \Lambda_U, \alpha)
&\sim& p(R | U_i, V, \alpha) p(U_i | \mu_U, \Lambda_U) \\
&=& \prod_{j=1}^M \bigg[\mathcal{N}(R_{ij} | U_i^T V_j, \alpha^{-1} )\bigg]^{I_{ij}} \mathcal{N}(U_i | \mu_U, \Lambda_U^{-1})
\end{eqnarray}
となります。平均が未知、分散が既知で共役事前分布に正規分布を使っているので、事後分布も正規分布になります。指数部分に着目して整理すると事後分布は
\begin{eqnarray}
p(U_i| R, V, \mu_U, \Lambda_U, \alpha) &=& \mathcal{N}(U_i | \mu_i^{\ast}, {\Lambda_i^{\ast}}^{-1}) \\
\Lambda_i^{\ast} &=& \Lambda_U + \alpha \sum_{j=1}^M [V_j V_j^T]^{I_{ij}} \\
\mu_i^{\ast} &=& {\Lambda_i^{\ast}}^{-1} \bigg( \alpha \sum_{j=1}^M [V_j R_{ij}]^{I_{ij}} + \Lambda_U \mu_U \bigg)
\end{eqnarray}
となることがわかります。\(V_j\)についても同様なので、他のパラメータを固定すれば、\(U\)、\(V\)は正規分布を使ってサンプリングできます。
ハイパーパラメータについては共役事前分布が設定されているので、Normal-Wishart分布を共役事前分布としたときの事後分布の公式が使えて
\begin{eqnarray}
p(\mu_U, \Lambda_U | U, \mu_0, \nu_0, W_0) &=&
\mathcal{N}(\mu_U | \mu_0^{\ast}, {\beta_0^{\ast} \Lambda_U}^{-1})
\mathcal{W}(\Lambda_U | W_0^{\ast}, \nu_0^{\ast}) \\
\mu_0^{\ast} &=& \frac{\beta_0 \mu_0 + N \overline{U}}{\beta_0 + N} \\
\beta_0^{\ast} &=& \beta_0 + N \\
{[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\\
\overline{U} &=& \frac{1}{N} \sum_{i=1}{N} U_i\\
\overline{S} &=& \frac{1}{N} \sum_{i=1}^N U_I U_i^T
\end{eqnarray}
となります。\(\mu_V\)、\(\Lambda_V\)も同様です。これで\(\mu_U\)、\(\Lambda_U\)、\(\mu_V\)、\(\Lambda_V\)も他のパラメータを固定すればサンプリングできることがわかりました。
必要な条件付き分布がわかったので、これらを使ってギブスサンプリングを実装します。実装は以下のようになります。
def bpmf_gibbs_sampling(R, U, V, N, M, D, n_sample):
beta0 = 2
mu0 = 0
nu0 = D
W0 = np.identity(D)
alpha = 2
for t_ in range(n_sample - 1):
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)
mu0_ast =(beta0 * mu0 + N * U_bar) / (beta0 + N)
mu_u = multivariate_normal(mu0_ast, inv((beta0 + N) * lam_u))
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)
mu0_ast = (beta0 * mu0 + M * V_bar) / (beta0 + M)
mu_v = multivariate_normal(mu0_ast, inv((beta0 + M) * lam_v))
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)
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割程度向上していました。思ったより効果があったようです。
