Bayesian Probabilistic Matrix FactorizationをJuliaで実装
概要
リブセンスさんの以下記事で実装されている、Bayesian Probabilistic Matrix Factorizationに興味が湧いたので実装してみました。
とはいえ全く同じことをやっても面白くないので、プラスαとしてJuliaで実装しPythonと速度を比較してみます。
Bayesian Probabilistic Matrix Factorizationに関しては、ここでは詳しく数式は見ませんが、Probabilistic Matrix Factorizationという手法を拡張したものになります*1。
論文は以下になります。
Bayesian Probabilistic Matrix Factorization using Markov Chain Monte Carlo
検証
コードは以下です*2。
Bayesian_Probabilistic_Matrix_Factorization.ipynb · GitHub
データはおなじみのMovie Lensを使います。
細かいところは実際のコードを見てもらうとして*3、Juliaのコードは以下のように実行しました。
また今回は簡易のため、学習データと精度検証用データを分ける等の処理はしていないです*4。
#時間計測 function main() #hyper hyper parameter beta_0 = 2 mu_0 = 0 D = 10 #分解後の潜在変数の次元 nu_0 = D tmp_matrix = rand(D, D) W_0 = Float64.(Int.(inv(tmp_matrix) * tmp_matrix .>= 0.9)) #逆行列かけて型変換をして単位行列に変換 alpha = 2 T = 50 #サンプリング回数 N = user #ユーザー数 M = item #アイテム数 R = df#今回は簡易的な実験のため学習と精度検証用のデータはわけない #学習 BPMF = BayesianProbabilisticMatrixFactorization.BayesianProbabilisticMatrixFactorizationModel(D, alpha, beta_0, mu_0, nu_0, W_0, N, M, T, R, [], []) BayesianProbabilisticMatrixFactorization.fit(BPMF) #予測 R_P = BayesianProbabilisticMatrixFactorization.predict_all(BPMF, 10) return R_P end @time R_P = main()
計算時間は、
209.223477 seconds (212.96 M allocations: 304.652 GiB, 9.71% gc time)
という結果です。
一方でPythonは、
np.random.seed(1) import time start = time.time() D = 10 n_sample = 50 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 = 10 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) time.time() - start
こんな感じで、リブセンスさんと全く一緒のコードです。
比較のためパラメータ等条件はJuliaでの実行と同じです。
計算時間はというと、
4599.483886003494
という結果でした。
Juliaのほうが20倍くらい早いですね、Juliaやばいです。
最後に一応RMSEも見ておくとJuliaが、
sqrt(mean((df[:, 3] - BayesianProbabilisticMatrixFactorization.predict(df, R_P)) .^ 2))
で実行し、
0.7190330856378937
でした。
一方Pythonは、
#精度確認 new_data = np.array(df) ret = np.zeros(new_data.shape[0]) for i in range(len(ret)): ret[i] = p[new_data[i, :][0], new_data[i, :][1]] np.sqrt(np.mean(pow((new_data[:, 2] - ret), 2)))
で実行、
0.71530658114895218
とほぼ同等な結果が出ました*5。
まとめ
Juliaはやっぱ早いですね！
もはやデータサイエンスに実務で関わる人で使わない理由はあまり無い気がします。
また今回数式を自分で読んで自分で実装する力はかなり大事だなと言うこと改めて感じました*6。
次はRobust Bayesian Matrix Factorisationという、Bayesian Matrix Factorisationを変分推論（variational inference）でやっている論文を見つけたのでこちらをJuliaで実装していきたいと思います*7。
それでは間違い等ありました指摘お願いいたします。
