LIVESENSE Data Analytics Blog

リブセンスのデータ分析、機械学習、分析基盤に関する取り組みをご紹介するブログです。

Factorization MachinesによるBayesian Personalized Rankingの実装

 こんにちは、リブセンスでデータサイエンティストをしている北原です。前回に続きコンテキストを扱えるFactorization Machines(FM)をモデルとした、Bayesian Personalized Ranking(BPR)(以下ではBPR-FMと略)を紹介します。今回はBPR-FMのモデルパラメータ推定の実装の話をします。実装にはJuliaを使います。モデルやアルゴリズムの詳細については下記の記事をご参照ください。

analytics.livesense.co.jp

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}$(1サンプルのコンテキストは$\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を

$$ \begin{eqnarray} \hat{y}({\bf x}^u, {\bf x}^i) &=& \sum^{n^i}_{j=1}{w^i_j x^i_j} + \sum^{n^i}_{l=1}\sum^{n^i}_{j=l+1}{\hat{w}^{i}_{l, j} x^i_l x^i_j} + \sum^{n^u}_{l=1}\sum^{n^i}_{j=1}{\hat{w}^{u,i}_{l, j} x^u_l x^i_j} \label{eq:fm_bpr} \\ \hat{w}^i_{l, j} &=& \sum^k_{f=1}{v^i_{l,f} v^i_{j,f}} \\ \hat{w}^{u,i}_{l, j} &=& \sum^k_{f=1}{v^u_{l,f} v^i_{j,f}} \end{eqnarray} $$

とします。

 モデルパラメータは

$$ \begin{eqnarray} \hat{y}^{u, i, j} &=& \sum^{n^i}_{j=1}{w^i_j (x^i_j - x^j_j)} \nonumber \\ &+& \frac{1}{2} \sum^{k}_{f = 1} \biggl\{ \biggl(\sum^{n^i}_{j = 1} v^{i}_{j,f} x^{i}_{j}\biggr)^2 - \biggl(\sum^{n^i}_{j = 1} v^{i}_{j,f} x^{j}_{j} \biggr)^2 - \sum^{n^i}_{j=1} (v^{i}_{j, f})^2 \bigl( (x^{i}_{j})^2 - (x^{j}_{j})^2 \bigr) \biggr\} \nonumber \\ &+& \sum^{k}_{f=1} \biggl( \sum^{n^u}_{j=1} v^{u}_{j,f} x^{u}_{j} \biggr) \biggl( \sum^{n^i}_{j = 1} v^{i}_{j, f} x^{i}_{j} - \sum^{n^i}_{j = 1} v^{i}_{j, f} x^{j}_{j} \biggr) \label{eq:bpr_fm} \\ \frac{\partial \hat{y}^{u, i, j}}{\partial w^{i}_{j}} &=& x^{i}_{j} - x^{j}_{j} \label{eq:dydw} \\ \frac{\partial \hat{y}^{u, i, j}}{\partial v^i_{j, f}} &=& x^{i}_{j} \biggl( \sum^{n^i}_{l = 1} v^{i}_{l, f} x^{i}_{l} + \sum^{n^u}_{l = 1} v^{u}_{l, f} x^{u}_{l} \biggr) - x^{j}_{j} \biggl( \sum^{n^i}_{l = 1} v^{i}_{l, f} x^{j}_{l} + \sum^{n^u}_{l = 1} v^{u}_{l, f} x^{u}_{l} \biggr) - v^{i}_{j, f} \bigl( (x^{i}_{j})^2 - (x^{j}_{j})^2 \bigr) \label{eq:dydv_i} \\ \frac{\partial \hat{y}^{u, i, j}}{\partial v^u_{j, f}} &=& x^{u}_{j} \biggl( \sum^{n^i}_{l = 1} v^{i}_{l, f} x^{i}_{l} - \sum^{n^i}_{l = 1} v^{i}_{l, f} x^{j}_{l} \biggr) \label{eq:dydv_u} \end{eqnarray} $$

を使って、アルゴリズム

f:id:livesense-analytics:20200720093512p:plain
BPR-FMのアルゴリズム
で計算できます。

データ

 まず学習データについて考えてみましょう。

 サンプリングをして$(u, i, j)$が得られたときにモデルパラメータの更新計算に使うデータは$({\bf x}^u, {\bf x}^i, {\bf x}^j)$という組み合わせになります。同一ユーザが2回以上評価することもあれば、同一アイテムが複数回評価されることもあるため、コンテキストデータを評価数分保持するのは冗長です。そこで、ユーザーコンテキストについてはユニークユーザー数$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ごとにアクセスしやすくしておきます。

 そうすると、学習データは以下のようなものになります。

f:id:livesense-analytics:20200720093836p:plain
BPR-FMの学習データ

例ではコンテキストデータはすべて0か1の2値になっていますが、実数であればそれ以外の数値も使えます。ただし、コンテキストを使わない時にMatrix Factorizationと同じ挙動をさせたいのであれば、ユーザーIDやアイテムIDに対応するカラムの値は0もしくは1にします。なお、複数段階評価のデータを使う場合、ユーザーIDと相対評価の高いアイテムID、相対評価の低いアイテムIDの三つ組データを使います。弊社で運用しているコードも三つ組データですが、サンプルコードが煩雑になるので本記事ではユーザーと評価アイテムのペアのみを扱います。

 モデルパラメータ更新ごとに1評価分のデータから学習データ$({\bf x}^u, {\bf x}^i, {\bf x}^j)$を作成します。学習される可能性のある評価アイテムと未評価アイテムの組み合わせを全てデータとして保持しておくのは冗長で不必要にメモリを消費してしまうので、1更新ごとに学習データを作成します。例えば、評価データから$u=1, i=2$が得られ未評価アイテム$j=4$がサンプリングされたときの学習データ$({\bf x}^u, {\bf x}^i, {\bf x}^j)$は以下のように作成されます。コンテキストを利用しないBPRでは評価データが学習対象ですが、本記事の方法でFMを利用する場合は評価データから得られたコンテキストデータの部分が学習対象になります。

f:id:livesense-analytics:20200720101656p:plain
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}$)をそのまま計算するのは冗長なので少し置き換えをします。頻出する積和演算を

