こんにちは、リブセンスでデータサイエンティストをしている北原です。前回に続きコンテキストを扱えるFactorization Machines(FM)をモデルとした、Bayesian Personalized Ranking(BPR)(以下ではBPR-FMと略)を紹介します。今回はBPR-FMのモデルパラメータ推定の実装と実務適用の話をします。実装にはJuliaを使います。モデルやアルゴリズムの詳細については下記の記事をご参照ください。
Factorization MachinesによるBayesian Personalized Ranking
まずはおさらいとしてBPR-FMについて簡単に紹介します。
BPRはユーザーごとの好みの順位を学習するアルゴリズムのフレームワークで、暗黙的評価データを利用したレコメンデーションで使われます。微分可能であれば好みを表す関数は扱う問題に合わせて選択することができるため、この関数にコンテキストを扱えるFMを使うことができます。本記事ではこれをBPR-FMと呼んでいます。
BPR-FMを使うことで、暗黙的評価データでもコンテキストを利用したレコメンデーションが可能になります。例えば、弊社の求人レコメンデーションであれば、勤務地や職種、年収などの情報を利用できるようになるため、ユーザーによりマッチした求人をおすすめできるようになることが期待できます。
評価データ数を$s$、IDなども含むユーザーコンテキスト数（特徴量数）を$n^u$、アイテムコンテキスト数を$n^i$、因子数を$k$とし、評価値$\mathbf{y} \in {\mathbb R}^{s}$、評価と対応するアイテムコンテキストデータ${\bf X}^i \in{\mathbb R}^{s \times n^i}$（１サンプルのコンテキストは$\bf{x}^i$と表記）、ユーザーコンテキストデータ${\bf X}^u \in{\mathbb R}^{s \times n^u}$、アイテムコンテキストに対応するモデルパラメータ${\bf w}^i \in {\mathbb R}^{n^i}$、${\bf V}^{i} \in {\mathbb R}^{n^i \times k}$、ユーザーコンテキストのモデルパラメータ${\bf V}^{u} \in {\mathbb R}^{n^u \times k}$を使って、BPRで使うFMを
とします。
モデルパラメータは
を使って、アルゴリズム
データ
まず学習データについて考えてみましょう。
サンプリングをして$(u, i, j)$が得られたときにモデルパラメータの更新計算に使うデータは$({\bf x}^u, {\bf x}^i, {\bf x}^j)$という組み合わせになります。同一ユーザが２回以上評価することもあれば、同一アイテムが複数回評価されることもあるため、コンテキストデータを評価数分保持するのは冗長です。そこで、ユーザーコンテキストについてはユニークユーザー数$s^u$、アイテムコンテキストについてはユニークアイテム数$s^i$だけ保持します。ユーザーやアイテムを表す$u, i, j$は連番IDとし、ユーザーコンテキスト${\bf X}^u \in {\mathbb R}^{s^u \times n^u}$、アイテムコンテキスト${\bf X}^i \in {\mathbb R}^{s^i \times n^i}$はそのID順に格納します。評価済みアイテム集合$I^{+}$はリストなどに格納しユーザーIDごとにアクセスしやすくしておきます。
そうすると、学習データは以下のようなものになります。
例ではコンテキストデータはすべて0か1の２値になっていますが、実数であればそれ以外の数値も使えます。ただし、コンテキストを使わない時にMatrix Factorizationと同じ挙動をさせたいのであれば、ユーザーIDやアイテムIDに対応するカラムの値は0もしくは1にします。なお、複数段階評価のデータを使う場合、ユーザーIDと相対評価の高いアイテムID、相対評価の低いアイテムIDの三つ組データを使います。弊社で運用しているコードも三つ組データですが、サンプルコードが煩雑になるので本記事ではユーザーと評価アイテムのペアのみを扱います。
モデルパラメータ更新ごとに１評価分のデータから学習データ$({\bf x}^u, {\bf x}^i, {\bf x}^j)$を作成します。学習される可能性のある評価アイテムと未評価アイテムの組み合わせを全てデータとして保持しておくのは冗長で不必要にメモリを消費してしまうので、１更新ごとに学習データを作成します。例えば、評価データから$u=1, i=2$が得られ未評価アイテム$j=4$がサンプリングされたときの学習データ$({\bf x}^u, {\bf x}^i, {\bf x}^j)$は以下のように作成されます。コンテキストを利用しないBPRでは評価データが学習対象ですが、本記事の方法でFMを利用する場合は評価データから得られたコンテキストデータの部分が学習対象になります。
ここからコードを使った説明をします。まず本記事で使うパッケージを事前に読み込んでおきます。
# 本記事で使うパッケージ using BenchmarkTools using DataFrames using Distributions using LinearAlgebra using Random using SparseArrays
本記事で使うサンプルデータは以下のようになります。
# 評価データ S = [1 1; 1 2; 1 3; 2 1; 2 2; 3 5; 3 6; 4 4; 4 5; 5 3; 5 4; 5 6; 6 1; 6 3] # ユーザー集合 U = collect(1:maximum(S[:, 1])) # 全アイテム集合 Iᵃ = collect(1:maximum(S[:, 2])) # ユーザー別評価済みアイテム集合 I⁺ = [d[:, :i] for d in groupby(DataFrame(u = S[:,1], i = S[:, 2]), :u)] # ユーザーコンテキスト Xᵘ = [1 0; 1 0; 0 1; 0 1; 0 1; 1 0] Xᵘ = [Matrix(I, maximum(U), maximum(U)) Xᵘ] # アイテムコンテキスト Xⁱ = [1 0 1; 1 0 0; 1 0 1; 0 1 0; 0 1 0; 0 1 0] Xⁱ = [Matrix(I, maximum(Iᵃ), maximum(Iᵃ)) Xⁱ] # スパース行列に変換 Xᵘ = sparse(Xᵘ) Xⁱ = sparse(Xⁱ)
レコメンデーションで使うコンテキストデータは非常にスパースなので通常はスパース行列として扱います。本記事でもスパース行列として扱いますが、慣れないうちは通常の行列の形で確認しておいたほうが理解が捗るかもしれません。
julia> convert(Matrix, sparse(Xᵘ)) 6×8 Array{Int64,2}: 1 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 1 1 0 julia> convert(Matrix, sparse(Xⁱ)) 6×9 Array{Int64,2}: 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 1 0
実装
本節ではBPR-FMのシンプルな実装を紹介します。最初から計算効率を重視したコードにすると計算式やアルゴリズムとの対応がわかりにくくなるので、まずはアルゴリズムをそのままコードに落とします。
計算効率を重視しないとはいえ式($\ref{eq:bpr_fm}$)〜($\ref{eq:dydv_u}$)をそのまま計算するのは冗長なので少し置き換えをします。頻出する積和演算を
と置きます。そうすると式($\ref{eq:bpr_fm}$)〜($\ref{eq:dydv_u}$)は
となります。事前に$q^i_f, q^j_f, q^u_f$を計算しておくことで計算量が削減されますし、計算式も簡潔になり実装しやすくなります。
BPR-FMの実装は以下のようになります。評価ユーザー$u$とアイテム$i$のペアをランダムサンプリングする部分は、評価データの学習順序のシャッフルに置き換えています。また、収束判定はせず指定回数後に終了します。式($\ref{eq:bpr_fm_q}$)〜($\ref{eq:dydv_u_q}$)の計算には行列演算を利用していますが、計算式ほぼそのままの実装になっていることがわかると思います。なお、このコードは$X^i, X^u$がスパース行列のときだけでなく通常の行列でも動作します。ただし、このような書き方だと現在のJulia(version 1.4.2)ではあまり最適化が効かないので遅いしメモリも消費するので注意が必要です。
function bpr_fm1(S, Xⁱ, Xᵘ, Iᵃ, I⁺, k, α, λ, n_itr) s = size(S)[1] nᵘ = size(Xᵘ)[2] nⁱ = size(Xⁱ)[2] wⁱ = rand(Normal(0.0, 0.001), nⁱ) Vⁱ = rand(Normal(0.0, 0.001), nⁱ, k) Vᵘ = rand(Normal(0.0, 0.001), nᵘ, k) for it in 1:n_itr for s_ in shuffle(1:s) u, i = S[s_, :] j = sample(setdiff(Iᵃ, I⁺[u])) xⁱ = Xⁱ[i, :] xʲ = Xⁱ[j, :] xᵘ = Xᵘ[u, :] qⁱ = Vⁱ' * xⁱ qʲ = Vⁱ' * xʲ qᵘ = Vᵘ' * xᵘ yᵘⁱʲ = wⁱ ⋅ (xⁱ .- xʲ) + 0.5 * (sum(qⁱ.^2 .- qʲ.^2) - sum((xⁱ.^2 .- xʲ.^2)' * Vⁱ.^2)) + qᵘ' * (qⁱ .- qʲ) ∂lnσ = 1 / (exp(yᵘⁱʲ) + 1) wⁱ .+= α .* (∂lnσ .* (xⁱ .- xʲ) - λ .* wⁱ) Vⁱ .+= α .* (∂lnσ .* (xⁱ * (qⁱ .+ qᵘ)' - xʲ * (qʲ + qᵘ)' - Vⁱ .* (xⁱ.^2 .- xʲ.^2)) - λ .* Vⁱ) Vᵘ .+= α .* (∂lnσ .* (xᵘ * (qⁱ .- qʲ)') - λ .* Vᵘ) end end return (wⁱ, Vⁱ, Vᵘ) end
高速化
前節の実装はあまり速くないので、前節のコードを書き換えて高速化を行います。高速化は２種類紹介します。一つはスパースデータに特化させることで無駄な計算をなくす方法です。もう一つは行列演算をループに置き換えることで冗長なメモリ確保を回避する方法です。ここではほぼ確実に有効な方法を紹介しますが、Juliaはちょっとした書き方の違いで最適化が効くことがあるので今回紹介するものよりよい方法があるかもしれません。
まず、スパースデータに特化させる方法から紹介します。前節の書き方だとゼロ成分も計算してしまうので遅くなります。そこで非ゼロ成分のみを計算することで効率化を図ります。データの書き換えが生じなければゼロ成分はずっとゼロなので、うまい書き方をすればコンパイラの最適化でゼロ成分の計算を回避できるのかもしれませんが、ここでは自分で非ゼロ成分を取り出して計算します。
非ゼロ成分のみ計算することで高速化を図ったコードは以下のようになります。アイテムコンテキストの非ゼロ成分はアイテム$i$と$j$とでは異なるため、個別に成分を取り出して計算しています。また、転置してカラム単位でユーザーあるいはアイテムにアクセスできるようにしています。Juliaのスパース行列はCompressed Sparse Column(CSC)なのでカラム方向にデータが配置されており、カラム単位で行列にアクセスするほうが高速なためです。BPR-FMではユーザーあるいはアイテム単位でコンテキストデータにアクセスしますが、本記事で扱っているコンテキストデータでは各行に１ユーザー分あるいは１アイテム分のデータを入れているのでそのままだと遅くなります。
function bpr_fm2(S, Xⁱ, Xᵘ, Iᵃ, I⁺, k, α, λ, n_itr) s = size(S)[1] nᵘ = size(Xᵘ)[2] nⁱ = size(Xⁱ)[2] wⁱ = rand(Normal(0.0, 0.001), nⁱ) Vⁱ = rand(Normal(0.0, 0.001), nⁱ, k) Vᵘ = rand(Normal(0.0, 0.001), nᵘ, k) # カラム単位でレコードにアクセスできるよう転置 Xⁱᵗ = sparse(Xⁱ') Xᵘᵗ = sparse(Xᵘ') for it in 1:n_itr for s_ in shuffle(1:s) u, i = S[s_, :] j = sample(setdiff(Iᵃ, I⁺[u])) # 非ゼロ成分のインデックスと値 (iᵘ, xᵘ) = findnz(Xᵘᵗ[:, u]) (iⁱ, xⁱ) = findnz(Xⁱᵗ[:, i]) (iʲ, xʲ) = findnz(Xⁱᵗ[:, j]) qⁱ = sum(Vⁱ[iⁱ, :] .* xⁱ, dims = 1) qʲ = sum(Vⁱ[iʲ, :] .* xʲ, dims = 1) qᵘ = sum(Vᵘ[iᵘ, :] .* xᵘ, dims = 1) yᵘⁱʲ = wⁱ[iⁱ] ⋅ xⁱ - wⁱ[iʲ] ⋅ xʲ + 0.5 * (sum(qⁱ.^2 .- qʲ.^2) - (sum((xⁱ.^2)' * Vⁱ[iⁱ, :].^2) - sum((xʲ.^2)' * Vⁱ[iʲ, :].^2))) + qᵘ ⋅ (qⁱ .- qʲ) ∂lnσ = 1 / (exp(yᵘⁱʲ) + 1) wⁱ[iⁱ] .+= α * ∂lnσ .* xⁱ wⁱ[iʲ] .-= α * ∂lnσ .* xʲ wⁱ .-= α * λ .* wⁱ Vⁱ[iⁱ, :] .+= α * ∂lnσ .* (xⁱ * (qⁱ .+ qᵘ) - Vⁱ[iⁱ, :] .* xⁱ.^2) Vⁱ[iʲ, :] .-= α * ∂lnσ .* (xʲ * (qʲ .+ qᵘ) - Vⁱ[iʲ, :] .* xʲ.^2) Vⁱ .-= α * λ .* Vⁱ Vᵘ[iᵘ, :] .+= α * ∂lnσ .* (xᵘ * (qⁱ .- qʲ)) Vᵘ .-= α * λ .* Vᵘ end end return (wⁱ, Vⁱ, Vᵘ) end
次に、冗長なメモリ確保を回避する方法を紹介します。ここまでの実装では行列演算を多用していましたが、この行列演算が問題になります。Juliaでは行列演算時に暗黙的にメモリが確保されるため、RやPythonと違って行列演算を利用すると実行速度が遅くなることが多いです。view()関数やLinearAlgebraパッケージのmul!()関数などを使って冗長なメモリ確保を回避する方法はあるのですが、シンプルかつ確実な方法は行列演算をループに置き換えて要素単位で計算する方法なので今回はこれを使います。
行列演算をループに置き換えることで高速化を図ったコードは以下のようになります。繰り返し使われる配列や行列は事前に初期化してメモリ領域を確保しておき、要素への代入を行うことで暗黙的なメモリ確保を回避しつつメモリ領域を使い回すことができます。後はひたすらループに置き換えます。
function vx!(q, V, i, x, k) for f in 1:k q[f] = 0.0 for (i_, x_) in zip(i, x) q[f] += V[i_, f] * x_ end end end function y_uij(qⁱ, qʲ, qᵘ, wⁱ, Vⁱ, iⁱ, iʲ, xⁱ, xʲ, k) y = 0.0 for (i_, x_) in zip(iⁱ, xⁱ) v² = 0.0 for v in Vⁱ[i_, :] v² += v^2 end y += wⁱ[i_] * x_ - 0.5 * x_^2 * v² end for (i_, x_) in zip(iʲ, xʲ) v² = 0.0 for v in Vⁱ[i_, :] v² += v^2 end y -= wⁱ[i_] * x_ - 0.5 * x_^2 * v² end y_ = 0.0 for f in 1:k y_ += qⁱ[f]^2 - qʲ[f]^2 end y += 0.5 * y_ for f in 1:k y += qᵘ[f] * (qⁱ[f] - qʲ[f]) end return y end function bpr_fm3(S, Xⁱ, Xᵘ, Iᵃ, I⁺, k, α, λ, n_itr) s = size(S)[1] nᵘ = size(Xᵘ)[2] nⁱ = size(Xⁱ)[2] wⁱ = rand(Normal(0.0, 0.001), nⁱ) Vⁱ = rand(Normal(0.0, 0.001), nⁱ, k) Vᵘ = rand(Normal(0.0, 0.001), nᵘ, k) # メモリ確保 qⁱ = zeros(k) qʲ = zeros(k) qᵘ = zeros(k) αλ = α * λ Xⁱᵗ = sparse(Xⁱ') Xᵘᵗ = sparse(Xᵘ') for it in 1:n_itr for s_ in shuffle(1:s) u, i = S[s_, :] j = sample(setdiff(Iᵃ, I⁺[u])) (iᵘ, xᵘ) = findnz(Xᵘᵗ[:, u]) (iⁱ, xⁱ) = findnz(Xⁱᵗ[:, i]) (iʲ, xʲ) = findnz(Xⁱᵗ[:, j]) vx!(qⁱ, Vⁱ, iⁱ, xⁱ, k) vx!(qʲ, Vⁱ, iʲ, xʲ, k) vx!(qᵘ, Vᵘ, iᵘ, xᵘ, k) yᵘⁱʲ = y_uij(qⁱ, qʲ, qᵘ, wⁱ, Vⁱ, iⁱ, iʲ, xⁱ, xʲ, k) ∂lnσ = 1 / (exp(yᵘⁱʲ) + 1) α∂lnσ = α * ∂lnσ for (i_, x_) in zip(iⁱ, xⁱ) wⁱ[i_] += α∂lnσ * x_ for f in 1:k Vⁱ[i_, f] += α∂lnσ * x_ * (qⁱ[f] + qᵘ[f] - Vⁱ[i_, f] * x_) end end for (i_, x_) in zip(iʲ, xʲ) wⁱ[i_] -= α∂lnσ * x_ for f in 1:k Vⁱ[i_, f] -= α∂lnσ .* x_ * (qʲ[f] + qᵘ[f] - Vⁱ[i_, f] * x_) end end for i_ in eachindex(wⁱ) wⁱ[i_] -= αλ .* wⁱ[i_] end for i_ in eachindex(Vⁱ) Vⁱ[i_] -= αλ * Vⁱ[i_] end for (i_, x_) in zip(iᵘ, xᵘ) for f in 1:k Vᵘ[i_, f] += α∂lnσ * x_ * (qⁱ[f] - qʲ[f]) end end for i_ in eachindex(Vᵘ) Vᵘ[i_] -= αλ * Vᵘ[i_] end end end return (wⁱ, Vⁱ, Vᵘ) end
最後に実行速度を計測してみましょう。サンプルデータを読み込んで、下記の調整パラメータで実行します。
k = 7 α = 0.01 λ = 0.01 n_itr = 30 @benchmark bpr_fm1(S, Xⁱ, Xᵘ, Iᵃ, I⁺, k, α, λ, n_itr) @benchmark bpr_fm2(S, Xⁱ, Xᵘ, Iᵃ, I⁺, k, α, λ, n_itr) @benchmark bpr_fm3(S, Xⁱ, Xᵘ, Iᵃ, I⁺, k, α, λ, n_itr)
私の環境では以下のような結果になりました。今回のサンプルデータではスパースデータに特化させることで1.6倍ほど速くなり、行列演算をループにすることでさらに3.3倍ほど速くなっているようです。どちらの高速化もそれほど大きな効果はないように見えますが、実際のレコメンデーションで使うデータはサンプルデータよりはるかにスパースかつサイズも大きいので効果が高いです。私の場合、bpr_fm1()のようなコードでモデルやアルゴリズムの確認をし、実データで動かす段階でbpr_fm2()のようなコードに書き換え、実行時間が苦しくなるとbpr_fm3()のようなコードに書き換えています。
julia> @benchmark bpr_fm1(S, Xⁱ, Xᵘ, Iᵃ, I⁺, k, α, λ, n_itr) BenchmarkTools.Trial: memory estimate: 8.91 MiB allocs estimate: 68500 -------------- minimum time: 3.675 ms (0.00% GC) median time: 4.371 ms (0.00% GC) mean time: 5.451 ms (20.10% GC) maximum time: 13.909 ms (54.20% GC) -------------- samples: 916 evals/sample: 1 julia> @benchmark bpr_fm2(S, Xⁱ, Xᵘ, Iᵃ, I⁺, k, α, λ, n_itr) BenchmarkTools.Trial: memory estimate: 4.44 MiB allocs estimate: 51710 -------------- minimum time: 2.151 ms (0.00% GC) median time: 2.541 ms (0.00% GC) mean time: 3.298 ms (22.75% GC) maximum time: 17.467 ms (77.42% GC) -------------- samples: 1514 evals/sample: 1 julia> @benchmark bpr_fm3(S, Xⁱ, Xᵘ, Iᵃ, I⁺, k, α, λ, n_itr) BenchmarkTools.Trial: memory estimate: 1.09 MiB allocs estimate: 12071 -------------- minimum time: 675.084 μs (0.00% GC) median time: 778.361 μs (0.00% GC) mean time: 982.032 μs (18.45% GC) maximum time: 10.302 ms (89.67% GC) -------------- samples: 5076 evals/sample: 1
実務での適用事例
今回BPR-FMの適用対象となったサービスでは、各ユーザーに合った求人を見つけてもらいやすくするために、応募や求人詳細ページの閲覧を学習データとして使った求人のレコメンデーションを行っています。会員登録時に勤務地や職種に関する希望条件を登録してもらう仕組みになっているので、レコメンデーションにはこれらの希望条件を反映させる必要があります。１ユーザーあたりの評価データ数が少ないため、今まで使っていたコンテンツベースフィルタリングや通常のFactorization Machinesだけではユーザーの希望条件を適切に反映させたレコメンデーションができませんでした。そこで、ルールベースのフィルタリングを併用して希望条件を満たす求人のみをレコメンドする方法をとっていました。
ところが、ルールベースで希望条件を満たす求人のみをレコメンドすると機会損失が発生する可能性があることがわかりました。希望条件外でも希望勤務地の近隣であれば勤務に支障のないユーザーもいますし、業務内容と職種の関係が明瞭でないこともあるので、希望条件に指定していない求人に応募するケースがあります。調べてみると、勤務地についてはおよそ１割、職種についてはおよそ３割が希望条件を満たしていない求人への応募でした。希望条件を満たす求人のみしかレコメンドしないとこれらの応募を取りこぼしやすくなります。
希望条件に近い求人をレコメンドできれば機会損失を抑制できるので、これを実現するために今回はBPR-FMを使うことにしました。希望条件外でも応募されやすい勤務地や職種を調べ、それをルールベースのフィルタリングに組み込む方法でも希望条件に近い求人をレコメンドできます。しかし、ルールが複雑になり管理や今後の拡張が難しくなるので今回は避けました。
BPR-FMならば、コンテキストを扱えるので希望条件をモデルに入れることができ、未評価アイテムのサンプリングによって希望条件を学習できます。高評価扱いされる応募や閲覧が多い求人は希望条件に合うもしくは希望条件外であっても許容範囲のものが多くなります。一方で、サンプリングされて低評価扱いで学習されるのは希望条件に合わない求人が多くなります。そのため、両者の関係が学習されて、希望条件に合うもしくは希望条件に近い求人がレコメンドされやすくなります。ただし、職種は希望条件を満たしていない応募・閲覧が多く、これだけではこの仕組が機能しにくいため、希望条件に合わない未評価求人をサンプリングすることで希望条件に合う求人をレコメンドしやすくしています。
このBPR-FMを実際のサービスに導入したところ、アルゴリズムの差し替えだけで、従来と比較してレコメンデーション経由の応募や閲覧はおよそ３割増えました。導入直後のA/Bテストのときだけでなくその後も安定的に高いKPIが維持されています。
まとめ
今回はBPR-FMのモデルパラメータ推定の実装と実務での適用事例を紹介しました。数式をほぼそのままコードに落とした実装に加えて、スパースデータへの特化と冗長なメモリ確保回避によって高速化した実装についても紹介しました。また、１ユーザーあたりの評価数が少ない場合でもユーザーの希望条件に近い求人をレコメンドできるようにするためにBPR-FMを導入したことや、それによって実サービスにおいてKPIを向上させることができたことについても紹介しました。
(2020年8月31日加筆)
