こんにちは、リブセンスでデータサイエンティストをしている北原です。今回は5年ぶりにBPMF（Bayesian Probabilistic Matrix Factorization）を扱います。5年前は論文の内容をそのままギブスサンプリングで実装しましたが、今回は同じモデルをStanで実装します。BPMFのポイントは因子行列の分散の扱いにあります。今回もBPMFの特徴がわかりやすくなるようにPMF（Probabilistic Matrix Factorization）の説明も行います。利用言語はRとStanです。なお、元論文や過去記事では分散ではなく分散の逆数である精度でモデルが表現されていましたが、本記事ではStanでの実装を想定し分散を使って表現しています。
BPMFのギブスサンプリングによる実装については下記の過去記事をご覧ください。
PMF
まず、PMF（Probabilistic Matrix Factorization）の簡単な説明から行います。
PMFはMF（Matrix Factorization）を確率モデルで表現したものです。MFと同じく協調フィルタリングで利用されます。PMFは、ユーザー数を$N$、アイテム数を$M$、因子数を$D$とし、評価行列を$\mathbf{R} \in \mathbb{R}^{N \times M}$、分解後のユーザー因子行列を$\mathbf{U} \in \mathbb{R}^{D \times N}$、アイテム因子行列を$\mathbf{V} \in \mathbb{R}^{D \times M}$としたとき
と表現されるモデルです。ただし、協調フィルタリングでの利用を想定するとデータはスパースなので、欠損データは学習に利用しません。
PMFのポイントは、以下のように、分解後の因子行列をシンプルなものに限定しているところにあります。
- 因子同士が無相関
- いずれの因子も分散が等しい
要するに、因子が直交した理想的な分解ができ、かつ、因子のばらつきも同程度となることを仮定していることになります。この条件を満たしやすいデータであればよいですが、そうでなければデータにフィットしにくくなる可能性があります。
PMFのStan実装
それでは、まずPMFを実装してみましょう。
Mnih et al. (2007)のPMFでは分散に適切な値を指定することになっていますが、Stanであれば分散に弱情報事前分布を設定することで容易に推定ができます。そこで、本記事でも分散は推定することにします。実用上も、分散パラメータが推定できると、パラメータ調整の手間が必要なくなるので便利です。
まず、今回学習させるサンプルデータを作成します。Rコードは以下の通りです。因子間に弱い相関(0.3)があり、因子の平均や標準偏差にはややばらつきがあるデータとしています。
library(tidyverse) library(cmdstanr) library(mvtnorm) library(SparseM) library(doParallel) registerDoParallel(detectCores()) options(mc.cores = parallel::detectCores()) cmdstan_version() # サンプルデータ生成のパラメータ set.seed(1) N <- 60 # ユーザー数 M <- 50 # アイテム数 D <- 10 # 因子数 mu <- rnorm(n = D, mean = 0, sd = 1) # 因子の平均 sig <- rnorm(n = D, mean = 1, sd = 0.5) # 因子の標準偏差 rho <- 0.3 # 因子間の相関 missing_rate <- 0.9 # 欠損させる割合 make_UVR <- function(N, M, D, mu, rho, sig, missing_rate) { S <- diag(sig) L <- matrix(rep(rho, D^2), nrow = D) diag(L) <- 1 Sig <- S %*% L %*% S U <- t(rmvnorm(n = N, mean = mu, sigma = Sig)) V <- t(rmvnorm(n = M, mean = mu, sigma = Sig)) R <- t(U) %*% V # スパースにするため0を入れる for (i in 1:N) { R[i, sample(1:M, M * missing_rate)] <- 0 } list(U = U, V = V, R = R) } uvr <- make_UVR(N, M, D, mu, rho, sig, missing_rate) U <- uvr$U V <- uvr$V R <- uvr$R
PMFのStanコードとRの呼び出しコード、出力結果は以下のとおりです。推定結果の量が多いので、推定値と実際の値とのRMSE(Root Mean Squared Error)を計算して出力するようにしています。
// pmf.stan
data{
int<lower=0> N;
int<lower=0> M;
int<lower=0> D;
int<lower=0> N_ra; // 非ゼロ成分の数
vector[N_ra] RA; // 値
array[N_ra] int<lower=1, upper=N> IA; // 行インデックス
array[N_ra] int<lower=1, upper=M> JA; // 列インデックス
}
parameters {
matrix[D, N] U;
matrix[D, M] V;
real<lower=0> sig;
real<lower=0> sig_u;
real<lower=0> sig_v;
}
model {
for (i in 1:N_ra) {
RA[i] ~ normal(U[:, IA[i]]' * V[:, JA[i]], sig);
}
to_vector(U) ~ normal(0, sig_u);
to_vector(V) ~ normal(0, sig_v);
sig ~ cauchy(0, 5);
sig_u ~ cauchy(0, 5);
sig_v ~ cauchy(0, 5);
}
generated quantities {
vector[N_ra] RA_est;
for (i in 1:N_ra) {
RA_est[i] = normal_rng(U[:, IA[i]]' * V[:, JA[i]], sig);
}
}
mod_pmf <- cmdstan_model("pmf.stan")
fit_pmf <- mod_pmf$sample(
data = list(N = N,
M = M,
D = D,
N_ra = length(R_coo@ra),
RA = R_coo@ra,
IA = R_coo@ia,
JA = R_coo@ja)
)
fit_pmf$summary("RA_est") %>%
mutate(v = R_coo@ra) %>%
rename(est = mean) %>%
summarise(rmse = sqrt(mean((est - v)^2)))
# A tibble: 1 × 1 rmse <dbl> 1 1.29
BPMF
BPMF（Bayesian Probabilistic Matrix Factorization）は、分解後の因子行列をより一般的なものにしています(Salakhutdinov et al. 2008)。
- 因子の相関を分散共分散行列$\Sigma$で表現
- 因子の平均$\mu$を表現
- 因子の平均と分散共分散行列の事前分布にNormal-Wishart分布を設定
要するに、因子行列にPMFほど強い仮定を置かない代わりに、事前分布を設定しています。仮定を緩和するとパラメータ推定に使える情報が少なくなり推定が不安定になるので、それを補うため事前分布を設定しています。
今回はStanを使うので、Stanコードにしやすい形でモデルを表現します。具体的には、自由度パラメータ$\nu$、対称正定値行列のパラメータ$\mathbf{W}$をもつ逆Wishart分布を$\mathcal{W}^{-1}(\nu, \mathbf{W})$で表し、因子行列の平均ベクトル$\mathbf{\mu}_U$、$\mathbf{\mu}_V$と、分散共分散行列$\mathbf{\Sigma}_U$、$\mathbf{\Sigma}_V$を導入して
とします。PMFの因子行列は単変量正規分布で表現されていましたが、BPMFでは多変量正規分布で表現されいます。オリジナルのBPMFとほぼ同じですが、平均と分散共分散それぞれに事前分布を設定しているところや、逆Wishart分布を使っているところが違います。Wishart分布ではなく逆Wishart分布を使っているのは、精度行列ではなく分散共分散行列を使っているためです。本記事では、モデルの性質に大きな影響を与えない定数パラメータを$\mathbf{\mu}_0 = \mathbf{0}$、$\beta_0 = 1$、$\nu_0 = D + 1$、$\mathbf{W}_0 = \mathbf{I}_D$として扱います。なお$\mathbf{I}_D$は$D$次元の単位行列です。
PMFと比較してかなり複雑なモデルになっているので、大規模データには向いていません。特に、Stan実装では計算時間やメモリ使用量が膨大になります。そのため、実運用で使う場合は、フルスクラッチで実装する必要があります。しかし、毎回フルスクラッチで実装しているとモデルの検証に非常に時間がかかってしまいます。Stan実装でデータや課題に応じたモデルを検証し、フルスクラッチ実装を実運用に使うといったように、用途に応じて使い分けるのがよいと思われます。
BPMFのStan実装
それではStanでBPMFを実装してみましょう。サンプルデータはPMFと同じものを使います。また、呼び出しコードもPMFとほぼ同じなので省略します。
BPMFのStanコードと出力結果は以下のとおりです。
// bpmf.stan
data{
int<lower=0> N;
int<lower=0> M;
int<lower=0> D;
int<lower=0> N_ra; // 非ゼロ成分の数
vector[N_ra] RA; // 値
array[N_ra] int<lower=1, upper=N> IA; // 行インデックス
array[N_ra] int<lower=1, upper=M> JA; // 列インデックス
}
parameters {
matrix[D, N] U;
matrix[D, M] V;
vector[D] mu_u;
vector[D] mu_v;
cov_matrix[D] Sig_u;
cov_matrix[D] Sig_v;
real<lower=0> sig;
}
model {
for (i in 1:N_ra) {
RA[i] ~ normal(U[:, IA[i]]' * V[:, JA[i]], sig);
}
for (i in 1:N) {
U[, i] ~ multi_normal(mu_u, Sig_u);
}
for (j in 1:M) {
V[, j] ~ multi_normal(mu_v, Sig_v);
}
sig ~ cauchy(0, 5);
mu_u ~ multi_normal(rep_vector(0, D), Sig_u);
mu_v ~ multi_normal(rep_vector(0, D), Sig_v);
Sig_u ~ inv_wishart(D + 1, identity_matrix(D));
Sig_v ~ inv_wishart(D + 1, identity_matrix(D));
}
generated quantities {
vector[N_ra] RA_est;
for (i in 1:N_ra) {
RA_est[i] = normal_rng(U[:, IA[i]]' * V[:, JA[i]], sig);
}
}
# A tibble: 1 × 1 rmse <dbl> 1 0.718
まとめ
今回は、PMFとBPMFのStan実装を紹介しました。PMFでは因子行列に強めの条件を課していますが、BPMFでは条件を緩和することで表現能力を高めています。一方で、条件を緩和すると推定が不安定になるので、BPMFでは(逆)Wishart分布を事前分布とすることで安定化を図っています。今回はStanを利用することで、PMFとBPMFのいずれもシンプルなコードで容易に実装できることを示しました。
なお、Stanのようにギブスサンプリングを使わない場合は、必ずしも共役事前分布を使う必要はありません。実際に、分散共分散行列の事前分布に(逆)Wishart分布を使うと問題があることを過去記事にて紹介しました。LKJ相関行列を使ったBPMFについては別記事にて紹介します。