$$ \begin{eqnarray} q^{i}_{f} &=& \sum^{n^i}_{l = 1} v^{i}_{l, f} x^{i}_{l} \label{eq:q_i} \\ q^{j}_{f} &=& \sum^{n^j}_{l = 1} v^{i}_{l, f} x^{j}_{l} \label{eq:q_j} \\ q^{u}_{f} &=& \sum^{n^u}_{l = 1} v^{u}_{l, f} x^{u}_{l} \label{eq:q_u} \end{eqnarray} $$

と置きます。そうすると式($\ref{eq:bpr_fm}$)〜($\ref{eq:dydv_u}$)は

$$ \begin{eqnarray} \hat{y}^{u, i, j} &=& \sum^{n^i}_{j=1}{w^i_j (x^i_j - x^j_j)} \nonumber \\ &+& \frac{1}{2} \sum^{k}_{f = 1} \bigl\{ (q^{i}_{f})^2 - (q^{j}_{f})^2 \bigr\} - \frac{1}{2}\sum^{n^i}_{j=1} \bigl( (x^{i}_{j})^2 - (x^{j}_{j})^2 \bigr) \sum^{k}_{f = 1} (v^{i}_{j, f})^2 \nonumber \\ &+& \sum^{k}_{f=1} q^{u}_{f} ( q^{i}_{f} - q^{j}_{f}) \label{eq:bpr_fm_q} \\ \frac{\partial \hat{y}^{u, i, j}}{\partial w^{i}_{j}} &=& x^{i}_{j} - x^{j}_{j} \label{eq:dydw_q} \\ \frac{\partial \hat{y}^{u, i, j}}{\partial v^i_{j, f}} &=& x^{i}_{j} ( q^{i}_{f} + q^{u}_{f}) - x^{j}_{j} ( q^{j}_{f} + q^{u}_{f} ) - v^{i}_{j, f} \bigl( (x^{i}_{j})^2 - (x^{j}_{j})^2 \bigr) \label{eq:dydv_i_q} \\ \frac{\partial \hat{y}^{u, i, j}}{\partial v^u_{j, f}} &=& x^{u}_{j} ( q^{i}_{f} - q^{u}_{f} ) \label{eq:dydv_u_q} \end{eqnarray} $$

となります。事前に$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

高速化

 前節の実装はあまり速くないので、前節のコードを書き換えて高速化を行います。高速化は2種類紹介します。一つはスパースデータに特化させることで無駄な計算をなくす方法です。もう一つは行列演算をループに置き換えることで冗長なメモリ確保を回避する方法です。ここではほぼ確実に有効な方法を紹介しますが、Juliaはちょっとした書き方の違いで最適化が効くことがあるので今回紹介するものよりよい方法があるかもしれません。

 まず、スパースデータに特化させる方法から紹介します。前節の書き方だとゼロ成分も計算してしまうので遅くなります。そこで非ゼロ成分のみを計算することで効率化を図ります。データの書き換えが生じなければゼロ成分はずっとゼロなので、うまい書き方をすればコンパイラの最適化でゼロ成分の計算を回避できるのかもしれませんが、ここでは自分で非ゼロ成分を取り出して計算します。

 非ゼロ成分のみ計算することで高速化を図ったコードは以下のようになります。アイテムコンテキストの非ゼロ成分はアイテム$i$と$j$とでは異なるため、個別に成分を取り出して計算しています。また、転置してカラム単位でユーザーあるいはアイテムにアクセスできるようにしています。Juliaのスパース行列はCompressed Sparse Column(CSC)なのでカラム方向にデータが配置されており、カラム単位で行列にアクセスするほうが高速なためです。BPR-FMではユーザーあるいはアイテム単位でコンテキストデータにアクセスしますが、本記事で扱っているコンテキストデータでは各行に1ユーザー分あるいは1アイテム分のデータを入れているのでそのままだと遅くなります。

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のモデルパラメータ推定の実装を紹介しました。数式をほぼそのままコードに落とした実装に加えて、スパースデータへの特化と冗長なメモリ確保回避によって高速化した実装についても紹介しました。

 機械学習を実際のサービスで利用して成果を上げるためには問題やデータの特徴に合わせた対応が必要です。BPR-FMを使うことでコンテキストを利用したレコメンデーションができるようになりますが、何も考えずに使うだけではうまくいきません。実際にBPR-FM導入対象となったサービスのデータでは、そのままBPR-FMを使っただけでは期待したレコメンド結果は得られません。

 BPR-FMは運用中のレコメンデーションで生じていたある課題を解決するために導入したものです。問題やデータの特徴に合わせてBPR-FMのアルゴリズムを少し拡張することでその課題を解決でき、結果的にKPIを向上させることができました。次回はどのような課題をどのように解決したのかを紹介しようと思います。