こんにちは、リブセンスでデータサイエンティストをしている北原です。今回は事前分布にLKJ相関分布を利用したBPMF（Bayesian Probalibistic Matrix Factorization）を扱います。元のBPMF(Salakhutdinov et al. 2008)では因子行列の分散共分散行列の事前分布にWishart分布を使っています。しかし、Wishart分布を利用すると推定値にバイアスが生じるなど問題があることが知られています。一方で、このバイアスはLKJ相関分布を利用するなど、いくつかの方法で緩和できることも知られています。Stanを利用すると容易にLKJ相関分布を利用したBPMFを実装できるので、本記事で紹介します。利用言語はRとStanです。
おさらい
本節では、過去記事を引用しながらBPMFや逆Wishart分布の問題、コレスキー分解、LKJ相関分布について説明します。
BPMFは協調フィルタリングでよく用いられるMF（Matrix Factorization）をベイズモデルにしたものです。BPMFの特徴やギブスサンプリング実装、Stan実装については下記記事をご参照ください。
BPMFでは、因子行列の事前分布に（逆）Wishart分布を利用していますが、分散共分散の事前分布に逆Wishart分布を利用すると推定バイアスなどの問題が生じることが知られています。推定バイアスがやっかいなのは、収束状況がよさそうでも推定結果が誤っているところです。推定に失敗していることに気づきにくいので対処が遅れてしまいます。この問題については下記記事をご参照ください。
しかし、分散共分散を分散の対角行列と相関行列に分解し、相関行列にLKJ相関分布を利用すると推定バイアスの問題が緩和されます。LKJ相関分布やLKJ相関分布とともに使われるコレスキー分解については下記記事をご参照ください。
BPMFでも、（逆）Wishart分布を利用することによって問題が生じている可能性が高いため、LKJ相関分布を利用することで推定結果が改善されることが期待できます。Stanを利用すれば、通常のBPMFをLKJ相関分布を利用したモデルにするのは容易です。
BPMF
LKJ相関分布を利用した時との比較のために、通常のBPMFのモデルについて簡単に説明します。
Stanで実装しやすい形にしたBPMFは、ユーザー数を$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}$、因子行列の平均ベクトルを$\mathbf{\mu}_U$、$\mathbf{\mu}_V$、分散共分散行列を$\mathbf{\Sigma}_U$、$\mathbf{\Sigma}_V$とし、自由度パラメータ$\nu$、対称正定値行列のパラメータ$\mathbf{W}$をもつ逆Wishart分布を$\mathcal{W}^{-1}(\nu, \mathbf{W})$で表したとき
と表現できます。LKJ相関分布を利用するときは、逆Wishart分布の部分を変更します。
LKJ相関分布を利用したBPMF
LKJ相関分布をBPMFに導入するやり方はいろいろありますが、ここでは次のようなモデルにします。分散共分散行列$\mathbf{\Sigma}$を、分散を対角成分にもつ対角行列$\mathbf{S}$と相関行列$\mathbf{R}$とに$\Sigma = \mathbf{S} \mathbf{R} \mathbf{S}$と分解し、$\mathbf{R}$のコレスキー分解を$\mathbf{L}$とし、パラメータ$\eta$をもつLKJ相関分布を$\mathcal{LKJ}(\eta)$として
と表現されるモデルとします。違いは$\mathbf{\Sigma}_U$、$\mathbf{\Sigma}_V$の逆Wishart分布の事前分布を、相関行列のコレスキー分解$\mathbf{L}_U$、$\mathbf{L}_V$にLKJ相関分布を設定するように、($\ref{Sig_U}$)、($\ref{Sig_V}$)を($\ref{L_U}$)、($\ref{L_V}$)へと変更したところのみです。
LKJ相関分布を利用したBPMFのStan実装
それでは、実際にLKJ相関分布を利用したBPMFを実装してみましょう。
サンプルデータは以前の記事と同じものを使いますが、Rコードを再掲します。因子行列については、因子同士に弱い相関(0.3)があり、個々の因子の分散は異なるものの平均1となるようにしてデータを生成しています。
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
Stanコードと、Stanの呼び出しおよび推定結果のRMSE(Root Mean Squared Error)を計算するRコードは次の通りです。多変量正規分布やLKJ相関分布ではコレスキー分解に対応した方法を使っています。過去記事でLKJ相関分布を扱った時と同じやり方です。
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;
cholesky_factor_corr[D] L_u; // コレスキー因子(U)
cholesky_factor_corr[D] L_v; // コレスキー因子(V)
vector<lower=0>[D] sig_u;
vector<lower=0>[D] sig_v;
real<lower=0> sig;
}
transformed parameters {
corr_matrix[D] R_u; // 相関行列(U)
corr_matrix[D] R_v; // 相関行列(V)
R_u = multiply_lower_tri_self_transpose(L_u);
R_v = multiply_lower_tri_self_transpose(L_v);
}
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_cholesky(mu_u, diag_matrix(sig_u) * L_u);
}
for (j in 1:M) {
V[, j] ~ multi_normal_cholesky(mu_v, diag_matrix(sig_v) * L_v);
}
mu_u ~ multi_normal_cholesky(rep_vector(0, D), diag_matrix(sig_u) * L_u);
mu_v ~ multi_normal_cholesky(rep_vector(0, D), diag_matrix(sig_v) * L_v);
sig ~ cauchy(0, 5);
sig_u ~ cauchy(0, 5); // 標準偏差(U)
sig_v ~ cauchy(0, 5); // 標準偏差(V)
L_u ~ lkj_corr_cholesky(1); // コレスキー因子のLKJ相関分布(U)
L_v ~ lkj_corr_cholesky(1); // コレスキー因子のLKJ相関分布(V)
}
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_bpmf_lkj <- cmdstan_model("bpmf_lkj.stan") fit_bpmf_lkj <- mod_bpmf_lkj$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_bpmf_lkj$summary("RA_est") %>% mutate(v = R_coo@ra) %>% rename(est = mean) %>% summarise(rmse = sqrt(mean((est - v)^2)))
出力結果は以下の通りです。以前の記事で紹介した逆Wishart分布を使ったモデルと比較すると、表現能力が上がっているのか、RMSEが小さくなっていることがわかります。
# A tibble: 1 × 1 rmse <dbl> 1 0.397
LKJ相関分布を利用することで逆Wishart分布を使った時のような推定バイアスの懸念は小さくなりますが、デメリットもいくつかあります。まず、逆Wishart分布のときより収束しにくくなります。また、計算時間もかかることが多いです。さらに、モデルがうまく作成できても高速化の余地が大きくないので、実運用には使いにくいという問題もあります。レコメンデーション用のモデルとしては、あくまで検証用と割り切って、実運用向けのモデルに推定バイアスなどが生じていないかを確認するのに用いるのがよいと思われます。
まとめ
今回は、LKJ相関分布を利用したBPMFのStan実装を紹介しました。分散共分散行列の事前分布に逆Wishart分布を使うと推定バイアスなどの問題が生じることが知られていますが、LKJ相関分布を利用することでこの問題が緩和されることが期待されます。Stanで実装し簡単なサンプルデータで確認したところ、RMSEが小さくなっており表現能力が向上している可能性があることがわかりました。高速化の余地が乏しいのでレコメンデーションにおける利用用途は限られますが、推定バイアスが生じていないかを確認するのによいモデルだと思われます。推定バイアスの問題が生じにくく、かつ、高速化も可能なモデルについては別記事にて紹介します。
